text:
# How relevant is the torus activity/geometry for the TeV gamma-rays emitted in the jet of M87 ?

## 1. Introduction

We have already discussed in Donea & Protheroe (2002) some aspects related to the IR radiation from ubiquitous dusty tori in quasars and some blazars. We have shown that the torus radiation could present a serious problem for the escape of TeV $`\gamma `$-rays from emission regions which are not well above hot tori.

It is hard to accept that M87, a mis-aligned BL Lac object, that posses a huge black hole (Ford et al., 1995), a very active jet (Wilson et al., 2002) and BLR (Sankrit et al., 1999) and possibly NLR (Zavala & Taylor, 2002), completely lacks the dusty component, be it very cold. M87 has a low luminosity nucleus and the fact that strong mid-infrared emission was not detected from the centre of M87 (Perlman et al., 2001) suggests either that there is no torus in M87, or if there is one it should have a low radiative efficiency making its detection difficult. Using argument related to the age of M87, the fact that M87 could undergo on-off activity cycles (Di Matteo et al., 2002) and the possibility that M87 could harbour a double black hole system (Zier & Biermann, 2001) we assume that M87 could still have sufficient dust around its nucleus (Fig. 1a).

We calculate the optical depth for $`\gamma `$-$`\gamma `$ pair-production $`\tau _{\gamma \gamma }`$ of GeV-TeV photons emitted by the jet interacting with IR photons from the torus (Fig. 1b). Bai and Lee (2002) have suggested that M87 could emit detectable TeV photons with an inverse-Compton emission peak at only 0.1 TeV. Protheroe et al. (2002) have predicted that $`\gamma `$-ray emission involving hadronic interactions cuts off at 10 GeV - 1 TeV, depending on the luminosity of the low hump of SED. We find that $`\gamma `$-rays with energies around 0.1 TeV are not absorbed by the IR cold photons and therefore could be detectable. Although the high temperature of $`1000`$ K is unrealistic for the present near-dormant state of M87, it is interesting to see what the difference would be in absorption when M87 would power on the torus.

We note that $`\gamma `$-rays with energies above 10 TeV (if produced) could be absorbed even when the torus is cold ($`T=100`$ K) and thin, and the emission region lies within the torus. As one moves away from the black hole, the TeV emission if produced in the jet, escapes the $`\gamma `$-$`\gamma `$ absorption and could be detected by the Cherenkov telescopes. We conclude that the absorption of GeV to 10 TeV photons by $`\gamma `$-$`\gamma `$ interactions with the infrared photons produced by a cold torus is negligible. Only a strong and rapid ignition in the nucleus of M87 could heat up the dust and increase the absorption of $`\gamma `$-rays produced nearby the infrared torus.

## References

Bai J. M., Lee M. G., 2001, ApJ, 549, L173.

Di Matteo T., et al., astro-ph 0202238

A.-C. Donea, R.J. Protheroe, 2002, Astroph. Phys, to appear.

R.J. Protheroe, A.-C. Donea, Reimer A, 2002, Astroph. Phys, submitted.

Perlman, E. S., et al., 2001, ApJ, 561, L51.

Sankrit, R., Sembach, K. R., & Canizares, C. R. 1999, ApJ, 527, 733

Zavala R. T., Taylor G. B., 2002, ApJ, 566, 9.

Zier C., Biermann P. L., 2001, A&A, 377, 23.

Wilson A.S., Yang, Y., ApJ, 568,133, 2002.


format: text

text:
# Nanoscale phase separation in manganites

## 1 Introduction

Manganites, the $`\mathrm{Mn}`$-based magnetic oxide materials such as $`\mathrm{LaMnO}_3`$, have been known for more than 50 years. Jonker and van Santen , and in more details Wollan and Koehler  investigated a rich magnetic structure of $`\mathrm{Ca}`$-doped $`\mathrm{La}_{1x}\mathrm{Ca}_x\mathrm{MnO}_3`$ and conductivity in these materials. Specifically, a strong correlation between magnetic and transport properties in manganites was observed. These materials exhibit a metal-like resistance in a ferromagnetic phase and an insulating behavior in an antiferromagnetic phase.

While the variety of magnetic structures in manganites was explained by Goodenough  based on the theory of semicovalent exchange, the correlation between transport and magnetic properties was for a first time qualitatively explained by Zener . He suggested that the conduction electrons travel in $`\mathrm{La}_{1x}\mathrm{Ca}_x\mathrm{MnO}_3`$ through the Mn4+ ions and each ion carries a fixed magnetic moment which is strongly ferromagnetically coupled with the spin of a carrier by a generalized Hund rule. Since a spin of a conduction electron should be aligned parallel to the local spin, in the classical picture a conduction electron cannot move in antiferromagnetically ordered environment. Anderson and Hasegawa  solved a problem of two local spins and one conduction electron, providing a strong mathematical background to Zener ideas. They found that if the bare hopping amplitude was $`t`$, than a hopping amplitude of an electron moving between the two local spins is $`t_{eff}=t\mathrm{cos}\theta `$, where $`\theta `$ is an angle between the directions of the local spins. Thus, the kinetic energy of conduction electrons directly depends on an angle between the sublattice moments. Moreover, the conduction electrons tend to align ferromagnetically the local spins, which surround them. Hence, a competition between ferromagnetic coupling via a conduction electron (so called double exchange mechanism) and an antiferromagnetic coupling via superexchange of two neighboring local spins takes place. De Gennes  suggested that this competition results in a homogeneous canted state, namely an angle $`\theta `$ is uniform and monotonically changes from $`\pi `$ (collinear antiferromagnetic order) to $`0`$ (collinear ferromagnetic order) with increasing carriers concentration $`x`$. Based on these considerations de Gennes plotted the first phase diagram of the manganites.

Soon after this work, Nagaev  improved de Gennes results by considering quantum fluctuations associated with the local spins. Nagaev proved that electron can move even in antiferromagnetically ordered phase with a small hopping amplitude $`t/\sqrt{2S+1}`$. Nagaev, Kasuya, and Mott also proposed  that for small electron concentrations it is more favorable for conduction electrons to form a self-trapped state (ferromagnetic polaron or ferron) in antiferromagnetic matrix by creating a small ferromagnetic bubble, rather than to form a homogeneous canted state in the whole sample. Thus, it was one of the first hints for nanoscale phase separation in a double-exchange model made by that time.

Recent growth of interest in manganites was initiated in 1993 by the discovery of the colossal magnetoresistance (CMR) effect in doped LaMnO3. The CMR phenomena implies a drastic decrease of resistivity in manganites in the presence of a magnetic field . Soon after discovery of CMR in manganites a phase diagram of $`\mathrm{La}_{1x}\mathrm{Ca}_x\mathrm{MnO}_3`$ was revised . It was found that a variety of phases appears in manganites in addition to those predicted by de Gennes.

There are plenty of experimental and theoretical studies of manganites nowadays. They were initiated first of all by the potential technological applications of colossal magnetoresistance phenomena and also by the interesting physics of strong correlations, which manifests itself in these materials. In particular, the interaction of spin, charge, and orbital degrees of freedom in manganites as well as their rich phase diagrams draw much attention of theorists and experimentalists in recent years.

The important question that has to be answered is about a leading mechanism responsible for CMR in the optimum doping region ($`x0.3`$). Some authors argue that CMR could be explained in framework of the double-exchange mechanism alone , others  claim that it is necessary to take into account a lattice interaction (Jahn-Teller polarons), some insist that percolation-type arguments could explain CMR . However, as it was pointed out by Dagotto et al and Arovas and Guinea  both analytical and numerical calculations in various models related to manganites exhibit a strong tendency toward phase separation in the wide range of temperatures and concentrations. Thus, it is believed that CMR phenomena could be understood as a competition and co-existence of different phases in manganites as well as a phenomena related to the proximity of the optimum doping region to various phase-separation thresholds. Note that at higher concentrations close to half-filling there appears another threshold of phase separation in the system corresponding again to the formation of ferromagnetic droplets, but now in a charged ordered insulating matrix .

One of the authors  demonstrated that the double-exchange model at low doping is unstable toward phase separation and the energy of a homogeneous canted state is higher than an energy of a self-trapped state corresponding to a ferromagnetic polaron. Hence a legitimate question arises whether the stability of a polaronic state is preserved when quantum fluctuations of spins and tails of wave function of conduction electrons are taken into account.

## 2 Basic theoretical model

The general chemical formula of the most popular class of manganites is $`\mathrm{Ln}_{1x}\mathrm{A}_x\mathrm{MnO}_3`$, where Ln is a trivalent cation Ln3+ (La, Pr, Nb, Sm, $`\mathrm{}`$), and A2+ is a divalent cation (Ca, Sr, Ba, Mg, $`\mathrm{}`$). Oxygen is in a O2- state, and the relative fraction of Mn4+ and Mn3+ is controlled by a chemical doping $`x`$. This class of manganites has the perovskite structure. In the cubic lattice environment, the five-fold degenerate $`3d`$-orbitals of Mn-ions are split into three lower-energy levels ($`d_{xy}`$, $`d_{yz}`$, and $`d_{zx}`$), usually referred to as $`t_{2g}`$, and two higher-energy $`e_g`$ states ($`d_{x^2y^2}`$ and $`d_{3z^2r^2}`$). The $`t_{2g}`$ levels with three electrons form a state with a local spin $`S=3/2`$, whereas delocalized $`e_g`$ states contain an electron or are unoccupied depending on the chemical doping $`x`$. The $`e_g`$ states are further split by the static Jahn \- Teller effect and for the simplicity we will treat here only the lowest $`e_g`$ state, assuming a Jahn-Teller gap to be large enough and neglecting any orbital effects in our consideration.

The simplest theoretical model suggested for the explanation of the properties of manganites is the ferromagnetic Kondo lattice model ($`sd`$ model):

$$\widehat{H}=J_H\underset{i}{}𝐒_i𝝈_it\underset{i,j}{}Pc_{i\sigma }^+c_{j\sigma }P+J_{ff}\underset{i,j}{}𝐒_i𝐒_j$$
(1)
The first term in equation (1) represents a strong on-site Hund’s ferromagnetic coupling ($`J_H>0`$) between the local spin $`S=3/2`$ and the spin of a conduction electron. In real manganites, the Hund’s interaction $`J_H`$ is of the order of 1 eV. The second term in equation (1) is the kinetic energy of the conduction electrons. The projection operator $`P`$ corresponds to the case of singly occupied $`e_g`$ orbitals (a strong Hund’s interaction prevents two conduction electrons with different spin projections from occupying the same site). Note that a strong electron-lattice interaction significantly reduces the effective width $`W`$ of the conduction band ($`W=2zt`$) resulting in a rather small hopping amplitude $`t0.3`$ eV. The third term in equation (1) is a weak antiferromagnetic coupling between local spins on neighboring sites, with $`J_{ff}0.001`$ eV. In equation (1), symbols $`i,j`$ mean the summation over $`z`$ nearest neighbor sites.

In the case of a strong on-site Hund’s coupling ($`J_HWJ_{ff}`$) the model described by the first two terms of the Hamiltonian (1) is referred to as the double exchange model. Note that if all local spins are ferromagnetically aligned, the conduction electrons will move freely in their surrounding. Thus, model (1) describes the competition between the direct antiferromagnetic coupling of local spins and the double exchange via conduction electrons, which tends to order local spins ferromagnetically. In the strong-coupling limit, Hamiltonian (1) can be simplified:

$$\widehat{H}=\underset{i,j}{}t(\theta _{ij})a_i^+a_j+J_{ff}S^2\underset{i,j}{}\mathrm{cos}(\theta _{ij})$$
(2)
where $`a_i^+`$ and $`a_j`$ are creation and annihilation operators of spinless fermions (conduction electrons whose spins $`\sigma `$ are aligned parallel to the local spins), $`t(\theta _{ij})`$ is an effective hopping amplitude, and $`\theta _{ij}`$ is an angle between sublattice moments, as we already discussed. The hopping amplitude in the case of classical spins ($`S1`$) reads:

