text:
# Fluctuations in Stationary non Equilibrium States

## Abstract

In this paper we formulate a dynamical fluctuation theory for stationary non equilibrium states (SNS) which covers situations in a nonlinear hydrodynamic regime and is verified explicitly in stochastic models of interacting particles. In our theory a crucial role is played by the time reversed dynamics. Our results include the modification of the Onsager–Machlup theory in the SNS, a general Hamilton–Jacobi equation for the macroscopic entropy and a non equilibrium, non linear fluctuation dissipation relation valid for a wide class of systems.

The Boltzmann–Einstein theory of equilibrium thermodynamic fluctuations, as described for example in Landau–Lifshitz LL , states that the probability for a fluctuation from equilibrium in a macroscopic region of volume $`V`$ is proportional to $`\mathrm{exp}\{V\mathrm{\Delta }S/k\}`$ where $`\mathrm{\Delta }S`$ is the variation of entropy density calculated along a reversible transformation creating the fluctuation and $`k`$ is the Boltzmann constant. This theory is well established and has received a rigorous mathematical formulation in classical equilibrium statistical mechanics via the so called large deviation theory La . The rigorous study of large deviations has been extended to hydrodynamic evolutions of stochastic interacting particle systems KOV . In a dynamical setting one may asks new questions, for example what is the most probable trajectory followed by the system in the spontaneous emergence of a fluctuation or in its relaxation to equilibrium. The Onsager–Machlup approach ON2  answers precisely to this question: in the situation of a linear hydrodynamic equation, that is, close to equilibrium, the most probable emergence and relaxation trajectories are one the time reversal of the other. Developing the methods of KOV , this theory has been extended to nonlinear regimes JLV . Onsager–Machlup assume the reversibility of the microscopic dynamics; however microscopically non reversible models were constructed where the above results still hold, GJL1 ; GJL2 .

Emergence of large fluctuations, including Onsager–Machlup symmetry, has been observed in stochastically perturbed gradient type electronic devices LM . In their work, these authors study also non gradient (i.e. non reversible) systems and observe violation of Onsager–Machlup symmetry.

In the present paper we formulate a general theory of large deviations for irreversible processes, i.e. when detailed balance condition does not hold. This question was previously addressed in EY . Natural examples are boundary driven stationary non equilibrium states (SNS), e.g. a thermodynamic system in contact with two reservoirs, but our theory covers, as a special case, also the model systems considered in LM . In our approach a crucial role is played by the time reversed dynamics with respect to the stationary non equilibrium ensemble.

Our results are:

1. The Onsager–Machlup relationship has to be modified: the emergence of a fluctuation takes place along a trajectory which is determined by the time reversed process.

2. We show that the macroscopic entropy solves a Hamilton–Jacobi equation generalizing to a thermodynamic context known results for finite dimensional Langevin equations FW  as those studied in LM . The Hamilton–Jacobi equation can be solved perturbatively if we consider not too large fluctuations.

3. We test our theory in a stochastic model of interacting particles system, the boundary driven zero range process, in which we perform the computations explicitly. In particular it is possible to construct the microscopic time reversed process and to write the macroscopic entropy in a closed form.

4. For a large class of systems we obtain a non equilibrium non linear fluctuation dissipation relationship which links the macroscopic evolution of the system and of its time reversal to the thermodynamic force.

We are interested in many body system in the limit of infinitely many degrees of freedom. The basic assumptions of our theory are the following.

1) The microscopic evolution is given by a Markov process $`X_t`$ which represents the configuration of the system at time $`t`$. This hypothesis probably is not so restrictive because also the Hamiltonian case discussed in EPR  in the end is reduced to the analysis of a Markov process. The stationary non equilibrium state (SNS) is described by a stationary, i.e. invariant with respect to time shifts, probability distribution $`P_{st}`$ over the trajectories of $`X_t`$.

2) The macroscopic behavior of the system is described by diffusion type hydrodynamical equations of the form

$$_t\rho =\underset{i,j}{\overset{d}{}}_{u_i}\{D_{i,j}(\rho )_{u_j}\rho \}=F(\rho )$$
(1)
where $`\rho =\rho _t(u)`$ represents in general a vector of thermodynamic variables, e.g. the densities of different species of particles, and $`D_{i,j}`$ is a matrix acting on this vector. The interaction with the reservoirs appears as boundary condition to be imposed on solutions of (1). We assume that there exists a unique stationary solution $`\overline{\rho }`$ of (1), i.e. a profile $`\overline{\rho }(u)`$ which satisfies the appropriate boundary conditions such that $`F(\overline{\rho })=0`$.

These equations are derived from the underlying microscopic dynamics through an appropriate scaling limit. The hydrodynamic equations represent laws of large numbers with respect the probability measure $`P_{st}`$ conditioned on an initial state $`X_0`$. The initial conditions for (1) are determined by $`X_0`$. Of course many microscopic configurations give rise to the same value of $`\rho _0(u)`$. In general $`\rho _t(u)`$ is an appropriate limit of a $`\rho _N(X_t)`$ as the number $`N`$ of degrees of freedom diverges.

3) The stationary measure $`P_{st}`$ admits a principle of large deviations describing the fluctuations of the thermodynamic variables appearing in the hydrodynamic equations. This means the following. The probability that in a macroscopic volume $`V`$ containing $`N`$ particles the evolution of the variable $`\rho _N`$ deviates from the solutions of the hydrodynamic equations and is close to some trajectory $`\widehat{\rho }_t`$, is exponentially small and of the form

$`P_{st}\left(\rho _N(X_t)\widehat{\rho }_t(u),t[t_1,t_2]\right)`$ $``$ $`e^{N[S(\widehat{\rho }_{t_1})+J(\widehat{\rho })]}`$ (2)
$`=`$ $`e^{NI(\widehat{\rho })}`$
where $`J(\widehat{\rho })`$ is a functional which vanishes if $`\widehat{\rho }_t`$ is a solution of (1) and $`S(\widehat{\rho }_{t_1})`$ is the entropy cost to produce the initial value $`\widehat{\rho }_{t_1}`$. We adopt the convention for the entropy sign opposite to the usual one, so that it takes the minimum value in the equilibrium state. We also normalize it so that $`S(\overline{\rho })=0`$. Therefore $`J(\widehat{\rho })`$ represents the extra cost necessary in order that the system follow the trajectory $`\widehat{\rho }_t`$. Finally $`\rho _N(X_t)\widehat{\rho }_t(u)`$ means closeness in some metric. This formula is a generalization of the Boltzmann–Einstein formula in which we set the Boltzmann constant $`k=1`$.

4) Let us denote by $`\theta `$ the time inversion operator defined by $`\theta X_t=X_t=X_t^{}`$. The probability measure $`P_{st}^{}`$ describing the evolution of the time reversed process $`X_t^{}`$ is given by the composition of $`P_{st}`$ and $`\theta ^1`$ that is

$$P_{st}^{}(X_t^{}=\varphi _t,t[t_1,t_2])=P_{st}(X_t=\varphi _t,t[t_2,t_1])$$
(3)
We assume that the time reversed SNS also admits a hydrodynamic level of description.

Let $`L`$ be the generator of the microscopic dynamics. We remind that $`L`$ induces the evolution of functions of the process according to the equation $`_tE[f(X_t)]=E[(Lf)(X_t)]`$, where $`E`$ stands for the expectation with respect to $`P_{st}`$ conditioned on the initial state $`X_0`$ FE . The time reversed dynamics is generated by the adjoint $`L^{}`$ of $`L`$ with respect to the invariant measure $`P_{in}`$, that is $`E^{P_{in}}[fLg]=E^{P_{in}}[(L^{}f)g]`$.The measure $`P_{in}`$, which is the same for both processes, is a distribution over the configurations of the system and formally satisfies $`L^{}P_{in}=0`$. $`E^{P_{in}}`$ is the expectation with respect to $`P_{in}`$ and $`f,g`$ are functions of the configuration. We require that also the evolution generated by $`L^{}`$ admits a hydrodynamic description, that we call the adjoint hydrodynamics, which, however, is not necessarily of the same form as (1). In fact the adjoint hydrodynamics can be non local in space.

In order to avoid confusion we emphasize that what it is usually called an equilibrium state, as distinguished from a SNS, corresponds to the special case $`L^{}=L`$, i.e. the detailed balance principle holds. In such a case $`P_{st}`$ is invariant under time reversal and the two hydrodynamics coincide.

We now derive a fundamental consequence of our assumptions, that is the relationship between the action functionals $`I`$ and $`I^{}`$ associated to the dynamics $`L`$ and $`L^{}`$. From equation (3) and our assumptions it follows that also $`P_{st}^{}`$ admits a large deviation principle with functional $`I^{}`$ given by

$$I_{[t_1,t_2]}^{}(\widehat{\rho })=I_{[t_2,t_1]}(\theta \widehat{\rho })$$
(4)
with obvious notations. More explicitly this equation reads

$$S(\widehat{\rho }_{t_1})+J_{[t_1,t_2]}^{}(\widehat{\rho })=S(\widehat{\rho }_{t_2})+J_{[t_2,t_1]}(\theta \widehat{\rho })$$
(5)
where $`\widehat{\rho }_{t_1},\widehat{\rho }_{t_2}`$ are the initial and final points of the trajectory and $`S(\widehat{\rho }_{t_i})`$ the entropies associated with the creation of the fluctuations $`\widehat{\rho }_{t_i}`$ starting from the SNS. The functional $`J^{}`$ vanishes on the solutions of the adjoint hydrodynamics. To compute $`J^{}`$ it is necessary to know the macroscopic entropy $`S`$, which is determined by the microscopic invariant measure $`P_{in}`$. We shall see that $`S`$ can be also obtained as the solution of a Hamilton–Jacobi equation involving only macroscopic quantities. From (5) we can already obtain the generalization of Onsager-Machlup relationship for SNS.

The physical situation we are considering is the following. The system is macroscopically in the stationary state $`\overline{\rho }`$ at $`t=\mathrm{}`$ but at $`t=0`$ we find it in the state $`\widehat{\rho }_0`$. We want to determine the most probable trajectory followed in the spontaneous creation of this fluctuation. According to (2) this trajectory is the one that minimizes $`J`$ among all trajectories connecting $`\overline{\rho }`$ to $`\widehat{\rho }_0`$ in the time interval $`[\mathrm{},0]`$. From (5) we have

$$J_{[\mathrm{},0]}(\widehat{\rho })=J_{[0,\mathrm{}]}^{}(\theta \widehat{\rho })+S(\widehat{\rho }_0)$$
(6)
The right hand side is minimal if $`J_{[0,\mathrm{}]}^{}(\theta \widehat{\rho })=0`$ that is if $`\theta \widehat{\rho }`$ is a solution of the adjoint hydrodynamics. The existence of such a relaxational solution is due to the fact that the stationary solution $`\overline{\rho }`$ is attractive also for the adjoint hydrodynamics. We have therefore the following generalization of Onsager-Machlup to SNS

“In a SNS the spontaneous emergence of a macroscopic fluctuation takes place most likely following a trajectory which is the time reversal of the relaxation path according to the adjoint hydrodynamics”

Let us assume $`J`$ of the form

$$J_{[t_1,t_2]}(\widehat{\rho })=\frac{1}{2}_{t_1}^{t_2}𝑑tW,K(\widehat{\rho })W$$
(7)
where $`W=_t\widehat{\rho }F(\widehat{\rho })`$, $`F(\widehat{\rho })`$ has been defined in (1), $`,`$ denotes integration in the space variable and $`K(\widehat{\rho })`$ is an appropriate positive kernel which we assume to be known and reflects, at the macroscopic level, the stochasticity of the system. This form of $`J`$ is typical for diffusion processes described by finite dimensional Langevin equations FW  and is what we expect if the hydrodynamic equations are of the diffusion type (1). This is also what we find for the model discussed later.

From (5) or (6) we have that the entropy is related to $`J`$ by

$$S(\rho )=\underset{\widehat{\rho }}{inf}J_{[\mathrm{},0]}(\widehat{\rho })$$
(8)
where the infimum is taken over all trajectories connecting $`\overline{\rho }`$ to $`\rho `$. Therefore $`S`$ must satisfy the Hamilton–Jacobi equation associated to the action functional $`J`$. A simple calculation gives

