text:
# On quantum Griffiths effects in metallic systems

## I  Introduction

The interplay of disorder and quantum criticality is a long-standing and still open problem in condensed matter theory. One aspect of this problem which has received considerable recent attention is the ’quantum Griffiths’ behavior which has been shown to occur near quantum critical points in certain model systems Thill99 ; McCoy69 ; Guo96 ; Rieger96 ; Fisher98 . The model systems in which quantum Griffiths behavior has been unambiguously demonstrated all possess a crucial common feature, namely that in the absence of disorder the critical degrees of freedom exhibit dissipationless, Hamiltonian spin dynamics (indeed typically characterized by dynamical exponent $`z=1`$). However, many systems of experimental importance involve magnetic degrees of freedom coupled to conduction electrons Lohneysen96 ; Mackenzie01 ; Lonzarich99 ; Maple98 , and therefore overdamped dynamics implying a pure-system critical behavior characterized by $`z>1`$. Extension of the theory of quantum Griffiths behavior to this case is therefore an important issue. In a series of papers Neto98 ; Neto00 ; Neto01  Castro-Neto and Jones have argued from various points of view that such overdamped systems exhibit quantum Griffiths behavior similar to that exhibited by undamped systems, and they and others have further argued that this phenomenon is at the heart of the ’non-fermi-liquid’ behavior observed in many heavy fermion materials Maple98 ; Stewart01 .

In this paper we examine the issue of quantum Griffiths behavior in nearly critical systems exhibiting overdamped dynamics, finding that it is essentially nonexistent, being replaced instead by ’superparamagnetic’ behavior. The essence of our analysis is this: in undamped models quantum Griffiths effects arise from an interplay between the low probability of nucleating magnetic ’droplets’ in the paramagnetic state and a low but non-negligible quantum tunnelling of these droplets. In a metallic, dissipative environment there is a strong suppression of tunnelling by dissipation, so that the droplets which dominate the susceptibility behave more or less classically, leading to superparamagnetic behavior rather than quantum Griffiths behavior.

Our results amount to an implementation of ideas outlined in Sherrington75  and to a generalization, to a non-vanishing density of defects, of a previously reported analysis Millis01a  of the ’magnetic droplet’ produced by a single, spatially localized defect, and rely heavily on the results of this previous work. The method used to analyse the dynamics of a distribution of defects is similar to that used in Neto00  and the broad qualitative features of the results we obtain are very similar to those obtained in that work. However, the specifics and the physical implications seem different. The issue is discussed in more detail in the conclusion.

The outline of this paper is as follows. In section II we present the model and the method used in our analysis. In section III we show that the approach reproduces results previously obtained in the dissipationless case. Section IV presents our new results concerning Griffiths-like behavior in systems with overdamped dynamics. Section V is a summary, comparison to other work, and conclusion, and is written so that readers uninterested in the details of the derivations may obtain from it the essence of our results.

## II Model and method of solution

### II.1 Model

The canonical quantum Griffiths problem concerns the effect of weak disorder added to a ’pure’ (non-disordered) system which possesses an Ising symmetry and is tuned to be near a quantum critical point. We consider a system in imaginary time and $`3`$ spatial dimensions (differences occurring for two spatial dimensions warrant a separate treatment, which will be presented elsewhere). The model is described by the action

$$S=S_{\mathrm{static}}+S_{\mathrm{dyn}}+S_{\mathrm{disorder}}$$
(1)
with

$`S_{\mathrm{static}}`$ $`={\displaystyle \frac{E_0}{8\pi }}{\displaystyle _0^\beta }d\tau {\displaystyle }{\displaystyle \frac{d^3x}{\xi _0^3}}[{\displaystyle \frac{\xi _0^2}{\xi ^2}}\varphi ^2(x,\tau )`$
$`+\xi _0^2[\varphi (x,\tau )]^2+{\displaystyle \frac{1}{2}}\varphi ^4(x,\tau )]`$ (2)
Here $`\varphi `$ is a dimensionless scalar order parameter, $`E_0`$ is the basic energy scale of the theory (perhaps of the order of the mean Kondo temperature for a heavy fermion system) and is fixed by normalizing the coefficient of the $`\varphi ^4`$ term to unity, $`\xi _0`$ is the basic length scale (typically of the order of a lattice constant), $`\xi `$ is the magnetic correlation length, and $`\beta `$ is the inverse temperature. It is convenient to define a parameter $`r=(\xi _0/\xi )^2>0`$ which measures distance from criticality. We consider only parameters such that the pure system is in the paramagnetic phase.

We take the disorder to couple to the square of the order parameter via

$$S_{\mathrm{disorder}}=\frac{E_0}{8\pi }_0^\beta 𝑑\tau \frac{d^3x}{\xi _0^3}V(x)\varphi ^2(x,\tau )$$
(3)
and assume it to be Gaussian distributed with correlator ($`\mathrm{}`$ represents average over configurations of the disorder)

$$V(x)V(y)=V_0^2K\left(\frac{xy}{\xi _0}\right)$$
(4)
where the kernel $`K(u)`$ decays on the scale $`u1`$ and satisfies $`d^3uK(u)=1`$. Because we are interested only in length scales $`xy>\xi _0`$ we will take $`K`$ to be a $`\delta `$ function. The dimensionless quantity $`V_0`$ parameterizes the strength of the disorder. Weak disorder corresponds to $`V_01`$.

The dynamic term $`S_{\mathrm{dyn}}`$ is crucial to the quantum criticality described by Eq. 1 and to our subsequent discussions. We consider two cases: (i) dissipationless, $`z=1`$ dynamics, as is the usually assumed in studies of quantum Griffiths behavior, with

$$S_{\mathrm{dyn}}^{(z=1)}=\frac{E_0}{8\pi }_0^\beta 𝑑\tau \frac{d^3x}{\xi _0^3}\left(\frac{\xi _0}{c}\right)^2\left(\frac{\varphi (x,\tau )}{\tau }\right)^2$$
(5)
Here $`c`$ is a characteristic velocity of the undamped excitations, such that $`c/\xi _0`$ is an energy presumably of the order of $`E_0`$.
(ii) Hertz antiferromagnet, $`z=2`$ dynamics, corresponding to the generic antiferromagnetic transition in a fermi liquid:

$$S_{\mathrm{dyn}}^{(z=2)}=S_{\mathrm{dyn}}^{(z=1)}+\frac{T}{8\pi E_0}\underset{\omega _n}{}\frac{\left|\omega _n\right|}{\mathrm{\Gamma }}\frac{d^3x}{\xi _0^3}\left|\varphi (r,\omega )\right|^2$$
(6)
where

$$\varphi (r,\omega _n)=E_0_0^\beta 𝑑\tau \varphi (r,\tau )e^{i\omega _n\tau }$$
(7)
In these conventions the dynamics are dissipative (i.e. dominated by the $`\mathrm{\Gamma }`$ term) if $`\omega <\omega ^{}c^2/(\xi ^2\mathrm{\Gamma })`$, and non-dissipative at higher frequencies. One expects in most systems (and finds for example in a weakly coupled fermi liquid or in the slave boson theory of the Kondo lattice) that all scales are of roughly the same order, i.e. that $`E_0c/\xi _0\mathrm{\Gamma }`$.

### II.2 Method

#### II.2.1 Overview

The dissipative term in Eq.6 corresponds to a long ranged interaction in time and renders available numerical methods prohibitively difficult to apply. To analyze the model defined by Eq. 1 we use simple analytical arguments modelled on those of Ref. Thill99 . We note that the effective dimensionality of the model defined by Eq. 1 is $`d_{\mathrm{eff}}=d+z`$. In this paper we consider only the spatial dimension $`d=3`$ so we are concerned only with models at and above the upper critical dimension $`d_c=4`$, so that quantal and thermal fluctuations of the order parameter in a fixed disorder configuration can be treated by an essentially mean field approximation. The usual fluctuation analysis which justifies the mean-field approximation for $`d_{\mathrm{eff}}>d_c`$ involves a translation-invariant model and fluctuations for which momentum is a good quantum number. Here we must deal with fluctuations in a system whose translation invariance is broken. These were investigated in Refs. Nayaranan01  and Millis01a  and were found not to affect the structure of the static mean field solution when $`d_{\mathrm{eff}}d_c`$ (except for some insignificant changes in some constants).

As noted for example by Sherrington75 , in the presence of the random potential, the crucial feature of the mean field solution is the presence of droplets: regions in which the order parameter is locally non-vanishing. Quantum Griffiths effects then arise from dynamical fluctuations of these droplets; to study them one must estimate the droplet density and tunnelling rate. We use statistical arguments and mean field analysis to estimate the density and an adaptation to the present case of the analysis presented for a droplet produced by a single point defect in Ref. Millis01a  to estimate the tunnelling rate.

#### II.2.2 Probability for the existence of a droplet

The assumption that the model is at or above its upper critical dimension means that mean-field theory is a good starting point Millis01a ; Nayaranan01 . We therefore consider static configurations, $`\varphi (x)`$, which minimize the combination of Eqs.2 and (3). These satisfy

$$\xi _0^2^2\varphi (x)+r\varphi (x)+\varphi (x)^3=V(x)\varphi (x)$$
(8)
If $`V(x)=0`$, then because we assume $`r>0`$ the minimum corresponds to $`\varphi (x)=0`$; however regions in which $`V(x)<0`$ can lead to $`\varphi (x)0`$. In the regions where $`V(x)=const.<0`$, $`\varphi (x)`$ is roughly constant whereas in between these regions $`\varphi (x)`$ decays exponentially. We refer to the regions where $`\varphi `$ is not exponentially small as ’droplets’. If the droplets are reasonably dilute, one may set $`\varphi =0`$ in the exponential tail regions Nayaranan01  and estimate the density of droplets of a given size and mean amplitude.

To motivate our estimate we first consider solving Eq 8 if $`V(x)=\overline{V}`$ for $`\left|x\right|<R`$ and $`V_0=0`$ otherwise. A previous paper Millis01a  considered a special case of this equation, with $`V_0(x)=V\delta ^{(d)}(x)`$ and the solutions found in that work may easily be modified for the present case. In $`d=3`$ one finds that the solution is, roughly (and neglecting unimportant logarithmic factors in the $`x`$-dependence)