$$t(\theta )=t\mathrm{cos}(\theta /2)$$
(3)
In Nagaev’s quantum approach the local spins at empty sites have the maximum projection, $`+S`$, on the magnetization vector of the corresponding sublattice. At occupied sites, however, local spin $`𝐒`$ and spin $`𝝈`$ of a conduction electron form a state with the total spin $`S+1/2`$, but with two possible projections of it $`S\pm 1/2`$. So there are two effective bands in quantum canting corresponding to the two different projections of the total spin. Their bandwidths read :

$$t_\pm (\theta )=\frac{t}{2S+1}\left[\sqrt{2S+1+S^2\mathrm{cos}^2(\theta /2)}\pm S\mathrm{cos}(\theta /2)\right]$$
(4)
The quantum hopping amplitude drastically differs from the classical de Gennes one. In contrast to the de Gennes picture, an electron can still move through antiferromagnetic matrix creating a state with $`S_{tot}^z=S+1/2`$ at one site and a state with $`S_{tot}^z=S1/2`$ at a neighboring site (a string-like motion introduced by Zaanen and Oleś ):

$$|S_{tot}^z=S+\frac{1}{2}|S_{tot}^z=S\frac{1}{2}|S_{tot}^z=S+\frac{1}{2}\mathrm{}$$

Hence, there are two equal hopping amplitudes in the case of electron motion through antiferromagnetic background: $`t_+=t_{}=t/\sqrt{2S+1}`$. On the other hand, for ferromagnetic ordering one gets from equation (4) $`t_+=t`$ and $`t_{}=t/(2S+1)`$. The proportionality of an effective bandwidth in the case of quantum canting to $`1/\sqrt{S}`$ is just an implication of its quantum nature.

## 3 Homogeneous canted state

An energy of the classical de Gennes state with an account taken for the antiferromagnetic interaction between the local spins reads:

$$E=ztx\mathrm{cos}(\theta /2)+\frac{1}{2}zJ_{ff}S^2\mathrm{cos}(\theta ),$$
(5)
where $`z`$ is the number of nearest neighbors, and $`x`$ is a carriers concentration. The first term in this equation is the gain in the kinetic energy, and the second term is the loss in the energy of antiferromagnetic interaction between local spins. Minimization of the energy (5) with respect to the parameter $`\mathrm{cos}(\theta /2)`$ yields:

$$E=\frac{zt^2x^2}{4J_{ff}S^2}zJ_{ff}S^2,\mathrm{cos}(\theta /2)=\frac{tx}{2J_{ff}S^2}$$
(6)
Thus, we have a transition from a collinear antiferromagnetic state for $`x=0`$ to a collinear ferromagnetic state for $`x=x_{c4}=2J_{ff}S^2/t`$. For $`0<x<x_{c4}`$, the canting angle ($`\theta \pi `$), and a homogeneous canted state takes place.

Previously, various homogeneous states were considered taking into account the quantum hopping amplitudes . In contrast to the classical case, a collinear antiferromagnetic state remains energetically favorable up to the critical value of the carrier concentration $`x_{c1}`$, which is given by:

$$x_{c1}=\frac{\pi ^4}{3}\left[\frac{8J_{ff}S^2}{zt}\frac{1}{\sqrt{2S+1}}\right]^3$$
(7)
Thus, in a quantum case, the canted state occurs for $`x>x_{c1}`$, whereas in the classical case the canted state appears for arbitrarily low doping levels. At higher doping levels ($`x>x_{c1}`$), the two-band quantum canted state arises, namely conduction electrons are in the two bands with total spin projections $`|S_{tot}^z=S+1/2`$ and $`|S_{tot}^z=S1/2`$. However at $`x>x_{c2}(27/2)x_{c1}`$ the bottom of the second band lies above the chemical potential level and a one-band state of quantum canting becomes favorable . Finally, at much higher carriers concentration ($`x>x_{c3}=4J_{ff}S^2/t\sqrt{2S+1}`$), a transition to the classical canted state of de Gennes (5) occurs. Note that for $`x>x_{c4}=2J_{ff}S^2/t`$ the canted state transforms into a collinear ferromagnetic state with the angle $`\theta =0`$.

To test the stability of the homogeneous state the electronic compressibility was calculated  according to the standard formula $`\kappa ^1=d^2E/dx^2`$. At very low concentrations $`x<x_{c1}`$ the compressibility is positive and a homogeneous collinear antiferromagnetic state corresponds to at least a local minimum of energy. However in the range of concentrations $`x_{c1}<x<x_{c4}`$ the compressibly $`\kappa <0`$ and the canted state is unstable. For example, we can calculate the compressibility of the classical de Gennes canted state and get:

$$\kappa ^1=\frac{zt^2}{2J_{ff}S^2}<0.$$

Hence compressibility of this state is negative, i.e. de Gennes classical canted state is also unstable. The negative sign of the compressibility indicates the instability of homogeneous state toward the phase separation. The simplest case of nanoscale phase separation corresponds to the formation of small ferromagnetic polarons inside an antiferromagnetic matrix. This state was considered in Ref.  and it was shown that a polaronic state is more favorable energetically than all the homogeneous states in the total range of concentrations $`0<x<x_{c4}`$. Note that magnetic polarons, in this case, correspond to the electron in the self-trapped ferromagnetic state of a finite radius inside an antiferromagnetic insulating matrix.

## 4 Polaronic state

As we already discussed, in the case of the classical hopping amplitude (3) a conduction electron may be self-trapped and form ferromagnetic droplets (magnetic polarons) inside an antiferromagnetic matrix. The simplest assumption is to consider that the boundary between the ferromagnetic region and the antiferromagnetic matrix is abrupt without an extended region of inhomogeneous canting. Then an energy of a polaronic state reads:

$$E=tx\left(z\frac{\pi ^2a^2}{R^2}\right)+\frac{1}{2}zJ_{ff}S^2\frac{4\pi }{3}x\left(\frac{R}{a}\right)^3\frac{1}{2}zJ_{ff}S^2\left[1\frac{4\pi }{3}x\left(\frac{R}{a}\right)^3\right].$$
(8)
In equation (8), $`R`$ is a radius of a polaron and $`a`$ is a lattice constant. The first term in equation (8) describes the kinetic energy gain due to the formation of a ferromagnetic region. The corrections to this term proportional to $`ta^2/R^2`$ correspond to the localization energy of a conduction electron inside a ferromagnetic droplet of radius $`R`$. The second term in equation (8) is a loss in the Heisenberg antiferromagnetic energy of local spins inside the droplet. Finally, the third term describes the energy of an antiferromagnetic interaction between local spins in a region outside the ferromagnetic polarons. The polaron radius is obtained from the condition of energy minimization $`dE/dR=0`$. So we have the following expressions for an energy and a polaron radius:

$`E_{pol}=ztx+{\displaystyle \frac{5}{3}}\pi ^2tx\left({\displaystyle \frac{2zJ_{ff}S^2}{\pi t}}\right)^{2/5}{\displaystyle \frac{1}{2}}zJ_{ff}S^2,`$ (9)
$`R_{pol}=a\left({\displaystyle \frac{\pi t}{2zJ_{ff}S^2}}\right)^{1/5}`$ (10)
Note that in this case the transition from a polaronic to a ferromagnetic state occurs when ferromagnetic polarons start to overlap. The critical concentration for ferromagnetic transition reads:

$$x_{c5}=\frac{3}{4\pi }\left(\frac{a}{R}\right)^3=\frac{3}{4\pi }\left(\frac{2zJ_{ff}S^2}{\pi t}\right)^{3/5}$$
(11)
Now, let us consider the quantum corrections to these results. In the case corresponding to the quantum hopping amplitude (4), a conduction electron can move through antiferromagnetic background with a heavy mass $`m^{}\sqrt{S}`$, so it is interesting to study for this case a problem concerning the stability of a magnetic polaron. Since the analysis of the discrete model (2) is rather complicated, we consider the continuum limit assuming that a radius of a polaron is much larger than the lattice constant $`a`$ (further on, we put $`a=1`$). A total energy (2) can now be written in the following form :

$`E={\displaystyle \left[z|\psi |^2+\psi ^{}\mathrm{}\psi \right]t(\theta )𝑑V}+zJ_{ff}S^2{\displaystyle \mathrm{cos}^2(\theta )𝑑V},`$ (12)
$`t(\theta )={\displaystyle \frac{t}{2S+1}}\left[\sqrt{2S+1+S^2\mathrm{cos}^2(\theta /2)}+S\mathrm{cos}(\theta /2)\right]`$ (13)
As one can see from equation (12), the total energy will lie between the two limiting values corresponding respectively to the motion of the conduction electron via ferromagnetic or antiferromagnetic backgrounds:

$$E_{FM}=zt<E<E_{AFM}=\frac{zt}{\sqrt{2S+1}}$$

Since the electron wave function should be normalized $`|\psi |^2𝑑V=1`$, we minimize the functional $`F=Et\beta |\psi |^2𝑑V`$ with respect to parameters $`\theta `$ and $`\psi `$, where $`\beta `$ is a Lagrange multiplier. The corresponding Euler-Lagrange equations have the following form:

$`\left[2z\psi +\mathrm{}\psi \right]t(\theta )+\mathrm{}\left[t(\theta )\psi \right]2\beta t\psi =0`$ (14)
$`\left[(z|\psi |^2+\psi ^{}\mathrm{}\psi ){\displaystyle \frac{t(\theta )}{\mathrm{cos}(\theta /2)}}2zJ_{ff}S^2\mathrm{cos}(\theta /2)\right]\mathrm{sin}(\theta /2)=0`$ (15)
We solve these two coupled differential equations by the following iterative procedure : (a) we choose a trial function for the canting angle $`\theta (r)`$; (b) we solve the first differential equation (14) to obtain an electron wave function $`\psi (r)`$; (c) using an obtained value for $`\psi (r)`$ we solve equation (15) to get a canting angle function $`\theta (r)`$; (d) we return to the step (a) until our iteration process converges.

Functions $`\psi (r)`$ and $`\theta (r)`$ obtained by the numerical solution of equations (14) and (15) are shown in figure 1 for a broad range of the values of parameter $`\alpha =t/J_{ff}S^2`$. Note that exact numerical solution which takes into account both the effects of quantum canting and tales of the wave function coincides with the classical Nagaev-Mott solution for $`\alpha \mathrm{}`$. One can see that the magnetic polaron represents a very good localized object, and the transition region from the ferromagnetic ordering ($`\theta =0`$) to an antiferromagnetic matrix ($`\theta =\pi `$) is narrow enough. Nevertheless, a polaronic state can disappear at relatively small value of the parameter $`\alpha _c75`$. Indeed, as one can see from figure 2, there is a transition from the polaronic state to a collinear antiferromagnetic state at small values of parameter $`\alpha <\alpha _c`$. For this case, the total energy of a magnetic polaronic state is equal to the bottom of the band for electron motion through the antiferromagnetic background, and as a result for $`\alpha <\alpha _c`$ an electron can move freely through the antiferromagnetic matrix. Note that to get a more precise value of $`\alpha _c`$ we should solve a variational problem for the functional $`F`$ on the discrete lattice since for small values of $`\alpha `$ a continuous approximation is not accurate enough. The work along these lines is in progress now.

## 5 Conclusion

