text:
# A critique of q-entropy for thermal statistics

## Introduction

In 1988 a relation between entropy and probability for thermal statistics was proposed by C. Tsallis tsallis1  which is non-additive and depends on a parameter $`q`$ that is presumably determined by the nature of the thermodynamic system under consideration. For the special case $`q=1`$ this relation reduces to the Boltzmann-Gibbs entropy, but for other values of $`q`$ it is claimed that this q-entropy leads to an alternative formalism that is “entirely consistent” with the laws of thermodynamics tsallis1 ; curado1 ; tsallis2 ; tsallis3 ; baranger1 ; cohen1 . It will be shown, however, that this claim is not valid, because the total entropy of thermodynamic systems that are in thermal contact is not even defined in this formalism for the case where systems have q-entropies with different values of $`q`$. As a consequence, the fundamental concepts of temperature and heat exchange between such systems can not be introduced in this formulation of thermodynamics unless $`q`$ is a universal constant. But even in this case, maximizing the q-entropy of the combined systems leads to unphysical results unless q=1. In particular, for $`q1`$ the resulting joint probability for states of weakly coupled systems is not the product of the individual probabilities, and as a consequence the total energy of the combined systems is not the sum of the mean energies of each system, although the additivity property is assumed to be satisfied by the individual micro-states. This non-additive property of the mean energy in the q-entropy formalism is manifestly incorrect, and not surprisingly leads to unphysical results, but the consequences have been ignored in various applications of q-entropy in the literature, as will be illustrated in some of the examples discussed here.

The definition of q-entropy for a thermodynamic system with micro-states labelled by an index $`i`$ is given by tsallis1 ; curado1 ; tsallis2 ; tsallis3

$`S_q=k{\displaystyle \frac{(1_ip_i^q)}{(q1)}},`$ (1)
where $`k`$ is a constant, $`q`$ is an undetermined parameter, and the quantities $`p_i`$ are positive numbers which satisfy the condition $`_ip_i=1`$. In the limit that $`q1`$ one recovers the Boltzmann-Gibbs form of the entropy

$`S_1=k_B{\displaystyle \underset{i}{}}p_iln(p_i),`$ (2)
where $`k=k_B`$ is the Boltzmann constant. In this special case, $`p_i`$ is the probability for the occurrence of the $`i`$-th micro-state, and this identification has been extended to the case $`q1`$. This extension, however, is not valid as can be seen from the definition of mean values for physical quantities associated with the q-entropy. For example, the internal energy $`U_q`$ is given by the form

$`U_q={\displaystyle P_iϵ_i},`$ (3)
where $`ϵ_i`$ corresponds to the $`i`$-th energy eigenvalue of the system, and $`P_i=p_i^q/_jp_j^q`$. The quantities $`P_i`$ are called “escort” probabilities in the q-entropy literature, but according to the conventional definition of mean value in statistics, these $`P_i^{}s`$ are the actual probabilities for the states of the system. Hence for $`q1`$ the $`p_i^{}s`$ introduced in the definition of q-entropy, Eq. 1, are devoid of any physical meaning, and are just functions of the probabilities $`P_j`$ according to the relation

$`p_i=P_i^{1/q}/({\displaystyle \underset{j}{}}P_j^{1/q}).`$ (4)
For values of $`q1`$, the q-entropy expression introduced in Eq. 1 is shown to be non-additive by the following arguments tsallis1 -tsallis3 . Suppose that two thermodynamic systems $`A`$ and $`B`$ are weakly coupled or are the subsystems of a larger system, and assume that the joint probabilities for the states of the combined system are the product of the probabilities for the states of the individual systems. As we shall see later on, this fundamental factorization property is not satisfied by the q-entropy formalism, but surprisingly this fact has been ignored in the literature. According to Eh. 4, we then have

$`p_{i,j}^{AB}=p_i^Ap_j^B.`$ (5)
Substituting this form into the expression for q-entropy, Eq. 1, one obtains the relation

$`S_q(AB)=S_q(A)+S_q(B)+(1q)S_q(A)S_q(B)/k.`$ (6)
However, this non-additive relation for the q-entropy leads immediately to a difficulty in the interpretation of $`S_q`$ as an expression for the thermodynamic entropy guerb1 . Since weak coupling means that the energy eigenvalues of the combined systems are essentially additive, we have

$`ϵ_{i,j}^{AB}=ϵ_i^A+ϵ_j^B,`$ (7)
and, according to Eqs.3 and 5, the total mean energy $`U_q(AB)`$ of the combined system is also additive:

$`U_q(AB)=U_q(A)+U_q(B).`$ (8)
Assume now that the combined system is isolated while there is an infinitesimal exchange of energy between systems $`A`$ and $`B`$. Then the variations $`\delta S_q(AB)=0`$ and $`\delta U_q(AB)=0`$, which implies that

$`{\displaystyle \frac{\delta S_q(A)}{1+(1q)S_q(A)/k}}={\displaystyle \frac{\delta S_q(B)}{1+(1q)S_q(B)/k}},`$ (9)
and

$`\delta U_q(A)=\delta U_q(B).`$ (10)
Combining these two equations, one finds that

$`[1+(1q)S_q(A)/k]T(A)=[1+(1q)S_q(B)/k]T(B),`$ (11)
where $`T(A)`$ and $`T(B)`$ are the absolute temperatures defined by the standard thermodynamic relation

$`{\displaystyle \frac{S}{U}}={\displaystyle \frac{1}{T}}.`$ (12)
For two systems in thermal contact these two temperatures should be equal, but according to Eq.11 this applies only if $`q=1`$. This problem is not unexpected, because temperature should be an an intensive quantity, but this is not possible in a formalism where the energy is additive, Eq.10, while the entropy does not satisfy this property. To avoid this problem, it has been proposed tsallis3 ; martinez2 ; casas1  to re-define absolute temperature as the quantity

$`T_q=(1+(1q)S_q/k)T`$ (13)
in which case the condition for thermal equilibrium, $`T_q(A)=T_q(B)`$, is satisfied by Eq. 11. But this definition of temperature, which must be universal, cannot be extended to systems $`A`$ and $`B`$ which are described by q-entropies with different values $`q_A`$ and $`q_B`$. In this case the q-entropy of the combined system, which is characterized by the quantities $`p_{i,j}^{AB}`$, Eq. 5, is not defined in term of these two $`q`$ parameters. Following the q-entropy formalism, one would have to introduce a new parameter $`q^{}`$ for the the q-entropy and energy of the combined system, but then these thermodynamic variables cannot be expanded in terms of the corresponding variables for the component systems $`A`$ and $`B`$ sasaki1 , 111 In reference sasaki1  the authors define the entropy of the joint system $`AB`$ in terms of an ad-hoc function of the q-entropies $`S_{q_A}(A)`$ and $`S_{q_B}(B)`$ of systems $`A`$ and $`B`$. This procedure, however, ignores the consistency requirement that this joint system must also be described by the same q-entropy relation, Eq.1. Instead, according to Eq.6, the q-entropy of the combined system would be given in terms of pseudo q-entropies for systems $`A`$ and $`B`$ with the same new parameter $`q^{}`$, and a similar problem would occur with the expansion of the total energy, Eq.8, i.e.

$`S_q^{}(AB)=S_q^{}(A)+S_q^{}(B)(1q^{})S_q^{}(A)S_q^{}(B)/k,`$ (14)
and

$`U_q^{}(AB)=U_q^{}(A)+U_q^{}(B).`$ (15)
Consequently, the concepts of thermal equilibrium, temperature, and heat exchange cannot be formulated for such systems. For example, no meaning can be attached to the statement that a system described by the Boltzmann-Gibbs entropy is in thermal equilibrium with a system described by q-entropy with $`q1`$. In other words, a Boltzmann-Gibbs thermometer would not be able to measure the temperature of a q-entropic system, and the laws of thermodynamics would therefore fail to have general validity 222In reference tsallis3 , Tsallis states that “ if we have in thermal contact systems with different entropic indices, say $`q_A`$ and $`q_B`$, it seems plausible \[my italics\] that at equilibrium $`T_{q_A}(A)=T_{q_B}(B)`$.” But as we have shown here, such a relation is not justified by the q-entropy formalism..

It follows that the parameter $`q`$ must be a universal constant, just like the Boltzmann constant $`k`$, which is applicable to all systems in thermodynamic equilibrium. If q is universal, a thermodynamic formulation for an infinitesimal reversible transfer of heat $`dQ`$ can be given between systems $`A`$ and $`B`$, with

$`dQ=T(A)dS_q(A)=T(B)dS_q(B),`$ (16)
corresponding to an exchange $`dU_q(A)=dU_q(B)`$ in the internal energy of these systems. But one is faced with the problem that in this case the temperature $`T(A)`$ is not equal to $`T(B)`$, which violates a fundamental principle of thermodynamics for systems in thermal equilibrium. Moreover, the corresponding differentials $`T_q(A)dS_q(A)`$ and $`T_q(B)dS_q(B)`$ associated with the proposed re-definition of absolute temperature, Eq.13, do not have any physical significance. In principle, this problem can be solved abe3  by introducing a different form $`S_q^R`$ for the q-entropy in the thermodynamic relation for the temperature $`T_q`$, so that

