text:
# 1 Introduction

## 1 Introduction

The observational data do not severely restrict lepton asymmetry of the Universe which, if present, resides almost entirely in the neutrino sector. In particular, the chemical potentials and the temperature of relic neutrinos can be comparable in magnitude. The nucleosynthesis bounds on the neutrino degeneracy parameters, $`\xi _{\nu _f}=\mu _{\nu _f}/T_\nu `$, are $`0.06<\xi _{\nu _e}<1.1`$, $`|\xi _{\nu _{\mu ,\tau }}|<6.9`$ . Current measurements of the cosmic microwave background anisotropy put an upper limit of $`35`$ on degeneracy parameters for all neutrino species .

The large lepton asymmetry could have profound consequences in the early history of the Universe, when the temperature was of order of the electroweak scale. The lepton charge alters the pattern of electroweak symmetry breaking, because the lepton number density forces the Higgs field to Bose condense and therefore can impede restoration of electroweak symmetry at high temperature . If the lepton number density is sufficiently large, the Universe has always been in the symmetry-broken phase and there was no electroweak phase transition . This observation provides a possible resolution of the monopole and the domain wall problems .

Although it seems natural to assume that the lepton and the baryon asymmetries of the Universe are of the same order of magnitude, mechanisms that could generate primordial lepton asymmetry of order one without generating large baryon asymmetry are known , and degeneracy of relic neutrinos remains at least a logical possibility. The assumption of large neutrino degeneracy was invoked, in particular, to explain an origin of the ultra-high energy cosmic rays .

In this paper I will study the consequences of the large lepton asymmetry at temperatures above the quark-hadron transition but below the electroweak scale. It appears that, when quarks are deconfined, the lepton asymmetry indirectly affects hadronic (strongly interacting) sector in such a way that all light quark species acquire large chemical potentials. A deviation of the lepton number density from zero and an overall electric neutrality cause an asymmetric distribution of the electric charge between the lepton and the hadron sectors. As a consequence, the quark-gluon plasma is electrically charged when lepton asymmetry is of order one.

## 2 Chemical potentials

There are five charges that are conserved or almost conserved when the temperature drops below the electroweak scale: electric charge, baryon number and lepton numbers for each of the three generations of fermions. Accordingly, one can introduce five chemical potentials to describe a state with non-zero densities of these charges. On average, two of the conserved charge densities are zero, because the Universe must be electrically neural and the baryon asymmetry of the Universe is negligibly small. But densities of lepton numbers can potentially be large enough to require chemical potentials comparable to the temperature.

The chemical potential for each of the particle species is a linear combination of the chemical potentials for conserved charges:

$$\mu _i=\underset{\alpha }{}q_i^\alpha \mu _\alpha .$$
(1)
Here $`q_i^\alpha =(Q_i,B_i,L_i^e,L_i^\mu ,L_i^\tau )`$ are charges of $`i`$th species. In assumption of large lepton asymmetry, chemical potentials for lepton numbers are non-zero at present and cause asymmetries in equilibrium distributions of the relic neutrinos. When the temperature was higher than the mass of the electron, the lepton number chemical potential alone would lead to an asymmetric distribution of electrons and would induce a non-zero density of electric charge. The electric neutrality thus requires to set $`\mu _Q=\mu _{L^e}`$, which makes the chemical potential for electrons equal to zero.

The situation drastically changes when quarks are deconfined. Light quarks can be considered essentially massless at temperatures above the quark-hadron transition and, since they carry both the electric charge and the baryon number, non-zero $`\mu _Q`$ would induce large baryon asymmetry unless $`\mu _B`$ is also non-zero. The conditions of electric neutrality and baryon symmetry then become more complicated and all chemical potentials generically become of the same order of magnitude.

I will consider temperatures of order of few Gev, when the temperature suppresses strong interactions due to the asymptotic freedom, and the ideal gas approximation for the quark-gluon plasma is more or less accurate. The active species at such temperatures are neutrinos, electrons, muons, photons, gluons and u, d and s quarks. Weak processes at the QCD epoch are rapid enough to maintain thermal equilibrium, and I will also assume that strong interactions do not considerably distort equilibrium distributions and will neglect masses of all active particle species.

The excess of particles over anti-particles per unit volume for two-component non-interacting fermions with chemical potential $`\mu `$ at temperature $`T`$ is

$$n_+n_{}=\frac{T^3}{6}\left(\kappa +\frac{\kappa ^3}{\pi ^2}\right).$$
(2)
The degeneracy parameters at the QCD epoch will be denoted by $`\kappa `$: $`\kappa =\mu /T`$, to distinguish them from the degeneracy parameters now, which are denoted by $`\xi `$. With the help of this equality, the electric neutrality, absence of considerable baryon number and requirement that there is a given lepton asymmetry yield a set of five equations for five chemical potentials:

$`{\displaystyle \frac{T^3}{6}}{\displaystyle \underset{i}{}}Q_i\left(\kappa _i+{\displaystyle \frac{\kappa _i^3}{\pi ^2}}\right)=`$ $`0`$
$`{\displaystyle \frac{T^3}{6}}{\displaystyle \underset{i}{}}B_i\left(\kappa _i+{\displaystyle \frac{\kappa _i^3}{\pi ^2}}\right)`$ $`=`$ $`0`$
$`{\displaystyle \frac{T^3}{6}}{\displaystyle \underset{i}{}}L_i^f\left(\kappa _i+{\displaystyle \frac{\kappa _i^3}{\pi ^2}}\right)`$ $`=`$ $`n_{L^f}n_{\overline{L}^f},`$ (3)
where $`n_{L^f}`$ ($`n_{\overline{L}^f}`$) are number densities of leptons (anti-leptons) of flavor $`f`$.

Though the lepton numbers are not conserved exactly, and some lepton charge could be produced by neutrino oscillations during nucleosynthesis , the asymmetry generated in this way is generically much smaller than one . So, if large density of lepton number was present at the QCD epoch, later it was just diluted by the expansion of the Universe. Since the entropy density,

$$s=\frac{P+\rho \mu (n_+n_{})}{T},$$
(4)
scales in the same way***This assumption implies an entropy conservation and, in particular, an absence of considerable entropy production at the QCD phase transition., the ratio of the lepton number density to the entropy is time-independent:

$$n_{L^f}n_{\overline{L}^f}=s\left(\frac{n_{L^f}n_{\overline{L}^f}}{s}\right)_{\mathrm{now}}.$$
(5)
Taking into account that the entropy of a gas of non-interacting two-component fermions is

$$s_{\mathrm{ferm}}=T^3\left(\frac{7\pi ^2}{90}+\frac{\kappa ^2}{6}\right),$$
(6)
and that each bosonic degree of freedom contributes $`2\pi ^2T^3/45`$, we find:

$$s=T^3\left(\frac{247\pi ^2}{90}+\underset{i}{}\frac{\kappa _i^2}{6}\right)$$
(7)
at the QCD epoch. The entropy at present is carried by photons and by neutrinos:

$$s_{\mathrm{now}}=\frac{4\pi ^2}{45}T_\gamma ^3+T_\nu ^3\underset{f}{}\left(\frac{7\pi ^2}{90}+\frac{\xi _{\nu _f}^2}{6}\right)=T_\nu ^3\left(\frac{43\pi ^2}{90}+\underset{f}{}\frac{\xi _{\nu _f}^2}{6}\right),$$
(8)
where the relation $`T_\gamma ^3=11T_\nu ^3/4`$ was used in the last equality. The present value of the lepton asymmetry depends on the neutrino chemical potentials as follows:

$$\left(n_{L^f}n_{\overline{L}^f}\right)_{\mathrm{now}}=\frac{T_\nu ^3}{6}\left(\xi _{\nu _f}+\frac{\xi _{\nu _f}^3}{\pi ^2}\right).$$
(9)
Finally, we get:

$$n_{L^f}n_{\overline{L}^f}=\frac{T^3}{6}\left(\frac{247\pi ^2}{15}+\underset{i}{}\kappa _i^2\right)\frac{\xi _{\nu _f}+\frac{\xi _{\nu _f}^3}{\pi ^2}}{\frac{43\pi ^2}{15}+\xi _{\nu _e}^2+\xi _{\nu _\mu }^2+\xi _{\nu _\tau }^2}.$$
(10)
Plugging quantum numbersI normalize the baryon charge so that quarks have baryon number $`1/3`$. of active species in eq. (1), using eq. (10) and the identities for electric charges of light quarks:

$$\underset{u,d,s}{}Q_i=0,\underset{u,d,s}{}Q_i^2=\frac{2}{3},\underset{u,d,s}{}Q_i^3=\frac{2}{9},\underset{u,d,s}{}Q_i^4=\frac{2}{9},$$
(11)
we get a system of five equations for dimensionless degeneracy parameters $`\kappa _\alpha =\mu _\alpha /T`$:

$`{\displaystyle \frac{8}{3}}\kappa _Q\kappa _{L^e}\kappa _{L^\mu }+{\displaystyle \frac{2}{9\pi ^2}}(\kappa _Q^3+\kappa _Q^2\kappa _B+\kappa _Q\kappa _B^2){\displaystyle \frac{1}{\pi ^2}}(\kappa _{L^e}\kappa _Q)^3{\displaystyle \frac{1}{\pi ^2}}(\kappa _{L^\mu }\kappa _Q)^3=0`$
$`\kappa _B^3+(6\kappa _Q^2+9\pi ^2)\kappa _B+2\kappa _Q^3=0`$
$`2(\kappa _{L^e}\kappa _Q)+{\displaystyle \frac{2}{\pi ^2}}(\kappa _{L^e}\kappa _Q)^3+\kappa _{L^e}+{\displaystyle \frac{1}{\pi ^2}}\kappa _{L^e}^3=\eta _e`$
$`2(\kappa _{L^\mu }\kappa _Q)+{\displaystyle \frac{2}{\pi ^2}}(\kappa _{L^\mu }\kappa _Q)^3+\kappa _{L^\mu }+{\displaystyle \frac{1}{\pi ^2}}\kappa _{L^\mu }^3=\eta _\mu `$
$`\kappa _{L^\tau }+{\displaystyle \frac{1}{\pi ^2}}\kappa _{L^\tau }^3=\eta _\tau ,`$ (12)
where

$`\eta _f`$ $`=`$ $`\left[{\displaystyle \frac{247\pi ^2}{15}}+{\displaystyle \frac{4}{3}}\kappa _Q^2+{\displaystyle \frac{2}{3}}\kappa _B^2+2(\kappa _{L^e}\kappa _Q)^2+2(\kappa _{L^\mu }\kappa _Q)^2+\kappa _{L^e}^2+\kappa _{L^\mu }^2+\kappa _{L^\tau }^2\right]`$ (13)
$`\times {\displaystyle \frac{\xi _{\nu _f}+\frac{\xi _{\nu _f}^3}{\pi ^2}}{\frac{43\pi ^2}{15}+\xi _{\nu _e}^2+\xi _{\nu _\mu }^2+\xi _{\nu _\tau }^2}}.`$
The system of equations (2) allows to expresses the degeneracy parameters at the time in the history of the Universe when the temperature was of order of few Gev in terms of neutrino degeneracy parameters now. This system, upon closer look, does not contain any large numerical coefficients. So, if $`\xi _{\nu _f}`$ are of order one, all $`\kappa _\alpha `$, including $`\kappa _B`$ and $`\kappa _Q`$, will also be of order one. For instance, in the linear approximation, $`\kappa _B=0`$ and $`\kappa _Q=1.44(\xi _{\nu _e}+\xi _{\nu _f})`$.