The tendency toward phase-separation is very strong in the double-exchange model. We have shown that in the wide range of parameter $`\alpha =t/J_{ff}S^2`$ a conduction electron forms a self-trapped state. In this state, an electron is localized in the ferromagnetic droplet of finite radius embedded in the antiferromagnetic matrix. This construction seems rather natural in the de Gennes classical approximation of the double exchange, where hopping amplitude $`t_{eff}=t\mathrm{cos}(\theta /2)`$ and electron cannot move through an antiferromagnetic background since $`t_{eff}=0`$ for $`\theta =\pi `$. However, we have proved that even in the quantum case, when a conduction electron can travel slowly through antiferromagnetic matrix (since $`t_{eff}=t/\sqrt{2S+1}`$ for $`\theta =\pi `$), the polaronic state remains well defined and stable. Our approach to the one electron problem in an antiferromagnetic matrix corresponds to a very small doping region in real materials, where the concentration of charge carriers (e.g. holes in LaxCa1-xMnO3) is low. However we believe that even for higher concentrations the polaronic picture remains qualitatively correct. Moreover, at very low concentrations magnetic polarons should be localized on impurity sites  whereas at concentrations higher than the Mott threshold they are depinned from the impurities. Preliminary estimates show that the Mott threshold in manganites corresponds to $`x_M0.010.04`$, which is significantly lower than the critical concentration $`x_{c5}0.15`$ for the overlap of magnetic polarons.

Our model of polarons embedded in ferromagnetic matrix allows one also to calculate a magnetoresistance and a noise spectrum of manganites, if we suggest that an electron transport takes place via the hopping of electrons from one polaron to another, while a polaron itself is immobile. These calculations were carried out in Ref.  and found experimental support in the recent paper of Babushkina et al.

If we proceed now to the experimental confirmation of small-scale phase separated picture we should mention that there already exist many evidences in favor of nanoscale phase separation in low and moderatly doped manganites. The confirmation of an inhomogeneous state in manganites comes from the nuclear magnetic resonance experiments of Allodi et al, where the two different hyperfine lines corresponding to ferromagnetic and antiferromagnetic regions were observed. Experiments on neutron scattering of Biotteau et al support the idea of small ferromagnetic droplets embedded in an antiferromagnetic or a canted matrix. And if we turn to transport properties, a very natural picture of electron percolation, which is in agreement with our model, was experimentally confirmed by Babushkina group in Ref. .

Thus the experimental and theoretical picture strongly confirms a phase-separated state in manganites in the region of low doping. Moreover, we believe that this picture remains qualitatively correct for the concentrations optimal for the CMR effect in the high temperature region $`T>T_C`$ ($`T_C`$ is a Curie temperature), where the ferromagnetic fluctuations of the short range (the temperature polarons) are present . Hence, a combination of very intuitive picture of polarons and the ideas of the percolation theory could provide a correct description for the behavior of manganites in the wide range of temperatures and carrier concentrations.

The authors acknowledge helpful discussions with D.I. Khomskii, V.I. Marchenko, I.A. Fomin, B.E. Meierovich, I. González, M. Hennion, and E. Pchelkin. This work was supported by the Russian Foundation for Basic Research (grant 02-02-16708 and 00-15-96570), INTAS grant 01-2008, CRDF grant RP2-2355-MO-02 and Russian President Program for Science Support (grant 00-15-9694).

## References


text:
# A Search for Eclipsing Binaries in Galactic Globular Clusters

binaries: eclipsing, blue stragglers, color-magnitude diagram, dust, extinction, globular clusters: individual (NGC 3201, M10, M12)

Dissertation Summary

Electronic mail: kaspar@dtm.ciw.edu

Thesis work conducted at: Dept of Astronomy, University of Michigan

Ph.D. Thesis directed by: Mario L. Mateo ; Ph.D. Degree awarded: July 2002

Variable stars have historically served as tools and laboratories in our understanding of stellar formation and evolution, the formation of star clusters, the calibration of distance determination methods, and a variety of other areas. In particular, the study of eclipsing binaries (EBs) in globular clusters (GCs) can provide direct distance estimates to clusters as well as constraints for the turnoff masses of GC stars. Knowing precise distances to GCs would constitute an independent check of widely used distance determination methods in astronomy. Obtaining masses of GC turnoff stars provides a fundamental test of low-metallicity stellar models which calculate ages for these stars - and thus the GCs themselves - and thereby allow a lower limit estimate of the age of the universe.

My dissertation research with M. Mateo at the University of Michigan consists of a monitoring survey of 10 Galactic GCs with the aim of identifying photometrically variable EBs around or below the main-sequence turnoff (MSTO). My thesis work comprises the results of our research of the GCs NGC 3201, M10, and M12. Our observing strategy is aimed at detecting variables with periods between 0.2 to around 5 days and $`16.5<V<20`$, thereby optimizing our chance to find the valuable detached EBs in the target clusters. In addition, we correct for differential reddening variations across the cluster fields of up to several tenths of magnitudes in $`V`$ by internally creating a differential extinction map with a arcmin resolution from our cluster photometry data, and calculate the additional reddening zero point using isochrone fitting. The resulting improvement in the appearance of the color-magnitude diagrams (CMDs) of the clusters is considerable (see von Braun & Mateo 2001, AJ, 121, 1522).

We obtained approximately 200 $`VI`$ epochs with about 20000 stars per image for each of the three clusters. Analysis of these data revealed the existence of 14 variable stars (11 EBs) in the field of NGC 3201, three variables (1 EB) in the field of M10, and 2 EBs in the field of M12. Spectroscopic follow-up work showed that only one variable (a blue straggler W Ursa Majoris contact EB) in the field of NGC 3201 is associated with the cluster. Another W UMa EB is most likely a member of M12, based on its location in M12’s CMD and its empirically calculated absolute magnitude. The rest of the variable stars we detected are members of the Galactic disk (see von Braun & Mateo 2002, AJ, 123, 279; and von Braun et al. 2002, AJ, 124, 2067 for details). We thus calculate a ratio of observable, short-period, main-sequence binaries to main-sequence stars of around 1/500 for the Galactic disk (consistent with literature estimates) and 1/9000 for GC members. The latter value is considerably smaller than the commonly quoted value of 1/1000. While the discrepancy may be due to small-number statistics, we attribute it to the fact that we lose about 20% of our fields to crowding toward the cluster center where the binary fraction should be highest after a few cluster relaxation times.

In my dissertation, I show the cluster fields and CMDs (before and after dereddening) with the locations of the variable stars, our differential extinction maps, as well as the phased lightcurves and spectra (whereever applicable) of the variable stars. For additional information, or to obtain a copy of the dissertation, please contact me or visit http://www.ciw.edu/vonbraun.

I would like to express my sincere gratitude to the U of Michigan Dept of Astronomy faculty (especially M. Mateo), postdocs (past and present), and particularly my fellow graduate students (especially K. Chiboucas) for their help, advice, and support.


text:
# A weakly non-adiabatic one-zone model of stellar pulsations: application to Mira stars

## 1 Introduction

The region in the Hertzsprung–Russell diagram which provides the majority and, probably, the most interesting classes of pulsating stars is the Asymptotic Giant Branch (AGB). As the stars of low and intermediate mass (from say $`1`$ to $`11M_{}`$) evolve along the AGB phase, they experience recurrent thermal instabilities and substantial mass loss. Depending on the phase in the thermal pulse cycle, they may spend intermittent time intervals as pulsating stars (Groenewegen & Jong 1994). As a consequence, they eject freshly synthesized material into the interstellar medium. Thus, they also play a crucial role for our understanding of the chemical evolution of galaxies. Last but not least, the corresponding stellar models are basic tools to study how planetary nebulae and their central objects form — see Habbing (1996) and Willson (2000) for reviews on the subject.

A substantial number of pulsating stars in the Galaxy show irregular behaviour and have the characteristics of AGB stars (Wood 2003). This makes them an interesting research target for theoretical modeling. The work done so far on non-linear stellar pulsations can be grouped into two different categories: the full numerical hydrodynamical approach and the more straightforward approach of the theory of dynamical systems. The former is directly based on an extensive knowledge of the physical processes in the star — equation of state, opacities, and many more — and requires the implementation of a hydrodynamical code coupled with a detailed treatment of the transfer of radiation. Hence, this approach definitely provides the most detailed and accurate description of pulsations. However, within this framework it is sometimes difficult to interpret the origin and shapes of the resulting light curves as the stellar parameters are varied. On its hand, the second approach is complementary to the former in the sense that it gives a qualitative framework in which the general features of the pulsations are easily understood and, consequently, allows to develop intuitive explanations in terms of a few very basic and relatively simple physical processes (Buchler 1993).

The hydrodynamic simulations of pulsating stars succeeded in reproducing the phenomenon of period-doubling (Buchler & Kóvacs 1987; Kóvacs & Buchler 1988; Aikawa 1990) and tangent bifurcations (Buchler, Goupil & Kóvacs 1987; Aikawa 1987) found in the light curves of some specific stars by changing the surface temperature or surface gravity. In particular, it has been found that there is strong evidence of underlying low-dimensional chaos from the non-linear analysis of the light curves of the irregular stars R Sct, AC Her, and SX Her (Buchler, Kolláth & Cadmus 2001). This paved the road to construct simple models within the framework of the theory of dynamical systems which turned out to be capable of reproducing the complex behaviour. Within this approach perhaps the simplest models of stellar pulsations are the so-called one-zone models. In this family of models the star is treated as a rigid core surrounded by a homogeneous gas shell. The one-zone models are not intended to be a substitute for fully hydrodynamical and sophisticated models. Instead, they merely seek a global and qualitative explanation of what process might be at the root of the chaotic behaviour. In the present paper, we propose a one-zone model of non-linear non-adiabatic stellar oscillations. We present the parametric study of the system and draw the attention on a particular set of numerical solutions which may have implications in the study of the variability of Long Period Variables (LPVs).

## 2 The one-zone model

Baker et al. (1966) were the first to introduce the one-zone model as tool to study the nonlinear behaviour of stellar pulsations. Buchler & Regev (1982) and Auvergne & Baglin (1985) later studied one-zone models under the assumption that the nonlinearity of the adiabatic coefficient $`\mathrm{\Gamma }_1`$ is the main trigger for nonlinear pulsations, whereas Tanaka & Takeuti (1988) pointed out that dynamical stability might be necessary for realistic models of pulsating stars, which is the approach followed here.

The one-zone model approach is especially suited for the AGB stars as the density difference between the central core and the outer layers is so large that these two regions can be considered as effectively decoupled. Following Icke et al. (1992) and Munteanu et al. (2002), we consider the stellar pulsation to be simulated by a variable inner boundary located well beneath the photosphere and moving with constant frequency (piston approximation). This sinusoidal driving consists of pressure waves originating in the interior and propagating through a transition zone until they dissipate in the outer layers. We denominate the driving oscillator “the interior” and the driven oscillator “the mantle”. While the radius of the former region, $`R_0`$, can be estimated as the radius of the hydrogen burning shell, the radius of the latter, $`R`$ is determined by the dissipation of the aforementioned pressure waves.

The variation of the interior radius, $`R_\mathrm{c}`$, around the equilibrium value, $`R_0`$, is given by $`R_\mathrm{c}=R_0+\alpha R_0\mathrm{sin}\omega _\mathrm{c}\tau `$, where $`\alpha `$ and $`\omega _\mathrm{c}`$ are, respectively, the fractional amplitude and the frequency of the driving. As in Stellingwerf (1972), we use the equation of motion and the energy equation without energy sources and in absence of any driving force in order to determine the final equation of motion describing the dynamics of the mantle. In the present work, we will study the particular case in which the luminosity at the base of the mantle equals the (constant) equilibrium luminosity of the star, $`L_{}`$. For convenience, we use as variables the non-dimensional radius $`rR/R_{}`$, pressure $`pP/P_{}`$ and time $`t\omega _\mathrm{m}\tau `$, where the stellar radius and pressure were normalized to their equilibrium values, while

