text:
# Spatio-temporal patterns in the Hantavirus infection

## I Introduction

Hantaviruses are infectious agents carried by rodents throughout the whole world schmaljohn97 ; mills99a ; mills99b . Some of them are able to cause severe disease in humans, with a mortality rate of around 50%, as in the case of the Hantavirus Pulmonary Syndrome (HPS) caused by the Sin Nombre Virus in the North American Southwest, or the Andes Virus in Patagonia. With few exceptions, each hantavirus is associated with a single primary rodent host species, chronically infected, and infects humans that come into contact with it or its excreta. Sin Nombre Virus is primarily carried by the deer mouse, *Peromyscus maniculatus*, the most numerous mammal in North America, prevalent throughout the region. It was the cause of an outbreak of fatal pulmonary disease in the Southwest of the United States in 1993, when the virus was first isolated and described. Since then, a great effort has been devoted to understand the nature of the virus reservoir, its temporal and spatial dynamics, and its relation to the human population, in an attempt to ultimately identify and predict the risk of the disease.

Needless to say, a complete mathematical description of the dynamics of the biological system, comprising the virus, the mice, the humans and the environment, is a daunting task. The goal of the present investigation is much less ambitious. From the biological complexities we extract a few major components centered on the basic ecological and epidemiological features of the mice population. As the motivation for our analysis we choose two observed characteristics of the disease. Both arise from the fact that environmental conditions strongly affect the dynamics and persistence of the infection. One of them, a temporal characteristic, is the reported observation that the infection can completely disappear from a population of mice if environmental conditions are inadequate, only to reappear sporadically or when conditions change mills99b ; calisher99 ; parmenter99 . The other, a spatial characteristic, is that there are indications of “focality” of the infection in “reservoir” populations mills99b ; kuenzi99 ; as environmental changes occur, these “refugia” yates2001  of the reservoir can expand or contract, carrying the infection to other places.

The model we introduce incorporates the decay by death of the mice population, the spread of the infection through their interaction, the increase by birth and effect of the environment to stabilize the population, and also their movement as a process of diffusion. We begin in Section II by first omitting the last feature (the movement from one location to another), motivating the different dynamical mechanisms, and obtaining some basic results including the observed temporal behavior. We proceed in Section III to perform a spatial extension of the model to include movement and obtain results relating to the refugia. A summary is given in the final section.

## II Basic model of mouse population

We can incorporate the basic ingredients of the biological system in a model of the mouse population only anderson . We suppose that the whole population is composed of two classes of mice, susceptible and infected, represented by $`M_S`$ and $`M_I`$ respectively. Sex and age composition of the population are disregarded in this basic model. The temporal evolution of $`M_S`$ and $`M_I`$ contains two basic ingredients: the contagion of the infection, that converts susceptible into infected, and a population dynamics independent of the infection:

$`{\displaystyle \frac{dM_S}{dt}}`$ $`=`$ $`bMcM_S{\displaystyle \frac{M_SM}{K}}aM_SM_I,`$ (1)
$`{\displaystyle \frac{dM_I}{dt}}`$ $`=`$ $`cM_I{\displaystyle \frac{M_IM}{K}}+aM_SM_I,`$ (2)
where $`M_S`$ and $`M_I`$ are the populations (or densities) of susceptible and infected mice, respectively, and $`M(t)=M_S(t)+M_I(t)`$ is the total population of mice. The motivation for the terms in Eqs. (1,2) follows.

*Births:* $`bM`$ represents births of mice, all of them born susceptible, at a rate proportional to the total density, since all mice contribute equally to the procreation mills99b .

*Deaths:* $`c`$ represents the rate of depletion by death for natural reasons, proportional to the corresponding density. If necessary, separate rates $`c_S`$ and $`c_I`$ could be introduced for the susceptible and infected populations respectively.

*Competition:* $`M_{S,I}M/K`$ represent a limitation process in the population growth, due to competition for shared resources. Each is proportional to the probability of an encounter of a pair formed by one mouse of the corresponding class, susceptible or infected, and one mouse of any class (since every mouse, either susceptible or infected, has to compete with the whole population). $`K`$ is a “carrying capacity,” characterizing in a simplified way the capacity of the medium to maintain a population of mice. Higher values of carrying capacity represent a higher availability of water, food, shelter and other resources that mice can use to thrive murray .

*Infection:* $`aM_IM_S`$ represents the number of susceptible mice that get infected, due to an encounter with an infected (and consequently infectious) mouse, at a rate $`a`$ that we assume constant. More elaborate models could incorporate a density dependence on $`a`$, for example due to an increased frequency of fights, during which contagion occurs through bites, when the density is too high and the population feels overcrowded calisher99 . The infection is chronic, infected mice do not die of it, and infected mice do not lose there infectiousness probably for their whole life mills99b ; kuenzi99 . For these reasons, this single term adequately describes the infection dynamics of the two subpopulations.

The sum of the two equations (1,2) reduces to a single equation for the whole population of logistic form:

$$\frac{dM}{dt}=(bc)M\left(1\frac{M}{(bc)K}\right).$$
(3)
Logistic growth has been observed in laboratory populations of *Peromyscus* terman68 , and is a well established metaphor of the dynamics of a self limitating population murray .

There are four parameters that characterize the system (1,2), viz. $`a`$, $`b`$, $`c`$ and $`K`$. Of these, we will choose $`K`$ as a control parameter of the dynamics, since it is the one that best represents the influence of the environment.

The system (1,2) has four equilibria. Two of them are irrelevant to the present analysis (the null state, which is always unstable, and a state with $`M_I<0`$ for any parameters). The other two equilibria interchange their stability character at a critical value of the carrying capacity, a result that we show in Fig. 1 as a bifurcation diagram. The critical value of the carrying capacity is

$$K_c=\frac{1}{a}\left(\frac{b}{bc}\right).$$
(4)
We can see that the prevalence of the infection can be correlated, through $`K`$, with the diversity of habitats and other ecological conditions. Thus, a scarcity of resources—that is to say, a low value of $`K`$—is accompanied by a lower number of infected mice, as found in field studies such as mills99b ; kuenzi99 ; abbot99 . Moreover, for values of $`K`$ below the threshold $`K_c`$ the number of infected animals is effectively zero, a fact that has also been observed in the field (see for example calisher99 ; parmenter99 ; mills99b ). That is, if temporarily the ecological conditions at a place in the landscape get adverse for the mice (because of a drought, for example) the infection can drop to zero. Correspondingly, when conditions improve again the infection reappears. The density of infected mice can even display a dramatic increase with respect to previous years, if a rare climatic event such as El Niño Southern Oscillation brings enhanced precipitation and the consequent increase in edible resources for the mice. An El Niño event in 1991-1992, precisely, preceded the outbreak of HPS in 1993 in the Southwest glass00 .

Figure 2 shows a simulation of such events, within the context of the present model. A time-dependent carrying capacity is shown in Fig. 2 (top), and the corresponding values of the susceptible and infected mice populations, $`M_S(t)`$ and $`M_I(t)`$ respectively, are displayed in Fig. 2 (bottom). We model the carrying capacity with a yearly sinusoidal behavior to emulate seasonal variations. A period of 20 years is shown, during which the carrying capacity oscillates around a value, sometimes above $`K_c`$ (shown as a horizontal line), sometimes below it. Discontinuities in the carrying capacity, some of which are present in Fig. 2 (top), do not necessarily occur in nature, and appear here because we keep the modeling of $`K(t)`$ at an elementary level, to illustrate the main features of the system. The period marked “a” in Fig. 2 (from years 6 to 8) is characterized by values of $`K`$ below $`K_c`$, and corresponds to very adverse environmental conditions. During these “bad years” the infection level effectively drops to zero, while the population of healthy mice, even if reduced, subsists. A return to “normal” carrying capacities after year 8 produces a very slow recovery of the infected population, which attains again appreciable values after year 11. An extraordinary event on year 17 is marked as “b” in Fig. 2. It corresponds to an increase in the carrying capacity (top), perhaps following an event such as El Niño the year before. These improved environmental conditions are followed by an immediate (if moderate) increase in the population of susceptible mice (bottom, dotted line), and by a slightly delayed outbreak of infection (bottom, full line). An event such as this would appreciably increase the risk for the human population to become infected.

## III Spatially extended model

The range of the deer mice is wide, comprising a diverse landscape with a variety of habitats. This spatial extension and the inhomogeneous way in which it affects local populations can be included in a spatially extended version of the model, where $`M_S`$, $`M_I`$ and $`K`$ become functions of a space variable $`𝐱`$. Diffusive movement of the mice provide an adequate mechanism of transport, since mice of the genus *Peromyscus* are known to hold a home range during most of their adult life, occasionally shifting it to nearby locations, in particular if these are vacant stickel68 ; vessey87 . In principle, different diffusion coefficients should be used for susceptible and infected mice. The observation that juvenile animals are the most mobile calisher99  and that the infection affects mainly adult males mills99a  certainly supports this. We will choose later, however, for the sake of simplicity of the model, to keep both diffusivities equal. The extended model can be written as:

$`{\displaystyle \frac{M_S}{t}}`$ $`=`$ $`f(M_S,M_I)+D_S^2M_S,`$ (5)
$`{\displaystyle \frac{M_I}{t}}`$ $`=`$ $`g(M_S,M_I)+D_I^2M_I,`$ (6)
where $`f`$ and $`g`$ are the r.h.s. of Eqs. (1) and (2) respectively (and contain the specific form of the spatial dependence $`K(𝐱)`$), and we include separate diffusion coefficients $`D_S`$ and $`D_I`$ for the two classes of mice.

The solution of the system (5,6), and even its stationary solution, may impossible to find, analytically, for an arbitrary function $`K(𝐱)`$. We describe below some general considerations about stability, followed by numerical solution for $`𝐱`$-dependent $`K`$.

### III.1 Stability of the extended solutions

Suppose that $`M_S^{}(𝐱)`$ and $`M_I^{}(𝐱)`$ are stationary solutions of Eqs. (5,6), i.e. they are solutions of a Laplace equation with nonlinear, space-dependent sources:

$`^2M_S`$ $`=`$ $`f(M_S,M_I)/D_S,`$ (7)
$`^2M_I`$ $`=`$ $`g(M_S,M_I)/D_I,`$ (8)
found by setting the time derivative of Eqs. (5,6) equal to zero. A perturbation around this equilibrium can be written as:

$`M_S(𝐱,t)`$ $`=`$ $`M_S^{}(𝐱)+u_S(𝐱,t),`$ (9)
$`M_I(𝐱,t)`$ $`=`$ $`M_I^{}(𝐱)+u_I(𝐱,t).`$ (10)
When the two-component vector $`u=(u_S,u_I)`$ describing the perturbation is inserted into the differential equations (5,6), a linearization around the equilibrium solutions yields

$$\frac{u(𝐱,t)}{t}=A(𝐱)u(𝐱,t)+D^2u(𝐱,t),$$
(11)
where $`A(𝐱)`$ is the linearization of the nonlinear terms of Eqs. (5,6) around the equilibrium, viz.,

