format: text

text:
# Magnetic Fields and the Cosmic Microwave Background

## Abstract

Observations of the high degree of isotropy of the cosmic microwave background are commonly believed to indicate that the Universe is ‘almost’ Friedmann-Lemaître-Robertson-Walker (at least since the time of last scattering). Theoretical support for this belief comes from the so-called Ehlers-Geren-Sachs Theorem. We show that a generalization of this theorem rules out any strong magnetic fields in the Universe. Our theoretical result is model-independent and includes the case of inhomogeneous magnetic fields, complementing previous results. We thus prove that cosmic microwave background observations severely constrain all types of primordial and protogalactic magnetic fields in the universe.

Magnetic fields seem to be abundant in the universe, and it is an open and important question as to whether the origin of these fields are either primordial (i.e. originating in the early universe and already present at the onset of structure formation), or protogalactic (i.e. generated by battery mechanisms during the initial stages of structure formation). One way to distinguish between these two possibilities would be to detect or rule out the presence of fields coherent on cosmological scales during recombination via their imprint on the cosmic microwave background (CMB) radiation. Dynamically significant magnetic fields present during recombination must be primordial. Indeed, any large-scale primordial magnetic field with a strength comparable to that inferred from the lowest measured intergalactic fields and close to the observational upper limits via Faraday rotation measurements  may well be of cosmological origin. The detected magnetic field strengths in high redshift galaxies  and in damped Lyman alpha clouds  are consistent with the observed Faraday rotation measurements. Primordial magnetic fields can be created in the early universe through phase transitions or via inflation. However, even the invocation of protogalactic dynamos to explain the magnitude of the field involves many uncertain assumptions and still requires a small primordial (pregalactic) seed field . Any seed magnetic field is then amplified adiabatically during gravitational collapse. Hence the possibility of a primordial field merits serious consideration.

It is consequently of primary interest to determine the effects these magnetic fields have on the CMB and, in particular, to consider the limits on any large-scale primordial field from the CMB isotropy measurements. Indeed, the CMB measurements may lead to severe constraints on *homogeneous* magnetic fields in particular classes of cosmological models . Here we demonstrate quite generally that severe limits may be put on *any* large scale magnetic field – be it homogeneous, as considered before, or inhomogeneous – from consideration of the high isotropy of the CMB. That is, we show that any universe in which its observers measure an exactly isotropic CMB cannot have a large-scale magnetic field. We show this by demonstrating the magnetic field evolves in a manner which is irreconcilable with the spacetime geometry required by the (exact) isotropy of the CMB. Contrary to popular belief, this has not been done before.

We wish to take a form of the Copernican principle – which is the fundamental assumption in most of cosmology – combined with our well established knowledge of the high isotropy of the CMB as our starting point. Thus we wish to consider the class of spacetimes which allow *every* observer to see an isotropic CMB, and determine whether any of these are compatible with a large-scale magnetic field. Depending on additional assumptions that can be made, a surprising amount can be inferred from this starting point; for example, if the matter in the spacetime is purely baryonic cold dark matter (‘dust’) with a contribution from radiation then the resulting spacetime must be Friedmann-Lemaître-Robertson-Walker (FLRW). This profound result is known as the Ehlers-Geren-Sachs theorem (EGS, ) and underpins the standard model of cosmology, not to mention many of the results of CMB physics. However, the Universe is not made up only of baryons but also dark matter, most of which could well be non-baryonic . In addition, in light of the supernovae Ia data , there could be a large contribution from a dynamic scalar field (e.g., quintessence ), and there is no reason to believe that this would be homogeneous. Thus the EGS theorem is not the last word on the validity of the homogeneity of the standard model, and recent work has investigated certain generalizations of the EGS theorem  to inhomogeneous models. Most importantly for this work is the role a magnetic field will play in relation to these theorems.

In its most general form, the ‘isotropic radiation field’ theorem identifies spacetimes in which all observers see an ‘isotropic radiation field’. This may be derived from the Einstein-Boltzmann equations for photons in a curved spacetime . It is then easy to show from the multipole expansions of  that in a spacetime in which all observers on some timelike congruence $`u^a`$ see an exactly isotropic radiation field, then this velocity field has two important properties : First, the expansion $`\theta `$ and acceleration $`\dot{u}_a`$ (the motion of the observers under non-gravitational forces) of these observers are related by

$$\dot{u}_a=\stackrel{~}{}_aQ,\theta =3\dot{Q},$$
(1)
where $`Q`$ is related to the energy density of the radiation field. Here ‘$`\stackrel{~}{}_a`$’ refers to the gradient in the instantaneous rest space of the observers (a generalization of the familiar ‘grad’ of Euclidean 3-dimensional vector calculus), and ‘$`^{}`$’ is a proper time derivative along the fluid flow. Second, the congruence must be shearfree, which means that the congruence cannot have any distortion (note that this condition is *not* satisfied in the spacetimes of ). In many cases of physical interest it can then be proven that these observers have zero rotation  For example, if part of the matter consists of a conserved comoving barotropic perfect fluid other than radiation, or for geodesic motion with any matter source, it follows from (1) that the expansion or the rotation must be zero. (For a conserved barotropic perfect fluid, we have $`\dot{u}_a=\stackrel{~}{}_a\varphi `$, and $`p^{}\theta =\dot{\varphi }`$, where $`\varphi 𝑑p/(\mu (p)+p)`$, and $`p^{}=dp/d\mu `$; so, $`\eta _{abc}\stackrel{~}{}^b\stackrel{~}{}^c(Q\varphi )=2(\frac{1}{3}p^{})\theta \omega _a=0`$. For geodesic motion, $`\eta _{abc}\stackrel{~}{}^b\stackrel{~}{}^cQ=\frac{2}{3}\theta \omega _a=0`$ .), a condition supported in part by isotropy of number-counts . For simplicity of exposition we shall examine only the irrotational case here. The general case of non-zero rotation is complicated, and will be considered elsewhere .

Now, if we have a pure magnetic field $`B^a`$, it will contribute an anisotropic pressure, $`\mathrm{\Pi }_{ab}=B_aB_b`$Angled brackets denote the ‘projected, symmetric, trace-free’ part of a tensor, where the projection is into the observers’ instantaneous rest space, and are also used to denote the projected symmetric and trace-free part of time derivatives of spatial tensors (i.e., the angled brackets are applied after the time derivative) . to the energy-momentum tensor of the spacetime, as well as an energy density $`\mu _B=\frac{1}{2}B^2`$, and an isotropic pressure, $`p_B=\frac{1}{3}\mu _B`$ . The energy of the magnetic field must be conserved, as required by one of the Einstein-Maxwell equations (the induction equation), which is an evolution equation for the magnetic field vector;

$$\dot{B}_a=\frac{2}{3}\theta B_a.$$
(2)
This implies that the anisotropic pressure evolves as

$$\dot{\mathrm{\Pi }}_{ab}=\frac{4}{3}\theta \mathrm{\Pi }_{ab}:$$
(3)
it is this evolution equation – which is a direct consequence of Maxwell’s equations – which is inconsistent with Einstein’s field equations when an isotropic radiation field is present, as we shall now show.

In the usual $`1+3`$ covariant approach, Einstein’s field equations are broken into a set of evolution and constraint equations , which we specialize to our particular case, to determine whether the magnetic field may form part of the source. The evolution equation for the shear, in our case, becomes an algebraic equation between the anisotropic pressures, $`\mathrm{\Pi }_{ab}`$, and the electric part of the Weyl tensor, $`E_{ab}`$, which governs tidal forces and gravitational waves:

$$E_{ab}\frac{1}{2}\mathrm{\Pi }_{ab}=A_{ab},$$
(4)
where $`A_{ab}=\stackrel{~}{}_a\dot{u}_b+\dot{u}_a\dot{u}_b`$ is a contribution from acceleration terms. The electric Weyl tensor must satisfy the evolution equation

$$\dot{E}_{ab}+\theta E_{ab}=\frac{1}{2}\left(\dot{\mathrm{\Pi }}_{ab}+\frac{1}{3}\theta \mathrm{\Pi }_{ab}\right).$$
(5)
(These equations may be found in e.g.,  for a general spacetime.) It may be shown that the acceleration terms $`A_{ab}`$ evolve as

$$\dot{A}_{ab}=\frac{1}{3}\theta A_{ab}.$$
(6)
\[This may be shown using (1), the Ricci identities, and Eqns. (5) and (14) from : i.e., the fact that there is no energy flux from the magnetic field and the matter implies that the expansion is homogeneous when the rotation is zero, which in turn implies that the acceleration vector evolves parallel to the velocity vector.\] Now, combining (4), (5) and (6) yields the following evolution equation for the anisotropic pressure:

$$\dot{\mathrm{\Pi }}_{ab}=\frac{2}{3}\theta \mathrm{\Pi }_{ab}\frac{2}{3}\theta A_{ab}.$$
(7)
This is the general evolution equation for the anisotropic pressure for an irrotational spacetime without energy flux which allows an isotropic radiation field. It is now easy to see that a magnetic field is inconsistent with this evolution rate: a magnetic field in these spacetimes must satisfy (3), which, upon substitution into (7), and using (6) implies that

$$\theta \mathrm{\Pi }_{ab}=0:$$
(8)
i.e., the spacetime must be static (which is not relevant to cosmology), or the magnetic field must vanish (in which case the inhomogeneous models will be those found in , and are not FLRW in general). Note that we have made no assumptions on the other matter present, other than it have a perfect fluid form; it needn’t be barotropic, or homogeneous (both assumed in ), or geodesic; it may also consist of many different components, e.g., a perfect fluid and a scalar field (e.g., quintessence). No additional assumption on the spacetime geometry has been made; this result is consequently model-independent. *Thus we see that, in quite general terms, a magnetic field in the Universe is not compatible with an exactly isotropic CMB.*

However, it should be noted that we chose the magnetic field to be ‘pure’ in the frame of the radiation – i.e., ‘frozen into’ the matter; if it were pure in another frame (i.e., with respect to some other observers), the energy flux (the Poynting vector arising from the electric field due to motion through the magnetic field ) would not be zero and the result would not follow directly from the above argument. Also, in principle, models which contain other anisotropic matter stresses need to be considered separately (e.g., if there were anisotropic pressures from some other source which happened to cancel with those of the magnetic field, then this result would not necessarily apply from the above calculation).