$$\omega _\mathrm{m}\sqrt{\frac{GM}{R_{}^3}}$$
(1)
is the characteristic frequency of the star.

Considering that the additional perturbative acceleration is proportional to the driving acceleration with a transmission coefficient $`Q`$ as in Icke et al. (1992) and following the reasoning of Saitou, Takeuti & Tanaka (1989) for the energy equation, we obtain the following equation of motion

$`{\displaystyle \frac{dr}{dt}}`$ $`=`$ $`v`$
$`{\displaystyle \frac{dv}{dt}}`$ $`=`$ $`pr^2r^2`$
$`Q\alpha \omega ^{4/3}\mathrm{sin}(\omega r\omega t\alpha \omega ^{1/3}\mathrm{sin}\omega t)`$ (2)
$`{\displaystyle \frac{dp}{dt}}`$ $`=`$ $`3\mathrm{\Gamma }_1r^1vp\xi r^3(r^\beta p^\delta 1),`$
where

$$\beta =a(r^3p1.2)+21.6,\delta =3.6r^3p(r^3p0.2)$$
(3)
are the coefficients introduced in Saitou et al. (1989) as a saturation effect of the $`\kappa `$-mechanism with $`a`$ being the control parameter. The characteristic frequency of the perturbation, $`\omega `$ results from the use of the dimensionless time unit and from the assumption that $`R_0`$ encompasses almost the entire stellar mass. It is defined as

$$\omega \frac{\omega _\mathrm{c}}{\omega _\mathrm{m}}=\left(\frac{R_0}{R_{}}\right)^{3/2}.$$
(4)
In its simplest form, the study of stellar pulsations can be considered as a thermo-mechanical, coupled oscillator problem (Gautschy & Glatzel 1990). The coupling constant is given by the ratio of the dynamical to the thermal time scale in the outer layers of the star. Whenever the thermal time scale ($`\tau _{\mathrm{th}}4\pi r^2\rho \mathrm{\Delta }rc_\mathrm{V}T/L_\mathrm{r}`$) of an outer region of the star of radial extension $`\mathrm{\Delta }r`$ happens to become comparable to the sound-traveling time through that region ($`\tau _{\mathrm{dyn}}\mathrm{\Delta }r/c_\mathrm{s}`$), non-adiabatic effects are relevant. This implies that there is an efficient exchange of mechanical and thermal energies in that region. The ratio of the time scales is much smaller than one throughout most of the envelope and is close to unity only in the outermost regions. Significant non-adiabatic effects are relevant for helium stars, very massive AGB stars, some post-AGB stars and in the ionization zones of hot stars (Stellingwerf 1986; Gautschy & Saio 1995). In the context of our system, the parameter $`\xi `$ in Eq.(2) is a measure of the non-adiabaticity and is given by the ratio of the dynamical timescale to the thermal timescale of the shell

$$\xi =\frac{L_{}}{\omega _mc_\mathrm{V}mT}.$$
(5)
The system of Eq.(2) constitutes the final set of relations for the unknowns $`r`$, $`v`$ and $`p`$. The parameters that must be specified are $`Q`$, $`\alpha `$, $`\omega `$, $`\mathrm{\Gamma }_1`$, $`a`$ and $`\xi `$. We consider an ideal gas with an adiabatic coefficient $`\mathrm{\Gamma }_1=5/3`$. The parameter fixing the evolutionary status of the star is $`\omega `$ as it provides a measure of the contrast between the core and the envelope. For AGB stars, $`R_0/R_{}`$ is of the order of 15%. This led Icke et al. (1992) to consider a value of $`\omega =20.1`$ which is the value adopted here. We have performed a parametric study involving the parameters of the perturbation, $`Q`$ and $`\alpha `$. The results of the numerical integrations show that the pair $`(Q,\alpha )`$ — and not the individual specific values of $`Q`$ and $`\alpha `$ — determines the dynamics of the system. In other words, given a fixed value of $`\alpha `$, a certain range of values of $`Q`$ can be found for which the same peculiar dynamics discussed below develops. Therefore, we present here only the case of small internal perturbation ($`\alpha 34\%`$) and amplified transmission through the envelope $`(Q>1)`$. It is worth noticing that for higher (lower) values of the parameter $`\alpha `$, the same behaviour is encountered if lower (higher) values of the coupling coefficient $`Q`$ are chosen.

## 3 Numerical results

### 3.1 Zero perturbation

The system studied by Saitou et al. (1989) can be obtained from the system given in (2) by considering the case of zero perturbation ($`Q=0`$). In spite of the inexistence of internal perturbation, this system is the prototype of non-linear self-excited oscillators due the $`\kappa `$-mechanism provided that $`\xi 0`$. The behaviour of any dynamical system $`\dot{𝒙}=𝑭(𝒙)`$, with $`𝒙=(x_1,x_2,\mathrm{},x_n)`$ and $`𝑭=(F_1,F_2,\mathrm{},F_n)`$ is critically determined by its fixed points $`𝒙_0`$ given by $`𝑭(𝒙_0)=0`$. The associated eigenvalues of the Jacobian matrix $`(J_{𝒙_0})_{ij}=(F_i/x_j)_{𝒙_0}`$ determine the nature of these fixed points. For instance, for a fixed point to be stable, it is required that all eigenvalues have negative real parts. In our case, the system has three fixed points: a trivial one, $`(r_0,v_0,p_0)=(1,0,1)`$, with mainly adiabatic origin and two other fixed points, $`(r_0^+,v_0^+,p_0^+)(0.68,0,4.75)`$ and $`(r_0^{},v_0^{},p_0^{})(8.8558,0,0002)`$, entirely due to non-adiabatic effects (that is, they exist only for $`\xi 0`$). For initial conditions close to the trivial fixed point, the period-doubling route to chaos was obtained by Saitou et al. (1989) by decreasing the surface temperature, more precisely by varying the control parameter $`a`$ in the range $`a[14,20]`$, while $`\xi `$ was kept constant. No investigation of an equivalent effect produced by varying $`\xi `$ was carried on. Concerning the dynamics near the trivial fixed point, our analysis of the dependence of the Jacobian matrix on $`a`$ and $`\xi `$ reveals that for $`a36`$, either the increase of $`\xi `$ or the decrease of $`a`$ leads to the same effect on the eigenvalues consisting in chaotic behaviour through the increase of pulsational instability. In order to prove the importance of $`\xi `$ in the dynamics of the system and to complete the study of Saitou et al. (1989), in Figure 1 we present a period-doubling route to chaos with the increase of the parameter $`\xi `$ as it appears in the space $`(r,v,p)`$. In order to ease the exploration of the dynamic details, we use mainly the Poincaré map. As we deal with a periodically driven system, the Poincaré map reduces to a stroboscopic sampling of the $`r`$, $`v`$, and $`p`$ values at multiples of $`T=2\pi /\omega `$.

### 3.2 Non-zero perturbation

We present now the properties of the non-linear oscillator that are due to the presence of the time-dependent perturbation. We fix $`a=20`$ as corresponding to the regular pulsation found by Saitou et al. (1989) and, moreover, we choose $`\xi =0.06`$ in order to compare with their work. Additionally we fix $`\omega =20.1`$ and $`\alpha 0.037`$. Our study considers the coupling coefficient $`Q`$ as the primary control parameter for the strength of the perturbation, while $`\omega `$ and $`\alpha `$ are kept constant at the values previously mentioned. In the top panels of Figure 2 we present the stroboscopic map of the system for increasing values of $`Q`$. We notice the successive creation of loops, finally leading to a knot-like structure for strong perturbation ($`Q=1.16`$). For the sake of comparison with real astronomical data, we also plot in the bottom panels of Figure 2 the light curves for the corresponding cases. The luminosity $`L`$ is obtained from the condition of radiative transfer together with the perfect gas law (Stellingwerf 1972):

$$\frac{L}{L_{}}=r^\beta p^\delta ,$$
(6)
where the normalizing constant $`L_{}`$ is the equilibrium stellar luminosity. The temporal scale is expressed in years and for this task we have used the stellar parameters associated to a typical Mira of $`1M_{}`$ as they result from the work of Vassiliadis & Wood (1993).

The most noticeable feature of the light curves of Figure 2 consists in a highly energetic sporadic burst followed by a series of smaller peaks. The number of small peaks which appear between the major ones depends on $`Q`$. Moreover, the creation of every new inward loop in the stroboscopic map is equivalent to the appearance of a new small peak in the light curve. We also found that each change in the number of loops is accompanied by a chaotic regime. That is, the transition from a regular regime to a chaotic one occurs through a sequence of period-doubling bifurcations similar to Figure 1.

For a clear exemplification of the chaotic regimes, in Figure 3 we show the light curves and the associated stroboscopic maps for a case of regular dynamics and for two cases of different degrees of irregularity. The central panels correspond to the case of regular dynamics prior to the development of the successive loops in the stroboscopic map. For a slightly lower value of the parameter $`Q`$ (top panels), a completely irregular light curve results and the stroboscopic map clearly shows it. Finally, in the lower panels, we illustrate the dynamics of the system when $`Q`$ takes a value corresponding to the accumulation of a sequence of period-doubling bifurcations.

Another important feature of the dynamics is the fact that the time interval $`\tau _\mathrm{B}`$ between major bursts increases with the strength of the perturbation, as it can be seen in the bottom panels of Figure 2. This is illustrated in Figure 4a, where the time interval between the successive major bursts is shown as a function of the coupling coefficient $`Q`$. As a visual guide we also plot the shape of the stroboscopic map at the fixed values of $`Q`$ where an additional inner loop appears. As it can be seen in this panel and in panel 4b, $`\tau _\mathrm{B}`$ significantly increases for $`Q[1.32,1.39]`$ whereas for $`Q>1.39`$ (Figure 4c) the separation between successive bursts increases drammatically. Also, in Figure 4b we show that the above mentioned increase of $`\tau _\mathrm{B}`$ can also be obtained by decreasing the non-adiabaticity of the system (that is, decreasing $`\xi `$). Nevertheless the main parameter for tuning the time interval between major bursts turns out to be $`Q`$. Note as well that these major bursts favor mass loss at exceptionally high rates and, moreover, the time intervals between them are long. Hence, it is tantalizing to directly connect them with the periodicities observed in the circumstellar shells that can be found surrounding some planetary nebulae (Van Horn et al. 2002).

### 3.3 Mathematical interpretation of the results

In order to validate the numerical results discussed above, we briefly present the mathematical characteristics of our system. Given that our system is non-autonomous — that is, the Hamiltonian is explicitly time-dependent — the typical methods of analysis of the theory dynamic systems cannot be used. To overcome this drawback, we used an averaging method (Sanders & Verhulst 1985) to transform our system into an autonomous one. The high value of the characteristic frequency ($`\omega 1`$) assures us that this method is applicable to our case. In the time-averaged framework, the time-dependent perturbation from Eq. (2) becomes $`F(r)=Q\alpha \omega ^{4/3}A\mathrm{sin}\omega r`$, where $`A=0.04993`$ is a constant evaluated numerically. The fixed points of the averaged system in the case of nonzero perturbation must satisfy the conditions:

$`p_0=r_0^4+Q\alpha \omega ^{4/3}Ar_0^2\mathrm{sin}\omega r_0`$ (7)
$`r_0^\beta p_0^\delta 1=0.`$ (8)
or, equivalently,