$`{\displaystyle \frac{S_q^R}{U_q}}={\displaystyle \frac{1}{T_q}}.`$ (17)
From the definition of $`T_q`$ given by Eq. 13 it follows that

$`S_q^R={\displaystyle \frac{k}{(1q)}}ln(1+(1q)S_q),`$ (18)
This alternative expression for the q-entropy was introduced by A. Réyni reyni1  in the form

$`S_q^R={\displaystyle \frac{k}{(1q)}}ln({\displaystyle \underset{i}{}}p_i^q),`$ (19)
which is additive 333 In reference abe3  the author concludes “that at the level of macroscopic thermodynamics, transmutation occurs from Tsallis theory to Réyni-entropy based theory”. This appears to be a convoluted admission that a non-additive form of q-entropy such as Eq. 1 is incompatible with the laws of thermodynamics, as can be verified by substituting for the $`p_i^{}s`$ the form in Eq. 5.

However, there are additional problems with either of these two definitions of q-entropy, even when $`q_A=q_B=q`$. Maximizing such entropy functions subject to the constraint of constant energy, Eq. 3 yields the q-probability distribution tsallis1 ; curado1 ; tsallis2 ; tsallis3

$`p_i[1(1q)\beta \delta ϵ_i]^{1/(1q)},`$ (20)
where $`\delta ϵ_i=ϵ_iU_q`$, $`U_q`$is the mean energy, Eq. 3, and $`\beta `$ is a parameter related to the inverse temperature. It is clear that the corresponding distribution $`p_{i,j}^{AB}`$ for the combined system $`AB`$ does not factor into the product $`p_i^Ap_j^B`$ even when the energies of the micro-states are additive, Eq. 7, unless $`q=1`$. For example, to first order in $`q1`$ 444A similar expression is obtained in casas1  which is, however, not quite correct, because in the second term a factor $`p(ϵ_i^A)p(ϵ_j^B)`$ is left out  one finds that

$`p(ϵ_i^A+ϵ_j^B)=p(ϵ_i^A)p(ϵ_j^B)[1+(q1)\beta ^2(\delta ϵ_i^A\delta ϵ_j^B)]`$ (21)
In the limit $`q=1`$ Eq. 20 reduces to the Boltzmann-Gibbs exponential form for the probability gibbs1 ,

$`p_iexp(ϵ_i/kT),`$ (22)
where $`\beta =1/kT`$. As is well known, this canonical distribution follows uniquely from the factorization requirement that

$`p(ϵ_i^A)p(ϵ_j^B)=p(ϵ_i^A+ϵ_j^B).`$ (23)
The Boltzmann-Gibbs form for the entropy $`S`$, Eq. 2, then follows from the assumption that

$`S={\displaystyle \underset{i}{}}f(p_i).`$ (24)
where the function $`f(p)`$ is determined uniquely. Applying the thermodynamic definition of temperature, $`S/U=1/T`$ and Boltzmann-Gibbs definition of the probability distribution, Eq. 22, one obtains the relation

$`{\displaystyle \underset{i}{}}p_i(ϵ_iU)[{\displaystyle \frac{df(p_i)}{dp_i}}+kln(p_i)]=0.`$ (25)
Hence

$`f(p_i)=kp_iln(p_i)`$ (26)
provided $`f(1)=f(0)=0`$, which corresponds to the requirement that the entropy vanishes at $`T=0`$.

## Unphysical properties resulting from applications of q-entropy to thermodynamic systems

To illustrate the consequences of disregarding such basic considerations, I would like to call attention to some unphysical results, left unmentioned in the literature, that follow from recent applications of q-entropy to some well-known thermodynamic systems. For example, many papers have been published on the application of q-entropy to black-body radiation black1 ; black2 ; black3 ; black4 ; black5 ; black6 ; black7 ; black8 . After laborious analysis and numerical computations, the authors in references black5 ; black6 ; black7  find that for $`q1`$ there are deviations from the well-known Stefan-Boltzmann law which states that the radiation energy depends on the fourth power of the temperature 555 The calculations in references black1 ; black2 ; black3 ; black4 ; black5  do not exhibit deviations from the Stefan-Boltzmann law because these calculations have left out an anomalous dependence on $`V/V(T)`$, where $`V`$ is the volume of the cavity and $`V(T)=(\mathrm{}c/kT)^3`$. It is clear, that such a dependence must appear, because both the energy and the entropy density are non-extensive for $`q1`$ (except to lowest order in an expansion in $`q1`$, in which case the energy is extensive).. But since Boltzmann derived this result from purely thermodynamical reasoning, without any statistical assumptions about the form of the entropy, it seems at first sight strange that such a deviation can occur in a formalism which is supposed to satisfy the laws of thermodynamics. The explanation is that as a consequence of the non-extensivity property of the q-entropy formalism, the black-body energy as well as the entropy do not dependent linearly on the volume of the cavity, as was originally assumed by Boltzmann. Hence, the q-energy density and the q-entropy density depend on the volume $`V`$ of the cavity, although there is no comment concerning this unphysical property in any of the q-entropy calculations in references black1 ; black2 ; black3 ; black4 ; black5 ; black6 ; black7 ; black8 . As follows from simple dimensional arguments, the volume dependence must be given by a dimensionless parameter $`V/a(T)^3`$, where $`a(T)`$ is a characteristic length which can depend only on the temperature $`T`$. From statistical mechanics we learn that $`a(T)=(\mathrm{}c/kT)`$ in fact corresponds to the mean thermal wavelength of the black-body photons. Moreover, expansions in a power series of this parameter which have been applied to fit the cosmic background radiation black1 ; black2 ; black3 ; black4 ; black5 ; black6 ; black7 ; black8  are nonsensical,because in this case $`a(T)`$ is of order $`1/10`$ $`cm`$ while the cavity volume $`V`$ has cosmological dimensions! On purely thermodynamic grounds it can also be shown that if the temperature dependence of the black-body energy density were to have the form $`uT^{4+\delta }`$, then the Maxwell-Boltzmann relation $`p=(1/3)u`$ for the thermal radiation pressure $`p`$ would lead to a power law volume dependence $`uV^{\delta /3}`$, and correspondingly $`UT(T^3V)^{1+\delta /3}`$, in accordance with our previous dimensional argument. Similarly, one finds that $`sV^{\delta /3}`$ for the entropy density. But unless $`\delta =0`$, such a volume dependence is incompatible with Kirchhoff’s law which states that the ratio of emissivity to absorption of radiation in the walls of the cavity must be a universal function of the temperature, and the frequency of the radiation. These properties are required in order that such a cavity reach thermal equilibrium. Historically, this law was the original basis for the universality properties of black-body radiation that culminated in Planck’s famous derivation.

Actually, the application of the q-entropy formalism to any system which consists of non-interacting components leads to unphysical properties. For example, for an ideal gas the energy is found to depend nonlinearly on the number of particles, the gas pressure $`p`$ is not equal to $`2/3`$ of the energy density, and there are correlations between the energies of any pair of particlesideal1 ; ideal2 ; ideal3 , contrary to very well-known result of kinetic gas theory and statistical mechanics. Similarly, for a system of non-interacting magnetic moments anomalous dependences are found for the magnetization and susceptibility on the number of spins spin1 , although claims have been also made that the magnetization reproduces experimental results in certain manganites spin2 ; spin3  The reason for the failure of this formalism to give physically sensible results is not hard to see - it is due to the fact that the probability distribution, Eq. 20, for any two components $`A_1`$ and $`A_2`$ is not the product of the probabilities for the separate components as would be expected if these components are non-interacting. Hence, contrary to basic physical principles, this formalism gives rise to spurious correlations among these components of such systems.

Some proponents of q-entropy have argued that this formalism should be considered only for systems for which the Boltzmann-Gibbs thermodynamic formalism supposedly “fails”. Frequently mentioned as candidates are systems with components that interact primarily through long range forces such as gravitational forces for which the total energy and entropy are non-extensive. In Nature, such system correspond to astrophysical objects such as stars and galaxies, but these objects are generally not found at maximum entropy and (correspondingly) uniform temperature. For example, stars are either evolving slowly in time, like main sequence stars, which have a large temperature gradient from the interior to the surface, and emit thermal radiation (black-body energy and entropy), or else have reached a degenerate state such as dwarf stars or neutron stars provided certain mass limits are satisfied chapline . Otherwise, stars eventually explode into supernovas, sometimes leaving remnants which collapse into such degenerate states or into a black hole (which are states having maximum entropy). The equilibrium properties of stars are obtained by hydrostatic equations supplemented by local thermodynamic equations for matter and radiation based on the Boltzmann-Gibbs entropy weisskopf  666 According to Boltzmann-Gibbs statistics, a gas of particles interacting via gravitational forces in thermal equilibrium (maximum entropy) would have to have infinite size and mass For the solution of the isothermal gas polytrope see chandra1 . Rather than being a ”failure”, the Boltzmann-Gibbs statistics has been applied to stellar structure with enormous success. A q-entropy formalism, however, predicts the existence of finite isothermal polytropes as the end products of stellar evolution ( states of maximum entropy) which fails completely to account for the observed property of stars in our universe. Such q-polytrope solutions have also been discussed as models for galaxies plastinograv ; taruya ; chavanis .