Of course, this result is a theoretical result based on perfect isotropy of the CMB, and this magnetic-EGS theorem is not *directly* applicable to the real Universe, since the CMB temperature is not *exactly* isotropic. The original EGS theorem has been generalized to the ‘almost EGS theorem’ , which states that if all fundamental observers measure the CMB temperature to be almost isotropic during some time interval in an expanding dust§§§Note that this assumption is crucial to the theorem  universe, then the universe is described by an almost FLRW model during this time interval. Therefore, there will be an accompanying approximate result (an ‘almost’ magnetic result) from which the observed isotropy of the CMB will lead to severe constraints on the magnetic field. This may be done by simply following the above proof, but now including terms which are zero here but keeping them small: thus we may derive a specific value for the maximum strength of any magnetic field from our knowledge of the CMB anisotropies.There will, however, be some additional assumptions required in such a calculation : for example, in an exact result such as ours, a quantity which is zero (such as the shear) will have zero derivatives; whereas, in the ‘almost’ case where a quantity is small, it may not necessarily have small derivatives; these terms must be kept small by assumption in order to complete the calculation. Indeed, Barrow et al.  have shown that the constraints from the CMB isotropy measurements may provide very strong limits on the strength of a homogeneous component of a primordial magnetic field, and stronger than those imposed by primordial nucleosythesis constraints. However, the results of  are not generic and are model dependent (see  for detailsIn particular, we note that in general Bianchi VIIh models, considered in , *cannot* have a magnetic field as a source .); furthermore, they only apply to homogeneous magnetic fields. In particular, as the models in  have shear, it is not clear as to whether the derived constraints are really characteristic of the presence of the magnetic field.\**\**\**The authors consider Bianchi VIIh models, which are homogeneous flat spacetimes with shear. Therefore the models will have an anisotropic CMB, regardless of how one chooses the matter content, as a consequence of the isotropic radiation field theorem, because it prohibits spacetimes with shear. In their work, the matter was chosen to be a magnetic field and hence was related to the shear in a particular fashion. It is not clear to what extent the anisotropy of the CMB in these models constrain the magnetic field through its particular coupling to the shear, and to what extent it restricts magnetic fields in general. Hence our theoretical result (which is *not* model dependent) complements their results very nicely and provides the necessary theoretical background to such limits. Indeed, a simple calculation shows that the constraints on magnetic fields from this analysis are of a similar order of magnitude obtained in , but an accurate limit requires a separate detailed analysis . We have thus demonstrated that *severe constraints may be placed on magnetic fields in the Universe as a consequence of the extremely high isotropy of the CMB.* In particular, our results also apply to inhomogeneous magnetic fields, and we note that any inflationary scenario leading to significant large-scale primordial magnetic fields would presumably result in magnetic inhomogeneities.

So far we have considered how constraints on *large-scale* magnetic fields may be derived from the CMB anisotropies. We have not, however, given consideration to other types of magnetic field which may be present in the Universe. Specifically, ‘tangled’ magnetic fields – i.e., strong magnetic fields which are inhomogeneous on small scales, but ‘average’ to zero on large scales – may also be constrained in a similar way. For example, consider a bundle of light rays coming from the CMB surface to us of some angular size. We have seen that a magnetic field will distort an isotropic radiation field, so that if our bundle were to pass through a strong small-scale magnetic field, it would produce temperature anisotropies in the CMB multipoles corresponding to size of the tangled field. This may be used to place limits on the strength and size of these inhomogeneous small-scale fields . With the coming launch of the next generation MAP and Planck satellites to measure the high-order multipoles of the CMB spectrum, we may be able to place stringent limits on these types of fields.

We have considered how constraints on magnetic fields may be derived from the temperature anisotropies of the CMB. What, then, would be the implication of discovering, by some other means, a large scale magnetic field in the Universe? It would follow immediately that the real Universe does not satisfy some or all of the assumptions used in this theorem. Most interestingly, it could be indicating that all observers in the Universe do not measure such high isotropy of the CMB; that is, that the Copernican principle does not hold. Thus we see an unexpected route for a physical test of the Copernican principle: looking for magnetic fields.

All of the discussion here has been confined to general relativity. It should be noted that the situation might be different in alternative theories of gravity. For example, in scalar-tensor theories of gravity and low-energy effective theories derived from string theory, the high isotropy of the CMB when combined with the Copernican principle implies that the Universe is isotropic and homogeneous in the case of geodesic matter . However, in string theory the electromagnetic field is coupled to the dilaton (unlike in general relativity in which the electromagnetic field is governed by the conformally invariant Einstein-Maxwell equations), and so the vacuum fluctuations of the electromagnetic field can be significantly amplified by accelerated growth of the dilaton in the pre-big-bang phase; in this case the primordial magnetic field might be strong enough to seed galactic dynamo effects and explain the origin of cosmic magnetic fields observed on galactic and intergalactic scales . Therefore the CMB constraints will likely be even more restrictive in alternative theories of gravity.

We would like to thank Richard Barrett, Roy Maartens and Christos Tsagas for helpful discussions and comments. This work was supported, in part, by the Natural Sciences and Engineering Research Council of Canada and The Royal Society.


text:
# A New Solution to Ginsparg-Wilson Relation from Generalized Staggered Fermion

## I Introduction

A long-lasting folklore that massless fermions are incompatible with lattice formalism of field theories has been almost brought to an end due to the breakthroughs in 90’s characterized in physical language by the introduction of infinitely many flavors . Mathematically, the problem of chiral fermions on lattices has been addressed to Ginsparg-Wilson relation(GWR) ; and overlap solution devised by Neuberger et al provides a promising candidate for the continuing numerical work . Intuitively speaking, overlap construction can be characterized by an operator function including a “compactification” of the argument operator. Dirac-Wilson operator with a heavy negative mass is chosen by Neuberger for the argument operator. In fact, we will present a generic algebraic analysis on the construction of this operator function in this letter, such that a large class of formal solutions to GWR is available.

On the other hand, staggered fermion formalism has dominated the simulations of LQCD for a long time when chiral properties are concerned ; we found an intriguing relation between staggered Dirac operator and noncommutative geometry(NCG) recently . However, in our understanding, the traditional definition of chirality for staggered fermions which reserves a $`U(1)`$ chiral symmetry is unreasonable from a very simple physical argument. The main contribution of this letter is that a new solution to GWR can be found out, if a generalized staggered Dirac operator(GSDO, we will refer conventional staggered Dirac operator as SDO) with a heavy negative mass is taken as argument operator in the operator function mentioned above and a more reasoning alternative definition of chirality is adopted.

This paper is organized as following. Anti-hermitian GSDO is defined in section II, then a brief review of the relation of this operator and NCG is given. Alternative definition of chirality can be fitted into GWR, unless the anti-hermitian GSDO is modified to be “$`\gamma ^5`$-hermitian”. In seciton III, an abstract algebraic analysis is forwarded to characterized the overlap solution to GWR, and GSDO can be matched into the solution class. Some open discussions are put in section IV.

## II Generalized Staggered Fermions, Noncommutative Geometry and “Wrong” Definition of Chirality

Define a GSDO on a 2m-dimensional lattice $`𝒵^{2m}`$ as

$$𝒟_{GS}=\underset{\mu =1}{\overset{2m}{}}(\gamma ^\mu \mathrm{𝟏}_\mu +\frac{ia}{2}\gamma ^{2m+1}\gamma ^\mu \mathrm{\Delta }_\mu )$$
(1)
in which we introduce notations as following: let $`a`$ be lattice constant, $`T_\mu ^\pm `$ be shift operators in conventional meaning; $`_\mu =\frac{1}{2a}(T_\mu ^+T_\mu ^{})`$, $`\mathrm{\Delta }_\mu =\frac{1}{a^2}(2T_\mu ^+T_\mu ^{})`$; $`\gamma ^\mu `$ consist of a set of hermitian gamma matrices in 2m-dimensional Euclidean space and $`\gamma ^{2m+1}=()^m\gamma ^1\gamma ^2\mathrm{}\gamma ^{2m}`$. We hope to point out that $`\{\gamma ^\mu \mathrm{𝟏}`$, $`\gamma ^{2m+1}\gamma ^\mu :\mu =1,2,\mathrm{},2m\}`$ generate 4m-dimensional Clifford algebra $`Cl(4m)`$; in another word, let $`\mathrm{\Gamma }^\mu =\gamma ^{2m+1}\gamma ^\mu `$, $`\mathrm{\Gamma }^{2m+\mu }=\gamma ^\mu \mathrm{𝟏}`$, and $`\mathrm{\Gamma }_\pm ^\mu =\frac{1}{2}(\mathrm{\Gamma }^\mu \pm i\mathrm{\Gamma }^{2m+\mu }),\mu =1,2,\mathrm{},2m`$, then GSDO can be put in a complex coordinate form

$$i𝒟_{GS}=\underset{\mu =1}{\overset{2m}{}}(\mathrm{\Gamma }_+^\mu _\mu ^++\mathrm{\Gamma }_{}^\mu _\mu ^{})$$
(2)
where $`_\mu ^\pm =\frac{1}{a}(T_\mu ^\pm 1)`$. It is easy to verify that $`𝒟_{GS}`$ satisfies the “square-root” condition $`𝒟_{GS}^2=\mathrm{\Delta }`$ where lattice laplacian $`\mathrm{\Delta }=_\mu ^+_\mu ^{}=_\mu \mathrm{\Delta }_\mu `$, and that $`𝒟_{GS}^{}=𝒟_{GS}`$.

It has been shown in  that definition (1) is equivalent to SDO adopting staggered phases $`\eta _\mu (x)`$ as gamma matrices when 2m=4 and Clifford algebra $`Cl(8)`$ stands in a specific representation. We have given a more explicit proof to this equivalency in cases 2m=2,4 in . Here two points shall be clarified. The first one is that when 2m=2,4, the difference between staggered phase formalism and GSDO is that no peculiar representation is specified in definition (1) or (2), so that generically GSDO may have no “stagger” interpretation as a geometric picture. The second one is that the equivalence of SDO and GSDO in higher even-dimensional cases has not yet been proved.

A deep relation between GSDO and NCG has been discovered in our recent work and we give a outline of this relation in this paragraph. Slightly modify Eq.(2) to make the following definition of a “geometric Dirac operator”

$`F:={\displaystyle \underset{\mu =1}{\overset{2m}{}}}(\mathrm{\Gamma }_+^\mu T_\mu ^++\mathrm{\Gamma }_{}^\mu T_\mu ^{})`$$`F`$ possess follow properties:
i) Connes’ distance formula $`d_F(x,y)=sup\{|f(x)f(y)|:[F,f]1\}`$, $`x,y𝒵^{2m}`$ endows an induced metric on lattices. We proved that when 2m=2, $`d_F(x,y)=_{i=1}^2(x^iy^i)^2,x,y𝒵^2`$ in , which recovers Euclidean geometry on 2D lattices. Note that since the difference part of $`F`$ and $`𝒟_{GS}`$ commutes with any $`f`$, $`d_F(,)=d_{𝒟_{GS}}(,)`$.

ii) $`i[F,f]`$ provides an involutive representation of differential form $`df`$ on lattices under Connes’ NCG construction. This representation can be extended to higher forms and due to that $`F`$ is an index-zero Fredholm operator subjected to the so-called “geometric square-root” condition $`F^21`$, this extension is “Junk-free” .

iii) Gauge coupling is added by twisting $`F`$ with link variables $`\omega `$

$`F(\omega ):={\displaystyle \underset{\mu =1}{\overset{2m}{}}}(\mathrm{\Gamma }_+^\mu \omega _\mu ^{}T_\mu ^++\mathrm{\Gamma }_{}^\mu T_\mu ^{}\omega _\mu )`$Then Wilson action for lattice gauge fields can be expressed as

$$S[\omega ]=Tr((F(\omega )F(\omega ))(F(\omega )F(\omega )))$$
(3)
with a properly-defined wedge product, providing that $`\omega `$ is unitary .