$$\varphi (x)=\{\begin{array}{cc}\varphi _0\hfill & \mathrm{for}x<R\hfill \\ \frac{\varphi _0R}{x}e^{(xR)/\xi }\hfill & \mathrm{for}R<x\hfill \end{array}$$
(9)
In other words, the magnetic order induced by the region of attractive $`V`$ is roughly constant inside the region and decays outside it, initially as $`1/x`$ and exponentially for distances larger than a correlation length from the boundary of the attractive potential region. Inserting the above ansatz, Eq. 9, into Eqs. 2 and 3, and minimizing the resulting action with respect to $`\varphi _0`$ yields

$$\overline{V}=\frac{\xi _0^2}{\xi ^2}a(R/\xi )+\varphi _0^2b(R/\xi )$$
(10)
with $`a(x)=1+3/x+3/x^2`$ and $`b(x)=1+3/x4e^{4x}\mathrm{\Gamma }[0,x]`$. These particular forms for $`a`$, $`b`$ depend on the specific potential configuration studied (here $`\overline{V}=const`$ for $`x<R`$ and $`\overline{V}=0`$ otherwise) and on the variational approximation used; but we argue that a generic droplet is described by a similar equation with $`a,b`$ functions which vary on the scale $`R/\xi 1`$ and which tend to unity as $`R/\xi \mathrm{}`$. Also in $`3`$ dimensions $`a(x)1/x^2`$ as $`x0`$ while $`b(x)`$ tends to a constant for $`x<1`$ Millis01a . The precise forms of $`a,b`$ affect only nonuniversal details such as widths of crossover regions. In this paper we shall assume

$`a(x)`$ $`=1+3x^2`$ (11a)
$`b(x)`$ $`=1`$ (11b)
where the $`3`$ arises from the difference in integrating a constant or $`1/r^2`$ over $`r^2dr`$.
One sees from Eq 10 that in order to obtain a solution at all the average potential, $`\overline{V}`$, must be smaller than a (negative) $`R`$dependent critical value,

$$V_c=\frac{\xi _0^2}{\xi ^2}a(R/\xi )$$
(12)
which tends to $`\xi _0^2/\xi ^2`$ as $`R\mathrm{}`$ and to a number of order $`1`$ as $`R\xi _0`$. As is evident from these formulae, the natural scale of the droplets is the correlation length $`\xi `$ which diverges as the quantum critical point is approached.

Eqs. 10 and 12 thus imply that one obtains a droplet in a region of linear dimension $`R`$ only if the average value $`\overline{V}`$ of the potential in that region is larger than a value of the order of $`V_c(R/\xi )`$ ($`V_c`$ is not an exact estimate because it pertains to the idealized disorder configuration discussed above). The standard estimate of the probability of a region of linear dimension $`R`$ with mean potential $`\overline{V}`$ is

$$P(R^3,\overline{V})\frac{\left(R/\xi _0\right)^{3/2}}{\sqrt{\pi }V_0}\mathrm{exp}\left(\left(\frac{R}{\xi _0}\right)^3\left(\frac{\overline{V}}{V_0}\right)^2\right)$$
(13)
and we therefore argue that the density $`N(R^3,\varphi _0^2)`$ of droplets of amplitude $`\varphi _0^2`$ and core size $`R`$, must be proportional to $`\frac{1}{V_0}\mathrm{exp}\left(\frac{R^3\left(\varphi _0^2+V_c(R/\xi )\right)}{V_0^2}^2\right)`$. This argument does not determine the preexponential factors (which involve, e.g. the issue of whether the region of size $`R`$ considered in Eq. 13 is part of a larger region which can sustain a droplet and numerical factors arising from the difference between idealized disorder configuration and typical one, which we have absorbed into $`V_0`$ and $`V_c`$ ). Because some of our subsequent considerations will require an estimate of the preexponential factors, we present the following arguments to fix them.

We begin by making a rough estimate of the fraction of sites contained in droplets (i.e. of the fraction of sites having a $`\varphi _0^2>0`$), as a function of distance from criticality. As noted above, in principle within mean field theory $`\varphi _0^2`$ is non-vanishing everywhere, but we neglect the regions where it is exponentially small, in other words we set $`\varphi _0=0`$ in the ’inter-droplet’ regions. To perform the estimate we coarse-grain the theory to the scale $`\xi `$. A given correlation volume $`\xi ^3`$ will have a non-vanishing $`\varphi _0^2`$ if the potential averaged over the droplet volume, $`V_{\mathrm{ave}},`$ is larger than $`V_c(\xi )a(1)\xi _0^2/\xi ^2`$. From Eq 13 we see that the probability $`P_\varphi `$ that a given correlation volume will have a non-vanishing $`\varphi _0`$, i.e. a $`V<V_c(\xi )`$, is (recall $`V_c<0`$)

$$P_\varphi =\frac{1}{2}\left(1\mathrm{erf}\left[\left(\frac{\xi }{\xi _0}\right)^{3/2}\frac{V_c(\xi )}{V_0}\right]\right)$$
(14)
where $`erf`$ is the error function.

Clearly, a picture of independent droplets must break down if $`P_\varphi `$ exceeds the percolation probability $`P_{\mathrm{perc}}`$ at which the set of correlation volumes with non-vanishing $`\varphi _0`$ percolate. Use of Eq. 14 and the estimate for three dimensional cubic lattices $`P_{\mathrm{perc}}0.2`$ shows that percolation will have occurred by the time $`\xi `$ exceeds $`\xi _{\mathrm{perc}}2.8a(1)^2/V_0^2`$. ($`\xi _{\mathrm{perc}}`$ is an underestimate because droplets larger than $`\xi `$ may occur). These estimates also show that the natural scale for $`\xi `$ is $`V_0^2`$ and strongly suggest that the probability that a given site is in a droplet (of any size) is a function only of the combination $`\xi V_0^2`$.

We therefore argue that the prefactors in the droplet density must be such that the total probability of finding a site in a droplet, $`P_{\mathrm{tot}}=\xi _0^6𝑑R^3𝑑\varphi _0^2R^3N(R^3,\varphi _0^2)`$ must be a function only of $`\xi V_0^2`$ and must be of the order of $`P_{\mathrm{perc}}`$ when $`\xi `$ is of the order of $`\xi _{\mathrm{perc}}`$. This implies

$$N[R^3,\varphi _0^2]=\frac{R^{9/2}}{C_{V_0}V_0}\mathrm{exp}\left(\frac{R^3\left(\varphi _0^2+V_c(R/\xi )\right)}{V_0^2}^2\right)$$
(15)
where the factor of $`R^{9/2}`$ ensures the correct scaling with $`\xi `$ and the numerical factor $`C_{V_0}11.25`$ ensures that when $`\xi =\xi _{\mathrm{perc}}`$, we have $`P=P_{\mathrm{perc}}`$. We emphasize that these formulae are phenomenological and must in particular break down when $`\xi `$ approaches $`\xi _{\mathrm{perc}}`$.

It is convenient to adopt a dimensionless system of units in which

$`R`$ $`=y\xi `$ (16)
$`\varphi _0`$ $`=f{\displaystyle \frac{\xi _0}{\xi }}`$ (17)
$`\xi `$ $`=u{\displaystyle \frac{\xi _0}{V_0^2}}`$ (18)
for which

$`N(y^3,f^2)dy^3df^2={\displaystyle \frac{y^{9/2}}{C_{V_0}u^{1/2}\xi ^3}}`$
$`\times \mathrm{exp}\left({\displaystyle \frac{y^3\left(f^2+a(y)\right)^2}{u}}\right)dy^3df^2`$ (19)
The factor of $`\xi ^3`$ expresses the fact that if the probability of a given site being in a droplet is a function only of $`u`$, then the density of droplets must be smaller by an extra factor of the typical droplet volume $`\xi ^3`$.

#### II.2.3 Tunnelling of the droplet for undamped, $`z=1`$, dynamics

We now estimate the rate $`\omega _{\mathrm{tun}}`$ at which a droplet characterized by the mean amplitude $`\varphi _0`$ and length scale $`R`$ tunnels in the case of undamped, $`z=1`$ dynamics by performing a variational instanton calculation using Eqs.2, 5 and the solution Eq. 9. In the simplest estimate one assumes that the droplet maintains its shape while collapsing and re-forming. To estimate the action associated with this process we write the droplet solution as

$$\varphi (x,\tau )=\varphi (x)\eta (\tau )$$
(20)
Substitution into Eqs. 2, 5 leads to

$$S_{\mathrm{inst}}=S_{\mathrm{kin}}+S_{\mathrm{barrier}}$$
(21)
$`S_{\mathrm{kin}}`$ involves the integral of $`\left(_t\varphi \right)^2`$ over the droplet and as noted in Ref. Millis01a  involves the $`1/r`$ ’tail’ of the droplet in a crucial manner; in contrast, the cost $`S_{\mathrm{barrier}}`$ of creating the instanton does not. One obtains

$`S_{\mathrm{kin}}^{(z=1)}`$ $`=C_{\mathrm{kin}}\xi f^2y^3a^{}(y){\displaystyle \frac{d\tau }{E_0}\left(\frac{\eta }{\tau }\right)^2}`$ (22)
$`S_{\mathrm{barrier}}`$ $`=C_{\mathrm{barrier}}\xi ^1f^4y^3b^{}(y)`$
$`\times {\displaystyle }E_0d\tau (2\eta (\tau )^2+\eta (\tau )^4)`$ (23)
Here $`C_{\mathrm{kin}}`$ and $`C_{\mathrm{barrier}}`$ are nonuniversal constants. $`C_{\mathrm{kin}}`$ involves the square of the ratio $`E_0/(c/\xi _0)`$ of the basic energy scale to the kinetic (or zone boundary magnon) energy while $`C_{\mathrm{barrier}}`$ is just a number. In the approximation we have employed $`C_{\mathrm{kin}}=E_0^2\xi _0^2/c^2`$ and $`C_{\mathrm{barrier}}=1`$. The functions $`a^{}=x^2𝑑x\varphi _0(x)^2`$ and $`b^{}=x^2𝑑x\varphi _0(x)^4`$ are functions with behavior similar to $`a,b`$; in our explicit calculations we set $`a^{}=a/3=1/3+y^2`$ and $`b^{}=b/3=1/3`$ for simplicity; again different choices affect only nonuniversal details.

The action associated with one instanton may now be determined by a standard minimization of Eqs 22,23 and is

$$S_{\mathrm{inst}}^{(z=1)}=S_1d(y)f^3y^3$$
(24)
For the present model in the present approximation the nonuniversal constant $`S_1=\sqrt{C_{\mathrm{kin}}C_{\mathrm{barrier}}/3}`$and $`d(y)=3\sqrt{a^{}(y)b^{}(y)}=\sqrt{a(y)}`$ where the last equality follows from our simplifying assumptions $`a^{}=3a`$ and $`b^{}=3b`$. The value of $`S_1`$controls the width of the crossover regime before the universal behavior is reached, and is linearly proportional to $`E_0\xi /c`$

The tunnelling rate is then given by

$$\omega _{\mathrm{tun},z=1}=\omega _0e^{S_{\mathrm{inst}}(\tau _0)}.$$
(25)
Here, $`\omega _0`$ is an attempt frequency presumably of order $`E_0`$ whose value is beyond the scope of this theory.

To conclude this section we briefly estimate the action associated with a different tunnelling mechanism, namely nucleation of a domain wall. For small droplets (’core size’ $`R`$ less than $`\xi `$) the important process was shown to be collapse and reformation of the entire droplet Millis01a . We therefore need consider only the case $`R\xi .`$ We observe that by expanding about the static uniform solution one obtains a domain wall with width $`W\varphi _0^1.`$ The kinetic term associated with the domain wall motion therefore has one fewer factor of the small quantity $`\varphi _0f/\xi V_0^2f`$, leading to a larger action and hence a smaller rate, in the weak disorder, near criticality limit. We note in passing that for $`\varphi _01`$ the powers of $`R`$ will be the same as we have considered but the extra factor of $`\varphi _0^1`$ will work in the other direction, favoring domain wall motion.

#### II.2.4 Tunnelling of the droplet for overdamped, $`z=2`$, dynamics

For overdamped dynamics two important differences occur. First, as shown in Refs. Leggett87 ; Neto00 ; Millis01a  the damping changes the action associated with a single instanton, strongly suppressing the bare tunnelling rate relative to that found for undamped dynamics. Essentially, the tunnelling is limited by the droplet’s ability to move through a viscous medium rather than by its ability to climb over a barrier. Second, and much more important, the overdamped dynamics leads to a long-ranged (in time) instanton-instanton interaction, which reduces the tunnelling rate further and indeed drives it to zero if the damping exceeds a critical value, as noted by previous authors Neto00 ; Millis01a .

To calculate the effects of damping we insert the ansatz, Eq. 20 into Eq 6. The new term arising from the overdamped dynamics is

$$S_{\mathrm{diss}}=\frac{\gamma }{4}𝑑\tau 𝑑\tau ^{}\frac{d\eta }{d\tau }\frac{d\eta }{d\tau ^{}}\mathrm{ln}(\frac{(\tau \tau ^{})^2+\tau _m^2}{\tau _m^2})$$
(26)
with $`\tau _m`$ a ’microscopic’ time of the order of $`\omega ^{}=c^2/(\xi _0^2\mathrm{\Gamma }`$. The net dissipative coefficient $`\gamma `$ is given for the Hertz antiferromagnet by Millis01a

$$\gamma =\frac{E_0}{4\pi \mathrm{\Gamma }}\frac{d^3x}{\xi _0^3}\varphi _0(x)^2=c_\gamma f^2y^3a^{\prime \prime }(y)\xi /\xi _0$$
(27)
The approximations employed in the previous section imply that the nonuniversal constant $`c_\gamma =E_0/\mathrm{\Gamma }`$ and $`a^{\prime \prime }=a^{}`$. In a generic system one expects all scales to be of roughly the same order, so that in particular $`c_\gamma `$ is expected to be of order unity.

The estimate of $`\gamma `$ is subject to the important caveat that the electron bath which causes the dissipation can penetrate the entire droplet. A reasonable estimate of the penetration depth, $`L_p`$, may be obtained by dividing the electron velocity, $`v_\mathrm{F}`$, by the magnitude of the order parameter; in rescaled units $`L_p/\xi v_\mathrm{F}/(E_0f)`$. We shall see below Eq. 48 that the parameters are such that the electrons can penetrate the entire droplet.

We have not been able to solve analytically for the instanton; instead we estimate the action by inserting the variational ansatz

$$\frac{d\eta }{d\tau }=2\frac{\mathrm{\Theta }(\tau _0^24\tau ^2)}{\tau _0}$$
(28)
into Eqs 6,23 obtaining $`S=S_{kin}+S_{diss}+S_{barrier}`$ with

$$S_{kin}=\frac{2C_{kin}\xi y^3a^{}(y)f^2}{E_0\tau _0\xi _0}$$
(29)
$$S_{diss}=2c_\gamma \frac{\xi }{\xi _0}y^3a^{}(y)f^2ln(c_d\tau _0/\tau _m)$$
(30)
$$S_{barrier}=\frac{2}{15}C_{barrier}\frac{\xi _0}{\xi }y^3b^{}(y)f^4E_0\tau _0$$
(31)
where $`ln(c_d)=_{1/2}1/2𝑑x𝑑yln(1+(xy)^2)0.1152\mathrm{}`$.

Minimization over the instanton duration then leads to

$$1=\frac{c_\gamma }{C_{kin}}\tau _0E_0+\frac{C_{barrier}}{15C_{kin}}\frac{b^{}f^2\xi _0^2}{a^{}\xi ^2}(\xi _0/\xi )^2(E_0\tau _0)^2$$
(32)
As previously remarked, we expect the ratios of the various dimensional parameters to be of the order of unity; also as we shall see below, in this problem the important droplets have $`f\xi ^{1/2}`$, so that provided $`\mathrm{\Gamma }`$ is less than a number of the order unity times $`\frac{\xi c}{f\xi _0}\xi ^{3/2}`$ (within our approximations the precise numerical factor is $`\sqrt{15}`$) the $`\tau _0^2`$ term is negligible and one has

$$\tau _0=\frac{\mathrm{\Gamma }\xi _0^2}{c^2}$$
(33)
and thus

$$S_{\mathrm{inst}}^{(z=2)}=c_\gamma C_2f^2\xi \xi _0^1y^3a(y)$$
(34)
where $`C_2`$ (=2.283… in the present approximations) is a numerical factor of the order of unity arising from combining the factors in Eqs 29 and 30).