I have solved the system of equations for the degeneracy parameters numerically in two representative cases. In the first case, neutrino degeneracies were chosen ’democratically’ – equal for all three neutrino species: $`\xi _e=\xi _\mu =\xi _\tau \xi `$. The dependence of the baryon chemical potential and the chemical potential for electric charge at the QCD epoch on $`\xi `$ are shown in Fig. 1. Fig. 2 shows the same quantities in the second case, when I took $`\xi _e=0`$ and $`\xi _\tau =0`$.

## 3 Discussion

The main result of this paper is that an overall charge neutrality, equilibrium with respect to the weak interactions and a non-zero density of the lepton number imply that the quark-gluon sector of the primordial plasma in the early Universe carries large (comparable to the entropy) electric charge, which compensates the opposite charge in the lepton sector. It is necessary to mention that, in explicit calculations of the chemical potentials carried out in the previous section, the quark-gluon plasma was treated as an ideal gas, which is an accurate approximation only at sufficiently high temperature. In practice, the ideal gas approximation should work well only for extensive quantities (entropy, charge density, etc.) and at temperature of order of few Gev or higher. The most interesting phenomena happen at lower temperatures when the ideal gas approximation definitely breaks down. Still, the qualitative conclusion that the quark-gluon plasma is charged in the presence of lepton asymmetry should remain unchanged, because it relies solely on the fact that quarks carry electric charges and can be thermally excited above the QCD phase transition.

The physics of charged quark-gluon plasma can differ considerably from that of the better studied neutral plasma. The charge density potentially can affect the fate of various topological defects  and can change the nature of the quark-hadron transition. At zero chemical potentials, the quark-gluon and the hadron phases are likely to be smoothly connected, so that there is a smooth crossover instead of the phase transition . It is difficult to say without detailed study whether large chemical potentials can make the transition stronger, but if this happens, it can have important cosmological consequences.

### Acknowledgments

I am grateful to N. Agasian, R. Brandenberger, J. Ng and A. Zhitnitsky for discussions. This work was supported by NSERC of Canada, by Pacific Institute for the Mathematical Sciences and in part by RFBR grant 00-15-96557 for the promotion of scientific schools.


format: text

text:
# Destruction of CD4 T Lymphocytes Alone Cannot Account for their Long-term Decrease in AIDS

## 1 Introduction

The last few years have witnessed an ongoing debate over the kinetics of AIDS. Ho et al. and Wei et al. showed that virions and infected CD4 T cells (TL) have a very short life span, and that the HIV and the immune system are in a QSS during the latent stage. Many mechanisms were proposed to account for the slow decrease in the CD4 TL concentration during the latent stage; most of them involve the direct or the indirect destruction of the CD4 TL by the virus. For a review and a simplified network representation see Atlan 1996 . We here employ a simple mathematical formalism to show that none of these mechanisms can explain the long term decrease in the CD4 TL concentration, and we propose other mechanisms to account for this decrease.

## 2 QSS T Cell Concentration

In order to describe qualitatively the change in the overall CD4 TL concentration ($`T`$) during HIV infection, we will use a QSS formalism. The QSS formalism assumes that all the mechanisms affecting $`T`$ can be averaged as a function of $`T`$ 111This assumption is true for the virions and for the infected CD4 TL  and for all the cell types reacting directly to the virions or the infected cells.. Before infection, $`T`$ is determined by a production rate ($`a`$), a proliferation rate ($`b`$) and a death rate ($`\delta _T`$). (We will address the possibility of a homeostatic mechanism and saturation in the next section.)

$$\dot{T}=a+(b\delta _T)T=ayT;y=b+\delta _T$$
(1)
The steady state is $`T=a/y`$, and the TL concentration reaches this steady state within a time of $`O(1/y)`$ (Figure 1a). 222The precise solution to this equation is $`T=\frac{a}{y}+(T(0)\frac{a}{y})e^{yt}`$. During the latent stage of the infection, $`T`$ is lower than in the healthy state. If the decrease in $`T`$ QSS is due to a constant destruction of CD4 TL (directly by the virions or by anti CD4 CTL), the equation we get is:

$$\dot{T}=a_1yT\gamma T$$
(2)
Where $`\gamma T`$ represents the destruction rate. The steady state of this equation is $`T=\frac{a}{y+\gamma }`$, and $`T`$ reaches the steady state in a time $`O(\frac{1}{y+\gamma })`$, which is faster than $`O(\frac{1}{y})`$ 333The precise solution in this case is $`T=\frac{a}{y+\gamma }+(T(0)\frac{a}{y+\gamma })e^{(y+\gamma )t}`$. (Figure 1b). Thus, adding a constant destruction rate can lower the QSS concentration, but this alone cannot produce the long-term decrease in $`T`$ characteristic of AIDS; $`T`$ only stabilizes around a new steady state.

One might reason that AIDS could emerge if one assumes that the CD4 TL destruction rate changes as a function of $`T`$. In other words imagine that the concentration of the agent (cell or virion) destroying the CD4 TL depends on $`T`$. For example, if we have an agent $`D`$ with a proliferation rate proportional to $`T`$ and a constant destruction rate, the dynamics of $`D`$ would be described by:

$$\dot{D}=xT\delta _DD$$
(3)
The QSS concentration of D then would be $`D=\frac{x}{\delta _D}T`$. The evolution of $`T`$ would then be described as :

$$\dot{T}=ayT\frac{x}{\delta _D}T^2$$
(4)
This case can be generalized to:

$$\dot{T}=ayT\gamma T^n;n>1$$
(5)
Note that $`T`$ will still reach a steady state, but more quickly. The steady state value of $`T`$ will be in this case: $`T=\frac{a+(n1)\gamma (a/y)^n}{y+n\gamma (a/y)^{n1}}`$ which converges to the linear steady state for $`n=1`$. The time at which $`T`$ reaches a steady state is to first order $`O(\frac{1}{y+\gamma n(\frac{a}{y})^{n1}})`$ (Figure 1c). Thus, neither a constant rate of CD4 TL destruction nor a $`T`$-dependent rate can explain the long-term fall in TL concentration during HIV infection. This is true even if TL destruction is caused by some indirect mechanism, for example, apoptosis, destruction of CD4 TL by CTL, anergy or other means.

## 3 Saturation

The same reasoning applies even if the CD4 TL steady state is regulated by a homeostatic mechanism. In this case the proliferation rate can be higher than the death rate, and the total TL population is limited by a higher order mechanism; we thus replace eq 1 by:

$$\dot{T}=a+yT\gamma T^2$$
(6)
This equation can be solved analytically444The precise solution is: $`T(t)=\frac{y+\sqrt{y^2+4a\gamma }\mathrm{tanh}\frac{\gamma \sqrt{y^2+4a\gamma }t\sqrt{y^2+4a\gamma }C}{2\gamma }}{2\gamma }`$ and has 2 types of behaviors:

1. If the T cells are produced (in the thymus or elsewhere) at a constant rate $`a>yT`$, we can neglect $`yT`$ and the dynamics of $`T`$ will be determined by the putative regulatory factor $`\gamma T^2`$. This factor could be the target of a process that raises the value of $`\gamma `$ (Figure 2b), or lowers the value of $`a`$ (Figure 2c); but in either case the rate with which $`T`$ reaches a steady state is faster than in the healthy state (Figure 2a,2b)555The steady state is $`T=\sqrt{\frac{a}{\gamma }}`$. Starting from an initial value of $`T(0)`$ the initial evolution is: $`T(0)(T(0)^2\gamma a)t`$. Thus lowering $`a`$ or rising $`\gamma `$ will both lead to a lower steady state, but this steady state will be aproched faster than in the healthy case. Indeed if ones lower $`a`$ the evolution rate of $`T`$ during the last stages of the disease is slower than in the healthy case, but this slow evolution occurs only very close to the steady state, and has no effect in most stages of $`T`$ evolution.
2. If, in contrast, the main source of T cells is by way of the proliferation of existing T cells ($`a`$ is negligible), then the dynamics of $`T`$ will be dominated by the net proliferation rate $`y`$. In this case a lower steady state of $`T`$ can be achieved either by raising $`\gamma `$, or by lowering $`y`$ (that is, by adding a constant destruction rate.). As in the cases described above, a new steady state is achieved, but it is achieved faster than in the healthy case (Figure 3 b-e) 666The steady state is $`T=\frac{y}{\gamma }`$. Starting from an initial value of $`T(0)`$ the initial evolution is: $`T(0)(T(0)\gamma y)t`$. Thus lowering $`y`$ or rising $`\gamma `$ will both lead to a lower steady state, but this steady state will be achieved faster than in the healthy case..

To summarize, we can say that any direct or indirect destruction of CD4 TL will lower the value of the TL concentration at the steady state, but such destruction alone will not delay the time it takes to reach this steady state. Thus a destructive mechanism alone cannot account for the very slow decrease in T actually observed in AIDS.

Moreover, any mechanism describing $`T`$ as a function of time generates a concave curve. In other words, the computed rate of decrease of T as a function of time decreases. However, the typical observed curve of evolution of AIDS is actually convex (Figure 3f).

## 4 Slow Change in QSS

The QSS analysis presented in the previous sections shows that any factor we add to the destruction rate of TL will reduce their concentration at the steady state, but such a factor will also reduce the time it takes for $`T`$ to reach this steady state. Thus, although many mechanisms can explain the decrease in the TL concentration, none of those proposed until now are able to explain why the decrease evolves so slowly.

The previous analysis fails only if the QSS assumptions are wrong. These assumptions can be wrong only if the process explaining the decrease in TL acts on cells whose lifespan is very long, and only if this process is enhanced as the value of $`T`$ decreases.

Thus the CD4 TL concentration can decrease slowly only if all the following conditions are fulfilled (Figure 4):

* The mechanism of destruction does not act directly on TL.
* The destruction rate is very slow.
* The rate of destruction increases when the TL concentration decreases.

These requirements are foreign to most known mechanisms. For example, if TL lose is due to a constant destruction of TL, the number of T cells destroyed will merely decrease when their concentration decreases. Accordingly, the destruction rate can only decrease as the T cell concentration decreases. Thus, for the destruction rate to increase when the concentration of TL decreases, the destruction of TL must feature some positive feedback on itself, and the positive feedback must be stronger than the negative feedback due to the decrease in TL concentration.

## 5 Proposed Mechanisms

We propose here some simple mechanisms that can fulfill the requirements that we showed above to be necessary to explain the long-term evolution of $`T`$ during HIV infection. These include:

1. A slow evolution of the virus toward more destructive variants .
2. A reduction in the CTL proliferation rate by cytokine regulation.
3. A competition between humoral and cellular immunity.
4. A reduction in the rate by which B cells produce virus-neutralizing antibodies due to destruction of the functional lymph node architecture.

Here are some details:

* Cytokines - The effect of cytokines on HIV evolution is still controversial . The observed decrease in TH1 type cytokines and the increase in TH2 type cytokines might be either a marker of the decrease in the number of CTLs or its cause. The tendency of cytokines to positively influence the proliferation of the T cells that produce them can serve as a positive feedback loop which might balance the effect of the negative feedback caused by the loss of infected cells.Infected T cells secrete TH2 type cytokines and CTLs secrete TH1 type cytokines . During the latent stage, the CTL concentration decreases faster than does that of the infected cell population. Hence, the TH1 concentration decreases and the TH2 concentration increases. At some stage, there is an inversion in the cytokine profile and a collapse of CTL production. This leads to a rapid rise in the viral load and a fast decrease in CD4 due to accelerated infection.
* Competition between humoral and cellular immune responses. - The decrease in the concentration of TH1 cytokines and the increase in TH2 cytokine concentration can also be explained by an extrapolation of classical TH1/TH2 competition. The humoral and cellular immune systems act in parallel against HIV. Antibodies destroy the virus, while, CTLs destroy the infected cells. An increase in the B-cell concentration will lead to a decrease in the virus concentration leading, in turn, to a decrease in the concentration of infected cells. But the CTLs become activated by the infected cells. Consequently the CTL concentration will decrease when the infected cell concentration decreases. This mechanism leads to a competition in which raising the humoral response will decrease the cellular response and vice versa. This simple competition can be enlarged and strengthened if we assume that CTLs, or the CD4 TL that activate CTLs, secrete TH1 type cytokines and that B cells secrete TH2 type cytokines. Both cytokine types enhance their own production by T helper cells and limit the secretion of the opposite type. This competition can enlarge an existing competition between B and T cells leading to the collapse of the cellular immune response.
* B-Cell depletion - The destruction of lymph nodes and of their germinal centers can limit the proliferation both of B cells and of T cells. If B-cell antibody production is reduced, the number of free virions can increase despite a reduction in the number of infected cells. The action of such a mechanism might be deduced from the negative correlation between virus and CTL concentration. The CTL concentration is proportional to the infected cell concentration. Thus, a negative correlation between CTL and virus concentrations can be due to a negative correlation between the concentrations of virus and infected cells. The virus concentration rises while the infected cell and CTL concentration is either constant or decreases. This can be achieved only if the viral lifespan increases. This increase in lifespan could be due to a reduction in the capacity of B cells to make antibodies that neutralize HIV.

## 6 Conclusions

We have shown that no simple direct or indirect mechanism of TL destruction can explain the long-term decrease in CD4 TL concentration leading to AIDS. The long term decrease in CD4 TL concentration can be explained only by some positive feedback of the destructive mechanism on itself. We propose four mechanisms, which separately or together, can produce the required positive feedback. Understanding the precise mechanism that is responsible for the long-term decrease of CD4 TL is of vital importance. Treating this mechanism immunologically might help halt the progression to AIDS, in synergy with drug therapy.

## 7 Acknowledgment

Irun R.Cohen is the incumbent of the Mauerberg Chair in Immunology, Director of the Robrt Koch-Minerva Center for Research in Autoimmune Disease, and Director of the Center for the Study of Emerging Diseases.

Henri Atlan is the incumbent of the Frank Taper chair in medical diagnosis at the Hebrew University in Jerusalem


text:
# Phase diagram of a generalized Hubbard model applied to orbital order in manganites

## Abstract

The magnetic phase diagram of a two-dimensional generalized Hubbard model proposed for manganites is studied within Hartree-Fock approximation. In this model the hopping matrix includes anisotropic diagonal hopping matrix elements as well as off-diagonal elements. The antiferromagnetic (AF), ferromagnetic (F), canted (C) and paramagnetic (P) states are included in the analysis as possible phases. It is found that away from half-filling only the canted and F states may exist and AF and P states which are possible for the usual Hubbard model do not appear. This is because the F order has already developed for on-site repulsion $`U=0`$ due to the hopping matrix of the generalized model. When applied for manganites the orbital degree is described by a pseudospin. Thus our “magnetic” phase diagram obtained physically describes how orbital order changes with $`U`$ and with doping for manganites. Part of our results are consistent with other numerical calculations and some experiments.

The three-dimensional (3D) cubic manganites R1-xAxMnO3 (R is a rare earth element such as La or Nd, A is a divalent alkali such as Sr or Ca) and the layered ones Lan-nxSr1+nxMnnO3n+1 ($`n=1,2`$) have attracted intensive interest due to their rich material properties . Most prominently they show a colossal magnetoresistance which may be qualitatively understood on the basis of double exchange (DE) model or the more general ferromagnetic Kondo lattice model. However, many recent experimental findings such as the doping dependence of charge, orbital, spin ordered phases  suggest that additional physical mechanisms beyond the DE model should be involved. The possible ingredients are: $`e_g`$ orbital degeneracy , electron correlation  as well as electron-lattice interaction, e.g., by Jahn-Teller coupling . Though there exist some works to study the complex interplay of spin, orbital and lattice degrees of freedom , it is still very helpful to clarify the role of each individual mechanism separately with the assumption that other degrees of freedom are frozen out. Actually, experiments showed that the spins in many 3D manganites R1-xAxMnO3 are ferromagnetically ordered or A-type antiferromagnetic ordered (ferromagnetic layers coupled antiferromagnetically) in a large region of the doped phases, and A-type antiferromagnetic ordered in the undoped phase . This means that a perfect spin-polarized two-dimensional (2D) plane is always retained over large region of doping. Therefore, the spin degree of freedom may be discarded and the following effective 2D Hamiltonian with only orbital degree of freedom is proposed :

$``$ $`=`$ $`{\displaystyle \underset{<ij>}{}}{\displaystyle \underset{\sigma \sigma ^{}}{}}t_{ij}^{\sigma \sigma ^{}}c_{i\sigma }^{}c_{j\sigma ^{}}+U{\displaystyle \underset{i}{}}(n_i1/2)(n_i1/2),`$ (1)
where $`<ij>`$ indicates summation over nearest-neighbor sites, $`t_{ij}^{\sigma \sigma ^{}}`$ denotes the hopping integral and $`U`$ is the effective inter-orbital Coulomb interaction. The two orbitals $`|x^2y^2`$ and $`|3z^2r^2`$ have been assumed as pesudospin $``$ and $``$, respectively. The hopping matrix is explicitly given by $`t_{ij}^{}=t_1=3t/4,t_{ij}^{}=t_2=t/4,t_{ij}^{}=t_{ij}^{}=(+)t_3=(+)\sqrt{3}t/4`$ along the $`x(y)`$-direction . The most important feature here is the orbital dependence of the integral $`t_{ij}^{\sigma \sigma ^{}}`$, which distinguishes the present model from the usual Hubbard model where we have $`t_{ij}^{\sigma \sigma ^{}}=t\delta _{\sigma \sigma ^{}}`$.

Although the above 2D model ignores some important effects of three dimensionality as we discuss in the later, it may still contain important physics for 3D cubic manganites. For this the above model or its strong coupling version, the so-called orbital $`t`$-$`J`$ model has been used to study the unusual optical conductivity observed in doped LaMnO3 . At the same time it should be also noted that the above model may be of direct relevance to layered manganites, which consist of stacking of single or double MnO2 layers. Thus a full understanding of the properties of the above model is necessary. While its optical conductivity was studied by others, we discuss a different aspect here. It is natural to ask in which properties the present model deviates from the usual Hubbard model. For example, as a basic problem, what are the differences between their magnetic (it actually means “orbital” for the real manganites) phase diagrams? This will be the main topic addressed in this work. Though the orbital ordering for manganites based on the current model has been numerically studied at some points , the full phase diagram has not been discussed. So the purpose of this paper is twofold: On one hand we want to know how its magnetic phase diagram is modified by the generalized hopping matrix, in analogy to the recent investigation of the influence of next-nearest-neighbor hopping  on the phase diagram. On the other hand, as a real model for manganites, our “orbital” phase diagram obtained below is certainly useful for the nontrivial question of how the orbital order is changed upon doping.

In the following we use the Hartree-Fock (HF) approximation to study the magnetic phases of the Hamiltonian (1). Although the HF approximation may overestimate the appearance of the ordered states due to neglect of quantum fluctuations, this simple approximation may still give many useful information. So it should be kept in mind that the quantitative aspects of the results given here should not be taken literally but that these results provide a basis for further more elaborated theoretical studies. In this respect the so called slave-boson mean field theory is an improved method, but requires much more numerical efforts . Besides the usual ferromagnetic (F), antiferromagnetic (AF) ordered and paramagnetic (P) states considered here the canted ordered (C) state is also taken into account as a possible phase. Although it has rarely been considered in previous studies of the Hubbard model, we include it here motivated by the fact that experiments on LaMnO3  have confirmed the type of $`|3x^2r^2/|3y^2r^2`$ orbital ordering which corresponds to a canted spin order in the pesudo-spin representation , as illustrated in Fig. 1. Other complicated spin ordered states which have been studied for the usual Hubbard model such as spiral phase  etc. will not be considered.

In HF approximation the Hamiltonian (1) can be rewritten as  (irrelevant constants are ignored):

$`_{\mathrm{HF}}`$ $`=`$ $`{\displaystyle \underset{<ij>}{}}{\displaystyle \underset{\sigma \sigma ^{}}{}}t_{ij}^{\sigma \sigma ^{}}c_{i\sigma }^{}c_{j\sigma ^{}}{\displaystyle \frac{2}{3}}U{\displaystyle \underset{i}{}}(2\stackrel{}{S}_i\stackrel{}{S}_i\stackrel{}{S}_i\stackrel{}{S}_i),`$ (2)
where $`\stackrel{}{S}_i=\frac{1}{2}_{\sigma \sigma ^{}}c_{i\sigma }^{}\stackrel{}{\sigma }_{\sigma \sigma ^{}}c_{i\sigma ^{}}`$ and $`\stackrel{}{\sigma }`$ are Pauli matrices. We restrict to solutions with uniform electron density and canted order of the magnetic moments characterized by

$`\stackrel{}{S}_i=\{\begin{array}{cc}m(\mathrm{sin}\theta ,0,\mathrm{cos}\theta )\hfill & i\mathrm{A},\hfill \\ m(\mathrm{sin}\theta ,0,\mathrm{cos}\theta )\hfill & i\mathrm{B}.\hfill \end{array}`$ (5)
Here we have assumed a bipartite lattice with sublattices A and B. For later discussion it is indispensable to understand the nature of all possible ordered states characterized by the parameter $`\theta `$ (with $`m`$ describing the size of the order parameter), as well as their corresponding orbital ordering .

1. $`\theta =0(\pi )`$: F order with spin $``$ ($``$) or homogenous $`|x^2y^2`$ ($`|3z^2r^2`$) orbital ordering.
2. $`\theta =\pi /2`$: AF order or staggered ($`|x^2y^2+|3z^2r^2`$) /($`|x^2y^2|3z^2r^2`$) orbital ordering.
3. other $`\theta `$ values: canted-F or AF order. For example, $`\theta =\pi /3(2\pi /3)`$ corresponds to staggered $`|3x^2r^2/|3y^2r^2`$ ($`|y^2z^2/|x^2z^2`$) orbital ordering.

In addition we note that due to broken SU(2) symmetry by the hopping matrix in Eqs. (1,2) the direction of spin polarization in the ordered states is no longer rotationally invariant, contrary to the case for the usual Hubbard model. For example, the F phases with spin $``$ and $``$ are not equivalent now. We also remark that the spins are lying within the plane in the AF ordered state described by Eq. (5).