Let $`ϵ=\gamma ^{2m+1}\gamma ^{2m+1}`$, then it is obvious that $`\{𝒟_{GS},ϵ\}=0`$. Till now chirality for SDO is defined to be $`ϵ`$; accordingly, staggered formalism has a $`U(1)_A`$ chiral symmetry which is the main advance of this theory. However, we understand that $`\gamma :=\gamma ^{2m+1}\mathrm{𝟏}`$ is a more meaningful choice for chirality in physics. In fact, consider continuum limit(C.L.) in Eq.(1), $`𝒟_{GS}\stackrel{C.L.}{}\gamma ^\mu \mathrm{𝟏}\frac{}{x^\mu }`$, and $`\gamma `$ is nothing but the chirality in continuum theory. Of course on a lattice, $`\{𝒟_{GS},\gamma \}0`$. But it has been a general belief that GWR provides a final answer to the problem of chiral fermions on lattices. Hence, what we plan to do next is to define a lattice Dirac operator $`\stackrel{~}{𝒟}_{GS}`$ such that GWR holds for $`\stackrel{~}{𝒟}_{GS}`$ and $`\gamma `$. As the first step, the anti-hermitian $`𝒟_{GS}`$ must be modified to be a $`\gamma `$-hermitian $`D_{GS}`$ defined as

$$D_{GS}=\underset{\mu =1}{\overset{2m}{}}(\gamma ^\mu \mathrm{𝟏}_\mu +\frac{a}{2}\gamma ^{2m+1}\gamma ^\mu \mathrm{\Delta }_\mu )$$
(4)
such that $`\gamma D_{GS}\gamma =D_{GS}^{}`$.

## III Solution to Ginsparg-Wilson Relation from generalized staggered Dirac Operator

The most popular solution to GWR now is overlap operator devised by Neuberger. We will describe Neuberber’s construction in a pure algebraic manner so that generalizations are easy to be searched out. Note that as in most physical literature, mathematical rigidity on convergency is not cared in detail.

GWR can be abstracted as a operator pair $`(D,\mathrm{\Gamma })`$ and a constraint relation

$$\{D,\mathrm{\Gamma }\}=aD\mathrm{\Gamma }D$$
(5)
where $`\mathrm{\Gamma }`$ satisfies $`\mathrm{\Gamma }^2=1,\mathrm{\Gamma }^{}=\mathrm{\Gamma }`$ and $`D`$ is $`\mathrm{\Gamma }`$-hermitian $`\mathrm{\Gamma }D\mathrm{\Gamma }=D^{}`$. Any solution can be written in the form $`aD=1+V`$; thus it is easy to verify that Eq.(5), together with the properties of $`D,\mathrm{\Gamma }`$, imply $`V`$ is $`\mathrm{\Gamma }`$-hermitian and unitary. On the other hand, any $`\mathrm{\Gamma }`$-hermitian unitary operator $`V`$ will provide a formal solution to Eq.(5).

Define a compactification transformation $`𝒯`$ of operators to be $`𝒯(A)=A(A^{}A)^{\frac{1}{2}}`$ for a generic operator $`A`$. We claim these three statement hold:

$`𝒯(A)`$ is unitary; (6)
$`𝒯(A)=(AA^{})^{\frac{1}{2}}A;`$ (7)
$`\mathrm{\Gamma }A\mathrm{\Gamma }=A^{}\mathrm{\Gamma }𝒯(A)\mathrm{\Gamma }=𝒯(A)^{}`$ (8)
In fact, from knowledge of linear algebra, $`A`$ can be decomposed as $`A=U\lambda V^{}`$ where $`U,V`$ are two unitary operators and $`\lambda `$ is a diagonal operator, so $`A^{}=V\overline{\lambda }U^{}`$. Based on these facts, any $`\alpha `$-order power of $`A^{}A`$ and $`AA^{}`$ are defined to be $`(A^{}A)^\alpha =V|\lambda |^{2\alpha }V^{}`$, $`(AA^{})^\alpha =U|\lambda |^{2\lambda }U^{}`$. Hence, both $`(A^{}A)^\alpha `$ and $`(AA^{})^\alpha `$ are hermitian. Then $`𝒯(A)𝒯(A)^{}=A(A^{}A)^{\frac{1}{2}}(A^{}A)^{\frac{1}{2}}A^{}=A(A^{}A)^1A^{}=1`$, and $`𝒯(A)^{}𝒯(A)=(A^{}A)^{\frac{1}{2}}A^{}A(A^{}A)^{\frac{1}{2}}=1`$, from which Eq.(6) follows. $`(AA^{})^{\frac{1}{2}}A=U\lambda |\lambda |^1V^{}=A(A^{}A)^{\frac{1}{2}}`$, i.e. Eq.(7). Eq.(8) can be shown by $`(\mathrm{\Gamma }A\mathrm{\Gamma }=A^{})(\mathrm{\Gamma }A^{}A\mathrm{\Gamma }=AA^{})([U^{}\mathrm{\Gamma }V,|\lambda |^2]=0)([U^{}\mathrm{\Gamma }V,|\lambda |^{2\alpha }]=0)(\mathrm{\Gamma }(A^{}A)^\alpha \mathrm{\Gamma }=(AA^{})^\alpha )\stackrel{(\text{7})}{}\mathrm{\Gamma }𝒯(A)\mathrm{\Gamma }=𝒯(A)^{}`$. Consequently, an operator function $`\text{D}(A):=(1+𝒯(A))/a`$ provides a formal solution to Eq.(5), if $`A`$ is $`\mathrm{\Gamma }`$-hermitian.

$`A`$ is chosen by Neuberger to be Dirac-Wilson operator $`𝒟_W`$ with a negative mass $`1/a`$, and $`\mathrm{\Gamma }=\gamma ^5`$. We can continue to abstract and to generalize this construction to a finite summation decomposition $`A=\frac{1}{a}+A_0+aA_1+a^2A_2+\mathrm{}+a^NA_N`$ for some natural number $`N`$ in which $`A_i,i=1,2,\mathrm{},N`$ are $`𝒪(1)`$ operators under C.L. and $`A_{even}^{}=A_{even},\mathrm{\Gamma }A_{even}\mathrm{\Gamma }=A_{even}`$, $`A_{odd}^{}=A_{odd},\mathrm{\Gamma }A_{odd}\mathrm{\Gamma }=A_{odd}`$. Hence the $`\mathrm{\Gamma }`$-hermiticity of $`A`$ follows, and $`\text{D}(A)\stackrel{C.L.}{}A_0`$. As a physical requirement, $`A_0`$ has to approach $`i\gamma ^\mu p_\mu `$ in momentum representation under C.L.. More subtleties on the choices of other $`A_i`$ have to be taken care to make $`\text{D}(A)`$ meaningful in physics. In the case of $`A=𝒟_W`$, $`A_0=\gamma ^\mu _\mu `$, $`A_1=\mathrm{\Delta }/2`$.

Now let $`A=D_{GS}\frac{1}{a}`$, and above abstract algebraic construction can be specified as $`\mathrm{\Gamma }=\gamma ^{2m+1}\mathrm{𝟏}=\gamma `$, $`N=1`$, $`A_0=\gamma ^\mu \mathrm{𝟏}_\mu `$, $`A_1=\frac{1}{2}\gamma ^{2m+1}\gamma ^\mu \mathrm{\Delta }_\mu `$. Then a solution to GWR for a 2m-dimensional lattice is

$`\stackrel{~}{𝒟}_{GS}=\text{D}(D_{GS})={\displaystyle \frac{1}{a}}(1+𝒯(D_{GS}))`$It is easy to check that $`\stackrel{~}{𝒟}_{GS}`$ has continuum limit $`i\gamma ^\mu \mathrm{𝟏}\frac{}{x^\mu }`$.

## IV Discussions

i) The old problem concerning staggered formalism of the $`2^m`$ flavor doubling is still presented in our formulation as being shown in the C.L. of $`\stackrel{~}{𝒟}_{GS}`$.

ii) An physically-admissible massless Dirac operator $`𝒟`$ on a 4-dimensional lattice shall fulfill the following conditions:

1. (Locality) in moment space, $`\stackrel{~}{𝒟}(p)`$ is an analytic function of momentum $`p`$;
2. (No doubling) $`\stackrel{~}{𝒟}(p)`$ is invertible for any $`p`$, except $`p=0`$;
3. (Continuum limit) $`\stackrel{~}{𝒟}(p)i\gamma ^\mu p_\mu +𝒪(ap^2)`$, when $`p0`$;
4. (Chirality) GWR holds for $`𝒟`$

$`\{𝒟,\gamma ^5\}=a𝒟\gamma ^5𝒟`$

Only condition 3,4 are considered in this letter, and it is in proceeding whether condition 1,2 are satisfied by $`\stackrel{~}{𝒟}_{GS}`$.

iii) It can be referred as a “Wick rotation” for the modification from $`𝒟_{GS}`$ to $`D_{GS}`$ and it is a profound change since the relations of the former to NCG that we outlined in Sec.II can not be expected to be hold for the latter. However, LQCD action can be put in this form

$`S[\omega ,\psi ,\overline{\psi }]=Tr((F(\omega )F(\omega ))(F(\omega )F(\omega )))+(\psi ,\stackrel{~}{𝒟}_{GS}\psi )`$compared with Eq.(3).

Acknowledgements
This work was supported by Climb-Up (Pan Deng) Project of Department of Science and Technology in China, Chinese National Science Foundation and Doctoral Programme Foundation of Institution of Higher Education in China. J.D. expects to express his gratitude to Professor Chuan Liu for the helpful and inspiring discussion on GWR.


text:
# The Influence of Hydrodynamical Winds on Hot Accretion Disk Solutions

## 1 INTRODUCTION

Hot accretion disks are a natural site where externally or internally produced soft photons can be Comptonized to higher energies by disk electrons. It is this property that makes hot accretion disk models attractive in explaining the spectra of compact X-ray binaries and active galactic nuclei. Alternative models to this picture have also been proposed involving shocks (Chakrabarti & Titarchuk 1995), patchy or uniform coronae (Haardt & Maraschi 1993; Liang & Price 1977), or spatially extended regions in a transition disk model (Misra, Chitinis & Melia 1998). Detailed comparisons of the observed X-ray spectra with predictions of these different models can be used to determine the applicability of each model to a given system. In addition, constraints imposed by X-ray timing analysis can also be used to further distinguish between these models (Nowak et al 1999). Apart from these observational tests, theoretical studies of the stability and self-consistencies of these models are important for determining the circumstances under which they are physically realizable.

In a seminal study, Shapiro, Lightman & Eardley (1976) proposed a hot accretion disk model to explain the X-ray spectra of the black hole system Cygnus X-1. Here, a hot inner disk was assumed to be cooled by external soft photons originating from an outer cold disk. In this model, the gravitational energy dissipated in the hot disk heats the protons, and due to the inefficiency of the Coulomb interaction rate in these environments, leads to a large temperature differences between the protons and electrons. The outer cold disk was assumed to be transformed into a hot inner disk for radii less than a critical value due to a thermal instability associated with the dominance of radiation pressure (Lightman & Eardley 1974). Despite the success of this model in explaining the X-ray observations of Cygnus X-1, there remains several outstanding theoretical issues. Specifically, it was not clear whether an unstable cold disk would in reality evolve to a hot disk (Chen & Taam 1994) or whether the narrow transition between the hot and cold disk was energetically possible. Moreover, the hot accretion disk solution, itself, was thermally unstable.