$$\frac{1}{2}\frac{\delta S}{\delta \rho },K^1(\rho )\frac{\delta S}{\delta \rho }+\frac{\delta S}{\delta \rho },F=0$$
(9)
One can try to solve this functional derivative equation by successive approximations. Let $`\overline{\rho }(u)`$ be the stationary profile in the SNS. If we expand $`S`$ as a Volterra functional series in the argument $`\rho \overline{\rho }`$, equation (9) reduces to a system of equations for the kernels of the expansion which can be solved by iteration. Applications of (9) to variuos models will be discussed in NOI .

We consider now the so called zero–range process which models a nonlinear diffusion of a lattice gas, see e.g. KL . The model is described by a positive integer variable $`\eta _\tau (x)`$ representing the number of particles at site $`x`$ and time $`\tau `$ of a finite lattice which for simplicity we assume one–dimensional. The particles jump with rates $`g(\eta (x))`$ to one of the nearest-neighbor sites $`x+1,x1`$ with probability $`1/2`$. The function $`g(k)`$ is non decreasing and $`g(0)=0`$. We assume that our system interacts with two reservoirs of particles in positions $`N`$ and $`N`$ with rates $`p_+`$ and $`p_{}`$ respectively. The microscopic dynamics is then defined by the generator DF

$`(L_Nf)(\eta )`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{x=N}{\overset{N}{}}}g(\eta (x))[^{x,x+1}f+^{x,x1}f]`$
$`+`$ $`{\displaystyle \frac{1}{2}}p_+[f(\eta ^N)f(\eta )]+{\displaystyle \frac{1}{2}}p_{}[f(\eta ^N)f(\eta )]`$
where $`^{x,x\pm 1}f=f(\eta ^{x,x\pm 1})f(\eta )`$ and

$$\eta ^{x,y}(z)=\{\begin{array}{ccc}\eta (z)& \text{if}& zx,y\hfill \\ \eta (z)1& \text{if}& z=x\hfill \\ \eta (z)+1& \text{if}& z=y,\hfill \end{array}$$

that is, $`\eta ^{x,y}`$ is the configuration obtained from $`\eta `$ when a particle jumps from $`x`$ to $`y`$. Similarly $`\eta ^{N,N+1}(N)=\eta (N)1`$, $`\eta ^{N,N1}(N)=\eta (N)1`$ and $`\eta ^N(\eta ^N)`$ is the configuration $`\eta `$ after addition of a particle at the point $`N(N)`$.

It is remarkable that the invariant measure for this process can be constructed explicitly DF . It is the gran–canonical measure obtained by the product of the marginal distributions

$$P_N(\eta (x)=k)=\frac{\lambda _N^k(x)}{g(1)\mathrm{}g(k)}Z_N^1(x)$$
(11)
where

$$\lambda _N(x)=\frac{p_+p_{}}{2(N+1)}x+\frac{p_++p_{}}{2}$$
(12)
and

$$Z_N(x)=1+\underset{k=1}{\overset{\mathrm{}}{}}\frac{\lambda _N^k(x)}{g(1)\mathrm{}g(k)}$$
(13)
We emphasize that, if $`p_{}p_+`$, the generator $`L_N`$ is not self–adjoint with respect to the invariant measure so that the process is different from its time reversal and detailed balance does not hold.

Let us introduce now the macroscopic time $`t=\tau /N^2`$ and space $`u=x/N`$ and the empirical density

$$\rho _N(t,u)=\frac{1}{N}\underset{x=N}{\overset{N}{}}\eta _{N^2t}(x)\delta (ux/N)$$
(14)
Then one can prove that in the limit $`N\mathrm{}`$ the empirical density (14) tends in probability to a continuous function $`\rho _t(u)`$ which satisfies the following hydrodynamic equation

$$_t\rho =\frac{1}{2}_u\{D(\rho )_u\rho \}=\frac{1}{2}_u^2\varphi (\rho )=F(\rho )$$
(15)
where $`\varphi (\rho )`$ is the inverse function of

$$\rho (\varphi )=\frac{1}{Z(\varphi )}\underset{k=0}{\overset{\mathrm{}}{}}\frac{\varphi ^k}{g(1)\mathrm{}g(k)}k$$
(16)
where $`Z(\varphi )`$ is the normalization constant as in (13). The function $`\varphi `$ is well defined because the right hand side of (16) is strictly increasing in $`\varphi 0`$. Then $`\rho (\varphi )`$ is the equilibrium density corresponding to the activity $`\varphi `$. The boundary conditions for (15) are $`\varphi (\rho _t(\pm 1))=p_\pm `$.

From the knowledge of the invariant measure one can calculate the adjoint generator $`L_N^{}`$ given by

$`(L_N^{}f)(\eta )={\displaystyle \frac{1}{2}}{\displaystyle \underset{x=N}{\overset{N}{}}}g(\eta (x)){\displaystyle \frac{\lambda _N(x+1)}{\lambda _N(x)}}^{x,x+1}f`$ (17)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{x=N}{\overset{N}{}}}g(\eta (x)){\displaystyle \frac{\lambda _N(x1)}{\lambda _N(x)}}^{x,x1}f`$
$`+{\displaystyle \frac{1}{2}}\lambda _N(N)[f(\eta ^N)f(\eta )]+{\displaystyle \frac{1}{2}}\lambda _N(N)[f(\eta ^N)f(\eta )].`$
Notice that the form of (17) is the same as (Fluctuations in Stationary non Equilibrium States) with the rates modified in such a way to invert the particle flux. From (17) one can derive the adjoint hydrodynamics

$$_t\rho =\frac{1}{2}_u^2\varphi (\rho )\alpha _u\{\frac{\varphi (\rho )}{\lambda (u)}\}=F^{}(\rho )$$
(18)
with $`\lambda (u)=\frac{p_+p_{}}{2}u+\frac{p_++p_{}}{2}`$ and $`\alpha =\frac{p_+p_{}}{2}`$. The boundary condition for (18) are the same as for (15). For the detailed computations we refer the reader to NOI .

The second term on the right hand side of (18) is produced by the new rates in (17). As expected, it is proportional to the difference of the chemical potentials.

In this model one can compute the action functionals $`J(\widehat{\rho })`$ and $`J^{}(\widehat{\rho })`$. The result is

$`J_{[t_1,t_2]}(\widehat{\rho })`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{t_1}^{t_2}}𝑑t_u^1W,{\displaystyle \frac{1}{\varphi (\widehat{\rho })}}_u^1W`$ (19)
$`J_{[t_1,t_2]}^{}(\widehat{\rho })`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{t_1}^{t_2}}𝑑t_u^1W^{},{\displaystyle \frac{1}{\varphi (\widehat{\rho })}}_u^1W^{}`$ (20)
where, $`W^{}=_t\widehat{\rho }F^{}(\widehat{\rho })`$. Moreover $`J`$ and $`J^{}`$ are defined to be infinite if $`\widehat{\rho }_t`$ does not satisfy the boundary conditions stated above.

The entropy $`S(\rho )`$ can be easily computed from the expression (11) for the invariant measure

$$S(\rho )=𝑑u\left[\rho (u)\mathrm{log}\frac{\varphi (\rho (u))}{\lambda (u)}\mathrm{log}\frac{Z(\varphi (\rho (u)))}{Z(\lambda (u))}\right]$$
(21)
From equations (19), (20) and (21), one can verify explicitly equation (5) and the generalized Onsager-Machlup relationship.

We deduce now a twofold generalization of the celebrated fluctuation–dissipation relationship: it is valid in non equilibrium states and in non linear regimes. Such a relationship will hold provided the rate function $`J_{[t_1,t_2]}^{}`$ of the time reversed process is of the form (7) with the same metric $`K(\rho )`$ but a different vector field $`F^{}`$ describing the hydrodynamic of the adjoint process, namely

$$J_{[t_1,t_2]}^{}(\widehat{\rho })=\frac{1}{2}_{t_1}^{t_2}𝑑tW^{},K(\widehat{\rho })W^{}$$
(22)
By taking the variation of the equation (5), a simple computation gives

$$F(\rho )+F^{}(\rho )=K(\rho )^1\frac{\delta S}{\delta \rho }$$
(23)
This relation holds for the non–equilibrium zero–range model discussed before. We also note that it holds for the equilibrium reversible models for which the large deviation principle has been rigorously proven such as the simple exclusion process KOV  and the Landau–Ginzburg model DV  and its non–gradient version Q . It is also easy to check that the linearization of (23) around the stationary profile $`\overline{\rho }`$ yields a fluctuation–dissipation relationship which in equilibrium reduces to the usual one.

We are grateful to E. Presutti for insisting that the results obtained for special models should actually exemplify general principles. We also thank A. De Masi, P. Ferrari, G. Giacomin and M. E. Vares for very illuminating discussions. D. G. was supported by FAPESP 98/11899-2 at IME-USP, Brazil and is very grateful to A. Galves and P. Ferrari for the warm hospitality. He wishes to thank P. Markowich and the University of Vienna for support during the final stage of the work through the Wittgenstein award. L. B. and G. J-L. acknowledge the support of Cofinanziamento MURST.


format: text

text:
# Angle–Resolved Loss of Landau Quasiparticles in 2D Hubbard Model

## Abstract

The problem of weakly correlated electrons on a square lattice is studied theoretically. A simple renormalization group scheme for the angle–resolved weight $`Z(\theta )`$ of the quasiparticles at the Fermi surface is presented and applied to the Hubbard model. Upon reduction of the cutoff the Fermi surface is progressively destroyed from the van Hove points toward the zone diagonals. Due to the renormalized $`Z(\theta )`$, divergences of both antiferromagnetic and superconducting correlation functions are suppressed at the critical scale, where the interactions diverge.

LPTHE/01-17 April 2001

Understanding of the one–particle spectrum of strongly correlated systems near the metal–insulator transition is an extremely difficult task, particularly if one wants to construct a microscopic theory. One standard example is the pseudogap regime of the HTC superconductors. ARPES measurements ARPES\_psgap  showed that the Fermi surface is destroyed by correlations. This happens first near the van Hove points where the one-particle spectrum develops the characteristic 2-peak structure. Remaining parts of the Fermi surface, often called Fermi patches, get progressively narrower around Brillouin zone diagonals as the temperature decreases. Regions around van Hove points contain non–Fermi liquid with a pseudogap and other signatures of strong correlations such as flat bands.flat\_band  In other words, the pseudogap has the form similar to the absolute value of the $`d_{x^2y^2}`$–superconducting (SC) order–parameter. This vision of the pseudogap is in agreement with STM results tunel\_psgap  as well, regardless of details on how these results are interpreted. The above experiments still do not reveal much on the origins of the pseudogap, namely whether it is simply the signature of a liquid of pre–formed pairs or something much richer in fluctuations. In fact it is known that the antiferromagnetic (AF) fluctuations are also strong in the pseudogap regime.neutrons

If we assume that AF and SC fluctuations are somehow together the major reason for the strong renormalization of the one–particle selfenergy, and for the consequent partial destruction of the Fermi surface, then a many–body analysis of the one–particle propagator can be done in a controlled way. In fact the weak coupling theory easily reproduces AF and SC fluctuations from particle–hole (p-h) and particle–particle (p-p) loop–logarithms. Even if the coupling in realistic HTC systems is of the order of Fermi energy (i.e. intermediate–to–strong), already the weak–coupling theory contains the observed two–particle correlations. In this paper I answer the question of how the angle–resolved quasiparticle weight $`Z(\theta )`$ is renormalized by strong and coupled AF and SC fluctuations, and of the main consequences of the renormalized $`Z(\theta )`$ to the characteristic angle–resolved two–particle correlation functions. The renormalization of the quasiparticle weight in 2D was recently studied by Kishine and Yonemitsu.Kishine99  To calculate the renormalization of $`Z`$ resolved in the position on a flat Fermi surface they used a two–loop selfenergy expansion with the two–loop–renormalized vertex. The results show clearly that the flatness of the Fermi surface induces the suppression of the quasiparticle residue and that this effect is anisotropic. In the present work I consider the whole square Fermi surface of the Hubbard model. For this purpose I employ the N–patch renormalization group theory.