$$G(r)r^\beta (r^4+Q\alpha \omega ^{4/3}r^2A\mathrm{sin}\omega r)^\delta 1=0.$$
(9)
The roots of Eq. (9) are to be found numerically. The new fixed points of interest are $`r_00.65`$ and $`r_01.03`$ and represent only slight displacements of the fixed points mentioned before in the case of zero perturbation. Moreover, their associated eigenvalues maintain the form from the previous case (unstable fixed points of saddle-focus type). Increasing the parameter $`Q`$, a new fixed point is created at about $`Q_01.3855`$ at $`r0.86`$. This new fixed point is stable as the associated eigenvalues have all negative real parts. With the creation of this fixed point, the looping behaviour dissappears. For higher values of $`Q`$, the fixed point is replaced by two unstable fixed points again of saddle-focus type. The distance between them increases with $`(QQ_0)^{1/2}`$. Note, however, that these very high values of $`Q`$ result in transmissions through the mantle which are not very realistic and, moreover, the time intervals between successive bursts for larger values of $`Q`$ become very large.

### 3.4 The role of $`\omega `$

Before interpreting the behaviour of the system in the context of real data of stellar variability, we consider worth of interest to explore the role of the parameter $`\omega `$ and to justify the particular value we have used throughout this work. Icke et al. (1992) concluded that in the case of complete adiabaticity ($`\xi =0`$), a decrease of $`\omega `$ leads to stronger chaotic pulsations. The values of $`\omega `$ used in their work were equivalent to adopting stellar models in the family of low-mass stars ($`M58M_{}`$) reaching the AGB. In a recent publication (Munteanu et al., 2002), we extended their conclusion to intermediate-mass stars ($`8M_{}M11M_{}`$) also in the AGB phase, more precisely to values of $`\omega `$ around 3. In the previous sections, we have shown that in the case of $`\omega =20.1`$ a peculiar behaviour is born from the interplay between non-adiabaticity and internal perturbation. Our analysis of the system corresponding to the parametric interval $`5\omega 25`$ has revealed that such a behaviour is found for values of $`\omega `$ close to 20, that is for low-mass stars. Mathematically, we attribute this fact entirely to the creation of new fixed points mentioned in the previous section which critically alter the dynamics of the system. They exist for values of $`\omega `$ higher than about 18. For values slightly lower than 18, the dynamics resembles the one encountered for $`\omega =20.1`$ (middle panels of Figure 5), but it does not present the successive creation of new loops. Instead, the change of the control parameter $`Q`$ leads to a mixture of chaotic regimes and uncorrelated creation and dissapearance of new loops. For completeness, we present in Figure 5 the light curves and the stroboscopic maps for three values of $`\omega `$. For each case, the temporal scaling factor was computed according to the work of Vassiliadis & Wood (1993) which provides a complete set of stellar parameters for AGB stars with initial masses in the range $`0.89M/M_{}5.0`$. Due to the richness of the dynamics in the case of $`\omega =20.1`$, we use it for a comparison with observations.

## 4 Comparison with observations

The values of the parameters used throughout the numerical integrations were intended to locate the stellar models we are dealing with in the family of Long Period Variables (LPVs) and, specifically, in the families of semiregular and Mira variables. Perhaps the most important of the Mira stars is $`o`$ Ceti. Its light curve shows a peculiar variability consisting in an exceptional peak occurring every two, three or five “cycles” (Barthés & Mattei 1997). Our simple model naturally recovers this behaviour by tuning the strength of the perturbation or the coupling coefficient. Moreover, we have shown that within our simple model large peaks in the light curve are associated to large inter-pulse intervals.

The Mira stars belonging to the Large Magellanic Cloud constitute the best sample of Miras concerning both periodicities and luminosities. In Figure 6a we show the observational data for some Mira stars in the LMC (Feast 1989) and the best fit to the observational data. Among them, the ones having periods longer than about 400 days clearly appear to be over-luminous with respect to the period-luminosity relationship found for Miras with relatively short periods (Zijlstra et al. 1996; Bedding et al. 1998). With this in mind we have obtained a period-luminosity relationship from our theoretical models using the time interval $`\tau _\mathrm{B}`$ and averaged values of the major peaks in luminosity. More precisely, we have varied the parameter $`Q`$ in the range $`[1,1.26]`$ resulting in light curves whose major peaks have periodicities (in days) within the range $`2.4\mathrm{log}\tau _\mathrm{B}3.0`$. The corresponding bolometric magnitudes were computed using a reference value for the equilibrium stellar luminosity of $`\mathrm{log}(L_{}/L_{})=3.5`$ (Vassiliadis & Wood 1993) which is typical for Mira stars. Our theoretical period-luminosity relationship is shown in Figure 6b. For the sake of comparison we also show in this panel all those Miras with $`\mathrm{log}P>2.4`$ and $`M_{\mathrm{bol}}<5.0`$, that is all the stars which are found to be over-luminous in Figure 6a. We have also included two other interesting Miras: R Hya and V Hya. As it can be seen our theoretical period-luminosity relationship fits very well the observational data. Moreover, note as well that the observations tend to cluster around fixed regions of the period-luminosity relationship. These regions are coincident with the regions where we find regular oscillations with a fixed number of loops in the stroboscopic map and its sequence of period-doubling cascade leading to chaos. The gaps between the theoretical data are the consequence of the drastic change in the characteristics of our light curves when a new loop is created, and correspond to a very small change of the coupling coefficient. Hence, according to our analysis we should not find almost any star in these regions, which is exactly what it is found.

Another peculiarity of Mira stars is that some of them show alternating deep and shallow minima, giving the appearance of double maxima. Some examples are R Cen, R Nor, U CMi, RZ Cyg and RU Cyg — see Hawkins (2001) and references therein. Among these, R Cen has the most persistent and stable double maxima in the light curve while for the rest of the cases the second maximum is often weak and the light curve sometimes reverts to that of a normal Mira. Our simplistic model provides such light curves for small values of the parameter $`Q`$. As for the general “chaotic connection” for the Mira variables, the first (and unique) case of evidence of chaotic pulsation in a Mira star (R Cyg) comes from the study by Kiss & Szatmáry (2002). They associate the long sub-segments of alternating maxima in R Cyg to a period-doubling event, supporting therefore the well-known scenario of period-doubling to chaos, which we also find in our model. Buchler et al. (2001) present an overview of observational examples of chaotic behaviour in some semiregular variables (SX Her, R UMi, RS Cyg, and V CVn). They argue that AGB stars are prone to chaotic pulsations due to the fact that relative growth rates of the lowest frequency modes are of the order of unity. Higher relative growth rates are a consequence of higher luminosity–mass ratios, that is more non-adiabatic stars. Hence, not only semiregular variables, but also Mira stars should also be candidates for chaotic pulsators. In spite of its simplicity, our toy-model presents chaotic pulsations for certain intervals of the parameters characterizing the strength of the internal driving and, thus, could provide some support to the conjecture that the evolution of semiregulars and of Mira stars is strongly connected.

## 5 Conclusions and caveats

We have introduced a weakly non-adiabatic one-zone model driven by sinusoidal pressure waves intended to reproduce the irregular pulsations of Mira-like variables for which other simple models already exist (Buchler & Regev 1982; Auvergne & Baglin 1985; Reid & Goldston 2002). Our model extends the works of Icke et al. (1992) and Saitou et al. (1989). In particular, Icke et al. (1992) proposed an adiabatic model driven by pressure waves (the piston approximation) whereas Saitou et al. (1989) studied a simple non-adiabatic model without driving (the self-excited pulsation model). Our approach is justified by the large density contrast between the interior and the outer layers. One-zone models have limitations since some Mira stars are suspected to have at least two frequencies (Mantegazza 1996) which may vary independently, suggesting that more than one mode is involved. Clearly, one-zone models are unable to reproduce this behaviour. Nevertheless, the majority of Miras do not show this behaviour. Another interesting approach would have been to use the modal coupling. However, although this approach has been used to model the pulsation of Miras (Buchler & Goupil 1988), it is more appropriate for classical Cepheids, RR Lyrae and W Vir stars — see Buchler et al. (1993) and references therein. It is also worth noting here that the piston approximation — first introduced by Bowen (1988) — used in this paper has been used since then by several authors (Fleischer et al. 1995; Höfner et al. 2003) even if it has never been deeply scrutinized for validity. Nevertheless this approximation appears to correctly reproduce the velocities and mass loss rates typical of AGB stars. Hence, this approximation can be regarded as a reasonable first-order approximation of the dynamical effect of the pulsation on the atmosphere. To summarize, our model should be considered as a toy model that qualitatively reproduces the general features of the light curves.

We have thoroughly explored the interesting particularities of the system from both the astrophysical and from the mathematical points of view. We have found that the degree of non-adiabaticity turns out to be a determining factor in the development of the period-doubling route to chaos. As far as the transition to chaos is concerned, an increase in the degree of non-adiabaticity plays the same role as a decrease of the effective temperature in the model studied by Saitou et al. (1989). Moreover, due to the time-dependent perturbation, a knot-like structure is created in the phase space when varying the transmission coefficient ($`Q`$) while the strength of the perturbation is kept fixed. The resulting periodic light curves are characterized by a repetitive pattern consisting in a major peak followed by $`n`$ minor peaks, where $`n`$ is given by the number of loops in the phase space. Our theoretical light curves are in qualitative agreement with those of several well-known Mira stars. In particular the prototype of Mira stars, $`o`$ Ceti, shows this pattern of alternating major peaks followed by several minor peaks. Furthermore, for a given choice of our input parameters the resulting light curves also resemble those of some peculiar Miras (R Cen, R Nor, U CMi, RZ Cyg, or RU Cyg) which appear to have double maxima due to alternating deep and shallow minima. We have found as well that our dynamical system presents both chaotic regimes and patterns of periodicity. The chaotic regions occur for ranges of $`Q`$ between those corresponding to the creation of a new luminosity peak between major bursts. We have also noticed that for increasing strengths of the perturbation the time interval between major bursts increases. This interval increases with $`Q`$ as well. At high values of the strength of the perturbation within the range yielding the knot-like structure, we obtained a peculiar behaviour: the creation of new loops stops and the resulting light curves show a time interval during which the luminosity preceeding every major peak remains constant.

We have also obtained a theoretical period-luminosity relationship and compared it with the observational data of Miras in the LMC. In particular, we have focused on those peculiar Miras with long periods which are known to be over-luminous with respect to the best fit of Feast (1989). We have found that our model provides a reasonable fit to the period-luminosity relationship of these stars and also to the observed clustering of the over-luminous stars around certain regions in the period-luminosity diagram. The ultimate reason for this fact is closely related to the creation of a new loop in the stroboscopic map, and, consequently, of a new luminosity peak in the corresponding time series.

Finally, we would like to stress that although our simple model succeeds in furnishing reasonable comparisons with real data, a detailed quantitative interpretation is beyond its capabilities because of its crude physical assumptions. In particular, next steps towards improving the one-zone approach could eventually include the treatment of convection which is supposed to play an important role in the envelopes of these stars. This improvement will affect also the functional form of the luminosity which was rather simplified in the present work.

Acknowledgements. This work has been supported by the MCYT grant AYA2000–1785, by the MCYT/DAAD grant HA2000–0038 and by the CIRIT grants 1995SGR-0602 and 2000ACES-00017. We would also like to acknowledge our anonymous referee for valuable criticisms and suggestions.


text:
# Insulator, conductor and commensurability: a topological approach

(December 15th, 2002)

## Abstract

I discuss a topological relation of the conduction property of a many-particle system on a periodic lattice at zero temperature to the energy spectrum. When the particle number per unit cell is an irreducible fraction $`p/q`$, an insulator must have $`q`$ low-lying states of energy $`o(1/L)`$ in one dimension and of energy $`o(1)`$ in two dimensions, where $`L`$ is the linear system size.