Another example which has been cited as a so-called failure of Boltzmann-Gibbs entropy because it involves long-range electromagnetic forces, is the divergence of the partition function $`Z`$ calculated for the bound states of the hydrogen atom hydrogen1 . In this case

$`Z={\displaystyle \underset{n}{}}exp(ϵ_n/T),`$ (27)
where $`ϵ_n=R/n^2`$ are the bound-state energy levels, and $`R`$ is the Rydberg constant. But the divergence of $`Z`$, which occurs in this case because the terms of the series approach unity as $`n`$ becomes large, is related to the growth of the the mean radius of the hydrogen atom which increases as $`n^2`$. Obviously, in a gas of hydrogen atoms in thermal equilibrium, this radius can not become larger than the mean distance between atoms. Therefore this distance provides an effective cutoff for the applicability of the hydrogen bound-state energy eigenvalues in the partition sum, because for larger values of $`n`$ these atoms can no longer be treated even approximately as a gas of non-interacting particles. Instead, for these states the gas must be viewed as a neutral plasma of electrons and protons interacting via long range electromagnetic forces. Hence the q-entropy formalism, which supposedly gives a finite partition function hydrogen1 , actually fails to account for the correct physics of this problem. These and other failures neutrinos1 , 777 It is frequently mentioned tsallis2  that the solar neutrino problem discussed in reference neutrinos1  is an example of “plain failures” of the Boltzmann-Gibbs entropy, because the calculated rates of nuclear reactions in the Sun supposedly do not account for the neutrino observations. Recently, however, the SNO data (Q.R. Ahmad et. al., Physical Review Letters 89 (2002)) show that there is now excellent agreement between theory and neutrino observations. Among other things, this result is a spectacular new confirmation of the applicability of Boltzmann-Gibbs entropy in the interior of a star. in the application of q-entropy to well-known physical systems mirror the inconsistencies which are inherent in a formulation of thermodynamics based on q-entropy, Eq. 1.

## .1 Conclusion

We have shown that a prerequisite to have a q-entropy formalism which is consistent with the laws of thermodynamics is that the parameter $`q`$ must be a universal constant, as is the case also with the constant $`k`$ that corresponds to Boltzmann’s constant for $`q=1`$. Moreover, for weakly coupled systems the entropy as well as the energy must be additive, which is a condition explicitly violated by the Tsallis entropy, Eq. 1. Although this condition appears to be satisfied by the Réyni q-entropy, Eq. 19, this is actually not the case, because the probability distribution obtained by maximizing this entropy for $`q1`$, Eq. 20, does not satisfy the required factorization condition, Eq. 5, which is required for weakly coupled systems. Thus, we have shown that the only value of $`q`$ consistent with the laws of thermodynamics is $`q=1`$, which corresponds to the familiar Boltzmann-Gibbs form for the entropy, Eq. 2. Indeed, it has been explicitly demonstrated here that the application of a q-entropy formalism to black-body radiation, and to other systems with weakly interacting components, leads to unphysical results when $`q1`$.

It has been suggested in the literature that q-entropy calculations are useful because they provide an additional parameter $`q`$ for comparing theory with observations, but this rationale fails to take into account the fact that such calculations would be inconsistent with the fundamental principles of thermodynamics and statistical mechanics. Hence, if a small departure of $`q`$ from unity is found in a fit to data, as is claimed for example in various analyses of the cosmic black-body radiation black1 ; black2 ; black3 ; black4 ; black5 ; black6 ; black7 ; black8 , such a fit can not provide any physical insight whatsoever into the source or meaning of the deviations.

Finally, we remark that the q-probability distribution, Eq. 20, which is obtained by maximizing either of the two q-entropy relations, Eqs. 1 and 19, constrained to a fixed pseudo-energy function, has also been applied to non-equilibrium problems. For example, data on turbulence in a pure-electron plasma column bogo1  and the velocity distribution in turbulent flow bodenschatz  swinney  have been fitted by such a probability function beck1 ; beck2 . But absolutely no physical justification has been given for applying to non-equilibrium systems a fundamental condition - maximum entropy - which is associated in statistical mechanics with systems in thermal equilibrium. While such a probability distribution can also be obtained from other ad-hoc assumptions beck1 ; beck2 , its connection to q-entropy and non-extensivity is completely unfounded. Although the q-probability distribution appears to be a good phenomenological parameterization for some turbulence data, its deduction from physical principles has not been established 888 In references beck1 ; beck2  the author expressed the q-probability distribution, Eq. 20, as a Laplace integral, which is the basis for his suggestion that this distribution is due to temperature fluctuations in ordinary Boltzmann-Gibbs statistics. But no connection based on physics has been made between these assumed temperature fluctuations and q-entropy. .

## Acknowledgments

I would like to thank S. Abe, C.Beck, E.G. Cohen, A. Plastino, J. Riera and C. Tsallis for comments, and for providing references to the extensive literature on non-extensive q-entropy.


format: text

text:
# On the singularities of gravity in the presence of non-minimally coupled scalar fields

## I Introduction

Recently, we have investigated the dynamics of homogeneous and isotropic solutions of the cosmological model described by the action:

$$S=d^4x\sqrt{g}\left\{F(\varphi )R_a\varphi ^a\varphi 2V(\varphi )\right\},$$
(1)
with $`F(\varphi )=1\frac{1}{6}\varphi ^2`$, corresponding to the so-called conformal coupling, and $`V(\varphi )=\frac{m}{2}\varphi ^2\frac{\mathrm{\Omega }}{4}\varphi ^4`$. Some novel dynamical behaviors were identified: superinflation regimes, a possible avoidance of big-bang singularities through classical birth of the universe from empty Minkowski space, spontaneous entry into and exit from inflation, and a cosmological history suitable for describing quintessence. Since one of the proposals of inflationary models is to describe our universe without finely-tuned parameters, these results would be viable only if they are robust against small perturbations in initial conditions and in the model itself. This is the question to be addressed here.

We study the robustness of our previous results by taking two generalizations of the model considered previously: we relax the isotropy requirements (perturbations in the initial conditions) and we consider a general coupling $`F(\varphi )`$ (perturbations in the model parameters). Models with more general $`F(\varphi )`$ have been considered recently. Our results show that the model is not robust. Its main properties are radically changed, even for small disturbances in initial conditions and in the model itself, due to the appearance of real, gravitational singularities that are dynamically unavoidable in general. The singularities are, essentially, of two types. The first one corresponds to the hypersurfaces $`F(\varphi )=0`$. It is not present in the isotropic case, and it implies that all previous homogeneous and isotropic solutions passing from the $`F(\varphi )>0`$ to the $`F(\varphi )<0`$ region are extremely unstable against anisotropic perturbations. The second type of singularity corresponds to $`F_1(\varphi )=0`$, with

$$F_1(\varphi )=F(\varphi )+\frac{3}{2}\left(F^{}(\varphi )\right)^2,$$
(2)
and it is present even for the homogeneous and isotropic cases. Although for small deviations of the conformal coupling the latter singularities are typically very far from the region of interest, in the general case they can alter qualitatively the global dynamics of the model due to restrictions that it imposes on the phase space. Again, the persistence of some of our previously described results, in particular the ones concerning heteroclinic and homoclinic solutions, are challenged.

Both kinds of singularities have already been described before. To the best of our knowledge, Starobinski was the first to identify the singularity corresponding to the hypersurfaces $`F(\varphi )=0`$, for the case of conformally coupled anisotropic solutions. Futamase and co-workers identified both singularities in the context of chaotic inflation in $`F(\varphi )=1\xi \varphi ^2`$ theories (See also ). The first singularity is always present for $`\xi >0`$ and the second one for $`0<\xi <1/6`$. Our conclusions are, however, more general since we treat the case of general $`F(\varphi )`$ and our results are based on the analysis of true geometrical invariants. Our main result is that the system governed by (1) is generically singular on both hypersurfaces $`F(\varphi )=0`$ and $`F_1(\varphi )=0`$. Here, generically means that it is possible to construct non-singular models if one fine-tunes $`F(\varphi )`$ and $`V(\varphi )`$, as we will show below.

One can advance that there are some geometrically special regions on the phase space of the model in question by a very simple analysis of the equations derived from the action (1). They are the Klein-Gordon equation

$$\mathrm{}\varphi V^{}(\varphi )+\frac{1}{2}F^{}(\varphi )R=0,$$
(3)
and the Einstein equations

