text:
# Dynamic Jahn-Teller Mechanism of Superconductivity in MgB2

## Abstract

We propose a novel mechanism of superconductivity in MgB2 based on the dynamic electronic structure of the $`p\sigma `$-orbitals coupled with $`e_{2g}`$ phonons. A nonconventional superconducting state is found to arise from electron-phonon interactions in the presence of additional pairing channels made available by the dynamic Jahn-Teller effects. A partially broken pseudo-spin symmetry in this Jahn-Teller system, together with two-phonon exchange pairing, naturally gives rise to two distinct gaps both of which are basically isotropic in the ($`k_x`$, $`k_y`$) space. Important experimental observations including high $`T_c`$ and the anomalous specific heat are explained using this theory.

Following the discovery of high $`T_c`$ superconductivity near 40 K in MgB2 akimi01 , a lot of efforts have been devoted to revealing the underlying mechanism of this novel superconductor buzea01 . Now it is widely agreed upon that a relatively strong electron-phonon coupling is responsible for the superconductivity of this material and there exit, at least phenomenologically, two superconducting gaps shin01 ; chen ; szabo ; bouquet01a ; junod01 ; yang01 ; bouquet01b . It has been proposed kortus01b ; choi01a  that the complex Fermi surfaces having both 2 dimensional (2D) cylindrical sheets and 3 dimensional (3D) tubular networks produce two gaps. However, it is still an open question whether one gap is mainly associated with the 2D cylindrical Fermi surfaces and the other is associated with the 3D Fermi surfaces.

In this letter, we propose a mechanism of superconductivity in MgB2 based on a dynamic Jahn-Teller effect arising from the interplay between the doubly degenerate $`p\sigma `$ electronic states and the $`e_{2g}`$ phonon modes. The hopping motion of holes of the $`p\sigma `$ character on the Boron layers is constrained by the ‘pseudo-spin’, which will be defined later, and partial breakdown of the pseudo-spin symmetry is also important in the pairing mechanism. Unlike existing theories, we will explain both high $`T_c`$ and two distinct gaps in MgB2 without invoking the additional 3D Fermi surface of $`p\pi `$ character.

According to the band structure calculations kortus01a ; pickett01a , the $`p\sigma `$ states on the Boron layer form doubly degenerate $`E_g`$ bands at $`\mathrm{\Gamma }`$ (to be precise, along the $`\mathrm{\Gamma }`$A line on which the component of $`𝐤`$ parallel with the Boron plane, $`𝐤_{||}`$, is zero), where we have small hole pockets presumably responsible for the superconductivity in MgB2. (Splitting of the band degeneracy away from $`\mathrm{\Gamma }`$ will be taken into account later.) These hole states are known to be strongly coupled with $`e_{2g}`$ phonon modes with $`\mathrm{}\omega 70meV`$ cava01b ; kong01 . The $`E_g`$ hole doublet coupled with $`e_{2g}`$ phonon doublet states on the Boron layer can be represented by the so-called $`Ee`$ Jahn-Teller model. A typical adiabatic potential energy surface of the $`Ee`$ Jahn-Teller model is illustrated in Fig. 1, where $`q_\theta `$ and $`q_\epsilon `$ correspond to the two phonon coordinates.

It is well known that one must take account of the dynamical phases of the holes and phonons in order to describe the $`Ee`$ system properly. When the ionic configuration traverses around the potential well in Fig. 1 once, the corresponding orbital wavefunction of the hole (for simplicity of discussion, we leave aside the real spin wavefunction which simply follows the elementary BCS theory) acquires a phase of $`1`$, as if the orbital part had spin 1/2. We call this a pseudo-spin of the hole and denote its projected value by $`\alpha =\pm \frac{1}{2}`$. This subject was extensively discussed in the literature bersuker  in the context of the Jahn-Teller effect or Berry’s phase. We also note that the dynamic Jahn-Teller effect in doped C60 superconductors was studied before, although the methods and results there were different from our work auer .

We introduce a single site Jahn-Teller Hamiltonian per unit cell using mass-weighted coordinate for the symmetrized displacement $`Q_\gamma `$($`\gamma =\theta ,\epsilon `$), i.e. $`q_\gamma =\sqrt{M_\gamma }Q_\gamma `$,

$$h_o=\frac{1}{2}(p_\theta ^2+p_\epsilon ^2)+\frac{1}{2}\omega _0^2(q_\theta ^2+q_\epsilon ^2)+A(q_\theta \sigma _x+q_\epsilon \sigma _y),$$
(1)
where $`A`$ is the hole-phonon coupling strength and $`\sigma _x`$, $`\sigma _y`$ are Pauli spin matrices for the pseudo-spin representation of $`E_g`$ hole doublet states. In an analogy to the 2D harmonic oscillator, the vibrational mode in the $`(q_\theta ,q_\epsilon )`$-space belongs to the eigenmode of the “phonon” angular momentum operator $`L_z`$ (=$`q_\theta p_\epsilon q_\epsilon p_\theta `$) with the eigenvalue of $`l`$, and the principal quantum number $`\nu `$ becomes $`\nu =l+2k`$ ($`k`$ is a non-negative integer and $`\nu l\nu `$). In the present theory of superconductivity, only $`\nu =`$ 0, 1, 2 and $`l=0,\pm 1`$ are relevant. Defining the total pseudo-spin of the system by $`J_z=L_z+\frac{1}{2}\sigma _z`$, it is straightforward to prove $`[h_0,J_z]=0`$, that is, the total pseudo-spin of the system is conserved, while individual $`L_z`$ and $`\sigma _z`$ are not. (The phonon angular momentum component, $`L_z`$, can be regarded as a part of the total pseudo-spin in this context.) These pseudo-spin states are derived from the $`Ee`$ Jahn-Teller coupling between $`p\sigma `$ holes and $`e_{2g}`$ phonons, and reflect the non-adiabatic nature of the coupled states. As a consequence, relevant electron and phonon states are expected to be strongly renormalized and contribute to anomalous behaviors (e.g., broad line-widths) observed in recent photoemission spectroscopy and Raman experiments shin01 ; takahashi ; szabo .

We generalize the theory to the Bloch state in the $`𝐤`$-space of a crystal by linear combinations of the local holes of the same pseudo-spin, but mathematical details will be reserved for a longer paper YWS . We first consider the one-phonon ($`\nu =1`$) exchange for the Cooper pair formation. Each phonon has $`l=\pm 1`$ in this case. In the one-phonon absorption or emission, the pseudo-spin of the hole should flip to conserve the total pseudo-spin. The only allowed one-phonon exchange mechanism is the opposite pseudo-spin Cooper pair as illustrated in Fig. 2(a). Pairing of the same pseudo-spin holes violates the conservation law in either one of the two (emission and absorption) processes. Such a selection rule has already been known in the self-energy correction to the Jahn-Teller system takada .

We next consider the two-phonon ($`\nu =2`$) case. In MgB2, the two-phonon exchange process has been thought to be important in producing high $`T_c`$ and included in some previous works kortus01b . When we properly consider the pseudo-spin symmetry of the Jahn-Teller system such as MgB2, incorporating the two-phonon exchange turns out to be essential in order to obtain two distinct gaps. Unlike the one-phonon exchange process, both opposite- and same- pseudo-spin pairs are allowed here as illustrated in Fig. 2(b) and (c). The total $`l`$ (sum of two individual $`l`$’s) of intermediate phonons is restricted to zero (i.e., the hole pseudo-spin change of $`\pm 2`$ is not possible) to conserve the total pseudo-spin. A new model Hamiltonian with different pseudo-spin channels explicitly labelled is, in the $`𝐤`$-space representation,

$$=_0+_1+_2+_3.$$
(2)
$`_0`$ $`=`$ $`{\displaystyle \underset{\alpha 𝐤\sigma }{}}\epsilon _𝐤c_{\alpha 𝐤\sigma }^{}c_{\alpha 𝐤\sigma },`$ (3a)
$`_1`$ $`=`$ $`{\displaystyle \underset{\alpha \mathrm{𝐤𝐤}^{}}{}}V_{\mathrm{𝐤𝐤}^{}}^0c_{\alpha 𝐤^{}}^{}c_{\alpha 𝐤^{}}^{}c_{\alpha 𝐤}c_{\alpha 𝐤},`$ (3b)
$`_2`$ $`=`$ $`{\displaystyle \underset{\alpha \mathrm{𝐤𝐤}^{}}{}}V_{\mathrm{𝐤𝐤}^{}}^1c_{\alpha 𝐤^{}}^{}c_{\alpha 𝐤^{}}^{}c_{\alpha 𝐤}c_{\alpha 𝐤},`$ (3c)
$`_3`$ $`=`$ $`{\displaystyle \underset{\alpha \mathrm{𝐤𝐤}^{}}{}}V_{\mathrm{𝐤𝐤}^{}}^2c_{\alpha 𝐤^{}}^{}c_{\alpha 𝐤^{}}^{}c_{\alpha 𝐤}c_{\alpha 𝐤}.`$ (3d)
$`_1`$, $`_2`$, and $`_3`$ correspond to the diagrams in Fig. 2(a), (b), and (c), respectively. The explicit form of the derived interactions, $`V_{\mathrm{𝐤𝐤}^{}}^i`$’s ($`i=`$0, 1, 2), can be calculated and will be given elsewhere YWS . To discuss the nature of superconducting states, it suffices here to treat them as given parameters. The additional pairing channels due to the presence of pseudo-spin symmetry can give rise to a non-trivial superconducting ground state. Although the pairing is still based on electron-phonon interactions, the pairing order parameter can now attain more than one symmetry representation leading to the ground state with multi-dimensional order parameters while retaining the conventional $`s`$-wave type superconductivity.

In the pairing Hamiltonian described above, however, important features of a real crystal are missing. First, the pseudo-spin symmetry is exact only at $`𝐤_{||}=0`$ and there occurs small but finite deviation from the exact symmetry for $`𝐤_{||}0`$ where Fermi surfaces actually exist. The band splitting at $`𝐤_{||}0`$ is a clear manifestation of the breakdown of the rigorous Jahn-Teller theory for Bloch wavefunctions of $`𝐤_{||}0`$ in a crystal as opposed to the case of a highly symmetric molecule. Second, higher-order Jahn-Teller coupling terms (not considered so far) do not commute with $`J_z`$ and break the pseudo-spin conservation pseudo . A simple way to incorporate the symmetry breaking effect is to introduce a small symmetry breaking parameter $`\varphi _𝐤`$ ($`\varphi _{𝐤_{||}=0}=0`$ and $`\varphi _{𝐤_{||}0}0`$) in the hopping Hamiltonian (Eq. (3a)) as follows,

$$\stackrel{~}{_0}=\underset{\alpha 𝐤\sigma }{}\epsilon _𝐤c_{\alpha 𝐤\sigma }^{}c_{\alpha 𝐤\sigma }\underset{\alpha 𝐤\sigma }{}\varphi _𝐤\epsilon _𝐤c_{\alpha 𝐤\sigma }^{}c_{\alpha 𝐤\sigma }.$$
(4)
This replacement naturally introduces the symmetry-breaking hopping process which gives the band splitting away from $`𝐤_{||}=0`$. (If the pseudo-spin $`\frac{1}{2}`$ were continuous symmetry ($`SU(2)`$) for all $`𝐤`$’s, there would be a single finite gap plus a zero-gap Goldstone mode, unlike the experimental observation in MgB2.) Now we choose the following two order parameters, one for the opposite pseudo-spin pair and the other for the same pseudo-spin pair,

$`A_𝐤=c_{\frac{1}{2}𝐤}c_{\frac{1}{2}𝐤}=c_{\frac{1}{2}𝐤}c_{\frac{1}{2}𝐤},`$ (5a)
$`B_𝐤=c_{\frac{1}{2}𝐤}c_{\frac{1}{2}𝐤}=c_{\frac{1}{2}𝐤}c_{\frac{1}{2}𝐤}.`$ (5b)
It is to be noted that the pseudo-spin pairing assumed above is a triplet state (even for exchange of two particles) to preserve the singlet (odd) state of the real spin pair and the s-wave (even) state of the orbital wavefunction. Following the standard mean field approximation, we can derive the effective Hamiltonian $`^{}`$ ($`=\epsilon _F𝒩`$),

$`^{}`$ $``$ $`{\displaystyle \underset{\alpha 𝐤}{}}\xi _𝐤(c_{\alpha 𝐤}^{}c_{\alpha 𝐤}c_{\alpha 𝐤}c_{\alpha 𝐤}^{})`$ (6)
$`{\displaystyle \underset{\alpha 𝐤}{}}\delta _𝐤(c_{\alpha 𝐤}^{}c_{\alpha 𝐤}c_{\alpha 𝐤}c_{\alpha 𝐤}^{})`$
$`{\displaystyle \underset{\alpha 𝐤}{}}\mathrm{\Delta }_{1,𝐤}(c_{\alpha 𝐤}^{}c_{\alpha 𝐤}^{}+c_{\alpha 𝐤}c_{\alpha 𝐤})`$
$`{\displaystyle \underset{\alpha 𝐤}{}}\mathrm{\Delta }_{2,𝐤}(c_{\alpha 𝐤}^{}c_{\alpha 𝐤}^{}+c_{\alpha 𝐤}c_{\alpha 𝐤})`$
$`+\text{(constant terms)},`$
where $`\xi _𝐤=\epsilon _𝐤\epsilon _F`$, $`\delta _𝐤=\varphi _𝐤\epsilon _𝐤`$, and the gap-like paremeters are expressed as