After recent theoretical studies from several groups it emerged that the N–patch model describes in a systematic and controlled way weakly correlated electrons near half–filling, and explains the major aspects of HTC-s.these ; doucot\_97 ; ZS\_prb\_00 ; Halboth\_00 ; Honerkamp\_01 ; Tsai\_Marston  Until now the RG analysis of the N–patch model has been done only on the level of the two–particle scattering amplitudes or, in field–theoretical jargon, of the four–point vertex $`U(K_1,K_2,K_3)`$. The analysis of the renormalization group flow of $`U`$ as a function of three patch indices $`i_1,i_2,i_3`$ gave several important results. Typically the amplitudes $`U(i_1,i_2,i_3)`$ diverge at some interaction– and doping–dependent critical energy scale $`\mathrm{\Lambda }_c`$. For the case of the Hubbard model, we distinguish two main renormalization regimes ZS\_prb\_00 , the parquet regime and the BCS regime.

In the parquet regime ($`|\mu |<\mathrm{\Lambda }`$) both particle-particle and particle-hole propagators have strong contributions to the beta–function due to the van Hove singularities and the Fermi surface nesting. In this regime, provided $`\mathrm{\Lambda }\mathrm{\Lambda }_c`$ and neglecting the selfenrgy corrections, both SC and AF tendencies are strong and build divergent correlation functions $`\chi ^{SC}`$ and $`\chi ^{AF}`$. The dominant component of the antiferromagnetic susceptibility is of the $`s`$–type and the dominant component of the superconducting one is of the $`d_{x^2y^2}`$–type. Both static compressibility $`\chi _c`$ and homogeneous magnetic susceptibility $`\chi _s`$ go to zero as the cutoff approaches to its critical value $`\mathrm{\Lambda }_c`$.these ; Halboth\_00 ; Honerkamp\_01  Consequently $`\mathrm{\Lambda }_c`$ is the energy scale of the crossover between the strange metal and the strongly correlated regime with gap or pseudogap. A question arising from already existing results on the Hubbard model these ; ZS\_prb\_00  and on its extensions Halboth\_00 ; Honerkamp\_01 ; Tsai\_Marston  is the following: if one wants to interpret the critical scale $`\mathrm{\Lambda }_c`$ (or temperature) as the energy (temperature) $`T^{}`$ for the onset of the pseudogap, why then are all signatures of the pseudogap not seen? This means in particular that $`\chi ^{AF}`$ and $`\chi ^{SC}`$ should be finite and not diverging. Emergent is the necessity to calculate the correlation functions with the corrections due to the one–particle selfenergy.

At stronger doping ($`|\mu |>\mathrm{\Lambda }`$) nesting properties get weaker so that eventually only remaining renormalization channel is superconducting (p-p). This is the BCS regime. There the superconductivity is simply BCS–like with the coupling constants and the angular profile of the order parameter determined at higher scales by a parquet–like flow, where $`\mathrm{\Lambda }`$ was larger than the chemical potential.

Equivalent to the RG approach is the fast parquet theory DzYak ; Zheleznyak , where $`\theta `$ variable (continuum version of the patch index) is called the fast variable in addition to the cutoff logarithm called the slow variable. In two dimensions the parquet integro–differential equations always have mobile pole solutions, i.e. the AF and SC fluctuations decouple one from the other. This seems to be in disagreement with the RG results, where we detected only immobile poles, the type of solution in which all scattering amplitudes develop the pole at the same scale $`\mathrm{\Lambda }_c`$. The question of the consilience between the two theories is still controversial. However, the results of De Abreu and Douçot doucot\_00  indicate that the mobile pole solution is dominant only in the very vicinity of $`\mathrm{\Lambda }_c`$ and that the fixed pole solution is an intermediate solution, valid over several decades of energy scale. The width of the “very vicinity” characterized by the mobile poles depends on the coupling constant so that for reasonable and not too weak $`U_0`$ the final regime is so close to $`\mathrm{\Lambda }_c`$, that the couplings are already too strong and out of reach of a weak coupling theory. Consequently the real physical interpretation can be given only to the immobile pole regime. We will concentrate on this “intermediate” regime in which all fluctuations are coupled and at least behave as if having an immobile pole.

We suppose that the electronic Green function has the form

$$G_l(K)=\frac{Z_l(\theta )}{i\omega \xi (𝐤)}.$$
(1)
$`Z_l(\theta )`$ is the angle–resolved scale dependent quasiparticle weight and $`\xi (𝐤)`$ is the tight–binding dispersion. The formalism keeps the notation introduced in the reference ZS\_prb\_00 . The form (1) contains two main approximations. The first one is to keep trace only of the renormalization of the coherent part of the propagator. The second approximation is to assume that the spectrum $`\xi (𝐤)`$ remains non–renormalized. This assumption implies that we ignore the flow of the Fermi surface (FS) and of the Fermi velocity. Because of the particle–hole symmetry the flow of the FS is zero at half–filling, where we can expect that our form of the Green function is closer to reality than in the imperfectly nested (non–half filled) case.

The flow equation for $`Z(\theta )`$ is derived from the general and exact one–loop RG equation for the complete selfenergy $`\mathrm{\Sigma }(K)`$ given in ZS\_prb\_00 . Let us suppose that we are at some scale $`l`$, and that we know the propagator (1). We integrate $`dl`$ further and look what is the effective two–point vertex $`\mathrm{\Gamma }_2`$ in the effective action $`S(l+dl)`$; it is

$$\mathrm{\Gamma }_2(l+dl)=Z_l^1(\theta )(i\omega \xi (𝐤))+d\mathrm{\Sigma }_l(K).$$
(2)
To find $`Z_{l+dl}`$ we expand $`d\mathrm{\Sigma }_l(K)`$ in first order in $`i\omega `$ to obtain

$$Z_{l+dl}=Z_l(1Z_l_{i\omega }d\mathrm{\Sigma }).$$
(3)
The differential equation for $`Z`$ follows immediately

$$_lZ_l(\theta )=Z_l^2(\theta )_{i\omega }\left[_l\mathrm{\Sigma }(K)\right]|_{\xi =i\omega =0}.$$
(4)
Only the terms of $`_l\mathrm{\Sigma }`$ which are linear in energy contribute. These are just the terms which are marginal upon zeroth order scaling in Shankar’s sense.Shankar  We will look for these terms.

The equation for $`_l\mathrm{\Sigma }`$ can be written as

$$_l\mathrm{\Sigma }(\theta ,ϵ,\omega )=\frac{\mathrm{\Lambda }}{(2\pi )^2}\frac{d\omega ^{}}{2\pi }\underset{\nu }{}𝒥_\nu (\theta ^{},\mathrm{\Lambda })𝑑\theta ^{}G_l(\theta ^{},\omega ^{},\nu \mathrm{\Lambda })𝒟_l(K,K_\nu ^{}),$$
(5)
where $`𝒟_l=2F_l\stackrel{~}{F}_l`$; $`F_l`$ and $`\stackrel{~}{F}_l`$ are energy–momenta dependent forward and backward scattering processes at the scale $`l`$, related to the effective interaction in a way that $`F_l(K_1,K_2)=U(K_1,K_2,K_1)`$ and $`\stackrel{~}{F}_l(K_1,K_2)=U(K_1,K_2,K_2)`$. $`𝒥_\nu (\theta ,\mathrm{\Lambda })`$ is the angle–resolved density of states at energy $`\xi =\nu \mathrm{\Lambda }`$. There is another, approximate but physically justified way to decompose $`𝒟`$. In fact, we will suppose that $`D_l`$ can be written as a sum of p-p and p-h terms:

$$𝒟_l(K,K^{})=𝒟_l^{pp}(K+K^{})+𝒟_l^{ph}(KK^{}).$$
(6)
The p-p part of the propagator $`𝒟`$ depends only on the total energy–momentum $`Q_{pp}=(\omega _{pp},𝐪_{pp})K+K^{}`$ while the p-h part depends only on the energy–momentum transfer $`Q_{ph}=(\omega _{ph},𝐪_{ph})KK^{}`$. As usual we skip the marginal part of the dependence on $`𝐪_{pp}`$ and $`𝐪_{ph}`$. For that purpose we note that both momenta can be written in the form

$$𝐪=𝐪^{(0)}(\theta ,\theta ^{})+𝐪^{(1)}(\theta ,\theta ^{},\xi ,\xi ^{}),$$
(7)
where $`𝐪`$ stands either for $`𝐪_{pp}`$ or for $`𝐪_{ph}`$, $`𝐪^{(0)}(\theta ,\theta ^{})`$ is the value of q when both momenta $`𝐤`$ and $`𝐤^{}`$ are at the Fermi surface, while $`𝐪^{(1)}`$ is the correction due to non–zero energies $`\xi `$ and $`\xi ^{}`$. Using standard scaling arguments we can skip $`𝐪^{(1)}`$ in the limit $`\mathrm{\Lambda }/ϵ_F0`$ because the cutoff is imposed to momenta. Similar argument cannot be used for $`\omega _{pp}`$ and $`\omega _{ph}`$ because the integral in (5) runs over all frequencies independently of the actual cutoff $`\mathrm{\Lambda }(l)`$. We are therefore left with the locally dispersionless phononic propagators

$$𝒟_l^{pp}(\theta ,\theta ^{},Q_{pp})𝒟_l^{pp}(\theta ,\theta ^{},i\omega +i\omega ^{})=2F_l^{pp}(\theta ,\theta ^{},i\omega +i\omega ^{})\stackrel{~}{F}_l^{pp}(\theta ,\theta ^{},i\omega +i\omega ^{})$$
(8)
and

$$𝒟_l^{ph}(\theta ,\theta ^{},Q_{ph})𝒟_l^{ph}(\theta ,\theta ^{},i\omega i\omega ^{})=2F_l^{ph}(\theta ,\theta ^{},i\omega i\omega ^{})\stackrel{~}{F}_l^{ph}(\theta ,\theta ^{},i\omega i\omega ^{}).$$
(9)
In these expressions we made the same pp-ph decomposition of the forward and backward amplitudes as we did with $`𝒟`$ in eq.6. The following step is to re-constitute the $`i\omega `$ dependence from the cutoff dependence. This can be done with logarithmic precision simply replacing $`\mathrm{\Lambda }`$ with $`i\omega `$. The derivatives over frequency in eq.(4) (acting only on $`F`$-parts) are then readily calculated

$$_{i\omega }F_l^{pp,ph}(\theta ,\theta ^{},i\omega \pm i\omega ^{})|_{i\omega =0}=\pm \frac{1}{i\omega ^{}}_lF_l^{pp,ph}(\theta ,\theta ^{}),$$
(10)
and equivalently for backward amplitudes $`\stackrel{~}{F}_l^{pp,ph}`$. The frequency–independent quantities $`_lF_l^{pp}(\theta ,\theta ^{})`$ and $`_lF_l^{ph}(\theta ,\theta ^{})`$ are the p-p and p-h parts of the $`\beta `$–function of the N–patch model, with appropriate configurations of the external momenta:

$$\begin{array}{cccc}_lF_l^{pp}(\theta ,\theta ^{})=\beta _{pp}\{U,U\}(\theta ,\theta ^{},\theta )\hfill & & & \\ _l\stackrel{~}{F}_l^{pp}(\theta ,\theta ^{})=\beta _{pp}\{U,U\}(\theta ,\theta ^{},\theta ^{})\hfill & & & \\ _lF_l^{ph}(\theta ,\theta ^{})=[X\beta _{ph}\{XU,XU\}](\theta ,\theta ^{},\theta )\hfill & & & \\ _l\stackrel{~}{F}_l^{ph}(\theta ,\theta ^{})=[2\beta _{ph}\{U,U\}\beta _{ph}\{U,XU\}\beta _{ph}\{XU,U\}](\theta ,\theta ^{},\theta ^{}),\hfill & & & \end{array}$$
(11)
where all $`\beta `$–functions are given in ref.ZS\_prb\_00 , but for the moment with dressed Green functions (1). Notice that the forward scattering has finite p-h contributions only from the ZS’ channel (1=3) while only the ZS channel (1=4) contributes to the backward scattering. This means that we forget about the contributions at zero momentum transfer. They are somewhat tricky, but don’t have any logarithmic part so that we can forget them.

We can also get rid of the $`Z`$ factors in beta–functions of the eq.(11) by rescaling the fermions at every step of the RG in a way that

$$\overline{\mathrm{\Psi }}(\theta )Z_l^1(\theta )\mathrm{\Psi }(\theta )\overline{\mathrm{\Psi }}(\theta )\mathrm{\Psi }(\theta )$$
(12)
and re–defining the effective interaction