$`F(\varphi )G_{ab}`$ $`=`$ $`(1+F^{\prime \prime }(\varphi ))_a\varphi _b\varphi `$ (4)
$``$ $`{\displaystyle \frac{1}{2}}g_{ab}\left[(1+2F^{\prime \prime }(\varphi ))_c\varphi ^c\varphi +2V(\varphi )\right]F^{}(\varphi )\left(g_{ab}\mathrm{}\varphi _a\varphi _b\varphi \right).`$ (5)
We will consider here the simplest anisotropic homogeneous cosmological model, the Bianchi type I, whose spatially flat metric is given by

$$ds^2=dt^2+a^2(t)dx^2+b^2(t)dy^2+c^2(t)dz^2.$$
(6)
The dynamically relevant quantities here are

$$H_1=\frac{\dot{a}}{a},H_2=\frac{\dot{b}}{b},\mathrm{and}H_3=\frac{\dot{c}}{c}.$$
(7)
For such a metric and a homogeneous scalar field $`\varphi =\varphi (t)`$, after using the Klein-Gordon Eq. (3), Eq. (4) can be written as

$`F(\varphi )G_{00}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{\varphi }^2+V(\varphi )F^{}(\varphi )\left(H_1+H_2+H_3\right)\dot{\varphi },`$ (8)
$`{\displaystyle \frac{1}{a^2}}F(\varphi )G_{11}`$ $`=`$ $`{\displaystyle \frac{1+2F^{\prime \prime }(\varphi )}{2}}\dot{\varphi }^2V(\varphi )F^{}(\varphi )\left(H_1\dot{\varphi }+V^{}(\varphi ){\displaystyle \frac{F^{}(\varphi )}{2}}R\right),`$ (9)
$`{\displaystyle \frac{1}{b^2}}F(\varphi )G_{22}`$ $`=`$ $`{\displaystyle \frac{1+2F^{\prime \prime }(\varphi )}{2}}\dot{\varphi }^2V(\varphi )F^{}(\varphi )\left(H_2\dot{\varphi }+V^{}(\varphi ){\displaystyle \frac{F^{}(\varphi )}{2}}R\right),`$ (10)
$`{\displaystyle \frac{1}{c^2}}F(\varphi )G_{33}`$ $`=`$ $`{\displaystyle \frac{1+2F^{\prime \prime }(\varphi )}{2}}\dot{\varphi }^2V(\varphi )F^{}(\varphi )\left(H_3\dot{\varphi }+V^{}(\varphi ){\displaystyle \frac{F^{}(\varphi )}{2}}R\right).`$ (11)
It is quite simple to show that Eqs. (9)-(11) are not compatible, in general, on the hypersurface $`F(\varphi )=0`$. Subtracting (10) and (11) from (9) we have, on such hypersurface, respectively,

$$F^{}(\varphi )(H_1H_2)\dot{\varphi }=0,\mathrm{and}F^{}(\varphi )(H_1H_3)\dot{\varphi }=0.$$
(12)
Hence, they cannot be fulfilled in general for anisotropic metrics. As it will be shown, it indeed corresponds to an unmovable (in the Painlevé sense) geometrical singularity which cannot be prevented in general by requiring that $`F^{}(\varphi )=0`$ or $`\dot{\varphi }=0`$ on the hypersurface.

As to the second singularity we have, after taking the trace of the Einstein equations, that:

$$R=R(\varphi ,\dot{\varphi })=\frac{1}{F_1(\varphi )}\left(4V(\varphi )+3V^{}(\varphi )F^{}(\varphi )(1+F^{\prime \prime }(\varphi ))\dot{\varphi }^2\right).$$
(13)
Inserting Eq. (13) in the Klein-Gordon Eq. (3), one can see that it contains terms which are singular for $`F_1(\varphi )=0`$. Again, as we will see, this corresponds to an unmovable geometrical singularity, and it cannot be eliminated, in general, by demanding that $`F^{}(\varphi )=0`$ on the hypersurface $`F_1(\varphi )=0`$. In both the hypersurfaces $`F(\varphi )=0`$ and $`F_1(\varphi )=0`$ the Cauchy problem is ill-posed, since one cannot choose general initial conditions.

The hypersurfaces $`F(\varphi )=0`$ and $`F_1(\varphi )=0`$ also prevent the global definition of an Einstein frame for the action (1), defined by the transformations

$`\stackrel{~}{g}_{ab}`$ $`=`$ $`F(\varphi )g_{ab},`$ (14)
$`\left({\displaystyle \frac{d\stackrel{~}{\varphi }}{d\varphi }}\right)^2`$ $`=`$ $`{\displaystyle \frac{F_1(\varphi )}{2F(\varphi )^2}}.`$ (15)
It is well known that in the Einstein frame the Cauchy problem is well posed. Again, the impossibility of defining a global Einstein frame shed some doubts about the general Cauchy problem. Moreover, the standard perturbation theory for helicity-2 and helicity-0 excitations, derived directly from Eqs. (14)-(15), fails in both hypersurfaces.

## II The singularities

In order to check the geometrical nature of these singular hypersurfaces, let us consider the Einstein Eqs. (8)-(11) in detail. For the metric (6), we have the following identities

$`G_{00}`$ $`=`$ $`H_1H_2+H_2H_3+H_1H_3,`$ (16)
$`G_{11}`$ $`=`$ $`a^2\left(\dot{H}_1+H_1(H_1+H_2+H_3){\displaystyle \frac{1}{2}}R\right),`$ (17)
$`G_{22}`$ $`=`$ $`b^2\left(\dot{H}_2+H_2(H_1+H_2+H_3){\displaystyle \frac{1}{2}}R\right),`$ (18)
$`G_{33}`$ $`=`$ $`c^2\left(\dot{H}_3+H_3(H_1+H_2+H_3){\displaystyle \frac{1}{2}}R\right),`$ (19)
$`R`$ $`=`$ $`2\left(\dot{H}_1+\dot{H}_2+\dot{H}_3+H_1^2+H_2^2+H_3^2+H_1H_2+H_2H_3+H_1H_3\right).`$ (20)
After using expressions (16) and introducing the new dynamical variables $`p=H_1+H_2+H_3`$, $`q=H_1H_2`$, and $`r=H_1H_3`$, Einstein Eqs. can be cast in the form

$`E(\varphi ,\dot{\varphi },p,q,r)`$ $`=`$ $`{\displaystyle \frac{1}{3}}F(\varphi )\left(p^2+qrq^2r^2\right)+{\displaystyle \frac{\dot{\varphi }^2}{2}}+V(\varphi )pF^{}(\varphi )\dot{\varphi }=0,`$ (21)
$`\dot{q}`$ $`=`$ $`\left(p+{\displaystyle \frac{F^{}(\varphi )}{F(\varphi )}}\dot{\varphi }\right)q,`$ (22)
$`\dot{r}`$ $`=`$ $`\left(p+{\displaystyle \frac{F^{}(\varphi )}{F(\varphi )}}\dot{\varphi }\right)r,`$ (23)
$`2F_1(\varphi )\dot{p}`$ $`=`$ $`(F(\varphi )+2F^{}(\varphi )^2)p^2+{\displaystyle \frac{3}{2}}(1+2F^{\prime \prime }(\varphi ))\dot{\varphi }^23V(\varphi )3F^{}(\varphi )V^{}(\varphi )`$ (25)
$`p\dot{\varphi }F^{}(\varphi )+(F(\varphi )+F^{}(\varphi )^2)(q^2+p^2qr)`$
Using the energy constraint (21), this last equation can be put in the form

$`\dot{p}`$ $`=`$ $`\left(p+{\displaystyle \frac{F^{}(\varphi )}{F(\varphi )}}\dot{\varphi }\right)p+3{\displaystyle \frac{F(\varphi )+\frac{1}{2}(F^{}(\varphi ))^2}{F(\varphi )F_1(\varphi )}}V(\varphi )`$ (27)
$`{\displaystyle \frac{3}{4}}\left({\displaystyle \frac{2F(\varphi )F^{\prime \prime }(\varphi )(F^{}(\varphi ))^2}{F(\varphi )F_1(\varphi )}}\right)\dot{\varphi }^2+{\displaystyle \frac{3}{2}}{\displaystyle \frac{F^{}(\varphi )}{F_1(\varphi )}}V^{}(\varphi ).`$
The Klein-Gordon equation (3) reads simply

$$\ddot{\varphi }+p\dot{\varphi }+V^{}(\varphi )\frac{F^{}(\varphi )}{2}R(\varphi ,\dot{\varphi })=0,$$
(28)
with $`R(\varphi ,\dot{\varphi })`$ given by Eq. (13). The energy constraint, equation (21), is evidently compatible with the other ones. Indeed $`E(\varphi ,\dot{\varphi },p,q,r)=0`$ is an invariant surface since one has that