$$\mathrm{\Delta }_{1,𝐤}=\underset{𝐤^{}}{}(V_{\mathrm{𝐤𝐤}^{}}^0+V_{\mathrm{𝐤𝐤}^{}}^1)A_𝐤^{},\mathrm{\Delta }_{2,𝐤}=\underset{𝐤^{}}{}V_{\mathrm{𝐤𝐤}^{}}^2B_𝐤^{}.$$
(7)
Using the basis $`(c_{\frac{1}{2}𝐤}^{},c_{\frac{1}{2}𝐤}^{},c_{\frac{1}{2}𝐤},c_{\frac{1}{2}𝐤})`$, it is convenient to cast the Hamiltonian into a matrix form (by omitting of the subscript $`𝐤`$ for notational simplicity),

$$^{}=\underset{𝐤}{}\left(\begin{array}{cccc}\xi & \delta & \mathrm{\Delta }_2& \mathrm{\Delta }_1\\ \delta & \xi & \mathrm{\Delta }_1& \mathrm{\Delta }_2\\ \mathrm{\Delta }_2& \mathrm{\Delta }_1& \xi & \delta \\ \mathrm{\Delta }_1& \mathrm{\Delta }_2& \delta & \xi \end{array}\right).$$
(8)
This Hamiltonian is easily diagonalized, for each $`𝐤`$, yielding two positive eigenvalues,

$$E_\pm =\sqrt{(\xi \delta )^2+(\mathrm{\Delta }_2\pm \mathrm{\Delta }_1)^2}.$$
(9)
The larger gap $`\mathrm{\Delta }_+(=\mathrm{\Delta }_2+\mathrm{\Delta }_1)`$ is associated with the band of $`E_{\text{normal}}=\xi \delta `$ and the smaller gap $`\mathrm{\Delta }_{}(=|\mathrm{\Delta }_2\mathrm{\Delta }_1|)`$ is associated with the band of $`E_{\text{normal}}=\xi +\delta `$. Because of this band splitting, two gaps are associated with two different normal-state (above $`T_c`$) densities of states (DOSs). If we define $`N`$ to be the normal-state DOS of unsplitted $`p\sigma `$ band at the Fermi level before the symmetry-breaking $`\varphi `$ is introduced, normal-state DOSs of two bands are $`N_\pm =\frac{N}{1\varphi }`$ with the assumption of $`\varphi _𝐤`$’s being independent of $`𝐤`$ ($`\varphi _𝐤=\varphi `$) in the vicinity of the Fermi level. Experimentally, it is observed that the state involved in the larger gap has the larger normal-state DOS than that of the smaller one bouquet01b . In this situation, we naturally assign a positive value for $`\varphi `$ ($`N_+=\frac{N}{1\varphi }>N_{}=\frac{N}{1+\varphi }`$).

We first follow, for heuristic purposes, the approximation procedure of the weak-coupling BCS theory. We assume some effective cutoff energy $`\omega `$ and define $`\lambda _1=N(V^0+V^1)`$ and $`\lambda _2=NV^2`$. Only $`\lambda _1`$, $`\lambda _2`$ and $`\omega `$ are independent variables for the gap equations. The coupled gap equations for finite temperature read as

$`\mathrm{\Delta }_+`$ $`=`$ $`{\displaystyle \frac{\lambda _2+\lambda _1}{1\varphi }}{\displaystyle _0^\omega }𝑑\xi {\displaystyle \frac{\mathrm{\Delta }_+}{2E_+}}\mathrm{tanh}{\displaystyle \frac{\beta E_+}{2}}`$
$`+{\displaystyle \frac{\lambda _2\lambda _1}{1+\varphi }}{\displaystyle _0^\omega }𝑑\xi {\displaystyle \frac{\mathrm{\Delta }_{}}{2E_{}}}\mathrm{tanh}{\displaystyle \frac{\beta E_{}}{2}},`$
$`\mathrm{\Delta }_{}`$ $`=`$ $`{\displaystyle \frac{\lambda _2+\lambda _1}{1+\varphi }}{\displaystyle _0^\omega }𝑑\xi {\displaystyle \frac{\mathrm{\Delta }_{}}{2E_{}}}\mathrm{tanh}{\displaystyle \frac{\beta E_{}}{2}}`$ (10)
$`+{\displaystyle \frac{\lambda _2\lambda _1}{1\varphi }}{\displaystyle _0^\omega }𝑑\xi {\displaystyle \frac{\mathrm{\Delta }_+}{2E_+}}\mathrm{tanh}{\displaystyle \frac{\beta E_+}{2}}.`$
By fitting to $`T_c`$ and the specific heat behavior, we are going to obtain below $`\lambda _1=0.235`$, $`\lambda _2=0.25`$ and $`\varphi =0.2`$ for the assumed value of $`\mathrm{}\omega =74.5meV`$. The normalized solutions of the above gap equations (with $`\frac{\lambda _1}{\lambda _2}0.94`$ and $`\varphi =0.2`$) are shown in Fig. 3(a). Our theory gives two distinct gaps of the identical $`T_c`$ originated from two split-off $`p\sigma `$ bands alone, in contrast to other theories which invoke the $`p\pi `$ band to obtain multi-gaps. Furthermore, we predict that there exist two gaps for the same $`𝐤`$ vector and both gaps are basically isotropic in the Boron plane (except for a small anisotropic influence from the $`p\pi `$ bands located at different regions of the $`𝐤`$ space) unlike other theories proposing generically anisotropic gaps. Without pseudo-spins, there would be only one s-wave gap solution giving the lowest free energy. It should be possible, with angle-resolved experimental probes, to test the symmetry or angular dependence of the gap order parameters. Also, it is important in the present theory to include the symmetry-breaking parameter $`\varphi _𝐤`$; otherwise, two gaps should be exactly degenerate and only one gap would manifest itself. With the same degree of approximation as in the BCS formalism, we can derive the expression for $`T_c`$ (weak coupling limit),

$$k_BT_c1.13\mathrm{}\omega \mathrm{exp}\left[\frac{(\lambda _1+\lambda _2)\sqrt{(\lambda _1\lambda _2)^2+4\varphi ^2\lambda _1\lambda _2}}{2\lambda _1\lambda _2}\right].$$
(11)
We then calculate the specific heat with this gap structure. The resulting specific heat is plotted in a usual manner in Fig. 3(b). To obtain both $`T_c=39K`$ and the specific heat data in agreement with experiment, we require $`\varphi =0.2`$, $`\lambda _1=0.235`$, $`\lambda _2=0.25`$. The symmetry breaking parameter $`\varphi =0.2`$ gives about the same band splitting as ab initio calculations for the two $`p\sigma `$ bands near the Fermi level. Needless to say, since $`T_c`$ is high and the material parameters are not in the weak-coupling limit, Eq. (11) is not to be applied directly to MgB2. We have to use a more realistic (Eliashberg) equation which gives, for example, $`T_c\frac{\omega _{\text{ln}}}{1.20}\mathrm{exp}\left[\frac{1.04(1+\lambda )}{\lambda (1+0.62\lambda )\mu ^{}}\right]`$ allen . It is beyond the scope of the present paper to calculate $`\lambda _1`$ and $`\lambda _2`$ from the first principles. The effective coupling constant $`\lambda `$ in the exponent in Eq. (11), namely, $`\frac{2\lambda _1\lambda _2}{(\lambda _1+\lambda _2)\sqrt{(\lambda _1\lambda _2)^2+4\varphi ^2\lambda _1\lambda _2}}`$ should be about 0.95, as is well known pickett01a ; kong01 , in order to give $`T_c=39`$ K with $`\mu ^{}=0.13`$. Then we obtain $`\lambda _1=0.71`$, $`\lambda _2=0.76`$ by simple scaling of the previous values. These values are not to be taken literally because the two-phonon exchange can modify the empirical expression for $`T_c`$ of the Eliashberg equation. In any cases, two-phonon exchange together with extra pairing channels makes an important contribution to $`T_c`$ and it is the source of larger effective $`\lambda `$ in MgB2 than in other materials.

In summary, we have proposed a pairing mechanism in superconducting MgB2 based on the dynamic Jahn-Teller effect of the doubly degenerate $`p\sigma `$ hole states coupled to degenerate $`e_{2g}`$ phonons. The pseudo-spin contraint imposed on the motion of holes by the Jahn-Teller interactions is an important factor in understanding the dynamical properties of the system. The one-phonon and two-phonon couplings lead to effective pairing interactions with additional pseudo-spin channels. As results, the high $`T_c`$ observed in MgB2 is attributed mainly to the two-phonon exchange mechanism, and further the symmetry breaking hopping term plays a crucial role in producing two gaps which are responsible for the anomalous behaviors observed in experiment including the specific heat measurement.

We are grateful to Prof. T.W. Noh for helpful discussions. This work was supported by the BK21 Project of KRF. J.Y. acknowledges the support by the KOSEF through CSCMR. Y.-W.S and J.I. acknowledge the support by the KOSEF through the CNNC of the Sungkyunkwan University.


format: text

text:
# Violation of the Wiedemann-Franz Law in a Large-N Solution of the t-J Model

## Abstract

We show that the Wiedemann-Franz law, which holds for Landau Fermi liquids, breaks down in a large-n treatment of the t-J model. The calculated ratio of the in-plane thermal and electrical conductivities agrees quantitatively with experiments on the normal state of the electron-doped Pr2-xCexCuO4 ($`x=0.15`$) cuprate superconductor. The violation of the Wiedemann-Franz law in the uniform phase contrasts with other properties of the phase that are Fermi liquid like.

A recent experiment that measured the electrical and thermal conductivities of the copper-oxide superconductor Pr2-xCexCuO4 in its normal state found striking deviations from the Wiedemann-Franz lawhill . Simply stated, the Wiedemann-Franz law says that fermion quasiparticles transport both electrical and heat currents, with the ratio of the heat conductivity $`\kappa `$ to the electrical conductivity $`\sigma `$ given by:

$$\frac{\kappa }{\sigma T}=\frac{\pi ^2}{3}\left(\frac{k_B}{e}\right)^2$$
(1)
Ordinary Landau Fermi liquids respect the Wiedemann-Franz law, so deviations from it indicate the presence of non-Fermi liquid physics. The results of Ref. hill,  are therefore broadly consistent with other experimental evidence that points to non-Fermi liquid behavior in the cuprate phase diagram.

One possible interpretation of the breakdown of the Wiedemann-Franz law is that the quasiparticle fractionalizes into separate carriers of charge and spin. To see what effects such fractionalization might have, consider first heat transport in a system of weakly-interacting electrons, which could be viewed as a crude approximation to electrons in the highly overdoped region of the cuprate phase diagram. The electrons transport both charge and heat. The electrical conductivity is given by the Drude formula $`\sigma =ne^2\tau /m`$ and the thermal conductivity $`\kappa T`$ because while each quasiparticle carries fixed charge $`e`$, it only carries an energy of order the temperature. Next consider a model for the undoped cuprates: the Néel ordered antiferromagnetic insulator with zero electrical conductivity. Phonons and spinwaves transport heat, and as both excitations are bosonic in character with linear dispersions at low energy, each contributes similarly, yielding a thermal conductivity $`\kappa T^3`$.

Leaving aside these ordinary states of matter, consider the case in which the spins in the insulator, instead of ordering, fractionalize into spinons with an extended Fermi surface and a non-zero density of statesPWA ; BZA . Now spins transport heat in much the same way as charges do in the non-interacting electron system, with $`\kappa T`$. The Lorenz ratio is infinite, and remains significantly larger than one upon doping with holes or electrons.

Just this scenario is predicted in a large-n treatment of the t-J model on the square lattice. We follow the approach of Refs. AM, ; MA, ; Ted, . Implementing the single-occupancy constraint by introducing slave-boson operators $`b_i`$, the t-J model may be written:

$`H`$ $`=`$ $`t{\displaystyle \underset{<i,j>}{}}(c_i^\alpha b_ic_{j\alpha }b_j^{}+H.c.)`$ (2)
$`+`$ $`J{\displaystyle \underset{<i,j>}{}}(\stackrel{}{S}_i\stackrel{}{S}_j{\displaystyle \frac{1}{4}}n_in_j)`$
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}V(|\stackrel{}{r}_i\stackrel{}{r}_j|)n_in_j.`$
Number and spin operators are related to the electron creation and annihilation operators by $`n_ic_i^\alpha c_{i\alpha }`$ and $`\stackrel{}{S}_i(1/2)c_i^\alpha \stackrel{}{\sigma }_\alpha ^\beta c_{i\beta }`$ where there is an implicit sum over pairs of raised and lowered Greek indices. The single-occupancy constraint is now holonomic, $`b_i^{}b_i+c_i^\alpha c_{i\alpha }=1`$, with the physical meaning that only a hole, or a single electron, may occupy each site of the lattice. (In the following we generally refer to hole doping with the understanding that our calculations apply equally well to electron doped systems.) Included in Eq. 2 is the off-site Coulomb repulsion $`V(r)`$; however, in the uniform and staggered flux phases discussed below it plays no role other than to contribute an additive constant to the total energy.

Because spin-exchange involves no net flow of charge, electrical current only arises from the hopping term. The continuity equation relates the time rate of change of the occupancy on a given site to the lattice divergence of the current flowing into the site:

$$e\frac{dn_j(t)}{dt}=ie[n_j,H]=\underset{\widehat{e}=\widehat{x},\widehat{y}}{}\frac{j_{j,j+\widehat{e}}^ej_{j\widehat{e},j}^e}{a}$$
(3)
where $`a`$ is the lattice spacing between copper atoms. Thus

$$j_{j,j+\widehat{e}}^e=ieta(c_j^\alpha b_jc_{j+\widehat{e}\alpha }b_{j+\widehat{e}}^{}H.c.)$$
(4)
is the electric current flowing from site $`j`$ into site $`j+\widehat{e}`$ along the link connecting the two sites. We emphasize that neither the spin-spin exchange interaction $`J`$ nor the Coulomb interaction $`V(r)`$ appear in the expression for the electrical current. This result, which is a direct consequence of the gauge invariance of the spin-exchange and Coulomb interactions in our microscopic calculation, contrasts with that obtained recentlyyang  within the more phenomenological “d-density wave” picturenayak .

The heat current $`j_{j,j+\widehat{e}}^q`$ can be found by taking the time derivative of the Hamiltonian density:

$$\frac{dh_j(t)}{dt}=\underset{\widehat{e}=\widehat{x},\widehat{y}}{}\frac{j_{j,j+\widehat{e}}^qj_{j\widehat{e},j}^q}{a}$$
(5)
where

$$h_j=\underset{i=j\pm \widehat{x},\widehat{y}}{}(tc_i^\alpha b_ic_{j\alpha }b_j^{}+H.c.)+J\stackrel{}{S}_i\stackrel{}{S}_j$$
(6)
and the sum is only over the four sites $`i`$ that are nearest-neighbors of site $`j`$. We have dropped the $`n_in_j`$ interaction terms. These do not contribute to the DC thermal conductivity because $`n_i`$ remains unchanged in the presence of currents. However, in contrast to the electrical current, the heat current has contributions both from hopping, and from spin-exchange:

$`j_{j,j+\widehat{e}}^q`$ $`=`$ $`ta[c_j^\alpha b_j_t(c_{j+\widehat{e}\alpha }b_{j+\widehat{e}}^{})+H.c.]`$ (7)
$`+`$ $`{\displaystyle \frac{Ja}{2}}[c_j^\alpha (_tc_{j+\widehat{e}\alpha })c_{j+\widehat{e}}^\beta c_{j\beta }`$
$`+`$ $`c_j^\alpha c_{j+\widehat{e}\alpha }(_tc_{j+\widehat{e}}^\beta )c_{j\beta }].`$
The model generalizes from the physical case of $`n=2`$ (up and down spins) to arbitrary (even integer) values of $`n`$ by letting the Greek indices run over $`1,\mathrm{},n`$. In the large-n limit the spin-spin interaction factorizes in the particle-hole channel. Formally this factorization is implemented via a Hubbard-Stratonovich transformation within the functional integral approach. Complex-valued mean-fields along the bonds are introduced:

$$\chi _{ij}=\frac{J}{n}c_i^\alpha c_{j\alpha }.$$
(8)
The $`\chi _{ij}`$ fields function as the order parameter, and as they are spin-rotation invariant, there is no possibility of Néel or other spin order, and the mean-field Hamiltonian respects spin-rotational symmetry. Furthermore, at sufficiently low temperatures the (holon) boson fields condensecondense , and we may make the replacement $`b_i=b_i^{}=\sqrt{x}`$ where $`x`$ is the hole density. The mean-field Hamiltonian, which is exact in the $`n\mathrm{}`$ limit, then takes the form:

$$H_{MF}=\underset{}{}[\frac{n}{J}|\chi _{ij}|^2(tx+\chi _{ij})(c_j^\alpha c_{i\alpha }+H.c.)].$$
(9)
For parameters appropriate to the cuprates, $`t=0.44`$ eV (following Hybertsen et al.hybertsen ) and $`J=0.13`$ eV (obtained by Singh et al.singh ), the minimum energy configuration has $`\chi _{ij}`$ both real and constant when the doping exceeds $`x>0.12`$. There are no broken symmetries in this “uniform” phase. Upon suppressing dimerization with the addition of a biquadratic spin exchange interactionMA ; Vojta , a staggered flux (SF) phase with counter-circulating orbital currentsAM ; MA ; Ted ; john  occurs for $`x<0.12`$, that is, in the underdoped region of the phase diagram. (The biquadratic interaction which simultaneously exchanges four fermions disappears in the physical $`n=2`$ limit of up and down spins, and thus does not alter the physics of the t-J modelMA .) More realistic models include non-zero next- and third-nearest neighbor hopping amplitudesOKA  but these terms do not qualitatively affect the phase diagram or transport behavior.

The expressions for the two currents simplify in the large-n limit. As the boson operators may be replaced by the c-number $`\sqrt{x}`$, the electrical current becomes:

$$j_{j,j+\widehat{e}}^e=ietxa(c_j^\alpha c_{j+\widehat{e}\alpha }H.c.).$$
(10)
Upon further replacing the fermion bilinear operator $`c_i^\alpha c_{j\alpha }`$ with $`\frac{n}{J}\chi `$, the heat current also simplifieswhy :

$$j_{j,j+\widehat{e}}^q=(tx+\chi )a(c_j^\alpha _tc_{j+\widehat{e}\alpha }+H.c.).$$
(11)
As the heat current differs only by the $`(tx+\chi )`$ prefactor from that of a non-interacting tight-binding system, in the low-temperature limit the Lorenz ratio for in-plane transport is simply:

$$\frac{\kappa }{\sigma T}=\frac{\pi ^2}{3}\left(\frac{k_B}{e}\right)^2\left(\frac{tx+\chi }{tx}\right)^2.$$
(12)
Thus for any $`\chi 0`$ the ratio differs from unity, indicating a breakdown of Fermi liquid theory. Note that all details of the scattering mechanisms cancel out in the ratio. Direct calculation of the two conductivities in linear response shows that the integrals over momentum have identical form. Only the frequency integrals differ; hence for static impurities in the weak-scattering limit the Lorenz ratio is given by Eq. 12. We note that while the order parameter $`\chi `$ is perturbed by the application of external electric fields or thermal gradients, the perturbation in $`\chi `$ does not alter the DC responsecaroli .

In the SF phase the $`\chi `$-fields are complex numbersAM ; MA ; Ted , and the prefactor $`(tx+\chi )^2`$ should be replaced by $`|tx+\chi |^2`$. Specifically, for $`\chi =|\chi |\mathrm{exp}(i\theta )`$, with phase $`\theta `$, the Lorenz ratio generalizes toin\_prep :

$$\frac{L}{L_0}=\frac{(tx)^2+|\chi |^2+2tx|\chi |\mathrm{cos}(\theta )}{(tx)^2}$$
(13)
where $`L\kappa /(\sigma T)`$ and $`L_0(\pi ^2/3)(k_B/e)^22.45\times 10^8\mathrm{W}\mathrm{\Omega }\mathrm{K}^2`$ is the Lorenz number.

We turn now to a comparison of our predictions, Eqs. 12 and 13, with experiments. In Fig. 1 we plot the Lorenz ratio as a function of the doping. As expected, the ratio diverges as the insulating limit $`x0`$ is approached because the spinons transport only heat, not charge. In the opposite limit of large doping $`\chi 0`$, and the ratio approaches unity. Landau Fermi liquid theory is recovered in the dilute limit of widely spaced electrons. We emphasize that the uniform phase with $`\chi _{ij}`$ constant and real does not break any symmetries. It exhibits weak pseudogap behavior because, according to the mean-field equations, $`|\chi |`$ increases slightly in size at low temperatures, which in turn increases the quasiparticle bandwidth (see Eq. 9) and decreases the density of states (DOS)chung . For example, at a hole concentration of $`x=0.15`$, $`|\chi |=0.024`$ eV at $`T=500`$ K rising to $`|\chi |=0.026`$ eV at zero temperature; consequently the DOS drops by 2%. This contrasts with the strong pseudogap behavior of the SF phase which has a gap along most of the Fermi surface and which breaks time reversal invariance. In either phase, however, the fermionic quasiparticles are non-interacting in the $`n\mathrm{}`$ limit and hence behave as long-lived Landau quasiparticles such as those found in ordinary Fermi liquids.

We note that the Wiedemann-Franz law is strongly violated in s-wave superconducting states because while Cooper pairs carry charge, the condensate has no entropy. In a d-wave superconductor, quasiparticle excitations at the nodes result in a modified Wiedemann-Franz lawdurst . The violation that we find occurs in the normal state, and is a consequence of the spin-charge separation inherent in the large-n solution of the t-J model, and not of any incipient superconducting tendencies.

In Fig. 2 we plot the temperature dependence of the Lorenz ratio for three dopings at which transport experiments have been conducted: $`x=0.06`$, $`0.15`$, and $`0.26`$. For a single crystal of the La2-xSrxCuO4 material with hole doping $`x=0.06`$, the resistivity was measured upon suppressing the superconductivity by application of a 18 T magnetic field along the c-axistakeya . The thermal conductivity was, however, measured in the superconducting state, so it is not possible to extract a real Lorenz ratio. Nevertheless it is intriguing that $`L/L_05`$ at low temperatures, based on the numbers appearing in the inset to Fig. 3 of Ref. takeya, . This compares reasonably well with the theoretical value of $`4.1`$ seen in Fig. 2.

At optimal doping, experimentally available magnetic fields can only eliminate superconductivity in electron-doped compounds. In Ref. hill,  a 13 T field was applied to Pr2-xCexCuO4 at $`x=0.15`$ to access the normal state. The measured ratio of $`L/L_02`$ found at temperatures above $`0.3`$ K is again in reasonable quantitative agreement with the theoretical value of $`1.95`$. At temperatures below $`0.18`$ K, however, the experimentally determined ratio drops rapidly below one. We have no explanation for the observed behavior at the lowest temperaturessi .

Finally, in the highly overdoped regime Proust et al.proust  have studied the Tl2Ba2CuO6+δ material at a hole concentration of $`x0.26`$. Superconductivity was suppressed in a 13 T field, and $`L/L_0=0.99\pm 0.01`$ in good agreement with the Wiedemann-Franz law for Fermi liquids. The theoretical value of the ratio is $`1.5`$. We speculate that the persistence of non-Fermi liquid behavior at large doping in the mean-field theory is an artifact of the large-n approximation. Finite-n corrections could possibly restore Fermi liquid behavior in the overdoped region. At large doping $`|\chi |`$ is small compared to the effective hole hopping amplitude $`tx`$, so fluctuations in $`\chi `$ may be expected to be relatively more important than at low doping.

In summary we have shown that the Wiedemann-Franz law is violated in a mean-field treatment of the t-J model. Our analysis, which holds for weak scattering, is exact in the $`n\mathrm{}`$ limit. The Lorenz ratio is significantly larger than one both in the uniform phase ($`x>0.12`$) and in the SF phase ($`x<0.12`$). The theoretical prediction is in reasonably good quantitative agreement with existing experimental measurements on the cuprate materials.

Note added:  After this work was completed a paper by Kim and Carbotte (KC) appearedkim  that examined the Wiedemann-Franz law within the context of the phenomenological d-density wave picture. There are several differences between their work and ours. The main difference is that we study both the uniform phase which has no broken symmetries, and the SF or d-density phase with time-reversal breaking counter-circulating currents. We find that the Lorenz ratio is significantly larger than one in both phases. Furthermore, at low temperatures KC find only small deviations from the Wiedemann-Franz law. This is due in part to the fact that their d-density order parameter (the analog of our $`\chi _{ij}`$) was chosen to be purely imaginary (equivalent to setting $`\theta =\pi /2`$ in our Eq. 13) and also because their kinetic energy is not rescaled by the slave-boson doping factor, $`x`$, as it is in our microscopic analysis of the t-J model. KC also find a large temperature variation in the Lorenz ratio for the case of strong scattering because the quasiparticle lifetime has a strong frequency dependence.

We thank John Fjærestad and Louis Taillefer for helpful comments. This work was supported in part by the NSF under grant No. DMR-9712391.


text:
# Traveling waves of infection in the Hantavirus epidemics

## I Introduction

The Hantavirus zoonosis is a persistent problem in many regions of the world mills99a . Each species of Hantavirus is almost exclusively associated with a single rodent species from which, eventually, humans may result infected. The disease can range from mild to very severe, with a mortality in excess of 50%. Such is the case of the Hantavirus Pulmonary Syndrome (HPS) produced by the virus associated with the deer mouse (*Peromyscus maniculatus*). In the North American Southwest, indeed, a serious outbreak of this disease in 1993 led to the identification of the virus and its association with the deer mouse parmenter99 . Since then, enormous effort and resources have been devoted to the understanding of the ecology and the epidemiology of the virus-mouse association, with the ultimate goal of being able to predict and prevent the risk for the human population mills99b ; cparmenter98 .

As observed by Yates et al. yates02 , the role of the environment seems to be determinant in the prevalence of the infection within the mouse population, and the related incidence of the HPS. Its has been observed that the disease can disappear completely from a local population during times of adverse environmental conditions, only to reappear sporadically mills99b ; cparmenter99 . Besides, there are indications of focality of the infection in “refugia” of the mouse population. Both phenomena are most probably related, the refugia acting as reservoirs of the virus during times when the infection has disappeared from most of the landscape. When environmental conditions change, the infection spreads again from these refugia.

In a recent work abramson2002 , Abramson and Kenkre have shown that a simple epidemic model is able to display qualitatively these behaviors, as the result of a bifurcation of the equilibrium states of the system as controlled by the carrying capacity of the medium. The purpose of the present paper is to analyze the dynamics of simple traveling waves, as a model of the mechanisms in which an epidemic wave of Hantavirus might propagate into a previously uninfected region.

## II Spatially extended model

The model of Ref. abramson2002  is a mean-field continuous model which has been intentionally kept simple enough to facilitate the comprehension of the basic mechanisms, and at the same time to incorporate as many known facts about the ecology and epidemiology of the biological system as possible. The reader is referred to abramson2002  for a detailed discussion, which we summarize here. It is known that the virus does not affect properties such as the mortality of the mice, so that no appreciable difference in death rate, for example, is to be expected between susceptible and infected mice. It is also not transmitted to newborns, so that all mice are born susceptible. The infection is transmitted from mouse to mouse through individual contacts, presumably during fights. More general facts of the ecology of *Peromyscus* indicate that adults occasionally shift their home range to nearby locations, in particular if these are vacant stickel68 ; vessey87 . This enables us to model the transport of mice as a diffusion process. Finally, intra-species competition for resources indicate a saturated population growth, which has been observed to be of a logistic form in the laboratory terman68 . Logistic growth is also a well established metaphor of the dynamics of a self-limitating population murray .

For the sake of simplicity, assume further that the only population structure is the division of the whole population into susceptible and infected mice, denoted by $`M_S`$ and $`M_I`$ respectively. With these ingredients, the model is described by the following equations:

$`{\displaystyle \frac{M_S}{t}}`$ $`=`$ $`bMcM_S{\displaystyle \frac{M_SM}{K(x,t)}}aM_SM_I+D^2M_S,`$ (1)
$`{\displaystyle \frac{M_I}{t}}`$ $`=`$ $`cM_I{\displaystyle \frac{M_IM}{K(x,t)}}+aM_SM_I+D^2M_S.`$ (2)
Observe that the carrying capacity $`K(x,t)`$, containing the most direct effect of the environment on the mouse population, is allowed, in our model, a spatial and a temporal variation, to accommodate for a diversity of habitats and temporal phenomena. The latter comprise the yearly variation due to seasonality, but also sporadic fluctuations such as droughts and El Niño effects.

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

$$\frac{M}{t}=(bc)M\left(1\frac{M}{(bc)K}\right)+D^2M.$$
(3)
This is Fisher’s equation, originally proposed as a deterministic model of the spread of a favored gene in a population fisher36 , and which eventually became a standard model for a self regulated field in a diversity of situations murray ; abramson2001 .

In Ref. abramson2002  it was shown that, as a function of $`K`$, the system undergoes a bifurcation between a stable state with only susceptible mice (and $`M_I=0`$) to a stable state with both positive populations. The value of the critical carrying capacity is a function of the parameters in the following way:

$$K_c=\frac{b}{a(bc)}.$$
(4)
This critical value does not depend on $`D`$, and the same bifurcation is observed either in a space-independent system ($`D=0`$) or in a homogeneous extended one in the presence of diffusion. In an inhomogeneous situation, for moderate values of the diffusion coefficient, the infected subpopulation remains restricted to those places where $`K>K_c`$, becoming extinct in the rest. During times of adverse environmental conditions, these regions become isolated and constitute the observed refugia. Figure 1 shows a typical situation of this phenomenon. A simulated, albeit realistic, landscape of $`K(x)`$ has been derived from satellite images in Northern Patagonia. The carrying capacity is supposed proportional to the vegetation cover, and results highest along a river, clearly inferred from the density plots. These show the distribution of the populations of susceptible and infected mice. It can be seen that susceptible mice cover most of the range. Meanwhile, the infected population has become extinct in most of it, and persists only in some of the places of highest $`K`$. The distributions shown in Fig. 1 constitute a stable equilibrium of the system, found by numerical resolution of Eqs. (1)-(2) from an arbitrary initial condition, and with zero-current boundary conditions.

## III Traveling waves

When conditions in the landscape change, how do the infected phase evolve from the refugia, retracting from or, more importantly, invading previously uninfected regions? this is the primary question we address in the present paper. Fisher’s equation (3) has found wide applicability for the description of the formation and propagation of spatio-temporal patterns. Traveling wave solutions of Fisher’s equation paradigmatically show how a stable phase invades an unstable one in the form of a front propagating at a certain speed.

There is no reason to suppose *a priori* that the two waves, susceptible and infective, will travel at the same speed. Accordingly, we use an ansatz which incorporates two independent traveling waves. In one dimension, $`z_1=xv_St`$ in Eq. (1) and $`z_2=xv_It`$ in Eq. (2). This gives the following second-order system of ordinary differential equations:

$`D{\displaystyle \frac{d^2M_S(z_1)}{dz_1^2}}+v_S{\displaystyle \frac{dM_S(z_1)}{dz_1}}+f(M_S,M_I)`$ $`=`$ $`0,`$ (5)
$`D{\displaystyle \frac{d^2M_I(z_2)}{dz_2^2}}+v_I{\displaystyle \frac{dM_I(z_2)}{dz_2}}+g(M_S,M_I)`$ $`=`$ $`0,`$ (6)
where $`v_S`$ and $`v_I`$ are the speeds of the susceptible and infected waves respectively, and $`f`$ and $`g`$ are the non-diffusive term in (1)-(2).

There are two interesting scenarios for these waves. In the first one, a large part of the system is initially at a state of low carrying capacity, below $`K_c`$, and consequently the population consists of uninfected mice only, at the stable equilibrium. Let us suppose that this region is in contact with a refugium. If environmental changes occur, and the whole region finds itself at a value of the carrying capacity $`K>K_c`$, the population will be out of equilibrium. Two processes will occur simultaneously: the population of susceptible mice will evolve towards a new equilibrium, and a wave of infected mice will advance from the refugium, invading the susceptible population. The speed of this wave can be calculated from the stability analysis of the equilibrium states, requiring that the unstable infected mice density does not oscillate below zero. This unstable equilibrium is $`M_S^{}=K(bc)`$, $`M_I^{}=0`$, and a linear stability analysis of the system (5)-(6) provides the following four eigenvalues:

$`\lambda _{1,2}`$ $`=`$ $`{\displaystyle \frac{v\pm \sqrt{v^2+4D(bc)}}{2D}},`$ (7)
$`\lambda _{3,4}`$ $`=`$ $`{\displaystyle \frac{v\pm \sqrt{v^2+4D[baK(bc)]}}{2D}}.`$ (8)
The requirement that $`M_I(z)`$ does not oscillate below 0 imposes a restriction to the radical in Eq. (8), from which we find the following expression for the speed of the traveling wave:

$$v2\sqrt{D\left[b+aK(bc)\right]}.$$
(9)
An example of such a wave is shown in Fig. 2a, found by a numerical integration of Eqs. (1) and (2) in 1 dimension.

The second interesting scenario corresponds to a system which is initially empty of mice, both susceptible and infected. This situation is always unstable within the validity of our simple model, but it is certainly a biological possibility. Consider a system with $`K>K_c`$, and with $`M_S=M_I=0`$ in almost all of its range, but in contact with a refugium in equilibrium. A wave of both mice populations will develop invading the empty region. In fact, the total population will display just the behavior of a traveling wave of Fisher’s equation. This wave will be composed of two fronts, susceptible and infected respectively, with a delay of the latter with respect to the former. As before, we can find the speeds of these two waves from a stability analysis around the corresponding equilibria. The leading wave propagates into the null equilibrium, $`M_S^{}=M_I^{}=0`$, to which the following eigenvalues correspond:

$`\mu _{1,2}`$ $`=`$ $`{\displaystyle \frac{v_I\pm \sqrt{v_I^2+4Dc}}{2D}},`$ (10)
$`\mu _{3,4}`$ $`=`$ $`{\displaystyle \frac{v_S\pm \sqrt{v_S^24D(bc)}}{2D}}.`$ (11)
In this situations, we require that $`M_S(z)`$ does not oscillate below 0, and Eq. (11) provides the restriction on the speed of the susceptible front:

$$v_S2\sqrt{D(bc)},$$
(12)
which is, naturally, the same result as for Fisher’s equation. The second front, developed when part of the quasi-stable population of susceptible mice is converted into infected, evolves from the equilibrium $`M_S^{}=K(bc)`$, $`M_I^{}=0`$, as in the previous scenario. Consequently, the same linear stability analysis apply, and from Eq. (8) we find a speed analogous to Eq. (9). The front of infected move behind the susceptible one at a speed:

$$v_I2\sqrt{D[b+aK(bc)]},$$
(13)
which, unlike $`v_S`$, does depend on the contagion rate $`a`$ and the carrying capacity $`K`$. Figure 2b shows such a situation. The density of susceptible mice rises from zero and lingers near the positive unstable equilibrium before tending to the stable one. It is remarkable that a delay exists between the two fronts, even when no such effect was explicitly considered in the dynamics (such as incubation time, or age dependence). Such delays have been observed in some populations and rationalized in different ways (see morgan00  or mills99b  for a synthesis).

Even though Eqs. (9), (12) and (13) give only a lower bound to the speed of propagation of the fronts, and allow a continuous of speeds above this, in real situations only the lower bound is expected to be observed as a stationary solution. Higher speeds may, however, play a role in transient situations whose relevance in far from equilibrium systems such as real mice in the wild, subject to a fluctuating environment, cannot be underestimated.

The different functional dependence of $`v_S`$ and $`v_I`$ on the parameters of the system (Eqs. (12)-(13)) indicates that two regimes are possible. Indeed, when $`v_I<v_S`$ the front of infection lags behind the front of susceptible at a delay $`\mathrm{\Delta }`$ that increases linearly in time: $`\mathrm{\Delta }(t)=(v_Sv_I)t`$. Elementary manipulation of Eqs. (12) and (13) shows that this occurs whenever the carrying capacity satisfies:

$$K_c<K<K_0\frac{2bc}{a(bc)}$$
(14)
where $`K_0`$ is a new critical carrying capacity. At $`K=K_0`$ the delay becomes effectively constant. For values of $`K`$ greater than $`K_0`$, the velocities $`v_I`$ and $`v_S`$, calculated from linear considerations around the equilibria, satisfy $`v_I>v_S`$. This regime is clearly unphysical in a stationary situation, since the front of susceptible necessarily precedes the infected one. It could be realizable and relevant in transient situations, that will be analyzed elsewhere. From numerical resolution of the system, we can observe that $`v_Iv_S`$ and the delay tends to a constant value, whenever $`K>K_0`$. Figure 3 shows the temporal evolution of the delay in the two regimes $`K<K_0`$ and $`K>K_0`$, as well as in the critical case $`K=K_0`$, where it is seen to increase as $`\mathrm{\Delta }t^{1/4}`$. It can be seen that there is a transient time in both regimes, that gets progressively longer as $`K`$ approaches $`K_0`$ either from above or from below.

With the assumption that $`v_S=v_I=v`$, it is possible to perform a piecewise linearization of Eqs. (5)-(6) and find an approximate analytical expression for the front shapes shown in Fig. 2b and, consequently, for the delay $`\mathrm{\Delta }`$ in the stationary state. The details of the calculation can be found in the Appendix. The main result is the following expression for the delay, for the case of equal speeds of the two classes of mice:

$$\mathrm{\Delta }=\frac{\sqrt{D}}{i\sqrt{(bc)a(KK_0)}}\mathrm{log}(w_1w_2),$$
(15)
where $`w_1`$ and $`w_2`$ are complex numbers of unit modulus that depend on $`a`$, $`b`$ and $`c`$, so that the logarithm is effectively twice the phase difference between them. When $`KK_0^+`$, the arguments of $`w_1`$ and $`w_2`$ tend to $`\pi `$ and 0 respectively, so that the following critical behavior is predicted by the linear approximation:

$$\mathrm{\Delta }\sqrt{D}[(bc)a(KK_0)]^{1/2},\text{ when }KK_0^+.$$
(16)
We have analyzed the dependence of $`\mathrm{\Delta }`$ on $`K`$ by means of numerical resolution of the full system. In Fig. 4 we plot the asymptotic value of $`\mathrm{\Delta }`$ as a function of $`KK_0`$ for a variety of system parameters. The behavior is found to be

$$\mathrm{\Delta }\sqrt{D}[(bc)a(KK_0)]^\alpha $$
(17)
for values of $`K`$ immediately above $`K_0`$. The exponent $`\alpha `$ is parameter-independent, and its value, calculated from the simulations is, approximately $`0.388`$. There is a discrepancy in the exponent found numerically in the fully nonlinear system and the one found in the linearized approximation, which shows the limitation of the linearized solution.

This critical behavior of the delay persists while the wavefronts are of the kind shown in Fig. 2b. However, if $`K`$ continues to grow the equilibrium value of the infected population turns greater than the equilibrium value of the susceptible. There is a gradual crossover to a situation where most of the population becomes infected. The value of the carrying capacity at which this happens can be estimated from the known equilibria, $`M_S^{}=M_I^{}`$, giving:

$$K_1=\frac{2b}{a(bc)}.$$
(18)
Indeed, the numerical results show that the power law decay of the delay as a function of $`K`$ starts to flatten for values of $`K>K_1`$ and reaches a minimum which shows the same $`a^1`$ dependence that $`K_1`$ does. Greater values of $`K`$ are probably unrealistic in the arid and semi-arid habitats of the Southwest. It is nevertheless interesting to point out that the interaction of the two fronts in this regime results in an *increase* of the delay as a function of $`K`$. This increase is also algebraic:

$$\mathrm{\Delta }(KK_0)^\beta ,$$
(19)
with $`\beta 0.16`$, as found in the numerical calculations.

## IV Conclusion

We have analyzed the propagation of traveling fronts in 1 dimension in a simple model of the ecology and epidemiology of the Hantavirus in deer mouse. We have found that, when a mouse-free region is in contact with an infected region in equilibrium, two waves propagate into the empty region. The first one is a wave of susceptible mice. A wave of infected mice propagates behind it with a certain delay. Two regimes of propagation exist, controlled by the environmental parameter $`K`$. If $`K_c<K<K_0`$, the lag between the two fronts increase linearly in time. Conversely, if $`K>K_0`$, the two fronts propagate at the same speed and the delay depends critically on the difference $`KK_0`$.

The occurrence of this double regime may be of relevance for the control of the propagation of an epidemic wave. Indeed, the controlled reduction of $`K`$ ahead of a propagating wave seems the most effective mean of stopping or reducing its advance. Ideally, the carrying capacity should be reduced below $`K_0`$, to ensure the complete extinction of the infection. However, if such a reduction is not feasible, the fact that $`K_0>K_c`$ provides an alternative: a reduction of the carrying capacity below $`K_0`$ would make the wave of infection start to lag more and more behind the wave of healthy mice. Possible implementations of these strategies, based on the propagation of waves in the presence of “barriers,” will be analyzed in detail elsewhere.

The existence of dynamical phenomena such as these traveling fronts also opens the interesting possibility of subjecting our predictions to experimental verification. Controlled experiments of front propagation could be possible in the Sevilleta LTER facility, that the University of New Mexico operates near Socorro, NM lter . Measurements of uncontrolled mice populations along lines radiating from the refugia of infection will also provide evidence of the propagation mechanisms. The observation of these in real mice population will provide a valuable source of data to assign realistic values to the parameters of the mathematical model.

###### Acknowledgements.

V. M. K. and G. A. acknowledge many discussions with Fred Koster and Jorge Salazar from which we learnt much regarding the peculiarities of the Hantavirus. We also thank Greg Glass, Karl Johnson, Luca Giuggioli and María Alejandra Aguirre for discussions. 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. A part of the numerically intensive computations was carried out on the Albuquerque High Performance Computing Center facilities.

## V Appendix: linearized solutions

For the purpose of finding approximate solutions of the waves profiles, it is better to replace the system (5)-(6) with one involving $`M`$ and $`M_I`$ instead:

$`DM^{\prime \prime }+vM^{}+(bc)M{\displaystyle \frac{M^2}{K}}`$ $`=`$ $`0,`$ (20)
$`DM_I^{\prime \prime }+vM_I^{}+q(M)M_IaM_I^2`$ $`=`$ $`0,`$ (21)
where $`q(M(z))=c+aM(z)M(z)/K`$ and both speeds are assumed equal, as observed in numerical resolutions for the regime $`K>K_0`$. Primes denote differentiation with respect to $`z`$. The reason for using this system is that the equation for $`M(z)`$, Eq. (20), being closed, can be solved independently of $`M_I`$. Its solution can be used then into Eq. (21) as a $`z`$-dependent parameter and solve for $`M_I(z)`$.

Linearized solutions for the traveling waves of Fisher’s equation (20) are well known, and essentially consist of two exponentials, matched smoothly at $`z=0`$, representing a front that travels to the right at a speed $`v2\sqrt{D(bc)}`$. Such a function needs a further simplification in order to solve Eq. (21). We approximate $`M(z)`$ with a step function discontinuous at $`z=0`$: $`M(z)=M^{}=K(bc)`$ if $`z<0`$ and $`M(z)=0`$ if $`z>0`$. Consequently, we have that $`q(z)=aM_I^{}=b+aK(bc)`$ if $`z<0`$ and $`q(z)=c`$ if $`z>0`$. We also linearize on both sides of the center value $`M_I(\mathrm{\Delta })=M_I^{}/2`$, which defines $`\mathrm{\Delta }`$ as the delay between the two waves. Figure 5 shows the geometry of the procedure. The linearized equation for $`M_I`$ breaks into three regimes:

$`DM_I^{\prime \prime }+vM_I^{}a(M_IM_I^{})`$ $`=`$ $`0\text{ if }z<\mathrm{\Delta },`$ (22)
$`DM_I^{\prime \prime }+vM_I^{}+aM_I^{}M_I`$ $`=`$ $`0\text{ if }z[\mathrm{\Delta },0],`$ (23)
$`DM_I^{\prime \prime }+vM_I^{}CM_I`$ $`=`$ $`0\text{ if }0<z.`$ (24)
The solutions of these are respectively:

$`M_1(z)`$ $`=`$ $`{\displaystyle \frac{b}{a}}+K(bc)+a_1e^{\lambda z}\text{ if }z<\mathrm{\Delta },`$ (25)
$`M_2(z)`$ $`=`$ $`b_1e^{\mu _+z}+b_2e^{\mu _{}z}\text{if }z[\mathrm{\Delta },0],`$ (26)
$`M_3(z)`$ $`=`$ $`c_1e^{\nu z}\text{if }0<z,`$ (27)
where:

$`\lambda `$ $`=`$ $`{\displaystyle \frac{\sqrt{bc+a}\sqrt{bc}}{\sqrt{D}}},`$ (28)
$`\mu _\pm `$ $`=`$ $`\sqrt{{\displaystyle \frac{bc}{D}}}\pm i\sqrt{{\displaystyle \frac{bc}{D}}a(KK_0)},`$ (29)
$`\nu `$ $`=`$ $`{\displaystyle \frac{\sqrt{bc}+\sqrt{b}}{\sqrt{D}}},`$ (30)
and the terms that diverge in infinity have already been cancelled. It is clear that $`\lambda >0`$ and $`\nu <0`$, as required for the solution not to diverge in $`\mathrm{}`$. The real part of $`\mu _\pm `$ is always positive, and there is an imaginary part when $`K>K_0`$.

The solutions $`M_1`$, $`M_2`$ and $`M_3`$ need to be matched smoothly at $`z=\mathrm{\Delta }`$ and $`z=0`$, in the following way ($`r=bc`$ from now on):

$`M_1(\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{b}{2a}}+{\displaystyle \frac{kr}{2}},`$ (31)
$`M_2(\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{b}{2a}}+{\displaystyle \frac{kr}{2}},`$ (32)
$`M_1^{}(\mathrm{\Delta })`$ $`=`$ $`M_2^{}(\mathrm{\Delta }),`$ (33)
$`M_2(0)`$ $`=`$ $`M_3(0),`$ (34)
$`M_2^{}(0)`$ $`=`$ $`M_3^{}(0).`$ (35)
These are five equations with five unknowns that can be solved without difficulty:

$`a_1`$ $`=`$ $`{\displaystyle \frac{b/a+kr}{2}}e^{\lambda \mathrm{\Delta }},`$ (36)
$`b_1`$ $`=`$ $`{\displaystyle \frac{(\lambda +\mu _{})}{2(\mu _+\mu _{})}}(krb/a)e^{\mu _+\mathrm{\Delta }},`$ (37)
$`b_2`$ $`=`$ $`{\displaystyle \frac{(\lambda +\mu _+)}{2(\mu _+\mu _{})}}(krb/a)e^{\mu _{}\mathrm{\Delta }},`$ (38)
$`c_1`$ $`=`$ $`b_1+b_2,`$ (39)
and:

$$\mathrm{\Delta }=\frac{1}{\mu _+\mu _{}}\mathrm{log}\left[\frac{(\lambda +\mu _+)(\mu _{}\nu )}{(\lambda +\mu _{})(\mu _+\nu )}\right].$$
(40)
Expressions (25)-(27), (28)-(30) and (36)-(40) define the linearized solution for $`M_I`$. An example of such a front wave is shown in Fig. 5.

Using expressions (28)-(30) in (40), it easy to see find:

$$\mathrm{\Delta }=\frac{\sqrt{D}}{i\sqrt{(bc)a(KK_0)}}\mathrm{log}(w_1w_2),$$
(41)
where

$`w_1`$ $`=`$ $`{\displaystyle \frac{\lambda +\mu _+}{\lambda +\mu _{}}}e^{i\varphi _1},`$ (42)
$`w_2`$ $`=`$ $`{\displaystyle \frac{\mu _{}\nu }{\mu _+\nu }}e^{i\varphi _2},`$ (43)
where $`\varphi _1`$ and $`\varphi _2`$ are respectively the arguments of the complex numbers $`\lambda +\mu _+`$ and $`\mu _+\nu `$. It is easy to see that $`\varphi _1`$ and $`\varphi _2`$ do not depend on $`D`$, so that the logarithm only corrects the dependence of $`\mathrm{\Delta }`$ on $`(bc)a(KK_0)`$. It can be also easily observed that $`w_1`$ lies in the second quadrant, and $`w_2`$ in the fourth, and that $`\text{Im}w_1=\text{Im}w_2`$. When $`KK_0^+`$, the imaginary parts tend to zero and the phase difference between both tend to $`\pi `$. Consequently the leading behavior of $`\mathrm{\Delta }`$ results:

$$\mathrm{\Delta }\sqrt{D}[(bc)a(KK_0)]^{1/2}.$$
(44)
This result from the piecewise linearization does not agree with that found numerically in the nonlinear system.


text:
# Untitled Document

Nature, Vol 416, 59 (7 March 2002)

Massive star formation in 100,000 years from turbulent and pressurized molecular clouds

Christopher F. McKee1,2 and Jonathan C. Tan2,3

1. Department of Physics, University of California, Berkeley, CA 94720, USA.

2. Department of Astronomy, University of California, Berkeley, CA 94720, USA.

3. Princeton University Observatory, Peyton Hall, Princeton, NJ 08544, USA.

Massive stars (with mass $`m_{}>8\mathrm{M}_{}`$) are fundamental to the evolution of galaxies, because they produce heavy elements, inject energy into the interstellar medium, and possibly regulate the star formation rate. The individual star formation time, $`t_f`$, determines the accretion rate of the star; the value of the former quantity is currently uncertain by many orders of magnitude1,2,3,4,5,6, leading to other astrophysical questions. For example, the variation of $`t_f`$ with stellar mass dictates whether massive stars can form simultaneously with low-mass stars in clusters. Here we show that $`t_f`$ is determined by conditions in the star’s natal cloud, and is typically $`10^5\mathrm{yr}`$. The corresponding mass accretion rate depends on the pressure within the cloud—which we relate to the gas surface density—and on both the instantaneous and final stellar masses. Characteristic accretion rates are sufficient to overcome radiation pressure from $`100\mathrm{M}_{}`$ protostars, while simultaneously driving intense bipolar gas outflows. The weak dependence of $`t_f`$ on the final mass of the star allows high- and low-mass star formation to occur nearly simultaneously in clusters.

Massive stars form in dense molecular clumps inside molecular clouds7. These regions are highly turbulent and are in approximate virial equilibrium, with comparable values of the gravitational energy and the kinetic energy. Observed star forming regions have a virial parameter $`\alpha _{\mathrm{vir}}=5\sigma ^2R_{\mathrm{cl}}/GM_{\mathrm{cl}}`$ of order unity8, where $`\sigma ^2`$ is the mass-averaged one-dimensional velocity dispersion, $`R_{\mathrm{cl}}`$ the radius, and $`M_{\mathrm{cl}}`$ the mass of the clump. The regions of high-mass star formation studied in ref. 7 are characterized by masses $`M_{\mathrm{cl}}3800\mathrm{M}_{}`$ and radii $`R_{\mathrm{cl}}0.5\mathrm{pc}`$. The corresponding mean column density and visual extinction are $`\mathrm{\Sigma }1\mathrm{g}\mathrm{cm}^2`$ and $`A_V=(N_\mathrm{H}/2\times 10^{21}\mathrm{cm}^2)\mathrm{mag}=214\mathrm{\Sigma }\mathrm{mag}`$, with a dispersion of a factor of a few. These column densities are far greater than those associated with regions of low-mass star formation, but are comparable to the central surface density in the Orion Nebula Cluster9, which is about $`1\mathrm{g}\mathrm{cm}^2`$. These regions have very high mean pressures $`\overline{P}`$:

$$\frac{\overline{P}}{k_B}=\frac{3M_{\mathrm{cl}}\sigma ^2}{4\pi R_{\mathrm{cl}}^3k_B}=\left(\frac{3\pi f_{\mathrm{gas}}\alpha _{\mathrm{vir}}}{20k_B}\right)G\mathrm{\Sigma }^22.28\times 10^8f_{\mathrm{gas}}\alpha _{\mathrm{vir}}\mathrm{\Sigma }^2\mathrm{K}\mathrm{cm}^3,$$
$`(1)`$
where $`k_B`$ is Boltzmann’s constant and $`f_{\mathrm{gas}}`$ is the fraction of the cloud’s mass that is in gas, as opposed to stars. They have power-law density profiles10 $`\rho r^{k_\rho }`$, with $`k_\rho 1.5\pm 0.5`$.

Molecular clouds are observed to be inhomogeneous on a wide range of scales, and numerical simulations show that this is a natural outcome of supersonic turbulence11. Dense, self-gravitating inhomogeneities that are destined to become stars are termed cores. We assume that the distribution of core masses determines the resulting distribution of stellar masses (the initial mass function) and we take this distribution as given. Because of the disruptive effects of protostellar outflows, only a fraction $`ϵ_{\mathrm{core}}`$ of the total core mass $`M_{\mathrm{core}}`$ ends up in the star: $`m_f=ϵ_{\mathrm{core}}M_{\mathrm{core}}`$, where $`m_f`$ is the final stellar mass (for binary and other multiple star systems, $`m_f`$ is the total mass of stars in the system). To estimate $`ϵ_{\mathrm{core}}`$ for massive stars, we assume that their outflows are scaled versions of the outflows from low-mass stars, which is qualitatively consistent with observation12. For low-mass stars, protostellar outflows typically eject 25-75% of the core mass13. A similar analysis applied to high-mass stars yields comparable results (J.C.T. & C.F.M., in prep.), so we adopt $`ϵ_{\mathrm{core}}=0.5`$ as a typical value.

Our fundamental assumption is that star-forming clumps and the cores embedded within them are each part of a self-similar, self-gravitating turbulent structure that is virialized—that is, in approximate hydrostatic equilibrium—on all scales above the thermal Jeans mass. The density and pressure are then power laws in radius (for the pressure, $`Pr^{k_P})`$, so that the cloud is a polytrope with $`P\rho ^{\gamma _p}`$. In hydrostatic equilibrium14, $`k_\rho =2/(2\gamma _p)`$ and $`k_P=\gamma _pk_\rho =2\gamma _p/(2\gamma _p)`$. Let $`c(P/\rho )^{1/2}r^{(1\gamma _p)/(2\gamma _p)}`$ be the effective sound speed; if the pressure is dominated by turbulent motions, then $`c`$ is proportional to the line width. Molecular clouds and cloud cores satisfy a line width-size relation in which $`c`$ increases outwards15, corresponding to $`\gamma _p<1`$. The equation of hydrostatic equilibrium gives $`M=k_Pc^2r/G`$ for the mass inside $`r`$, and $`\rho =Ac^2/(2\pi Gr^2)`$ for the density at $`r`$, where $`A=\gamma _p(43\gamma _p)/(2\gamma _p)^2`$. It is then immediately possible to determine the properties of a core in terms of the pressure at its surface and the mass of the star that will form within it (see Fig. 1). The radius of a core is then $`r_s=0.074(m_f/30\mathrm{M}_{})^{1/2}\mathrm{\Sigma }^{1/2}\mathrm{pc}`$; recall that the typical clump observed in ref. 7 has a radius of 0.5 pc, which is set by the angular resolution of those observations. The mean density in a core is $`\overline{n}_\mathrm{H}=1.0\times 10^6(m_f/30\mathrm{M}_{})^{1/2}\mathrm{\Sigma }^{3/2}\mathrm{cm}^3`$, and the r.m.s. velocity dispersion is $`1.65(m_f/30\mathrm{M}_{})^{1/4}\mathrm{\Sigma }^{1/4}\mathrm{km}\mathrm{s}^1`$. The thermal sound speed at temperature $`T`$ is $`0.3(T/30\mathrm{K})^{1/2}\mathrm{km}\mathrm{s}^1`$, and so a high-mass core is dominated by supersonic turbulence and should be very clumpy.

We now consider the timescale for a star to form in such a core. On dimensional grounds, we expect the protostellar accretion rate to be16:

$$\dot{m}_{}=\varphi _{}\frac{m_{}}{t_{\mathrm{ff}}},$$
$`(2)`$
so long as radiation pressure is not important. Here $`m_{}`$ is the instantaneous stellar mass, $`t_{\mathrm{ff}}=[3\pi /(32G\rho )]^{1/2}`$ is the free-fall time and $`\varphi _{}`$ is a dimensionless constant of order unity. This equation could be violated in the sense that $`\varphi _{}1`$ only in the unlikely case that the star forms from a coherent spherical implosion; if the star formation is triggered by an external increase in pressure, $`\varphi _{}`$ could be increased somewhat, but deviations from spherical symmetry in the triggering impulse and in the protostellar core will generally prevent $`\varphi _{}`$ from becoming too large. It could be violated in the opposite sense that $`\varphi _{}1`$ if the core is magnetically dominated, so that collapse could not begin until the magnetic field diffused out of the core. However, magnetic fields are not observed to be dominant in molecular cores17. Equation (2) implies that the accretion rate, and thus the star formation time $`t_fm_{}/\dot{m}_{}`$, depends weakly on the properties of the ambient medium, $`\dot{m}_{}\rho ^{1/2}`$.

We now show that if the collapse is spherical and self-similar, then $`\varphi _{}`$ is quite close to unity provided that radiation pressure does not disrupt the flow. Although we assume the collapse is spherical far from the star, it will naturally proceed via a disk close to the star owing to the angular momentum of the accreting material; we assume this does not limit the flow of matter onto the star as otherwise the disk would become very massive and gravitationally unstable. Shielding by the disk reduces the importance of radiation pressure on the accretion flow, and allows the formation of massive stars provided the accretion rate is sufficiently large18. Under the assumption of a polytropic structure in hydrostatic equilibrium, equation (2) implies:

$$\dot{m}_{}=\varphi _{}ϵ_{\mathrm{core}}\frac{4}{\pi \sqrt{3}}k_PA^{1/2}\frac{c^3}{G},$$
$`(3)`$
where the value of $`\rho `$ entering $`t_{\mathrm{ff}}`$ is evaluated at the radius in the initial cloud that just encloses the gas that goes into the star when its mass is $`m_{}`$. For the isothermal case ($`\gamma _p=1`$) and for $`ϵ_{\mathrm{core}}=1`$, this reduces to the classic result of Shu19 with $`\varphi _{}=0.975\pi \sqrt{3}/8=0.663`$. McLaughlin & Pudritz2 show that the accretion rate $`\dot{m}_{}`$ is proportional to $`t^{33\gamma _p}`$, so that for $`\gamma _p<1`$, as observed, the accretion rate accelerates2,4. As discussed in ref. 2, termination of the accretion breaks the self-similarity once the expansion wave encloses $`m_f`$, which occurs at $`t0.45t_f`$. Thereafter, the relation $`\dot{m}_{}t^{33\gamma _p}`$ becomes approximate, but the approximation should be reasonable good2. The star-formation time is then given by:

$$t_f=\frac{(43\gamma _p)}{\varphi _{}}t_{\mathrm{ff}}.$$
$`(4)`$
Evaluating $`\varphi _{}`$ (see Fig. 1), we conclude that $`\dot{m}_{}m_{}/t_{\mathrm{ff}}`$ to within a factor 1.5 for spherical cores in which the effective sound speed increases outward.

We can express the accretion rate in terms of the mean pressure of the clump in which the star forms, $`\overline{P}`$ (eq. 1), and the final mass of the star, $`m_f`$,

$$\dot{m}_{}=4.75\times 10^4ϵ_{\mathrm{core}}^{1/4}(f_{\mathrm{gas}}\varphi _P\alpha _{\mathrm{vir}})^{3/8}\left(\frac{m_f}{30\mathrm{M}_{}}\right)^{3/4}\mathrm{\Sigma }^{3/4}\left(\frac{m_{}}{m_f}\right)^j\mathrm{M}_{}\mathrm{yr}^1,$$
$`(5)`$
where $`\varphi _PP_s/\overline{P}`$ is the ratio of the core’s surface pressure to the mean pressure in the clump and $`j3(2k_\rho )/[2(3k_\rho )]`$; for $`k_\rho =3/2`$, $`j=1/2`$. This typical accretion rate is large because the core is embedded in a high-pressure environment, so that in hydrostatic equilibrium it has a high density and short free-fall time. Because massive cores are turbulent and clumpy, we expect the accretion rate to exhibit large fluctuations. Whereas clumpiness is an attribute of massive star-forming regions in our model, it is not a prerequisite, as suggested in ref. 3. The star formation time is:

$$t_f=1.26\times 10^5ϵ_{\mathrm{core}}^{1/4}(f_{\mathrm{gas}}\varphi _P\alpha _{\mathrm{vir}})^{3/8}\left(\frac{m_f}{30\mathrm{M}_{}}\right)^{1/4}\mathrm{\Sigma }^{3/4}\mathrm{yr}.$$
$`(6)`$
The weak dependence of the star formation time on the final stellar mass means that low-mass and high-mass stars can form approximately at the same time. Furthermore, $`t_f`$ is small compared to the estimated timescale $`10^6\mathrm{yr}`$ for cluster formation20. With the fiducial values of the parameters, a $`100\mathrm{M}_{}`$ star forming in a clump with $`\mathrm{\Sigma }=1\mathrm{g}\mathrm{cm}^2`$ has a final accretion rate of $`1.1\times 10^3\mathrm{M}_{}\mathrm{yr}^1`$ and a star formation time of $`1.8\times 10^5\mathrm{yr}`$.

A necessary condition for massive star formation is that the ram pressure associated with accretion exceed the radiation pressure at the point where the dust in the infalling gas is destroyed21. For a 100 $`M_{}`$ star, this requires18 $`\dot{m}_{}>6\times 10^4M_{}`$ yr-1, which is satisfied for $`\mathrm{\Sigma }>0.5`$ g cm-2. Once the core has formed a star, the star can continue to grow by Bondi-Hoyle accretion22 provided $`m_{}<10\mathrm{M}_{}`$ so that radiation pressure does not prevent the focussing of gas streamlines in the wake of the star. The Bondi-Hoyle accretion rate $`\dot{m}_{}=2.2\times 10^7(M_{\mathrm{cl}}/4\times 10^3\mathrm{M}_{})^{5/4}\mathrm{\Sigma }^{3/4}(m_{}/10\mathrm{M}_{})^2\mathrm{M}_{}\mathrm{yr}^1`$ (CFM & JCT, in prep.) is so low that this process does not significantly alter the stellar mass under the conditions observed in the clumps of ref. 7.

Direct comparison of our results with observation is difficult because the actual masses of massive protostars are poorly determined. Our approach is to predict the properties of some well-studied massive protostars in terms of their bolometric luminosities. The bolometric luminosity $`L_{\mathrm{bol}}`$ has contributions from main-sequence nuclear burning $`L_{\mathrm{ms}}`$, deuterium burning $`L_D`$, and accretion $`L_{\mathrm{acc}}`$. The accretion luminosity $`L_{\mathrm{acc}}=f_{\mathrm{acc}}Gm_{}\dot{m}_{}/r_{}`$, where $`f_{\mathrm{acc}}`$ is a factor of order unity accounting for energy radiated by an accretion disk, advected into the star or converted into kinetic energy of outflows, and where the stellar radius $`r_{}`$ may depend sensitively on the accretion rate $`\dot{m}_{}`$. Massive stars join the main sequence during their accretion phase at a mass that also depends on the accretion rate23. To treat accelerating accretion rates, we have developed a simple model for protostellar evolution based on that of Nakano et al.6,24 The model accounts for the total energy of the protostar as it accretes and dissociates matter and, if the central temperature $`T_c>10^6\mathrm{K}`$, burns deuterium. We have modified this model to include additional processes, such as deuterium shell burning, and we have calibrated these modifications against the more detailed calculations of Palla & Stahler25,23.

Our model allows us to make predictions for the masses and accretion rates of embedded protostars that are thought to power hot molecular cores (CFM & JCT, in prep.). Figure 2 compares our theoretical tracks with the observed bolometric luminosities of several sources. We find that uncertainties in the value of the pressure create only small uncertainties in $`m_{}`$ for $`L_{\mathrm{bol}}>\mathrm{few}\times 10^4\mathrm{L}_{}`$.

The infrared and submillimeter spectra of accreting protostars and their surrounding envelopes have been modelled in ref. 5., modelling the same sources shown in Fig. 2. We note that uncertainties in the structure of the gas envelope and the possible contributions from additional surrounding gas cores or diffuse gas will affect the observed spectrum. Comparing results, our inferred stellar masses are similar, but our accretion rates are systematically smaller by factors of $`25`$. The modelled5 high accretion rates of $`10^3\mathrm{M}_{}\mathrm{yr}^1`$ for stars with $`m_{}10\mathrm{M}_{}`$ would be difficult to achieve unless the pressure was increased substantially; for example, if the stars are destined to reach $`m_f30\mathrm{M}_{}`$, pressure increases of a factor $`40`$ are required.

Received 17 September 2001; accepted 2 January 2002

References

1. Bernasconi, P. A. & Maeder, A. About the absence of a proper zero age main sequence for massive stars. Astron. Astrophys. 307, 829-839 (1996).

2. McLaughlin, D. E. & Pudritz, R. E. Gravitational Collapse and Star Formation in Logotropic and Nonisothermal Spheres. Astrophys. J. 476, 750-765 (1997).

3. Stahler, S. W., Palla, F. & Ho, P. T. P. The Formation of Massive Stars. Protostars & Planets IV, eds. Mannings, V., Boss, A. P. & Russell, S. S. (University of Arizona Press, Tucson), 327-351 (2000).

4. Behrend, R. & Maeder, A. Formation of massive stars by growing accretion rate. Astron. Astrophys. 373, 190-198 (2001).

5. Osorio, M., Lizano, S. & D’Alessio, P. Hot Molecular Cores and the Formation of Massive Stars. Astrophys. J. 525, 808-820 (1999).

6. Nakano, T., Hasegawa, T., Morino, J.-I. & Yamashita, T. Evolution of Protostars Accreting Mass at Very High Rates: Is Orion IRc2 a Huge Protostar? Astrophys. J. 534, 976-983 (2000).

7. Plume, R., Jaffe, D. T., Evans, N. J. (II), Martin-Pintado, J. & Gomez-Gonzalez, J. Dense Gas and Star Formation: Characteristics of Cloud Cores Associated with Water Masers. it Astrophys. J. 476, 730-749 (1997).

8. Bertoldi, F. & McKee, C. F. Pressure-confined clumps in magnetized molecular clouds. Astrophys. J. 395, 140-157 (1992).

9. Hillenbrand, L. A. & Hartmann, L. W. A Preliminary Study of the Orion Nebula Cluster Structure and Dynamics. Astrophys. J. 492, 540-553 (1998).

10. van der Tak, F. F. S., van Dishoeck, E. F., Evans, N. J. (II) & Blake, G. A. Structure and Evolution of the Envelopes of Deeply Embedded Massive Young Stars. Astrophys. J. 537, 283-303 (2000).

11. Vázquez-Semadeni, E., Ostriker, E. C., Passot, T., Gammie, C. F. & Stone, J. M. Compressible MHD Turbulence: Implications for Molecular Cloud and Star Formation. Protostars & Planets IV, eds. Mannings, V., Boss, A. P. & Russell, S. S. (University of Arizona Press, Tucson), 3-28 (2000).

12. Richer, J. S., Shepherd, D. S., Cabrit, S., Bachiller, R., & Churchwell, E. Molecular Outflows from Young Stellar Objects. Protostars & Planets IV, eds. Mannings, V., Boss, A. P. & Russell, S. S. (University of Arizona Press, Tucson), 867-894 (2000).

13. Matzner, C. D. & McKee, C. F. Efficiencies of Low-Mass Star and Star Cluster Formation. Astrophys. J. 545, 364-378 (2000).

14. McLaughlin, D. E. & Pudritz, R. E. A Model for the Internal Structure of Molecular Cloud Cores. Astrophys. J. 469, 194-208 (1996).

15. Larson, R. B. Turbulence and star formation in molecular clouds. Mon. Not. R. Astron. Soc. 194, 809-826 (1981).

16. Stahler, S. W., Shu, F. H. & Taam, R. E. The evolution of protostars. I - Global formulation and results. Astrophys. J. 241, 637-654 (1980).

17. Crutcher, R. M. Magnetic Fields in Molecular Clouds: Observations Confront Theory. Astrophys. J. 520, 706-713 (1999).

18. Jijina, J. & Adams, F. C. Infall Collapse Solutions in the Inner Limit: Radiation Pressure and Its Effects on Star Formation. Astrophys. J. 462, 874-887 (1996).

19. Shu, F. H. Self-similar collapse of isothermal spheres and star formation. Astrophys. J. 214, 488-497 (1977).

20. Palla, F. & Stahler, S. W. Star Formation in the Orion Nebula Cluster. Astrophys. J. 525, 772-783 (1999).

21. Wolfire, M. G. & Cassinelli, J. Conditions for the formation of massive stars. Astrophys. J. 319, 850-867 (1987).

22. Bonnell, I. A., Bate, M. R., Clarke, C. J. & Pringle, J. E. Competitive accretion in embedded stellar clusters. Mon. Not. R. Astron. Soc. 323, 785-794 (2001).

23. Palla, F. & Stahler, S. W. The evolution of intermediate-mass protostars. II - Influence of the accretion flow. Astrophys. J. 392, 667-677 (1992).

24. Nakano, T., Hasegawa, T. & Norman, C. The Mass of a Star Formed in a Cloud Core: Theory and Its Application to the Orion A Cloud. Astrophys. J. 450, 183-195 (1995).

25. Palla, F. & Stahler, S. W. The evolution of intermediate-mass protostars. I - Basic results. Astrophys. J. 375, 288-299 (1991).

26. Schaller, G., Schaerer, D., Meynet, G. & Maeder, A. New grids of stellar models from 0.8 to 120 solar masses at Z = 0.020 and Z = 0.001. Astron. Astrophys. Supp. 96, 269-331 (1992).

27. Hunter, T. R. et al. G34.24+0.13MM: A Deeply Embedded Proto-B-Star. Astrophys. J. 493, L97-L100 (1998).

28. Molinari, S., Testi, L., Brand, J., Cesaroni, R. & Palla, F. IRAS 23385+6053: A Prototype Massive Class 0 Object. Astrophys. J. 505, L39-L42 (1998).

29. Kaufman, M. J., Hollenbach, D. J. & Tielens, A. G. G. M. High- Temperature Molecular Cores near Massive Stars and Application to the Orion Hot Core. Astrophys. J. 497, 276-287 (1998).

30. Wyrowski, F., Schilke, P., Walmsley, C. M. & Menten, K. M. Hot Gas and Dust in a Protostellar Cluster near W3(OH). Astrophys. J. 514, L43-L46 (1999).

ACKNOWLEDGEMENTS. We thank Steve Stahler, Ralph Pudritz, Malcolm Walmsley and Mark Krumholz for discussions. This work was supported by the NSF, by NASA (which supports the Center for Star Formation Studies) and (for J.C.T.) by a Spitzer-Cotsen fellowship from Princeton University.

Correspondence should be addressed to CFM (e-mail: cmckee@mckee.berkeley.edu) or to JCT (e-mail: jt@astro.princeton.edu).

Fig. 1. Variation of model parameters and results with $`k_\rho `$. Let $`P_s=\varphi _P\overline{P}`$ be the surface pressure of a core. The properties at the surface of the core are given by $`r_s=[AGm_f^2/(2\pi k_P^2ϵ_{\mathrm{core}}^2P_s)]^{1/4}`$, $`\rho _s=[Ak_P^2ϵ_{\mathrm{core}}^2P_s^3/(2\pi G^3m_f^2)]^{1/4}`$, and $`c_s=[2\pi G^3m_f^2P_s/(Ak_P^2ϵ_{\mathrm{core}}^2)]^{1/8}`$. We anticipate that the overall star formation efficiency in the clump will be relatively high, so in equation (1) we adopt $`f_{\mathrm{gas}}=2/3`$ as a fiducial value. We estimate $`\varphi _P=P_s/\overline{P}2`$ (C.F.M. & J.C.T., in prep.), and we set $`\alpha _{\mathrm{vir}}`$ equal to unity7. We take $`k_\rho =1.5`$ as a typical value, corresponding to $`\gamma _p=2/3`$, $`k_P=1`$ and $`A=3/4`$. Following equation (3), we evaluate $`\varphi _{}`$ using the results of McLaughlin & Pudritz2, $`\varphi _{}=\pi \sqrt{3}[(2\gamma _p)^2(43\gamma _p)^{(76\gamma _p)/2}8^{(3\gamma _p5)/2}m_0]^{1/(43\gamma _p)}`$, where $`m_0`$ is tabulated in their paper. Over the entire range of $`\gamma _p`$ and $`k_\rho `$ relevant to molecular clouds ($`0\gamma _p1`$, $`1k_\rho 2`$), $`\varphi _{}1.620.48k_\rho `$ to within about 1%. The star formation time decreases from $`3.5t_{\mathrm{ff}}`$ to $`1.5t_{\mathrm{ff}}`$ as $`\gamma _p`$ varies from 0 to 1. The variation of $`t_f`$ and $`\dot{m}_f`$ relative to the $`k_\rho =1.5`$ case is also shown. Note that the singular polytropic model in hydrostatic equilibrium breaks down for $`k_\rho =1`$, $`\gamma _p=0`$, since then the pressure gradient vanishes ($`k_P=0`$).

Fig. 2. Derived properties of nearby massive protostars. Solid lines show the predicted evolution in luminosity of protostars (including their accretion disks, but allowing for the powering of protostellar outflows so that $`f_{\mathrm{acc}}=0.5`$; J.C.T. & C.F.M., in prep.) of final mass 7.5, 15, 30, 60 and 120 $`\mathrm{M}_{}`$ accreting from cores with $`k_\rho =1.5`$ embedded in a $`\mathrm{\Sigma }=1\mathrm{g}\mathrm{cm}^2`$ clump, typical of Galactic regions7. The luminosity step occurring at around 5 to 7 $`\mathrm{M}_{}`$, depending on the model, corresponds to the onset of deuterium shell burning, which swells the protostellar radius by a factor of about two and thus reduces the accretion luminosity by the same factor. The dashed and long dashed lines show a 30 $`\mathrm{M}_{}`$ star forming in a clump with mean pressure ten times smaller and larger than the fiducial value, respectively. The dotted line shows the zero age main sequence luminosity from ref. 26. Four observed hot molecular cores are shown: G34.24+0.13MM27, IRAS 23385+605328, Orion Hot Core29 and W3($`\mathrm{H}_2\mathrm{O}`$)—the Turner-Welch object30. The vertical error bar illustrates the uncertainty in the bolometric luminosity. The horizontal error bar shows the corresponding range of allowed values of $`m_{}`$ for the $`\mathrm{\Sigma }=1\mathrm{g}\mathrm{cm}^2`$ models. These values and the constraints on $`\dot{m}_{}`$ are listed here for each source.


text:
# Constraints on tan𝛽 in SUSY SU(5) GUT

## I Introduction

* \- What are unified ?
As one of the general features of the SU(5) GUT models, both Yukawa coupling constants of the bottom-quark and the tau-lepton are unified at GUT scale. Additionally, in several SU(5) models including the SUSY SU(5) GUT modelNilles , the evidence for the unification of gauge coupling constants is numerically providedAmaldi  .
* \- What is $`M_{SUSY}`$ ?
In the conventional analysisAmaldi , the gauge unification in the SUSY SU(5) GUT model is achieved by introducing an energy-scale constant $`M_{SUSY}`$ by-hand to decouple the lower energy-scale region from the effects of the supersymmetric particles. Thus, $`M_{SUSY}`$ is regarded as the threshold between the standard model(SM) and SUSY. The value of $`M_{SUSY}`$ depends on the strong coupling constant $`\alpha _3`$ very sensitivelyAmaldi . For instance, while $`M_{SUSY}`$ 1(TeV) for $`\alpha _3=0.106`$ is givenAmaldi  in 1991, $`M_{SUSY}`$ 100(GeV) is obtained for the recent valueParticle  of $`\alpha _3=0.1185`$. If $`M_{SUSY}`$ were really less than 100(GeV), we should have already detected such supersymmetric particles in the region of $``$ 100(GeV) experimentally, however we have not yet. This fact may allow us to think of $`M_{SUSY}`$ as a negligible parameter, or, at least, the authors should discuss on the possibility of the models without $`M_{SUSY}`$. Note that such model without $`M_{SUSY}`$ means the pure supersymmetric theory without including the conventional standard model.
* How much $`\mathrm{tan}\beta `$ is optimal ?
$`\mathrm{tan}\beta `$ is evaluated by the two different ways of the likelihood analysis at two-loop level in this paper. One of the two ways is on the unification of the gauge coupling constants, and another is of the bottom-tau’s Yukawa coupling constants. The sequence to optimize $`\mathrm{tan}\beta `$ is as follows:
In the first, $`M_{SUSY}`$ and the GUT scale $`M_X`$ are fixed by the $`\chi ^2`$-fitting on the unification of the gauge$`(\chi _g^2)`$ coupling constants. Once $`M_{SUSY}`$ and $`M_X`$ are fixed, corresponding $`\mathrm{tan}\beta `$ is derived numerically to realize the bottom-tau unification at $`M_X`$ scale.
Additionally, the analysis on the model without $`M_{SUSY}`$ is also carried out in the same manner.

## II Analysis of $`\mathrm{tan}\beta `$ with $`M_{SUSY}`$

On the optimization of $`\mathrm{tan}\beta `$ is discussedBarger  by the renormalization group equations(RGEs) of the gauge and Yukawa coupling constants at two-loop levelArason ; Castano  in the SM and SUSY. As the Yukawa coupling constants, heavier three of the ordinal Fermions (top-quark, bottom-quark, and tau-lepton) are taken into consideration.
First of all, the explicit values of the parameters in this paper are referred from the latest issue of the Particle Data GroupParticle .
On the boundary conditions of the RGEs, they are described further below. The conditions of the gauge coupling constants are given on Z-boson’s mass shell $`M_Z=91.184`$(GeV) in the $`\overline{\text{MS}}`$-scheme as:

$`\alpha _1`$ $`=`$ $`{\displaystyle \frac{5}{3}}{\displaystyle \frac{\alpha _{\overline{MS}}}{1\mathrm{sin}^2\theta _{\overline{MS}}}},`$ (1)
$`\alpha _2`$ $`=`$ $`{\displaystyle \frac{\alpha _{\overline{MS}}}{\mathrm{sin}^2\theta _{\overline{MS}}}},`$ (2)
as their explicit values are:

$`\alpha _{\overline{MS}}^1`$ $`=`$ $`127.934\pm 0.027,`$ (3)
$`\mathrm{sin}^2\theta _{\overline{MS}}`$ $`=`$ $`0.23117\pm 0.00016,`$ (4)
$`\alpha _3`$ $`=`$ $`0.1185\pm 0.0020.`$ (5)
Moreover, the definition of $`\mathrm{tan}\beta `$ at $`M_Z`$ is given as usual:

$$\mathrm{tan}\beta (M_Z)v_2(M_Z)/v_1(M_Z),$$
(6)
where, $`v_1`$ and $`v_2`$ imply the vacuum expectation values(VEVs) of two Higgs doublets, and as their explicit value:

$$v(M_Z)=(v_1^2(M_Z)+v_2^2(M_Z))^{1/2}=246\text{(GeV)},$$
(7)
is fixed. For simplicity, mere $`\mathrm{tan}\beta `$ means $`\mathrm{tan}\beta (M_Z)`$ in this paper.
On the Higgs quartic coupling constant $`\lambda `$, this is precisely defined at $`M_{SUSY}`$Barger  as:

$$\lambda (M_{SUSY})=\pi \left[\frac{3}{5}\alpha _1(M_{SUSY})+\alpha _2(M_{SUSY})\right]\mathrm{cos}^22\beta (M_{SUSY}),$$
(8)
however, it is very difficult to derive $`\lambda `$ from this definition. As $`\lambda `$, the authors make use of an approximate, and easier alternative definition at $`M_Z`$:

$$\lambda \pi \left[\frac{3}{5}\alpha _1(M_Z)+\alpha _2(M_Z)\right]\mathrm{cos}^22\beta (M_Z).$$
(9)
Note that the numerical difference between the definitions Eq. (8) and Eq. (9) is very small as about $`0.1\%`$, and $`\lambda `$ is decoupled to the Yukawa coupling constants among the RGEs at one-loop level, therefore, all the numerical results in this paper are effectually independent of the choice of the definition of $`\lambda `$.
Boundary conditions of the Yukawa coupling constants are given at each on-shell energy-scales as:

$$\alpha _q(m_q)=\frac{m_q^2(m_q)}{2\pi v^2(m_q)},$$
(10)
where, the suffix $`q`$ implies the flavor of Fermions as $`q=t`$, $`b`$, or $`\tau `$, respectively. Their explicit values are as follows:

$`m_t(m_t)`$ $`=`$ $`174.3\pm 5.1\text{(GeV)},`$ (11)
$`m_b(m_b)`$ $`=`$ $`4.2\pm 0.2\text{(GeV)},`$ (12)
$`m_\tau (m_\tau )`$ $`=`$ $`1.77699\pm 0.00029\text{(GeV)}.`$ (13)
When assuming the existence of $`M_{SUSY}`$, corresponding boundary conditions are given as:

$`\alpha _i^{SUSY}(M_{SUSY})`$ $`=`$ $`\alpha _i^{SM}(M_{SUSY}),(i=1,2,3),`$ (14)
$`\alpha _t^{SUSY}(M_{SUSY})`$ $`=`$ $`{\displaystyle \frac{m_t^2(M_{SUSY})}{2\pi v^2(M_{SUSY})\mathrm{sin}^2\beta (M_{SUSY})}}`$
$`=`$ $`\alpha _t^{SM}(M_{SUSY}){\displaystyle \frac{1}{\mathrm{sin}^2\beta (M_{SUSY})}}`$
$`=`$ $`\alpha _t^{SM}(M_{SUSY}),`$
$`\alpha _{b,\tau }^{SUSY}(M_{SUSY})`$ $`=`$ $`\alpha _{b,\tau }^{SM}(M_{SUSY}){\displaystyle \frac{1}{\mathrm{cos}^2\beta (M_{SUSY})}}`$
$`=`$ $`\alpha _{b,\tau }^{SM}(M_{SUSY}),`$
where, the variables with a prime$`(^{})`$ mean the redefinition of them in the SM region to keep their continuity between the regions of the SM and SUSY. With such redefinitions, the graphs in this paper is continuous even if crossing on the $`M_{SUSY}`$ border.
The conditions at GUT energy-scale $`M_X`$ are described in the following subsections.

### II.1 On the gauge coupling constants

The average value $`\alpha _g`$ among $`\alpha _1,\alpha _2`$, and $`\alpha _3`$ is introduced around the GUT scale $`M_X`$. $`\chi _g^2`$-fitting function is defined as the squared sum of each differences between $`\alpha _g`$ and $`\alpha _i(i=13)`$ as follows:

$$\chi _g^2\underset{i=1}{\overset{3}{}}\left(\frac{\alpha _i\alpha _g}{\delta \alpha _i}\right)^2,$$
(17)
where, $`\delta \alpha _i`$ means the uncertainty width of $`\alpha _i`$, respectively. the GUT scale $`M_X`$ is defined as the energy-scale where $`\chi _g^2`$ becomes the minimum. In order to make $`\chi _g^2`$ the minimum, $`M_{SUSY}`$ is fine-tuned. As the results, we find $`\chi _g^20`$ can be achieved by tuning $`M_{SUSY}`$ at any $`\mathrm{tan}\beta `$ in the ranges of the errorbars of the other parameters. Fine-tuned values of $`M_{SUSY}`$ are shown as functions of $`\mathrm{tan}\beta `$ with the uncertainty width of $`\alpha _3`$ in Fig. 1 . We find the derived $`M_{SUSY}`$’s value depends on only the strong coupling constant $`\alpha _3`$ rather than the other parameters.

### II.2 On the Yukawa coupling constants

When the SU(5) GUT scenario is presumed, Yukawa coupling constants of the bottom-quark and tau-lepton are unified from each other at $`M_X`$. Their $`\chi _Y^2`$-fitting function is defined as the following equations:

$$\chi _Y^2\underset{q=b,\tau }{}\left(\frac{\alpha _q\alpha _Y}{\delta \alpha _q}\right)^2,$$
(18)
$$\alpha _Y=\left(\alpha _b+\alpha _\tau \right)/2,$$
(19)
where, $`q`$ implies $`b`$ or $`\tau `$, and $`\delta \alpha _q`$ means the uncertainty width of $`\alpha _q`$, respectively. With the uncertainty width $`\delta =0.2`$(GeV) of $`m_b`$, $`\chi _Y^2`$ at $`M_X`$ scale as functions of $`\mathrm{tan}\beta `$ are shown in Fig. 2. The behavior of $`\chi _Y^2`$ is not sensitive on the variance of $`m_t`$ or $`\alpha _3`$. Also shown in Fig. 2, only small $`\mathrm{tan}\beta 2`$ is allowed by the likelihood analysis with the $`\chi _Y^2`$ function, and this minimal point of $`\chi _Y^2`$ is stable despite of the uncertainty of the parameters like $`m_b`$, $`m_t`$, or $`\alpha _3`$.

## III Analysis of $`\mathrm{tan}\beta `$ without $`M_{SUSY}`$

As mentioned above, on the model without $`M_{SUSY}`$ is discussed in this section. Therefore, the RGEsCastano  in this section are pure SUSY. On the boundary conditions of the gauge coupling constants, the VEVs of Higgs bosons, and $`\mathrm{tan}\beta `$, they are given at $`M_Z`$ in the same way of the model with including $`M_{SUSY}`$ as described in the previous section as from Eq. (3) to Eq. (6). The boundary conditions of the Yukawa coupling constants are given at the energy-scales of each on-shell masses as follows:

$$\alpha _t(m_t)=\frac{m_t^2(m_t)}{2\pi v^2(m_t)\mathrm{sin}^2\beta (m_t)},$$
(20)
and

$$\alpha _{b,\tau }(m_{b,\tau })=\frac{m_{b,\tau }^2(m_{b,\tau })}{2\pi v^2(m_{b,\tau })\mathrm{cos}^2\beta (m_{b,\tau })}.$$
(21)
The definitions of $`\chi _g^2`$ and $`\chi _Y^2`$ at $`M_X`$ are the same ones Eq. (17) and Eq. (18) in the previous section, however, in this section, $`\chi _g^2`$ depends on $`\mathrm{tan}\beta `$ or other parameters directly in contrast to the previous section, because the decoupling boundary $`M_{SUSY}`$ is vanished away from the model.
The dependence of $`\mathrm{tan}\beta `$ on $`\chi _g^2`$ with the uncertainty width of $`\alpha _3`$ is shown in Fig. 3 . As Fig. 3 shows, two minimal points exist at small and large $`\mathrm{tan}\beta `$. The variance of $`\chi _Y^2`$ at $`M_X`$ is shown with each uncertainty width of $`m_b`$, $`m_\tau `$, or $`\alpha _3`$ in Fig. 4a, 4b, or 4c, respectively. As shown in these graphs, usually $`\chi _Y^2`$ has three local minima.

## IV Discussions

### IV.1 On the model with $`M_{SUSY}`$

Since $`\chi ^2`$ of the gauge coupling constants($`\chi _g^2`$) can be always equal to zero by fine-tuning of $`M_{SUSY}`$, there exists no effective constraint on $`\mathrm{tan}\beta `$ from the unification of the gauge coupling constants. $`\mathrm{tan}\beta `$ is evaluated as 1.6 by the criteria of the bottom-tau unification in this model with $`M_{SUSY}`$. $`M_{SUSY}`$ is estimated about 2(TeV) to unify the gauge coupling constants at two-loop level. (As a reference, $`M_{SUSY}M_Z`$ is derived at one-loop level. ) Nevertheless, such small $`\mathrm{tan}\beta 2`$ is difficult to be allowed by $`g2`$ experimentsChattopadhyay , because small $`\mathrm{tan}\beta `$ makes the masses of sparticles quite smaller than expected. The evolutions of the gauge and Yukawa coupling constants at $`\mathrm{tan}\beta =1.6`$ are shown in Fig. 6a and 6b, respectively. The summary of the arguments in this subsection is shown in the Table 1.

### IV.2 On the model without $`M_{SUSY}`$

The value of $`\mathrm{tan}\beta `$ is possibly settled at 2, $`625`$, or 50 in accordance with three minima of $`\chi _Y^2`$. However, large $`\mathrm{tan}\beta 50`$ is denied theoretically, because $`\chi _g^2`$ and $`\chi _Y^2`$ cannot be minimal simultaneously, i.e., they are exclusive from each other as shown in Fig. 5. On the small $`\mathrm{tan}\beta 2`$ choice, both of $`\chi _g^2`$ and $`\chi _Y^2`$ are enough small, however, this is excluded by the recent $`g2`$ experimentsChattopadhyay  as mentioned previously. Large $`\mathrm{tan}\beta 50`$ is also not probable if $`bs\gamma `$ experimentb-sr  is true. Moreover, the authors are dubious whether these two local minima of $`\chi _Y^2`$ at small and large $`\mathrm{tan}\beta `$ are physical ones or not, because they make the RGEs divergent at one-loop level, therefore, the evolution of the Yukawa coupling constants seems strange even at two-loop level, for instance, as shown in Fig. 7 at $`\mathrm{tan}\beta =1.9`$ . Medium $`\mathrm{tan}\beta `$ as $`625`$ is good for $`\chi _Y^2`$ of the bottom-tau unification, however, $`\chi _g^2`$ of the gauge unification cannot be enough small at two-loop level. Thus, precise unification of the gauge coupling constants cannot occur with medium $`\mathrm{tan}\beta `$ as $`625`$. The evolutions of the gauge or Yukawa coupling constants at $`\mathrm{tan}\beta =6.0`$ are shown in Fig. 8a or 8b, respectively. These discussions are summarized in the Table 2.

## V Summary and Conclusion

In the SU(5) GUT scenario, $`\mathrm{tan}\beta `$ is evaluated by making use of the RGEs of the gauge and Yukawa coupling constants. Only small $`\mathrm{tan}\beta 2`$ is allowed in the model including $`M_{SUSY}`$, and $`\mathrm{tan}\beta `$ is estimated as 2 or 6 to 25 in the model without $`M_{SUSY}`$, theoretically. However, small $`\mathrm{tan}\beta 2`$ is difficult to be allowed by the analysis of the recent $`g2`$ experimentChattopadhyay . As the result, only medium $`\mathrm{tan}\beta `$ as $`625`$ is acceptable in the SUSY SU(5) GUT model without $`M_{SUSY}`$ if the bottom-tau unification is prior to the unification among the gauge coupling constants.

###### Acknowledgements.

The authors thank Y. Okada, M. Kakizaki and T. Nihei for useful advice. We are grateful to D. W. Hertzog for helpful comments on our preliminary stage.