We observe that for the value of $`\tau _0`$ given in Eq 33, the term written in Eq 34 is larger than $`S_{\mathrm{barrier}}`$ (Eq 23) by two powers of the correlation length (provided that the quantity $`f`$ is of order unity or less, as is the case for the situations considered here.) Thus, in the metallic case and near to criticality, the difficulty in tunnelling arises from moving through the viscous medium, not climbing over the barrier. This result was noted previously Millis01a .

The bare tunnelling amplitude is thus

$$\omega _{\mathrm{bare}}^{(z=2)}=\omega _0e^{S_{\mathrm{inst}}^{(z=2)}}$$
(35)
and is much smaller than in the dissipationless case, because of the factor $`f\xi `$ in the argument of the exponential.

The standard macroscopic quantum tunnelling arguments Leggett87 ; Neto00 ; Millis01a  imply that the instanton-instanton interaction renormalizes the bare tunnelling rate so that if $`\gamma <1`$ then the $`T=0`$ tunnelling rate is

$$\omega _{\mathrm{tun}}=\omega _0\left(\frac{\omega _{\mathrm{bare}}}{\omega _0}\right)^{\frac{1}{1\gamma }}$$
(36)
whereas if $`\gamma >1`$ tunnelling stops at $`T=0`$. We see from Eq. 27 that $`\gamma `$ is a strong function of the droplet size and amplitude; droplets which may tunnel (i.e have $`\gamma <1`$) have a very weak amplitude even in rescaled units: $`f\xi ^{1/2}`$.

Eq 36 is a zero temperature result. At $`T>0`$ the ’Caldeira-Leggett’ renormalization is temperature dependent. The key question for this paper is the temperature at which $`\omega _{tun}(T)<T`$. If $`\gamma >1`$ then $`\omega _{tun}(T)<T`$ at all $`T<E_0`$, implying that the droplet behaves classically at all $`T`$. If $`\gamma <1`$ then the usual arguments shows that $`\omega _{tun}(T)`$ drops below $`T`$ when $`T`$ becomes greater than $`\omega _{tun}(T=0)`$, so that Eq 36 gives the temperature scale separating a high-T region, in which the droplet behaves classically, from the low-T region, in which it behaves quantum mechanically.

## III Estimate of Quantum Griffiths Behavior

### III.1 Overview

The standard Griffiths estimate is that a droplet of magnetic moment $`M_d=d^3re^{i\stackrel{}{Q}\stackrel{}{r}}\varphi _0(r)`$ ($`\stackrel{}{Q}`$ is the ordering vector) and tunnelling frequency $`\omega _{\mathrm{tun}}[R,\varphi _0^2]`$ gives rise to a susceptibility $`\chi _d`$ proportional to $`M_d^2/\left(\omega _{\mathrm{tun}}+T\right)`$. The susceptibility of a system with a distribution of droplets is then given by

$$\chi (T)=d^3Rd^2\varphi _0\frac{N(R^3,\varphi _0^2)M_d^2[R,\varphi _0]}{\omega _{\mathrm{tun}}[R,\varphi _0^2]+T}$$
(37)
For a droplet in an antiferromagnetic system, we find $`M_d`$ is a random function with magnitude $`\varphi _0R`$–the term proportional to $`R`$ comes from the boundary of the droplet, where the order parameter amplitude is dropping and the cancellation over one unit cell of the antiferromagnetic order is not complete. A different dependence would change prefactors but not affect our results crucially.

It is convenient to introduce an explicit integral over frequency, writing

$$\chi (T)=\xi ^3𝑑\omega \frac{I(\omega )}{\omega +T}$$
(38)
so that after conversion to dimensionless units we have

$$I(\omega )=𝑑y^3𝑑f^2(\xi ^3N(y^3,f^2))f^2y^2\delta (\omega \omega _{\mathrm{tun}}(y,f))$$
(39)
The prefactor $`\xi ^3`$ in $`\chi `$ arises because each droplet has magnetic moment of the order of unity and the density of droplets is $`\xi ^3`$. The quantity $`\xi ^3N`$ has no explicit dependence on $`\xi `$ (see Eq. 19).

We will use the delta function to eliminate the $`f`$ integral in $`I`$ and perform the integration over y either numerically or via an extremal value argument.

### III.2 z=1

Using Eq. 25 yields

$$f(\omega ,y)=\left(\frac{\mathrm{ln}\left(\frac{\omega _0}{\omega }\right)}{S_1d(y)}\right)^{1/3}\frac{1}{y}$$
(40)
Substituting this result into Eq.39 yields

$$I(\omega )=\frac{2\mathrm{ln}^{1/3}(\omega _0/\omega )}{\omega S_1^{4/3}}_0^{\mathrm{}}\frac{\xi ^3N(y^3,f^2(\omega ,y))dy}{d(y)^{4/3}}$$
(41)
where $`N(y,f(\omega ,y))`$ is $`N(y,f)`$ (Eq 15) with $`f`$ given by Eq 40.

In the limit of very low frequency one may use asymptotic methods to analyse the integral in Eq 41; the extremum is at

$$y_{\mathrm{max}}=\frac{\mathrm{ln}^{1/3}\left(\frac{\omega _0}{\omega }\right)}{\sqrt{3}S_1^{1/3}}$$
(42)
Substitution leads to

$$\chi (T)\frac{1}{\xi ^3C_{V_0}}\frac{1}{T^{1d_{\mathrm{asympt}}}}$$
(43)
with (restoring units)

$$d_{\mathrm{asymp}}(\xi )=\frac{16}{3\sqrt{3}S_1}\frac{1}{\xi V_0^2}$$
(44)
This is the familiar quantum Griffiths result: if one is sufficiently close to the pure system critical point ($`d(\xi )<1`$) then the susceptibility diverges, with degree of divergence characterized by an exponent which approaches unity proportional to one power of the inverse correlation length.

Note that the prefactor in Eq 43 rapidly vanishes as criticality is approached, so although the susceptibility diverges more strongly, the amplitude of the divergence decreases. Note further that in the asymptotic limit, $`f1`$ so that the mean order parameter density (integrated order parameter divided by droplet volume) is of the order of $`\xi ^1`$. Thus the picture that emerges is of large, weak droplets.

We have evaluated

$$d_{\mathrm{eff}}(\omega )=1+\frac{d\mathrm{ln}\left[I(\omega )\right]}{d\mathrm{ln}(\omega )}$$
(45)
via a numerical computation of Eq.41. Fig. 1 shows $`d_{\mathrm{eff}}(\xi ,\omega )`$ as a function of $`\xi `$ for $`\omega =0.001\omega _0`$ and several different values of the non-universal parameter $`S_1`$ (solid lines) along with the asymptotic limit estimates from Eq 44. We observe that for these low frequencies and not too long $`\xi `$ the asymptotic limit provides a reasonable (but not perfect) estimate of the effective exponent: relative corrections are of the order of $`(\xi /ln(\omega _0/\omega )^{2/3})`$. We see also that depending on the value of the non-universal parameter $`S_1`$, the effective exponent may remain above the critical value of unity (corresponding to a non-divergent susceptibility) until $`\xi `$ becomes of the order of $`\xi _{perc}`$. For $`\xi `$ of the order of $`\xi _{perc}`$ the standard quantum griffiths approximation (independent droplets) breaks down, and one must deal instead with the critical singularities appropriate to a phase transition in a disordered system; in other words with the still unsolved problem of the mixing of quantum critical and quantum griffiths singularities.

### III.3 z=2

For overdamped dynamics, some droplets will have $`\gamma >1`$ and therefore will not tunnel at all at $`T=0.`$ The function $`I(\omega )`$ will thus have a contribution proportional to $`\delta (\omega )`$ leading to the $`1/T`$ behavior expected of classical droplets. For those droplets which do tunnel we must use Eq 36 in Eq 39. We write

$$I(\omega )=I_0\delta (\omega )+I_{rest}(\omega )$$
(46)
with $`I_{rest}`$ given by Eq 39 and $`I_0`$ by

$$I_0=𝑑y^3𝑑f^2f^2y^2N(y^3,f^2)\mathrm{\Theta }(\gamma (y,f)1)$$
(47)
From Eq 38 we see that if $`I_0`$ is appreciable, then $`\chi 1/T`$: this is the superparamagnetism expected from essentially classical droplets.

We begin by estimating $`I_0`$. The $`\mathrm{\Theta }`$ function limits the $`f`$ integration to

$$f^2>f_{\mathrm{min}}^2(y)=\frac{\xi _0}{c_\gamma \xi y^3a(y)}$$
(48)
Note that for large $`\xi `$, $`f_{\mathrm{min}}a(y)`$. Further, the typical scale for $`f`$ is $`\xi ^{1/2}`$ so that the penetration depth $`L_p`$ of electrons into the droplet is large: $`L_p/\xi \xi ^{1/2}`$ so the assumption that electrons penetrate the droplet is indeed self consistent.

Use of Eq 48 in Eq 47 gives

$`I_0(\xi )`$ $`={\displaystyle \frac{3\sqrt{\pi }}{2C_{V_0}}}{\displaystyle _0^{\mathrm{}}}y^2dy[{\displaystyle \frac{\sqrt{u}e^{y^3(f_{\mathrm{min}}^2(y)+a(y))^2/u}}{\sqrt{\pi }y^{3/2}}}`$
$`+a(y)(\mathrm{erf}\left({\displaystyle \frac{y^{3/2}\left(f_{\mathrm{min}}^2(y)+a(y)\right)}{u^{1/2}}}\right)1)]`$ (49)
$`I_0`$, normalized to the total weight in $`I`$, $`𝑑\omega I(\omega )`$ is plotted in Fig. 2 as a function of $`\xi `$ for different values of the disorder strength $`V_0`$. We see that the factor of $`\xi ^1`$ in Eq 48 means that as criticality is approached, almost all of the weight in the droplet probability distribution is in droplets which do not tunnel.

For the droplets which are able to tunnel at frequency $`\omega `$, we find from Eqs. 34 and 36 that

$$f_{z=2}^2(y)=\frac{\xi _0}{c_\gamma \xi y^3a(y)}\frac{\mathrm{ln}\left(\frac{\omega _0}{\omega }\right)}{C_2+\mathrm{ln}\left(\frac{\omega _0}{\omega }\right)}$$
(50)
Note that in contrast the expression for $`f`$ in the $`z=1`$-case shown in Eq. 40, in the $`z=2`$ case, $`f`$ does not diverge as $`\omega 0`$. As in the $`z=1`$ case considered above, one obtains an expression for $`I_{\mathrm{rest}}(\omega )`$ by substituting the result for $`f`$ into Eq 39 yielding