$$\frac{d}{dt}E(\varphi ,\dot{\varphi },p,q,r)=\left(2p+\frac{F^{}(\varphi )}{F(\varphi )}\dot{\varphi }\right)E(\varphi ,\dot{\varphi },p,q,r)$$
(29)
along solutions of Eqs. (22), (23), (27), and (28). Note that Eqs. (27) and (28) are decoupled from the equations for $`\dot{q}`$ and $`\dot{r}`$. Equations (22) and (23) are, hence, linear first order equations, and they could be easily integrated after the solutions of (27) and (28) have been found. Moreover, since one has $`r\dot{q}q\dot{r}=0`$, $`q(t)/r(t)`$ is a constant of motion fixed only by the initial conditions. Suppose the initial ratio is $`q(0)/r(0)=\gamma `$: this would imply that $`(H_1H_2)=\gamma (H_1H_3)`$ for all $`t`$, leading to, for instance, $`c^\gamma (t)a^{\gamma 1}(t)b(t)`$ in the metric (6). This simplification is a consequence of the scalar character of our source field, and it does not suppose any loss of generality in our arguments.

A closer analysis of Eqs. (22)-(27) reveals the presence of the singularities. In general, the right-hand side of these equations diverge for $`F(\varphi )=0`$ and for $`F_1(\varphi )=0`$. One can check that these divergences are indeed related to real geometrical singularities by considering the Kretschman invariant $`I=R_{abcd}R^{abcd}`$, which for the metric (6) is given by

$$I=4\left(\left(\dot{H}_1+H_1^2\right)^2+\left(\dot{H}_2+H_2^2\right)^2+\left(\dot{H}_3+H_3^2\right)^2+H_1^2H_2^2+H_1^2H_3^2+H_2^2H_3^2\right).$$
(30)
As one can see, $`I`$ is the sum of non negative terms. Moreover, any divergence of the variables $`H_1`$, $`H_2`$, $`H_3`$, or of their time derivatives, would suppose a divergence in $`I`$, characterizing a real geometrical singularity. Since the relation between the variables $`p`$, $`q`$, $`r`$, and $`H_1`$, $`H_2`$, $`H_3`$ is linear, any divergence of the first, or of their time derivative, will suppose a divergence in $`I`$.

Suppose, now, that $`F(\varphi _0)=0`$, and that $`F(\varphi )`$ is (real) analytical for $`\varphi =\varphi _0`$. In this case, $`F^{}(\varphi )/F(\varphi )`$ diverges as $`(\varphi \varphi _0)^1`$ near $`\varphi _0`$ for nonvanishing $`q`$ and $`r`$, rendering $`\dot{q}`$ and $`\dot{r}`$ divergent by (22) and (23), unless $`p`$ also diverges in order to keep $`\left(p+\frac{F^{}(\varphi )}{F(\varphi )}\dot{\varphi }\right)`$ finite on $`\varphi _0`$. In both cases, $`I`$ diverges. There is no dynamical restriction to ensure that $`\dot{\varphi }`$ vanishes on $`\varphi _0`$ \- it can take any value compatible with the energy constraint (21). Indeed, the latter implies that on $`\varphi _0`$

$$\frac{\dot{\varphi }^2}{2}pF^{}(\varphi _0)\dot{\varphi }+V(\varphi _0)=0.$$
(31)
There is no way of having $`\dot{\varphi }=0`$ on $`\varphi _0`$, unless $`V(\varphi _0)=0`$, and even in this case, $`\dot{\varphi }=pF^{}(\varphi _0)`$ is also possible. Note that the hypothesis of $`F(\varphi )`$ analytical at $`\varphi _0`$ is not a necessary one. For any differentiable function $`F(\varphi )`$ with a zero in $`\varphi _0`$ one has $`|F(\varphi )|=|_{\varphi _0}^\varphi F^{}(s)𝑑s|k|\varphi \varphi _0|`$, with $`k=\mathrm{max}_{s[\varphi _0,\varphi ]}|F^{}(s)|`$, implying that $`|F^{}(\varphi )/F(\varphi )||\varphi \varphi _0|^1|F^{}(\varphi )|/k`$. Since $`F^{}(\varphi )`$ is assumed to be continuous, the last ratio tends to 1 when $`\varphi \varphi _0`$, implying the divergence of $`F^{}(\varphi )/F(\varphi )`$ in that limit.

Now, let us suppose $`F_1(\varphi _1)=0`$. If $`F^{}(\varphi _1)0`$, the right-hand side of Eq. (27) diverges. The vanishing of $`F^{}(\varphi _1)`$ implies, by (2), that $`F(\varphi _1)=0`$, and the arguments of the last paragraph can be repeated.

A singularity-free model can be constructed by demanding that $`F(\varphi _0)=F^{}(\varphi _0)=0`$, by choosing a $`V(\varphi )`$ that goes to 0 at a proper rate when $`\varphi \varphi _0`$, and by demanding that $`F_1(\varphi )`$ have no other zeros than the ones of $`F(\varphi )`$. Models for which $`F(\varphi )=\zeta \varphi ^{2n}`$ and $`V(\varphi )=\alpha \varphi ^{2(2n1)}+\mathrm{high}\mathrm{order}\mathrm{terms}`$, for instance, fulfill these requirements. However, such a highly fine-tuned class of model is of no physical interest here, since it does not contain $`F(\varphi )>0`$ and $`F(\varphi )<0`$ regions and consequently has no solution for which the effective gravitational constant $`G_{\mathrm{eff}}`$ changes it sign along the cosmological history. The stability of such solutions were the starting point of the analyses of the pioneering work  and of the present one as well. Note that by Eq. (2), models with an $`F(\varphi )<0`$ region will allways have singularities of the type $`F_1(\varphi )=0`$. This fact shall be taken into account to better understand the recently proposed dynamical stability of the $`F(\varphi )<0`$ region.

## III Conclusion

The singularities described in the precedent section imply that the model presented in  is not robust, since our main conclusions were a consequence of very especial initial conditions. For instance, all homogeneous and isotropic solutions crossing the $`F(\varphi )=0`$ hypersurface are extremely unstable against anisotropic perturbations. By Eqs. (22) and (23), any deviation from perfect isotropy (expressed by nonvanishing $`q`$ and $`r`$ variables) for these solutions, however small, will lead catastrophically to a geometrical singularity. Many of the novel dynamical behaviors presented in  depend on these solutions. This is the case, for instance, of some solutions exhibiting superinflation regimes. The heteroclinic and homoclinic solutions identified in  can cross the $`F(\varphi )=0`$ hypersurface and, hence, they also suffer the same instability against anisotropic perturbations. The homoclinic solutions were considered as candidates to describe a non-singular cosmological history, with the big-bang singularity being avoided through a classical birth of the universe from empty Minkowski space. Apart from $`F(\varphi )=0`$ singularities, these solutions are also affected by the singularities of type $`F_1(\varphi )=0`$. Suppose that the conformal coupling is disturbed by a very small negative term: $`F(\varphi )=1(\frac{1}{6}ϵ)\varphi ^2`$. The $`F_1(\varphi )=0`$ singularities will be near the $`\varphi =\pm 1/\sqrt{ϵ}`$ hypersurfaces. Although they are located far from the $`F(\varphi )=0`$ regions, they alter the global structure of the phase-space. In this case, they restrict the existence of homoclinics, rendering a non-singular cosmological history more improbable.

The singularities do not affect the conclusions obtained by considering solutions inside the $`F(\varphi )>0`$ region. The asymptotic solutions presented in , for instance, are still valid. The conclusion that for large $`t`$ the dynamics of any solution (inside $`F(\varphi )>0`$) tends to an infinite diluted matter dominated universe remains valid. Moreover, for small anisotropic deviations ($`q`$ and $`r`$ small in comparision with $`p`$), Eqs. (22) and (23) allow us to conclude that solutions inside $`F(\varphi )>0`$, for large $`t`$, approach exponentially isotropic matter-dominated universe.

###### Acknowledgements.

The authors wish to thank G. Esposito-Farese for previous discussions concerning the Cauchy problem for models with $`F(\varphi )=0`$, and Profs. Albert and Anny Sanfeld for the warm hospitality in Mallemort, France, where this work was initiated. They also acknowledge the financial support from the EEC (project HPHA-CT-2000-00015), from OLAM - Fondation pour la Recherche Fondamentale (Belgium), from Fondation Science et Environnement (France), and from CNPq (Brazil).


text:
# A simple model for drag reduction

## Abstract

Direct Numerical Simulations established that the FENE-P model of viscoelastic flows exhibits the phenomenon of turbulent drag reduction which is caused in experiments by dilute polymeric additives. To gain analytic understanding of the phenomenon we introduce in this Letter a simple 1-dimensional model of the FENE-P equations. We demonstrate drag reduction in the simple model, and explain analytically the main observations which include (i) reduction of velocity gradients for fixed throughput and (ii) increase of throughput for fixed dissipation.