With new fermion operators $`a`$ and $`b`$ corresponding to sublattice A and B, respectively, we may write the Hamiltonian (2) in the momentum space as $`_{\mathrm{HF}}=2NUm^2/3+(a_k^{}b_k^{}a_k^{}b_k^{})^\mathrm{M}(a_kb_ka_kb_k)^\mathrm{T}`$ with $`^\mathrm{M}`$ denoting a $`4\times 4`$ Hermitian matrix. Here $`N`$ is the total number of lattice sites and the reduced Brillouin zone is $`|k_x|+|k_y|\pi `$ or $`\pi /2<k_1,k_2\pi /2`$ ($`k_{x,y}=\pm k_1+k_2`$). In principle $`^\mathrm{M}`$ leads to four energy bands $`\epsilon _k^\alpha `$ ($`\alpha =1,2,3,4`$). The ground state energy is given by $`E_{\mathrm{GS}}=_{k,\alpha }\epsilon _k^\alpha +2NUm^2/3`$, where the first contribution depends on the band filling $`n=N_e/N`$ where $`N_e`$ is the total electron number. Next we calculate numerically the values $`m`$, $`\theta `$ by minimizing the $`E_{\mathrm{GS}}`$ for given repulsion $`U`$ and band filling $`n`$ or doping $`\delta `$ ($`=1n`$). For this purpose we have used a mesh of $`10^2\times 10^2`$ $`k`$ points in the Brillouin zone. In the following discussion of results $`t`$ is set as energy unit.

With values of $`m`$, $`\theta `$ obtained from the minimization of EGS we may plot the magnetic phase diagrams in the plane ($`\delta ,U`$), as shown by Figs. 2 and 3. Although there are a lot of studies for the usual Hubbard model (i.e., $`t_1=t_2=1,t_3=0`$) , we give its phase diagram again in Fig. 2 since the canted ordered state is added here. We also want to compare it with the result for the current generalized Hubbard model whose phase diagram is displayed in Fig. 3 (solid line). In the present model, the hopping matrix involves anisotropic diagonal elements and off-diagonal ones. To see more clearly how the phase diagram is influenced by them, we plot it for a simplified model with a hopping matrix which only includes anisotropic diagonal elements (i.e., $`t_1=3/4,t_2=1/4,t_3=0`$); it is denoted by a dashed line in Fig. 3.

¿From Figs. 2 and 3, a number of differences can be seen. Firstly, the phases for the usual Hubbard model include AF, F, P states as well as a narrow region with a canted phase. It appears in a crossover region from AF to F phase, which is understandable. The phase diagram is then qualitatively modified by introduction of anisotropic diagonal elements where only the canted and F ordered phases exist and the original AF and P phases are excluded away from half-filling (i.e., $`\delta >0`$). Furthermore when the off-diagonal elements are introduced the phase diagram is mainly quantitatively modified, i.e., the stability region of the F ordered phase is enlarged. We point out that the F phase in Fig. 3 always means that spins align in the up direction ($``$) due to a larger hopping matrix element $`t_{ij}^{}`$ (i.e., $`t_1>t_2`$).

The above results may be understood if we consider the property of the kinetic term alone, i.e., assume that $`U=0`$. Then the Hamiltonian (1) may be exactly diagonalized with two bands:

$`ϵ_k^\pm `$ $`=`$ $`(t_1+t_2)(\mathrm{cos}k_x+\mathrm{cos}k_y)\pm \sqrt{(t_1t_2)^2(\mathrm{cos}k_x+\mathrm{cos}k_y)^2+4t_3^2(\mathrm{cos}k_x+\mathrm{cos}k_y)^2}.`$ (6)
¿From $`ϵ_k^\pm `$ the $`z`$-axis magnetization $`M_z=c_i^{}c_ic_i^{}c_i/2`$ may be calculated, as shown in Fig. 4 where the curves $`M_z/n`$ vs $`\delta `$ are given. For the usual Hubbard model $`M_z`$ is always zero, corresponding to P state for the abscissa in Fig. 2. However, a finite magnetization (with spin $``$) is always exhibited for the present model as shown by the solid line in Fig. 4; also for the simplified model (zero off-diagonal hopping) shown by the dashed line, except at the point $`\delta =0`$. Moreover, the $`M_z`$ value for the current model is always bigger than that for the simplified one, which means that the F order of the former is stronger. Because of the underlying F order already in the $`U=0`$ case only the canted (with nonzero F component involved) and F ordered phases may exist in finite $`U`$ region away from half-filling. It is also easily understood why the F phase region is enlarged after inclusion of the off-diagonal hopping matrix elements.

We consider the case of half-filling (i.e., $`\delta =0`$). It has been shown that the usual Hubbard model has an AF ground state for any finite $`U`$, corresponding to the ordinate in the phase diagram of Fig. 2. We speculate that this conclusion might also be true for our present generalized Hubbard model. Although the present ground state is F ordered at $`U=0`$ as given in Fig. 4, it is found to have AF order for $`U>1.2`$. This was concluded from a calculation employing a mesh of $`10^3\times 10^3`$ $`k`$ points. Whether or not there exists a crossover between F and AF for finite $`U`$ is not yet certain.

Our “magnetic” phases obtained from Hamiltonian (1) actually correspond to the orbital orderings relevant for the manganites as mentioned before, therefore we now discuss them in detail. To see clearly how the orbital order changes with doping $`\delta `$ we have plotted in Fig. 5 the relation of parameter $`\theta `$ with $`\delta `$ for fixed $`U`$. It is found $`\theta `$ decreases monotonically from $`\pi /2`$ to $`0`$ with increase of $`\delta `$ indicating the process AF$``$ C$``$ F order for the pseudospins as shown in Fig. 3. Let us take $`U=10`$ as an example in the following. (The real value for it in manganites is about $`410`$). In the undoped case, $`\theta =\pi /2`$ means a staggered ($`|x^2y^2+|3z^2r^2`$)/($`|x^2y^2|3z^2r^2`$) orbital ordering, which is consistent with the result by Motome and Imada from quantum Monte Carlo calculation, but slightly different from the experimentally observed $`|3x^2r^2/|3y^2r^2`$ . This may be due to some other physical mechanism, e.g, a Jahn-Teller effect which is present in real materials but not included here . On the other hand, it is found that the $`|3x^2r^2/|3y^2r^2`$ type orbital ordering may be also exhibited from our pure orbital model at very small doping $`\delta 0.06`$. With increase of doping the orbital order changes continuously until at some critical value $`\delta _c0.12`$ which is still small it becomes $`|x^2y^2`$ type. This behaviour for small $`\delta _c`$ has also been noted by Horsch et al. based on the orbital $`t`$-$`J`$ model (i.e., large $`U`$ version of ours) with the exact diagonalization method . Interestingly it is also seen that the value $`\delta _c`$ is not a monotonic function of $`U`$, which reaches a maximum $`0.2`$ at about $`U=5`$. Recently Akimoto et al. found an A-type antiferromagnetic metallic ground state for La1-xSrxMnO3 at doping $`x>0.5`$, and they proposed that $`|x^2y^2`$ orbitals should be occupied in this state. This is consistent with our result derived from the planar model (1). In addition, we also noticed the experimental findings that many cubic manganites are 3D ferromagnets in the intermediate doping regime nearby $`\delta =0.3`$ , where the hopping along the $`z`$-axis is allowed by the DE mechanism. Thus it is instructive to discuss the effect of three dimensionality neglected in our model in this doping region, see also Ref. . For a 3D model, the hopping along the $`z`$-axis is orbitally highly anisotropic; only the hopping between two $`|3z^2r^2`$ orbitals is permitted. So the 2D ferro-orbital order with $`|x^2y^2`$ type polarization, as found above, may be suppressed in order to gain the kinetic energy along the $`z`$-axis. Such a depolarization may be simulated by introducing a negative chemical potential $`\mu `$ on the $`|3z^2r^2`$ orbital . It is found that, with reasonable $`\mu `$ valid for 3D manganites the disordered orbital liquid state may appear at some doping nearby $`\delta =0.3`$. This orbital liquid was first discussed by Ishihara et al. , which may be enhanced by the quantum fluctuations beyond HF approximation . On the other hand, for the layered manganites with single or double layers of MnO2 where the inter-layer or bilayer hopping (or correspondingly $`\mu `$) is assumed very small, the orbital liquid state may be excluded and our ferro-orbital ordered state is still expected to exist in the doping region mentioned above.

In the above discussion we have shown the possible orbital ordering as function of doping $`\delta `$ and repulsion $`U`$. The results are partly consistent with previous numerical calculations, and also with some experimental results on manganites although only orbital degree of freedom is retained in the present model. For a full understanding of various phases observed in manganites experimentally, one needs to consider the spin and orbital orderings simultaneously , as well as include lattice degree of freedom in some cases . However, our orbitally ordered structures obtained above are very instructive to understand the role of orbital degree of freedom itself, i.e., how the interplay of the orbital kinetics (with attention to its peculiar matrix character) and orbital correlation induces the orbital ordering. Moreover, they may be also possibly observed under some conditions, e.g., at large external magnetic fields.

In summary, within HF approximation we have studied the magnetic phase diagram of a generalized Hubbard model proposed for manganites where the hopping matrix includes anisotropic diagonal elements and off-diagonal elements. The AF, F, canted ordered and P states are considered as possible phases. It is found that in the doped case only the canted and F states may exist, while the AF and P phases of the usual Hubbard model are not favored. This is because the F order has already been established at $`U=0`$ due to the peculiar hopping matrix. When applied for manganites, the spin index of this model actually represents the orbital degree of freedom. Thus our obtained magnetic phase diagram is able to describe how the orbital order changes upon doping.

Two of authors (Q. Yuan and T. Yamamoto) would like to thank the support of Visitor Program of MPI-PKS, Dresden. Q. Yuan also acknowledges the part support by Chinese NSF.

Note Added in Proof: Recently the scenario for phase separation has been extensively studied for various models in the manganites . Such a tendency of separation into orbitally ferro- and antiferro- ordered regions may also exist in our case of the canted phase in Fig. 3. This problem deserves further detailed investigation.


text:
# The Interaction of Type Ia Supernovae with their Surroundings: The Exponential profile in 2D

## 1 INTRODUCTION

The evolution of supernova remnants (SNRs) has been classically assumed to go through 4 stages (see for example Woltjer 1972). In the initial free-expansion or ejecta-dominated stage, the expansion of the remnant is determined by the properties of the supernova ejecta and the density of the surrounding medium. The mass of the ejected material, the kinetic energy of the explosion, the density profile of the ejected material and the density distribution in the ambient medium are all required to understand the expansion of the remnant in this stage. This phase is particularly important because many of the so-called “historical” supernovae (SNe) known in our galaxy may either be in this stage, or be making a transition from it to the Sedov or adiabatic phase, where the mass swept-up by the SN shock wave substantially exceeds the ejecta mass.

Determination of the density profile of the ejected material is not a simple task, and little observational evidence is available. The density depends on the structure of the progenitor star and the shock acceleration of the gas during the explosion (Chevalier & Fransson 1994). The best available information is for SN 1987A, where the density seems to be well constrained by models of the explosion (eg. Arnett  and Shigeyama & Nomoto ). In this case the density in the outer parts of the ejecta in the free-expansion phase varies approximately as a power-law with velocity. Similar models have been used to describe Type II SNe in general, and self-similar solutions for power-law ejecta interacting with a power law or constant density surrounding medium have been derived by Chevalier (1982a) and Nadozhin (1985). Recently, Truelove & McKee (1999) have presented more general approximate solutions for the evolution of non-radiative SNRs in the ejecta-dominated stage. The special case of interaction of a uniform density ejecta with the surrounding medium will also be addressed in Drury and Dwarkadas (2000).