$$U_l(1,2,3)[Z_l(1)Z_l(2)Z_l(3)Z_l(4)]^{1/2}U_l(1,2,3).$$
(13)
After transformations (12) and (12) the calculations of the $`\beta `$–functions to the one–loop order are identical to the case with $`Z=1`$. Performing $`\omega ^{}`$ integral to logarithmic precision, the flow equation for $`Z(\theta )`$ becomes

$$_l\mathrm{log}Z_l(\theta )=\frac{1}{(2\pi )^2}𝑑\theta ^{}𝒥_{}(\theta ^{},\mathrm{\Lambda })\eta _l(\theta ,\theta ^{})\eta _l(\theta ),$$
(14)
where

$$\eta _l(\theta ,\theta ^{})_l\{2[F_l^{pp}F_l^{ph}]\stackrel{~}{F}_l^{pp}+\stackrel{~}{F}_l^{ph}\}(\theta ,\theta ^{})$$

with $`_lF`$–terms given by eqs.(11) and calculated with the bare Green functions, as in ref.ZS\_prb\_00 . Generalization of eq.14 to finite temperatures can be done simply by replacing $`𝒥_{}`$ with $`_\nu (\nu )𝒥_\nu n_F(\nu \mathrm{\Lambda })`$.

Taking 1D “limit” of the eq.(14) is simple and instructive: instead of N patches we now have 2 patches: $`\theta =R`$ (right) and $`\theta =L`$ (left). Independent scattering amplitudes at non-rational filling are $`F(RL)=g_2=U(RLR)`$ and $`\stackrel{~}{F}(RL)=g_1=u(RLL)`$. We skip $`g_4=u(RRR)`$ from considerations because it has no logarithmic renormalization. The $`\theta `$–integrals reduce to summation over two points so that it is easy to reproduce the well–known result $`\eta =\frac{1}{4\pi ^2v_F^2}(g_1^2g_1g_2+g_2^2)`$ Q1D\_Bourbonnais . It is the Luttinger liquid exponent.

In two dimensions the angle resolved $`\eta _l(\theta )`$ can also be associated with some non-Landau (non-Fermi liquid) behavior. We see that $`\eta `$ becomes finite if the forward and backward amplitudes have some logarithmic flow over a wide range of $`\theta `$–space $`(\mathrm{\Delta }\theta 1)`$. Another possibility for having finite $`\eta `$ is when the Fermi surface is close to the van Hove singularities. “Close” means that the distance between the Fermi level and the van Hove singularity is comparable or inferior to the scale $`\mathrm{\Lambda }_c`$ at which interactions start to flow strongly.

We will now calculate the renormalization of $`Z_l(\theta )`$ in the 2D Hubbard model at half–filling, from the knowledge of the scale dependence of the patch-dependent interaction $`U_l(1,2,3)`$. In the RG equations of the previous section the discretization of $`\theta `$ is done in a way described in ref.ZS\_prb\_00 . In the present case the Fermi surface is square so that there are two mechanisms for the suppression of quasiparticle residues. Namely, both above mentioned conditions are fulfilled: (i) forward and backward amplitudes have logarithmic flows for any configuration $`(\theta ,\theta ^{})`$ if the two angles are at opposite sides of the Fermi surface, so that the available phase space is indeed large; (ii) van Hove singularities are at the Fermi surface and are nested. In fact, one can alternatively imagine a Fermi surface with non-nested van Hove singularities and nested parts elsewhere. Such a model would be even closer to the realistic situation in some HTC compounds. For the sake of rigor, we will however remain limited to the Hubbard model.

The result is shown in fig.1. The figure shows $`Z(\theta _i)`$ with $`0\theta _i\pi /2`$ on 9 equidistant points. Different lines correspond to different values of the scaling parameter as it approaches its critical point, i.e. when couplings diverge. Settings are the same as in ref.ZS\_prb\_00 : the Fermi surface was discretized into 32 patches and the initial interaction is $`U_0/(4t)=0.333`$. One sees that the Fermi surface is first destroyed at the van Hove points and than the regions of the FS destroyed by correlations grow larger and larger. This kind of flow is compatible with the interpretation that $`\mathrm{\Lambda }_c`$ is not the critical temperature for some symmetry braking, but merely the scale at which coherent quasiparticle cease to exist at the Fermi surface giving place to a gapped or pseudogapped liquid. The magnitude of the (pseudo-)gap is largest in van Hove points and smallest on the diagonals of the Brillouin zone. The whole gap function $`\mathrm{\Delta }(\theta )`$ then scales as $`\mathrm{\Lambda }_c\times f(\theta )`$, where $`f(\theta )`$ is a function with the symmetry of the absolute value of the $`d_{x^2y^2}`$–harmonic. This is the angle–resolved (pseudo-)gap, responsible for the correlation–induced angle–resolved localization. In other words the electrons near the van Hove points get much less mobile than those near diagonals. The similar scenario has been proposed by the Zürich group Honerkamp\_01  even without concrete calculations of the quasiparticle weight.

The question of antiferromagnetism and superconductivity remains to be clarified. Let’s discuss this problem taking into account the scale dependent $`Z(\theta )`$ in the flow equations for the susceptibilities $`\chi ^{AF}(\theta ,\theta ^{})`$ and $`\chi ^{SC}(\theta ,\theta ^{})`$. Following the procedure given in ref.ZS\_prb\_00  and dressing the electronic propagators with $`Z`$–factors we get

$$\dot{\chi }_l^\delta (\theta _1,\theta _2)=\frac{1}{Z_l(\theta _1)Z_l(\theta _2)}𝑑\theta \stackrel{~}{z}_l^\delta (\theta _1,\theta )D_l^\delta (\theta )\stackrel{~}{z}_l^\delta (\theta ,\theta _2).$$
(15)
This equation has the same structure as the one in ref.ZS\_prb\_00 , with two modifications. First, we skip the retardation effects, replacing $`l_\delta `$ simply by $`l`$, because we are at half–filling. Second, the quantity $`\stackrel{~}{z}_l^\delta (\theta _1,\theta )`$ that has the role of a triangular vertex is somewhat modified. Its flow writes:

$$\left[_l\eta (\theta _1)\eta (\theta _2)\right]\stackrel{~}{z}_l^\delta (\theta _1,\theta _2)=𝑑\theta \stackrel{~}{z}_l^\delta (\theta _1,\theta )D_l^\delta (\theta )V_l^\delta (\theta ,\theta _2).$$
(16)
The meaning of $`\stackrel{~}{z}_l^\delta (\theta _1,\theta _2)`$ is that

$$\stackrel{~}{z}_l^\delta (\theta _1,\theta _2)Z_l(\theta _1)z_l^\delta (\theta _1,\theta _2)Z_l(\theta _2)$$

so that the initial conditions for $`\stackrel{~}{z}_l^\delta `$ and for $`z_l^\delta `$ are the same. After discretization we integrate numerically equations (15) and (16). Fig.2 shows the flow of the dominant eigenvalues of susceptibilities $`\delta =AF`$ and $`\delta =SC`$ near the divergence of scattering amplitudes. The thin line represents both (degenerated) susceptibilities for $`U=0`$. Including only the one–loop vertex renormalization we get the strong enhancement and, as far as my numerics can say, even divergences of both AF and SC susceptibilities. If we now include also the one–particle–weight renormalization, both susceptibilities are radically reduced and lose their divergent behavior. On the other hand, the flow of the compressibility $`\chi _c`$ and magnetic susceptibility $`\chi _\sigma `$ is not affected by $`Z`$–renormalization because of the Ward identities.Metzner\_Ward  All above results support the statement that what happens at energy scale $`\mathrm{\Lambda }_c`$ is a flow towards a state with spin– and charge–gap or pseudogap, insulating and without AF or SC ordering.

We finish this discussion with a few words about the effects of doping. Two regimes, parquet and BCS, exist also in the scaling properties of $`Z_l(\theta )`$. In the whole parquet regime we expect the behavior governed by the proximity to the half–filling situation, so that the present results can be applied at all energies larger than the chemical potential $`|\mu |`$. Upon doping the nesting becomes more and more imperfect, the p-h logarithm loses its divergence, the critical scale is more and more suppressed and the p-h and p-p channels get progressively less coupled. Eventually at strong enough doping and low enough $`\mathrm{\Lambda }`$ one is in the BCS regime where the effective physics is described by the 2D BCS theory, with renormalized and $`\theta `$–dependent interaction and quasiparticle weight. The anomalous dimension $`\eta (\theta )`$ goes to zero because (i) the the range $`\mathrm{\Delta }\theta `$, over which the forward and the backward amplitude have strong flow due to BCS diagram, scales with the cutoff, and (ii) the van Hove singularities are outside of the cutoff. The critical temperature for the onset of the superconductivity is not affected by the renormalized $`Z`$, but the magnitude and angular dependence of the superconducting order parameter are dressed by $`Z(\theta )`$. The symmetry of the gap remains $`d_{x^2y^2}`$.

To summarize, I proposed a simple renormalization group theory for the angle–dependent destruction of the Fermi surface in the Hubbard model. The results offer a theoretical comprehension of the angle–dependent Fermi surface truncation in the cuprate superconductors in terms of the scattering processes of the electrons on the low–energy collective excitations of both particle–particle and particle–hole types. The theory, based on the N–patch model, is in its essence a controlled weak–coupling procedure that keeps trace of the dependence of the effective interaction and one–particle spectral weight on the position of the particles at the Fermi surface. As one approaches the critical scale $`\mathrm{\Lambda }_c`$ the quasiparticle weight goes to zero first near the van Hove points, and the effect progresses toward Brillouin zone diagonals as one lowers the temperature. Dressing the flow equations for AF and SC response functions with the one–particle weight factors results in dramatical reduction of correlations of both types. The strongly correlated state just below $`\mathrm{\Lambda }_c`$ is gapped or pseudogapped and without any long–range–order. Critical scale $`\mathrm{\Lambda }_c`$ is interpreted as the pseudogap temperature $`T^{}`$ found in cuprate superconductors.

###### Acknowledgements.

I am grateful to Serguei Brazovskii, Benoit Douçot, Benedikt Binz, Claude Bourbonnais, Nicolas Dupuis, and Sebastien Dusuel for important discussions and comments. Laboratoire de Physique Théorique et Hautes Energies, Universités Paris VI Pierre et Marie Currie – Paris VII Denis Diderot, is supported by CNRS as Unité Mixte de Recherche, UMR7589.


text:
# LATE AGB MAGNETIC CYCLES: MHD SOLUTIONS FOR THE HST PN RINGS

## 1 Introduction

The multiple concentric rings, or arcs, that were recently discovered by the Hubble Space Telescope around a handful of planetary nebulae (PNe) is one of the most puzzling and unexpected results delivered by the HST (see Kwok, Su, & Hrivnak 1998; Hrivnak, Kwok & Su 2001; Terzian & Hajian 2000 and references therein). The best documented case of these systems of faint concentric rings (hereafter called HST rings) is displayed by NGC 6543 (Balick, Wilson & Hajian 2000). In reality, they are regularly spaced concentric shells, indicating quasi-periodic events with time intervals, assuming typical expansion velocities of AGB winds, in the range of 500 to 1500 years. These time scales are too short for thermal pulses and too long for acoustic envelope pulsations and, thus, the origin of the rings cannot be ascribed to any of this type of events. Soker (2000) made a critical review of the mechanisms that have been proposed to explain them, and he indicates that mass-loss variations associated with solar-like magnetic cycles are perhaps the best alternative for their origin.

In a more recent paper, Simis, Icke & Dominik (2001) discuss a new non-magnetic alternative, and present detailed 1D hydrodynamical simulations for the acceleration of a dusty AGB wind. They include dust grain nucleation and growth in the AGB atmosphere, and the dust and gas components are treated as two separated fluids, without assumptions about grain drift or dust-gas coupling. Recently formed grains are accelerated first and then, as they drift throught the gas at relatively large speeds (sometimes well above the equilibrium drift velocity), the drag force transmits the momentum to the gas. For non-magnetic dusty flows, and depending on grain properties (i. e., chemical composition, size, and charge) and local conditions (i. e., ionization fraction and temperature), an efficient momentum transfer requires that the grains sweep throught gas column densities of the order of $`N10^{19}a_5`$ cm-2 (where $`a_5`$ is the grain radius in units of $`10^5`$ cm; Franco et al. 1991). Thus, the grain drift velocity can be larger than its equilibrium value for a certain distance, compressing the gas at some locations and creating regions with larger dust-to-gas mass ratios. These compressions appear in a cyclic manner in their 1D simulations and, if such a flow configuration is stable in two or three dimensions, they can drive a variable mass loss rate and may lead to the creation of dusty shells. The stability of this process, however, is difficult to explore at the moment, and it is unclear if the compressions can truly create shells or if they only lead to inhomogeneities in these radiatively driven outflows. Thus, at present, one can consider that dusty flow oscillations (if they really exist) or solar-like magnetic cycles are two possible candidates to generate the HST rings.