$`I_{\mathrm{rest}}(\omega )`$ $`={\displaystyle \frac{3\xi _0^3}{\omega C_Vc_\gamma ^3\xi ^3}}{\displaystyle \frac{C_2\mathrm{ln}\left(\frac{\omega _0}{\omega }\right)}{\left(C_2+\mathrm{ln}\left(\frac{\omega _0}{\omega }\right)\right)^3}}`$
$`\times {\displaystyle }{\displaystyle \frac{dyy^{1/2}e^{\frac{y^3\left(f_{z=2}^2(y)+a(y)\right)^2}{u}}}{\left(y^3a(y)\right)^2}}`$ (51)
The resulting expression is to good accuracy proportional to $`1/\omega `$ times logarithms. The physics is that even the average of droplets which are able to tunnel is dominated by those droplets on the verge of freezing, leading again to a superparamagnetic contribution to the susceptibility.

We have numerically evaluated the integral in Eq. 51 for parameters such that $`I_0`$ is not too large. Sample results are shown in Fig. 3, which plots the quantity $`J_{rest}=\omega I_{rest}`$ for a relatively small value of the damping. The frequency dependence is a consequence of the logarithmic factors in Eq 51; the non-vanishing intercept as $`\omega 0`$ means that up to logarithmic corrections the contribution to the susceptibility arising from this term is $`1/T`$.

## IV Conclusion

This paper presents an investigation of the possibility of quantum Griffiths effects in three dimensional metallic system which, when ’pure’ (non-disordered) is near an antiferromagnetic quantum critical point with Ising symmetry. For comparison we present also a parallel investigation of quantum Griffiths effects in a model of an insulating system near a similar critical point. The key feature of metallic systems is the dissipative dynamics arising from the particle-hole continuum of electrons; in the model insulating system the dynamics are undamped. Comparison of the two calculations shows that dissipation suppresses quantum Griffiths effects completely, leaving instead an effectively superparamagnetic behavior.

A simple precis of our results follows. Quantum Griffiths effects are a consequence of randomness: essentially, in a random system which is on average in the paramagnetic phase, regions (’droplets’) may occur in which the randomness pushes the system locally to the ordered side of the phase diagram, so that local formation of an order parameter is favored. In certain circumstances (first noted by McCoy McCoy69 ) these droplets may dominate the response. In this situation one may approximately write the susceptibility, $`\chi `$, as an average over droplets times a susceptibiltiy for each droplet, i.e.

$$\chi =_{\mathrm{droplets}}P(\mathrm{droplet})\chi _{\mathrm{droplet}}$$
(52)
We have used simple extremal statistics arguments (similar to those used by Thill and Huse Thill99 ) to estimate the droplet probability distribution $`P(\mathrm{droplet})`$ and an extension of earlier work which studied a particular class of droplets Millis01a  to obtain the susceptibility $`\chi _{\mathrm{droplet}}`$ of a given droplet. We were then able to perform the average over droplets and obtain an estimate for the susceptibility.

This method reproduces the essential features of the standard results for quantum Griffiths effects in undamped (insulating) systems, namely that the low $`T`$ behavior of the susceptibility is governed by a new exponent $`d_{\mathrm{eff}}`$ given by the product of the inverse correlation length $`\xi ^1`$ and inverse mean square disorder amplitude $`V_0^2`$ and a non-universal number (which we estimate for the particular model we consider). A divergent susceptibility results when $`d_{\mathrm{eff}}`$ becomes less than unity, and the results are functions only of $`\xi V_0^2`$. We note one additional interesting finding. The standard arguments which produce the standard quantum Griffiths results are based on a picture of dilute ’droplets’ and apply only if the $`\xi `$ is not too large (otherwise the droplets percolate, and an isolated droplets picture fails). For the model we consider we obtain an estimate for the critical value of $`\xi `$, and find that depending on the value of the non-universal factor in $`d_{\mathrm{eff}}`$, droplets may reach the percolation point before the Griffiths exponent drops below unity. In other words, in the models we consider the existence of a quantum Griffiths regime (which one may somewhat imprecisely define as a divergent susceptibilty arising from fluctuations of isolated droplets) is not guaranteed–it may or may not occur depending on the value of a non-universal coefficient. Sufficiently near a critical point a regime of divergent susceptibility does of course occur, but the proper theory of this regime would have to go beyond the model of isolated droplets and treat correctly the mixing of critical and griffiths singularities.

We also found that for undamped systems near antiferromagnetic critical points the amplitude of the divergent term in the susceptibility vanishes rapidly as criticality is approached, indeed as $`\xi ^3`$, essentially because each relevant droplet has a magnetic moment of the order of unity and as criticality is approached the droplets get larger in size but fewer in number.

The main new result of our work, however, pertains to metallic systems with overdamped (dissipative) dynamics. For these systems (i.e. for quantum critical phenomena in metals) the answer is entirely different. The physics in the undamped case is a balance between the probability of a droplet occurring (which vanishes rapidly as the droplet size or amplitude increases) and $`\chi _{\mathrm{droplet}}`$, which is of the order of the inverse of the quantum tunnelling rate of the droplet and diverges rapidly as the droplet size or amplitude increases. The effect of dissipation is to strongly decrease the tunnelling rate, and indeed to drive it to zero for droplets larger than a particular, amplitude-dependent, size. For relevant parameters we find that a non-vanishing density of droplets does not tunnel at $`T=0`$; these give rise to a ’superparamagnetic’ ($`\chi 1/T`$) susceptibility rather than a quantum griffiths (continuously varying exponent) behavior. For those droplets which do behave quantum mechanically, the effect of dissipation on the tunneling rate is found to change the balance between probability and $`\chi _{\mathrm{droplet}}`$ dramatically. We find that even considering only the droplets which can tunnel quantum mechanically, those which dominate the integral for $`\chi `$ are those which are right on the edge of classical (non-tunnelling) behavior, leading again to superparamagnetism rather than to quantum Griffiths behavior. We also find that the dependence on parameters is different: in the undamped case, apart from prefactors the mean square disorder strength $`V_0^2`$ and the correlation length enter via the combination $`\xi V_0^2`$. In the damped case additional factors of $`\xi `$ occur which drive the system more rapidly to classical behavior.

Our results raise questions about the claims Maple98 ; Neto98 ; Neto00 ; Neto01  that quantum Griffiths effects are important in heavy fermion materials, which are precisely three dimensional metals with Ising symmetry, typically near antiferromagnetic quantum critical points. Ref. Maple98  contains a phenomenological description of data. If the theoretical results presented here are accepted, then these data require a different, non-Griffiths interpretation. Ref. Neto98  argued that a disordered system near a quantum critical point could be mapped onto the disspationless Ising model in a transverse field; the results of the present paper and of Millis01a  indicate on the contrary that dissipation is essential.

Ref. Neto01  uses a novel variant of a technique introduced by Dotsenko Dotsenko99  to study essentially the same model as is studied here. A rather different result was obtained, namely that quantum Griffiths effects can be important in a reasonable range of the phase diagram even in the metallic case. We outline the differences between the results found here and those of Ref. Neto01 . The method introduced by Dotsenko Dotsenko99  and used by Ref Neto01  begins from a classical theory defined by a functional integral with action given by the static term in Eq. 1 and evaluates the disorder-average by the replica method. Whereas other workers Boyanovsky84 ; Nayaranan01  then used the replicated field theory to derive scaling equations for variables including the mean disorder strength, Dotsenko argued that one should look for spatially localized energetically unstable configurations of the replicated field theory, which correspond to local maxima of the replicated action and are to be identified with the ’droplets’ discussed above. Dotsenko shows that the leading nonanalytic contribution to the free energy in the vicinity of an assumed $`T>0`$ critical point comes from droplets with size of the order of the magnetic correlation length, $`\xi `$; we refer to these henceforth as ’typical droplets’. The authors of Ref. Neto01  assume that the $`T0`$ limit of this classical theory may be straightforwardly taken, and then add to this theory estimates of the dynamics of ’typical’ droplets. The results reported in Ref. Neto01  disagree in a number of specific details with the results presented here, including for example the way in which the bare tunnelling rate is estimated. The most important difference, however, is in the interpretation of the results. Ref. Neto01  argues that one should identify the boundary of the Griffiths region with the value of $`\xi ^2`$ at which a ’typical droplet’ ceases to tunnel. Our analysis, which involves averaging over all droplets, indicates that independent of whether the ’typical droplet’ (however defined) may tunnel, the susceptibility is dominated by droplets which are at or beyond the edge of ceasing to tunnel; these give an essentially ’superparamagnetic’ ($`\chi 1/T`$) behavior, instead of the continuously varying exponent characteristic of quantum Griffiths behavior.

Ref. Neto00  presented a detailed analysis of a different model in which spins are added to a pure system which itself is far from any critical point. In this model the phase transition is disorder-driven: it occurs when the density of added spins is high enough that these order; whereas our interest here has been in models in which even the non-disordered system is near a critical point. Further in the model studied in Neto00  the way the disorder is introduced means that the local spin amplitude $`\varphi _0`$ (c.f our Eq. 9) is always of order unity, whereas in our treatment the local spin amplitude may be considerably smaller. An approximate mapping between the model considered in Neto00  and the one considered here may be obtained by setting our parameters $`\varphi _0`$ and $`\xi `$ equal to unity and considering the behavior as the disorder strength $`V_0`$ is increased (whereas we consider a fixed $`V_0`$ and study the behavior as $`\xi `$ is increased).

Although specific details differ, in a broad qualitative sense results obtained in Ref Neto00  are similar to those obtained here. In particular, Ref Neto00  states that at sufficiently low temperature dissipation will suppress the quantum griffiths behavior. However, Ref Neto00  argued that in an extremely wide temperature regime could exist in which behavior characteristic of the undamped system occurs, whereas in the model we consider, for any reasonable parameters there is no such temperature regime. A crucial point is that Neto00  focussed on model parameters such that the damping coefficient was extremely weak (i.e in our notations (see below Eq. 27) they took $`c_\gamma <<1`$). In this limit, it is plausible that there is a temperature regime in which behavior characteristic of the undamped model may occur, before finally a crossover occurs to a regime (similar to the one we considered) in which damping is important. Important avenues for future investigation include more detailed studies of the crossovers between the weak-damping and order unity damping cases and between the disorder driven-criticality effects studied in Neto00  and the pure system criticality-driven effects studied here, as well as determination of the damping coefficient values appropriate to the heavy fermion materials of interest.

Our work has the following implications for experiment. First, the canonical quantum Griffiths effects are due to weak disorder added to a pure critical point. We have shown that in the limit of weak disorder and a pure critical point described by the Hertz theory Hertz74 ; Millis95 , the dissipation characteristic of metallic systems changes the quantum Griffiths singularities into a kind of superparamagnetic behavior. In other words, as a matter of principle the canonically defined quantum Griffiths behavior should not be observable in metals near magnetic quantum critical points. This suggests that claims Maple98  to have observed quantum Griffiths behavior in heavy fermion systems should be treated with caution (at least for systems with Ising symmetry). Further, we showed that the droplets which dominate the susceptibility can tunnel only when the system is not close to criticality, and in these cases the droplet size is not much larger than the basic scale of the theory. Thus if the susceptibility is dominated by the tunnelling of droplets, the picture which emerges is more similar to the ’Kondo disorder’ picture of Bernal95 ; Dobrosavljevic97  than it is to the conventionally-defined quantum Griffiths picture. Indeed, the experimental claims involve ’heavy fermion’ systems where the interaction which favors a non-magnetic phase is the Kondo effect. As noted by many authors Dobrosavljevic97 , the fact that Kondo temperatures are exponentially sensitive to system parameters means that a slight variation in system parameters can lead to a wide variation in Kondo temperatures. The canonical assumption of weak disorder which we and others Thill99 ; Rieger96 ; Neto01 ; McCoy69  have made may not be valid for these systems. The interplay between quantum criticality and a broad distribution of disorder should be treatable by the methods introduced here, and seems worth examining.