Power-law models (with an exponent of 7) have also been used to describe the ejecta density profile in Type Ia SNe (Chevalier 1982a, 1982b; Fabian, Brinkmann & Stewart , Band & Liang ), which are believed to result from the degenerate ignition of a C/O white dwarf in a binary system. However Dwarkadas & Chevalier (1998, hereafter DC98) showed that an exponential ejecta profile for Type Ia SN explosions gives the best fit to the density profile. Furthermore, they compared the interaction of Type Ia SNe described by power-law, exponential and constant density ejecta profiles with the surrounding medium, in the simplest case of spherical symmetry.

This paper expands on the work of DC98, presenting high-resolution two-dimensional simulations of ambient medium interaction in Type Ia’s using an exponential density profile to describe the ejecta. We compare the properties of the two-dimensional flow to that in one-dimension. The 2D grid allows us to investigate the multi-dimensional aspects of the flow, especially the formation of hydrodynamic instabilities in the free-expansion phase. In particular we study the evolution and growth of the Rayleigh-Taylor (R-T) instability that arises at the decelerating contact discontinuity.

The role of the R-T instability in young supernova remnants has been receiving great attention in recent years (see Chevalier 1996). The accumulation of observational evidence of substantial non-radial motions and mixing occurring in the explosion of SN 1987A led several groups to perform multi-dimensional simulations of supernova explosions (Arnett et al. 1989; Benz & Thielemann 1990; Fryxell et al. 1991; Herant & Benz 1991). Hydrodynamic instabilities in self-similar driven waves were investigated by Chevalier, Blondin & Emmering (1992, hereafter CBE92). The presence of instabilities can lead to fragmentation and the formation of clumps or blobs in the interface between the ejecta and the surrounding material. Mixing among various layers can destroy any early stratification and cause certain elements to be found at larger radii than would otherwise be expected. In this paper our aim is to characterise the nature and growth of these instabilities for the case of SN ejecta described by an exponential density distribution. Implications of mixing of the ejecta and circumstellar (CS) material for the observable supernova remnants, in particular Tycho’s SNR (SN1572), is also addressed. The growth of these instabilities under various conditions is investigated, such as a constant density medium (indicative of the interstellar medium) or a circumstellar medium formed by mass loss from the progenitor star. We compare and contrast the results in these two cases.

The rest of this paper proceeds as follows. In §2 we first review briefly the properties of an exponential ejecta density profile. This is followed by initial results of the 2D simulations in §3. In §4 we discuss the R-T instabilities in more detail. §5 focuses on the application of our results to Type Ia SNe, especially Tycho’s SNR, while §6 contains a brief discussion and summary.

## 2 The Exponential Ejecta Profile

The exponential profile used to describe the SN ejecta can be written as (DC98):

$$\rho _{SN}=A\mathrm{exp}(v/v_e)t^3$$
(1)
where, by integrating over the total ejecta distribution, taken to have mass $`M_e`$ and kinetic energy $`E`$, we find that

$$v_e=\left[\frac{E}{6M_e}\right]^{1/2}=2.44\times 10^8E_{51}^{1/2}\left[\frac{M_e}{M_{ch}}\right]^{1/2}\mathrm{cm}\mathrm{s}^1,$$
(2)
where $`M_{ch}1.4M_{}`$ is approximately the Chandrasekhar mass, and $`E_{51}`$ is the explosion energy in units of 1051 ergs. The value of $`E_{51}`$ quoted in the literature lies anywhere between about 0.3 and 1.8. The parameter $`A`$ is given by:

$$A=\frac{6^{3/2}}{8\pi }\frac{M_{e}^{}{}_{}{}^{5/2}}{E^{3/2}}=7.67\times 10^6\left[\frac{M_e}{M_{ch}}\right]^{5/2}E_{51}^{3/2}\mathrm{g}\mathrm{s}^3\mathrm{cm}^3.$$
(3)
The interaction of the SN ejecta with the ambient medium gives rise to a double shocked structure, consisting of a forward shock expanding into the surrounding medium, and a reverse shock that moves inward in a Lagrangian sense, separated by a contact discontinuity. We have carried out 2-dimensional hydrodynamic simulations of the supernova ejecta colliding with the ambient medium. The assumption of spherical symmetry with regard to the exponential ejecta density distribution is certainly an approximation. In reality the formation of instabilities during the explosion will probably destroy any such early symmetry, leading to large-scale, possibly turbulent convective motions (M$`\ddot{\mathrm{u}}`$ller 1998). Nevertheless, this assumption is a valid starting point for such a study, and the inclusion of clumps and non-radial flows would be the subject of future, more complicated studies (see Wang & Chevalier 2000)

Our lack of knowledge of the progenitors of Type Ia SNe (eg. Livio 1999) makes it unclear what the density structure is in the immediate surroundings. The lack of radio and X-ray emission, and of nebular excitation lines from circumstellar interaction, points to a medium of constant low density. However if Type Ia’s are formed from degenerate ignition in a white dwarf, then it is clear that the white dwarf must be in a binary system, and mass transfer must be occurring (Wheeler 1996). A circumstellar region due to mass loss from the secondary star would then be expected, although its extent may be small. Cumming et al. (1996) place limits on the properties of such a region. In this paper we consider both cases, that of a constant density surrounding medium and a circumstellar medium. If the mass loss rate and velocity of the wind from the companion star are constant, the density of the circumstellar medium (CSM) decreases inversely as the square of the radius. The supernova remnant will remain in such a medium only for a few decades at most. Nevertheless, as shown by DC98, the passage of the shock wave through such a CS region can leave an impression on the density structure of the interaction region even after the shock wave has exited the CS region. Thus it is both observationally relevant as well as academically interesting to study the evolution of the shock within a circumstellar medium. Furthermore, while an exponential profile has at present been postulated only in relation to Type Ia SNe, it must be emphasised that it is really not clear what the structure of the ejecta density profile is for any supernova other than possibly SN 1987A, despite recent progress in this direction (see for instance Matzner & McKee ). A detailed study of the exponential ejecta density expanding into various ambient media is thus warranted.

## 3 2D Simulations

The simulations were carried out using the VH-1 code, a 3-Dimensional finite-difference hydrodynamic code based on the Piecewise Parabolic Method (Colella and Woodward 1984), written by the Virginia Hydro group. The code, obtained from Dr. John Blondin at NCSU, solves the equations of conservation of mass, momentum and energy in Lagrangian co-ordinates, and then remaps the variables back to the original Eulerian grid. It has been extensively tested on various astrophysical hydrodynamical problems. The version employed here uses an expanding grid that tracks the outer shock front and expands along with it, making it ideally suited to problems where the dimensions change by many orders of magnitude over the course of the run.

The evolution of the outer shock front in 2D is similar to that in 1D, as described by DC98. One important quantity is the evolution of the expansion parameter with time. If the radius is proportional to time raised to some power (R $``$ tδ), then $`\delta `$ is the expansion parameter. This is a dimensionless quantity that can be used to discriminate between the various phases of remnant evolution. In Figure 1 we show the evolution of the expansion parameter with time for a Type Ia SN evolving in a constant density medium. The time axis is normalised to the quantity $`T^{}248E_{51}^{0.5}\left(\frac{M_e}{M_{ch}}\right)^{5/6}n_0^{1/3}\mathrm{yr}`$ using equation (4) in DC98. The number density $`n_0`$ is related to the ambient medium density $`\rho `$ by $`n_0`$ = $`\rho `$/2.34 $`\times 10^{24}`$ cm-3, appropriate for a medium with a H/He ratio of 10:1. In Figure 2 we show the same plot for a CS medium whose density decreases as $`r^2`$. In this case the time normalisation is taken to be $`T_2^{}=E^{0.5}M_{e}^{}{}_{}{}^{1.5}D^1`$, where $`D=\rho r^2`$. These plots show how the expansion parameter varies over the evolution of the remnant in the free expansion phase, until it reaches the Sedov value, which is 0.4 in the case of the constant density medium and 0.67 in the case of a circumstellar medium. The continuous evolution of the expansion parameter is consistent with the higher-resolution 1D calculations of DC98.

As the outer shock expands into the circumstellar medium, it begins to sweep up the ambient material. The Sedov value of the expansion parameter in either case is reached only once the swept-up mass significantly exceeds the ejecta mass, by a factor of 30 or more. However, once the mass of swept-up material becomes comparable to the ejecta mass, the pressure of the surrounding medium causes the SN ejecta to decelerate. As the mass of shocked ambient material continues to increase, the reverse shock begins to move inwards towards the center, and the contact discontinuity between the two shocks begins to decelerate. The high density shocked SN ejecta behind the contact discontinuity are decelerated by the swept-up material ahead of it. The pressure and density gradients lie in opposite directions, and the system becomes unstable to Rayleigh-Taylor (R-T) instabilities. Small-amplitude perturbations, with wavelength on the size of a few grid zones, then occur. The size of these perturbations appears to grow linearly, and soon the familiar R-T fingers spread out into the shocked surrounding medium. The tips of these fingers terminate in bulbous-shaped mushroom-like heads, a result of the Kelvin-Helmholtz instability as the shocked gas surrounding the fingers tries to slide past them. This evolution is depicted in Figures 3 and 4, which are discussed in detail below.

## 4 The Rayleigh-Taylor instability

The Rayleigh-Taylor instability arises at the interface between a low-density fluid supporting a high-density fluid (e.g. Chandrasekhar 1961). In our case the analogous situation arises at the contact discontinuity separating the shocked SN ejecta from the shocked ambient medium, as described above. Small ripples arise initially along the entire length of the contact discontinuity. No explicit perturbation is required to induce the instability, which can be triggered purely by the numerical noise in the system. We have also tested the case when explicit perturbations of various amplitudes are introduced (§4.2). The ripples in our higher resolution simulations appear to start off close to the symmetry axis and then move towards the equator, but this is not easily confirmed, since by the time the ripples acquire an appreciable amplitude they are spread over the entire contact discontinuity. It is however quite likely that approximations at the boundary lead to minor fluctuations that trigger the instability.

Figure 3 shows grey-scale images of the evolution of the instability with time, for a SN evolving in a constant density medium, which we refer to henceforth as Case 1. Figure 4 shows the same for a Type Ia SN evolving in a circumstellar medium, referred to as Case 2. The times listed for each plot are normalised in the same manner as for the expansion parameters (§3). The radius for Case 1 is normalised to the value $`R^{}=\left(\frac{M_e}{\frac{4}{3}\pi \rho }\right)^{1/3}2.19\left(\frac{M_e}{M_{ch}}\right)^{1/3}n_0^{1/3}\mathrm{pc}`$. For Case 2 the radius normalisation is $`R^{}=M_e/D`$, where $`D=\rho r^2`$. The growth of the finger size (the penetration of the denser layer into the less dense layer) with time for Case 1 is shown in Figure 5. The plot for Case 2 looks quite similar when normalized to the appropriate units and is not shown.