Against this backdrop, the importance of radial advection of energy in such disks was noted by Ichimaru (1977). Taking this effect into account, Narayan & Yi (1994) constructed a second type of hot disk solution with a self-similar structure called Advection Dominated Accretion Flows (ADAF). These solutions were, in general, radiatively inefficient (with most of the energy advected into a black hole) compared to the non-advective, radiatively cooled solutions found earlier. In contrast to the radiatively cooled solutions, ADAFs are thermally stable, and it was shown that a sharp transition between the hot inner regions and the cool outer regions could exist, provided an additional mechanism (e.g. conduction/turbulent diffusion) is assumed which transfers energy from the hot disk to the transition radius (Honma 1996). The topology of these solutions was studied by Chen et al. (1995) and Zdziarski (1998), and it was found that the advection dominated and radiatively cooled solutions only existed below a critical mass accretion rate. Observational support for the radiatively inefficient ADAF solutions was provided by the X-ray spectral fitting of low luminosity galactic centers and quiescent X-ray Novae (Narayan et al. 1998).

It was noted by Narayan & Yi (1994) that steady state rotating solutions do not exist for an adiabatic index of the accreting flow, $`\mathrm{\Gamma }_D=5/3`$. In this case the ADAF solutions become self-similar spherical Bondi accretion flows. This is not a general result, however, since Ogilvie (1999) has shown that time dependent self similar solutions for $`\mathrm{\Gamma }_D=5/3`$ do exist. If we consider only time independent flows, Narayan & Yi (1994) constructed self-similar rotating solutions in the case of a disk with a nearly equipartition magnetic field (corresponding to an effective adiabatic index of $`4/3`$). Since it is not clear that all accretion flows will contain such fields (e.g., Balbus & Hawley 1998), the applicability of these solutions could be limited. Moreover, Bisnovatyi-Kogan & Lovelace (2000) have recently argued that, in the presence of an equipartition field, the electrons in the disk will be preferentially heated compared to the protons and a two temperature disk may not exist. In contrast, Blandford & Begelman (1999) have recently shown that a rotating self-similar flow with $`\mathrm{\Gamma }_D=5/3`$ is possible, in the presence of a wind, if the accretion flow rate in the disk varies with radius as $`\dot{M}R^P`$. A strong hydrodynamic wind from hot disks is possible, in principle, since the ion temperature in the disk is close to the virial value as shown in the early work by Takahara, Rosner, & Kusunose (1989) and Kusunose (1991). Recently confirmation of these results has been obtained by Beckert (2000) for a variety of viscosity prescriptions. In a parallel development, Chakrabarti (1999) and Das (1999) have studied the possibility of outflows in the context of shock/centrifugal barrier models for the accreting flow. Comparisons with observations reveal that the X-ray spectra of such a wind driven self-similar flow called ADIOS (Advection dominated inflow-outflow solutions) can explain the observed spectra of black hole candidates in quiescence (Quataert & Narayan 1999). In addition, Lin, Misra, & Taam (2000) find that the spectral characteristics of high luminosity black hole systems suggest that winds may be important for these systems as well.

Theoretically, the conditions under which wind dominated accreting flows exist are unclear. For example, Abramowicz, Lasota, & Igumenshchev (2000) argue, based on the sign of the Bernoulli function and on detailed two dimensional hydrodynamical simulations that winds are not expected for disks characterized by low viscosities ($`\alpha <0.1`$). Although modest outflows have been seen in simulations of disks characterized by high viscosities (Igumenshchev & Abramowicz 1999; see also Stone et al. 1999), the relaxation of the assumption of identical adiabatic indices in the disk and wind (reflecting different mechanisms for driving accretion and outflow) can enhance the effectiveness of winds (see Chakrabarti 1999). The existence of winds can also be affected by the presence of convection which can modify the temperature profile (Narayan, Igumenshchev, & Abramowicz 2000). In particular, they find that the thermal disk structure is dependent on the magnitude of the $`\alpha `$ viscosity parameter and the direction in which angular momentum is transported by the convection. For outward transport of angular momentum the ADAF solutions are recovered, whereas for inward transport both non-accreting solutions with lowered temperatures ($`\alpha <0.05`$) and the standard ADAF solutions are found ($`\alpha >0.05`$). Thus the effect of convection on the existence of winds remains inconclusive especially since it is dependent on the unknown distribution of energy dissipation in the vertical extent. It is clear that a strong wind could exist under certain conditions, and it would be useful to study the effects of such a wind on the structure of hot accretion disks.

In ADIOS solutions it is generally assumed that the mass flow rate has a power-law dependence on radius with the power law index, $`P`$, treated as a parameter (Blandford & Begelman 1999; Beckert 2000). In principle, the variation of the accretion rate with radii should be obtained self-consistently from the surface wind velocity and the specific energy of the wind matter. This, in turn, depends on the disk structure, geometry and the thermal properties of the wind. A realistic computation of this disk/wind configuration is difficult due to the complexity of the problem and the unknown vertical distribution of the viscous dissipation in the disk. However, we show that significant progress can be attained if only the dependence of the surface velocity on the disk structure is estimated and its numerical value is parameterized to take into account the uncertainties involved.

In the next section, the basic equations of hot disks with winds are described. In §3 the results of the calculations are presented and in §4 the stability of the solutions is discussed. Finally, we summarize the results and conclude in the last section.

## 2 DISK STRUCTURE EQUATIONS

The standard equations for a Newtonian, keplerian hot accretion disk are modified to take account of mass, energy, and angular momentum loss due to a hydrodynamical wind. To facilitate the presentation of these equations in a compact form, the following dimensionless quantities are introduced: the ion temperature $`\theta kT_i/m_pc^2`$; the electron temperature $`\theta _ekT_e/m_ec^2`$; the radius $`r=R/r_g`$ where $`r_gGM/c^2`$ and $`M`$ is the mass of the black hole; the accretion rate $`\dot{m}\dot{M}/\dot{M}_{Edd}`$ where the Eddington accretion rate is defined to be $`\dot{M}_{Edd}=4\pi m_pcr_g/\sigma _T`$; Thomson optical depth $`\tau n\sigma _TH`$, where $`n`$ and $`H`$ are the number density and half-height of the disk respectively; the speed of the wind at the disk surface $`v_wV_w/c_s`$ where $`c_s=c\theta ^{1/2}`$ is the local isothermal sound speed of the disk and the Coulomb energy exchange frequency $`\nu _c^{}=(r_g/c)\nu _c`$.

For a non-magnetized two-temperature disk i.e. when the ion temperature ($`T_i`$) is much higher than the electron temperature ($`T_e`$), the local pressure is simply $`P_g=nkT_i`$ and the equation for hydrostatic equilibrium can be written as

$`(H/R)^2`$ $`=`$ $`\theta r`$ (1)
In the presence of a wind the accretion rate is no longer constant and hence mass conservation gives,

$`{\displaystyle \frac{d\dot{M}}{dR}}`$ $`=`$ $`4\pi Rm_pnV_w`$ (2)
which in dimensionless form can be re-written using equ. (1) as,

$`{\displaystyle \frac{d\dot{m}}{dr}}`$ $`=`$ $`\tau v_wr^{1/2}.`$ (3)
The torque equation (e.g. Beckert 2000) is,

$`{\displaystyle \frac{dT}{dR}}`$ $`=`$ $`\dot{M}(R){\displaystyle \frac{dl}{dR}}`$ (4)
where $`T`$ is the viscous torque and where it has been assumed that the specific angular momentum of the wind is same as that in the disk i.e. $`l=\sqrt{GMR}`$. The viscous integrated stress is described in terms of the $`\alpha `$-prescription of Shakura & Sunyaev (1973), and it scales as the pressure, i.e. $`T=\alpha P_gH`$ for a disk in hydrostatic equilibrium. The above equation can then be written in dimensionless form as

$`\alpha {\displaystyle \frac{d(\tau \theta r^2)}{dr}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{m}(r)r^{1/2}.`$ (5)
Note that since the accretion rate has a radial dependence (due to the presence of a wind) the above differential equation cannot be trivially integrated.

For a disk in thermal equilibrium, the local vertically integrated energy dissipation rate ($`Q_D`$) should be equal to the sum of the Coulomb energy loss ($`Q_C`$), the advective energy flux ($`Q_A`$) which is non-local radial energy transfer rate associated with the flow and the energy flux carried away by the wind ($`Q_W`$). The energy dissipation rate due to the viscous torque for each half of the disk is

$`Q_D`$ $`=`$ $`{\displaystyle \frac{3}{8\pi }}{\displaystyle \frac{\dot{M}_{Edd}c^2}{r_g^2}}\alpha \theta \tau r^{3/2}.`$ (6)
A fraction of this dissipated energy will be transfered by the protons to electrons via Coulomb interactions. The Coulomb energy exchange flux is

$`Q_C`$ $`=`$ $`{\displaystyle \frac{3}{8\pi }}{\displaystyle \frac{\dot{M}_{Edd}c^2}{r_g^2}}\tau \nu _c^{}(\theta {\displaystyle \frac{m_e}{m_p}}\theta _e)`$ (7)
$``$ $`{\displaystyle \frac{3}{8\pi }}{\displaystyle \frac{\dot{M}_{Edd}c^2}{r_g^2}}\tau \nu _c^{}\theta .`$ (8)
Here, the dimensionless Coulomb frequency is given by,

$`\nu _c^{}`$ $`=`$ $`{\displaystyle \frac{r_g}{c}}\nu _c=4.4\times 10^4\text{ln}\mathrm{\Delta }\tau \theta _e^{3/2}\theta ^{1/2}r^{3/2}`$ (9)
where ln$`\mathrm{\Delta }15`$ is the Coulomb logarithm (Spitzer 1962).

The advective energy flux is

$`Q_A`$ $`=`$ $`{\displaystyle \frac{\dot{M}}{4\pi R^2}}{\displaystyle \frac{P_g}{nm_p}}\xi _A={\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{\dot{M}_{Edd}c^2}{r_g^2}}\dot{m}\theta \xi _Ar^2`$ (10)
where $`\xi _A`$ the advective factor is defined to be (Chen et. al. 1995)

$`\xi _A`$ $`=`$ $`({\displaystyle \frac{3r}{2\theta }}{\displaystyle \frac{d\theta }{dr}}{\displaystyle \frac{r}{\tau }}{\displaystyle \frac{d\tau }{dr}}+{\displaystyle \frac{r}{H}}{\displaystyle \frac{dH}{dr}}).`$ (11)
It has been implicitly assumed in this definition that the specific energy of the gas is $`(3/2)kT_i`$ or in other words, the adiabatic index of the accretion flow $`\mathrm{\Gamma }_D=5/3`$. Note that self-similar (ADAF) solutions without winds, $`\theta r^1`$, $`Hr`$ and $`\tau r^{1/2}`$, leads to $`\xi _A=0`$. This is consistent and expected since, for $`\mathrm{\Gamma }_D=5/3`$ and in the absence of a wind, there are no steady self-similar (ADAF) rotating solutions and only spherical Bondi accretion flows are allowed (Narayan and Yi 1994).

The energy flux carried away by the wind is

$`Q_W`$ $`=`$ $`nV_WE_W={\displaystyle \frac{1}{8\pi }}{\displaystyle \frac{\dot{M}_{Edd}c^2}{r_g^2}}\tau v_we_wr^{5/2}`$ (12)
Here $`E_W`$ is the average energy imparted to a proton in the wind from the local disk and $`e_wE_W/(GMm_p/2R)`$ is its dimensionless form. The external energy required for a particle to escape to infinity is $`GMm_p/2R(3/2)(kT_ikT_f)`$ where $`T_i`$ and $`T_f`$ are the proton temperature of the wind at the disk surface and at infinity respectively. If the wind is isothermal or if $`kT_i<<GMm_p/2R`$ then this external energy input per particle is simply $`GM/2R`$ or the binding energy of the particle. Further, if this energy is obtained only from the local region of the disk from where the particle originates, then $`e_w1`$. Thus, $`e_w<1`$ corresponds to the case where only a fraction of the required energy is obtained from the local region of the disk and the rest presumably is due to heating of the wind from a different region of the disk. The $`e_w>1`$ corresponds to the unlikely case where more energy than required is deposited from the local disk to the wind. This additional energy may be radiated away by a radiatively cooled wind and/or appear as kinetic energy of the wind particles at infinity. The uncertainties involved in the energy coupling between the wind and the disk is factored into $`e_w`$ which is treated as a parameter in this work.