$$A(𝐱)=\left[\begin{array}{cc}\frac{f}{M_S}& \frac{f}{M_I}\\ \frac{g}{M_S}& \frac{g}{M_I}\end{array}\right]_{\{M_S^{},M_I^{}\}},$$
(12)
and $`D`$ is the $`2\times 2`$ diagonal matrix of the diffusivities.

Solutions of Eq. (11) can be looked for in the form of plane waves,

$$u(𝐱,t)e^{i𝐤𝐱+\lambda t},$$
(13)
which, in Eq. (11), satisfies:

$$[\lambda IA(𝐱)+k^2D]u(𝐱,t)=0,$$
(14)
where $`I`$ is the identity matrix. The nontrivial solutions of Eq. 14) will provide a dispersion relation $`\lambda (k^2)`$, implicitly:

$$det[\lambda IA(𝐱)+k^2D]=0.$$
(15)
In the general situation of $`𝐱`$-dependent $`K`$, it is not possible to proceed further without the knowledge of the equilibria. However, in a system where $`K`$ does not depend on the space variable, an analytic assessment of the stability of the *homogeneous* steady states is possible. We have again two relevant steady states: $`\{M_S^{}=(bc)K,M_I^{}=0\}`$ and $`\{M_S^{}=b/a,M_I^{}=b/a+(bc)K\}`$. The dispersion relations corresponding to each one of these are easily found from Eq. (15). Those corresponding to the first one (the equilibrium with $`M_I^{}=0`$) are shown in Fig. 3. They provide a direct stability criterion. The slopes of the two lines are determined by the diffusion coefficients only, and as such are always negative. It can be seen that one of the temporal eigenvalues is always negative, provided that $`b>c`$, which is, obviously, the sensible case in the biological context since otherwise no positive solutions are found. The other eigenvalue is negative provided that $`K<K_c`$, which is the same stability condition found in the nonextended case. Furthermore, when the state becomes unstable, the fastest growing mode of the perturbation (the one with larger $`\lambda `$) is that with $`k^2=0`$, an homogeneous perturbation. Under such conditions, the perturbation eventually drives the system to the other homogeneous steady state, having a nonzero infected population. In this simple model, hence, there are no spatially dependent instabilities to the homogeneous steady state.

### III.2 Refugia

Certainly, the most interesting situations arise when $`K`$ exhibits a spatial dependency. This is in fact the case in the field, where $`K`$ follows the diversity of the landscape. We have analyzed two cases of this situation, by means of a numerical solution of Eqs. (5,6). The first case is a 1-dimensional system, where the profile displayed by the stationary solutions of the populations is readily accessible. The second one is a 2-dimensional system, intended to provide a more realistic picture of the consequences of the bifurcation.

We consider first a 1-dimensional landscape, consisting of a spot of high carrying capacity ($`K>K_c`$) in the middle of a bigger region of low carrying capacity ($`K<K_c`$). A typical situation is shown in Fig. 4, where vertical lines represent the boundaries between the three zones. From an arbitrary initial condition of the populations, a steady state is attained in which the infected population is concentrated at the spot of higher $`K`$, that constitutes a “refugium.” A “leak” of infection is seen outside the high-$`K`$ region, due to the diffusion. Far from this, the mouse population remains effectively not infected.

In Fig. 5 we show the steady state of a 2-dimensional realization of the system (5,6) on a square grid which simulates a hypothetical landscape by assigning different values to $`K_{ij}`$, the carrying capacity at each site. This is supposed higher along a “river” as can be inferred from the density plots shown. The non-infected population occupies the whole landscape, with a non-homogeneous density. Moreover, as expected from the results of the homogeneous model, for small and moderate values of the diffusion coefficient, the infected population survives in a patchy pattern, only in the regions of high carrying capacity, becoming extinct in the rest. These “islands” of infection become reservoirs of the virus kuenzi99  or “refugia” yates2001 , which are the places of highest risk for human exposure and contagion of the virus. It is also from these refugia that the disease would spread (blurring the patchiness, as observed in mills99b ; abbot99 ) when environmental conditions change. While our model is qualitative at this stage, this is precisely what is observed in the field. We comment in passing that the steady state distribution of neither infected nor susceptible mice reproduces exactly the distribution of the carrying capacity. This is the result of the interaction of diffusion with the nonlinear interactions. Thus, notice in the 1-dimensional representation shown in Fig. 4 that, although the carrying capacity follows a step distribution, the mice populations are not steps. Both $`M_S`$ and $`M_I`$ have diffusive “leaking” , the former exhibiting a dip as one moves out of the region of large capacity. Similarly, in the 2-dimensional case shown in Fig. 5, we see that the peaks of the populations represented by pure white appear at different places for the susceptible and infected. They do not occupy the entire “river” region or follow precisely the peaks of the distribution of the carrying capacity.

## IV Concluding remarks

Two observed characteristics of Hantavirus infection have served as the focus of our present investigation: temporal patterns in the evolution of the population of infected mice, and emergence of spatial features in the landscape of infection, the so-called “refugia.” Our theoretical model, represented by (5,6), incorporates nonlinear terms describing infection transfer between mice populations, a logistic description of their interactions with the environment, and diffusive terms representing their motion over the terrain. We have shown that the combination of these various terms, while simple, naturally predicts the temporal and spatial patterns whose observations have motivated the analysis. Our tools of investigation comprise of analytic stability considerations which result in features such as bifurcation behavior (e.g., Fig. 1) as well as numerical procedures which yield the temporal evolution (e.g., Fig. 2). The spatial extension inherent in our model allows us to analyze the dispersion relation describing in a simplified case departures from stationary states (see Fig. 3) and to deduce more generally the existence of the “refugia” (see Figs. 45).

We are currently in the process of investigating a number of further features of the spread of infection on the basis of the model and techniques explained in the present paper. They include among others: traveling waves which can depict the spread of fronts of infection emanating from the refugia in periods favorable to the propagation of the infection; situations in which the mice are limited in their meanderings to more or less localized regions for territorial reasons but spread the infection when the localized regions overlap; non-diffusive effects in the motion of the mice over the terrain; the effect of stochastic disturbances in the environment; and relevant details of the infection process such as delay effects related to finite incubation periods. The results of these investigations will be reported elsewhere.

###### Acknowledgements.

We acknowledge many discussions with Terry Yates, Bob Parmenter, Fred Koster and Jorge Salazar from which we learnt much regarding the peculiarities of the hantavirus including the observation of refugia. We also thank Greg Glass, Karl Johnson and Luca Giuggioli for discusions. V. M. K. acknowledges a contract from the Los Alamos National Laboratory to the University of New Mexico and a grant from the National Science Foundation’s Division of Materials Research (DMR0097204). G. A. thanks the support of the Consortium of the Americas for Interdisciplinary Science and the hospitality of the University of New Mexico.


format: text

text:
# DEUTERIUM CHEMISTRY IN THE PRIMORDIAL GAS

## 1 Introduction

The formation of H2 and HD molecules in the post-recombination Universe plays a central role in the evolution of gas condensations. Even trace abundances of these molecules strongly affect the cooling properties of the primordial gas which would be otherwise an extremely poor radiator (cooling by Ly-$`\alpha `$ photons is ineffective at temperatures less than $`10^4`$ K).

The chemistry of the early Universe has been investigated in several studies starting with the seminal paper by Lepp & Shull (1983). We mention in particular the work of Puy et al. (1993), Galli & Palla (1998; hereafter GP) and the comprehensive analysis of deuterium chemistry by Stancil, Lepp & Dalgarno (1998; hereafter SLD). The abundances of H2 and HD predicted at low redshift in these studies are of the order of $`10^6`$ and $`10^9`$, respectively, depending on the cosmological model adopted (see Palla, Galli & Silk 1995 for the variation of the chemical abundances with the assumed baryon-to-photon ratio). A comparison of the abundance of HD obtained at $`z=10`$ by GP ($`n[\mathrm{HD}]/n[\mathrm{H}]=1.1\times 10^9`$), and SLD ($`n[\mathrm{HD}]/n[\mathrm{H}]=1.6\times 10^9`$) for the same cosmological model ($`h=0.5`$, $`\mathrm{\Omega }_0=1`$, $`\mathrm{\Omega }_b=0.0367`$) shows satisfactory agreement, although this might not be completely significant since both calculations were based on a compilation of reaction rates obtained basically from the same sources. To repeat the words of SLD, “further studies are needed to reduce the uncertainty in the HD abundance”. One of the goals of this paper is to review the progress made in our understanding of HD chemistry in the three years subsequent to the publication of these studies. We will outline the improvements which have occurred in the meantime, and discuss the remaining uncertainties.

From an observational standpoint, the abundance of atomic deuterium and molecular hydrogen has been measured in several cosmological clouds at redshift $`z2`$–3 in absorption along the line of sight to bright quasars. Deuterium in particular has attracted renewed attention because of the controversy about its abundance in high-redshift Ly-$`\alpha `$ clouds (see Hogan 1998 for an account and a resolution of the controversy). Molecular hydrogen at high redshift was first detected by Levshakov & Varshalovich (1985) in a damped Ly-$`\alpha `$ system at $`z=2.8`$. Since then, the presence of H2 has been confirmed in at least four additional systems (see e.g. Levshakov et al. 2000 and references therein).

In Fig. 1 we compare observational data on D and H2 at high redshift with the corresponding abundances calculated with the standard model of GP that follows the homogenous expansion of the universe. The agreement between the theoretical and observed deuterium abundance is excellent. In the case of H2, one should keep in mind the sensitivity of this molecule to ambient conditions in Ly-$`\alpha `$ clouds, such as the stronger ultraviolet radiation field and the lower dust content with respect to local interstellar medium. These fctors can account for the considerable spread of observed abundances. Despite the complex phenomenology associated with damped Ly-$`\alpha `$ systems and the resulting uncertainty in the interpretation of the results, it is encouraging, to say the least, to witness the emergence of the observational foundations of the highly theoretical discipline of primordial chemistry.

## 2 Chemical Reactions

The formation of HD in the primordial gas follows two main routes, involving a deuteron exchange with H2:

$$\mathrm{D}+\mathrm{H}_2\mathrm{HD}+\mathrm{H},$$
(1)
and

$$\mathrm{D}^++\mathrm{H}_2\mathrm{HD}+\mathrm{H}^+.$$
(2)
Being an isotopic modification of the the most fundamental three-electron interaction, namely, the $`\mathrm{H}+\mathrm{H}_2`$ reaction, reaction (1) has received considerable interest. Thermal rate constants for this reaction have been measured by Ridley, Schulz, & LeRoy (1966), Westenberg & de Haas (1967), Mitchell & LeRoy (1973), and Michael & Fisher (1990) over a wide range of temperatures. Theoretical calculations employing statistical, semiclassical and quantal method have been performed by several groups, and show very good agreement with each other and with the experimental data. The most recent studies are by Zhang & Miller (1989), Michael, Fisher, & Bowman (1990), Mielke et al. (1994), and Charutz, Last, & Baer (1997). These results are compared in the usual Arrhenius plot shown in Fig. 2 (the high-temperature experimental results of Michael & Fisher 1990 are not shown). In this paper we adopt the reaction rate computed by Mielke et al. (1994) with the DMBE surface (Varandas et al. 1987) which predicts much more accurate low-temperature kinetics than other surfaces. Their results agree with the values computed by Michael & Fisher (1990) (with the same surface) within $`20\%`$. We also note that Zhang & Miller (1989) computed state-by-state cross sections for $`vJv^{}J^{}`$ transitions and state-to-state rate constants in the temperature range 200–1000 K. The agreement with the experimental cross section values at $`E1`$ eV (Phillips, Levene, & Valentini 1989) is within a factor $`2`$.