The addition of few tens of parts per million (by weight) of long-chain polymers to turbulent fluids can bring about a reduction of the friction drag by up to 80% 00SW . This “drag reduction” phenomenon has important practical implications besides being interesting from the fundamental point of view, integrating turbulence research with polymer physics. In spite of intense interest for an extended period of time 69Lu ; 75Virk ; 90Ge , Sreenivasan and White 00SW  recently concluded that “it is fair to say that the extensive - and continuing - activity has not produced a firm grasp of the mechanisms of drag reduction”. Recently however it was shown that drag reduction is observed in Direct Numerical Simulation of model viscoelastic hydrodynamic equations 97THKN ; 98DSB ; 00ACP . From the theoretical viewpoint these observations are crucial, indicating that the phenomenon is included in the solutions of the model equations. Understanding drag reduction then becomes a usual challenge of theoretical physics. In this Letter we present a further simplification of the model equations and gain analytic insights into the phenomenon. The FENE-P equation for the fluid velocity $`𝒖(𝒓,t)`$ contains an additional stress tensor related to the polymer:

$$\frac{𝒖}{t}+(𝒖\mathbf{})𝒖=\mathbf{}p+\nu _s^2𝒖+\mathbf{}𝓣+𝑭,$$
(1)
where $`\nu _s`$ is the viscosity of the neat fluid, $`𝑭`$ is the forcing and the stress tensor $`𝓣`$ is determined by the polymer conformation tensor $`𝑹`$ according to

$$𝓣(𝒓,t)=\frac{\nu _p}{\tau _p}\left[\frac{f(𝒓,t)}{\rho _0^2}𝑹(𝒓,t)\mathrm{𝟏}\right].$$
(2)
here $`\nu _p`$ is a viscosity parameter, $`\tau _p`$ is a relaxation time for the polymer conformation tensor, $`\rho _0`$ is the rms extension of the polymers in equilibrium, and $`f(𝒓,t)`$ is a function that limits the growth of the trace of $`𝑹`$. The model is closed by the equation of motion for the conformation tensor which reads

$`{\displaystyle \frac{R_{\alpha \beta }}{t}}+(𝒖\mathbf{})R_{\alpha \beta }`$ $`={\displaystyle \frac{u_\alpha }{r_\gamma }}R_{\gamma \beta }+R_{\alpha \gamma }{\displaystyle \frac{u_\gamma }{r_\beta }}`$
$`{\displaystyle \frac{1}{\tau _p}}\left[f(𝒓,t)R_{\alpha \beta }\rho _0^2\delta _{\alpha \beta }\right]`$
These equations were simulated on the computer in a channel or pipe geometry. The main observations on the effect of the polymer on the turbulent flow that we need to focus on are the following:

(i) For a fixed pressure gradient at the wall the fluid throughput is increased (see Fig. 1). (ii) For a fixed throughput the gradient at the wall decreases (i.e. the dissipation decreases). (iii) The trace of the conformation tensor $`𝑹`$ follows qualitatively the rms streamwise velocity (see Fig. 2). We are particularly interested in point (iii) since in our opinion the space dependence of the amount of stretching (and with it of the effective viscosity) is crucial, and compare 01GLP ; 02GLP  for a discussion of this point in the context of the instability of laminar flows.

Obviously, Eqs. (1)-(A simple model for drag reduction) as they stand are not amenable to analytic investigation in the turbulent regime. To gain insight we therefore attempt to simplify them as much as possible without losing the main phenomena (i)-(iii). Consider therefore a model for the streamwise velocity which is the Burger’s equation ($`u`$ in the streamwise directions with gradients in the $`y`$ (wall-normal) direction), to which the effect of a scalar $`R`$ is added:

$`u_t+uu_y`$ $`=`$ $`\nu u_{yy}+sR_y+F,`$ (4)
$`R_t+uR_y`$ $`=`$ $`{\displaystyle \frac{1}{\tau }}R+Ru_y.`$ (5)
where a subscript $`y`$ stands for a partial derivative with respect to $`y`$. In the following we shall denote Eqs. (4)-(5) with the acronimous uR model. The parameter $`s`$ is related to the polymer concentration, and $`\tau `$ is the relaxation time of $`R`$. We will consider the model in the domain $`LyL`$, with boundary conditions chosen later. We will denote spatial averages by pointed brackets, $`A_L^LA(y)𝑑y`$. The simplicity of the uR model allows us to state the energy budget in simple terms. Multiplying (4) by $`u`$ and taking the spatial average of (4) and (5) we obtain:

$$\frac{1}{2}\frac{d}{dt}u^2=\nu u_y^2+suR_y+Fu$$
(6)
$$\frac{d}{dt}R=\frac{1}{\tau }R+2Ru_y$$
(7)
The term $`uR_y`$ measures the “energy” given by the velocity field $`u`$ to the polymer field $`R`$. Multiplying (7) by $`s/2`$ and summing (7) with (6) we obtain:

$$\frac{d}{dt}(\frac{1}{2}u^2+\frac{s}{2}R)=\nu u_y^2\frac{s}{2\tau }R+Fu$$
(8)
In the steady state the overall power $`Fu`$ is balanced by the overall energy dissipation per unit time $`D`$, $`D=\nu u_y^2+\frac{s}{2\tau }R`$. The term $`\frac{1}{2}u^2+\frac{s}{2}R`$ represents the sum of the kinetic energy of the flow plus the potential energy of the stretched polymers. We remark that already from these elementary consideration it becomes clear that for a fixed power input the existence of the term $`\frac{s}{2\tau }R`$ necessarily reduces the gradients of $`u`$ in agreement with point (ii) above. To address points (ii) and (iii) further we consider next the solution of the model with $`F=0`$ and with a fixed velocity $`u`$ and stretching $`R`$ at $`L`$ and $`L`$. In other words, we take as boundary conditions

$$u(L)=u_0,u(L)=u_0,R(L)=R(L)=R_0.$$
(9)
In Fig. 3 we compare the solution of the uR model to that of the pure Burger’s equation (i.e. Eq. (4) with $`R=F=0`$.) To focus our thinking we would like the reader to consider the solution in the left half space as a model of the streamwise velocity component in the lower half channel, with the solution in the right half space being simply an anti-symmetric copy. The position of the lateral “wall” is modeled by the point where $`u=0`$. Thinking this way points (ii) and (iii) are clearly demonstrated. We proceed now analytically to demonstrate drag reduction (point (ii)) and to understand the profile of $`R`$ (point (iii)).

First we consider the stationary solution of the pure Burger’s equation. Integrate equation (4) in $`y`$ to find

$$\frac{1}{2}u^2=\nu u_y+\frac{1}{2}u_0^2$$
(10)
where the constant of integration was fixed by noticing that for $`L`$ sufficiently large $`u_y`$ is expected to vanish at the boundaries. Multiplying (10) by $`u_y`$, and integrating between $`L`$ and $`L`$ using the boundary conditions, we find the viscous dissipation $`ϵ`$

$$ϵ\nu u_x^2=\frac{2}{3}u_0^3$$
(11)
Next we consider the solution of the uR model for the same boundary conditions (9) and $`R0`$. In the stationary state $`R_t=0`$, and by dividing Eq. (5) by $`Ru`$ we can integrate it formally in $`y`$ and obtain:

$$R=a|u|exp(\frac{1}{\tau }\frac{dy}{u})$$
(12)
where $`a`$ is a constant of integration. This equation is the explanation of point (iii). It says that for small velocity $`u0`$, i.e. for $`y0`$, $`R`$ necessarily goes to $`0`$. In particular, approximating $`u=my`$ near the point $`y=0`$, we obtain $`R|u|^b`$ where $`b=1+1/(m\tau )`$. Thus we should expect that at positions with small $`u`$ where the gradient of $`u`$ is large the generic behavior of $`R`$ is a cusp with $`R=0`$ for $`m\tau 1`$. To compute the dissipation analytically we consider the limit $`\tau \mathrm{}`$, i.e. we look for a solution at the zero order of the perturbation series in $`1/\tau `$. In this limit:

$$R=\frac{R_0}{u_0}|u|.$$
(13)
Returning to Eq. (4) we integrate it in $`y`$ to obtain

$$\nu u_y=\frac{1}{2}u^2SR+sR_0\frac{1}{2}u_0^2.$$
(14)
We can now substitute Eq. (13) in the domain $`Ly0`$ where $`|u|=u`$, and integrate between $`L`$ and $`0`$. Multiplying the result by a factor of 2 we find the viscous dissipation $`ϵ_R`$

$$ϵ_R=\frac{2}{3}u_0^33sR_0u_0$$
(15)
This result is an analytic demonstration of point (ii). We note that our analysis has been performed in the limit $`\tau \mathrm{}`$. For large but finite values of $`\tau `$ the qualitative picture we have drawn is unchanged. Needless to say, the above discussion can be reformulated by keeping constant the energy dissipation while increasing the value of $`u`$ at the boundary, to demonstrate point (i). We choose however to demonstrate point (i) next, using a forced solution. Point (i) is most clearly demonstrated in the uR model using periodic boundary conditions and constant forcing. We consider $`0y2\pi `$ and choose the external forcing $`F`$ to be:

$`F(y)`$ $`=`$ $`f_0sin(4y/3),\text{for}0y{\displaystyle \frac{3\pi }{2}},`$
$`F(y)`$ $`=`$ $`f_1sin(4y),\text{for}{\displaystyle \frac{3\pi }{2}}y2\pi .`$ (16)
We examined the solutions of the uR model for the set of parameters $`f_0=0.1`$, $`f_1=0.05`$, $`\nu =0.01`$, $`s=0.01`$ and initial conditions $`R(y)=2sin(y)`$. The remaining parameter is the relaxation time $`\tau `$. It turns out that for very small values of $`\tau `$, no effect of the polymer field is observed. (We remark that for periodic boundary conditions the limit $`\tau 0`$ corresponds to the case of no polymer.) For $`\tau \mathrm{}`$ no stationary solutions can be obtained. For $`\tau `$ smaller than some critical value $`\tau _c`$, the solution of the uR model shows stable stationary solutions with drag reduction. The typical situation is presented in Fig. 4, showing the numerical solutions for $`\tau =(0.15)^1\tau _c`$, compared against the solution of Burger’s equation. The uR model shows a larger amplitude near the strongest shock due to forcing at $`x=3\pi /4`$. It is worth noting that the gradient is maintained extremely close to the one obtained by the Burger equation, demonstrating nicely point (i).

Point (iii) is nicely demonstrated in Fig. 5 which presents the solution for $`R`$ together with $`u^2`$ for both the uR model and the Burger equation. As one can clearly see, the behavior of $`R`$ is similar to what observed in Fig. 1, namely there is a qualitative similarity between the space dependence of $`R`$ and $`u^2`$, here with sharp cusp in $`R`$ near the point of maximum gradient of $`u`$. On the other hand, the smallest shock present in the solution of Burger’s equation has been completely smoothed out by the uR model. This is an indication that when $`R`$ is not sufficiently suppressed where the gradient of $`u`$ is significant, there can be drag enhancement. This important point will be addressed again in the concluding remarks.

Again, the simplicity of the model affords an analytic explanation of why the solution near the biggest shock shows a larger velocity amplitude compared to the Burger equation. Let $`y_0`$ be the position of the maximal velocity near the shock. The position $`y_0`$ is unchanged in the two models. We can expand $`u`$, $`F`$ and $`R`$ as power series near $`y_0`$. Let $`\mathrm{\Delta }`$ be defined as $`\mathrm{\Delta }=y_sy_0`$, where $`y_s=3\pi /4`$ is where the velocity vanishes (the “wall” position in our model). We have:

$$u(y)=u_0u_2y^2,F(y)=F_0+F_1y,R(y)=R_0+R_1y$$
(17)
where $`F_0f_0\mathrm{\Delta }`$ and $`F_1f_0`$ because of Eq. (16). Inserting (17) into the uR model, we obtain:

$$u_0=\frac{\nu f_0+sR_0\tau ^1}{\mathrm{\Delta }f_0}=u_b(1+P),P=\frac{sR_0}{\nu \tau f_0}$$
(18)
where $`u_b=\nu /\mathrm{\Delta }`$ is the solution obtained without polymer, ($`P=0`$). The increment of the velocity, which is responsible for the drag reduction, is proportional to $`sR_0`$.

In summary, we have introduced a simple model of the effect of polymeric additives to Newtonian fluids, with the aim of understanding in simple mathematical terms some of the prominent features associated with the phenomenon of drag reduction. Needless to say, the model cannot be taken as quantitative; we were concerned with the qualitative features summarized above for convenience as point (i)-(iii). We demonstrated that our model reproduces these qualitative features, and provided straightforward analytic explanation to all those features. It appears that we can draw from the results of this model a few important conclusions: (i) Arguments concerning the turbulent cascade process do not appear essential. These arguments are the hallmark of the theory presented in 90Ge  which proposed that the main effect of the polymer is to introduce a dissipative cutoff at scales larger than the Kolmogorov scale, due to the polymer relaxation time matching there the hydrodynamic time scale. Looking at Eq. (5) one could think that the largest velocity gradients could be estimated as $`|u_y|1/\tau `$. But since $`R`$ can go to zero where the gradient is largest, such estimates cannot be made. Moreover, just an increase in the dissipative scale cannot account for drag reduction; a homogeneous increase in the effective viscosity should lead by itself to drag enhancement rather than reduction. (ii) Drag reduction is a phenomenon that appears on the scales of the system size, involving energy containing modes rather than dissipative, small scale modes 02ACLPP . (iii) The main point appears to be the space dependence of the stretching of the polymer, here modeled by the value of $`R(y)`$. It is crucial that $`R`$ is small where the velocity gradients are large. It is the space dependence of the effective viscosity which should be looked at as the source of drag reduction. A similar conclusion was arrived to in the context of the study of the stability of laminar flows in a channel geometry 01GLP ; 02GLP  accept that there the space dependence of the effective viscosity had been introduced by hand. In the FENE-P context as well as in our model (and presumably in actual experiments) the space dependence appears self consistently. It remains to understand this self consistent build up of differential effective viscosity in the context of the much more elaborate FENE-P model. In light of the present results this appears an extremely worthwhile endeavor that will shed important light on the phenomenon of drag reduction.

###### Acknowledgements.

This work was supported in part by the European Commission under a TMR grant, the German Israeli Foundation, and the Naftali and Anna Backenroth-Bronicki Fund for Research in Chaos and Complexity. We thank E. de Angelis for Figs. 1 and 2 taken from her PhD thesis.


text:
# Forward-Backward Multiplicity Correlations in Symmetric and Asymmetric High Energy Collisions

## Abstract

Forward-backward correlations are explored within the two-component clan model of multiparticle production. It is found that existing data are well described, and, in $`hh`$ collisions, that clans must be allowed to leak particles from one hemisphere to the other. General formulae given for the symmetric case are then extended to the asymmetric one, which is relevant for pA and AB collisions.

The analysis of multiplicity correlations between two different regions of phase-space provides important information on the dynamics of multiparticle production processes.

Due to space limitations, the reader is referred for all details and references to papers  and .

I will at first consider symmetric reactions (like $`e^+e^{}`$ annihilation, $`pp`$ and $`p\overline{p}`$ collisions) with a symmetric choice of the forward (F) and backward (B) regions: in $`e^+e^{}`$ annihilations the F regions is chosen randomly between those defined by a plane through the collision point and perpendicular to the thrust axis; in $`hh`$ collisions the collision axis is chosen as reference. B is always the symmetric region. The forward-backward (FB) correlation strength, $`b_{FB}`$, is defined as usual:

$$b_{FB}=\frac{(n_F\overline{n}_F)(n_B\overline{n}_B)}{\left[(n_F\overline{n}_F)^2(n_B\overline{n}_B)^2\right]^{1/2}},$$
(1)
where $`n_F`$ if the number of particles in F and $`n_B`$ the number of particles in B. A more precise study of the correlations investigates the dependence of the average number of particles in one hemisphere on the number of particles in the opposite one, say $`\overline{n}_F(n_B)`$.

In $`hh`$ collisions the correlation strength grows logarithmically with the c.m. energy $`\sqrt{s}`$ and is rather large (from 0.16 at 63 GeV to 0.7 at 1800 GeV); in $`e^+e^{}`$ annihilation the correlation strength is rather small (0.1 at 91 GeV). Furthermore, $`b_{FB}`$ was measured in 2-jet and 3-jet events separately at LEP1 and was found to be compatible with zero in both cases.

The mechanism of the weighted superposition of different classes of events, which was found to describe the shoulder structure of multiplicity distributions and the oscillations of their $`H_q`$ moments, can be tested on FB correlations. The classes of events are defined according to the reaction considered: in $`hh`$ collisions we define soft events (events without mini-jets) as class 1, and semi-hard events (events with mini-jets) as class 2; in $`e^+e^{}`$ annihilation, class 1 are two-jet events, class 2 are three-jet events.

The starting point is the weighted superposition equation for the joint distribution of $`n_F`$ and $`n_B`$ particles:

$$P(n_F,n_B)=\alpha P_1(n_F,n_B)+(1\alpha )P_2(n_F,n_B),$$
(2)
where $`\alpha `$ is the fraction of events of class 1. Based on the above equation and symmetry properties, it is possible to obtain explicitly a general formula for $`b_{FB}`$, which depends of course on the first two moments of the general multiplicity distribution (MD) in each class $`i=1,2`$ (namely the average multiplicity $`\overline{n}_i`$ and the dispersion $`D_{n,i}`$) and on the correlation strength in each class, $`b_i`$:

$$b_{FB}=\frac{\alpha b_1D_{n,1}^2(1+b_2)+(1\alpha )b_2D_{n,2}^2(1+b_1)+\frac{1}{2}\alpha (1\alpha )(\overline{n}_2\overline{n}_1)^2(1+b_1)(1+b_2)}{\alpha D_{n,1}^2(1+b_2)+(1\alpha )D_{n,2}^2(1+b_1)+\frac{1}{2}\alpha (1\alpha )(\overline{n}_2\overline{n}_1)^2(1+b_1)(1+b_2)}.$$
(3)
We apply first our result to OPAL data, setting $`b_1=b_2=0`$ as the data suggest. We take the values of the moments of the MD in each class from previous fits and obtain $`b_{FB}=0.101`$ to be compared with the OPAL result $`0.103\pm 0.007`$. Thus we conclude that the weighted superposition mechanism describes well FB correlations in $`e^+e^{}`$ annihilation.