The equation describing conservation of energy ($`Q_D=Q_C+Q_A+Q_W`$ ) can be written in a compact form using equs. (6), (8), (10) and (12) as

$`{\displaystyle \frac{\nu _c^{}r^{3/2}}{\alpha }}+{\displaystyle \frac{2}{3}}{\displaystyle \frac{\dot{m}\xi _A}{\alpha \tau r^{1/2}}}+{\displaystyle \frac{1}{3}}{\displaystyle \frac{v_we_w}{\alpha \theta r}}=1.`$ (13)
The value of the electron temperature ($`\theta _e`$) which is required to calculate the Coulomb frequency (equ. 9), depends on the dominant radiative cooling mechanism operating in the hot disk. Following Zdziarski (1998) it is assumed here that the disk cools primarily by unsaturated Comptonization of soft photons. This is motivated by the high energy spectral characteristics of black hole binaries in the low-state since the Comptonization process is characterized by the Compton y-parameter $`(y4kT_e/m_ec^2)`$ max$`(\tau ,\tau ^2)1`$. Thus

$`\theta _e{\displaystyle \frac{y}{4max(\tau ,\tau ^2)}}`$ (14)
can be used to estimate the electron temperature using $`y`$ as a parameter of order unity. Technically the electron temperature profile in the disk should be obtained by solving the geometry dependent radiative transfer equations in both vertical and radial directions. Thus it should be emphasized that equation (14) is a simplifying assumption which may be approximately valid only under certain conditions (Shapiro, Lightman & Eardley 1976).

Determining the wind velocity at the surface of the disk is complicated because it depends on the vertical structure of the disk which in turn depends on the detailed variation of the unknown viscosity with height. Furthermore, the opening angle for the wind flow needs to be determined using the geometry of the accreting flow and the pressure gradient in the wind in the radial direction. However, a qualitative estimate of the dependence of the wind velocity with respect to the average temperature and gravitational potential in the disk, can be obtained by studying the simpler geometry of a spherical flow from a radius $`R`$, temperature $`T_i`$ and potential $`GM/2R`$. The factor of two in the potential takes into account the angular momentum of the disk. The estimate obtained from such a spherical flow can be scaled by an unknown factor $`\eta `$ which would then take into account the uncertainties in the geometry and vertical structure of the disk,

$`v_w=\eta v_s(\theta r)`$ (15)
where $`v_s`$ is the velocity of the wind at the surface of the sphere divided by the isothermal sound speed ($`=\theta ^{1/2}c`$). Thus the motivation here is to obtain a qualitative dependence of the disk surface wind velocity with temperature and radius of the disk, rather than a quantitative estimation of its magnitude. The standard equations of the hydrodynamical theory of stellar winds have been used to compute the wind velocity at the surface of the sphere. Using Bernoulli’s equation and the fact that the Mach number at the sonic point is unity, the following relation for the sonic radius is obtained

$`{\displaystyle \frac{r_c}{r_s}}`$ $`=`$ $`{\displaystyle \frac{\lambda [\frac{53\mathrm{\Gamma }}{4(\mathrm{\Gamma }1)}]}{\frac{v_s^2}{2}+\frac{\mathrm{\Gamma }}{\mathrm{\Gamma }1}\lambda }}`$ (16)
where $`r_c`$ is the location of the sonic point, $`r_s`$ is the radius of the sphere, $`\mathrm{\Gamma }`$ is the adiabatic index of the wind and $`\lambda `$ is the ratio of the gravitational potential at the surface to the temperature of the sphere i.e. $`\lambda \varphi /kT_i=GMm_p/(2r_skT_s)`$. Note that the adiabatic index in the wind may differ from that in the disk (see Chakrabarti 1999). Conservation of mass and the equation of state, $`P_g\rho ^\mathrm{\Gamma }`$ (where $`\rho `$ is the mass density) gives,

$`v_s`$ $`=`$ $`({\displaystyle \frac{\lambda }{2}})^{\frac{1}{2}\frac{\mathrm{\Gamma }+1}{\mathrm{\Gamma }1}}\mathrm{\Gamma }^{\frac{1}{\mathrm{\Gamma }1}}({\displaystyle \frac{r_c}{r_s}})^{\frac{3\mathrm{\Gamma }5}{2(\mathrm{\Gamma }1)}}`$ (17)
In this work, equations (16) and (17) are solved simultaneously to yield $`v_s(\lambda ,\mathrm{\Gamma })`$ where $`\lambda `$ for an isothermal disk is $`1/(2\theta r)`$. If $`v_s<<\lambda `$, then

$`v_s(\lambda ,\mathrm{\Gamma })`$ $`=`$ $`({\displaystyle \frac{1}{2}})^{\mathrm{\Gamma }_2}\mathrm{\Gamma }^{\frac{3}{2}}\lambda ^2[(53\mathrm{\Gamma })(1{\displaystyle \frac{\lambda (\mathrm{\Gamma }1)}{\mathrm{\Gamma }}})]^{\mathrm{\Gamma }_1}`$ (18)
where $`\mathrm{\Gamma }_1=\frac{3\mathrm{\Gamma }5}{2(\mathrm{\Gamma }1)}`$ and $`\mathrm{\Gamma }_2=\frac{7\mathrm{\Gamma }9}{2(\mathrm{\Gamma }1)}`$ . For the isothermal case ($`\mathrm{\Gamma }1`$), equ. (18) reduces to