In Soker’s (2000) view, the magnetic field plays no direct role in the evolution of the AGB wind, however temporal variation in the number of magnetic spots would be able to modify the mass-loss rate. The number of cool spots over the AGB surface, which could be preferred sites for dust formation, is controlled by the magnetic cycle. Thus, given that the mass-loss is driven by radiation pressure on dust grains, the same cycle may also regulate periodic variations in the mass-loss rate. In this interpretation, the magnetic field is a passive player with no dynamical effects, and the mass-loss rate simply follows the spot cycle activity. In addition, to make a logical association with bipolar PNe, he suggests that a stronger magnetic activity could be expected from dynamo amplification in binary systems.

Here we develop a different point of view and explore some of the actual dynamical effects of a solar-like magnetic cycle. The possibility of a solar-like magnetic dynamo at the AGB phase has been recently discussed by Blackman et al. (2001), and they conclude that dynamo amplification is likely to operate in rotating AGB stars. A logical extension of this result is that a solar-like activity, including dynamo and cycles, is also expected in some AGB stars. We thus build a very simple model, without mass-loss variations, and perform 2(1/2)D magnetohydrodynamic computations considering a cyclic polarity inversion of the surface magnetic field of an AGB star. Our results show the importance of MHD effects in the formation of the HST rings and successfully reproduce their main features. This indicates that modulated mass-loss episodes are not really necessary to generate the rings. In Section 2 we describe the MHD models and results. A brief discussion is given in Section 3.

## 2 Numerical Models

### 2.1 The Method

The simulations have been performed using the magnetohydrodynamic code ZEUS-3D (version 3.4), developed by M. L. Norman and the Laboratory for Computational Astrophysics. This is a finite-difference, fully explicit, Eulerian code descended from the code described in Stone & Norman (1992). A method of characteristics is used to compute magnetic fields, as described in Clarke (1996), and flux freezing is assumed in all the runs. We have used spherical polar coordinates ($`r,\theta ,\mathrm{\Phi }`$), with reflecting boundary conditions at the equator and the polar axis. Rotational symmetry is assumed with respect to the rotational (polar) axis, and our models are effectively two-dimensional. The simulations are carried out in the meridional ($`r`$, $`\theta `$) plane, but three independent components of the velocity and magnetic field are computed (i. e., the simulations are “two and a half” dimensions). We are unable to include, as Simis et al. (2001) have done with their two-fluid code, the effects of dust in our simulations.

Our grids consist of $`200\times 180`$ equidistant zones in $`r`$ and $`\theta `$, respectively (with a radial extent of 0.1 pc, and an angular extent of $`90^{}`$), and the innermost radial zone lies at $`r=2.5\times 10^3`$pc from the central star. These values are used in all the simulations shown in Figures 1 – 3.

The equations for the stellar wind flow are simple, since we model an isothermal, slow AGB wind with constant mass-loss rate ($`10^6`$ M yr-1 ) and radial velocity (10  km s-1 ). These values represent the boundary conditions used in the first five innermost radial zones. Rotational effects are not important in these outflows, since we have used a small value for the stellar rotation velocity, 0.01  km s-1 (see García-Segura et al. 1999 for more details).

The novel aspect in this paper is the simple treatment of the stellar magnetic field ($`B_\mathrm{s}`$), which it is allowed to change sign in a cycle of the form:

$$B_\mathrm{s}(t)=B_{\mathrm{max}}\mathrm{cos}(2\pi \frac{t}{P}),$$
(1)
where $`B_{\mathrm{max}}`$ is the maximum average $`B`$-field at the AGB surface, and $`P`$ is the period of the magnetic cycle. Since we do not know the true variation form of the field, this functional form is just a first simple approximation. As in the case of the Sun, we assume that $`B_{\mathrm{max}}`$ has opposite signs at each hemisphere, with a neutral current sheet near the equatorial plane (e. g. Wilcox & Ness 1965; Smith, Tsurutani & Rosenberg 1978). Its average thickness in the solar case is of about $`10^8`$ cm, and its presence does not affect the field outside the equatorial sections. For simplicity, given that we compute only one hemisphere, we neglet the size of this current sheet. The equation for the inner boundary wind toroidal field on the computational domain, which is restricted to one hemisphere ($`0\theta <\pi /2`$), is:

$$B_\varphi (t)=B_\mathrm{s}(t)\frac{v_{\mathrm{rot}}}{v_{\mathrm{}}}\left(\frac{R_\mathrm{s}}{r}\right)^2\left(\frac{r}{R_\mathrm{s}}1\right)\mathrm{sin}\theta ,$$
(2)
where $`v_{\mathrm{rot}}`$ is the stellar rotation velocity, $`v_{\mathrm{}}`$ the wind velocity, and $`R_\mathrm{s}`$ the stellar radius. The function $`\mathrm{sin}\theta `$ cancels the toroidal component at the symmetry axis, $`\theta =0`$ (pole). The poloidal field component can be neglected at large distances, so that our field configuration naturally satisfies the condition $`B=0`$.

The parameter which identifies our models is the ratio of the magnetic field energy density to the kinetic energy density in the wind (Begelman & Li 1992)

$$\sigma (t)=\frac{B^2}{4\pi \rho v_{\mathrm{}}^2}=\frac{B_\mathrm{s}^2(t)R_\mathrm{s}^2}{\dot{M}v_{\mathrm{}}}\left(\frac{v_{\mathrm{rot}}}{v_{\mathrm{}}}\right)^2.$$
(3)
Obviously, this parameter is always positive and oscillates with twice the frequency of the stellar magnetic cycle.

### 2.2 Results

The first model, shown in Figure 1, corresponds to the case of a slow and dense wind ($`v_{\mathrm{}}=10\text{ km s}\text{-1}`$ and $`\dot{M}=10^6\text{M}\text{}\text{ yr}\text{-1}`$), in which the magnetic field at the surface of the AGB star cycles with a full period of 2000 yr. A peak value $`\sigma _{\mathrm{max}}=0.01`$ is used for a slowly rotating AGB star ($`v_{\mathrm{rot}}=0.01`$  km s-1 ), which is achieved, for example, when the star has $`B_\mathrm{s}=53`$ G and $`R_\mathrm{s}=1`$ AU, or $`B_\mathrm{s}=113`$ G and $`R_\mathrm{s}=100`$ R. For the cases discussed by Soker (2000), $`R_\mathrm{s}=2`$ AU and $`0.2\text{ km s}\text{-1}v_{\mathrm{rot}}2\text{ km s}\text{-1}`$, the required surface magnetic fields are within 1.3 G $`B_\mathrm{s}`$ 0.13 G. These field strenght values are modest compared withn the ones that could be achieved by dynamo amplification during AGB evolution (Blackman et al. 2001).

Figure 1 displays radial cuts, on the equatorial plane, of some of the wind variables at $`t=10,000`$ yr after the onset of the magnetized wind. The first (top) panel shows the existence of equidistant concentric shells, separated by 1,000 yr, a half of the magnetic cycle period. Each density peak corresponds to a minimum in the magnetic pressure (second panel), which also occurs at the moment when the magnetic field changes polarity (fourth panel). The last panel of Figure 1 displays the radial velocity, which clearly shows the gas response to the local variations in the total pressure (third panel) of the wind.

Thus, shell formation in a magnetized wind with variable field strength is a straightforward process. The outflowing plasma notices the magnetic pressure depressions, both upstream and downstream in the wind, and moves towards the low-pressure sites. The wind is then compressed at these locations, increasing the local density, to compensate for the low magnetic pressure values. Obviously, the particular values achieved at these density peaks, and the density constrast between ring and inter-ring zones, depend on the model assumptions. Thus, these values can be modified by changing the amplitude of the pressure fluctuations (i. e., by changing the maximum field intensities, or adding mass-loss or velocity variations, etc.) but, for our purpouses, it is sufficient to show that a reasonable density constrast is achieved by this simple model. The total pressure plot (third panel), shows that the local pressure fluctuations decrease with time (or position in the wind). As the plasma flows from the positions of the magnetic peaks and compresses the gas of the valleys, a series of MHD waves are continuosly driven that maintains small gas oscillations in the expanding wind.

Figure 2 shows the emission measure of four different models, with $`\sigma _{\mathrm{max}}=0.001,0.01,0.05,0.1`$. The plots nicely reproduce the typical spacing of the rings, and show the general drop-off expected for a constant mass-loss wind. The structures are in reality a set of spherical or quasi-spherical hemispheric shells, and their projected maximum column densities in the plane of the sky give the impression of concentric rings. For large enough values of $`\sigma _{\mathrm{max}}`$, the magnetized AGB wind is able to self-collimate towards the polar axis. This self-collimation is also present in the free streaming fast winds that form some PNe (see García-Segura et al. 1999) but, given the large densities of the slow AGB winds, the self-collimated structures are particularly prominent, and easily observable features in this case. The two right panels in Figure 2 clearly show the resulting structures, that are better defined and more conspicous in the $`\sigma _{\mathrm{max}}=0.1`$ case. The self-collimated gaseous structures also follow the periodic variations of the magnetic field, resulting in two strings of regularly spaced blobs located along the polar axis resembling those observed in He 2-90 (Sahai & Nyman 2000).

## 3 Discussion and Summary

A relevant point for the applicability of this model is whether or not the $`B`$-field can be considered frozen into the outflowing dusty AGB wind, and if the simple MHD effects that we have described are truly operative during wind evolution. This question has been previously explored by several authors and in much more restrictive environments. For instance, in collapsing molecular clouds with ionization fractions as low as $`10^7`$, where ion drift could lead to substantial magnetic flux leakage. Nontheless, even at these very low ionization fraction values, numerical models of cloud collapse show that the field remains nearly frozen into the cloud for long periods of time (e. g. Black & Scott 1982). For AGB winds the situation is certainly more favorable because the ions are provided by a number of species with low ionization potentials (Habing 1996; Cox 1997), and charged dust grains are also well coupled to the magnetic field of the outflow. Thus, flux leakage is not considered important during the AGB phase, which occurs during time scales of the order of 103 to 105 yr (for cases where ambipolar diffusion is important see Mac Low et al. 1995 and references therein). Later on, the AGB wind is photoionized by the central star of the PN (as in NGC 6543 or Hb 5), and the temperature and ionization fraction of the outflowing plasma is suddenly increased. Some pressure readjustments occur within the photoionized plasma at this time, but the basic shell structuring remains unchanged. Obviously, flux leakage becomes irrelevant after this moment, and the rings can survive for longer periods of time.

An important issue related with the validity of the present results is whether or not magnetic reconnection can occur in the flow, modifiying some of the features appearing in the models. As stated before, in analogy with the solar case, it is likely that the general AGB $`B`$-field has a dipolar structure, and a neutral current sheet should be present around the equatorial plane of their magnetized winds. The process creates a pair of half-shells at each cycle, one per hemisphere and separated by the neutral current sheet at the equatorial plane. Each half-shell, in turn, is compressed by magnetic intershell regions with opposite polarities. This implies that the compressed half shells also act as neutral current sheets, separating regions with oppositely directed field lines. Our simulations assume well ordered fields and the results do not show, even in the more magnetic models, any signs of reconnection within the shells. If one assumes that a twisted random field component could be also present, a faster reconnection mode could also occur (e. g. Lazarian & Vishniac 1999). Given the low field values in the shells, however, any random field reconnection that may occur within them is severily limited by the corresponding low Alvén speeds, and the energy and flux involved in such a process has to be neglegible.