An important question in the condensed matter physics is that under what conditions the system becomes a conductor. While we do not know, of course, the complete answer, there is a characteristic of the system which may be naturally related to the conduction. In a many-particle system on the periodic lattice, the known insulating states such as band and Mott insulators can be understood as a “locking” of the particles by the lattice. It requires the particle density to be commensurate with the lattice, namely the filling factor (the number of particles per unit cell of the lattice) to be an integer. When the filling factor is not integral but is a (simple) rational number, the particles may still be locked by spontaneously breaking the translation symmetry and thus effectively enlarging the unit cell. When a small number of particles are “doped” into the system, they become mobile carriers, and we expect the system to be gapless and conducting. While this argument looks natural, it is far from trivial in the strongly interacting systems. It is needed to clarify how far this physical expectation holds, and what are the precise consequences.

In the recent several years, it has been pointed out that the commensurability (in the sense of the particle number per unit cell) is related to the topological nature of the system. Extending the ideas by Lieb-Schultz-Mattis LSM  (LSM) and by Laughlin Laughlin , several statements on the energy spectrum related to the commensurability are derived OYA ; comme . On the other hand, the naive expectation of the spontaneous breaking of the translation symmetry apparently contradicts with the “Resonating-Valence-Bond” (RVB) type gapped quantum liquid states. They have $`1/2`$ particle per unit cell, but there is no long-range order which can be defined in terms of the local operator. (For recent examples, see e.g. Refs. MS-QDM ; MSP-QDM ; BFG .) However, this is still consistent Bonesteel ; Misg:degeneracy  with the topological argument as there is a “topological degeneracy” of the ground states under the periodic boundary condition. Moreover, the same underlying mechanism was used in a topological proof YOA ; Luttad  of the Luttinger’s theorem Luttinger . It was also applied to the Fermi sea volume of a Kondo lattice, which had been unamenable to rigorous analyses.

On the other hand, there still remain several fundamental questions, including (but not limited to) the followings. While the intuitive argument would relate the commensurability to the conduction property, the relation has not been discussed in the topological context. Moreover, argument in Ref. comme  for systems in more than one dimension relies on a rather nontrivial assumption that a finite gap does not collapse due to an Aharonov-Bohm (AB) flux AB . It is still necessary to clarify whether (or when) that assumption can be justified. In the present Letter, I shall address these questions by discussing the Drude weight at zero temperature in view of the topological argument.

The classification of the insulator and conductor has been discussed by several authors. Drude weight represents the “free acceleration” of the current by the electric field, and thus gives a criterion to define the (perfect) conductor. The Drude weight $`D`$ is defined by

$$\sigma (\omega )=i\frac{D}{\pi }\frac{1}{\omega +i\delta }+\sigma _{reg}(\omega ),$$
(1)
where $`\sigma (\omega )\sigma (\omega ,\stackrel{}{q}=0)`$ is the zero momentum component of the complex AC conductivity, $`\sigma _{reg}(\omega )`$ is assumed to be finite at $`\omega 0`$, and $`\delta `$ is a positive infinitesimal. Kohn related the dependence of the ground state energy on the AB flux to the Drude weight Kohn . His argument has been further developed by several authors, in particular by Scalapino, White and Zhang SWZ . In addition, Resta and Sorella RestaSorella  put forward an alternative but deeply related criterion of an insulating state which I will comment later.

Since the analyses on the Drude weight in Refs. Kohn ; SWZ  and the topological argument on the commensurability are both associated with the insertion of the AB flux, one may expect they can be related to each other. In order to establish a concrete connection, let me first reformulate the Drude weight in the real time. I consider a system defined on a finite $`d`$-dimensional periodic lattice of volume $`V=L_x\times L_y\times \mathrm{}\times L_d`$ with periodic boundary conditions, where $`L_j`$ denote the lengths in each directions. The lengths and the volume are measured in terms of the unit cell of the lattice. Let the filling factor (particle number per unit cell of the lattice) be $`\nu =p/q`$. I choose $`L_x`$ to be an integral multiple of $`q`$, and all other $`L_j`$’s to be mutually prime with $`q`$. In taking the thermodynamic limit $`L_j\mathrm{}`$ eventually, $`\nu `$ is kept constant.

Suppose, for the time $`t<0`$, no electric field is applied and the system was in the ground state. It can be chosen so that it is an eigenstate of the $`x`$ component of the total momentum operator $`P_x`$. For $`0<t<T`$ a small constant electric field $`E`$ is applied uniformly throughout the system. I then define the “time-dependent Drude weight”

$$D(t)=\frac{j(t)}{Et}$$
(2)
where $`j(t)`$ represents the current in $`x`$-direction. Since the Drude weight represents the free acceleration of the current by the electric field, $`D(t)`$ converges to $`D`$ in the limit of $`t\mathrm{}`$. Actually it can be proved, provided that the linear response theory holds. For a later convenience, I further define a weighted average $`\overline{D}`$ of $`D(t)`$ over the time domain $`0<t<T`$ by

$$\overline{D}=\frac{2}{T^2}_0^T𝑑ttD(t).$$
(3)
In the limit of large $`T`$, $`\overline{D}`$ is dominated by the long-time asymptotic of $`D(t)`$, and thus also should converge to $`D`$.

One still has to be careful in using the present formulation, however. For a given lattice model and a given (finite) electric field $`E`$, $`j(t)`$ cannot keep the $`t`$-linear behavior indefinitely. In particular, in a lattice model, $`j(t)`$ would not exceed some finite upper limit. This is a manifestation of the non-linear effects. Generally, when $`j(t)`$ becomes large, non-linear effects are important and the linear response theory breaks down. In order to apply the linear response theory (so that eqs. (2) and (3) gives the standard Drude weight), the current $`j(t)`$ should be sufficiently small.

Let us represent the electric field by the time-dependent vector potential $`A(t)`$, which is uniform throughout the system, in the $`x`$ direction. It is given by $`A(t)=Et`$. In fact, this $`A(t)`$ is related to the AB flux $`\mathrm{\Phi }(t)`$ which pierces through the “hole” of the torus. $`\mathrm{\Phi }(t)=L_xA(t)`$. The current $`j(t)`$ is given by

$$j(t)=\frac{1}{V}\frac{}{A}=\frac{L_x}{V}\frac{(\mathrm{\Phi })}{\mathrm{\Phi }},$$
(4)
where $`(\mathrm{\Phi })`$ is the Hamiltonian with the AB flux $`\mathrm{\Phi }`$. Since $`\mathrm{\Phi }`$ is the function of time, Hamiltonian is also a function of time.

$$\frac{d}{dt}=\frac{d\mathrm{\Phi }}{dt}\frac{(\mathrm{\Phi })}{\mathrm{\Phi }}=EL_x\frac{(\mathrm{\Phi })}{\mathrm{\Phi }}.$$
(5)
Combining eqs. (4) and (5), we obtain the interesting relation $`j(t)=(1/EV)d/dt`$. Now the averaged Drude weight $`\overline{D}`$ is given by

$`\overline{D}`$ $`=`$ $`{\displaystyle \frac{2}{T^2E^2V}}{\displaystyle _0^T}{\displaystyle \frac{d}{dt}}𝑑t`$ (6)
$`=`$ $`{\displaystyle \frac{2L_x^2}{\mathrm{\Phi }(T)^2V}}\left(_{t=T}_{t=0}\right).`$
In the Schrödinger picture, $`_{t=T}`$ should be understood as $`\mathrm{\Psi }(T)|(\mathrm{\Phi }(T))|\mathrm{\Psi }(T)`$ where $`|\mathrm{\Psi }(t)`$ is the state vector at time $`t`$. This formula is exact and does not rely on the adiabatic theorem.

If one takes $`T`$ to infinity and $`\mathrm{\Phi }(T)`$ to infinitesimal, the time evolution may become adiabatic (however, see discussions below). In this case, the above formula reduces to the well known Kohn formula Kohn

$$D=\overline{D}=\frac{L_x^2}{V}\frac{^2E_0(\mathrm{\Phi })}{\mathrm{\Phi }^2},$$
(7)
where $`E_0(\mathrm{\Phi })`$ is the ground state energy as a function of the AB flux ($`\mathrm{\Phi }`$).

Scalapino, White and Zhang SWZ  pointed out a subtlety in the Kohn formula in particular for a gapless system in $`d>1`$. Namely, many level crossing occur even for a very small AB flux $`\mathrm{\Phi }`$ in a large system, and the definition of $`^2E(\mathrm{\Phi })/\mathrm{\Phi }^2`$ becomes actually rather subtle. They offer some definitions which would give the correct Drude weight. An alternative solution proposed here is to use the real-time formulation (6) instead of the Kohn formula (7). I fix the time interval $`T`$ (and thus the “sweeping rate” of $`\mathrm{\Phi }`$) and take the thermodynamic limit $`L_j\mathrm{}`$ first. After taking the thermodynamic limit, I let $`T\mathrm{}`$. There could be more and more level crossings and the time evolution is not necessarily adiabatic, as the system size is increased with a fixed $`T`$. With huge numbers of crossings in a large system, the final state $`|\mathrm{\Psi }(T)`$ may actually be very complicated LZS  but $`\overline{D}`$ is still well-defined. By taking $`T\mathrm{}`$ limit after taking the thermodynamic limit, $`\overline{D}`$ should converge to the standard Drude weight $`D`$. This procedure in the real time corresponds to, in the frequency space, taking $`\omega 0`$ limit after taking the thermodynamic limit with $`\stackrel{}{q}=0`$.

Now let me choose $`\mathrm{\Phi }(T)`$ to be the unit flux quantum $`\mathrm{\Phi }_0=2\pi /e^2`$ where $`e`$ is the charge of the particle and I set $`\mathrm{}=c=1`$. I note that, this choice corresponds to applying a small electric field $`E=O(1/(TL_x))`$. Even at $`t=T`$, the current is $`j(T)DET=O(1/L_x)`$. Thus, it can be arbitrarily small as one takes the limit of the large system size $`L_x`$. Therefore the linear response theory should be valid in a sufficiently large system. I now invoke the argument used in Refs. AligiaOrtiz ; comme ; Luttad . Namely, as the Hamiltonian is translationally invariant for any value of $`\mathrm{\Phi }`$, the total momentum $`P_x`$ of the system does not change during the process. It can be most simply seen in the Heisenberg equation of motion $`dP_x/dt=i[,P_x]=0`$. However, at $`t=T`$, the Hamiltonian $`(\mathrm{\Phi }_0)`$ with a unit flux quantum inserted is not identical to the original one $`(0)`$, although they must have the same spectrum. The large gauge transformation

$$Ue^{2\pi i_\stackrel{}{r}xn_\stackrel{}{r}/L_x},$$
(8)
where $`n_\stackrel{}{r}`$ is the particle density operator at site $`\stackrel{}{r}`$, brings the Hamiltonian $`(\mathrm{\Phi }_0)`$ back to the original form $`(0)`$. In the course, the total momentum $`P_x`$ of the system is also changed by $`\mathrm{\Delta }P_x=2\pi \nu L_y\times \mathrm{}L_d`$ because of

$$U^1T_xU=T_xe^{i\mathrm{\Delta }P_x},$$
(9)
where $`T_x=e^{iP_x}`$ is the lattice translation operator to $`x`$ direction. The momentum change $`\mathrm{\Delta }P_x`$ is different from $`0`$ for an incommensurate filling ($`\nu `$ not an integer). I note that this argument does not require the process to be adiabatic. Even for a non-adiabatic time evolution and a complicated final state $`|\mathrm{\Psi }(T)`$, the momentum is still exactly conserved before making the large gauge transformation, thanks to the translation invariance of the Hamiltonian.

Let me denote $`E_{\mathrm{\Delta }P_x}`$ be the energy of the lowest excited state with total momentum $`\mathrm{\Delta }P_x`$. Then $`\overline{D}`$ is bounded from below as $`\overline{D}e^2L_x^2E_{\mathrm{\Delta }P_x}/(4\pi V)`$. Considering the case $`\mathrm{\Phi }(T)=2\pi n/e^2`$, we can obtain similar bounds on $`E_{n\mathrm{\Delta }P_x}`$ for $`n=1,2,\mathrm{}q1`$.