Theory predicts that the size of the R-T fingers should increase as the square of the time. However a simple quadratic fit was found to be inadequate to model the growth, while a second order least-squares polynomial fit was found to provide a reasonable solution. The deviations are much larger earlier on, which is partly due to our inability to determine the amplitude of the perturbation in the very early stages. The fact that a linear as well as a quadratic component is needed is consistent with the fact that the flow is not self-similar (eg. Dalziel et al. 1999).

A standard linear analysis of the R-T instability shows that the growth rate is proportional to the square-root of the wavenumber (eg. Chandrasekhar 1961; Sharp 1984). According to this the smallest wavelengths should be the most unstable. In our simulations we find that the initial ripples are on the scale of a few (5-10) grid zones. Smaller wavelengths are damped by any inherent viscosity in the numerical method. The ripples grow into spike-like protrusions, forming the characteristic Rayleigh-Taylor ‘fingers’. As the spiky protrusions expand outwards into the surrounding medium, the shocked ambient material tries to slide past them. This gives rise to shear flow at the edges, leading to the growth of Kelvin-Helmholtz (K-H) instabilities, visible as the bulbous, mushroom-shaped heads on the R-T fingers. The flow of the shocked ambient medium past the fingers also tends to limit their growth, more so in the constant density medium case. Over time the fingers appear to merge into one another and the fastest growing wavelength increases with time, or correspondingly the number of individual R-T fingers seen decreases with time. The scale of the fastest-growing wavelength is set by the system, independent of the initial perturbation that triggered the instability. The result is that as the reverse shock approaches the center, the unstable wavelength, or alternatively the number of R-T fingers seen, is the same no matter what the amplitude of the initial perturbation, or even if no initial perturbation was applied to the system.

All our simulations show a jet-like structure along the symmetry axis. This structure is even more prevalent in those simulations where an initial perturbation was introduced. When the jet breaks through the outer shock it tends to “fan out”, thus giving rise to a large mushroom-shaped cloud. These structures, prevalent in many other numerical simulations embodying the same symmetry, are most likely numerical artifacts and are not real. They are formed at the symmetry axis, where the geometry results in grid zones whose volume is very small. Any fluid component perpendicular to the axis of symmetry (i.e. flowing into the axis) has nowhere to go and is deflected along the axis, thus resulting in these spurious structures.

We elaborate later on the differences between the instabilities seen in the constant density and the 1/$`r^2`$ surrounding medium, after first discussing the similarities between the two cases. A point to note in all our simulations is that, irrespective of the initial amplitude of the perturbation or of the nature of the surrounding medium, the unstable region never occupies more than about 50% of the interaction region. The size of the R-T projections grows faster than the size of the interaction region, so that the width of the mixing region tends to increase with time. However in all cases the reverse shock reaches the center, and the instabilities ceased to grow, when the mixing width was less than or about half the size of the interaction region. One implication of this is that the instabilities never penetrate the outer shock, and in fact have no direct effect on it. The appearance of the outer shock is not changed by the considerable amount of mixing occurring in the interaction region close to the contact discontinuity. The same is not true of the reverse shock in the circumstellar case, as described later.

While the above results are consistent with those found by CBE92 in their study of the Rayleigh-Taylor instability for power-law ejecta, a major difference exists between the two cases. The latter describe the quasi-steady evolution of the instability in the self-similar driven case, and in fact show (their Figure 7) how snapshots of the evolution taken at times two orders of magnitude apart hardly reveal any difference. The same is not true in the exponential case. Introduction of the velocity scale length in the exponential profile means that the solution is no longer self-similar, even though it is scalable. The appearance of the unstable region tends to change, albeit slowly with time, and no two snapshots which are two decades apart in time resemble each other closely as in the models of CBE92.

A close look at Figures 3 and 4 shows that there is considerable difference between the growth of the instabilities in the two cases. In the constant density medium case the growth of the fingers appears stunted. The formation of the mushroom caps at the head of a finger tends to limit the growth of the projection. The projections grow to a much larger extent in the CS case, and occupy a greater fraction of the interaction region. We attribute these differences in large part to the difference in Atwood number between the two cases (eg. Sharp 1984). The Atwood number, or equivalently the density contrast between the shocked ambient and shocked ejecta material, governs the growth of the instabilities for small-amplitude perturbations of wavelength close to the critical wavelength required for unstable behaviour (Sharp 1984). The Kelvin-Helmholtz instability caused by the shocked fluid flowing past the spikes is more pronounced at lower values of the Atwood number, which is the case when the shock expands in a constant density medium. This can cause the caps to mushroom out, increasing the effect of the drag on the spike and slowing its growth. This difference in the two cases is similar to that seen by CBE92 in their study of R-T instabilities for self-similar driven waves.

The two cases are further differentiated by the effect of the instabilities on the reverse (inner) shock. In the constant density medium case, the contact discontinuity is far enough away from the reverse shock that the instabilities do not affect the shape of the shock, which stays smoothly spherical. However in the CS case the reverse shock is much closer to the contact discontinuity, and the pressure variations are higher, causing the reverse shock itself to get wrinkled. The corrugation of the inner shock (Figure 4, middle panels) occurs on a scale similar to the most unstable wavelength, and seems to get more pronounced as the reverse shock heads towards the origin. Again, a similar effect was noticed by CBE92.

Note also that the time taken for the reverse shock to reach the centre in Figure 3 is about two-thirds of that found in the work of DC98. This is to be expected. The instabilities cause a “leakage” of pressure from the shocked ambient medium to the shocked ejecta. Thus the pressure driving the reverse shock is higher than in the spherically symmetric work of DC98, and the reverse shock consequently reaches the origin earlier.

The most likely scenario for supernova shock evolution is that it traverses a region of circumstellar medium followed by a constant density interstellar medium. As the shock moves from the CS to the constant density medium, the nature of the instabilities will change. The details will depend strongly on the extent of the CS region. If the size of the initial CS region is small, such that the instabilities do not have time to grow, then the effect of the CS region on the later evolution of the remnant is small. However we do find that in most cases the reverse shock remains corrugated even after the remnant is expanding primarily in a constant density medium. Thus, in addition to the spike in the density profile at the contact discontinuity, and the corresponding temperature drop mentioned in DC98, the shape of the reverse shock itself is affected by passage through an initial CS region. The presence of a CS shell at the interface between the two media was not taken into account and would further affect the evolution of the instabilities.

### 4.1 Effect of Grid Resolution

The effects of grid resolution on the instabilities merits consideration. Increasing the resolution leads to the appearance of structure on smaller scales, as would be expected. In Case 2 the instabilities appear to clump together with a number of sub-fingers arising from each finger, instead of appearing as individual spikes. The turbulent mixing behaviour, with vortices being formed close to the projecting fingers, causes individual fingers to be bent at an angle, producing structures resembling question marks. The stems show the usual Kelvin-Helmholtz mushroom caps due to shear flows around the edges. The number of R-T clumps that exist when the reverse shock reaches the center remains the same in all simulations, although the substructure in each clump continually increases. The appearance somewhat resembles a tree-structure, with higher resolution revealing more branches in the tree and more stems arising from the branches (see top panels in Figure 6). The pattern is fractal in nature, although the code is unable to resolve individual stems very well. The mixing width (the width of the region occupied by the instabilities) shows a definite increase of a few percent with resolution in the earlier stages of the instability. This may be partly due to the fact that as the shock expands outwards, individual structures are more easily discernible in a higher resolution simulation as the increase in the size of a grid zone is not as large as in a lower-resolution simulation. As the instabilities evolve the mixing width occupied by the unstable region becomes approximately constant.

In Case 1 a similar overall effect is visible. Simulations with a higher grid resolution lead to more R-T fingers being visible. As the simulation evolves the amount of substructure also increases, as would be expected. Also visible is an increase in the size of the K-H caps at the heads of the stems in the higher resolution cases. As before however, the simulations appear quite similar once the reverse shock nears the inner boundary, and the number of R-T fingers is more or less the same. No fractal structure is seen, in contrast to case 2, although the fingers are much better defined with an increase in the number of grid zones, as would be expected. Another effect is that the instabilities appear to start earlier as the spatial resolution increases. An increase in grid resolution by 50% leads to a decrease in the time when the ripples are first seen by about 20%. Also important is the fact that, as the growth rate is limited, a minimum resolution is required to capture the growth of the instabilities. This is consistent with the results of Fryxell et al. (1991) for SN 1987A. We have found that even a 300 $`\times `$ 300 run does not generate the instabilities till the reverse shock is close in to the centre of the explosion, and that a grid of at least 500 azimuthal zones is necessary to capture the instabilities in detail.

### 4.2 Effect of initial perturbation

We have investigated the effect of adding an initial perturbation of various amplitudes. A sinusoidal perturbation of wavelength equal to the angular dimensions of the grid was added. Figure 6 shows snapshots at approximately the same epoch in time from various runs for Case 2. Two cases with no perturbation added but different resolutions are shown in the first row. The remaining frames show cases where a perturbation of the given amplitude was introduced, all with the resolution set to 500 by 500 zones. In general we find that for Case 2 the overall effect of small-amplitude perturbations ($`10\%`$) is tiny, and the instabilities develop in almost the same way as they do in the case of zero initial perturbation. Initially the differences appear large, with the no-perturbation simulations showing a more homogeneous growth of projections whereas simulations with initial perturbations show fewer and more well-defined fingers. However the differences decrease as the reverse shock moves inwards, and in most cases the appearance of the unstable region looks quite similar by the time the reverse shock nears the inner edge of the grid. The effects of perturbations of amplitudes $`15\%`$ are more pronounced. There is an increase in the amount of small-scale structure to be seen, and the larger perturbation simulation presents a more clumpy appearance. A small increase in the mixing width accompanies the increase in the amplitude of the perturbation. It is evident from Figure 6 that the effects of higher resolution and introduction of a perturbation are not independent, in the sense that increasing the resolution has a similar effect to adding a perturbation. The higher resolution case with no perturbation is quite similar to the runs with an initial perturbation.

For Case I the effect of adding an initial perturbation is more pronounced, as shown in Figure 7. The inhomogeneous appearance of very small wavelength projections is replaced by individually discernible fingers. The projections seem to stretch out further, but their number has decreased. Thus, although there is a definite increase in the mixing width, the volume occupied by the instabilities remains about the same. Again the differences decrease as the simulation proceeds, although even when the inner shock is close to the center the differences are immediately perceptible, unlike in Case 2. Another distinction is that with the addition of an initial perturbation, the instabilities form much earlier than in the case of no perturbation, and a lower grid resolution is adequate to show that the instabilities are present, although a higher resolution grid is still required to resolve them. A grid with 300 azimuthal zones is sufficient to at least demonstrate the growth of the instabilities (Fig 7).

## 5 Tycho’s SNR (SN 1572)