Another related question is if the shells can be formed only by modulated mass-loss rate episodes, without magnetic pressure in the wind, as proposed by Soker (2000) and Simis et al. (2001) (keeping in mind that the Simis et al. results could be unstable in 2D or 3D, precluding the actual formation of a shell). In these cases, the periodic formation of shells in the expanding wind is only due to the increase in mass-loss, and the higher shell densities can be maintained as long as the wind temperature decreases in these same locations. This is a perfectly reasonable possibility during the AGB phase. However, once the wind becomes photoionized, the plasma temperature becomes fairly homogeneous, and the rings tend to be washed away in a sound crossing time. Thus, ring survival is severely compromised in these purely hydrodynamic cases.

There is one additional observational element that suggests that the magnetohydrodynamic mechanism is indeed relevant in these cases. Recent HST images of He 2-90 show a series of brigth knots which have been interpreted as due to symmetric and well collimated fast jets (Sahai & Nyman 2000). New spectroscopic observations by Guerrero et al. (2001) have found that these jet-like features are moving away from He 2-90 at a remarkably constant radial velocity of only 26  km s-1 . Obviously, this radial speed value has to be corrected for the inclination angle with respect to the plane of the sky. This angle is unknown, but assuming that the axis of the jet-like features is located almost on the plane of the sky, they find an upper bound to the actual space velocity of about 290  km s-1 . However, the lack of any trace of collisional excitacion within the flow or interaction with the ambient medium (i. e., the sulfur and oxygen lines in the collimated structures are weak, and there are no signs of the bow-shocks expected for the leading parts of supersonic jets), makes it difficult to ascribe these highly collimated structures with a supersonic jet (this is reinforced by the observed narrow line profiles and lack of acceleration within the structures). In contrast, the constant and moderate velocitie values, and the similarity of the observed knotty structure with the largest $`\sigma _{\mathrm{max}}`$ model presented here provides a very atractive alternative for its origin. Figure 3 shows a qualitative comparison of one of the models ($`\sigma =0.1`$) with the structuring of the “jets” in He 2-90 (figure 3 in Sahai & Nyman 2000). As stated above, the peak densities depend on model assumptions and, despite the fact that we have not attempted to specifically reproduce this nebula, the similarity is outstanding. Thus, although the images actually show features that resemble supersonic, fast jets, their observed properties, and our model, indicate that these features could have a different nature. If this is the case, the velocity of the knots should be close to that of the AGB wind, and the low radial velocity observed by Guerrero et al. (2001) is compatible with our model.

Since jets are actually formed by magnetic collimation, the jet interpretation of Sahai & Nyman (2000) and Guerero et al. (2001) has a logical link with our results, but there is one important difference. This key difference is that the gas in our case is $`only`$ collimated but $`not`$ accelerated. The collimation of these outflows is solely due to the hoop stress of the toroidal field, but their acceleration require an additional process. For instance, the acceleration of a jet-like flow in some planetaries is achieved after the magnetized wind passes throught a reverse shock and then is subsequently pinched at the polar axis (Różyczka & Franco 1996). The case of He 2-90 is probably different since, although magnetic collimation is most likely present, the acceleration mechanism to create a “supersonic fast jet” seems to be absent. Note that collimation without acceleration can occur in a freely outflowing magnetized wind.

It is also relevant to ask if the shells formed in our models are stable against the 2D undular or interchange modes, or the 3D mixed mode, of the Parker instability (the 3D mixed mode is a combination of the interchange and undular modes; see Hughes & Cattaneo 1987, and Matsumoto et al. 1993). The growth rates of these instabilities depend on the field strenght, the acceleration of the magnetized plasma, and size of the magnetically supported regions, but the fastest growing mode corresponds to the interchange instability (in this case, the flux tubes are displaced in the radial direction, following the pressure gradient). The undular instability, on the other hand, grows at a slower pace and distorts the $`B`$-field lines, creating a wavy structure along the field lines. For the model configurations considered here, the simulations are blind to the undular instability (because this branch grows in a plane perpendicular to the computational domain; see Kim et al. 2000) but they should be able to resolve some modes (i. e., those with wavelengths larger than the grid size) of the interchange instability. The accelerations of the flow, however, are very small and the resulting growth time scales are much too large. Thus, for any practical purpouses, the models are effectively stable to these perturbances.

In addition, given that the higher pressure intershell regions become ligther than the shells, the growth of plain Rayleigh-Taylor instabilities in the radial direction could become relevant at some stage. Our simulations cannot resolve the fine scales of this instability, but a simple analysis indicates that the growth times are also large. One cannot define a sharp interface between the (magnetic) intershell and (non-magnetic) shell regions (their initial densities are equal), but the maximum pressure difference is simply equal to $`B_{\mathrm{max}}^2/8\pi `$. Thus, for a shell thickness $`L`$, an upper limit to the pressure gradient in the radial direction is simply $`B_{\mathrm{max}}^2/8\pi L`$, and the corresponding upper limit to the acceleration during shell formation is $`av_a^2/2L`$ (where $`v_a`$ is the average Alfvén speed inside the shell). The resulting growth time for R-T instabilities is then (Shore 1992),

$$\tau \left(\frac{\rho _2+\rho _1}{\rho _2\rho _1}\right)^{1/2}\left(\frac{L}{a}\right)^{1/2}\left(\frac{\rho _2+\rho _1}{\rho _2\rho _1}\right)^{1/2}\left(\frac{2^{1/2}L}{v_a}\right),$$
(4)
where $`\rho _2`$ is the shell density, and $`\rho _1`$ is the intershell density. Clearly, this time scale is always large: at the early stages because the density contrast is very small, and because the acceleration drops as the density difference grows at later stages. Thus, the growth time of RT instabilities is always larger than the time scale for shell formation.

In summary, the dynamical effects of a solar-like magnetic cycle operating in late AGB stars have been explored with MHD models. Although a number of simplifying assumptions have been made, these 2(1/2)D numerical simulations are able to succesfuly reproduce for the first time the main observed properties of the HST rings. One obvious limitation of our models is that we are unable to include the dynamical effects of dust grains, and cannot explore the magnetic counterpart of the Simis et al. (2001) models. Further studies with less restrictive assumptions may shed more light on the range of MHD effects that may operate in these type of objects.

Acknowledgments It is a pleasure to thank Yervant Terzian for his helpful discussions and encouragement about this topic during his visit to Ensenada. We also thank Jack Thomas for pointing out the recent work on dynamo amplification in AGB stars. We warmly thank M. Guerrero, L.F. Miranda, Y-H. Chu, M. Rodriguez and R.M. Williams for sharing their results with us prior to publication. The comments and criticisms of the referee are also gratefully acknowledged. As usual, we also thank Michael L. Norman and the Laboratory for Computational Astrophysics for the use of ZEUS-3D. The computations were performed at Instituto de Astronomía-UNAM. This work has been partially supported by grants from DGAPA-UNAM (IN130698, IN117799 & IN114199) and CONACyT (32214-E).

List of References

Balick, B., Wilson, J. & Hajian, A. R. 2000, AJ, 121, 354

Begelman, M. & Li, Z. Y. 1992, ApJ, 397, 187

Black, D. C. & Scott, E. H. 1982, ApJ, 263, 696

Blackman, E. G., Frank, A., Markiel, J. A., Thomas, J. H. & Van Horn, H. M. 2001, Nature, 409, 485

Clarke, D. A. 1996, ApJ, 457, 291

Cox, P. 1997, in IAU Symp. No 180 “Planetary Nebulae”, eds. H. J. Habing & H. J. G. L. M. Lamers, (Kluwer Academic Publishers, Dordrecht), 139

Franco, J., Ferrini, F., Ferrara, A. & Barsella, B. 1991, ApJ, 366, 443

Habing, H. J. 1996, A&A Rev., 7, 97

Hrivnak, B. J., Kwok, S. & Su, K. Y. L. 2001, AJ, 121, 2775

Hughes, D. W. & Cattaneo, F. 1987, Geophys. Astrophys. Fluid Dyn., 39, 65

García-Segura, G., Langer, N., Różyczka, M., & Franco, J. 1999, ApJ, 517, 767

Guerrero, M., Miranda, L.F., Chu, Y-H., Rodríguez, M. & Williams, R.M., 2001, pre-print.

Kim, J., Franco, J., Hong, S. S., Santillán, A. & Martos, M. 2000, ApJ, 531, 873

Kwok, S., Su, K. Y. L., Hrivnak, B. J. 1998, ApJ, 501, L117

Lazarian, A. & Vishniac, E. 1999, ApJ, 517, 700

Mac Low, M.-M., Norman, M. L., Konigl, A., & Wardle, M. 1995, ApJ, 442, 726

Matsumoto, R., Tajima, T., Shibata, K. & Kaising, M. 1993, ApJ, 414, 357

Różyczka, M. & Franco, J. 1996, ApJL, 469, L127

Sahai, R., & Nyman, L.- A. 2000, ApJ, 538, L145

Shore, S. N. 1992, An Introduction to Astrophysical Hydrodynamics, (Academic Press, London)

Smith, E.J., Tsurutani, B.T., & Rosenberg, R.L. 1978, J. Geophys. Res., 83, 717

Soker, N. 2000, ApJ, 540, 436

Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753

Terzian, Y., & Hajian, A. R. 2000, Asymmetrical Planetary Nebulae II: From Origins to Microstructures, ed. J. H. Kastner, N. Soker, S. A. Rappaport, A.S.P. Conf. Series, 199, 33

Wilcox, J., & Ness, N. 1965, J. Geophys. Res., 70, 1233

Figure Captions

Figure 1. Equatorial cut for the run with $`\sigma =0.01`$ showing the effects of a cyclic polarity inversion of the AGB stellar magnetic field. Each density peak correspond to 1,000 yr, half of the full period. Units are in c.g.s., except for the radial velocity ( km s-1 ).

Figure 2. Emission measure of four different runs projected with a tilt of $`45^{}`$ from the plane of the sky. The panels correspond to: (top-left) $`\sigma =0.001`$, (bottom-left) $`\sigma =0.01`$, (top-right) $`\sigma =0.05`$, (bottom-right) $`\sigma =0.1`$. Each concentric shell correspond to 1,000 yr. Note the piled-up blobs at the polar axis for the most magnetized runs, resulting from the self-collimation of the AGB wind.

Figure 3. Comparison between the polar region of the model with $`\sigma =0.1`$ (projected without tilt on the plane of the sky) (top), with figure 3 by Sahai & Nyman (2000) of He 2-90 (bottom). In their figure, the relative intensity of the jet is plotted as a function of radius showing the radial offsets of the knots from the center (thick curve: southeast jet, thin curve: northwest jet). Both plots are logarithmic in the Y axis, and linear in the X axis.


text:
# 1 Introduction

## 1 Introduction

The distance to the Large Magellanic Cloud is not agreed upon. A recent compilation of the values of distance modulus published in 1998 and 1999 (Gibson 1999) cover a range from 18.07 to 18.74, which corresponds to a factor 1.36 in the distance. An attempt to harmonize the distance scale based on RR Lyrae variables and red clump giants, the two best calibrated indicators, gives a range 18.24 - 18.44 (Popowski 2001). I think the current empirical calibration of the red clump giants is the most reliable, and this gives $`18.24\pm 0.08`$ (Udalski 2000), but the issue is far from being settled.

In the near future the most reliable results will be obtained with analysis of detached eclipsing binaries. A very good description of method, and many historical references dating from 1910, can be found in Kruszewski and Semeniuk (1999). The method was applied by Lacy (1979) to several dozens nearby systems, and the distances he obtained have been verified with the Hipparcos parallaxes (Ribas et al. 1998, Semeniuk 2000). A much better calibration will be provided by the future Full-sky Astrometric Mapping Explorer (FAME, Horner et al. 1999, Semeniuk 2001). However, even today this method provides the LMC distance which is competitive with the best alternatives. The first modern determination of the LMC distance modulus based on eclipsing binaries were by Bell et al. (1991) and Bell et al. (1993), who used HV2226 and HV5936 to obtain $`18.6\pm 0.3`$ and $`18.1\pm 0.3`$, respectively. More recently HV 2274 was observed by Guinan and his collaborators, and they obtained for the distance modulus values $`18.54\pm 0.08`$, $`18.42\pm 0.07`$, and $`18.30\pm 0.07`$ (Guinan et al. 1997, 1998a, 1998b, respectively). The same data but different analysis of the interstellar reddening provided somewhat different values: $`18.22\pm 0.13`$ (Udalski et al. 1998b) and $`18.40\pm 0.07`$ (Nelson et al. 2000). Finally, Fitzpatrick et al. (2000) obtained $`18.31\pm 0.09`$ using HV982. While the agreement between various determinations is not perfect, they all are in a reasonable agreement with the ‘short’ distance scale to the LMC.