$`v_s(\lambda ,\mathrm{\Gamma }1)=({\displaystyle \frac{\lambda }{2}})^2\text{exp}(3/2\lambda )`$ (19)
There are several points which should be noted here. First, equs. (18) and (19) seem to suggest that $`v_s`$ has a maxima for variations with $`\lambda `$. This is an artifact of the assumption $`v_s<<\lambda `$ and if equs. (16) and (17) are used instead (as is done in this work), then $`v_s`$ is a monotonic function of $`\lambda `$. Second, there is a maximum value of $`\lambda `$ beyond which there are no wind solutions or $`v_w=0`$ for $`\theta <\theta _{max}=(\mathrm{\Gamma }1)/(2\mathrm{\Gamma }r)`$. Third, inspection of equ. (18) reveals that for θ
>θmax
>𝜃subscript𝜃𝑚𝑎𝑥\theta\mathrel{\mathchoice{\lower 0.5pt\vbox{\halign{$\mathnew\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}\theta_{max}, $`v_w`$ is a rapidly varying function of $`\theta `$, a factor two increase in $`\theta `$ could lead to an increase of two orders of magnitude in $`v_w`$ (Lin, Misra, & Taam 2000).

The equations corresponding to conservation of mass (3), angular momentum (5) and energy (13) can be solved to obtain the disk structure using the estimated values for electron temperature (14) and wind velocity (15).

## 3 DISK SOLUTIONS WITH WINDS

### 3.1 Local Solutions

The differential equations governing the disk structure can be solved using appropriate boundary conditions. However, these equations can be reduced to algebraic ones, by treating the logarithmic derivatives of physical quantities as constant (i.e., independent of radii) parameters. The solutions for this set of algebraic equations depends only on the local radius and hence are local solutions. Such solutions generally provide important insight on the disk structure and its dependence on disk parameters.

To obtain local solutions it is assumed that the logarithmic derivative of temperature ($`\xi _T=d\text{log}T/d\text{log}R`$) and accretion rate ($`P=d\text{log}\dot{M}/d\text{log}R`$) with respect to the logarithm of the radius are constants. Integrating equation (5) with the constant of integration chosen such the viscous torque vanishes in the inner radius, leads to

$$\tau =\frac{\dot{m}(r)r^{\frac{3}{2}}\text{J}(r)}{(1+2P)\alpha \theta }$$
(20)
where J $`(r)1(6/r)^{P+1/2}`$ takes into account the torque free inner boundary condition at $`r=6`$. The advective fraction (equ. 11) is reduced to

$$\xi _A=3(\xi _T1)+P$$
(21)
where it has been assumed that J$`(r)1`$ is a constant. While it is convenient to consider $`\xi _T`$ as a parameter, $`P`$ can be obtained self-consistently by mass conservation (equ. 3),

$$P=\frac{\tau r^{1/2}}{\dot{m}}\eta v_s(\theta r)=\frac{1}{1+2P}\left(\frac{\eta }{\alpha }\right)\left(\frac{v_s(\theta r)}{\theta r}\right).$$
(22)
Hence, the local disk solution is obtained by solving the energy equ. (13) with equs. (20), (21) and (22). Figure 1 shows the variation of the normalized accretion rate, $`\dot{m}`$, with optical depth, $`\tau `$, for $`r=40`$. Comparison with Figure 1 of Zdziarski (1998) shows that the topology of local solutions does not change with the inclusion of a wind. For accretion rates less than a critical value ($`\dot{m}<\dot{m}_{crit}`$), there exists a radiatively cooled solution branch and an advection/wind dominated solution branch. No solutions exist for $`\dot{m}>\dot{m}_{crit}`$. Note, that for $`\mathrm{\Gamma }=5/3`$, comparisons cannot be made with respect to solutions without winds since steady state solutions do not exist.

For the advection/wind dominated solution branch, the Coulomb exchange rate can be neglected ($`Q_C0`$) and equ. (13) can be written for $`\xi _T=1`$ as

$`{\displaystyle \frac{3}{2}}`$ $`=`$ $`P(1+2P)\theta r+{\displaystyle \frac{e_W}{2}}({\displaystyle \frac{\eta }{\alpha }})({\displaystyle \frac{v_s(\theta r)}{\theta r}})`$ (23)
$`=`$ $`({\displaystyle \frac{\eta }{\alpha }})v_s(\theta r)\left(1+{\displaystyle \frac{e_W}{2\theta r}}\right)`$
The above with equs. (20) and (22) lead to self-similar solutions for an advection/wind dominated flow: $`\theta 1/r`$, $`\tau r^{P1/2}`$ and $`\dot{m}r^P`$. Thus, with the inclusion of a wind, self-similar flows are maintained which is a direct consequence of the form of the wind velocity chosen in equ. (15). For example, if the parameter $`\eta `$ is a function of radius rather than a constant as assumed here, self-similar solutions will no longer be allowed. The solution to equ (23) is shown in Figure 2, where $`\theta r`$ is plotted as a function of $`\eta /\alpha `$ for different values of $`\mathrm{\Gamma }`$ and $`E_W`$. For values of $`\eta /\alpha >2`$, $`\theta r`$ depends only logarithmically with $`\eta /\alpha `$ and is $`0.1`$. The ratio of the wind to advective cooling $`Q_W/Q_A=e_W/(2\theta r)5`$ making both cooling channels equally important for this solution branch. It is illustrative to solve for $`P`$ using equs. (22) and (23) yielding

$$P=\left(1+\sqrt{1+24/(e_W+2\theta r)}\right)/4$$
(24)
which implies that $`P0.9`$ for $`e_W1`$ and $`\theta r0.1`$. This is shown in Figure 3 where $`P`$ is plotted as a function of $`\eta /\alpha `$. Note that from equ. (24) a smaller value of $`e_W=0.5`$ corresponds to $`P1.3`$ and a larger value $`e_W=1.5`$ corresponds to $`P0.7`$. Thus P remains of order unity as long as $`e_W1`$ which is expected unless the wind is predominately heated non-locally or if the energy transferred by the disk to the wind greatly exceeds the energy required to generate the wind.

At the critical accretion rate, $`\dot{m}_{crit}`$, radiative cooling becomes comparable to the local dissipation ($`Q_CQ_D`$). Using the self-similar solution obtained above, $`\dot{m}_{crit}r^{3/10}`$ which is the same dependence obtained for hot solutions without winds (Zdziarski 1998). The value of $`\dot{m}_{crit}`$ is however, generally lower than for solutions without winds.

In summary, for a wide range of values for the viscosity parameter ($`\alpha `$) and the wind parameter ($`\eta `$), advection/wind dominated solutions are characterized by $`\theta 0.1/r`$ and $`P0.9`$, independent of the accretion rate or radius. We note that the result that the ion temperatures of wind dominated solutions are about an order of magnitude smaller than the no-wind ADAF solutions justifies the assumption of a keplerian disk.

### 3.2 Global Solutions

It is prudent to confirm the inferences obtained from local analysis by constructing global solutions of the disk. These global solutions involve solving the three differential equations for conservation of mass, angular momentum and energy (equs. 3, 5 and 13) with three independent boundary conditions. For convenience, the boundary conditions adopted in this work are that the viscous torque vanishes at the inner radius ($`r_i=6`$), and the accretion rate, $`\dot{m}(r_o)`$, and the temperature, $`\theta (r_o)`$, are specified at the outer radius, $`r_o`$. For a given $`\dot{m}(r_o)`$, $`\theta (r_o)`$ is computed from the local solution described in §3.1.

Figure 4 illustrates the variation of the mass accretion rate with respect to radius for different values of $`\dot{m}(r_o)`$. Here, $`\theta (r_o)`$ is chosen from the advection/wind dominated branch of the local solutions and hence, the global solution presented is also advection/wind dominated. The global solutions are, indeed, similar to the type inferred from local analysis with $`P0.9`$ except when r
>6
>𝑟6r\mathrel{\mathchoice{\lower 0.5pt\vbox{\halign{$\mathnew\displaystyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\textstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\scriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\scriptscriptstyle\hfil#\hfil$\cr>\crcr\sim\crcr}}}}6 where as expected the flow is no longer self-similar. The dotted line is the critical accretion rate (beyond which hot disk solutions no longer exist) computed using local analysis but with values of $`P`$ and $`\xi _T`$ taken from the global solution.

Figure 5 presents the global solutions for the case when $`\theta (r_o)`$ is chosen from the radiatively cooled branch of the local solutions. For $`\dot{m}(r_o)<<\dot{m}_{crit}(r_o)`$, the accretion rate is a constant since winds are not important (i.e., $`\theta <<1/r`$). However, for m˙(ro)
<m˙crit(ro)
<˙𝑚subscript𝑟𝑜subscript˙𝑚𝑐𝑟𝑖𝑡subscript𝑟𝑜\dot{m}(r_{o})\mathrel{\mathchoice{\lower 0.5pt\vbox{\halign{$\mathnew\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}\dot{m}_{crit}(r_{o}), winds do become important for certain regions of the disk. In such cases the wind regulates the flow such that m˙(r)
<m˙crit(r)
<˙𝑚𝑟subscript˙𝑚𝑐𝑟𝑖𝑡𝑟\dot{m}(r)\mathrel{\mathchoice{\lower 0.5pt\vbox{\halign{$\mathnew\displaystyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\textstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\scriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}{\lower 0.5pt\vbox{\halign{$\mathnew\scriptscriptstyle\hfil#\hfil$\cr<\crcr\sim\crcr}}}}\dot{m}_{crit}(r). These wind regulated solutions occur for a limited range of $`\dot{m}(r_o)`$ which is due to the weak dependence of $`\dot{m}_{crit}`$ with radius ($`\dot{m}_{crit}r^{3/10}`$). In contrast, the accretion rate at the innermost regions of the disk can be significantly reduced from $`\dot{m}(r_o)`$ for the advection dominated branch. For radiative cooling mechanisms other than soft photon Comptonization, $`\dot{m}_{crit}`$ may depend more steeply on radius, leading to a wider range of situations where a wind-regulated accretion flow would exist.

## 4 STABILITY OF DISK SOLUTIONS WITH WINDS

Radiatively cooled hot disks (without winds) are known to be thermally unstable. This result remains valid for the radiatively cooled branch described in §3.1, simply because winds are not important for these solutions. Advection dominated hot disks (without winds) are known to be thermally stable and again this result is not altered for the advection/wind dominated solution since the energy loss rate due to the wind increases with temperature.

Formally the criterion for thermal stability is that the derivative of the energy dissipation rate with respect to temperature should be smaller than the derivative cooling terms with respect to temperature, at a constant surface density. This leads to (Appendix: equ. A1)

$$S_T\frac{Q_C}{2}+2Q_A+\xi _{v\theta }Q_WQ_D>0$$
(25)
as the criterion for thermal stability provided that the advection timescale is longer than the thermal timescale. Here, $`\xi _{v\theta }d\text{log}v_s/d\text{log}\theta >1`$. For the wind/advection dominated solutions, $`Q_D=Q_A+Q_W`$ and $`Q_C<<Q_D`$, hence the disk will be locally stable to thermal perturbations. For the radiatively cooled branch, at high accretion rates, advection and wind cooling may be sufficiently important to stabilize the disk.

On timescales longer than the viscous timescale, the surface density (or optical depth) will respond to fluctuations in the disk and its temporal behavior is governed by,

$$\frac{\tau }{t}=\frac{1}{R}\frac{}{R}(R\tau V_r)\frac{\tau V_W}{H}$$
(26)
where $`V_r`$ is the radial velocity. From the torque equation (see equ. 5), it can be seen that fluctuations in $`V_r`$ are linearly proportional to variations in temperature in the long wavelength limit. On the other hand $`V_W/H(v_w)`$ is a rapidly increasing function of temperature. Thus, in thermal equilibrium, if $`\tau `$ is inversely proportional to temperature, the disk will be unstable to long wavelength perturbations. In other words, if a decrease in $`\tau `$ leads to an increase in the temperature, the mass loss rate by the wind will increase rapidly which will cause a further decrease in $`\tau `$ leading to a runaway situation.

Technically, for the local wind dominated solution described in §3.1 without radiative cooling ($`Q_C=0`$), the temperature of the disk is independent of $`\tau `$ (equ 23). However, even if there is only a small contribution due to radiative cooling (i.e. $`Q_C<<Q_D`$), an increase in $`\tau `$ will cause $`Q_C`$ to increase which leads to a decrease in the temperature. Thus, wind dominated solutions are subject to the instability described above in the long wavelength limit. A linear stability analysis leads to the dispersion relation for $`\lambda >>R`$ (Appendix: equ. A4)

$$i\omega =\frac{5}{2}(\frac{V_r}{R})(\frac{Q_C}{S_T})P(\xi _{v\theta }1).$$
(27)
For wind dominated solutions, $`\xi _{v\theta }>1`$ and the disk is thermally stable ($`S_T>0`$), hence $`i\omega `$ is real and positive indicating that the mode is unstable. For a radiatively cooled disk, $`S_T<0`$ and this mode is stable for such solutions. The growth time of the wind instability ($`1/\omega `$) can be significantly longer than the viscous (or accretion) timescale ($`R/V_r`$) if $`Q_C/S_TQ_C/Q_D<<1`$. This implies that the disk structure can be described in terms of local steady state solutions during the linear growth of the wind driven instability. Equ. (27) and the inferred arguments depend on the form of the radiative cooling mechanism (equ. 14) assumed in the analysis. In the general case, $`Q_C\tau ^\beta `$, and it can be shown that the disk will suffer this instability provided $`\beta >1`$ and if $`S_T>0`$. For the unsaturated Comptonization model (i.e. equ. 14) $`\beta =7/2`$ and hence the disk is found to be unstable. Note however, that $`\beta <1`$ corresponds to the situation where the electron temperature is proportional to the optical depth and varies faster than $`\theta _e\tau ^{2/3}`$. Since for most physical radiative processes, $`\theta _e`$ is expected to be inversely proportional to $`\tau `$, hot disks with winds would in general be unstable. The above results are based on local and linear stability analysis and require confirmation by global simulations.

## 5 SUMMARY AND DISCUSSIONS

The effect of a hydrodynamical wind on a hot accretion disk has been studied for disks characterized by an adiabatic index, $`\mathrm{\Gamma }_D=5/3`$. It is found that the topology of local ADAF solutions found earlier (for $`\mathrm{\Gamma }_D<5/3`$) remains the same for $`\mathrm{\Gamma }_D=5/3`$ with the inclusion of a wind. In particular, there exists two branches of solutions for accretion rates below a critical value, while there are no solutions for higher accretion rates. One solution corresponds to the radiative cooling branch where the effects of the wind are negligible whereas both winds and advection are important for the second branch.

We have shown that the local solutions for the wind/advection dominated branch may be described in a self-similar form with the temperature, $`\theta 1/r`$, optical depth $`\tau r^{P1/2}`$ and accretion rate $`\dot{m}r^P`$. For a wide range of viscosity and wind parameters, the temperature is about 10% of the virial value ($`\theta r0.1`$) and $`P0.9`$, independent of the accretion rate and radius. These results were confirmed by global analysis for which self-similar solutions were obtained for radii much greater than the inner disk radius, thereby providing a theoretical basis for the model advanced by Blandford & Begelman (1999). Global solutions also revealed that under certain boundary conditions (see §3.2), the accretion rate $`\dot{m}\dot{m}_{crit}r^{3/10}`$. For these “critical” solutions radiative cooling and energy loss due to the wind, were both important and hence they are not self-similar.

The wind/advection dominated solutions are thermally stable, but are locally unstable to a new type of instability which we denote as “wind driven instability.” The presence of a wind causes the disk to be unstable to long wavelength perturbations in the surface density. Since the growth rate of this instability is inversely proportional to the ratio of the radiative cooling to the gravitational energy dissipation rates, the timescale for the instability could be much larger than the viscous time scale in the disk for a low radiative cooling efficiency.

The accretion rate index, P, $`1`$ implies that, if the outer radius of the hot disk is much larger than the inner one, most of the accreting matter is lost in the form of a wind rather than accreted onto the star. If the compact object is a neutron star or white dwarf, it is this residual accretion which gives rise to the X-ray luminosity from the surface of the star. Loeb, Narayan & Raymond (2000) point out that in such a scenario (with nearly similar outer boundary conditions), the surface emission from the star will scale roughly as $`R_{}^{P1}`$, where $`R_{}`$ is the radius of the star. From the observed luminosities of soft X-ray transients and cataclysmic variables, they estimate $`P=0.9\pm 0.5`$, which is in agreement with the theoretical results presented here.

For this radial dependence of the mass flow rate in the disk, the total energy dissipated in the disk scales as log $`(R_o/R_i)`$, where $`R_o`$ and $`R_i`$ are the outer and inner radii of the hot disk, rather than $`R_i^1`$. Thus both the inner and outer regions of hot disk are equally important energetically. It is generally assumed that the hot disk is surrounded and connected to an outer cold disk by a narrow transition region. The energy required to make this transition cannot be generated locally if the transition zone is small (Misra & Melia 1996) and is probably due to global transfer from the inner hot regions to the transition zone (Honma 1992, see also Abramowicz et al. 2000). For the wind/advection dominated case this may not be possible since the energy generated in the inner regions is only comparable to that of the outer region. In this case, the transition zone would then have to be extended and could significantly contribute to the luminosity of the source (Misra & Melia 1996). However, the above speculation has to be verified by detailed multi-dimensional simulations of the transition zone taking into account the possibility that the hot disk may be wind dominated.

Since the radiative cooling rate (equ. 8) for the self-similar solution scales as $`R^{(14P15)/4}R^{1/2}`$ for $`P0.9`$, most of the luminosity arises from the outer regions of the hot disk. This result, however, depends on the radiative cooling mechanism operating in the hot disk. For example, if the electron temperature is chosen to be $`\theta _eR^{\alpha _e}`$ instead of equ. 14, then the radiative flux is $`R^{2P33\alpha _e/2}R^{(1+3\alpha _e/2)}`$. Hence if $`\alpha _e>2/3`$ then the luminosity will be dominated by the inner region.

The terminal velocity of the wind should be approximately equal to the sound speed at the surface of the disk. This implies that the bulk energy per proton of the wind flow at infinity is $`0.1GM/R`$ since the temperature of the disk $`\theta 0.1r`$ For $`P1`$ most of the wind arises from the outermost radius and hence the total bulk energy flux in the wind $`0.1\dot{M}_{wind}GM/R_o`$, where $`\dot{M}_{wind}`$ is the total mass loss rate and $`R_o`$ is the outer radius of the hot disk. Since the radiative efficiency of such disks is $`<<0.1`$, the bulk energy of the wind flow in some circumstances may be much larger than the radiative flux from the source.

The viscous timescale for the self-similar solution is

$$t_v=\frac{R}{V_R}\frac{R^{1/2}}{(1+2P)(\theta r)\alpha }(\frac{r_g}{c})$$
(28)
which numerically is $`t_v250(r_g/c)0.01`$ s for $`\alpha =0.1`$, $`R=50r_g`$ and a 10 $`M_{}`$black hole. Thus, a wind driven instability may occur on timescales $`110`$ s for $`Q_C/Q_D0.10.01`$. The instability may not lead to disruption of the disk, but instead give rise to large amplitude oscillations. For example, the surface density could increase, accompanied by a decreasing temperature. When the temperature decreases sufficiently, the wind is quenched and radiative cooling becomes important. Since the radiatively cooled disk will be thermally unstable, the system may return to the wind dominated state on thermal timescales, giving rise to cyclic behavior. Alternatively, the surface density may decrease and the disk may then be evacuated and finally filled again on the viscous timescale. In either case, the disk may evolve through epochs of high wind activities on timescale $``$ seconds. Such a disk/wind interaction may be important for explaining the spectral and temporal behavior of the regular bursting activity exhibited by the microquasar GRS 1915-105 (Taam et al 1997).

Future work involving time dependent numerical simulations of hot disk with winds are needed to confirm the above speculations and to quantify the rate of mass loss and varying spectral characteristics during such events. Such studies may also shed light on the relevance of wind instabilities as a mechanism for either oscillatory behavior or evaporative behavior in the hot inner disk regions. Besides these issues, further theoretical investigations of the transition from a cold to a hot wind dominated disk and the role of two dimensional effects along the lines suggested by Turolla & Dullemond (2000) would also be illuminating.

R.M. gratefully acknowledges support from the Lindheimer Fellowship at Northwestern University.

## Appendix A Linear stability analysis in the long wave-length approximation

The linear response of the disk structure is calculated for perturbation of the ion temperature $`\theta \theta (1+\delta \theta e^{i(\omega tkr)})`$ and optical depth $`\tau \tau (1+\delta \tau e^{i(\omega tkr)})`$, where $`\delta \theta `$ and $`\delta \tau `$ are small dimensionless quantities. In the long wavelength approximation, it is assumed that $`kr<<1`$ such that $`e^{ikr}1`$ which implies that the radial derivatives of perturbed quantities can be neglected.

In this formalism, the linear response of physical quantities (e.g., the energy flux) can be obtained directly from the dependence of the quantity on temperature and optical depth. The dissipative energy rate $`Q_D\theta \tau `$ (equ. 6) and hence $`\delta Q_D=\delta \theta +\delta \tau `$. The Coulomb energy exchange flux $`Q_C\tau ^{7/2}\theta ^{1/2}`$ (equs. 8, 9, 14) while the energy flux carried away by the wind $`Q_W\tau v_w(\theta )`$ (equ. 12 ). The advective energy flux (equ. 10), $`Q_A\dot{m}\theta \tau \theta ^2`$ since $`\dot{m}\tau \theta `$ (equ. 5) for a self-similar steady state solution with $`kr<<1`$. Note that in the long wavelength approximation the advective factor ($`\xi _A`$) is invariant to perturbations. Thus the normalized response of the total energy flux $`Q_TQ_DQ_CQ_WQ_A`$ is,

$`Q_T\delta Q_T`$ $`=`$ $`(Q_D{\displaystyle \frac{1}{2}}Q_C\xi _{v\theta }Q_W2Q_A)\delta \theta +(Q_D{\displaystyle \frac{7}{2}}Q_CQ_WQ_A)\delta \tau `$ (A1)
$`=`$ $`S_T\delta \theta {\displaystyle \frac{5}{2}}Q_C\delta \tau `$
where $`\xi _{v\theta }d\text{log}v_s/d\text{log}\theta `$. On thermal timescales $`\delta \tau <<\delta \theta `$ and $`Q_T\delta Q_T/\delta \theta =S_T>0`$ is the criterion for thermal stability (equ. 25).

The time dependent form of the mass conservation equation (3) governs the temporal behavior of the surface density or optical depth (equ. 26)

$$(\frac{r_g}{c})\frac{\tau }{t}=\frac{1}{r}\frac{\dot{m}}{r}\tau v_wr^{3/2}$$
(A2)
Noting that $`\dot{m}\tau \theta `$, the linearized form of the above equation is

$$i\omega \tau \delta \tau =(\frac{c}{r_g})(\frac{1}{r}\frac{\dot{m}}{r})\delta \theta (1\xi _{v\theta })$$
(A3)
For timescales much longer than the thermal timescale, $`\delta Q_T0`$ i.e. the system will be in thermal equilibrium. Combining equations (A1) and (A3), yields the dispersion relation in the long wavelength limit (equ. 27)

$$i\omega =\frac{5}{2}(\frac{V_r}{R})(\frac{Q_C}{S_T})P(\xi _{v\theta }1)$$
(A4)
where $`V_r`$ is the magnitude of the radial velocity at radius $`R`$.


text:
# Interplane Transport and Superfluid Density in Layered Superconductors

## Abstract

We report on generic trends in the behavior of the interlayer penetration depth $`\lambda _c`$ of several different classes of quasi two-dimensional superconductors including cuprates, Sr2RuO4, transition metal dichalcogenides and organic materials of the $`(BEDTTTF)_2X`$ \- series. Analysis of these trends reveals two distinct patterns in the scaling between the values of $`\lambda _c`$ and the magnitude of the DC conductivity: one realized in the systems with a Fermi liquid (FL) ground state and the other seen in systems with a marked deviation from the FL response. The latter pattern is found primarily in under-doped cuprates and indicates a dramatic enhancement (factor $`10^2`$) of the energy scale $`\mathrm{\Omega }_C`$ associated with the formation of the condensate compared to the data for the FL materials. We discuss implications of these results for the understanding of pairing in high-$`T_c`$ cuprates.

The formation of the superconducting condensate in elemental metals and their alloys is well understood within the theory of Bardeen, Cooper and Schrieffer (BCS) in terms of a pairing instability in the ensemble of Fermi liquid quasiparticles. Applicability of the FL description to high-$`T_c`$ cuprate superconductors is challenged by remarkable anomalies found in both the spin- and charge response of these compounds in the normal state . Because quasiparticles are not well defined at $`T>T_c`$ in most cuprates it is natural to inquire into the distinguishing characteristics of the superconducting condensate which appears to be built from entirely different ”raw material”. Infrared spectroscopy is perfectly suited for the task. Indeed, the analysis of the optical constants in far infrared (IR) unfolds the process of the formation of the condensate $`\delta (0)`$-peak in the dynamical conductivity and also gives insight into the single-particle excitations in a system both above and below $`T_c`$ .

In this paper we focus on the interplane properties of high-$`T_c`$ superconductors. We will show that the distinctions in the behavior of the condensate in the FL superconductors and high-$`T_c`$ cuprates are most radical in the case of the $`c`$-axis interplane response. The analysis of the generic trends seen in the behavior of the $`c`$-axis condensate (correlation between the penetration depth $`\lambda _c`$ and the DC conductivity) allows us to infer the energy scale $`\mathrm{\Omega }_C`$ associated with the development of the superfluid. This energy scale is of the order of the energy gap $`2\mathrm{\Delta }`$ in FL superconductors but may dramatically exceed the gap in systems lacking well defined quasiparticles at $`T>T_c`$ (primarily in under-doped cuprates). We discuss a connection between the magnitude of $`\mathrm{\Omega }_C`$ and the nature of the normal state response.

Experimentally, the inter-layer $`c`$-axis penetration depth of a layered superconductor can be measured using a variety of techniques such as microwave absorption, magnetization, vortex imaging  or IR spectroscopy . In the IR approach, the magnitude of $`\lambda _c`$ can be reliably extracted from the imaginary part of the complex conductivity $`\sigma _1(\omega )+i\sigma _2(\omega _2)`$ as $`1/\lambda ^2=\omega \times \sigma _2(\omega )`$. Regardless of the method employed, the penetration depth in several families of cuprates reveals a universal scaling behavior with the magnitude of $`\sigma _{DC}(T=T_c)`$ (Fig. 1): the absolute value of $`\lambda _c`$ is systematically suppressed with the increase of the normal state conductivity. The scaling is obeyed primarily in under-doped cuprates. The deviations from the scaling are also systematic and are most prominent in over-doped phases (red symbols in Fig. 1). Such deviations are a direct consequence of a well-established fact: on the over-doped side of the phase diagram the DC conductivity increases whereas $`\lambda _c`$ is either unchanged or may show a minor increase .

We find a similar scaling pattern in other classes of layered superconductors, including organic materials, transition metal dichalcogenides and Sr2RuO4 (Fig. 1). While the non-cuprate data set is not nearly as dense, the key trend is analogous to the one found for cuprates. The slope of the $`\lambda _c\sigma _{DC}`$ dependence is also close for both cuprates and non-cuprate materials. The principal difference is that the cuprates universal line is shifted down by approximately one order of magnitude in $`\lambda _c`$. The latter result shows that the superfluid density ($`1/\lambda ^2`$) is significantly enhanced in under-doped cuprates compared to non-cuprate materials with the same DC conductivity.

Possible origins of the $`\lambda _c\sigma _{DC}`$ correlation were recently discussed in the literature . A plausible qualitative account of this effect can be based on the Ferrel-Glover-Tinkham (FGT) sum rule:

$$\frac{c^2}{\lambda ^2}=_{0+}^{\mathrm{\Omega }_C}[\sigma _1^N(\omega )\sigma _1^{SC}(\omega )]𝑑\omega $$
(1)
where $`\sigma ^N(\omega )`$ is the normal state conductivity and $`\sigma ^{SC}(\omega )`$ is the conductivity due to un-paired carriers in the superconducting state at $`\omega >0`$ and $`TT_c`$. For a dirty limit superconductor $`\sigma _1^N(\omega )\sigma _{DC}`$, and Eq. 1 can be approximated as:

$$\frac{c^2}{\lambda ^2}2\mathrm{\Delta }\sigma _{DC}.$$
(2)
Such an approximation is possible because within the BCS model the energy scale $`\mathrm{\Omega }_C`$ from which the condensate is collected is of the order of magnitude of the gap: $`\mathrm{\Omega }_C2\mathrm{\Delta }35kT_c`$. A connection between $`1/\lambda ^2`$, $`\sigma _{DC}`$, and the energy gap is illustrated in the inset of Fig. 1. In the dirty limit the magnitude of $`\sigma _{DC}`$ sets the amount of spectral weight available in the normal state conductivity whereas the magnitude of $`\mathrm{\Omega }_C2\mathrm{\Delta }`$ defines the fraction of this weight which is transferred into condensate at $`T<T_c`$. Therefore, the magnitude of $`\lambda _c`$ can be expected to systematically decrease with the enhancement of the DC conductivity, in accord with the FGT sum rule. Notably, an approximate form (Eq. 2) yields the $`\lambda _c(\sigma _{DC})`$ scaling with the power law $`\alpha =1/2`$ which is close to $`\alpha =0.59`$ seen in Fig. 1.

Interestingly, the sum rule arguments discussed above successfully describe the universal $`\lambda _c\sigma _{DC}`$ scaling in the underdoped cuprates despite the fact that the superconducting energy gap is not well-defined in the interlayer conductivity of these materials. The gap-less response of cuprates is exemplified in left panel in Fig. 2 displaying $`\sigma _1(\omega )`$ data for Pr0.3Y0.7Ba2Cu3O6.95 with $`T_c60`$ K. The conductivity was determined through a Kramers-Kronig analysis of reflectivity measured over the frequency range from 10 cm-1 to 48,000 cm-1. Regardless of the low- and high-energy extrapolations, the $`\sigma _1(\omega )`$ data is unchanged over the energy interval displayed in Fig. 2. The response of the Pr-doped sample is similar to the features found in the conductivity of under-doped oxygen deficient YBa2Cu3O7-δ (YBCO) . Namely, the normal state conductivity is suppressed as the sample is cooled down to $`T_c`$, with a transfer of spectral weight to higher energies. Below $`T_c`$, one does not find any radical changes in the $`\sigma _1(\omega )`$ spectra of Pr-doped YBCO or in the oxygen-deficient YBCO with a similar $`T_c`$. Most importantly, this system, along with all other under-doped compounds, shows significant absorption in the superconducting state so that $`\sigma _1^{SC}>0`$. Therefore, in cuprates only a small fraction of the far-infrared spectral weight is contributing to the condensate. The latter result is in apparent conflict with the assumption $`\sigma _1^{SC}(\omega <2\mathrm{\Delta })0`$, which allows one to reduce Eq. 1 to an approximate form given by Eq. 2. However, the strong condensate density seen in the cuprates located on the universal line can be understood in terms of the dramatic enhancement of the energy scale $`\mathrm{\Omega }_C`$ over the magnitude of the energy gap. This latter conclusion also follows from the explicit sum rule analysis for samples of underdoped La2-xSrxCuO4 (La214) and YBCO materials, suggesting that $`\mathrm{\Omega }_C`$ in these compounds exceeds 0.1-0.2 eV .

The involvement of a broad energy scale into the formation of the condensate in underdoped cuprates is also supported through a comparison of the universal scaling patterns seen for these materials and of a similar pattern detected for non-cuprate superconductors. According to Eqs. 1 and 2, the energy scale associated with the condensate formation for materials on the upper line is of the order of 3 meV which is close to estimates of the gap for most ”conventional” materials in Fig. 1. If Eq. 2 is employed to describe the difference between the upper and the lower lines in Fig. 1, than one is forced to conclude that the corresponding scale for underdoped cuprates is of the order of $``$ 0.3 eV. As pointed out above, this assessment of $`\mathrm{\Omega }_C`$ is supported by the explicit sum rule analysis and also makes $`\mathrm{\Omega }_C`$ the largest energy scale in the problem of cuprate superconductivity. Data points in Fig. 1 for overdoped materials support the notion that the $`\lambda _c\sigma _{DC}`$ plot provides means to learn about the energy scale associated with the condensate formation. Deflection of the over-doped cuprates from the universal line implies that $`\mathrm{\Omega }_C`$ is gradually suppressed with the increased carrier density. This trend is common for Tl2Ba2CuO6+δ (Tl2201), La214 and YBCO materials (see Fig. 1). Integration of the conductivity for all these materials shows that the FGT sum rule is exhausted at energies as low as 0.08 eV .

A quick inspection of the materials in Fig. 1 suggests that the $`3`$ meV scale is observed in systems in which superconductivity emerges out of a FL state whereas the enhanced value of $`\mathrm{\Omega }_C`$ 0.3 eV is found in non-FL underdoped cuprates. The experiments which in our opinion are most relevant to FL versus non-FL classification include quantum oscillations of the low-$`T`$ inter-layer resistivity (and of other quantities) in high magnetic fields . Quantum oscillations is a direct testimony of long-lived quasiparticles which are capable of propagating coherently between the layers. Such quasiparticles are a prerequisite of the Fermi-liquid ground state. On the contrary, quantum oscillations have never been reported for under-doped cuprates. The lack of coherence in the $`c`$-axis transport in these materials indicates that the ground state of cuprates is distinct from the FL picture.

Signatures of FL versus non-FL behavior also can be recognized in the spectra of the $`c`$-axis conductivity. A hallmark of the FL response is the Drude peak seen in $`\sigma _1(\omega )`$ of metals. Notably, this feature is never found in underdoped compounds (forming the non-FL line in Fig. 1). The electronic contribution to $`\sigma _1(\omega )`$ in these materials is usually structureless which is commonly associated with the incoherent (diffusive) motion of charge carriers across the planes. On the contrary, many materials that belong to the upper (FL) line in Fig. 1 demonstrate a familiar Drude behavior. This behavior has been found in Sr2RuO4  and is also shown in our data for inter-plane response of NbSe2 (Fig. 2, right panel), reported here for the first time. In both cases, the width of the peak decreases at low temperatures, which is characteristic of the response of ordinary metals. Earlier measurements of the $`c`$-axis thermo reflectance proved applicability of conventional BCS electrodynamics to the NbSe2 data. As for the over-doped compounds (located in a cross-over region between non-FL and FL lines in Fig. 1) their conductivity is indicative of the formation of the Drude peak (Fig. 2, middle panel), which is becoming more pronounced with the increase of the carrier density.

To summarize the experimental results, we wish to stress the following points: i) two distinct patterns in $`\lambda _c\sigma _{DC}`$ correlation (Fig. 1) originate from a dramatic difference ($`10^2`$) in the energy scale $`\mathrm{\Omega }_C`$ from which the interlayer condensate is collected; ii) the pattern with the typical energy scale of $`3`$ meV is realized in the materials with the coherent FL- type transport between the planes, whereas the one with $`\mathrm{\Omega }300`$ meV is found in cuprate superconductors with the incoherent non-FL response; iii) over-doped cuprates reveal a cross-over between the two behaviors. These results allow us to draw several conclusions regarding the features of the superconducting condensate in different layered systems.
$``$ The symmetry of the order parameter seems to be unrelated to trends seen in the condensate response. Indeed, the upper line in Fig. 1 is formed by $`s`$-wave transition metal dichalcogenides, $`p`$-wave Sr2RuO4 and organic materials for which both $`s`$\- and $`d`$-wave states have been proposed , while $`d`$-wave high-$`T_c`$ materials form the lower line and the crossover region between the lines.
$``$ Electrodynamics of the FL systems at $`TT_c`$ is determined by the magnitude of the gap (and hence by $`T_c`$), in general agreement with the BCS theory. It is therefore hardly surprising that the trend initiated by layered FL materials is also followed in 1-dimensional organic conductors as well as by more conventional systems such as Nb Josephson junctions, bulk Nb, Pb or amorphous $`\alpha `$Mo1-xGex (see Fig. 1).
$``$ We find no obvious connection between the broad energy scale $`\mathrm{\Omega }_C`$ and the critical temperature of the studied superconductors. While scaling of $`\lambda _c`$ by the magnitude of $`T_c`$ does reduce the ”scattering” of the data points, the two distinct $`\lambda _c\sigma _{DC}`$ patterns persist even if such scaling is implemented. In particular, the critical temperature of under-doped La214 materials is nearly the same as that of the several ET-compounds ($`1215`$ K). Nevertheless, the penetration depth is dramatically enhanced in the latter system.
$``$ While the pseudogap state was shown to be responsible for the anomalous superfluid response of the underdoped cuprates, the characteristic pseudogap temperature $`T^{}=90350`$ K is still much lower than our estimate of $`\mathrm{\Omega }_C`$ for these materials.

Dramatic enhancement of the $`\mathrm{\Omega }_C`$ scale well beyond the gap value in underdoped cuprates may be connected with a well-established feature of the electronic spectral function $`A(\omega )`$ in underdoped cuprates. At $`T>T_c`$ the spectral function shows no quasiparticle peak in the antinodal direction; the spectral weight of $`A(\omega )`$ appears primarily in the incoherent channel and is spread out over the energy interval extending beyond 0.3-0.5 eV. The coherent contribution to $`A(\omega )`$ is observed in the superconducting state and the $`T`$-dependence of its weight is similar to that of the superfluid density. These experiments suggest that the emergence of the condensate is associated with the changes in the incoherent region of the spectral function. We stress a quantitative agreement between our estimate of the $`\mathrm{\Omega }_C`$ and the energy region affected at $`T<T_c`$ in the ARPES measurements. Because carrier condensation is connected with the perturbation of the entire spectral function, the scale set by the superconducting gap becomes irrelevant to the problem of the interplane electrodynamics of underdoped cuprates in accord with our findings. The development of coherence in the spectral function at $`T<T_c`$ may be consistent with the lowering of the electronic kinetic energy , a result which was inferred earlier from the sum rule analysis of the interlayer conductivity .

Close correspondence between the properties of the spectral function in the $`(\pi ,0)`$ region and the features of the interlayer conductivity may stem from the peculiarities of the interlayer transfer matrix elements which are maximized for this particular region . Therefore it is natural to look for a common microscopic origin of the anomalies in the $`c`$-axis response and characteristic features of the ARPES spectra at $`(\pi ,0)`$. The latter has been associated with the electron fractinalization at $`T>T_c`$ . An important attribute of scenarios of cuprate superconductivity discussed in Ref. is that at $`T<T_c`$ the total energy is lowered due to reduction of the kinetic energy, in qualitative agreement with the data. A quantitative account of the distinct energy scales associated with the condensate in non-FL materials is a challenge for models attempting to describe superconductivity in the cuprates. This research was supported by the U.S. DOE, NSF and Research Corporation.