Rayleigh-Taylor instabilities have been often cited as the cause of corrugations and rippled structure in images of SNRs. Velazquez et al. (1998) have studied the Rayleigh-Taylor instability in Tycho’s SNR, following up on VLA observations carried out by Reynoso et al. (1997). They conclude that wavy structure seen in the NE quadrant of the VLA image, formed by small, regularly spaced protrusions just behind the outer shock front, is due to Rayleigh-Taylor instabilities. It is unclear if the results of our work support such an interpretation. As emphasised in the previous section, in no case did the instabilities in our simulations cover more than half the interaction region, let alone penetrate the outer shock front, as do some of the protrusions shown in Figure 2 of Velazquez et al. (1998). If the radio emission is coming from just behind the outer shock our simulations imply that it is unlikely that any unstable structure should be visible close to to the outer rim. Moreover, the calculations of DC98 and our own 2D calculations show that a tenable model for Tycho requires that the reverse shock be quite far out, and the size of the interaction region be small. The instabilities will then still be in the incipient stages, as is also found by Velazquez et al. (1998). However that would mean that they would occupy an even smaller fraction of the interaction region, and would not likely be as prominent as they are in the radio images, unless some effect led to an increase in entropy within the mixing region, leading to a larger growth of the R-T instabilities than is observed in our simulations. The latter is not entirely improbable - it is possible perhaps that thermal instabilities in the surrounding medium, operating over evolutionary timescales, set up a very clumpy and non-linear initial condition for the growth of these instabilities, resulting in much larger growth than that studied herein. There exists also the possibility that if the NE portion is running into a high density H1 cloud as revealed by Reynoso et al. (1997), some of the corrugated shock structure may be due to interaction of the remnant with an uneven cloud boundary.

Based on the results of the simulations carried out herein, we may argue that emission which arises predominantly from just behind the outer shock is unlikely to reveal evidence of the R-T instabilities studied in this paper. However the instabilities might show up better in maps of emission from close to the interface between the shocked ambient and shocked ejecta gas. A case in point is again Tycho’s SNR. X-ray line observations by Hwang and Gotthelf (1997) have provided narrow-band images in various emission lines, such as Si, S and Fe. These reveal a shell which is generally circular, but many of the images (their Figure 4e) show a corrugated, wavy outer shell and the presence of clumps. An isolated knot is also seen in several images near the southeast bulge. The presence of Si, S and Fe lines indicates that much of the emission is coming from the reverse-shocked ejecta material (see also Hwang et al. 1998). The Si and S emission may in fact be arising from close to the contact discontinuity (DC98). If so, the rippled outer boundaries may be a result of R-T instabilities similar to those described in this paper. We have carried out simulations with parameters that may be reasonable for Tycho’s SNR. Although these simulations have not exhaustively explored the parameter space (and as such are not intended as a specific model for the remnant), they do indicate that instabilities would arise within a couple of hundred years after the explosion, and that, depending on the initial conditions, they would not yet have progressed into the non-linear phase at the age of Tycho’s SNR.

We could perhaps summarise by saying that for SNe with an exponential ejecta density profile in the free-expansion phase, the presence of instabilities near the contact discontinuity is more likely to be revealed in X-ray than in radio images. This implies a lack of correlation between the relative brightness of the radio and X-ray images, at various points around the circumference. Such a lack of correlation is seen in images of Tycho (Hwang and Gotthelf 1997).

## 6 Summary

We have carried out 2-dimensional simulations of Type Ia SNe interacting with the ambient medium, expanding on the work of DC98. An exponential density distribution was assumed for the ejecta. The simulations reveal that the decelerating contact discontinuity separating the shocked supernova ejecta from the shocked surrounding medium is unstable to Rayleigh-Taylor instabilities, as expected.

The nature and characteristics of the instabilities was examined. The instabilities tend to cover at most 50% of the interaction region. In no case is the outer shock affected by the turbulent mixing in the interior. However, if the surrounding medium is a circumstellar medium formed by mass loss from the companion star, then the reverse shock is found to be wrinkled due to the effect of the instabilities. The scale of the instabilities as they enter the non-linear regime is independent of the initial perturbations introduced into the system.

The effect of increasing resolution is to increase the small-scale structure visible in each R-T finger. A fractal-like pattern is seen for a surrounding CS medium, with multiple fingers forming an R-T clump. However the number of such clumps remains the same irrespective of the initial perturbation. If the ambient medium is of constant density, then the amount of substructure seen is greater, but no clear, well-defined fractal pattern is formed.

Our results suggest that corrugated outer shocks in SNRs cannot be attributed to instabilities of the type studied herein. We do caution however that this may be a result of the initial conditions assumed herein, and it is possible that different initial conditions, or a non-axisymmetric surrounding medium, could result in much larger growth of the R-T fingers. We feel that the instabilities are more likely to be visible in emission that arises from the reverse-shocked ejecta, such as X-ray emission from young supernova remnants. This is perhaps the case for Tycho’s SNR, and some of the clumpy structure seen in the images reported by Hwang and Gotthelf (1997) could possibly be a result of R-T instabilities.

The study of instabilities helps us to understand the multi-dimensional structure that is formed during the evolution of supernova remnants, and thereby aid in interpretation of the observational data. Further studies are required that take into account instabilities that may occur during or just after the explosion process, and determine how these affect the later evolution. Also, 3-dimensional studies often tend to show subtle differences from 2-dimensional ones, and a fuller understanding will require 3-dimensional simulations.

Acknowledgements I would like to thank Dr. Roger Chevalier and Dr. Lewis Ball for useful discussions and comments on the manuscript, Dr. Noella D’Cruz for a critical reading, and the anonymous referee for useful comments.

Figure Captions


text:
# I Energy of the DP relative to the SP core, in meV/Å per dislocation, for the 90∘ partial in Si. Four different approximations (LDA, TBTE and two different Keating parameterizations) are used to compare quadrupole and dipole boundary conditions. Cell size refers to the DP case.

Stability of the Period-Doubled Core of the 90 Partial in Silicon

In a recent Letter , Lehto and Öberg (LO) investigated the effects of strain fields on the core structure of the 90 partial dislocation in silicon, especially the influence of the choice of supercell periodic boundary conditions in theoretical simulations. Specifically, they addressed the relative stability of the traditionally accepted single-period (SP) geometry vs. our recently proposed double-period (DP) structure . Performing supercell calculations that employed the original Keating potential for Si , they reached two main conclusions: (1) that a dipole-cell configuration minimizes the overall cell strain, thus being more adequate for small-cell simulations; (2) that the relative stability of the SP and DP geometries depends on the choice of boundary conditions, with the SP and DP cores being favored for “dipole” and “quadrupole” configurations, respectively. The purpose of the present comment is not to dispute their first point, with which we agree. Rather, we wish to focus on the more important DP-versus-SP stability issue, because we believe their second conclusion to be incorrect. Below, we show that their results for the relative stability between the two structures are in disagreement with cell-size converged tight-binding total energy (TBTE) calculations, which suggest the DP core to be more stable, regardless of the choice of boundary condition. Moreover, we argue that this disagreement is due to their use of a Keating potential.

Clearly, in the limit of sufficiently large supercells, all results should be independent of the the choice of boundary conditions. Here, supercell convergence is investigated by performing TBTE and Keating-potential calculations for three different supercell sizes, using both the dipole and quadrupole cells, and two different sets of parameters for the Keating-potential calculations . In our larger cell, all dislocation separations are similar to those of the 2048-atom cells studied by LO. The first observation we can draw from the TBTE results in Table I is that, as pointed out by LO, for small cells the dipole boundary condition gives results which are closer to the converged value, while the quadrupole cell has a bias of $``$10 meV/Å in favor of the DP structure, which decreases by one order of magnitude as we approach cell-size convergence. However, Table I also shows that the Keating potential has a much stronger bias in favor of the SP structure. For the 192-atom cell, the results of LDA (from Ref. ), TETB, and the Keating potentials of Refs.  and  for $`E_{DP}E_{SP}`$ are $``$80, $``$70, $``$20, and $``$15, respectively (in meV/Å, with error bars of $`\pm 5`$ in each case). Thus, the Keating potential has a systematic bias in favor of SP of $``$50 meV/Å relative to the more accurate methods. As it happens, going to cell-size convergence shifts all the $`E_{DP}E_{SP}`$ values by roughly 20 meV/Å. As a result, the converged values for the Keating potentials are close to zero, making the variations with choice of boundary conditions (which, in absolute numbers, are similar to those seen in TBTE) take on an artificial importance. In fact, the $`E_{DP}E_{SP}`$ values are so small that even the qualitative conclusions about which structure is favored can be seen to depend on the choice of Keating parameters. On the other hand, it seems clear that the more accurate TBTE results would be quite immune from displaying a sign change in $`E_{DP}E_{SP}`$ as a result of details of the choice of cell size or boundary conditions.

Finally, LO also based their conclusions in part on a density-functional calculation for a finite cylindrical sample. However, we suggest that the systematic errors involved in choosing surface boundary conditions for such a geometry may also be severe, and that such results should not be trusted in the absence of careful tests of convergence with respect to sample size.

R. W. Nunes acknowledges financial support from CNPq-Brazil. R. W. Nunes1 and David Vanderbilt2

$`^1`$ Depto. de Física, UFMG

CP 702, Belo Horizonte, MG 30123-970, Brazil

$`^2`$ Dept. of Physics and Astronomy, Rutgers University

P.O. Box 849, Piscataway, NJ 08855-0849  N. Lehto and S. Öberg, Phys. Rev. Lett. 80, 5568 (1998).

 J. Bennetto, R. W. Nunes, and D. Vanderbilt, Phys. Rev. Lett. 79, 245 (1997).

 P. N. Keating, Phys. Rev. 145, 637 (1966).

 G.-X. Qian and D. J. Chadi, J. Vac. Sci. Technol. B 4, 1079 (86).


text:
# References

Solar Neutrino Solutions in Non-Abelian Flavor Symmetry

Morimitsu Tanimoto

Department of Physics, Niigata University, Niigata, JAPAN

( Talk at ICHEP2000, July 27-August 2, 2000, Osaka, Japan )

We have studied the large mixing angle MSW solution for the solar neutrinos in the non-abelian flavor symmetry. We predict the MNS mixing matrix taking account of the symmetry breakings.

1 LMA-MSW Solution in Non-Abelian Flavor Symmetry Recent data in S-Kam favor the large mixing angle MSW (LMA-MSW) solution. How does one get the LMA-MSW solution as well as the maximal mixing of the atmospheric neutrinos in theory? It is not easy to reproduce the nearly bi-maximal mixings with LMA-MSW solution in GUT models.

The non-abelian flavor symmetry $`S_{3L}\times S_{3R}`$ or $`O_{3L}\times O_{3R}`$ leads to the LMA-MSW solution naturally. The mass matrices are

$`M_E\left(\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right),M_\nu \left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right).`$
The orthogonal matrix diagonalizes $`M_E`$ is

$`F=\left(\begin{array}{ccc}1/\sqrt{2}& 1/\sqrt{6}& 1/\sqrt{3}\\ 1/\sqrt{2}& 1/\sqrt{6}& 1/\sqrt{3}\\ 0& 2/\sqrt{6}& 1/\sqrt{3}\end{array}\right).`$ (1)
The MNS mixing matrix is given as $`U_\nu F^T`$, and so we predict $`\mathrm{sin}^22\theta _{}=1`$ and $`\mathrm{sin}^22\theta _{\mathrm{atm}}=8/9`$ in the symmetric limit.

In this talk, we discuss masses and flavor mixings of quarks/leptons in the non-abelian flavor symmetry with the SU(5) GUT. We consider $`O(3)_\mathrm{𝟓}^{}\times O(3)_{\mathrm{𝟏𝟎}}\times Z_6`$ symmetry. Our scenario for fermion masses is

$``$ Neutrinos have degenerate masses.

$``$ Quarks/charged-leptons are massless.

$``$ Symmetry breakings give $`\mathrm{\Delta }m^2`$ and other fermion masses.