While I expect that there will be a major improvement in the distance determinations using detached eclipsing binaries, and the errors in the distance modulus will be reduced to 0.05 mag, or even less, it is far from clear that the proponents of a different distance value will be convinced. Ultimately, it would be important to obtain a 1% distance using some purely geometric method. This will not be easy, as even if the Space Interferometry Mission (SIM, Allen et al. 1997) achieves astrometric accuracy of $`1\mu s`$, this will translate into 5% accuracy for the distance determination. While in principle the accuracy could be improved by measuring parallaxes to many stars in the LMC, it is not likely that systematic errors will be reduced so much as to make it practical.

In the following section I present an outline of a very traditional and purely geometrical method for distance determination based on a combination of astrometry and spectroscopy of visual binaries. In the subsequent section I present a list of candidate objects, selected from the published OGLE catalogs of BVI photometry for over 7 million stars in the LMC (Udalski et al. 2000), and over 2 million stars in the SMC (Udalski et al. 1998a). Finally, I outline an observational approach to the selection of the best systems.

A similar approach to the distance determination was proposed by Massa and Endal (1987ab).

## 2 Geometric Distance to the Magellanic Clouds

The simplest and best known geometrical method for distance determination is the traditional parallax. Unfortunately, at the LMC distance even SIM will be able to perform 5% measurement, at best. Fortunately, there is another, almost equally traditional method, based on visual binaries. Combining astrometric and spectroscopic orbits provides 1% distances and masses. Two good examples are $`\iota `$ Pegasi (Boden et al. 1999a) and 64 Piscium (Boden et al. 1999b). This method is likely to settle the controversy over Pleiades distance by determining the orbit of Atlas (Pan et al. 1999).

A great advantage of visual binary method over the parallax is that only small angle relative astrometry is needed. Also, binary orbits may be much larger than 1 AU, making it possible to reach farther with excellent accuracy. One problem is the need to resolve the binary, at least in the conventional use of the binary method. This is a serious problem at the LMC distance of $`50`$ kpc. For a binary in a circular orbit the angular separation is given as

$$\theta =0.25mas\left(\frac{P_{orb}}{10yr}\right)^{2/3}\left(\frac{M_1+M_2}{20M_{}}\right)^{1/3},\mathrm{ford}=50\mathrm{kpc}.$$
$`(1)`$
An orbital period of $`10`$ years is about as long as acceptable, considering rapid progress of technology. The binary HV2274 has the total mass of $`23.5M_{}`$ Guinan et al. 1998b), and the apparent magnitude $`V=14.16`$ (Udalski et al. 1998b), hence the mass scaling adopted in the eq. (1) is reasonable. There are two problems with an application of visual binary method to LMC distance determination: 0.25 mas is too small an angle, and at $`V14`$ mag a binary is too faint for current optical and infrared interferometers. However, it is not unreasonable to expect that within the next 5 or 10 years such a binary will be within reach of future ground based instruments.

There is another possibility: to use SIM do perform 1 micro-arc-second astrometry on the light centroid of a visual binary. With the expected $`10`$ mas angular resolution SIM could not resolve the binary but it could determine the motion of the light centroid with an accuracy better than 1%, if the amplitude of centroid motion is 0.25 mas. To make this useful we have to select visual binaries made of two diverse components: one very hot and blue, the other very cool and red, and to take advantage of the broad spectral response of the SIM, with many filters covering wavelengths in the range $`0.40.8`$ microns, approximately. The location of the light centroid would be close to the hot star in the blue, and close to the cool star in the near infrared. Of course, to be feasible, this approach requires us to decompose the binary spectrum into two components, and to determine the luminosity ratio of the two stars as a function of wavelength. While this is not an easy task, there is no obvious reason for it to be impossible to accomplish. In fact it should be done from the ground prior to any attempt to put a candidate binary on the SIM target list.

Note, that Hipparcos discovered many new visual unresolved binaries by measuring periodic variations in the position of their light centroid (Mignard 1998, and references therein).

For any of the two proposed approaches: either the resolution of an LMC binary with a future ground based interferometer, or the astrometry of the light centroid variations with the future SIM, it is necessary to identify suitable binaries and to determine accurate spectroscopic orbits. With the target stars brighter than 14 mag the determination of accurate radial velocity curve should not be a problem. The candidate stars are selected from the OGLE catalogs of millions of stars in the next section.

## 3 Candidate Binaries

Udalski et al. (1998a, 2000) published BVI photometry for over 2 million SMC stars, and over 7 million LMC stars, respectively. The results for the SMC and LMC are available over the Internet at:

http://bulge.princeton.edu/~ogle/ogle2/bvi\_maps.html

http://bulge.princeton.edu/~ogle/ogle2/lmc\_maps.html

and at:

http://www.astrouw.edu.pl/~ogle/ogle2/bvi\_maps.html

http://www.astrouw.edu.pl/~ogle/ogle2/lmc\_maps.html

The stars selected for analysis were brighter than $`V=14`$ mag, and had good photometry in all three bands: B, V, and I. There were 810 such stars in the SMC, and 1127 in the LMC. They are shown in Fig. 1 and Fig. 2 in $`(VI)(BV)`$ color-color diagrams, and in Fig. 3 and Fig. 4 in color-magnitude diagrams. The diagonal lines in Fig. 1 and Fig. 2 correspond to the relation:

$$(BV)=\frac{7}{6}(VI)0.5.$$
$`(2)`$
The open circles located below this line are the stars which have an excess of light in the B and I bands, i.e. these are the candidates for unresolved blue-red binaries. They are all listed in Table 1, in which the second columns give the OGLE field number, and the third columns the star number within the corresponding fields. All other columns are self-explanatory. Finding charts for all the candidate binaries are shown in Fig. 5 and Fig. 6. The candidate is located at the center of each square panel, which is 1,000 pixels on its side, corresponding to 7’ in the sky.

## 4 Discussion

The candidates listed in Table 1 have to be verified spectroscopically for the presence of composite spectra. The next step is more time consuming: for every star with a composite spectrum the two radial velocity curves have to be determined. Some genuine binaries will have periods too short, and hence angular separations too small to be of interest. Others will have periods too long to be of interest. Hopefully, there will be several objects in the period range of 5 – 15 years, which might be considered to be optimal for distance determination.

The list of candidate binaries will increase as OGLE upgrades to a large mosaic CCD camera, and multi-band photometry will become available for full SMC and LMC, not just for their central regions. Also, photometry for the stars which are brighter than the saturation limit apparent in Fig. 3 and Fig. 4 should be searched for possible candidates.

Acknowledgments. This work was supported with the NSF grants AST-9819787 and AST-9820314. It is a great pleasure to acknowledge that all figures were made using the SM plotting package provided by R. H. Lupton.

## REFERENCES

* Allen, R., Shao, M., and Patterson, D. 1997, Proc. SPIE, 2871, 504.
* Bell, S. A. et al. 1991, MNRAS, 250, 119.
* Bell, S. A. et al. 1993, MNRAS, 265, 1047.
* Boden, A. F. et al. 1999a, Astrophys. J., 515, 356.
* Boden, A. F. et al. 1999b, Astrophys. J., 527, 360.
* Fitzpatrick, E. L. et al. 2000, astro-ph/0010526.
* Gibson, B. K. 1999, astro-ph/9910574.
* Guinan, E. F. et al 1997, Bull. Am. Astron. Soc., 29, 1209.
* Guinan, E. F. et al. 1998a, astro-ph/9809132v1.
* Guinan, E. F. et al. 1998b, Astrophys. J. Letters, 509, L21.
* Horner, S. D. et al. 1991, astro-ph/9907213.
* Lacy, C. H. 1979, Astrophys. J., 228, 817.
* Kruszewski, A., and Semeniuk, I. 1999, Acta Astron., 49, 561.
* Massa, D., and Endal, A. S. 1987a, AJ, 93, 579.
* Massa, D., and Endal, A. S. 1987b, AJ, 93, 760.
* Mignard, F. 1998, in: Highlights of Astronomy, Vol. 11A, JD 14 of the XXIIIrd General Assembly of the IAU, Ed. J. Andersen. Kluwer Academic Publishers, 1998, p. 539.
* Nelson, C. A. et al. 2000, Astron. J., 119, 1205.
* Pan, X. P., Shao, M., and Kulkarni, S. 1999, Am. Astron. Soc. Meeting, 195, 132.03.
* Popowski, P. 2001, MNRAS, 321, 502.
* Ribas, I. et al. 1998, Astron. Astrophys., 330, 600.
* Semeniuk, I. 2000, Acta Astron., 50, 381.
* Semeniuk, I. 2001, astro-ph/0103511.
* Udalski, A. 2000, Astrophys. J. Letters, 531, L25.
* Udalski, A. et al. 1998b, Astrophys. J. Letters, 509, L25.
* Udalski, A. et al. 1998a, Acta Astron., 48, 147.
* Udalski, A. et al. 2000, Acta Astron., 50, 307.


text:
# Effect of electron irradiation on vortex dynamics in YBa₂Cu₃O_{7-𝛿} single crystals

\[

## Abstract

We report on drastic change of vortex dynamics with increase of quenched disorder: for rather weak disorder we found a single vortex creep regime, which we attribute to a Bragg-glass phase, while for enhanced disorder we found an increase of both the depinning current and activation energy with magnetic field, which we attribute to entangled vortex phase. We also found that introduction of additional defects always increases the depinning current, but it increases activation energy only for elastic vortex creep, while it decreases activation energy for plastic vortex creep.

\]

The effect of random pinning on crystalline order and on dynamics of the flux-line-lattice (FLL) was a subject of numerous experimental and theoretical investigations. Neutron diffraction  and $`\mu `$SR  experiments gave experimental evidence for the existence of two vortex solid phases with different positional correlations. Magnetization measurements  of $`Bi_2Sr_2CaCu_2O_8`$ showed that the field, $`H_{on}`$, at which a steep increase of the measured current $`J_m(H)`$ begins, coincides approximately with the field above which the intensity of Bragg peaks sharply decreases . The field $`H_{on}`$ was interpreted as a phase boundary between low field ordered and high field disordered vortex phases. These experimental results are supported by theoretical studies. It was shown that in presence of rather weak disorder the FLL retains a quasi-long-range order resulting in the so-called Bragg glass phase . However, with the increase of random pinning or magnetic field a transition to strongly disordered entangled vortex phase (glass phase) is predicted . A sharp increase in magnetization below the fish-tail peak position $`H_p`$ was observed in $`Nd_{1.85}Ce_{0.15}CuO_{4d}`$  and non twinned $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$  single crystals.

Magnetic measurements showed that in optimally doped $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ crystals no fishtail behavior is observed both in detwinned  and twinned  samples, while decrease of the oxygen content always induced non monotonous $`J_m(H)`$-curves. Two distinctive peculiarities in the magnetization curves of oxygen deficient  and electron irradiated  crystals were observed: (1) the peaks $`H_{on}`$ and $`H_p`$ shift toward lower magnetic fields with increasing defect concentration, and (2) in magnetic fields $`H<H_p`$ the current $`J_m`$ increases, while in magnetic fields $`H>H_p`$ the current $`J_m`$ decreases with increasing defect concentration. It is believed  that the peak $`H_p`$ separates elastic vortex creep in low and plastic vortex creep, mediated by motion of the FLL dislocations, in high magnetic fields. Thus the introduction of additional defects leads to an increase of vortex pinning in the region of elastic creep and such behavior is expected. On the other hand, decrease of pinning force with increasing disorder observed in the region of plastic creep is non trivial, and reasons of such behavior are not known. The aim of this paper is to show the effect of point-like defects concentration on vortex dynamics and pinning parameters in $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ single crystals.

The investigated sample was $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ single crystal with $`T_c`$ $``$ 93.5 K and $`\mathrm{\Delta }T_c<`$ 0.5 K. Twin planes (TP’s) inside the measured part of the sample were aligned in one direction. The transport current was applied along the $`ab`$-plane and at an angle $`\alpha =45^{}`$ with respect to TP’s. Measurements were performed in magnetic fields applied parallel to the $`c`$-axis. Temperature stability of the sample during measurements was better than 5 mK, and measurements in the normal state showed that overheating of the sample at the highest dissipation level of 50 $`\mu `$W did not exceed 10 mK.