A second point is that the very slow dynamics of the droplets makes it much easier for them to order. Further, in a metallic system the droplet-droplet interactions are of long range (see, e.g. Narozhny01  for a discussion in the context of the two dimensional metal insulator transition). For this reason we expect that in the presence of disorder the actual phase transition at which long ranged order sets in is an essentially classical affair, in which droplets lock together when the temperature becomes lower than some droplet-droplet coupling.

A third point, perhaps relevant beyond the present context, is that (as seen for example in Eq 34) dissipation can have a crucial effect on bare tunnelling rates: in the metallic problem we considered the cruical impediment to tunnelling of a droplet was found to be the viscosity of the medium, not the energy barrier which had to be surmounted.

Acknowledgements: We thank A. Castro-Neto, B. A. Jones and H. Rieger for very helpful conversations. This work was supported by NSF-DMR-00081075 (AJM), the Ames Laboratory, operated for the U.S. Department of Energy by Iowa State University under Contract No. W-7405-Eng-82 and an award from the Research Corporation (J. S.). We acknowledge the hospitality of the Aspen Center for Physics, where part of this work was performed.


format: text

text:
# An additive variant of Tsallis generalized entropy

## I Introduction

The non-extensivity is one of the major features of Tsallis entropy and perhaps the reason for most of the criticism it receives. This property was recognized in the very first paper tsallis88  where the generalized entropy has been proposed in its discrete form

$$S_q(X)=\frac{1_{i=1}^N[p(X=x_i)]^q}{q1},$$
(1)
where $`q`$ is a real constant and $`p(X=x_i)`$ is the probability distribution of a discrete random variable $`X`$ that takes value on $`\{x_1,\mathrm{},x_N\}`$. Tsallis tsallis88  obtained the non-extensive (non-additive) relation

$$S_q(X,Y)=S_q(X)+S_q(Y)(q1)S_q(X)S_q(Y)$$
(2)
computing the entropy of the joint distribution of two independently distributed random variables $`X`$ and $`Y`$. Since then, the words “extensive” and “additive” have been loosely used as synonymous in the literature about Tsallis statistics touchette . However, following Touchette’s touchette  guidelines for the usage of these terms, the property expressed by Eq. (2) is called here “non-additivity”.

In a recent paper, Abe PhysA305\_62  has shown that, for composite systems resulting from the joint observation of two independently distributed sub-systems, the "physical entropy" is additive at the macroscopic level, even in the context of Tsallis statistics, where by “physical entropy” he means the function (R nyi entropy)

$$S^{}=\frac{\mathrm{ln}(1+(1q)S_q)}{1q}=\frac{\mathrm{ln}(_{i=1}^N[p(X=x_i)]^q)}{1q},$$
(3)
where $`S_q`$ is the Tsallis entropy.

Power-law (Tsallis) canonical distributions have been connected to the finitude of the dimension of the phase space of the heat bath in Plastino\_94 , Alm01  and Adib .

The issue of additivity within Tsallis statistics has been addressed in the scientific literature up to now assuming the joint observation of independent random variables as done in the derivation of Eq. (2). Regardless of the validity of this rule for non-interacting systems and for “weakly interacting systems” in the thermodynamic limit, the same can not be generally applied for finite thermodynamic systems. We say that two systems are in weak interaction when the Hamiltonian of the composite system is just the sum of the Hamiltonians of its sub-components, $`H_{1+2}=H_1+H_2`$. As a consequence, the internal energy $`(U=H\rho 𝑑\text{q}𝑑\text{p})`$ of such a system in an equilibrium state $`\rho `$ is an additive variable.

Consider, for instance, the observation of the energy. When one joins (through weak interaction) two independent sub-systems, $`A_1`$ and $`A_2`$, the distribution of the energy of the composite system $`A`$, ($`E_A=E_{A_1}+E_{A_2}`$), is given by the convolution of the sub-component’s distributions, $`\rho (E_A)=\rho (E_{A_1})\rho (E_{A_2})`$, which is generally different from the distribution of the joint observation of $`E_{A_1}`$ and $`E_{A_2}`$, $`\rho (E_{A_1},E_{A_2})=\rho (E_{A_1})\rho (E_{A_2})`$. In other words, when observing $`E_A`$, it is not necessary to know how much of it comes from system $`A_1`$ or $`A_2`$, but only that the sum of these contributions adds up to the observed quantity. In addition to this, even the range of the energy value of each sub-system changes upon the composition: If each separate sub-system initially has internal energy equal to $`E_0`$, after their merging the composite system will have energy equal to $`2E_0`$ and the energy of each sub-system will now fluctuate in the interval $`[0,2E_0]`$. Therefore, the joint distribution of the simultaneous observation of two independent systems $`A_1`$ and $`A_2`$ does not represent the composite thermodynamic system $`A_1+A_2`$. Since the distribution of a sum of independent random variables is a “convolution” of the distributions of the component variables, the additivity of the energy implies that the correct distribution for the composite system shall be obtained through a convolution operation. Only in the thermodynamic limit, when the sub-systems become infinite, the resulting distribution will factor into the product of the component distributions. Moreover, it can be rigorously proof that for systems in weak interaction the Boltzmann-Gibbs canonical distribution (exponential) is the only one for which the convolution of the components equals their product.

We will derive a physical entropy function $`S`$ that is additive at the equilibrium distribution even for finite systems. For a given system $`A`$ resulting from a weak combination of two sub-systems $`A_1`$ and $`A_2`$, (all of them in equilibrium), the entropy function $`S`$ satisfies the additive rule $`S(A)=S(A_1)+S(A_2)`$. We will explicitly evaluate $`S`$ for canonical ergodic systems. This analysis will be based on the canonical ensemble theory of ergodic systems as presented by Kinchin Kinchin , whose basic elements we recollect in Section II. The present work is a follow up of the study started in Alm01  where it was presented a one-to-one connection between the canonical distribution and the form of the entropy, and the Tsallis entropy was derived from first principles. This alternative approach for Tsallis statistics has been already used by Andrade et al. PRE65\_036121  to formulate a molecular dynamics scheme for Hamiltonian systems exhibiting power-law distributions, generalizing the extended Hamiltonian approach presented by Nos Nose . This setup was also used by Adib et al. Adib  to analyse the statistical behavior of systems with homogenous Hamiltonians.

## II Ergodic systems

For a generic system with Hamiltonian $`H=H(\text{q},\text{p})`$ and phase space $`\mathrm{\Xi }`$, let $`V(E)`$ denote the volume of the region of $`\mathrm{\Xi }`$ where $`H(\text{q},\text{p})E`$, $`\mathrm{\Omega }(E)=V(E)/E`$ be the structure function Kinchin , and $`\beta (E)`$ be the function

$$\beta (E)=\frac{d}{dE}\mathrm{ln}\mathrm{\Omega }(E).$$
(4)
The distribution on the surface of constant energy ($`\mathrm{\Sigma }(E)`$) that is invariant under the motion of the system is called microcanonical and has density given by (see Chapter II of Kinchin )

$$\rho ^{\text{(mc)}}(\text{q},\text{p})=\frac{1}{\mathrm{\Omega }(E)H(\text{q},\text{p})},$$
(5)
where $`H(\text{q},\text{p})`$ is the norm of the gradient of the Hamiltonian $`H(\text{q},\text{p})`$.

For a system $`A`$ which is the weak combination of two sub-systems $`A_1`$ and $`A_2`$ (denoted by $`A=A_1A_2`$), we have that $`H(\text{q},\text{p})=H_1(\text{q}_1,\text{p}_1)+H_2(\text{q}_2,\text{p}_2)`$ with no interaction term between $`A_1`$ and $`A_2`$ in the Hamiltonian, i.e., the energies of the sub-systems are additive. In that case the canonical distribution in $`\mathrm{\Xi }_1`$, when the total energy of the composite system $`A`$ is equal to $`E`$, is given by the density function Kinchin

$$\rho _1(\text{q},\text{p})=\frac{\mathrm{\Omega }_2(EH_1(\text{q},\text{p}))}{\mathrm{\Omega }(E)}.$$
(6)
The structure function of $`A`$ satisfies the convolution relation

$$\mathrm{\Omega }(E)=_0^E\mathrm{\Omega }_1(\xi )\mathrm{\Omega }_2(E\xi )𝑑\xi .$$
(7)
Observe that the microcanonical distribution on $`\mathrm{\Sigma }_1(E_1)`$ is the conditional distribution of $`\rho _1`$ obtained by fixing the value of $`H_1(\text{q},\text{p})=E_1`$, which results in

$$\rho _1^{\text{(mc)}}(\text{q},\text{p})=\frac{1}{\mathrm{\Omega }_1(E_1)H_1(\text{q},\text{p})}.$$
(8)
It has been shown in Alm01  that the form of the canonical $`\rho _1`$ depends on the value of the parameter $`q_2`$, which is related to $`\beta _2(E)`$ through the expression

$$\frac{d}{dE}\left(\frac{1}{\beta _2(E)}\right)=q_21.$$
(9)
When $`q_2=1`$, $`\rho _1`$ is an exponential (Boltzmann-Gibbs) and, when $`q_21`$, $`\rho _1`$ is a power-law (Tsallis) distribution. Identifying system $`A_2`$ with the heat bath $``$ and integrating Eq. (9) with the initial condition $`\mathrm{\Omega }_2(0)=0`$ we get that

$$\mathrm{\Omega }_2(E)=C_0E^{1/(q_21)},$$
(10)
where $`C_0`$ is a constant, and consequently

$$\beta _2(E)=\frac{1}{(q_21)E}.$$
(11)
We are going to consider that the elementary systems appearing from now on have structure functions in the form $`\mathrm{\Omega }(E)=CE^{1/(q1)}`$, with specific constants $`C`$ and $`q`$ for each one. The parameter $`q`$, which is a signature of the system’s phase space geometry, may vary from system to system. Indeed, Almeida et al. temperature  have showed that the weak composition of two systems $`𝒪`$ and $``$, with parameters $`q_𝒪`$ and $`q_{}`$, respectively, results in a system $`A`$ with parameter $`q_A`$ given by

$$\frac{q_A}{q_A1}=\frac{q_𝒪}{q_𝒪1}+\frac{q_{}}{q_{}1},$$
(12)
and that the temperature definition obeying the theorem of equipartition of energy leads to the expressions

$`kT`$ $`=`$ $`\left({\displaystyle \frac{q_A1}{q_A}}\right)E_A`$ (13)
$`=`$ $`\left({\displaystyle \frac{q_𝒪1}{q_𝒪}}\right)<H_𝒪>`$ (14)
$`=`$ $`\left({\displaystyle \frac{q_{}1}{q_{}}}\right)<H_{}>.`$ (15)
Let’s consider a weak combination of independent canonical sub-systems $`A_i`$, each one having an observable sub-system $`𝒪_i=𝒪(A_i)`$ and a heat bath $`_i=(A_i)`$. The observable and the heat bath of the composite system $`A`$ are, respectively, the weak combinations of the observables and of the heat baths of the sub-systems, i.e., $`A=_iA_i`$, $`𝒪(A)=_i𝒪(A_i)`$, and $`(A)=_i(A_i)`$. Identifying the observable part $`𝒪(=𝒪(A))`$ and the heat bath $`(=(A))`$, respectively, with the systems $`A_1`$ and $`A_2`$ in Eq. (6), we write the canonical distribution of $`A`$ as

$$\rho _𝒪(H_𝒪)=\frac{\mathrm{\Omega }_{}(E_AH_𝒪)}{\mathrm{\Omega }_A(E_A)},$$
(16)
where $`E_A`$ is the total energy of $`A`$ , $`H_𝒪`$ is the energy of its observable sub-system, and $`\mathrm{\Omega }_A`$ and $`\mathrm{\Omega }_{}`$ are the structure functions of the total system $`A`$ and of its heat bath $``$, respectively. Writing $`\rho _𝒪`$ in the usual Tsallis density form

$$\rho (H_𝒪)=\rho _𝒪(H_𝒪)=\rho _0\left[1(q_{}1)\beta H_𝒪\right]^{1/(q_{}1)},$$
(17)
we have that

$$\rho _0=\frac{\mathrm{\Omega }_{}(E_A)}{\mathrm{\Omega }_A(E_A)}$$
(18)
and

$$\beta =\beta _{}(E_A)=\frac{1}{(q_{}1)E_A}.$$
(19)
Moreover, the equivalent to Eq. (12) reads as