This expression suggests the nontrivial dependence of the transport property and/or the energy spectrum on the aspect ratio SWZ . In taking the thermodynamic limit, here I want to keep an isotropic limit. Namely all the length $`L_j`$ are taken to be of same order of magnitude $`L`$ and then $`L`$ is sent to the infinity. In that limit, we obtain

$$\overline{D}=D\text{const.}L^{2d}E_{n\mathrm{\Delta }P_x},$$
(10)
where $`d`$ is the dimensionality of the system. This is the main result of the present letter. Let me discuss its consequences for the systems in different dimensions $`d`$ separately.

$`𝐝=1`$

If the system has an incommensurate filling and still an insulator ($`D=0`$),

$$E_{n\mathrm{\Delta }P_x}=o(\frac{1}{L}),$$
(11)
where $`n=1,2,\mathrm{}q1`$ and $`o(1/L)`$ means vanishing faster than $`1/L`$ as $`L\mathrm{}`$. It means that, if there is a gap between the ground state(s) and the excited states, there are $`q`$-fold (nearly) degenerate ground states $`|\mathrm{\Psi }_0^{(k)}`$ with different momenta $`P_0+2\pi k\nu `$. ($`k=0,2,\mathrm{},q1`$). Intuitively, the system could acquire a gap by spontaneously breaking the translation symmetry to become an insulator. Then generally we expect $`q`$-fold degenerate ground states below the gap. This is consistent with the above bound. On the other hand, the result (11) does not rule out the possibility of a gapless incommensurate insulator. Still, the bound (11) requires that its finite-size spectrum of such a system should not scale (entirely) as $`1/L`$, although it is the case when the system is described by a conformal field theory.

On the other hand, Resta initiated the argument RestaSorella ; AligiaOrtiz  that the phase of the ground-state expectation value of eq. (8)

$$z=\mathrm{\Psi }_0|U|\mathrm{\Psi }_0,$$
(12)
can be interpreted as the polarization which is well defined only in an insulator. Namely, the system is a conductor if $`z=0`$ and an insulator if $`z0`$ in the thermodynamic limit. That this approach is rather closely related to the LSM argument is pointed out only recently by Nakamura and Voit NakamuraVoit . Because of the LSM argument, $`U|\mathrm{\Psi }_0`$ belongs to a different momentum eigenvalue than the ground state and thus $`z=0`$. This naively implies that the system is always a conductor, which is of course not true. Observing this, Aligia and Ortiz AligiaOrtiz  proposed to use $`z^{(q)}=\mathrm{\Psi }_0|U^q|\mathrm{\Psi }_0`$ instead of $`z`$, for a fractional filling $`\nu =p/q`$. This indeed fits with the picture of spontaneously broken translation symmetry, which implies the expectation value of the position operator to be multi-valued. In such a case, the physical (pure) ground state in the thermodynamic limit should be taken as $`|\mathrm{\Phi }^{(j)}_ke^{2\pi jk\nu }|\mathrm{\Psi }_0^{(k)}`$. The original order parameter $`z`$ can be non-vanishing with respect to $`|\mathrm{\Phi }^{(j)}`$. More analyses would be needed in this direction to examine the implications.

In any case, the present result implies that it would be more ’difficult’ to make an insulator out of a system with highly incommensurate filling (i.e. large $`q`$) as there are large number of restrictions. Even when the system is conducting, requiring the Drude weight to be finite we obtain $`E_{n\mathrm{\Delta }P_x}=O(1/L)`$ if the filling is incommensurate. While this is nontrivial, the standard LSM argument LSM ; OYA ; YOA  already gave this bound irrespective of the conduction property. On the other hand, eq. (11) is a stronger result which applies only to an insulator. The bound (11) certainly should not apply to a conductor, as there are many counterexamples.

$`𝐝=2`$

If the system has an incommensurate filling and still an insulator ($`D=0`$),

$$E_{n\mathrm{\Delta }P_x}=o(1).$$
(13)
Namely, they should vanish (at some rate, possibly very slowly) in the thermodynamic limit. In particular, when there is a finite gap between the ground states and the excited states, the ground states must be $`q`$-fold degenerate. On the other hand, a gapless insulator phase is also consistent with the above bound. This conclusion is essentially the same as obtained in Ref. comme  assuming that the gap does not close due to the AB flux. Here, that (questionable) hypothesis is replaced by the assumption that the system is an insulator. For a conductor, requiring the Drude weight to be finite only yields a trivial bound $`E_{n\mathrm{\Delta }P_x}=O(1)`$, Unfortunately, no useful result can be obtained for $`d3`$ with the present argument.

Straightforward generalization to the multi-species cases is possible in parallel to Ref. YOA . For example, let me define $`\nu _{}`$ and $`\nu _{}`$ respectively as up-spin and down-spin electron density, and $`\nu =\nu _{}+\nu _{}`$. One can consider the “real” electric field which couples equally to electrons with either spin, to derive similar constraints if $`\nu `$ is fractional. For the half-filled Hubbard model, this does not give any constraint as $`\nu =1`$ is an integer. However, in this case one can still consider a fictitious electric field which couples to up-spin electrons only. Then applying the same argument as before, one can deduce the bounds (11) and (13) respectively for one and two dimensions, if the system is “insulating” with respect to the up-spin electrons. For the standard Hubbard model in one dimension at half-filling, the exact Bethe Ansatz solution shows that all the low-lying states scale as $`1/L`$ for any value of the Hubbard interaction $`U`$. This is actually consistent with the above result, because for $`U>0`$ the system is insulating with respect to the (real) charge but is a “conductor” in terms of spin. For $`U<0`$, it is a “spin insulator” but is a charge conductor. Therefore, the system does conduct “number of up spin electrons” for either signs of $`U`$.

To summarize, in this Letter I have discussed a relation among the commensurability, the conduction property of the system, and the finite-size energy spectrum, combining the ideas from Refs. Kohn ; SWZ  and those from Refs. LSM ; Laughlin ; OYA ; comme . As a consequence, I have derived an upper bound for energies of some low-lying states for an insulator with an incommensurate filling in one and two dimensions. It is amusing to observe that several of very fundamental concepts LSM ; Laughlin ; Luttinger ; AB ; Kohn ; RestaSorella  in condensed matter physics, some of them put forward about 40 years ago, are rather deeply connected with each other, and to the simple notion of the commensurability.

As a final remark, I re-examine the validity of the “gap protection” hypothesis made in Ref. comme . Suppose the system has a finite gap (at zero AB flux) above the ground state(s). If the hypothesis is valid (i.e. the gap remains finite for any value of the AB flux), the present argument concludes that the system is an insulator ($`D=0`$). Therefore, if there is a conductor ($`D>0`$) with a finite gap, it contradicts with the hypothesis. According to Ref. SWZ , in the presence of the gap, the superfluid weight $`D_s`$ is identical to the Drude weight $`D`$, thus the conductor with a finite gap automatically means a superconductor.

In fact, in the presence of the gauge (Coulomb) interaction, a superconductor generally has a finite mass gap, due to the Anderson-Higgs mechanism AH . Thus, the gauge (Coulomb) interaction, being long-ranged, is enough to invalidate the hypothesis. If there is no long-range Coulomb interaction, a “superconductor” (or rather a superfluid) is generally associated to off-diagonal long-range or quasi-long-range order, and thus to gapless excitations. Hence the hypothesis may still be valid for a systems with short-range interactions (and hoppings).

Recently, Misguich et al. Misg:degeneracy  verified the “gap protection” hypothesis numerically for the RVB spin liquid state in $`d=2`$. They actually proposed an even stronger hypothesis (at least for the spin liquid state in $`d=2`$) that the ground-state energy $`E(\mathrm{\Phi })`$ is independent of the AB flux $`\mathrm{\Phi }`$ in the thermodynamic limit. This “flatness” hypothesis in Ref. Misg:degeneracy  may be more naturally understood with respect to the conduction property. Namely, if the groundstate is always below the gap, its energy $`E(\mathrm{\Phi })`$ must satisfy $`E(\mathrm{\Phi }_0)=E(0)`$. Thus, if $`E(\mathrm{\Phi })`$ varies with $`\mathrm{\Phi }`$, it must be convex at one point while it is concave at another. According to eq. (7) the system must then change its character from a superconductor to an insulator due to the insertion of the AB flux. Excluding such an unnatural behavior leads to the “flatness” hypothesis in $`d2`$. On the other hand, this argument does not support the “flatness” hypothesis for $`d3`$.

It is a pleasure to thank Leon Balents, Claudio Chamon, Gregoire Misguich and Masaaki Nakamura for very useful comments. The present work is supported in part by Grant-in-Aid for Scientific Research from MEXT of Japan.


text:
# Direct measurements of the spin and the cyclotron gaps in a 2D electron system in silicon

## Abstract

Using magnetocapacitance data in tilted magnetic fields, we directly determine the chemical potential jump in a strongly correlated two-dimensional electron system in silicon when the filling factor traverses the spin and the cyclotron gaps. The data yield an effective $`g`$ factor that is close to its value in bulk silicon and does not depend on filling factor. The cyclotron splitting corresponds to the effective mass that is strongly enhanced at low electron densities.

A two-dimensional (2D) electron system in silicon metal-oxide-semiconductor field-effect transistors (MOSFETs) is remarkable due to strong electron-electron interactions. The Coulomb energy overpowers both the Fermi energy and the cyclotron energy in accessible magnetic fields. The Landau-level-based considerations of many-body gaps ando ; yang , which are valid in the weakly interacting limit, cannot be directly applied to this strongly correlated electron system. In a perpendicular magnetic field, the gaps for charge-carrying excitations in the spectrum should originate from cyclotron, spin, and valley splittings and be related to a change of at least one of the following quantum numbers: Landau level, spin, and valley indices. However, the gap correspondence to a particular single-particle splitting is not obvious brener , and the origin of the excitations is unclear. In a recent theory iordan , the strongly interacting limit has been studied, and it has been predicted that in contrast to the single-particle picture, the many-body gap to create a charge-carrying (iso)spin texture excitation at integer filling factor is determined by the cyclotron energy. This is also in contrast to the square-root magnetic field dependence of the gap expected in the weakly interacting limit ando ; yang .

A standard experimental method for determining the gap value in the spectrum of the 2D electron system in a quantizing magnetic field is activation energy measurements at the minima of the longitudinal resistance englert ; klein ; dol88 ; usher . Its disadvantage is that it yields a mobility gap which may be different from the gap in the spectrum. In Si MOSFETs, the activation energy as a function of magnetic field was reported to be close to half of the single-particle cyclotron energy for filling factor $`\nu =4`$, while decreasing progressively for the higher $`\nu `$ cyclotron gaps englert ; klein ; dol88 . At low electron densities, an interplay was observed between the cyclotron and the spin gaps, manifested by the disappearance of the cyclotron ($`\nu =4`$, 8, and 12) minima of the longitudinal resistance krav00 . On the contrary, for the 2D electrons in GaAs/AlGaAs heterostructures, the activation energy at $`\nu =2`$ exceeded half the single-particle cyclotron energy by about 40% usher . Another, direct method for determining the gap in the spectrum is measurement of the chemical potential jump across the gap smith85 ; aristov ; valley . It was applied to the 2D electrons in GaAs smith85  and gave cyclotron gap values corresponding to the band electron mass aristov . Recently, the method has been used to study the valley gap at the lowest filling factors in the 2D electron system in silicon which has been found to be strongly enhanced and increase linearly with magnetic field valley ; rem .