Additional defects were introduced by irradiation with 2.5-MeV electrons, which are suitable for production of point-like defects. Irradiation was performed at temperatures $`T10K`$  and after irradiation the current-voltage characteristics (CVC) were measured without heating the sample above 110 K. This excludes diffusion, and therefore annihilation and clustering of the defects. Irradiation was performed at dose rate $`4.2\times 10^{13}cm^2sec^1`$ and at an angle $`5^{}`$ off the c-axis to avoid electron channeling. Homogeneity of electron beam was about 5%. Following procedure used in Ref.  we estimated that irradiation dose $`10^{18}cm^2`$ produces the averaged over all sublattices concentration of the defects $`10^4`$ dpa, and that the penetration range $``$ 1mm is at least two orders higher then thickness of the crystal $`7\mu `$m resulting in homogeneity of the defects along the c-axis better then 1%. Reduction of the $`T_c`$ after irradiation was $`\mathrm{\Delta }T_{irr}=`$ 0.65, 1 and 1.3 K after irradiation with doses $`10^{18},2\times 10^{18}`$ and $`3\times 10^{18}el/cm^2`$, respectively.

Fig.1 shows field variation of the measured ”critical” current $`J_m`$ determined at an electric field level $`10^6V/cm`$. Before irradiation the current $`J_m`$ continuously decreases with increasing field as it was previously observed by Zhukov et al. in crystals with very small oxygen deficiency, $`\delta 0.03`$ . After irradiation the $`J_m(H)`$ curves show fish-tail behavior and the peak position $`H_p`$ in the $`J_m(H)`$ curves shifts toward low fields with the dose. It is also seen that introduction of additional defects increases $`J_m`$ for $`H<H_p`$, while the current $`J_m`$ decreases with the dose for $`H>H_p`$. Such variation of the $`J_m`$ and $`H_p`$ with increasing of the defect concentration is analogous to the behavior previously observed in oxygen deficient  and in electron irradiated  $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ single crystals.

In low magnetic fields ($`H15kOe`$ for non irradiated sample, and $`H<H_p`$ for irradiated samples) the CVC data follow the dependence

$$E=E_0exp[(U_0/kT)(J_c/J)^\mu ],$$
(1)
where the exponent $`\mu 1`$, $`E_0`$ is a constant, $`U_0`$ is the activation energy, and $`k`$ is the Boltzmann constant. The depinning critical current $`J_c`$ can be determined by extrapolation of the ratio $`\rho _d(J)/\rho _{BS}`$ to unity  assuming that at current density $`J=J_c`$ the differential resistivity $`\rho _ddE/dJ`$ equals the flux flow resistivity in the Bardeen-Stephen model, $`\rho _{BS}=\rho _NB/B_{c2}`$ , where $`\rho _N`$ is the normal state resistivity. Induction of the upper critical field was estimated assuming that $`B_{c2}=(dB_{c2}/dT)(TT_c)`$ with $`dB_{c2}/dT`$ = -2.5T/K . Field variation of the $`J_c`$ is shown in the inset of Fig.2a. Substituting these values of $`J_c`$ in Eq.1, and fitting experimental $`E(J)`$-curves by this equation we derived $`U_0(H)`$ dependence shown in the inset of Fig.2a.

As one can see in Fig.1 and in the inset of Fig.2a, vortex dynamics, namely field variation of the pinning parameters, strongly depends on the strength of disorder. In presence of weak disorder (before irradiation) $`J_c`$, $`U_0`$, and vortex velocity $`v=cE/B`$ do not depend on magnetic field indicating a single vortex creep regime, as it is predicted by the collective pinning theory  for low magnetic fields. Also, for this regime of vortex creep the current $`J_m`$ decreases with increasing magnetic field due to increase of the flux density. Increase of the disorder induces field variation of the pinning parameters. After irradiation $`J_c`$ and $`U_0`$ increase with magnetic field that leads to reduction of vortex velocity and to increase of the current $`J_m`$ with increasing magnetic field. These observations indicate that before and after irradiation we test different vortex phases.

In presence of weak quenched disorder the ordered Bragg-glass phase is expected  and our experimental data probably indicates that dynamics of just this vortex phase is described by the single vortex creep regime. Also note, that derived exponent $`\mu 1\pm 0.15`$ is close to the value of $`\mu 0.70.8`$ predicted for the Bragg-glass phase for not too small currents . With increase of disorder a transition of the Bragg-glass to entangled vortex phase is expected . The derived after irradiation weak increase of the $`U_0`$ and rapid increase of the $`J_m`$ with magnetic field correlates with previous experimental findings  for magnetic fields $`B>B_{on}`$, or in magnetic fields where the entangled vortex solid is expected. Therefore we believe that after irradiation we test the entangled vortex phase and dynamics of this phase is characterized by the increase of both $`J_c`$ and $`U_0`$ with increasing magnetic field, which lead to rapid increase of the current $`J_m`$. Also, introduction of additional defects increases both the $`U_0`$ and $`J_c`$ resulting in steep increase of the current $`J_m`$ with irradiation dose, see Fig.1, in agreement with previous investigations .

Let us consider experimental data obtained in magnetic fields $`H20kOe`$ for nonirradiated sample, and in magnetic fields $`H>H_p`$ for irradeated samples. As one can see in Fig.2b, the CVC data follow the dependence predicted for motion of the FLL dislocations

$$E(J)=\rho _0Jexp\{(U_{pl}/kT)[1(J_{pl}/J)^\mu ]\},$$
(2)
where $`\mu `$ = -1/2, and $`\rho _0`$ is a constant. The main peculiarities of the curves presented in Fig.2b are that the slope of the curves measured in the same magnetic field decreases with increasing defects concentration and that these curves intersect one another. For vortex creep described by Eq.2 such behavior is possible only in that case when the critical current increases, but the activation energy decreases with increasing defect concentration. Extrapolating the ratio $`\rho _d(J)/\rho _{BS}`$ to unity we derived field and temperature variation of the current $`J_{pl}`$ shown in the inset of Fig.2b and in Fig.3a, respectively. Substituting these values of $`J_{pl}`$ in Eq.2, and fitting experimental $`E(J)`$-curves by this equation we determined the field and temperature variation of $`U_{pl}`$ shown in the inset of Fig.2b and in Fig.3a, respectively.

The current $`J_{pl}`$ is determined by the interaction of the dislocations with the FLL and random pinning centers. Interaction with the FLL results in a shear limited contribution  $`J_{sh}c_{66}/BdB^{1/2}`$, where $`c_{66}=\mathrm{\Phi }_0B/(8\pi \lambda )^2`$ is the shear modulus, $`\mathrm{\Phi }_0`$ is the flux quantum, $`\lambda `$ is the penetration depth, $`da_0`$ is the width of channel for moving dislocation, and $`a_0(\mathrm{\Phi }_0/B)^{1/2}`$ is the intervortex distance. As one can see in the inset of Fig.2b, in magnetic fields not very close to the melting point the current $`J_{pl}`$ really increases with the field in agreement with theoretical predictions. In the presence of random pinning the shear modulus $`c_{66}`$ decreases  and therefore $`J_{sh}`$ also decreases. However, the contribution of the core pinning increases with the defect concentration . The derived increase of the current $`J_{pl}`$ with irradiation dose indicates that increase of the core pinning dominates over reduction of the $`J_{sh}`$.

The activation energy $`U_{pl}`$ decreases with the increase of both irradiation dose and magnetic field. A decrease of the activation energy with increasing field agrees with theoretical calculations, $`U_{pl}\epsilon \epsilon _0a_0B^{1/2}`$ , where $`\epsilon _0=(\mathrm{\Phi }_0/8\pi \lambda )^2`$, and previous experimental findings . But, a decrease of the activation energy with increasing irradiation dose is not evident. Of course, the $`T_c`$ decreases after irradiation which leads to an increase of the $`\lambda (T)=\lambda (0)/[1(T/T_c)^2]^{1/2}`$, and hence to a decrease of the activation energy $`U_{pl}\lambda ^2`$. However, the reduction of the $`T_c`$ can not describe the fast decrease of $`U_{pl}`$ in our measurements. This is demonstrated in Fig.3a, which shows dependence $`U_{pl}`$ vs. $`tT/(T_c\mathrm{\Delta }T_{irr})`$. It is evident that $`U_0`$ decreases with the dose increasing even if we take into account the reduction of $`T_c`$.

To explain this behavior let us consider a displacement of a vortex segment $`L_0`$ over intervortex distance $`a_0`$ as shown in the left-hand inset of Fig.3b. In the absence of random pinning the energy of such nucleus can be written as $`U_{pl}=2E_{el}+E_{pl}`$, where $`E_{el}\epsilon \epsilon _0a_0`$ is elastic energy of the vortex segment $`L_{ab}`$, and $`E_{pl}`$ is the energy required for displacement of the vortex segment $`L_0`$ over the distance $`a_0`$. At small driving forces such a nucleus is stable when $`E_{pl}2E_{el}`$. Thus we obtain a minimal activation energy $`U_{pl}4\epsilon \epsilon _0a_0`$, which coincides within a factor 4 with previous estimates $`U_{pl}\epsilon \epsilon _0a_0`$ . In the absence of random pinning the equilibrium positions of the vortices are the straight lines as it is shown by two parallel straight lines. In the presence of random pinning the equilibrium positions become curved as shown by dashed lines in the right-hand inset of Fig.3b. Therefore, along the vortex lines a vortex fragments, for which an average distance $`a_p`$ between two neighboring equilibrium positions of vortex lines is smaller then $`a_0`$ appear. In this case the energy $`U_{pl}4\epsilon \epsilon _0a_p`$ is smaller then activation energy in the absence of pinning because $`a_p<a_0`$. Substantial displacements of vortex segments, $`\mathrm{\Delta }a_0=a_0a_p0.2a_08nm`$, due to core interaction with individual pinning centers probably unreliable, and we attribute them to fluctuations of the defects concentration, which naturally present in real crystals. Indeed, equating the work $`A=J_{pl}F_0L_0a_p/c`$, required for displacement of the segment $`L_0`$ over distance $`a_p`$, to the elastic energy, $`2E_{el}=U_{pl}/2`$, we estimated $`L_0`$ 500 nm for $`\phi =10^{18}`$ el/cm2, and derived that for coherence length $`\xi `$(85K) = 5 nm the segment $`L_0`$ interacts with about 200 point defects. Therefore, fluctuations of the defects concentration of 10% can give difference in core interaction with about 20 point defects in different positions of the segment $`L_0`$.

TP s strongly affect pinning and dynamics of vortices . In particular, being plane defects they form channels of easy vortex motion along the plane of twins . Due to suppression of superconducting order parameter within the TP s some part of vortices is trapped by the TP s , and pinning of these vortices along the TP s may be reduced as compared with pinning in the bulk of the crystal . Therefore, for the same driving force velocity of the trapped vortices can be higher compared with velocity of vortices placed in the bulk of crystal. In high magnetic fields contribution of the trapped vortices to dissipation of energy is small due to small fraction of these vortices. However, in a magnetic field 1 kOe the intervortex separation $`a_0140nm`$ becomes comparable with the distance between twins $`d300nm`$ in our sample, and a significant part of vortices can be trapped by the TP’s. Therefore, contribution of the trapped vortices increases with decreasing magnetic field, and deviation from the field scaling for non irradiated sample in low fields presented in the inset of Fig.1 probably reflects the reduced pinning of the trapped vortices.

In conclusion, we have shown that vortex dynamics dramatically depends on the strength of disorder. In presence of weak disorder a single vortex creep is realized, which we attribute to dynamical property of the ordered Bragg-glass phase. In presence of strong disorder and for elastic vortex creep we have found rapid increase of the depinning current and weak increase of the activation energy with increasing magnetic field, which we attribute to dynamical property of the entangled vortex solid. We also found that the introduction of additional defects always increases the depinning critical current, but it increases the activation energy only for the elastic vortex creep, while it decreases the activation energy for the plastic vortex creep.

We gratefully acknowledge support by German Federal Ministry for Research and Technology under the Project UKR-032-96.