$$\frac{q_A}{q_A1}=\frac{q_𝒪}{q_𝒪1}+\frac{q_{}}{q_{}1},$$
(20)
where $`q_A`$, $`q_𝒪`$ and $`q_{}`$ are the parameters $`q`$ of $`A`$, $`𝒪`$ and $``$, respectively.

Applying Eq. (12) recursively to $`A`$, $`𝒪`$ and $``$ we get the relations

$$\frac{q_A}{q_A1}=\underset{i}{}\frac{q_{A_i}}{q_{A_i}1}$$
(21)
$$\frac{q_𝒪}{q_𝒪1}=\underset{i}{}\frac{q_{𝒪_i}}{q_{𝒪_i}1}$$
(22)
and

$$\frac{q_{}}{q_{}1}=\underset{i}{}\frac{q__i}{q__i1}.$$
(23)
## III Physical entropy

For a system with Hamiltonian $`H`$ and probability density function $`\rho `$ let $`U`$ be the internal energy, $`W`$ be the work of the external forces and $`Q`$ be the heat supplied to the system. In terms of the microscopic functions we have that

$$U=_\mathrm{\Xi }H\rho 𝑑\text{q}𝑑\text{p}.$$
(24)
The first law of thermodynamics states that

$$\delta U=\delta W+\delta Q.$$
(25)
Computing the variation $`\delta U`$ from its microscopic definition, Eq. (24), we get that

$$\delta U=_\mathrm{\Xi }(\delta H)\rho 𝑑\text{q}𝑑\text{p}+_\mathrm{\Xi }H(\delta \rho )𝑑\text{q}𝑑\text{p}.$$
(26)
The first term on the right hand side of this equation is the variation of the work $`\delta W`$ and so we get that

$$\delta Q=_\mathrm{\Xi }H(\delta \rho )𝑑\text{q}𝑑\text{p}.$$
(27)
The second law of thermodynamics states that there are an exact differential functions $`S`$, called entropies, such that

$$\delta S=\frac{1}{T}\delta Q,$$
(28)
where $`T`$, the integrating factor, is a temperature scale.

As pointed out in a previous work Alm01 , the functional form of the physical entropy is determined by the canonical distribution: Given a canonical density function in $`\mathrm{\Xi }`$, $`\rho =f(H)`$, with a non-negative monotonic function $`f(x)`$ defined for $`x0`$, the functions (entropies) that satisfy the second law of thermodynamics are given by

$$S(\rho )=_\mathrm{\Xi }F(\rho )𝑑\text{q}𝑑\text{p}$$
(29)
with

$$F(\rho )=\frac{1}{T}_0^\rho f^1(\xi )𝑑\xi +\psi \rho ,$$
(30)
where $`\psi `$ is an arbitrary constant. As will be seen later, it is convenient to choose $`\psi =k`$, where $`k`$ is the Boltzmann constant.

Applying the above procedure to the Tsallis density function

$$\rho (H)=\rho _0[1(q_{}1)\beta H]^{1/(q_{}1)},$$
(31)
where $`(q_{}1)\beta =1/E`$ is the inverse of the total energy of the system, considering that

$$\rho _0=\rho (0)=\left\{_\mathrm{\Xi }[1(q_{}1)\beta H]^{1/(q_{}1)}𝑑\text{q}𝑑\text{p}\right\}^1,$$
(32)
is constant, we get the function

$$(\rho )=\frac{1}{T(q_{}1)\beta }\left[1\frac{\rho _0}{q_{}}_\mathrm{\Xi }\left(\frac{\rho }{\rho _0}\right)^q_{}𝑑\text{q}𝑑\text{p}\right]k,$$
(33)
which is not the entropy yet, because of the imposed condition that $`\rho _0`$ be constant. Collecting and rearranging terms we get the alternative expression

$`(\rho )=`$ $`{\displaystyle \frac{1}{q_{}T\beta }}\left\{{\displaystyle \frac{\rho _0}{q_{}1}}{\displaystyle _\mathrm{\Xi }}\left[\left({\displaystyle \frac{\rho }{\rho _0}}\right)\left({\displaystyle \frac{\rho }{\rho _0}}\right)^q_{}\right]𝑑\text{q}𝑑\text{p}\right\}`$ (34)
$`+\left({\displaystyle \frac{1}{q_{}T\beta }}k\right).`$
This expression may be written also as

$$(\rho )=\left(\frac{1}{q_{}T\beta }\right)\rho _0S_{\text{TS}}\left(\rho /\rho _0\right)+\left(\frac{1}{q_{}T\beta }k\right),$$
(35)
where

$$S_{\text{TS}}(\rho /\rho _0)=\frac{1}{q_{}1}_\mathrm{\Xi }\left[\left(\frac{\rho }{\rho _0}\right)\left(\frac{\rho }{\rho _0}\right)^q_{}\right]𝑑\text{q}𝑑\text{p}$$
(36)
is the Tsallis entropy evaluated at the ratio $`\rho /\rho _0`$. Observe in Eq. (35) the splitting of the roles played by the density function form, $`\rho /\rho _0`$, and by the normalization value $`\rho _0`$.

Now, considering that both $`\rho `$ and $`\rho _0`$ vary along the process, the evaluation of $`\delta `$ will produce the sought expression for the entropy. From Eq. (33) we have that

$`\delta =`$ $`{\displaystyle \frac{1}{T(q_{}1)\beta }}[{\displaystyle _\mathrm{\Xi }}\left({\displaystyle \frac{\rho }{\rho _0}}\right)^{q_{}1}(\delta \rho )d\text{q}d\text{p}+`$ (37)
$`(\delta \rho _0){\displaystyle \frac{q_{}1}{q_{}}}{\displaystyle _\mathrm{\Xi }}\left({\displaystyle \frac{\rho }{\rho _0}}\right)^q_{}d\text{q}d\text{p}].`$
Substituting the expressions for $`\rho `$ and $`\rho _0`$, and considering the temperature expressions Eqs. (13)-(15) we get that

$$\delta =\frac{1}{T}_\mathrm{\Xi }H(\delta \rho )𝑑\text{q}𝑑\text{p}+k\delta (\mathrm{ln}\rho _0)=\frac{1}{T}\delta Q+k\delta (\mathrm{ln}\rho _0),$$
(38)
and from it we conclude that

$$\delta S=\frac{1}{T}\delta Q=\delta k\delta (\mathrm{ln}\rho _0),$$
(39)
and so the generic entropies are

$$S(\rho )=(\rho )k\mathrm{ln}\rho _0+C,$$
(40)
where $`C`$ is an arbitrary constant.

For a fixed value of $`D=q_{}T\beta `$, the maximization of $`S(\rho )`$ under the conditions

$`{\displaystyle \rho 𝑑\text{q}𝑑\text{p}}=1,`$ (41)
$`{\displaystyle H\rho 𝑑\text{q}𝑑\text{p}}=U,`$ (42)
and

$$\rho _0=\rho (0),$$
(43)
results in the Tsallis density function $`\rho (H)`$ given by Eqs. (31) and (32).

### III.1 Additivity of $``$

Let’s compute the explicit value of $`_A(\rho )`$ for a system $`A(=𝒪)`$ whose canonical distribution on $`\mathrm{\Xi }_𝒪`$, $`\rho (H_𝒪)`$, is given in Eq. (17). From a direct integration we get

$$\rho _0_{\mathrm{\Xi }_𝒪}\left(\frac{\rho }{\rho _0}\right)^q_{}𝑑\text{q}𝑑\text{p}=\left[1\frac{<H_𝒪>}{E_A}\right],$$
(44)
where we adopt the notation $`<f>=_{\mathrm{\Xi }_𝒪}f\rho 𝑑\text{q}𝑑\text{p}`$ for the average of a generic function $`f`$.

From Eqs. (19) and (13) we get that

$$\frac{1}{T(q_{}1)\beta }=k\left(\frac{q_A}{q_A1}\right).$$
(45)
Substituting Eqs. (44) and (45) in Eq. (33) and using Eq. (20) one gets that

$$_A(\rho )=k\left(\frac{q_𝒪}{q_𝒪1}\right).$$
(46)
Using the combination rule for the $`q_𝒪`$’s, Eq. (22), we observe that under weak combination of sub-systems the function $`_A`$ is additive. In general, if $`A`$ is a (weak) combination of $`N`$ sub-systems $`A_i`$, $`i=1,\mathrm{},N`$, $`(A=A_1A_2\mathrm{}A_N)`$, all at the same temperature $`T`$, then

$$_A(\rho )=\underset{i=1}{\overset{N}{}}_{A_i}(\rho _i).$$
(47)
Notice that $`_A(\rho )`$ does not depend on the parameter $`q_{}`$, hence Eq. (46) is valid in the limit case $`\text{q}_{}1`$ (Boltzmann-Gibbs), which can be also verified through a direct evaluation of Eq. (56) with $`\rho =\rho _0\mathrm{exp}(\beta H)`$. Observe that Eq. (14) together with Eq. (46) leads to the relation between $`_A`$, the internal energy $`U=<H_𝒪>`$ and the temperature $`T`$,

$$_A(\rho )=\frac{U}{T}.$$
(48)
Substituting the alternative forms of $`(\rho )`$ into Eq. (40) we obtain

$$S(\rho )=k\left(\frac{q_𝒪}{q_𝒪1}\right)k\mathrm{ln}\rho _0+C$$
(49)
and

$$S(\rho )=\frac{U}{T}k\mathrm{ln}\rho _0+C.$$
(50)
This last expression was also obtained through a different route by Potiguar and Costa fabricio02a .

Notice that the constants $`q_{}`$ and $`\beta `$ appearing in Eq. (34) are associated with the heat bath and that

$`{\displaystyle \frac{1}{q_{}T\beta }}`$ $`=`$ $`k\left({\displaystyle \frac{q_A}{q_A1}}\right)\left({\displaystyle \frac{q_{}1}{q_{}}}\right)`$ (51)
$`=`$ $`k\left({\displaystyle \frac{q_𝒪}{q_𝒪1}}\right)\left({\displaystyle \frac{q_{}1}{q_{}}}\right)+k`$ (52)
is determined by the constants $`q_𝒪`$ and $`q_{}`$. From this last equation, we see that the commonly used relation $`kT\beta =1`$ is valid only in the limit $`q_{}1`$, i.e., in the Boltzmann-Gibbs statistics. This fact has been overseen by some authors who inadvertently adopted $`kT=1/\beta `$ within Tsallis statistics, which might have led to misleading results.

### III.2 Microcanonical ensemble

The microcanonical ensemble is represented by a constant energy $`H_𝒪`$ and its density function in this surface $`\mathrm{\Sigma }_𝒪`$ is

$$\rho ^{\text{(mc)}}(H_𝒪)=\frac{1}{\mathrm{\Omega }_𝒪(H_𝒪)H_𝒪}.$$
(53)
Hence, we may consider that in Eq. (53) $`\rho _0=(\mathrm{\Omega }_𝒪)^1`$ and obtain from Eq. (49) that

$$S(\rho ^{\text{(mc)}})=k\mathrm{ln}\left[\mathrm{\Omega }_𝒪(H_𝒪)\right]+\left[k\left(\frac{q_𝒪}{q_𝒪1}\right)+C\right],$$
(54)
which is the Boltzmann entropy for the microcanonical ensemble plus constants.

### III.3 Boltzmann-Gibbs canonical ensemble

The Boltzmann-Gibbs canonical ensemble considers an infinite heat bath, which is represented by the limit condition $`\text{q}_{}1`$ keeping $`\beta `$ fixed. Taking this limit in Eq. (17) and in either of the Eqs. (33) or (34) we arrive at

$$\rho (H_𝒪)=\rho _0\mathrm{exp}(\beta H_𝒪)$$
(55)
and

$$(\rho )=\frac{1}{T\beta }\left\{1_\mathrm{\Xi }\rho \mathrm{ln}(\rho )𝑑\text{q}𝑑\text{p}+\mathrm{ln}(\rho _0)\right\}k,$$
(56)
from which, considering also that in the limit $`q_{}1`$ we have $`1/(T\beta )=k`$ (see Eq. (52) below), we get that

$$S=(\rho )k\mathrm{ln}(\rho _0)+C=k_\mathrm{\Xi }\rho \mathrm{ln}(\rho )𝑑\text{q}𝑑\text{p}+C$$
(57)
is the Boltzmann-Gibbs-Shannon entropy plus an arbitrary constant $`C`$.

### III.4 Additivity of $`k\mathrm{ln}\rho _0`$