Reaction (2) represents the major source of HD in diffuse interstellar clouds (Dalgarno, Weisheit, & Black 1973). Its rate coefficient is almost constant and close to the Langevin value ($`2\times 10^9`$ cm3 s-1, see Fig. 3). The reaction rate has been measured in the laboratory by Fehsenfeld et al. (1973, at $`T=200`$ and 278 K) using a flowing afterglow technique, and by Henchman, Adams, & Smith (1981, at $`T=205`$ and 295 K) using a variable-temperature selected ion flow tube technique (see also Smith, Adams, & Alge 1982). Gerlich (1982) performed accurate quantum-mechanical calculations of the rate coefficient, which are in excellent agreement with the experimental results of Henchmann et al. (1981) at $`T=295`$ K, but less at $`T=205`$ K (still within a factor $`2`$). There is a discrepancy with the results of Fehsenfeld et al. (1973) at both temperatures.

The destruction of HD occurs mainly via the reverse reactions of (1) and (2),

$$\mathrm{HD}+\mathrm{H}\mathrm{D}+\mathrm{H}_2,$$
(3)
and

$$\mathrm{HD}+\mathrm{H}^+\mathrm{D}^++\mathrm{H}_2.$$
(4)
In general, the rate coefficients for the forwards and reverse chemical reactions are related by the standard thermodynamic expression (e.g. Berry, Rice & Ross 1980)

$$\mathrm{ln}\left(\frac{k_\mathrm{f}}{k_\mathrm{r}}\right)=\frac{\mathrm{\Delta }H^0}{RT}+\frac{\mathrm{\Delta }S^0}{R},$$
(5)
where $`R`$ is the universal gas constant, and $`\mathrm{\Delta }H^0`$, $`\mathrm{\Delta }S^0`$ are the enthalpy and entropy changes. Thus, if the rate $`k_\mathrm{f}`$ is knwon, the rate for the reverse reaction $`k_\mathrm{r}`$ can be obtained directly from eq. (5). From the differences in the zero point vibrational energies of H2 and HD and in the ionization potentials of H and D, one obtains the enthalpy changes for reactions (1) and (2): $`\mathrm{\Delta }H^0/R=412`$ K and $`462`$ K, respectively. The entropy change $`\mathrm{\Delta }S^0/R`$ can be calculated on statistical grounds (see e.g. Flower 2000) and is the same for both reactions: $`\mathrm{\Delta }S^0/R=\mathrm{ln}2`$ for para-H2, $`\mathrm{\Delta }S^0/R=\mathrm{ln}(2/3)`$ for ortho-H2. These values are for reactant and product molecules in their ground states. Both the entropy and the enthalpy changes are modified when rovibrationally excited molecular states are involved (Flower 2000).

Rather than using eq. (5), in this paper we prefer to provide independent fitting formulae for both the direct and reverse reactions that dominate the chemistry of deuterium in the primordial gas. For specific applications, the reader may adopt the reaction rate of the forward (or reverse, if better constrained) reaction, and compute the rate of the reverse (forward) reaction from the thermodynamic relation expressed by eq. (5).

The rate coefficient of reaction (3) has been calculated by Shavitt (1959) using a semiempirical H3 energy surface. Only sparse laboratory data exist for this reaction. At $`T=10^3`$ K, the experimental result of Boato et al. (1956) is within $`20\%`$ from the value calculated by Shavitt (1959). At lower temperatures (720–880 K), the theoretical rate is a factor $`2`$ lower than the experimental data by van Meersche (1951), but the extrapolation of the adopted H3 energy surface introduces significant uncertainty in the results at low temperatures (see Fig. 4). As a challenge to chemical physicists, we recall the words of Shavitt (1959): “it is unfortunate that for a reaction as basically important as the one considered here, the experimental data are so incomplete and inconclusive”.

Since reaction (4) is endothermic by 0.0398 eV (462 K), the removal of HD at low temperatures is reduced by a factor $`\mathrm{exp}(462/T)`$, and this can lead to significant enhancement of the HD/H2 ratio (fractionation). The rate coefficient of this reaction has been measured in the laboratory by Henchmann et al. (1981) at $`T=205`$ and 295 K, and the results are in good agreement with the values obtained by applying the principle of detailed balance to the reverse reaction (2). As in GP, we adopt the rate coefficient calculated by Gerlich (1982) over the temperature range 30–600 K (see Fig. 5).

Finally, the relative abundance of D and D+ is determined by the charge exchange reactions

$$\mathrm{H}^++\mathrm{D}\mathrm{H}+\mathrm{D}^+$$
(6)
and its reverse

$$\mathrm{H}+\mathrm{D}^+\mathrm{H}^++\mathrm{D}.$$
(7)
The cross section for reaction (6) was computed by Matveenko (1974) and Hunter & Kuriyan (1977) for energies from $`10^3`$ to 7.5 eV. The two results differ by a factor $`2`$ at low energies for reasons unclear. Subsequent calculations by Hodges & Breig (1993), and, more recently, by Igarashi & Lin (1999) and Zhao, Igarashi, & Lin (2000), confirm the validity of the results of Hunter & Kuriyan (1977), and improve significantly the accuracy of the calculations around $`10^3`$ eV. Good agreement was also found with the experimental measurements of Newman et al. (1982) and Esry et al. (2000). The cross section of reaction (7) was also calculated by Igarashi & Lin (1999) and Zhao et al. (2000). Unfortunately, no experimental data are available for this reaction.

Watson (1976) estimated the rate of reactions (6) and (7) and Watson, Christensen & Deissler (1978) calculated the rate coefficient on the basis of the cross section obtained by Hunter & Kuriyan (1977). These rates have been widely adopted in studies of deuterium chemistry. Since reaction (6) has a threshold of $`43`$ K, the rate coefficient for reaction (7) is usually obtained by multiplying that of reaction (6) by $`\mathrm{exp}(43/T)`$. The situation has been reanalyzed recently by Wolf Savin (2001) who computed accurate rates for reactions (6) and (7) from the cross sections of Igarashi & Lin (1999) and Zhao et al. (2000). These results are compared to those of Watson et al. (1978) in Figs. 6 and 7.

In Table 1 we collect accurate fitting formulae (computed by us or by the authors quoted in the references) for the chemical reactions discussed in this section. These expressions update and replace the corresponding formulae given in Table 2 of GP.

In addition to these reactions, additional contributions to the formation of HD in the early Universe come from the associative detachment reactions

$$\mathrm{D}+\mathrm{H}^{}\mathrm{HD}+e$$
(8)
and

$$\mathrm{D}^{}+\mathrm{H}\mathrm{HD}+e,$$
(9)
whose rate coefficients however are not known, and can only be estimated from the corresponding H reactions (see SLD). Finally, minor contributions to the formation of HD come from the radiative association reaction

$$\mathrm{H}+\mathrm{D}\mathrm{HD}+h\nu ,$$
(10)
whose rate coefficient was computed by Stancil & Dalgarno (1997), and from reactions involving less abundant deuterated species like HD+ and H2D+ (see SLD for details).

Generally, the chemistry of HD in the primordial gas is dominated by the ion–neutral reactions (2) and (4) in a gas of low density (e.g. the primordial gas before the formation of the first structures), whereas the neutral–neutral reactions (1) and (3) become more important in conditions of high density and temperature (cloud collapse, shocked gas).

## 3 Heating and Cooling

In order to calculate the cooling properties of HD, one should know the population of all rovibrational levels. In steady-state, these are obtained by solving the balance equations

$$x_J\underset{J^{}}{}[R_{JJ^{}}(T_{\mathrm{rad}})+C_{JJ^{}}(n,T_{\mathrm{gas}})]=\underset{J^{}}{}x_J^{}[R_{J^{}J}(T_{\mathrm{rad}})+C_{J^{}J}(n,T_{\mathrm{gas}})],$$
(11)
where $`J`$ and $`J^{}`$ indicate a generic couple of rovibrational levels. The collisional transition probabilities $`C_{JJ^{}}(n,T_{\mathrm{gas}})`$ and $`C_{J^{}J}(n,T_{\mathrm{gas}})`$ are obtained by multiplying the corresponding excitation coefficients, $`\gamma _{JJ^{}}(T_{\mathrm{gas}})`$ and $`\gamma _{J^{}J}(T_{\mathrm{gas}})`$, and the density of the colliding species. The terms $`R_{JJ^{}}`$ and $`R_{J^{}J}`$ are the radiative excitation and de-excitation rates, that can be expressed in terms of the Einstein coefficients $`A_{JJ^{}}`$ and $`B_{JJ^{}}`$,