The effective electron mass, $`m`$, and $`g`$ factor in Si MOSFETs have been determined lately from measurements of the parallel magnetic field of full spin polarization in this electron system and of the slope of the metallic temperature dependence of the conductivity in zero magnetic field gm . It is striking that the effective mass becomes strongly enhanced with decreasing electron density, $`n_s`$, while the $`g`$ factor remains nearly constant and close to its value in bulk silicon. This result is consistent with accurate measurements of $`m`$ at low $`n_s`$ by analyzing the temperature dependence of the Shubnikov-de Haas oscillations in weak magnetic fields in the low-temperature limit us ; rem1 . A priori it is unknown whether or not the so-determined values $`g`$ and $`m`$ correspond to the spin and the cyclotron splittings in strong perpendicular magnetic fields.

In this paper, we report the first measurements of the chemical potential jump across the spin and the cyclotron gaps in a 2D electron system in silicon in tilted magnetic fields using a magnetocapacitance technique. We find that (i) the $`g`$ factor is close to its value in bulk silicon and does not change with filling factor, in contrast to the strong dependence of the valley gap on $`\nu `$; and (ii) the cyclotron splitting is determined by the effective mass that is strongly enhanced at low electron densities. We also verify the systematics of the gaps in that the measured $`\nu =4`$, 8, and 12 cyclotron gap decreases with parallel magnetic field component by the same amount as the $`\nu =2`$, 6, and 10 spin gap increases.

Measurements were made in an Oxford dilution refrigerator with a base temperature of $`30`$ mK on high-mobility (100)-silicon MOSFETs (with a peak mobility close to 2 m2/Vs at 4.2 K) having the Corbino geometry with diameters 250 and 660 $`\mu `$m. The gate voltage was modulated with a small ac voltage 15 mV at frequencies in the range 2.5 – 25 Hz and the imaginary current component was measured with high precision using a current-voltage converter and a lock-in amplifier. Care was taken to reach the low frequency limit where the magnetocapacitance, $`C(B)`$, is not distorted by lateral transport effects. A dip in the magnetocapacitance at integer filling factor is directly related to a jump, $`\mathrm{\Delta }`$, of the chemical potential across a corresponding gap in the spectrum of the 2D electron system, and therefore we determine $`\mathrm{\Delta }`$ by integrating $`C(B)`$ over the dip in the low temperature limit where the magnetocapacitance saturates and becomes independent of temperature valley .

Typical magnetocapacitance traces taken at different electron densities, temperatures, and tilt angles are displayed in Fig. 1 near the filling factor $`\nu =hcn_s/eB_{}=4`$ and $`\nu =6`$. The magnetocapacitance shows narrow minima at integer $`\nu `$ which are separated by broad maxima, the oscillation pattern reflecting the modulation of the thermodynamic density of states, $`D`$, in quantizing magnetic fields: $`1/C=1/C_0+1/Ae^2D`$ (where $`C_0`$ is the geometric capacitance between the gate and the 2D electrons, and $`A`$ is the sample area) smith85 . As the magnetic field is increased, the maximum $`C`$ approaches the geometric capacitance indicated by the dashed lines in Fig. 1. Since the magnetocapacitance $`C(B)<C_0`$ around each maximum is almost independent of magnetic field, this results in asymmetric minima of $`C(B)`$, the asymmetry being more pronounced for $`\nu =4`$, 8, and 12. The chemical potential jump at integer $`\nu =\nu _0`$ is determined by the area of the dip in $`C(B)`$:

$$\mathrm{\Delta }=\frac{Ae^3\nu _0}{hcC_0}_{\text{dip}}\frac{C_{\text{ref}}C}{C}𝑑B_{},$$
(1)
where $`C_{\text{ref}}`$ is a step function that is defined by two reference levels corresponding to the capacitance values at the low and high field edges of the dip as shown by the dotted line in Fig. 1. The so-determined $`\mathrm{\Delta }`$ is smaller than the level splitting by the level width. The last is extracted from the data by substituting $`(C_0C_{\text{ref}})/C`$ for the integrand in Eq. (1) and integrating for the case of resolved levels between the magnetic fields $`B_1=hcn_s/e(\nu _0+1/2)`$ and $`B_2=hcn_s/e(\nu _01/2)`$.

Tilting the magnetic field allows us to verify the systematics of the gaps in the spectrum and probe the lowest-energy charge-carrying excitations. As the thickness of the 2D electron system in Si MOSFETs is small compared to the magnetic length in accessible fields, the parallel field couples largely to the electrons’ spins while the orbital effects are suppressed simonian . Therefore, the variation of a gap with $`B_{}`$ should reflect the change in the excitation energy as the Zeeman splitting, $`g\mu _BB`$, is increased: the excitation energy change is determined by the difference between the spin projections onto magnetic field for the ground and the lowest excited states. Within single-particle picture, e.g., one can expect that with increasing $`B_{}`$ at fixed $`B_{}`$, the spin gap will increase, the valley gap will stay constant, and the cyclotron gap, which is given by the difference between the cyclotron splitting and the sum of the spin and the valley splittings, will decrease. In contrast, for spin textures (so-called skyrmions), the dependence of the excitation energy on $`B_{}`$ should be much stronger compared to the single-particle Zeeman splitting yang .

In Fig. 2(a), we show the value of the chemical potential jump, $`\mathrm{\Delta }_s`$, across the $`\nu =2`$ and $`\nu =6`$ gaps as a function of magnetic field for different tilt angles. It is insensitive to both filling factor and tilt angle, as expected for spin gaps. The data are best described by a proportional increase of the gap with the magnetic field with a slope corresponding to an effective $`g`$ factor $`g1.75`$. The so-determined value obviously gives a lower boundary for the $`g`$ factor because both the valley splitting at odd $`\nu `$ and the level width are disregarded.

In Fig. 2(b), we show how the $`\nu =6`$ gap changes with $`B_{}`$ at different values of the perpendicular field component. It is noteworthy that the level width contribution, which is indicated by systematic error bars, depends weakly on parallel field, and the valley splitting has been verified to be independent of $`B_{}`$. This, therefore, allows more accurate determination of the $`g`$ factor as shown by the solid line in Fig. 2(b). Its slope yields $`g2.6`$, which is in agreement with the data obtained for the $`\nu =2`$ and $`\nu =10`$ gaps. The fact that this value is close to the $`g`$ factor $`g=2`$ in bulk silicon points to the single spin-flip origin of the excitations for the $`\nu =2`$, 6, and 10 gaps.

Unlike spin gaps, the chemical potential jump, $`\mathrm{\Delta }_c`$, across the $`\nu =4`$, 8, and 12 gaps decreases with parallel magnetic field component, as already seen from Fig. 1. In Fig. 3(a), we compare the behaviors of the $`\nu =6`$ and $`\nu =4`$ gaps with $`B_{}`$ at fixed perpendicular field component. For $`B_{}`$ between 2.7 and 6.6 T, the absolute values of the slopes of these dependences are equal, within experimental uncertainty, to each other so that the sum of the gaps is approximately constant even if the level width contribution is taken into account. These results lead to two important consequences: (i) the $`\nu =4`$, 8, and 12 gaps are cyclotron ones, the conventional systematics of the gaps remaining valid in the studied electron density range down to $`1.5\times 10^{11}`$ cm-2; and (ii) the $`g`$ factor does not vary with filling factor $`\nu `$. Although our value of $`g2.6`$ is in agreement with the previously measured ones englert ; gm ; us , we do not confirm the conclusion on oscillations of the $`g`$ factor with $`\nu `$ based on activation energy measurements and made in line with theoretical predictions ando ; afs  under the assumption of $`B_{}`$-independent level width englert .

In Fig. 3(b), we compare the data for the chemical potential jump across the $`\nu =4`$, 8, and 12 gaps in perpendicular and tilted magnetic fields including the term $`g\mu _B(BB_{})`$ that describes the increase of the spin gap with $`B_{}`$. The data coincidence confirms that the changing spin gap is the only cause for the dependence of the cyclotron gap on parallel field component. As is evident from the figure, $`\mathrm{\Delta }_c`$ is considerably smaller than the value ($`\mathrm{}\omega _c2\mu _BB_{}`$) expected within single-particle approach ignoring both valley splitting and level width.

To reduce experimental uncertainty related to the inaccurate determination of the level width, we plot in Fig. 4 the difference, $`(\mathrm{\Delta }_c\mathrm{\Delta }_s)/2\mu _BB`$, of the normalized values of the cyclotron and the spin gaps in a perpendicular magnetic field as a function of electron density. Assuming that the cyclotron splitting is determined by the effective mass $`m`$, this difference corresponds to ($`m_e/mg`$), where $`m_e`$ is the free electron mass. Using data for $`m`$ and $`g`$ obtained in both parallel gm  and weak us ; smith72  magnetic fields, we find that the value ($`m_e/mg`$) is indeed consistent with our data, see Fig. 4. The effective mass determined from our high-$`n_s`$ data using $`g=2.6`$ is equal to $`m0.23m_e`$, which is close to the band mass of $`0.19m_e`$. As long as our $`g`$ value is constant, the decrease of the normalized gap difference with decreasing $`n_s`$ reflects the behavior of the cyclotron splitting, which is in agreement with the conclusion of the strongly enhanced effective mass at low electron densities gm ; us .

We now discuss comparatively the results obtained for the valley and the spin gaps. According to Ref. valley , the enhanced valley gap at the lowest filling factors $`\nu =1`$ and $`\nu =3`$ in Si MOSFETs is comparable to the single-particle Zeeman splitting. As our data for the spin gap correspond to the single-particle Zeeman splitting, this may lead to a different systematics of the gaps in the spectrum compared to the single-particle picture. Such a possibility has been supported by a recent theory brener  which shows the importance of the Jahn-Teller effect for the ground state of a 2D electron system in bivalley (100)-Si MOSFETs in quantizing magnetic fields. At $`\nu =2`$, due to static and dynamic lattice deformations, the valley degeneracy is predicted to lift off giving rise to a complicated phase diagram including three phases: spin-singlet, canted antiferromagnet, and ferromagnet. In our experiment, over the studied range of magnetic fields down to 3 T, we observe at $`\nu =2`$ a spin-ferromagnetic ground state only. This gives an estimate of the strength of the suggested mechanism for the valley splitting enhancement.

The fact that we do not observe oscillations of the $`g`$ factor as a function of $`\nu `$ is not too surprising, because our value of $`g`$ is close to the $`g`$ factor in bulk silicon so that those oscillations may be small. At the same time, our data for the $`g`$ factor allow us to arrive at a conclusion that at $`\nu =2`$, the valley gap is small compared to the spin gap. Therefore, the valley splitting does oscillate with filling factor ando , the conclusion being valid, at least, for the strongly enhanced gaps at $`\nu =1`$ and $`\nu =3`$. We stress that this effect occurs in the strongly correlated electron system, which is beyond the conventional theory of exchange-enhanced gaps ando .

Let us finally discuss the results obtained for the cyclotron gap. The data of Fig. 4 indicate unequivocally that the origin of the small $`\mathrm{\Delta }_c`$ value in Fig. 3(b) is not related to valley splitting and level width. Instead, it is renormalization of the effective mass and $`g`$ factor due to electron-electron interactions: the observed decrease of the gap difference with decreasing $`n_s`$ in Fig. 4 as well as the systematics of the gaps are in agreement with both the decrease of the ratio of the cyclotron and the spin gaps with decreasing $`n_s`$ krav00  and the sharp increase of the effective mass at low electron densities gm ; us . Needless to say that the conventional theory ando  yields an opposite sign of the interaction effect on the cyclotron splitting.

We gratefully acknowledge discussions with I. L. Aleiner, S. V. Iordanskii, A. Kashuba, and S. V. Kravchenko. This work was supported by the RFBR, the Russian Ministry of Sciences, and the Programme “The State Support of Leading Scientific Schools”. V.T.D. acknowledges support of A. von Humboldt foundation via Forschungspreis.