Let’s consider that the system $`A`$ is the weak combination of $`n`$ identical sub-systems $`A_i`$, $`i=1,\mathrm{},n`$, which by their turn are the weak combination of an observable $`𝒪_i`$ and a heat bath $`_i`$. Let $`q`$, $`q_1`$ and $`q_2`$ be the $`q`$-parameters of the $`A_i`$, $`𝒪_i`$ and $`_i`$, respectively. Therefore, we have that the parameters for the system $`A`$, its observable part $`𝒪`$, and its heat bath $``$ are respectively $`q_A`$, $`q_𝒪`$ and $`q_{}`$, and they satisfy the relations of Eqs. (21)-(23), i.e.,

$$\frac{q_A}{q_A1}=\underset{i}{}\frac{q_{A_i}}{q_{A_i}1}=n\frac{q}{q1}$$
(58)
$$\frac{q_𝒪}{q_𝒪1}=\underset{i}{}\frac{q_{𝒪_i}}{q_{𝒪_i}1}=n\frac{q_1}{q_11}$$
(59)
and

$$\frac{q_{}}{q_{}1}=\underset{i}{}\frac{q__i}{q__i1}=n\frac{q_2}{q_21}.$$
(60)
Computing the structure function of the system $`A`$ using Eq. (7) we obtain

$$\mathrm{\Omega }_A(E_A)=\mathrm{\Omega }_𝒪(E_A)\mathrm{\Omega }_{}(E_A)E_A\frac{\mathrm{\Gamma }\left(\frac{q_𝒪}{q_𝒪1}\right)\mathrm{\Gamma }\left(\frac{q_{}}{q_{}1}\right)}{\mathrm{\Gamma }\left(\frac{q_A}{q_A1}\right)}$$
(61)
Therefore,

$$\rho _0=\frac{\mathrm{\Omega }_{}(E_A)}{\mathrm{\Omega }_A(E_A)}=\frac{\mathrm{\Gamma }\left(\frac{q_A}{q_A1}\right)}{\mathrm{\Omega }_𝒪(E_A)E_A\mathrm{\Gamma }\left(\frac{q_𝒪}{q_𝒪1}\right)\mathrm{\Gamma }\left(\frac{q_{}}{q_{}1}\right)}$$
(62)
Recursively applying Eq. (7) to compute $`\mathrm{\Omega }_𝒪(E_A)`$ we get that

$$\mathrm{\Omega }_𝒪(E_A)E_A\mathrm{\Gamma }\left(\frac{q_𝒪}{q_𝒪1}\right)=\left[C_0\mathrm{\Gamma }\left(\frac{q_1}{q_11}\right)E_A^{\frac{q_1}{q_11}}\right]^n,$$
(63)
where $`C_0`$ is the multiplicative constant appearing in the structure functions $`\mathrm{\Omega }_{𝒪_i}`$. Plugging (63) into (62) we obtain

$$\rho _0=\frac{\mathrm{\Gamma }\left(\frac{q_A}{q_A1}\right)}{\mathrm{\Gamma }\left(\frac{q_{}}{q_{}1}\right)}\left[C_0\mathrm{\Gamma }\left(\frac{q_1}{q_11}\right)E_A^{\frac{q_1}{q_11}}\right]^n.$$
(64)
\>From Eq. (13) the total energy $`E_A`$ in terms of the temperature and $`q`$ is

$$E_A=kT\left(\frac{q_A}{q_A1}\right)=nkT\left(\frac{q}{q1}\right),$$
(65)
\>From these two last equations we get that

$`k\mathrm{ln}\rho _0`$ $`=`$ $`k\mathrm{ln}\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{nq_2}{q_21}\right)}{\mathrm{\Gamma }\left(\frac{nq}{q1}\right)}}\right]`$ (69)
$`+kn\left({\displaystyle \frac{q_1}{q_11}}\right)\mathrm{ln}[n]`$
$`+kn\left({\displaystyle \frac{q_1}{q_11}}\right)\mathrm{ln}\left[kT\left({\displaystyle \frac{q}{q1}}\right)\right]`$
$`+kn\mathrm{ln}\left[C_0\mathrm{\Gamma }\left({\displaystyle \frac{q_1}{q_11}}\right)\right].`$
The terms in Eqs. (69) and (69) are additive since they are linear in $`n`$, but the terms in (69) and (69) are non-additive. We can eliminate this non-additivity by taking the constant $`C`$ appearing on Eq. (40) equal to the negative of the sum of all terms on the right hand side of Eqs. (69)- (69) that depend, besides of $`n`$, only on the parameters related to the smallest system sub-components, viz., $`q`$, $`q_1`$, $`q_2`$ and $`C_0`$. This adequate choice of the constant $`C`$ is analogous to the solution presented by Gibbs to the paradox bearing his name. In this way, we are left solely with the only term depending on $`T`$,

$$k\mathrm{ln}\rho _0+C=kn\left(\frac{q_1}{q_11}\right)\mathrm{ln}[T]=k\left(\frac{q_𝒪}{q_𝒪1}\right)\mathrm{ln}[T].$$
(70)
Combining Eqs. (40), (46) and (70) the resulting entropy function becomes

$$S=k\left(\frac{q_𝒪}{q_𝒪1}\right)\left[1+\mathrm{ln}[T]\right]=\frac{U}{T}\left[1+\mathrm{ln}[T]\right].$$
(71)
Since the temperature $`T`$ is an intensive variable and the internal energy is additive, this entropy form is additive.

## IV Conclusions

We have seen that both the microcanonical and the classical canonical (Boltzmann-Gibbs) distributions can be obtained from the Tsallis distribution. We have obtained also an expression for the “physical entropy” of weak-interaction systems with finite heat baths ($`q_{}1`$) and showed that, for systems with constant energy, it yields the classical microcanonical Boltzmann entropy and that when the heat bath gets infinite ($`q_H1`$), it becomes the classical Boltzmann-Gibbs-Shannon entropy. Moreover, we saw that this “physical entropy” may be written as a variant of Tsallis entropy plus the product of the Boltzmann constant and the logarithm of the partition function. In addition to this, we saw that it may be turned into an additive variable by adding appropriate constants in a procedure reminiscent to what was done by Gibbs to overcome the Gibbs paradox of the microcanonical entropy for ideal gases.

Our results strength the physical meaning of the generalized entropy by showing that the generalized thermostatistics fills the gap between the microcanonical and the Boltzmann-Gibbs canonical ensembles. These results are summarized in the following frame:

| Ensemble | microcanonical | Tsallis | BG canonical |
| --- | --- | --- | --- |
| heat bath |  | finite | infinite |
| constraint | $`H_𝒪=\text{constant}`$ | $`q_{}1`$ | $`q_{}1`$ |
| density $`\rho `$ | $`\rho _0[H]^1`$ | $`\rho _0[1(q_{}1)\beta H_𝒪]^{1/(q_{}1)}`$ | $`\rho _0\mathrm{exp}(\beta H_𝒪)`$ |
| Entropy | $`S=k\left(\frac{q_𝒪}{q_𝒪1}\right)k\mathrm{ln}[\rho _0]`$ | |  |

The analysis presented here aimed at the additivity of the entropy within thermodynamic systems and it was not intended to contest the validity of Eq. (2), which is applicable for the original Tsallis entropy $`S_{\text{TS}}`$ under joint observation of independent random variables. It was also seen that when weakly combining systems the parameter $`q`$ varies with the dimensionality of the system’s phase space following the additive rule for weak interactions (Eq. (12)).

## V Acknowledgments

This work was partially supported by CNPq (Brazil).


text:
# High energy photons from passage of jets through quark gluon plasma

## Abstract

We calculate the production of high energy photons from Compton scattering and annihilation of a quark jet passing through a quark gluon plasma produced in a relativistic heavy ion collision. The contributions are large and reflect the momentum distribution of the jets and the initial conditions of the plasma.

preprint: DUKE-TH-02-224

Relativistic heavy ion collisions are studied with the aim of producing a plasma of quarks and gluons (QGP). Photons are considered to be an important probe for the investigation of the formation and evolution of such a plasma due to their weak final-state interactions Fei:76 . Once produced they carry the information about the conditions of the environment in which they were created, encoded in their momentum distribution, thus providing a glimpse deep into the bulk of strongly interacting matter. Though most of the measured photons have their origin in the decay of hadrons after the QGP phase, it has become possible to isolate the direct photons produced in such collisions WA98:00 .

The sources of direct photons considered so far include quark annihilation, Compton scattering, and bremsstrahlung following the initial hard scattering of partons of the nuclei Owens:87 , as well as thermal photons from the QGP KapLiSei:91 ; BNNR:92 ; RPSSS:96 ; AGKZ:98 ; ArnMoYaf:01  and from hadronic interactions in the hot hadronic gas after the hadronization of the plasma KapLiSei:91 ; XioShuBro:92 . The pre-equilibrium production of photons has also been investigated by several authors Strickl:94 . Results are available for production from the entire history of the system SriSin:01 .

In this letter we study a new source of direct photons originating from the passage of the produced high energy quark jets through the QGP (jet-photon conversion). A fast quark passing through the plasma will produce photons by Compton scattering with the thermal gluons and annihilation with the thermal antiquarks. This process is higher order in $`\alpha _s`$ compared with photons from initial hard scatterings, but it is not a subleading contribution, since it corresponds to double scattering, which is enhanced by the size of the system. For cold nuclear matter this effect is encoded in multi-parton matrix elements which are enhanced by powers of $`A^{1/3}`$ LQS:94 . Below, we find that this source is at least comparable in strength to the other direct photon sources and even dominates in the range $`p_{}6`$ GeV for Au+Au collisions at the Relativistic Heavy Ion Collider (RHIC).

We also demonstrate that the $`p_{}`$-distribution of these photons is directly proportional to the momentum distribution of jets at an early stage after their production, before they have lost energy on their travel through the plasma. Since the measured high-$`p_{}`$ hadron spectrum is proportional to the spectrum of partons after they have left the plasma, a comparision of both spectra could provide a quantitative determination of the energy loss and help confirm the mechanism of jet quenching GyulWang:94 ; BDMPS:97 .

Furthermore, the photon yield depends on the integrated density of the matter traversed by the jets and thus can provide a measurement of this quantity. We emphasize that our mechanism is distinct from the back-to-back correlation among direct photons and leading hadrons, which was proposed as a direct measurement of the jet energy loss in dense matter WangHS:96 .

The idea of jet-photon conversion is based on the properties of the annihilation and Compton cross sections. The kinematics of the annihilation of a quark antiquark pair $`(q+\overline{q}\gamma +g)`$ is expressed in terms of the Mandelstam variables $`s=(p_q+p_{\overline{q}})^2`$, $`t=(p_qp_\gamma )^2`$ and $`u=(p_{\overline{q}}p_\gamma )^2`$, where the $`p_i`$ are the four-momenta of the particles. The differential cross section at Born level for massless partons has the form $`d\sigma /dt=8\pi \alpha \alpha _se_q^2\left(u/t+t/u\right)/9s^2`$, where $`e_q`$ is the fractional charge of the annihilating quarks. This implies that the largest contribution to the production of photons arises from small values of $`t`$ or $`u`$, corresponding to $`𝐩_\gamma 𝐩_q`$ or $`𝐩_\gamma 𝐩_{\overline{q}}`$. The process can be visualized as a conversion of one of the annihilating quarks to a photon. This picture immediately suggests that the photons from this process provide a direct measurement of the quark (or antiquark) momentum. For the Compton process $`(g+q\gamma +q)`$ we have $`d\sigma /dt=\pi \alpha \alpha _se_q^2\left(u/s+s/u\right)/3s^2`$ . The dominant contribution now comes from the region of small $`u`$, when $`𝐩_\gamma 𝐩_q`$, corresponding to the conversion of the quark (or the antiquark) into a photon.

Following Wong:94 ; Nadeau:93  we approximate the invariant photon differential cross sections for the annihilation and Compton processes as