2 $`𝐎(\mathrm{𝟑})_\mathrm{𝟓}^{}\times 𝐎(\mathrm{𝟑})_{\mathrm{𝟏𝟎}}\times 𝐙_\mathrm{𝟔}`$ Symmetry Quarks and leptons belong to $`\mathrm{𝟓}^{}`$ and $`\mathrm{𝟏𝟎}`$ of the $`SU(5)`$ GUT and $`\mathrm{𝟑}`$ of the $`O(3)`$ symmetry. Higgs $`H`$ ($`\overline{H}`$) belong to $`\mathrm{𝟓}`$ ($`\mathrm{𝟓}^{}`$) of the $`SU(5)`$ and $`\mathrm{𝟏}`$ of the $`O(3)`$. Then, neutrinos have the $`O(3)_\mathrm{𝟓}^{}\times O(3)_{\mathrm{𝟏𝟎}}`$ invariant mass term

$`{\displaystyle \frac{<H>^2}{\mathrm{\Lambda }}}\nu _L\nu _L.`$ (2)
The $`Z_6`$ symmetry forbids $`\psi _{\mathrm{𝟏𝟎}}(\mathrm{𝟑})\psi _{\mathrm{𝟏𝟎}}(\mathrm{𝟑})H`$, which gives degenerate up-quark masses.

The flavor symmetry is broken explicitly by $`\mathrm{\Sigma }_\mathrm{𝟓}^{}^{(i)}(\mathrm{𝟓},\mathrm{𝟏}),\mathrm{\Sigma }_{\mathrm{𝟏𝟎}}^{(i)}(\mathrm{𝟏},\mathrm{𝟓})(i=1,2)`$, which transform as the symmetric traceless tensor 5’s of $`O(3)`$. Dimentionless breaking parameters are given as

$`\sigma _{\mathrm{𝟏𝟎},\mathbf{5}^{}}^{(1)}{\displaystyle \frac{\mathrm{\Sigma }_{\mathrm{𝟏𝟎},\mathbf{5}^{}}^{(1)}}{M_f}}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 2\end{array}\right)\delta _{\mathrm{𝟏𝟎},\mathbf{5}^{}},`$
$`\sigma _{\mathrm{𝟏𝟎},\mathbf{5}^{}}^{(2)}{\displaystyle \frac{\mathrm{\Sigma }_{\mathrm{𝟏𝟎},\mathbf{5}^{}}^{(2)}}{M_f}}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right)ϵ_{\mathrm{𝟏𝟎},\mathbf{5}^{}}.`$
Neutrinos get Majorana masses from a superpotential

$`W={\displaystyle \frac{H^2}{\mathrm{\Lambda }}}\mathrm{}(\mathrm{𝟏}+\alpha _i\sigma _5^{}^{(i)})\mathrm{},`$ (3)
which yields a diagonal neutrino mass matrix. In order to get the charged lepton masses, we introduce $`O(3)_\mathrm{𝟓}^{}`$-triplet $`\varphi _\mathrm{𝟓}^{}(\mathrm{𝟑},\mathrm{𝟏})`$ and $`O(3)_{\mathrm{𝟏𝟎}}`$-triplet $`\varphi _{\mathrm{𝟏𝟎}}(\mathrm{𝟏},\mathrm{𝟑})`$. These VEV’s are determined by the superpotential

$`W=Z_\mathrm{𝟓}^{}(\varphi _\mathrm{𝟓}^{}^23v_\mathrm{𝟓}^{}^2)+Z_{\mathrm{𝟏𝟎}}(\varphi _{\mathrm{𝟏𝟎}}^23v_{\mathrm{𝟏𝟎}}^2)`$
$`+X_\mathrm{𝟓}^{}(a_{(i)}\varphi _\mathrm{𝟓}^{}\sigma _\mathrm{𝟓}^{}^{(i)}\varphi _\mathrm{𝟓}^{})+X_{\mathrm{𝟏𝟎}}(a_{(i)}^{}\varphi _{\mathrm{𝟏𝟎}}\sigma _{\mathrm{𝟏𝟎}}^{(i)}\varphi _{\mathrm{𝟏𝟎}})`$
$`+Y_\mathrm{𝟓}^{}(b_{(i)}\varphi _\mathrm{𝟓}^{}\sigma _\mathrm{𝟓}^{}^{(i)}\varphi _\mathrm{𝟓}^{})+Y_{\mathrm{𝟏𝟎}}(b_{(i)}^{}\varphi _{\mathrm{𝟏𝟎}}\sigma _{\mathrm{𝟏𝟎}}^{(i)}\varphi _{\mathrm{𝟏𝟎}})`$
where $`Z_{\mathrm{𝟏𝟎},\mathbf{5}^{}}`$, $`X_{\mathrm{𝟏𝟎},\mathbf{5}^{}}`$, $`Y_{\mathrm{𝟏𝟎},\mathbf{5}^{}}`$ are all singlets of $`O(3)_\mathrm{𝟓}^{}\times O(3)_{\mathrm{𝟏𝟎}}`$. Minimizing the potential, we get

$`<\varphi _\mathrm{𝟓}^{}>\left(\begin{array}{c}1\\ 1\\ 1\end{array}\right)v_\mathrm{𝟓}^{},<\varphi _{\mathrm{𝟏𝟎}}>\left(\begin{array}{c}1\\ 1\\ 1\end{array}\right)v_{\mathrm{𝟏𝟎}}.`$
Masses of charged leptons arise from a superpotential

$`W={\displaystyle \frac{\kappa _E}{M_f^2}}(\overline{e}\varphi _{\mathrm{𝟏𝟎}})(\varphi _\mathrm{𝟓}^{}\mathrm{})\overline{H},`$ (4)
which is the realization of ”Democratic Mass Matrix”,

$`M_E\left({\displaystyle \frac{v_\mathrm{𝟓}^{}v_{\mathrm{𝟏𝟎}}}{M_f^2}}\right)\left(\begin{array}{ccc}1& 1& 1\end{array}\right)\left(\begin{array}{c}1\\ 1\\ 1\end{array}\right).`$ (5)
Adding the superpotential containing the flavor symmetry breaking parameters $`\sigma _{\mathrm{𝟓}^{},\mathbf{10}}^{(i)}`$, we get the charged lepton mass matrix:

$`M_E^HF^TM_EF=\kappa _E\left({\displaystyle \frac{v_\mathrm{𝟓}^{}v_{\mathrm{𝟏𝟎}}}{M_f^2}}\right)<\overline{H}>`$
$`\times \left(\begin{array}{ccc}ϵ_\mathrm{𝟓}^{}ϵ_{\mathrm{𝟏𝟎}}& ϵ_{\mathrm{𝟏𝟎}}\delta _\mathrm{𝟓}^{}& ϵ_{\mathrm{𝟏𝟎}}\\ ϵ_\mathrm{𝟓}^{}\delta _{\mathrm{𝟏𝟎}}& \delta _\mathrm{𝟓}^{}\delta _{\mathrm{𝟏𝟎}}& \delta _{\mathrm{𝟏𝟎}}\\ ϵ_\mathrm{𝟓}^{}& \delta _\mathrm{𝟓}^{}& 3\end{array}\right),`$ (6)
in which order one coefficients are omitted. The mass ratios are given as

$`{\displaystyle \frac{m_\mu }{m_\tau }}𝒪(\delta _\mathrm{𝟓}^{}\delta _{\mathrm{𝟏𝟎}}),{\displaystyle \frac{m_e}{m_\tau }}𝒪(ϵ_\mathrm{𝟓}^{}ϵ_{\mathrm{𝟏𝟎}}).`$
The quark/lepton masses and mixings fix

$`\delta _{\mathrm{𝟏𝟎}}\lambda ^2,ϵ_{\mathrm{𝟏𝟎}}\lambda ^3\lambda ^4,\delta _\mathrm{𝟓}^{}\lambda ,ϵ_\mathrm{𝟓}^{}\lambda ^2.`$
3 Neutrino Masses and Mixings Neutrino masses are given as

$`m_1`$ $`c_\mu (1+\alpha _1\delta _\mathrm{𝟓}^{}+\alpha _2ϵ_\mathrm{𝟓}^{}),`$
$`m_2`$ $`c_\mu (1+\alpha _1\delta _\mathrm{𝟓}^{}\alpha _2ϵ_\mathrm{𝟓}^{}),`$
$`m_3`$ $`c_\mu (12\alpha _1\delta _\mathrm{𝟓}^{}),c_\mu ={\displaystyle \frac{<H>^2}{\mathrm{\Lambda }}}`$
which leads to ($`\mathrm{with}\delta _\mathrm{𝟓}^{}\lambda ,ϵ_\mathrm{𝟓}^{}\lambda ^2`$)

$`\left|{\displaystyle \frac{\mathrm{\Delta }m_{21}^2}{\mathrm{\Delta }m_{32}^2}}\right|=`$ $`{\displaystyle \frac{2}{3}}{\displaystyle \frac{\alpha _2ϵ_\mathrm{𝟓}^{}}{\alpha _1\delta _\mathrm{𝟓}^{}}}{\displaystyle \frac{1+\alpha _2ϵ_\mathrm{𝟓}^{}}{1\frac{1}{2}\alpha _1\delta _\mathrm{𝟓}^{}}}\lambda ^2\lambda .`$
Putting $`\mathrm{\Delta }m_{32}^2=3\times 10^3\mathrm{eV}^2`$, we predict $`\mathrm{\Delta }m_{21}^2(\mathrm{factor})\times 10^4\mathrm{eV}^2`$, which is consistent with the LMA-MSW solution. Flavor mixings come from the charge lepton mass matrix since the neutrino one is diagonal. The charged lepton mass matrix is diagonalized by $`V_R^{}M_E^HV_L`$, in which

$`V_L^{}\left(\begin{array}{ccc}1& \lambda & \lambda ^2\\ \lambda & 1& \lambda \\ \lambda ^2& \lambda & 1\end{array}\right).`$ (7)
The neutrino mixing matrix is given by $`V_L^{}F^T`$. We predict

$`\mathrm{sin}^22\theta _{}=(1{\displaystyle \frac{4}{3}}\lambda ^2)^20.87`$
$`\mathrm{sin}^22\theta _{\mathrm{atm}}={\displaystyle \frac{8}{9}}(1\lambda ^2)(1+{\displaystyle \frac{1}{\sqrt{2}}}\lambda 2\lambda ^2)^2`$
$`0.95`$
$`\left|U_{e3}\right|={\displaystyle \frac{2}{\sqrt{6}}}\lambda (1{\displaystyle \frac{1}{\sqrt{2}}}\lambda )0.14.`$ (8)
4 Summary It is remarked that:

$``$ The solar neutrino mixing $`\mathrm{sin}^22\theta _{}`$ deviates from the maximal mixing ($`0.87`$).

$``$ The atmospheric neutrino mixing $`\mathrm{sin}^22\theta _{\mathrm{atm}}`$ deviates from $`8/9`$ depending phase of $`\lambda `$.

$``$ $`𝐔_{\mathrm{𝐞𝟑}}`$ is near to the experimental bound of CHOOZ ($`\mathbf{0.16}`$) .

Neutrino masses are degenerated within a factor 2. For example, we get $`m_10.030\mathrm{eV},m_20.033\mathrm{eV},m_30.058\mathrm{eV}`$, which is consistent with $`\beta \beta _{0\nu }`$ decay bound.