$$R_{JJ^{}}=\{\begin{array}{cc}A_{JJ^{}}+B_{JJ^{}}u(\nu _{JJ^{}},T_{\mathrm{rad}}),\hfill & J^{}<J\text{,}\hfill \\ B_{JJ^{}}u(\nu _{JJ^{}},T_{\mathrm{rad}}),\hfill & J^{}>J\text{,}\hfill \end{array}$$
(12)
where $`u(\nu _{JJ^{}},T_{\mathrm{rad}})`$ is the energy density of the cosmic background radiation (CBR) per unit frequency at the temperature $`T_{\mathrm{rad}}`$,

$$u(\nu _{JJ^{}},T_{\mathrm{rad}})=\frac{8\pi h\nu _{JJ^{}}^3}{c^2}\left[\mathrm{exp}(h\nu _{JJ^{}}/kT_{\mathrm{rad}})1\right]^1.$$
(13)
In the primordial gas, collisional excitation of HD is dominated by collisions with H, and, to a less extent, with He. The coefficients for inelastic scattering of He–HD were computed by Green (1974) and Schaefer (1990) at temperatures $`T600`$ K, for $`0J3`$ and $`\mathrm{\Delta }J=+1,+2`$. Collisional coefficients for rotational excitation of the system H–HD were usually derived by scaling the He-HD values with the square root of the ratio of the reduced masses of the two systems, $`\gamma _{\mathrm{H}\mathrm{HD}}=(\mu _{\mathrm{He}\mathrm{HD}}/\mu _{\mathrm{H}\mathrm{HD}})^{1/2}\gamma _{\mathrm{He}\mathrm{HD}}`$, where $`(\mu _{\mathrm{He}\mathrm{HD}}/\mu _{\mathrm{H}\mathrm{HD}})^{1/2}=1.51`$ (see e.g. Timmerman 1996, GP). In the words of Timmermann (1996), “this assumption is however an educated guess, and data on H–HD are urgently needed”. These data were recently provided by Flower & Roueff (1999) (for H and H2) and Roueff & Zeippen (1999) (for He), using full quantum-mechanical methods and updated energy potential surfaces to evaluate rovibrational excitation coefficients for collisions of HD with H, H2 and He. The results for H–HD collisions computed by Flower & Roueff (1999) in the temperature range $`100T_{\mathrm{gas}}2000`$ K for $`v2`$ and $`J9`$ for a total of 30 rovibrational levels are available from the CCP7 server http:$`//`$ccp7.dur.ac.uk$`/`$.

The energy levels and the Einstein coefficients $`A_{JJ^{}}`$ of HD were calculated by Abgrall et al. (1982) who considered both dipole and quadrupole transitions and included a large number of rovibrational levels. Since the energy spacing of the rotational levels of HD is quite large, $`E_1/k=128`$ K, $`E_2/k=383`$ K, $`E_3/k=764`$ K, etc., GP computed the HD cooling function with a simple four-level system ($`J=0`$–3) adopting the collisional coefficients of Schaefer (1990). Flower et al. (2000) updated the calculations of GP adopting the collisional rate coefficients of Flower & Roueff (1999) and Roueff & Zeippen (1999). The HD cooling function computed by Flower et al. (2000) is also available from the CCP7 server. A useful approximation in the low-density limit is the expression

$$\mathrm{\Lambda }_{\mathrm{HD}}[n(\mathrm{H}0)]=2\gamma _{10}(T_{\mathrm{gas}})E_{10}e^{E_{10}/kT_{\mathrm{gas}}}+(5/3)\gamma _{21}(T_{\mathrm{gas}})E_{21}e^{E_{21}/kT_{\mathrm{gas}}},$$
(14)
where $`E_{10}/k=128`$ K, $`E_{21}/k=255`$ K and the collisional rate coefficients $`\gamma _{JJ^{}}`$ are given by Flower & Roueff (1999):

$$\gamma _{10}(T_{\mathrm{gas}})=4.4\times 10^{12}+3.6\times 10^{13}(\mathrm{log}T_{\mathrm{gas}})^{0.77}\mathrm{cm}^3\mathrm{s}^1,$$
(15)
and

$$\gamma _{21}(T_{\mathrm{gas}})=4.1\times 10^{12}+2.0\times 10^{13}(\mathrm{log}T_{\mathrm{gas}})^{0.92}\mathrm{cm}^3\mathrm{s}^1.$$
(16)
The comparison between the HD cooling rate calculated by GP and Flower et al. (2000), shown in Fig. 8 is instructive. A simple four-level molecule is able to predict quite accurately the cooling rate in a wide range of temperatures ($`T_{\mathrm{gas}}2000`$ K) and densities ($`n[\mathrm{H}]10^7`$ cm-3), but of course fails badly in the high-temperature, high-density regime. However, for cosmological applications, it is important to keep in mind that the HD cooling function of Flower et al. (2000) has been computed assuming that the temperature of the CBR is much smaller than the gas temperature. This approximation is valid at low redshifts (where the two temperatures differ by more than a factor $`10`$ for $`z20`$), but becomes increasingly inaccurate at higher redshifts. Since most cosmological scenarios of structure formation begin at $`z100`$, when $`T_{\mathrm{rad}}300`$ K, the level population of molecules is strongly affected by stimulated emission and absorption. In such conditions, molecules may become an effective heating source for the gas, because the rate of collisional de-excitation of the rovibrational levels is faster than their radiative decay (Khersonskii 1986; Puy et al. 1993).

As an illustration of this effect, we plot in Fig. 9 the net heat transfer function $`(\mathrm{\Gamma }\mathrm{\Lambda })_{\mathrm{HD}}`$ computed with the GP model for $`n(\mathrm{H})=1`$ cm-3 as function of gas temperature at three selected redshifts. Note that for $`T_{\mathrm{gas}}>T_{\mathrm{rad}}`$ the cooling function is significantly decreased from the value computed with $`T_{\mathrm{rad}}=0`$, because of the radiative depopulation of excited states. The sudden drop signals the condition $`T_{\mathrm{gas}}=T_{\mathrm{rad}}`$. Finally, for $`T_{\mathrm{gas}}<T_{\mathrm{rad}}`$ the function changes sign and becomes a net heating term for the gas. In cosmological simulations the heat transfer should be computed self-consistently at each redshift.

## 4 Application: cloud collapse at $`z100`$

In both cold dark matter and baryonic dark matter cosmological scenarios, the first objects predicted to enter in the nonlinear stage are the smallest ones. In cold dark matter models, it is expected that overdense regions with masses $`10^5`$$`10^7`$ $`M_{}`$ first collapse in the redshift range $`10z100`$ (see e.g. Cen, Ostriker, & Peebles 1993). The crucial question is whether molecular cooling allows the baryonic component to dissipate its kinetic energy and collapse on a timescale shorter than a Hubble time. This important question is fully addressed e.g. in the 1-D hydrodynamical calculations by Haiman, Thoul, & Loeb (1996) and in the 3-D numerical simulations presented by Abel, Anninos, & Norman (1997) and Bromm, Coppi, & Larson (1999). Here, we consider a particular issue related to the choice of the initial conditions for collapse calculations.

Following Tegmark et al. (1997), we consider the growth of a “top-hat” overdensity region, an isolated spherical perturbation in a uniform density Universe (see e.g. Padmanabhan 1993). The radius of this region increases at a slower rate than the scale factor of the Universe (but still obeys the Friedman’s equations) and, after reaching a maximum value (turnaround), recollapses to a point. The formation of a singularity is clearly an artifact of the simplified assumptions of the “top-hat” model. In a realistic situation, gas dynamical process (internal pressure, shocks) will eventually halt the collapse of the baryonic component at some finite value of the density. The resulting quasi-equilibrium configuration is a virialized “halo”, and its subsequent evolution depends on the ability of the gas to cool on a timescale shorter than a dynamical timescale.

We show in Fig. 10 the results of a sample run of the Tegmark et al. (1997) evolutionary model, obtained with our chemical network and molecular cooling prescriptions. For the case shown, the cloud reaches virialization at $`z_{\mathrm{vir}}=110`$, where the temperature is $`T_{\mathrm{vir}}=2000`$ K. After virialization, the gas density is assumed to remain constant and uniform for the rest of the run, and we follow the cooling of the gas due to H2 and HD molecules. This is clearly a poor approximation since the cooling of the gas will of course induce an increase in the the density, and the post-virialization evolution of the cloud must be followed with a full hydrodynamical calculation (Galli et al., in preparation). Nevertheless, it is instructive to consider how rapidly the baryons are able to dissipate the thermal energy of the cloud via molecular cooling. As we see in Fig. 10, recombination reduces the ionization fraction of the cloud to negligible values, weakening the Compton cooling of the gas (since $`T_{\mathrm{rad}}<T_{\mathrm{gas}}`$). The rapid rise of H2 (via the H- channel) around $`z_{\mathrm{vir}}`$ causes the sudden cooling of the gas from $`T_{\mathrm{vir}}`$ down to a few hundred degrees, and the subsequent increase in the HD abundance further reduces the gas temperature down to few tens of degrees, establishing again thermal coupling of gas and cosmic radiation at redshift $`z10`$.

Loosely speaking, a Hubble time corresponds to the redshift dropping by a factor $`2^{2/3}1.6`$. For the particular case considered in this illustrative example, at a resdhift $`z=z_{\mathrm{vir}}/1.670`$ the cloud temperature has dropped to $`T_{\mathrm{gas}}350`$, i.e. of a factor $`6`$ with respect to the virial temperature. Molecular cooling (mostly H2 in this range of redshift) therefore enables the cloud to cool significantly within a Hubble time, and eventually to collapse and form luminous objects. Cooling by HD appears to induce a further substantial reduction of the temperature at later times, but the subsequent evolution of the cloud can only be followed with a realistic calculation. However, it must be kept in mind that the mass scale of the fragments formed by the collapse of the cloud depends sensitively on the gas temperature ($`M_J\rho ^{1/2}T_{\mathrm{gas}}^{3/2}`$). The cooling induced by the formation of HD may in fact reduce the gas temperature below $`100`$ K and therefore enable the formation of primordial low-mass stars or even brown dwarfs, a possibility recently explored by Uheara & Inutsuka (2000) (see also Flower & Pineau des Forêts 2001). These exciting results represent a radical cange of perspective in the field of primordial star formation and deserve further investigation.

## 5 Conclusions

We have examined the most effective chemical reactions leading to the formation/destruction of HD molecules in the early Universe, and presented a list of accurate reaction rate coefficients for primordial chemistry calculations that update those adopted in previous studies. We have analyzed the heating/cooling properties of HD molecules, stressing the relevance of a self-consistent calculation of the energy transfer rate between gas and radiation for cosmological applications. As an illustration, we have presented a simplified model for the evolution of density perturbations in the expanding Universe, following previous investigations by Haiman et al. (1996) and Tegmark et al. (1997), but including our comprehensive and updated chemical network. These results, together with recent findings by Uheara & Inutsuka (2000) and Flower & Pineau des Forêts (2001) underline the substantial contribution of HD to gas cooling during the collapse of primordial clouds. A preliminary conclusion from these studies is that HD is at least as important as H2 in determining the thermal balance of zero-metal clouds.

This work is supported by the Italian Ministry for the University and for Scientific and Technological Research (MURST) through a COFIN-2000 grant. It is a pleasure to thank Guillaume Pineau des Forêts and Malcolm Walmsley for informative dicussions on deuterium chemistry. We also thank A. Dalgarno, S. N. Shore and an anonymous referee for useful comments on an earlier version of this paper.


text:
# Segregation in hard sphere mixtures under gravity. An extension of Edwards approach with two thermodynamical parameters

## Abstract

We study segregation patterns in a hard sphere binary model under gravity subject to sequences of taps. We discuss the appearance of the “Brazil nut” effect (where large particles move up) and the “reverse Brazil nut” effects in the stationary states reached by “tap” dynamics. In particular, we show that the stationary state depends only on two thermodynamical quantities: the gravitational energy of the first and of the second species and not on the sample history. To describe the properties of the system, we generalize Edwards’ approach by introducing a canonical distribution characterized by two configurational temperatures, conjugate to the energies of the two species. This is supported by Monte Carlo calculations showing that the average of several quantities over the tap dynamics and over such distribution coincide. The segregation problem can then be understood as an equilibrium statistical mechanics problem with two control parameters.

Segregation is an ubiquitous and intriguing phenomenon observed in vibrated granular mixtures such as powders or sand: in presence of shaking the system is not randomized, but its components tend to separate . An example is the so called “Brazil nut” effect (BNE) where, under shaking, large particles rise to the top and small particles move to the bottom of the container. Several mechanisms have been proposed to explain these phenomena, as for instance “percolation” , geometrical reorganization , convection , “condensation”  or inertia . Interestingly, by changing grains sizes or mass ratio a crossover towards a “reverse Brazil nut” effect (RBNE) was more recently discovered  where small particles segregates to the top and large particles to the bottom. On the whole the criteria to predict segregation in a mixture, an issue of large practical relevance, are still largely unknown .

In the perspective of the Statistical Mechanics approach originally proposed by Edwards to describe powders , we study here segregation patterns, in absence of convection mechanisms, of a binary hard sphere model under gravity subject to sequences of “taps”. We show that a Statistical Mechanics description of segregation is indeed possible by the introduction of two “configurational temperatures” which characterize the macroscopic status of the system.

In Edwards’ approach granular media at rest (i.e., when grains have no kinetic energy and, thus, the bath temperature is zero) can be described in terms of Statistical Mechanics where a single control parameter, called compactivity, plays the role of the temperature of usual thermal systems . We show that the stationary states of the present mixture are indeed independent on the sample history as in a “thermodynamics” system, but are characterized by two control parameters, such as the gravitational energy of the first and of the second species. We postulate that the probability distribution of the system in the stationary state can be obtained by the principle of maximum entropy under the constraint that the two energies are independently fixed. This leads to a Gibbs canonical distribution characterized by two configurational temperatures, conjugate to the energies of the two species, generalizing in this way Edwards approach  . The postulate is supported by Monte Carlo simulations showing for several quantities that statistical averages over the tap dynamics coincide with those over such an ensemble distribution. The segregation problem can then be studied as standard Statistical Mechanics problems using such a distribution and, for instance, a segregation “phase” diagram can be theoretically derived.

We consider a binary hard sphere system made of two species 1 (small) and 2 (large) with grain diameters $`a_0`$ and $`\sqrt{2}a_0`$, under gravity on a cubic lattice of spacing $`a_0=1`$ (see inset of Fig.1). The hard core potential $`_{hc}`$ is such that two nearest neighbor particles cannot overlap. This implies that only pairs of small particles can be nearest neighbors on the lattice. The system overall Hamiltonian is:

$$=_{hc}+m_1gH_1+m_2gH_2,$$
(1)
where $`H_1=_i^{(1)}z_i`$ and $`H_2=_i^{(2)}z_i`$, the height of site $`i`$ is $`z_i`$ and the two sums are over all particles of species 1 and 2 respectively. We set the units such that the two kinds of grain have masses $`m_1=1`$ and $`m_2=2`$, and the gravity acceleration is $`g=1`$. In the above units, the gravitational energies in a given configuration are thus $`E_1=H_1`$ and $`E_2=2H_2`$.

Grains are confined in a box of linear size $`L`$ between hard walls and periodic boundary conditions in the horizontal directions. The system is subject to a Monte Carlo dynamics made of sequences of “taps” . During each “tap” the system evolves for a time $`\tau _0`$ (the “tap duration”, in units of lattice sweeps) in presence of a finite thermal bath of temperature $`T_\mathrm{\Gamma }`$; afterwards it is suddenly frozen at zero temperature in one of its stable states. These states correspond to a local minimum or saddle of the energy, such that any particle movement does not decrease the energy, and have been called “inherent states”  in analogy with the terminology of glasses . After each tap, when the system is at rest, we record the quantity of interest, which are a function of the number of “taps”, $`t`$.

$`N_1`$ grains of type 1 and $`N_2`$ grains of type 2 are initially prepared in a random loose stable pack . Under “shaking” the system is let to approach a stationary state for each value of the tap parameters $`T_\mathrm{\Gamma }`$ and $`\tau _0`$, as shown in Fig.1 for the large grains height $`h_2(t)=H_2(t)/N_2`$. In Fig.2, we plot as function of $`T_\mathrm{\Gamma }`$ (for several values of $`\tau _0`$) the asymptotic value of the vertical segregation parameter, i.e., the difference of the average heights of the small and large grains at stationarity, $`\mathrm{\Delta }h(T_\mathrm{\Gamma },\tau _0)h_1h_2`$. Here $`h_1`$ and $`h_2`$ are the average of $`H_1/N_1`$ and $`H_2/N_2`$ over the tap dynamics in the stationary state.

Fig.2 shows that at high $`T_\mathrm{\Gamma }`$ the Brazil Nut Effect (BNE) of granular media is observed: “large” grains are found on the top of the system (i.e., $`\mathrm{\Delta }h<0`$). For small values of $`T_\mathrm{\Gamma }`$, below a crossover amplitude, $`T_\mathrm{\Gamma }^{}(\tau _0)`$, a smooth crossover to the Reverse Brazil Nut Effect (RBNE) is instead recorded, with $`\mathrm{\Delta }h>0`$. The RBNE was first discovered in recent Molecular Dynamics simulations . Our aim below is to show that the properties of these segregation patterns can be explained in terms of a statistical mechanics approach to the problem.

The results given in the main panel of Fig. 2 apparently show that $`T_\mathrm{\Gamma }`$ is not a right “thermodynamical” parameter, since sequences of “taps” with different $`\tau _0`$ give different values for the system observables. However, if the stationary states corresponding to different tap dynamics (i.e., different $`T_\mathrm{\Gamma }`$ and $`\tau _0`$), are indeed characterized by a single “thermodynamical” parameter á la Edwards, the curves of Fig.2 should collapse on an universal master function when $`\mathrm{\Delta }h(T_\mathrm{\Gamma },\tau _0)`$ is parametrically plotted as function of an other macroscopic observable such as the average energy, $`e(T_\mathrm{\Gamma },\tau _0)=(E_1+E_2)/N=\frac{N_1}{N}h_1+\frac{2N_2}{N}h_2`$ (where $`N=N_1+N_2`$). This collapse of data, actually found in other situations , is clearly not observed here, as it is apparent in the inset of Fig.2 (see also ).

We show, instead, that two macroscopic quantities, such as $`h_1`$ and $`h_2`$, may be sufficient to characterize uniquely the stationary state of the system. Namely, we show that any macroscopic quantity $`Q`$, averaged over the tap dynamics in the stationary state, is only dependent on $`h_1`$ and $`h_2`$, i.e., $`Q=Q(h_1,h_2)`$. We have checked that this is the case for several independent observables, such as the number of contacts between large particles, $`N_c`$, the density of small and large particles on the bottom layer, $`\rho _1^b`$ and $`\rho _2^b`$, and others. As shown in Fig.3, we find, in particular, with good approximation that: $`N_cN_c(e)=N_c(\frac{N_1}{N}h_1+\frac{2N_2}{N}h_2)`$, $`\rho _2^b\rho _2^b(h_2)`$, $`\rho _1^b\rho _1^b(h_1)`$. Therefore we need both $`h_1`$ and $`h_2`$ (or equivalently $`E_1`$ and $`E_2`$) to characterize unambiguously the state of the system independently on the previous history (i.e., in our case independently on the particular tapping parameters $`T_\mathrm{\Gamma }`$ and $`\tau _0`$).

In this respect the stationary state can be genuinely considered as a “thermodynamical state”, and consequently one can attempt to construct an equilibrium statistical mechanics, as originally suggested by Edwards. This will allow to calculate any quantity in the stationary state using an appropriate ensemble, irrespective of the particular history.

In the present paper we generalize Edwards’ approach by introducing a canonical distribution characterized by two configurational temperatures, conjugate to the energies of the two species, and verify by Monte Carlo calculations that the average of several quantities over the tap dynamics and over such distribution coincide. More precisely we have to find the probability distribution, $`P_r`$, of observing the system in the generic inherent state $`r`$ corresponding to an energy $`E_{1r}`$ for the small particles and $`E_{2r}`$ for the large particles (see ). To answer this, we assume that the microscopic distribution is given by the principle of maximum entropy $`S=_rP_r\mathrm{ln}P_r`$ under the condition that the average energy $`E_1=_rP_rE_{1r}`$ and $`E_2=_rP_rE_{2r}`$ are independently fixed. This can be done by introducing two Lagrange multipliers $`\beta _1`$ and $`\beta _2`$, which are determined by the constraint on $`E_1`$ and $`E_2`$, and can be considered as the “inverse configurational temperature” of species 1 and 2 . This procedure leads to the Gibbs result:

$`P_r={\displaystyle \frac{e^{\beta _1E_{1r}\beta _2E_{2r}}}{Z}}`$ (2)
where $`Z=_r\mathrm{exp}\left[\beta _1E_{1r}\beta _2E_{2r}\right]`$ and, in the thermodynamic limit, the entropy, $`S`$, is given by:

$`S=\mathrm{ln}\mathrm{\Omega }(E_1,E_2),`$ (3)
and $`\beta _1`$ and $`\beta _2`$:

$`\beta _1={\displaystyle \frac{\mathrm{ln}\mathrm{\Omega }(E_1,E_2)}{E_1}},\beta _2={\displaystyle \frac{\mathrm{ln}\mathrm{\Omega }(E_1,E_2)}{E_2}}.`$ (4)
Here $`\mathrm{\Omega }(E_1,E_2)`$ is the number of inherent states corresponding to energy $`E_1`$ and $`E_2`$.

The hypothesis that the ensemble distribution at stationarity is given by eq.(2) can be tested as follows. We have to check that the time average of any quantity, $`Q(h_1,h_2)`$, as recorded during the taps sequences in a stationary state characterized by given values $`h_1`$ and $`h_2`$, must coincide with the ensemble average, $`Q(h_1,h_2)`$, over the distribution eq.(2). To this aim, we have calculated the ensemble averages $`N_c`$, $`\rho _2^b`$, $`\rho _1^b`$, for different values of $`\beta _1`$ and $`\beta _2`$; then we have expressed parametrically $`N_c`$, $`\rho _2^b`$, $`\rho _1^b`$, as function of the average of $`h_1`$ and $`h_2`$, and compared them with the corresponding quantities, $`N_c`$, $`\rho _1^b`$ and $`\rho _2^b`$, averaged over the tap dynamics. The two sets of data are plotted in Fig.3 showing a rather good agreement (notice, there are no adjustable parameters). Technically, in order to calculate the ensemble averages we simulated, with standard Monte Carlo methods, the model with $``$ from eq.(1) where we impose that the only accessible states are the inherent states. This can be done by using a configurational Hamiltonian

$$_{conf}=_{hc}+\beta _1m_1gH_1+\beta _2m_2gH_2\mathrm{ln}\mathrm{\Pi },$$
(5)
where $`\mathrm{\Pi }=1`$ or $`0`$ depending whether the configuration is an inherent state or not. In this way the weight $`e^{_{conf}}`$ gives the ensamble distribution (2) (for further details see ). The discovery of two configurational temperatures may suggest the existence of more than one “dynamical temperature” in the off-equilibrium fluctuation-dissipation relation approach of .

We show now how the ensemble averages with the distribution (2) can be analytically calculated for a simple toy model of a mixture of two particles. The model can be visualized as a cell of the above hard sphere pack, made of one large and two small grains seating on the vertices of a square lattice, with three layers each of 2 sites. Interestingly, such a three levels system is able to reproduce the general properties found before. In particular, the resulting diagram for segregation is plotted in Fig.4 in the plane $`T_1=(m_1ga_0\beta _1)^1`$, $`T_2=(m_2ga_0\beta _2)^1`$. It shows that for a given value of $`T_1`$ by decreasing $`T_2`$ the system undergoes a crossover, at $`T_2^{}(T_1)`$, from a region where type-2 particles are on average above (BNE) to a region where they are below (RBNE) type-1 particles. The structure of the diagram of Fig.4, for small values of $`T_1`$ and $`T_2`$, can be easily understood by looking at eq.(2): when the configurational temperature of small particles, $`T_1`$, is much smaller than the other, $`T_2`$, the potential energy of small grains is much smaller than that associated to the large ones, therefore large grains segregate on top (BNE). The opposite happens when $`T_2`$ is close to zero and much smaller than $`T_1`$. In particular, in this toy model the RBNE region cannot be met if $`T_2`$ is too large.

The crossover observed in vertical segregation from BNE to RBNE corresponds to the drift of the minimum of the free energy, $`F=\mathrm{ln}Z`$, from negative to positive values of $`\mathrm{\Delta }h`$. In this scenario, the pack undergoing a taps dynamics with “amplitudes” $`T_\mathrm{\Gamma }`$ (and given “duration”, $`\tau _0`$) at stationarity, as for the data shown in Fig.2, can be thus thought to be following a given path, $`T_2=T_2(T_1)`$ (where $`T_1=T_1(T_\mathrm{\Gamma },\tau _0)`$), in the $`(T_1,T_2)`$ plane, crossing from the BNE to the RBNE region.

In conclusion, we have shown that in the present hard sphere binary mixture model Edwards’ approach to describe segregation holds if two thermodynamical parameters are introduced. In the case considered here these two thermodynamical parameters play the role of configurational temperatures related to the two types of particles. By applying all the standard techniques of statistical mechanics of the canonical distribution one can thus predict the complex segregation patterns as found above. In particular we have shown on a hard sphere mixture on a lattice that the segregation phenomena in granular media can be studied applying standard statistical mechanics to a new Hamiltonian (5), which contains the original Hamiltonian of the grains plus a term which constrains the system to be in a inherent state.

In general for a more complex system one might expect more constraints to be imposed, leading to more than two thermodynamical parameters. In practice, the criteria to determine a priori the required parameters can be not easily accessible. However more recently we have extended data for very low energy  and found that only one thermodynamical parameter is necessary to describe the stationary state. This seems to be a general feature . If this is the case, a Statistical Mechanics approach with only one thermodynamical variable may be feasible for low energy.

This work was partially supported by the TMR-ERBFMRXCT980183, INFM-PRA(HOP), MURST-PRIN 2000 and MURST-FIRB 2002. The allocation of computer resources from INFM Progetto Calcolo Parallelo is acknowledged.


text:
# Untitled Document

Better be third then second in a search for a majority opinion

Dietrich Stauffer

Institute for Theoretical Physics, Cologne University

D-50923 Köln, Euroland

Abstract: Monte Carlo simulations of a Sznajd model show that if a near-consensus is formed out of four initially equally widespread opinions, the one which at intermediate times is second in the number of adherents usually loses out against the third-placed opinion.

If several opinions compete against each other via mutual persuasion, and finally a consensus or very large majority emerges, then usually (and also in the present work) the final winner is that opinion which at some intermediate stage had already a narrow majority. If at the end a tiny minority of dissenters remains, it seems plausible that they are the remnants of the second-most-important opinion at some intermediate stage. However, we present here simulations in a four-opinion model, where usually the tiny minority of dissenters was on third and not on second place halfway through the process, while the opinion which then was on second place finally died out. This model is the Sznajd model on the dilute square lattice with diffusing agents of four possible opinions.

The Sznajd model  (see  for a recent review) is put onto a square lattice. Each lattice site initially is either empty, with probability 1/2, or has one of four possible opinions 1, 2, 3, 4, with probability 1/8 each. Then at each time step every occupied site tries to move to an empty neighbour. Afterwards randomly selected pairs of nearest neighbours, who share the same opinion, convince all those neighbours of the pair’s opinion, which differ by at most one unit . If this is done for each lattice site, one time step is completed, and we start again with diffusion and convincing. In this way the rigidity of the standard Sznajd lattices is avoided; in principle everybody can exchange opinions with everybody else. The opinions no longer can change if they have settled onto the choices 1 and 3, or 2 and 4, or 1 and 4, or if one opinion covers everybody.

In all ten simulations of $`301\times 301`$ sites, an opinion fixed point was reached after about 4,000 to 100,000 time steps. In one case, only opinion 3 survived; in all other cases, of the two opinions which survived at the end, one had only 12 to 332 adherents compared with the about 45,000 of the winner. But this tiny surviving minority was on third place half way through the process, while the opinion which at half time had much more (4676 to about 22500) adherents finally had none. (In one case the leading and second opinion at half time had about the same number of votes. The opinion ranked fourth at half time always died out.) So to be first or third is good, while the second place is dangerous. (With 10000 samples of $`101\times 101`$ sites, an opinion fixed point was always reached, and the at half time second-ranked opinion finally vanished in about 92 percent of the cases; for $`31\times 31`$ exception were less rare.)

The explanation is based on the discreteness of the four opinions and the fact that opinion 1 is not regarded as similar to opinion 4. The two extreme opinions 1 and 4 thus can convince only one neighbouring opinion each, while the two centrist opinions 2 and 3 have two neighbours each. After some time, most of the opinions will be centrist (2 or 3). These two centrist opinions then fight for a clear majority, one is winning and also will convince the extremist opinion close to it, while the other centrist opinion is losing out completely and thus leaves its neighbouring extremist opinion untouched. This mechanism should also work in other models like Potts spins at low temperatures, as long as the opinions are discrete: Being second means to lose completely; being third allows a small chance of survival.

In the final fixed configuration for the 4 possible states, in 10000 simulations of $`101\times 101`$ sites each, opinions 1 and 3 survived in 41 percent of the cases, opinions 2 and 4 in another 41 percent, while 9 percent each had only opinion 2 or only opinion 3 surviving. For three instead of four possible opinions, nearly always at the end everybody shared the centrist opinion 2; with five possible opinions, usually a small number of opinions 1 and 5 together with a big majority for opinion 3, and without any opinions 2 and 4, survived. For $`31\times 31`$ lattices more exceptions occur.

In summary, for survival of an opinion among four choices it may be better to hide on third place then to be nearly the winner.

 K. Sznajd-Weron and J. Sznajd, Int. J. Mod. Phys. C 11, 1157 (2000).

 D. Stauffer, Journal of Artificial Societies and Social Simulation 5, No.1, paper 4 (2002) (jasss.soc.surrey.ac.uk).

 G. Deffuant, D. Neau, F. Amblard and G. Weisbuch, Adv. Complex Syst. 3, 87 (2000); R. Hegselmann and M. Krause, for Journal of Artificial Societies and Social Simulation 5 (2002).


text:
# IC/2002/8; UPR-980-T hep-ph/0202090 Weak Mixing Angle and the 𝑆𝑈(3)_𝐶×𝑆𝑈(3) Model on 𝑀⁴×𝑆¹/(𝑍₂×𝑍₂')

\[

## Abstract

We show that the desirable weak mixing angle $`\mathrm{sin}^2\theta _W=0.2312`$ at $`m_Z`$ scale can be generated naturally in the $`SU(3)_C\times SU(3)`$ model on $`M^4\times S^1/(Z_2\times Z_2^{})`$ where the gauge symmetry $`SU(3)`$ is broken down to $`SU(2)_L\times U(1)_Y`$ by orbifold projection. For a supersymmetric model with a TeV scale extra dimension, the $`SU(3)`$ unification scale is about hundreds of TeVs at which the gauge couplings for $`SU(3)_C`$ and $`SU(3)`$ can also be equal in the mean time. For the non-supersymmetric model, $`SU(2)_L\times U(1)_Y`$ are unified at order of 10 TeV. These models may serve as good candidates for physics beyond the SM or MSSM. PACS: 11.25.Mj; 11.10.Kk; 04.65.+e; 11.30.Pb

\]

Grand Unification Theory (GUT) has long been an intriguing proposal for physics beyond the Standard Model (SM) or Minimal Supersymmetric Standard Model (MSSM). The impressive successes of supersymmetric (SUSY) GUTs include the successful prediction of the weak mixing angle at the electroweak (EW) scale, the unification of three running gauge couplings of the $`G_{SM}=SU(3)_C\times SU(2)_L\times U(1)_Y`$ group at GUT scale $`M_{GUT}2\times 10^{16}`$ GeV. Besides the proton decay problem, an unsatisfactory feature of this idea is that there exists a SUSY desert for new physics from TeV scale to the GUT scale. Several years ago, gauge coupling unification has been reconsidered in the large extra dimension scenarios. The point of this proposal is that beyond the compactification scale, the KK modes arising from the bulk fields would give radiative corrections to the gauge couplings which would behave as polynomial functions of the energy scale . Gauge coupling unification can be qualitatively achieved around $`10^6`$ GeV if the compactification scale is $`10^5`$ GeV, and the power law unification can be correct at quantitative level if the high energy fundamental theory is known.

On the other hand, the discrete symmetry on the extra space manifold gives another interesting proposal for GUT breaking or gauge symmetry breaking . The doublet-triplet splitting problem and proton decay problem can be solved neatly in this kind of scenario. In addition, the discrete symmetry may not act freely on the extra space manifold. When the discrete symmetry does not act freely on the extra space manifold, there exists a brane at each fixed point, line or hypersurface, where only part of the gauge symmetry and SUSY might be preserved and the SM fermions can be located at . However, the SUSY desert still remains since the compactification scale is close to the 4-dimensional (4D) GUT scale.

About thirty years ago, Weinberg tried to unify the $`SU(2)_L`$ weak and $`U(1)_Y`$ hypercharge interactions into the $`SU(3)`$ gauge interaction where the tree-level weak mixing angle ($`\mathrm{sin}^2\theta _W`$) arising from the breaking of a $`SU(3)`$ group is $`0.25`$ which is within $`10\%`$ of the present experimental value of $`0.2312`$ . However, the quark doublets can not be accommodated in the theory due to their small hypercharge quantum numbers. In a recent paper , Dimopoulos and Kaplan pointed out that the above tree-level prediction of weak mixing angle ($`\mathrm{sin}^2\theta _W=0.25`$) can be obtained even if the weak and hyperecharge interactions were not unified. And the new physics scale is expected to be around $`3.75`$ TeV using the SM gauge couplings and beta functions. The idea is quite interesting but the model they proposed is too complicated. Their model involved the $`SU(3)\times SU(2)\times U(1)`$ gauge group altogether to get the SM $`SU(2)_L\times U(1)_Y`$ gauge group by Higgs mechanism. One problem is that, as discussed in their paper, to achieve a good value of $`\mathrm{sin}^2\theta _W`$ around $`0.25`$, the gauge couplings of the extra $`SU(2)\times U(1)`$ have to be much larger than that of $`SU(3)`$ group, which determines the gauge couplings of the $`SU(2)_L\times U(1)_Y`$ group. Consequently the $`U(1)`$ gauge group may have Landau pole problem not far beyond the $`3.75`$ TeV.

How to combine above ideas to construct a minimal simple model is an interesting question. As we know, the gauge symmetry breaking by the discrete symmetry on extra space manifold will not reduce the rank of the bulk gauge group if the discrete symmetry is abelian. So, only the $`SU(5)`$, $`SU(4)_C\times SU(2)_L`$ and $`SU(3)_C\times SU(3)`$ gauge symmetries can be broken directly down to $`G_{SM}`$ without extra U(1). The $`SU(5)`$ model has been discussed extensively and the $`SU(4)_C\times SU(2)_L`$ have relatively higher unification scale comparing to TeV. In this letter, we would like to discuss the $`SU(3)_C\times SU(3)`$ model on the space-time $`M^4\times S^1/(Z_2\times Z_2^{})`$ in which the correct weak mixing angle $`\mathrm{sin}^2\theta _W=0.2312`$ at the $`m_Z`$ scale can be obtained naturally with a TeV scale extra dimension. By the way, the idea that part of the Standard Model particles live in large (TeV) extra dimensions and the orbifold models based on $`SU(3)\times SU(3)`$ were also discussed in Ref. .

We will concentrate on the SUSY model because the discussion of non-SUSY model is similar. Here is our convention. We consider the 5-dimensional (5D) space-time as the product of the 4D Minkowski space-time $`M^4`$, and the circle $`S^1`$ with radius $`R`$. The corresponding coordinates for the space-time are $`x^\mu `$, ($`\mu =0,1,2,3`$), $`yx^5`$. In addition, the orbifold $`S^1/(Z_2\times Z_2^{})`$ is obtained from $`S^1`$ moduled by two equivalent classes: $`yy`$ and $`y^{}y^{}`$ where $`y^{}=y+\pi R/2`$. $`y=0,\pi R/2`$ are two fixed points. We assume that the bulk gauge symmetry is $`G=SU(3)_C\times SU(3)`$, and two Higgs hypermultiplets $`\mathrm{\Psi }_u`$ and $`\mathrm{\Psi }_d`$ transforming as $`(1,3)`$ and $`(1,\overline{3})`$ are in the bulk.

The $`N=1`$ SUSY theory with gauge group $`G`$ in 5-dimension have 8 real supercharges, corresponding to $`N=2`$ SUSY in 4-dimension. The vector multiplet physically contains a vector boson $`A_M`$ where $`M=0,1,2,3,5`$, two Weyl gauginos $`\lambda _{1,2}`$, and a real scalar $`\sigma `$. In the language of 4D $`N=1`$ SUSY , it contains a vector multiplet $`V(A_\mu ,\lambda _1)`$ and a chiral multiplet $`\mathrm{\Sigma }((\sigma +iA_5)/\sqrt{2},\lambda _2)`$ in the adjoint representation. And the 5-dimensional hypermultiplet physically has two complex scalars $`\varphi `$ and $`\varphi ^c`$, two Weyl fermions $`\psi `$ and $`\psi ^c`$, and can be decomposed into two chiral mupltiplets $`\mathrm{\Phi }(\varphi ,\psi )`$ and $`\mathrm{\Phi }^c(\varphi ^c,\psi ^c)`$, which transform as conjugate representations of each other under $`G`$.

We define the $`P`$ and $`P^{}`$ as corresponding operators for the $`Z_2`$ symmetry transformations $`yy`$ and $`y^{}y^{}`$, respectively. With general action for the 5-dimensional $`N=1`$ SUSY gauge theory and its couplings to the bulk hypermultiplets in , we obtain the symmetry transformations for the gauge fields, and $`\mathrm{\Phi }_u`$, $`\mathrm{\Phi }_u^c`$ under $`P`$

$`V(x^\mu ,y)V(x^\mu ,y)=PV(x^\mu ,y)P^1,`$ (1)
$`\mathrm{\Sigma }(x^\mu ,y)\mathrm{\Sigma }(x^\mu ,y)=P\mathrm{\Sigma }(x^\mu ,y)P^1,`$ (2)
$`\mathrm{\Phi }_u(x^\mu ,y)\mathrm{\Phi }_u(x^\mu ,y)=P\mathrm{\Phi }_u(x^\mu ,y),`$ (3)
$`\mathrm{\Phi }_u^c(x^\mu ,y)\mathrm{\Phi }_u^c(x^\mu ,y)=P^{}\mathrm{\Phi }_u^c(x^\mu ,y),`$ (4)
where $`P^{}`$ is the complex conjugate of $`P`$. And under $`P`$, the transformation properties of $`\mathrm{\Phi }_d`$ and $`\mathrm{\Phi }_d^c`$ are similar to those of $`\mathrm{\Phi }_u`$ and $`\mathrm{\Phi }_u^c`$. Similar results hold for $`P^{}`$. Because of the $`Z_2`$ symmetry, the eigenvalues of $`P`$ and $`P^{}`$ must be $`\pm 1`$. Denoting the general bulk field $`\varphi `$ with parities ($`\pm ,\pm `$) under ($`P`$, $`P^{}`$) as $`\varphi _{\pm \pm }`$, the solutions under Fourier-expansion can be found in . Reducing to 4D fields, $`\varphi _{++}`$ and $`\varphi _{}`$ correspond respectively to the KK modes with masses $`2n/R`$ and $`(2n+2)/R`$, $`\varphi _+`$ and $`\varphi _+`$ are the KK modes with masses $`(2n+1)/R`$ (n is non-negative integer). We emphasize that the zero modes are contained only in $`\varphi _{++}`$ fields. $`\varphi _{}`$ and $`\varphi _+`$ would vanish at $`y=0`$ and $`\varphi _{}`$ and $`\varphi _+`$ would vanish at $`y=\pi R/2`$.

We choose the matrix representaions for $`P`$ and $`P^{}`$ as $`P=\mathrm{diag}(1,1,1)\mathrm{diag}(1,1,1)`$ and $`P^{}=\mathrm{diag}(1,1,1)\mathrm{diag}(1,1,1)`$, which are the group elements of $`SU(3)_C\times SU(3)`$ up to a $`Z_2`$ phase. With general discussions and results for the gauge symmetry breaking by the discrete symmetry on extra space manifold in , we obtain that under $`P`$, the $`SU(3)_C\times SU(3)`$ gauge symmetry is broken down to $`G_{SM}`$ on the 3-brane at $`y=0`$ for all the KK modes and in the bulk for zero modes. And under $`P^{}`$, the 4D $`N=2`$ SUSY is broken down to $`N=1`$ for zero modes. The 3-branes at $`y=0`$ and $`y=\pi R/2`$ preserve only 4D $`N=1`$ SUSY due to the $`Z_2`$ projections. These results can also be obtained by noticing the parities of $`V`$ and $`\mathrm{\Sigma }`$ for the $`SU(3)`$ group, and of $`\mathrm{\Phi }_u`$ and $`\mathrm{\Phi }_u^c`$ (similar for $`\mathrm{\Phi }_d`$ and $`\mathrm{\Phi }_d^c`$) under $`(P,P^{})`$:

$`V:\left(\begin{array}{ccc}(+,+)& (+,+)& (,+)\\ (+,+)& (+,+)& (,+)\\ (,+)& (,+)& (+,+)\end{array}\right),`$ (5)
$`\mathrm{\Sigma }:\left(\begin{array}{ccc}(,)& (,)& (+,)\\ (,)& (,)& (+,)\\ (+,)& (+,)& (,)\end{array}\right),`$ (6)
$`\mathrm{\Phi }_u:\left(\begin{array}{c}(+,+)\\ (+,+)\\ (,+)\end{array}\right),\mathrm{\Phi }_u^c:\left(\begin{array}{c}(,)\\ (,)\\ (+,)\end{array}\right).`$ (7)
For $`SU(3)_C`$ group, $`V`$ and $`\mathrm{\Sigma }`$ have parities $`(+,+)`$ and $`(,)`$, respectively. Then, the MSSM quarks and leptons can be put on the observable 3-brane at $`y=0`$, and the Higgs doublets $`H_u`$ and $`H_d`$ are obtained from the zero modes of $`\mathrm{\Phi }_u`$ and $`\mathrm{\Phi }_d`$. The Yukawa terms are similar to those in the MSSM for on the observable 3-brane, we only have Higgs doublets $`H_u`$ and $`H_d`$ for all the KK modes.

$`SU(2)_L`$ group is generated by $`T^a\sigma ^a/2`$ ($`a=1,2,3`$) of the corresponding $`SU(3)`$ group, where $`\sigma ^a`$ is the Pauli matrix. $`U(1)_Y`$ is generated by $`\sqrt{3}T^8`$, so, the hypercharges of the remaining two Higgs doublets $`H_u`$ and $`H_d`$ are $`1/2`$ and $`1/2`$, respectively. From the $`SU(3)`$ relation we would get $`g_1=g_2/\sqrt{3}`$. On the observable 3-brane at $`y=0`$, we can then put the left handed quark and lepton fields in the fundamental representation of the $`SU(2)_L`$ group while the right handed quarks and leptons are singlets. Knowing that the total $`U(1)_Y`$ charges for each Yukawa term must be zero, together with four anomaly-free conditions: $`(3_C)^2Y`$, $`(2_L)^2Y`$, $`Y^3`$ and $`(grav)^2Y`$, we can determine the correct hypercharges for $`Q`$, $`U`$, $`D`$, $`L`$ and $`E`$ fields. Since the hypercharges of two Higgs doublets are quantized as they are from triplets of $`SU(3)`$ group, an understanding of the charge quantization is achieved as a result of the consistency of our setup.

There exist two energy scales in our models: the compactification scale $`M_C=1/R`$, and the $`SU(3)`$ unification scale which is considered as the cutoff scale $`\mathrm{\Lambda }`$ in the theory. Because there might exist additional $`SU(2)_L`$ and $`U(1)_Y`$ kinetic terms localized on the observable 3-brane which may violate the $`SU(3)`$ gauge coupling relation, the cutoff scale might be much larger than the compactification scale if those localized kinetic terms were not small. We will neglect those extra kinetic terms for two reasons: (1) It is natural to set those localized kinetic terms to zero at tree level in the fundamental theory, so, they can only be generated at loop level as counter terms and are very small if the theory is weakly coupled at cutoff scale; (2) If the gauge interaction is strong coupled at $`\mathrm{\Lambda }`$ and $`M_C0.01\mathrm{\Lambda }`$, the effects of those kinetic terms are very small .

First, let us study the scenario with $`M_C\mathrm{\Lambda }`$. Assuming the sparticle mass scale, $`M_S`$, in the range $`2001000`$ GeV, using the MSSM beta function $`(11,1,3)`$ above this scale and the SM beta function $`(41/6,19/6,7)`$ down to the $`m_Z`$ scale, in the limit that $`\mathrm{sin}^2\theta _W=0.25`$ is realised at $`M_C`$ or $`\mathrm{\Lambda }`$, we find out that $`M_C7714`$ TeV.

Second, we discuss the scenario with $`M_C0.01\mathrm{\Lambda }`$. For $`M_C<\mu <\mathrm{\Lambda }`$, we should include the massive KK modes in counting the radiative corrections to gauge couplings. There are two kinds of beta functions, $`b_o`$ and $`b_e`$, from the KK modes of bulk fields with masses $`2n/R`$ and $`(2n+1)/R`$, respectively. The radiative corrected gauge couplings above the $`M_C`$ scale is

$`\alpha _i^1(\mu )`$ $`=`$ $`\alpha _i^1(M_C)+{\displaystyle \frac{b_i}{2\pi }}\mathrm{ln}{\displaystyle \frac{\mu }{M_C}}`$ (9)
$`+{\displaystyle \underset{n}{}}{\displaystyle \frac{1}{2\pi }}(b_{io}\mathrm{ln}{\displaystyle \frac{\mu }{M_{2n1}}}+b_{ie}\mathrm{ln}{\displaystyle \frac{\mu }{M_{2n}}}),`$
where $`\alpha _i=g_i^2/4\pi `$, $`b_e=(2,2,6)`$ and $`b_o=(14,2,0)`$. One obtains that the compactification scale is several TeVs in this kind of scenario.

Third, because we include the $`SU(3)_C`$ gauge symmetry in the bulk, it is interesting to study whether there exists the scenario in which the $`SU(3)`$ unification scale or $`\mathrm{\Lambda }`$ is the scale at which the gauge couplings for $`SU(3)_C`$ and $`SU(3)`$ are also equal in the mean time. And if the 5D $`SU(3)`$ theory is strong coupled at cutoff scale, it is reasonable to assume that the gauge coupling for $`SU(3)`$ is similar to that of $`SU(3)_C`$ at $`\mathrm{\Lambda }`$. We would like to emphasize that because we do not present an explicit Grand Unified Theory which unifies the $`SU(3)`$ and $`SU(3)_C`$ gauge interactions in this letter, we do not mean that there exists the gauge unification for $`SU(3)`$ and $`SU(3)_C`$ even if the gauge couplings for $`SU(3)_C`$ and $`SU(3)`$ are equal or similar at the cut-off scale. However, the quantitative study can show whether a further complete gauge coupling unification might be possible or not from the point of view of gauge coupling runnings.

To have an estimate, we extract the leading $`SU(3)`$ symmetric contributions given by the KK modes, that is $`_k\frac{b_{io}+b_{ie}}{4\pi }\mathrm{ln}\frac{\mu }{M_k}`$. Similar to the case in , for $`\mu M_C`$, we may resum this leading contribution to be a power law running, that is $`\frac{b_{io}+b_{ie}}{4\pi }(\frac{\mu }{M_C}1)`$. Using the beta function listed above, we get a rough estimate of the cutoff scale $`\mathrm{\Lambda }`$:

$`{\displaystyle \frac{1}{2\pi }}({\displaystyle \frac{\mathrm{\Lambda }}{M_C}}1)\alpha _2^1\alpha _3^1,`$ (10)
where $`\alpha _{2,3}`$ are the running gauge couplings for $`SU(2)_L`$ and $`SU(3)_C`$ at scale $`M_C`$. The difference in the power law running is given by the KK tower from the Higgs fields which we choose living in the bulk. $`\frac{\mathrm{\Lambda }}{M_C}`$ is about $`100130`$ for $`M_S2001000`$ GeV and $`M_C<10`$ TeV. At first sight this large gap between the cutoff scale and $`M_C`$ requires less than TeV $`M_C`$ to meet the purpose of unifying simultaneously $`SU(2)_L`$ and $`U(1)_Y`$, because the $`SU(3)`$ asymmetric relative running given by the zero mode is $`\frac{b_1/3b_2}{2\pi }ln\frac{\mathrm{\Lambda }}{M_C}`$ which is about $`2.03`$ for $`\mathrm{\Lambda }/M_C=120`$ ($`\alpha _1^1/3\alpha _2^13.21`$ at $`m_Z`$ scale). However, there are also $`SU(3)`$ asymmetric radiative corrections given by the massive KK modes. It is $`_k\frac{b_ob_e}{4\pi }ln\frac{M_{2k}}{M_{2k1}}=\frac{b_ob_e}{4\pi }ln\frac{2}{1}\frac{4}{3}\frac{6}{5}\mathrm{}`$. This relative running between $`3\alpha _1`$ and $`\alpha _2`$ is completely given by the $`U(1)_Y`$ part and is about $`1.11`$ for $`\mathrm{\Lambda }/M_C=120`$. Detailed calculation using logarithmic running including step by step massive KK states shows that in the whole range $`M_S2001000`$ GeV, $`3\alpha _1`$, $`\alpha _2`$ and $`\alpha _3`$ can meet simultaneously for $`M_C9.31.5`$ TeV with $`\alpha ^16066`$ and $`\mathrm{\Lambda }/M_C`$ around 110. And the requirement that $`3\alpha _1`$, $`\alpha _2`$ and $`\alpha _3`$ meet simultaneouly further constrains the model, i. e., $`M_C`$ and $`\mathrm{\Lambda }`$ will be correlated if the SUSY threshold is known. In Fig. 1, we give a plot for $`M_S=400`$ GeV showing the $`\alpha ^1`$ versus the energy scale. We may see clearly the running behaviours of the three gauge couplings. $`3\alpha _1`$, $`\alpha _2`$ and $`\alpha _3`$ meet with a value of about $`1/63`$ at $`\mathrm{\Lambda }494`$ TeV. In this case $`M_C4.5`$ TeV and $`\mathrm{\Lambda }/M_C110`$. We may see that for $`\mu 10`$ TeV, the difference between $`3\alpha _1`$ and $`\alpha _2`$ is quite small (the discrepancy is less than $`0.2\%`$ for $`\mu >100`$ TeV).

The SUSY breaking can be done by the Scherk-Schwarz mechanism or by gauge mediated SUSY breaking. For the first approach, one can solve the $`\mu `$ problem, and get the gaugino masses, the $`\mu B`$ terms and Higgs masses . The radiative electroweak (EW) symmetry breaking may also be realized. For the second approach, the 3-brane at the fixed point $`y=\pi R/2`$ can be considered as hidden sector for SUSY breaking. The SUSY breaking effects may be communicated to the observable sector via gauge-mediated mechanism. This may be realized if there is a SUSY breaking gauge singlet field $`S`$ living in the bulk. Interaction like $`Tr(SW^\alpha W_\alpha )`$ gives gaugino masses via $`F`$-term. $`\mu `$ term, $`\mu B`$ and Higgs mass terms can also be generated via the interaction of $`\mathrm{\Phi }`$ and $`S`$. Squarks and sleptons can then obtain their masses through the radiative corrections. EW symmetry breaking may be achieved in the usual way via loop corrections at the EW scale because of the large top quark Yukawa coupling.

Similarly, we can construct a non-SUSY model on $`M^4\times S^1/(Z_2\times Z_2^{})`$. The triplet scalar Higgs field can live in the bulk, too. The SM fermions are confined on the observable 3-brane. Because the relative running between $`3\alpha _1`$ and $`\alpha _2`$ are much faster than that in the SUSY case for $`\mu <M_C`$, and the relative power law running between $`3\alpha _1`$ and $`\alpha _3`$ or between $`\alpha _2`$ and $`\alpha _3`$ is $`1/12`$ times slower for $`\mu >M_C`$, it is impossible that the three gauge couplings meet simultaneously at an energy scale if $`M_C>1`$ TeV involved. So, it may be better to put $`SU(3)_C`$ on the boundary. Since $`\mathrm{sin}^2\theta _W=0.25`$ is expected to be reached using the SM beta functions at about $`3.75`$ TeV, we should have $`M_C<3.5`$ TeV for this case. The cutoff scale, which is understood as the energy scale where $`3\alpha _1`$ and $`\alpha _2`$ unified, can be tens of TeV by setting $`M_C>1`$ TeV. The problem of this case is that the masses of Higgs zero modes might be a problem. Studying on EW symmetry breaking, if it’s accounted for solely by the radiative corrections (for which $`1\%`$ fine tuning is clearly needed), would give us the strong constraints on this setup.

Furthermore, we can discuss the $`SU(3)_C\times SU(3)`$ model on $`M^4\times S^1/Z_2`$ with gauge-Higgs unification. However, the weak mixing angle $`\mathrm{sin}^2\theta _W`$ is expected about $`0.75`$ at the cutoff scale.

The constraints on $`M_C`$ from the current high energy experiment data are model dependent, and one need to study it in detail. It is possible that the present knowledges might give the strong constraints on the non-SUSY model for $`M_C<`$ 3.5 TeV. Moreover, the extra $`SU(2)_L`$ doublets in the off-diagonal $`\mathrm{\Sigma }`$ ($`A_5`$ for non-SUSY case), and the third components of the $`\mathrm{\Phi }_u^c`$ and $`\mathrm{\Phi }_d^c`$, may be produced at LHC via s-channel $`W^\pm `$, $`Z^0`$ and $`\gamma `$ exchange for $`M_C`$ is several TeVs, and the decay modes are model dependent. For the $`SU(2)_L`$ doublets in the off-diagonal $`V`$ ($`A_\mu `$ for non-SUSY case) and the third components of the $`\mathrm{\Phi }_u`$ and $`\mathrm{\Phi }_d`$, there might exist the derivative interactions with the leptons, so, they might also be produced at LHC and decay into leptons. Of course, the decay channels of the KK modes of bulk fields are dependent on the loop corrections to the masses of KK modes, which should be considered in the collider test of the models.

Conclusion. We have presented the $`SU(3)_C\times SU(3)`$ model on the space-time $`M^4\times S^1/(Z_2\times Z_2^{})`$ with natural explanation for the weak mixing angle which is a consequence of the $`SU(3)`$ symmetry breaking in our setup. For the SUSY model with a TeV scale extra dimension, the $`SU(3)`$ unification scale is about hundreds of TeVs at which the gauge couplings for $`SU(3)_C`$ and $`SU(3)`$ can be equal in the mean time. For the non-SUSY model, the $`SU(3)`$ unification of $`SU(2)_L\times U(1)_Y`$ may be realized at tens of TeVs at which the gauge couplings for $`SU(3)_C`$ and $`SU(3)`$ can not be equal simultaneously. Therefore, the supersymmetric desert does not appear in our models, which is an unsatisfactory feature in the usual GUTs. And our study shows that from the point of view of the gauge coupling runnings, a complete gauge coupling unification might occur for the SUSY case. We stress that both models presented in this letter, supersymmetric or not, are of strong predictive power because of their clean setup. Studying the non-supersymmetric model, one may soon reveal whether it’s realistic or not. For the supersymmetric model, we can break the SUSY by Scherk-Schwarz mechanism or gauge mediated SUSY breaking, and might test it at LHC. We also comment on the charge quantization, radiative EW symmetry breaking, and the models with gauge-Higgs unification. We conclude that this kind of interesting models may be good candidates for the minimal extensions of the SM or MSSM, and their phenomenology deserves further study.

After sending our letter, we noticed the papers  which discuss similar issue.

Acknowledgments. We would like to thank W. Y. Keung for reading the manuscript and useful discussions. This work was supported in part by the U.S. Department of Energy under Grant No. DOE-EY-76-02-3071.