$$E_\gamma \frac{d\sigma ^{(\mathrm{a})}}{d^3p_\gamma }\sigma ^{(\mathrm{a})}(s)\frac{1}{2}E_\gamma \left[\delta (𝐩_\gamma 𝐩_𝐪)+\delta (𝐩_\gamma 𝐩_{\overline{𝐪}})\right]$$
(1)
$$E_\gamma \frac{d\sigma ^{(\mathrm{C})}}{d^3p_\gamma }\sigma ^{(\mathrm{C})}(s)E_\gamma \delta (𝐩_\gamma 𝐩_𝐪).$$
(2)
Here $`\sigma ^{(\mathrm{a})}(s)`$ and $`\sigma ^{(\mathrm{C})}(s)`$ are the corresponding total cross sections. Now consider an ensemble of fast quarks passing through a hot medium. A certain fraction, determined by the total cross section, will undergo these annihilation or Compton processes initiated by the medium and will emit photons. Whenever this happens the energy of the emitted photon will reflect the initial energy of the converted quark.

Using (1,2) the rate of production of photons due to annihilation and Compton scattering are given by Wong:94

$`E_\gamma {\displaystyle \frac{dN^{(\mathrm{a})}}{d^4xd^3p_\gamma }}`$ $`=`$ $`{\displaystyle \frac{16E_\gamma }{2(2\pi )^6}}{\displaystyle \underset{q=1}{\overset{N_f}{}}}f_q(𝐩_\gamma ){\displaystyle }d^3pf_{\overline{q}}(𝐩)[1+f_g(𝐩)]\sigma ^{(\mathrm{a})}(s){\displaystyle \frac{\sqrt{s(s4m^2)}}{2E_\gamma E}}+(q\overline{q}),`$ (3)
$`E_\gamma {\displaystyle \frac{dN^{(\mathrm{C})}}{d^4xd^3p_\gamma }}`$ $`=`$ $`{\displaystyle \frac{16E_\gamma }{(2\pi )^6}}{\displaystyle \underset{q=1}{\overset{N_f}{}}}f_q(𝐩_\gamma ){\displaystyle d^3pf_g(𝐩)\left[1f_q(𝐩)\right]\sigma ^{(\mathrm{C})}(s)\frac{(sm^2)}{2EE_\gamma }}+(q\overline{q}).`$ (4)
The $`f_i`$ are distribution functions for the quarks, antiquarks, and gluons. Inserting thermal distributions for the gluons and quarks one can obtain an analytical expression for these emission rates for an equilibrated medium KapLiSei:91 ; BNNR:92 ; Wong:94 .

We propose that the phase-space distribution of the quarks and gluons produced in a nuclear collision can be approximately decomposed into two components, a thermal component $`f_{\mathrm{th}}`$ characterized by a temperature $`T`$ and a hard component $`f_{\mathrm{jet}}`$ given by hard scattering of partons and limited to transverse momenta $`p_{}1`$ GeV/$`c`$: $`f(𝐩)=f_{\mathrm{th}}(𝐩)+f_{\mathrm{jet}}(𝐩)`$. We note that $`f_{\mathrm{jet}}`$ dominates for large momenta, while at small momenta $`f`$ is completely given by the thermal part.

The phase-space distribution for the quark jets propagating through the QGP is given by the perturbative QCD result for the jet yield LinGyu:95 :

$$\begin{array}{c}f_{\mathrm{jet}}(𝐩)=\frac{1}{g_q}\frac{(2\pi )^3}{\pi R_{}^2\tau p_{}}\frac{dN_{\mathrm{jet}}}{d^2p_{}dy}R(r)\hfill \\ \hfill \times \delta (\eta y)\mathrm{\Theta }(\tau \tau _i)\mathrm{\Theta }(\tau _{\mathrm{max}}\tau )\mathrm{\Theta }(R_{}r).\end{array}$$
(5)
where $`g_q=2\times 3`$ is the spin and colour degeneracy of the quarks, $`R_{}`$ is the transverse dimension of the system, $`\tau _i1/p_{}`$ is the formation time for the jet and $`\eta `$ is the space time rapidity. $`R(r)`$ is a transverse profile function. We shall take $`\tau _{\mathrm{max}}`$ as the smaller of the life-time $`\tau _f`$ of the QGP and the time $`\tau _d`$ taken by the jet produced at position $`𝐫`$ to reach the surface of the plasma. The boost invariant correlation between the rapidity $`y`$ and $`\eta `$ is assumed here Bj:83 . Similar results are obtained when a Gaussian correlation between $`y`$ and $`\eta `$, characterized by a width $`\mathrm{\Delta }_p2/p_{}\mathrm{cosh}y`$ is assumed LinGyu:95 .

We calculate the jet production and the direct production of photons in lowest order pQCD Owens:87 . We have used the CTEQ5L parton distributions cteq5  and EKS98 nuclear modifications EKS98  with the scale set to $`p_{}`$. We further use a $`K`$-factor of 2.5 in the jet production to account for higher order corrections, cf. EskTuo:00 . We neglect the weak dependence of $`K`$ on $`p_{}`$ and $`\sqrt{S}`$ in this study BFLPZ:00 . Our results for the number of quarks, antiquarks and gluons at $`y=0`$ can be parameterized as

$$\frac{dN^{\mathrm{jet}}}{d^2p_{}dy}|_{y=0}=T_{AA}\frac{d\sigma ^{\mathrm{jet}}}{d^2p_{}dy}|_{y=0}=K\frac{a}{(1+p_{}/b)^c}.$$
(6)
Here $`T_{AA}=9A^2/8\pi R_{}^2`$ is the nuclear thickness for a head-on collision. Numerical values for the parameters $`a`$, $`b`$ and $`c`$ are listed in Table 1.

The integral over $`𝐩`$ in Eqs. (3) and (4) is dominated by small momenta. We can therefore drop the jet part in the decomposition of the distributions $`f(𝐩)`$ in the integrands and approximate them by the thermal part. After performing the integrals the results for annihilation and Compton scattering are identical and read Wong:94

$$\begin{array}{c}E_\gamma \frac{dN_\gamma ^{(\mathrm{a})}}{d^3p_\gamma d^4x}=E_\gamma \frac{dN_\gamma ^{(\mathrm{C})}}{d^3p_\gamma d^4x}=\frac{\alpha \alpha _s}{8\pi ^2}\underset{f=1}{\overset{N_f}{}}\left(\frac{e_{q_f}}{e}\right)^2\hfill \\ \hfill \times \left[f_q(p_\gamma )+f_{\overline{q}}(p_\gamma )\right]T^2\left[\mathrm{ln}\left\{\frac{4E_\gamma T}{m^2}\right\}+C\right]\end{array}$$
(7)
with $`C=1.916`$. We have included the three lightest quark flavours so that $`_fe_{q_f}^2/e^2=2/3`$.

Once again we apply the decomposition of the distributions to (7). Inserting the thermal part we recover the known rate of emission of high energy photons from interactions within the plasma. Since we are dealing with (nearly) massless partons we have to address the infrared divergence in (7) in the limit $`m0`$. It has been shown KapLiSei:91 ; BNNR:92  that the divergence can be eliminated by including higher-order effects from the interaction of the quarks in the thermal medium with the result that $`m^2`$ in the argument of the logarithm must be replaced with $`2m_{\mathrm{th}}^2=g^2T^2/3`$. We assume that these considerations remain valid for the emission of photons even when one of the partons is from a jet while the other one is from a thermal medium. This assumption remains to be verified. We add that the thermal photon production from the plasma has recently been calculated up to two loops AGKZ:98  and to complete leading order ArnMoYaf:01 , and is about a factor of two larger than the lowest order results given here. It will be of considerable interest to extend our present work to account for these effects.

Now the rate of high energy photon emission from the passage of quark jets through the QGP is obtained by substituting the jet contributions for $`f_q`$ and $`f_{\overline{q}}`$ in (7). In order to calculate the photon yields we specify our model for the medium. Here we are interested in the emission of photons from the QGP phase alone. We assume that a thermally and chemically equilibrated plasma is produced in the collision at time $`\tau _0`$ with temperature $`T_0`$, and we ignore the transverse expansion of the plasma. Further assuming an isentropic longitudinal expansion Bj:83 , we can relate $`T_0^3\tau _0`$ to the observed particle rapidity density $`dN/dy`$. We take $`dN/dy1260`$, based on the charged particle pseudorapidity density measured by the PHOBOS experiment Phobos:02  for central collisions of Au nuclei at $`\sqrt{S_{NN}}=`$ 200 GeV. For central collision of Pb nuclei at LHC energies we use $`dN/dy=5625`$ estimated by KaMcLSri:92 . Imposing a rapid thermalization limited by $`\tau _01/3T_0`$ KaMcLSri:92  we fix our initial conditions to be $`T_0=446`$ MeV and $`\tau _0=0.147`$ fm/$`c`$ for RHIC and $`T_0=897`$ MeV and $`\tau _0=0.073`$ fm/$`c`$ for LHC.

We model the nuclei as uniform spheres and assume a transverse profile for the initial temperature as $`T(r)=T_0[2(1r^2/R_{}^2)]^{1/4}`$ where $`R_{}=1.2A^{1/3}`$ fm. We use the same profile $`R(r)=2(1r^2/R_{}^2)`$ for the jet production in (5) while performing the space-time integration $`d^4x=\tau d\tau rdrd\eta d\varphi `$. The limits of the $`\tau `$-integration are $`[\tau _0,\tau _f]`$ where $`\tau _f`$ is fixed through the relation $`T^3\tau =`$ const, when the temperature reaches $`T_f=160`$ MeV.

We also give results for the competing photon processes, the direct production in primary hard annihilation or Compton processes and the production via bremsstrahlung from a produced jet parton. For both cases the leading order results can again be found in Owens:87 . No $`K`$-factor is used for direct photon production. Our results here obey the scaling for single photons suggested in Sri:01  from empirical considerations.

In Fig. 1 we plot the results for thermal photons, direct photons due to primary processes, bremsstrahlung photons and the photons coming from jets passing though the QGP in central collision of gold nuclei at RHIC energies. The corresponding results for LHC energies are shown in Fig. 2. We see that the quark jets passing through the QGP give rise to a large yield of high energy photons. This contribution should be absent in $`pp`$ collisions. For RHIC this contribution is the dominant source of photons up to $`p_{}6`$ GeV. Note that the jet-to-photon conversion falls more rapidly with $`p_{}`$ than the direct photon yiield, similar to a higher twist correction.

A suppression of the bremsstrahlung contribution due to the multiple scattering suffered by the fragmenting partons ina  will further enhance the importance of the jet-photon conversion process. The effects of a possible chemical non-equilibrium may be estimated as follows. A reduced quark abundance will enter at least linearly into the thermal rate but can reduce the jet-photon conversion by at most a factor of two, as the Compton contribution will remain unaffected. Therefore the thermal rate will be more strongly suppressed than the jet-photon conversion. Similarly, the jet-photon conversion is only moderately sensitive to the initial conditions and the transit time of the jet in the medium. We have found that increasing $`\tau _0`$ to 1 fm/$`c`$, keeping $`T_0^3\tau _0`$ fixed, reduces our results for RHIC only by about 40 percent.

One may ask how our results will be modified by the energy loss experienced by the quark jets passing through the QGP. This effect will be relatively small for two reasons. First, we are concerned here with quark jets which lose energy at less than half the rate as gluons, which form the dominant fraction of jet events. Second, as the energy loss scales as the square of the distance traveled through the medium BDMPS:97 , the average energy loss at the moment of jet-photon conversion is only one-third of the total normal energy loss. This argument implies that, while the high-$`p_{}`$ hadron spectra encode information about the energy loss in the plasma, the photons from medium induced emission, on average, carry information about jets at an earlier stage. Obviously, the comparison of the high-$`p_{}`$ hadron and direct photon spectra can shed light on the evolution of the jet spectrum. We believe that this could be a useful method for the analysis of jet quenching and medium effects on hard particles in the QGP. More details of the calculation together with a careful analysis of the assumptions made here will be presented in a forthcoming paper.

In summary, we have discussed a novel mechanism for the production of high energy photons in relativistic heavy ion collisions. The photons arise from the passage of quark jets through the medium, when a quark is converted into a photon due to Compton scattering from a gluon or annihilation by an antiquark. Our calculations indicate that this mechanism is the leading source of directly produced photons at RHIC in the region $`p_{}6`$ GeV/$`c`$. The spectrum of the resulting photons reflects the initial distribution of the hard scattered quarks and is sensitive to the initial conditions and the traversal time of the jet in the plasma.

###### Acknowledgements.

This work was supported in part by DOE grants DE-FG02-96ER40945 and DE-AC02-98CH10886 and the Natural Sciences and Engineering Research Council of Canada. RJF is supported by the Feodor Lynen program of the Alexander von Humboldt Foundation.