We now turn to $`p\overline{p}`$ collisions at 546 GeV (UA5); the experimental result is $`0.58\pm 0.01`$. Taking the shoulder description from a fit to the overall MD, one finds that Eq. (3) with $`b_1=b_2=0`$ gives a too small result, indicating that FB correlations are needed in each substructure. On the other hand, if particles in each substructure were totally uncorrelated, one would get a result which is much larger than experimentally found. We conclude that particles are emitted with partial correlation, as already suggested by UA5.

Since we have succesfully used a two-step production process to describe multiparticle dynamics, leading to the negative binomial (NB) shape of the MD, and which naturally accommodates partial correlation, we apply it to the study of FB correlations in a generalised way:

* In the first step, $`N`$ independent objects (generalised clans) are produced according to a given MD, $`𝒫(N)`$.
* In the second step, particles are produced within each clan according to a MD $`Q(n_c)`$.

This two-step process is used to describe each component (with possibly different parameters in each step): from now on, all our formulae refer to each component separately, but we drop the component index in the formulae. We further assume that clans are produced independently in rapidity thus the distribution in the number of forward clans, $`N_F`$, at fixed total number of clan, $`N`$, is binomial. The successful NB parametrisation of the data is obtained when $`𝒫(N)`$ is the Poisson distribution and $`Q(n_c)`$ is the logarithmic distribution.

Consider now the situation in which particles produced by an F (B) clan all remain in the F (B) region: if clans follow a Poisson MD, as is implied by the very good fit to the data obtained with the NB distribution (NBD), one obtains $`b=0`$, independently of the details of the MD inside clans. It is therefore necessary to allow clans to leak particles from one hemisphere to the other. In order to do this, a new parameter is introduced, $`p`$, which is defined as the average fraction of particles in a clan which remain in the same region and do not leak away to the opposite one. As clans are classified F or B according to where the majority of the particles fall, one has $`0.5p1`$. We also define $`q=1p`$.

The general result is now

$$b=\frac{D_n^2/\overline{n}D_c^2/\overline{n}_c(pq)^2\overline{n}_c+4\gamma /\overline{n}_c}{D_n^2/\overline{n}+D_c^2/\overline{n}_c+(pq)^2\overline{n}_c4\gamma /\overline{n}_c},$$
(4)
where $`\gamma `$ is the covariance between F and B multiplicities within a clan. The above formula simplifies considerably in the case of the NBD, and even further if we assume that particles within a clan fall in F or B independently, since then $`\gamma =(D_c^2\overline{n}_c)pq`$:

$$b=\frac{2\beta pq}{12\beta pq}$$
(5)
where $`\beta `$ is related to the average number of particles per clan: $`\overline{n}_c=\beta \left[(\beta 1)\mathrm{ln}(1\beta )\right]^1`$.

We can now use experimental data to extract the values of the leakage parameter $`p`$ which we do not know, then we use the NBD structure to reproduce the almost linear relation between $`\overline{n}_F`$ and $`n_B`$. Since at 63 GeV the mini-jet sample can be considered negligible, we use only one component and obtain $`p_{\text{soft}}=0.78`$. Because the number of particles per clans in the soft component varies very little from 63 to 900 GeV, we assume now that also the leakage parameter is essentially constant in the GeV region. Using this value for $`p_{\text{soft}}`$ and the measured value of the correlation strength at 900 GeV, we can now calculate $`p_{\text{semi-hard}}=0.77`$. The just mentioned values can then be used to predict the behaviour of the correlation strength with energy: e.g., it can be shown that it will not continue to grow linearly with $`\mathrm{ln}\sqrt{s}`$ in the multi-TeV region if the leakage parameter $`p_{\text{semi-hard}}`$ remains constant .

Within the same two-step process framework, we can extend our formulae to the case in which either the reaction is asymmetric (e.g. consider proton-nucleus or nucleus A-nucleus B with A $``$ B) or the F and B regions are asymmetrically defined, or both cases are realised . In order to proceed, we must allow for the possibility that the average leakage from F to B may be different from the average leakage from B to F, and thus introduce different parameters: $`p_Fp_B`$; we must also allow for the possibility that, on average, more clans are produced in one region than in the other; we call $`r`$ the average fraction of clans produced in the F region: in general $`r1/2`$. The symmetric case is recovered for $`p_F=p_B`$ and $`r=1/2`$.

Using the formalism of generating functions, we arrive at formulae for the joint MD in the general case. The joint generating function (GF) for each component in Eq. (2) is

$$G(z_F,z_B)=\underset{n_F,n_B}{}z_F^{n_F}z_B^{n_B}P(n_F,n_B).$$
(6)
Since clans are Poisson distributed and independent of each other, the joint GF for $`N_F`$ forward clans and $`N_B`$ backward clans is

$$𝒢(z_F,z_B)=\mathrm{exp}\left\{\overline{N}[rz_F+(1r)z_B1]\right\},$$
(7)
where $`\overline{N}`$ is the overall average number of clans. Calling $`g_{c,I}(z_F,z_B)`$ the GF for the joint MD of particles from one forward ($`I=F`$) or one backward ($`I=B`$) clan, we make use again of the independence among clans and obtain the general formula:

$$G(z_F,z_B)=𝒢(g_{c,F}(z_F,z_B),g_{c,B}(z_F,z_B)).$$
(8)
The special case in which the MD is of NB type is of particular interest. In that case in fact one has simply (we assume that particles within a clan fall in F or B independently)

$$g_{c,F}(z_F,z_B)=g_{\text{log}}(z_Fp_F+z_Bq_F)\text{and}g_{c,B}(z_F,z_B)=g_{\text{log}}(z_Fq_B+z_Bp_B),$$
(9)
where $`g_{\text{log}}`$ is the GF for the logarithmic distribution of parameter $`\beta `$:

$$g_{\text{log}}(z)\mathrm{log}(1z\beta )/\mathrm{log}(1\beta ).$$
(10)
The generating function for the joint probability of $`n_F`$ and $`n_B`$ particles is then given by

$$G(z_F,z_B)=\mathrm{exp}\left\{r\overline{N}[g_{\text{log}}(z_Fp_F+z_Bq_F)1]\right\}\mathrm{exp}\left\{(1r)\overline{N}[g_{\text{log}}(z_Fq_B+z_Bp_B)1]\right\}.$$
(11)
This formalism can be used to describe FB multiplicity correlations in very general terms. It is interesting to point out a few properties:

* The marginal distribution is the convolution of two NBD’s, thus it can be shown that it is an infinitely divisible distribution.
* When $`r=1/2`$, one obtain for the GF of the marginal distribution a simple correction to the NBD generating function; e.g., for the forward marginal GF:

$$g_F(z)G(z,1)=\left[1+\frac{\overline{n}}{k}(1z)+\frac{\overline{n}^2}{k^2}p_Fq_B(1z)\right]^{k/2}.$$
(12)
* $`\overline{n}_F(n_B)`$ is not a linear function of $`n_B`$, and vice-versa, except for the particular cases $`p_i=1/2`$ and $`p_i=1`$, ($`i`$=F,B); examples are shown in .

In conclusion, general formulae for the FB multiplicity correlation strength and F and B joint distributions generating functions have been given in the framework of the superposition mechanism of two weighted MD’s for different classes of events. Assuming NB regularity behaviour for 2- and 3-jet samples of events in $`e^+e^{}`$ annihilation at LEP energy, results obtained by OPAL collaborations are correctly reproduced within experimental errors. In $`pp`$ collisions the FB multiplicity correlation strengths in the two substructures (soft and semi-hard events) turn out to be quite important: in particular an interesting connection is found between the particle populations within clans, particle leakage from clans in one hemisphere to the opposite hemisphere and superposition effect between different substructure of the collision. This finding favours structures with larger particle populations per clans and the decrease of the average number of clans.

It has been shown that assuming different particle leakage percentages ($`p_Bp_F`$) for binomially generated particles from clans in one hemisphere to the opposite one and asymmetric ($`r1/2`$) distribution in the two hemispheres of binomially generated clans, a general formula for the generating function of the joint $`(n_F,n_B)`$-charged particle multiplicity distribution for each class of events (or substructure) in asymmetric reactions can be obtained when the total multiplicity distribution is of negative binomial type, and can be easily generalised to any discrete infinitely divisible multiplicity distribution. This search is relevant for the study of forward-backward multiplicity correlations in non-identical heavy ion collisions and in proton-nucleus collisions. Accordingly, the newly introduced particle leakage and asymmetry parameters can be used for classifying different classes of collisions.