text:
# GENERALIZED UNCERTAINTY PRINCIPLE AND DARK MATTER

## I Introduction

It is by now widely accepted that dark matter (DM) constitutes a substantial fraction of the present critical energy density in the universe. However, the nature of DM remains an open problem. There exist many DM candidates, among which a contending category is weakly interacting massive particles, or WIMPs. It has been suggested that primordial black holes (PBHs) are a natural candidate for WIMPs. More recent studies based on the PBH production from the “blue spectrum” of inflation demand that the spectral index $`n1.3`$, but this possibility may be ruled out by the recent WMAP experiment.

In the standard view of black hole thermodynamics, based on the entropy expression of Bekenstein and the temperature expression of Hawking, a small black hole should emit blackbody radiation, thereby becoming lighter and hotter, leading to an explosive end when the mass approaches zero. However Hawking’s calculation assumes a classical background metric and ignores the radiation reaction, assumptions which must break down as the black hole becomes very small and light. Thus it does not provide an answer as to whether a small black hole should evaporate entirely, or leave something else behind, which we refer to as a black hole remnant (BHR).

Numerous calculations of black hole radiation properties have been made from different points of view, and some hint at the existence of remnants, but none appears to give a definitive answer. A cogent argument against the existence of BHRs can be made: since there is no evident symmetry or quantum number preventing it, a black hole should radiate entirely away to photons and other ordinary stable particles and vacuum, just like any unstable quantum system.

In a series of recent papers, a generalized uncertainty principle (GUP) was invoked to argue the contrary, that the total collapse of a black hole may be prevented by dynamics and not by symmetry, just like the prevention of hydrogen atom from collapse by the uncertainty principle. These arguments then lead to a modified black hole entropy and temperature, and as a consequence the existence of a BHR at around the Planck mass. This notion was then combined with hybrid inflation model and it was shown that primordial BHRs might in principle be the primary source for dark matte. In this paper we briefly reproduce these arguments, and include additional discussion on the reheating temperature as a result of black hole evaporation.

## II Generalized Uncertainty Principle

As a result of string theory or general considerations of quantum mechanics and gravity, the GUP gives the position uncertainty as

$$\mathrm{\Delta }x\frac{\mathrm{}}{\mathrm{\Delta }p}+l_p^2\frac{\mathrm{\Delta }p}{\mathrm{}},$$
(1)
where $`l_p=(G\mathrm{}/c^3)^{1/2}1.6\times 10^{33}`$cm is the Planck length. A heuristic derivation may also be made on dimensional grounds. We think of a particle such as an electron being observed by means of a photon with momentum $`p`$. The usual Heisenberg argument leads to an electron position uncertainty given by the first term in Eq.(1). But we should add to this a term due to the gravitational interaction of the electron with the photon, and that term must be proportional to $`G`$ times the photon energy, or $`Gpc`$. Since the electron momentum uncertainty $`\mathrm{\Delta }p`$ will be of order of $`p`$, we see that on dimensional grounds the extra term must be of order $`G\mathrm{\Delta }p/c^3`$, as given in Eq.(1). Note that there is no $`\mathrm{}`$ in the extra term when expressed in this way. The position uncertainty has a minimum value of $`\mathrm{\Delta }x=2l_p`$, so the Planck distance, $`l_p`$, plays the role of a fundamental length.

## III Black Hole Remnant

The characteristic energy $`E`$ of the emitted photons may be estimated from the uncertainty principle. In the vicinity of the black hole surface there is an intrinsic uncertainty in the position of any particle of about the Schwarzschild radius, $`\mathrm{\Delta }xr_s`$, due to the behavior of its field lines \- as well as on dimensional grounds. This leads to a momentum uncertainty

$$\mathrm{\Delta }p\frac{\mathrm{}}{\mathrm{\Delta }x}=\frac{\mathrm{}}{r_s}=\frac{\mathrm{}c^2}{2GM_{\mathrm{BH}}},$$
(2)
and hence to an energy uncertainty of $`\mathrm{\Delta }pc\mathrm{}c^3/2GM_{\mathrm{BH}}`$. We identify this as the characteristic energy of the emitted photon, and thus as a characteristic temperature; it agrees with the Hawking temperature up to a factor $`4\pi `$, which we will henceforth include as a “calibration factor” and write (with $`k_B=1`$),

$$T_\mathrm{H}\frac{\mathrm{}c^3}{8\pi GM_{\mathrm{BH}}}=\frac{M_p^2c^2}{8\pi M_{\mathrm{BH}}},$$
(3)
where $`M_p=(\mathrm{}c/G)^{1/2}1.2\times 10^{19}`$GeV is the Planck mass.

The blackbody energy output rate of BH is given by

$$\dot{x}=\frac{1}{t_{ch}(x_i^33t/t_{ch})^{2/3}},$$
(4)
where $`x=M_{\mathrm{BH}}/M_p`$ and $`x_i`$ refers to the initial mass of the hole. $`t_{ch}=60(16)^2\pi t_p4.8\times 10^4t_p`$ is a characteristic time for BH evaporation, and $`t_p=(\mathrm{}G/c^5)^{1/2}0.54\times 10^{43}`$sec is the Planck time. The black hole thus evaporates to zero mass in a time given by $`t/t_{ch}=x_i^3/3`$, and the rate of radiation has an infinite spike at the end of the process.

The momentum uncertainty according to the GUP is

$$\frac{\mathrm{\Delta }p}{\mathrm{}}\frac{\mathrm{\Delta }x}{2l_p^2}\left[1\sqrt{14l_p^2/(\mathrm{\Delta }x)^2}\right].$$
(5)
Therefore the modified black hole temperature becomes

$$T_{\mathrm{GUP}}=\frac{M_pc^2}{4\pi }x\left[1\sqrt{11/x^2}\right].$$
(6)
This agrees with the Hawking result for large mass if the negative sign is chosen, whereas the positive sign has no evident physical meaning. Note that the temperature becomes complex and unphysical for mass less than the Planck mass and Schwarzschild radius less than $`2l_p`$. At the Planck mass the slope is infinite, which corresponds to zero heat capacity of the black hole, and the evaporation comes to a stop.

If there are $`g`$ species of relativistic particles, then the BH evaporation rate is

$$\dot{x}=\frac{16g}{t_{ch}}x^6\left[1\sqrt{11/x^2}\right]^4.$$
(7)
Thus the hole with an initial mass $`x_i`$ evaporates to a Planck mass remnant in a time given by

$`\tau `$ $`=`$ $`{\displaystyle \frac{t_{ch}}{16g}}\left[{\displaystyle \frac{8}{3}}x_i^38x_i{\displaystyle \frac{1}{x_i}}+{\displaystyle \frac{8}{3}}(x_i^21)^{3/2}4\sqrt{x_i^21}+4\mathrm{cos}^1{\displaystyle \frac{1}{x_i}}+{\displaystyle \frac{19}{3}}\right]`$ (8)
$``$ $`{\displaystyle \frac{x_i^3}{3g}}t_{ch},x_i1.`$ (9)
The energy output given by Eq.(7) is finite at the end point where $`x=1`$, i.e., $`dx/dt|_{x=1}=16g/t_{ch}`$, whereas for the Hawking case it is infinite at the endpoint where $`x=0`$. The present result thus appears to be more physically reasonable. The evaporation time in the $`x_i1`$ limit agrees with the standard Hawking picture.

## IV Hybrid Inflation and Black Hole Production

The hybrid inflation, first proposed by A. Linde, can naturally induce large number of small PBHs. In the hybrid inflation model two inflaton fields, $`(\varphi ,\psi )`$, are invoked. Governed by the inflation potential, $`\varphi `$ first executes a “slow-roll” down the potential, and is responsible for the more than 60 e-folds expansion while $`\psi `$ remains zero. When $`\varphi `$ eventually reduces to a critical value, it triggers a phase transition that results in a “rapid-fall” of the energy density of the $`\psi `$ field, which lasts only for a few e-folds, that ends the inflation.

The evolution of the $`\psi `$ field during the second stage inflation, measured backward from the end, is

$$\psi (N[t])=\psi _e\mathrm{exp}(sN[t]),$$
(10)
where $`N(t)=H_{}(t_et)`$ is the number of e-folds from $`t`$ to $`t_e`$, $`H_{}`$ is the Hubble parameter during inflation, and $`s`$ is a numerical factor of the order unity.

Quantum fluctuations of $`\psi `$ induce variations of the starting time of the second stage inflation, i.e., $`\delta t=\delta \psi /\dot{\psi }`$. This translates into perturbations on the number of e-folds, $`\delta N=H_{}\delta \psi /\dot{\psi }`$, and therefore the curvature contrasts, $`\delta \rho /\rho \delta `$. With an initial density contrast $`\delta (m)\delta \rho /\rho |_m`$, the probability that a region of mass $`m`$ becomes a PBH is

$$P(m)\delta (m)e^{w^2/2\delta ^2}.$$
(11)
Let us assume that the universe had inflated $`e^{N_c}`$ times during the second stage of inflation. It can be shown that

$$e^{N_c}\left(\frac{2M_p}{sH_{}}\right)^{1/s},$$
(12)
and the curvature perturbations reentered the horizon at time

$$tt_h=H_{}^1e^{3N_c}.$$
(13)
At this time if the density contrast was $`\delta 1`$, then BHs with size $`r_sH_{}^1e^{3N_c}`$ would form with an initial mass

$$M_{\mathrm{BH}i}\frac{M_p^2}{H_{}}e^{3N_c}.$$
(14)
Following the numerical example given in Ref.21, we let $`H_{}5\times 10^{13}`$ GeV and $`s3`$. Then the density contrast can be shown to be $`\delta 1/7`$, and the fraction of matter in the BH is thus $`P(m)10^2`$. From Eq.(12), $`e^{N_c}54`$. So the total number of e-folds is $`N_c4`$. The black holes were produced at the moment $`t_h2\times 10^{33}`$ sec, and had a typical mass $`M_{\mathrm{BH}i}4\times 10^{10}M_p`$. Let $`g100`$. Then the time it took for the BHs to reduce to remnants, according to Eq.(9), is

$$\tau \frac{x_i^3}{3g}t_{ch}5\times 10^{10}\mathrm{sec}.$$
(15)
The “black hole epoch” thus ended in time for baryogenesis and other subsequent epochs in the standard cosmology. As suggested in Ref.21, such a post-inflation PBH evaporation provides an interesting mechanism for reheating.

## V Black Hole Remnants as Dark Matter

This process also provides a natural way to create cold dark matter. Although in our example $`P(m)10^2`$, PBHs would soon dominate the energy density by the time $`tP(m)^2t_h2\times 10^{29}`$s, because the original relativistic particles would be diluted much faster than non-relativistic PBHs. By the time $`t\tau `$, all the initial BH mass ($`x_i`$) had turned into radiation except one unit of $`M_p`$ preserved by each BHR. As BH evaporation rate rises sharply towards the end, the universe at $`t\tau `$ was dominated by the BH evaporated radiation.

Roughly, $`\mathrm{\Omega }_{\mathrm{BHR},\tau }1/x_i`$ and $`\mathrm{\Omega }_{\gamma ,\tau }1`$ at $`t\tau `$, and since the universe resumed its standard evolution after the black hole epoch ($`t>\tau `$), we find the density parameter for the BHR at present to be

$$\mathrm{\Omega }_{\mathrm{BHR},0}\left(\frac{t_{eq}}{\tau }\right)^{1/2}\left(\frac{t_0}{t_{eq}}\right)^{2/3}\frac{1}{x_i}\mathrm{\Omega }_{\gamma ,0},$$
(16)
where $`t_04\times 10^{17}`$s is the present time, and $`t_{eq}`$ is the time when the density contributions from radiation and matter were equal. It is clear from our construction that $`(t_{eq}/\tau )^{1/2}x_i`$. So $`t_{eq}10^{12}`$ sec, which is close to what the standard cosmology assumes, and Eq.(16) is reduced to a simple and interesting relationship:

$$\mathrm{\Omega }_{\mathrm{BHR},0}\left(\frac{t_0}{t_{eq}}\right)^{2/3}\mathrm{\Omega }_{\gamma ,0}10^4\mathrm{\Omega }_{\gamma ,0}.$$
(17)
In the present epoch, $`\mathrm{\Omega }_{\gamma ,0}10^4`$. So we find $`\mathrm{\Omega }_{\mathrm{BHR},0}𝒪(1)`$, about the right amount for dark matter!

## VI Black Hole Epoch and Reheating Temperature

As discussed above, shortly after PBHs were produced the density of the universe was dominated by the BHs. Eventually the universe was reheated through their continuous evaporation. To simplify the discussion we ignore BH accretions of the radiation as well as BH mergers. Then under Hubble expansion the effective reheating temperature at the end of the black hole epoch, or $`t\tau `$, can be expressed as

$$T_r(\tau [x_i])=\frac{1}{x_i1}_1^{x_i}𝑑xT_{\mathrm{GUP}}(x)\frac{a(t[x])}{a(\tau [x_i])},$$
(18)
where $`a(t)`$ is the scale factor. Since $`x_i1`$, the evaporation only became effective near the late times during this black hole epoch, when the energy density was dominated by the BH radiation. As a further approximation we assume radiation dominance throughout the BH epoch so that $`a(t)t^{1/2}`$. Expressing $`t`$ in terms of $`x`$ using Eq.(7), we find

$$T_r(\tau [x_i])\frac{M_pc^2}{16\pi (x_i1)}\left[2\mathrm{log}(2x_i)1\right]+𝒪\left(\frac{1}{x_i^3}\right).$$
(19)
In our model $`x_i4\times 10^{10}`$. So $`T_r(\tau [x_i])1.3\times 10^8`$ GeV, which is sufficiently lower than the Planck and the GUT scales, but higher than the baryogenesis scale.

## VII Acknowledgements

I deeply appreciate my early collaborations and fruitful discussions with Ronald J. Adler. I also thank S. Dimopoulos, A. Green, A. Linde, M. Shmakova, and K. Thompson for helpful discussions. This work is supported by the Department of Energy under Contract No. DE-AC03-76SF00515.


format: text

text:
# A note on non-thermodynamical applications of non-extensive statistics

## Abstract

It is pointed out that the constraint to be imposed to the maximization of the entropy for processes outside the class of thermodynamical systems, is generally not well defined. In fact, any probability distribution can be derived from Jaynes’s principle with a suitable choice of the constraint. In the case of Tsallis’s non-extensive formalism, this implies that it is not possible to establish any connection between specific non-thermodynamical processes and non-extensive mechanisms and, in particular, to assign any unambiguous non-extensivity index $`q`$ to those processes.

Boltzmann-Gibbs thermostatistics bridges the microscopic description of physical systems that obey the laws of mechanics with the macroscopic picture drawn from the principles of thermodynamics. From a historical perspective, it reconciled the physics of thermal processes, fully congenial with our everyday-life experience, and the mechanistic interpretation of the Universe as a huge ensemble of interacting particles –two views that, in the middle of the nineteenth century, were far from being perceived as compatible, both mathematically and philosophically .

Physical systems whose microscopic degrees of freedom are governed by the equations of (classical or quantum) mechanics, and which have reached a state of macroscopic equilibrium, will hereafter be called thermodynamical systems. The Boltzmann-Gibbs formulation makes it possible to derive a theory for the collective equilibrium state variables of a thermodynamical system from its microscopic dynamics. At the same time, it gives origin to the mesoscopic (statistical) description level. An overwhelming corpus of experimental work validates the results of this procedure, both at mesoscopic and macroscopic levels. It is known, however, that Boltzmann-Gibbs thermostatistics fails to give a mathematically consistent description of certain physical systems –notably, those driven by gravitational and other long-range interactions– since the predicted values of some of their state variables diverge . Though it is not clear whether such systems attain at all a state identifiable with macroscopic equilibrium , their thermodynamical properties would be characterized by non-extensive macroscopic variables.

In 1988, C. Tsallis proposed an extension of Boltzmann-Gibbs equilibrium thermostatistics based on a variational principle for a generalized form of the entropy ,

$$S_q=\frac{_ip_i^q1}{q1},$$
(1)
where $`p_i`$ is the probability that the thermodynamical system under study is found in its $`i`$-th quantum state and $`q(\mathrm{},\mathrm{})`$. The fact that, for $`q1`$, the entropy $`S_q`$ is not additive with respect to the factorization of the probabilities, has led to the assumption that this formalism may provide consistent thermostatistical description of non-extensive systems. The parameter $`q`$ has been called non-extensivity index, since it measures the deviation from additivity of $`S_q`$.

In the canonical scenario, the generalized entropy (1) is maximized with respect to $`p_i`$, with the constraint of probability normalization,

$$\underset{i}{}p_i=1,$$
(2)
and fixing the value of a generalized form of the mean energy ,

$$\frac{_ip_i^qϵ_i}{_ip_i^q}=E_q,$$
(3)
where $`ϵ_i`$ is the energy of the $`i`$-th state. This canonical maximization procedure yields the energy-dependent probability distribution

$$p_i=Z^1[1+\beta (q1)ϵ_i]^{1/(q1)},$$
(4)
where $`Z`$ is a normalization constant analogous to the partition function and $`\beta `$ is an auxiliary parameter. The variational formulation of canonical Boltzmann-Gibbs thermostatistics is fully recovered in the limit $`q1`$, where $`\beta `$ reduces to the inverse temperature.

A remarkable property of Tsallis’s generalization is that it preserves the mathematical structure of standard thermostatistics for any value of $`q`$ . This noticeable feature justifies the rather unexpected form of the constraint (3), which replaces the usual definition of the mean energy $`E=_ip_iϵ_i`$. In a long series of publications , it has been shown that most of the theorems of equilibrium statistical mechanics, as well as many results concerning linear and nonlinear non-equilibrium properties, can be formally generalized in the frame of the extended formalism. Apart from this formal equivalence with Boltzmann-Gibbs thermostatistics, the relevance of the non-extensive formulation should be validated by the observation of actual thermodynamical systems with an energy distribution of the form (4). At the same time, the role of the quantity $`E_q`$ as a macroscopic property of the system in question should be assessed, in comparison with the role of the mean energy $`E`$ as a state variable of extensive systems. Until now, however, there is no conclusive evidence that non-extensive thermostatistics might correctly describe any thermodynamical system .

On the other hand, many real systems have been identified where the statistical distribution of their relevant variables –not of the energy, however– are well fitted with functions of the type of Eq. (4) . To be specific, empirical distributions for quantities $`x`$ defined over the semi-infinite range $`(0,\mathrm{})`$ have been systematically fitted with the two-parameter function

$$p(x)=N(1+ax)^b,$$
(5)
\[cf. Eq. (4)\] where $`N(a,b)`$ is chosen in such a way that $`p(x)`$ is normalized to unity. Through identification of $`p(x)`$ with the distribution of Eq. (4), the fitting parameters $`a`$ and $`b`$ are used to assign an “inverse temperature” $`\beta `$ and a non-extensivity index $`q`$ to the empirical distribution under consideration. For quantities $`x(\mathrm{},\mathrm{})`$ the chosen function is, instead,

$$p(x)=N(1+ax^2)^b.$$
(6)
This approach has been applied, for instance, to momentum distributions in elementary particle interactions , velocity distributions in diffusing biological systems , and volume and return distributions in financial processes . In all the reported cases the result of the fitting seems to be quite good, a circumstance that has invariably led to the claim that the systems under study are governed by mechanisms characterized by non-extensivity. Frequently, moreover, connections have been established with presumably related concepts, such as self-similarity, scale invariance, non-ergodicity, meta- and quasi-equilibria, criticality, algorithmic complexity, et cætera .

Now, the fact that an empirical distribution is well fitted by a function derived from a variational principle analogous to that of canonical thermostatistics, does not necessarily mean that the nature of the underlying processes is the same as in a thermodynamical system, as described by equilibrium statistical mechanics, nor that the state of the system can be identified with canonical thermal equilibrium. The existence of a mechanical Hamiltonian formulation for processes such as, say, elementary interactions may be a matter of controversy, but a Hamiltonian-like realistic description of a biological population or the stock market should be out of question. However obvious, this remark raises a significant question associated with the maximization procedure that yields fitting functions such as those of Eqs. (5) and (6), and with Jaynes’s principle in general: Besides probability normalization, which is the “correct” constraint to be used in the canonical maximization of the entropy?

For extensive thermodynamical systems, numberless instances of experimental validation show that the mean energy $`E`$ is to be fixed. The formal equivalence of non-extensive themostatistics with the Boltzmann-Gibbs formulation suggests in turn that constraint (3) may be necessary to deal with non-extensive thermodynamical systems (see, however, Ref. ). On the other hand, no rigorous justification can generally hold for any of the infinitely many constraints that can be imposed in the case of non-thermodynamical systems . Assuming that the states of the system are well defined, any function $`\varphi (x)`$ of the relevant variable $`x`$ may be used to introduce the average

$$\mathrm{\Phi }_q=\frac{_ip_i^q\varphi (x_i)}{_ip_i^q},$$
(7)
in full analogy with Eq. (3). Maximization of the entropy (1) under constraints (2) and (7) leads to

$$p(x)=Z^1[1+\beta (q1)\varphi (x)]^{1/(q1)}.$$
(8)
This freedom in the choice of the variational constraint seems to have been systematically overlooked by those authors who applied non-extensive thermostatistics to the description of empirical data from non-thermodynamical systems. In fact, choosing fitting functions as in Eqs. (5) and (6) amounts to restricting $`\varphi (x)`$ to $`x`$ and $`x^2`$, respectively.

An exception to this rule, which dramatically points out the necessity of considering more general constraints when studying non-thermodynamical systems, is given by the fittings of velocity-difference distributions in turbulent flows . In this case, reasonable fittings with functions as in Eq. (8) are obtained only when $`\varphi (x)`$ is allowed to take rather complicated forms, typically, $`\varphi (x)=x^{2\alpha }/2c\text{sgn}(x)(|x|^\alpha |x|^{3\alpha }/3)`$ with $`c0`$ and $`0<\alpha <1`$. The circumstance that such form of $`p(x)`$ implies an unusual choice of the function being averaged in the variational constraint does not seem to have been discussed in the relevant literature, though.

The arbitrariness of the variational constraint for non-thermodynamical systems makes it possible to finely tune the fitting function by a suitable choice of $`\varphi (x)`$. As an illustration, we show in Fig. 1 an empirical distribution of high-frequency stock returns in the New York Stock Exchange , along with two fittings. The dashed curve corresponds to the fitting function of Eq. (6), with $`b=2.5`$ ($`q=1.4`$) . While the overall quality of the approximation is good, a systematic deviation from the empirical data is apparent for intermediate values of the distribution ($`10^4p(x)10^2`$). This deviation is considerably reduced if, as shown by the full curve, the maximization of entropy is subject to constraint (7) with $`\varphi (x)=|x|^{1.6}`$, which gives $`q=1.3`$. As discussed above, this constraint is as valid as any other, and has the advantage of yielding a better fitting for the empirical data.

It immediately results from Eq. (8) that introducing the function

$$\varphi (x)=\frac{1}{q1}[Ap(x)^{1q}+\varphi _0]$$
(9)
in constraint (7) leads to a variational principle which exactly yields any given distribution $`p(x)`$. The constants $`A`$ and $`\varphi _0`$ fix the origin and units of measure for the average $`\mathrm{\Phi }_q`$, while $`q`$ establishes the connection between the shapes of $`\varphi (x)`$ and $`p(x)`$. Note that the index $`q`$ can be chosen arbitrarily, and that a different form of $`\varphi (x)`$ is obtained for each value of $`q`$. For $`q1`$, we get $`\varphi (x)=A^{}\mathrm{ln}p(x)+\varphi _0^{}`$.

The simple observation that any distribution can be derived from a variational principle for the entropy if a suitable constraint is chosen, has far-reaching consequences when interpreting the fitting of non-thermodynamical empirical data in the frame of non-extensive statistics. In particular, it voids of meaning any claim of connection between the fitted data and possible non-extensive mechanisms underlying the system in question. A quantitative proof of this assertion is provided by the fact that, for a given system, the non-extensivity index $`q`$ is not uniquely defined, and can in fact be given any value by an appropriate choice of $`\varphi (x)`$. Any system, in fact, could be made “extensive” by simply using the constraint that yields $`q=1`$!

At this point, it could be argued that the kind of constraints arising from a choice of $`\varphi (x)`$ as in Eq. (9) will generally be unconventional, typically involving complicated functions of the relevant variables. We have seen that this is in fact the case if the distributions of velocity-differences in turbulent flows are forced to fit non-extensive thermostatistics. For other physical systems, it has already been remarked that insisting to stick to a variational principle may necessarily lead to consider “non-traditional” constraints . For non-physical systems, unfortunately, we can hardly discern between “traditional” and “non-traditional” constraints. Without a rigorous argument to decide on this point, no choice can possibly be dismissed on such basis.

In summary, we have pointed out that any probability distribution can be derived from the entropic variational formalism that underlies Tsallis’s non-extensive thermostatistics, if a suitable “state function” is chosen to define a constraint for the maximization procedure. The same remark should hold for any variational formalism based on Jaynes’s principle. While for macroscopic Hamiltonian systems in thermodynamical equilibrium it is well established that the mean energy is to be fixed, for non-thermodynamical systems it is generally not possible to argue for or against any choice. Therefore, any claim of connection between non-thermodynamical processes and non-extensive mechanisms, based on the fitting of empirical probability distributions with the functions derived from Tsallis’s variational formalism, is essentially insubstantial. Such kind of fittings provide at most a phenomenological description of the systems in question, and bear little information on their true nature.


text:
# Phonon thermal conductivity in doped 𝐋𝐚_𝟐𝐂𝐮𝐎_𝟒: relevant scattering mechanisms

## I Introduction

The segregation of spins and holes into stripelike arrangements appears to be a common feature of doped Mott insulators.Tranquada et al. (1994); Tranquada et al. (1995a, b); Lee and Cheong (1997); Mori et al. (1998) These so-called stripe correlations seem to be of particular importance in understanding the electronic phase diagram of high-temperature superconductors, where a competition between a static stripe phase and the superconducting phase is widely discussed.White and Scalapino (1999); Chakravarty et al. (2001); Zaanen (2003) Such a competition is expected to be reflected by the dynamics of stripes, i.e., static stripes should reduce the superconducting order parameter while stripes should be fluctuating in the presence of fully developed superconductivity. Indeed, there is growing experimental evidence for such a scenario.Büchner et al. (1994); Tranquada et al. (1995a); Klauss et al. (2000); Lake et al. (2001, 2002) While many experiments give evidence towards static stripes of holes and spins,Tranquada et al. (1995a); Tranquada et al. (1997, 1999) signatures of stripe fluctuations currently only comprise magnetic correlations.Cheong et al. (1991); Yamada et al. (1998); Mook et al. (1998) A direct observation of fluctuating charge stripes, however, is still missing. A promising alternative approach to study stripe dynamics involves the phonon heat transport which is an indirect probing method. Since charge stripes lead to lattice distortions, a sensitivity of the phonon heat transport to the dynamics as well as the degree of periodicity of stripes can be expected.

The thermal conductivity $`\kappa `$ of doped $`\mathrm{La}_2\mathrm{CuO}_4`$ has repeatedly been the subject of experimental research, yet no detailed understanding of the rather complicated and strong changes of the temperature dependence of $`\kappa `$ upon partial substituting Sr and/or small Rare Earths (RE) like Nd or Eu for La has been achieved. For example, in the antiferromagnetic insulators $`\mathrm{La}_{2\mathrm{y}}(\mathrm{RE})_\mathrm{y}\mathrm{CuO}_4`$ a huge peak at room temperature which arises due to magnetic heat transport is found in the thermal conductivity parallel to the $`\mathrm{CuO}_2`$-planes ($`\kappa _{ab}`$), while the thermal conductivity perpendicular to the $`\mathrm{CuO}_2`$-planes ($`\kappa _c`$) is purely phononic without a high temperature peak.Hess et al. (2003) Similarly intriguing is a strong suppression of the phononic low temperature peak in both $`\kappa _{ab}`$ and $`\kappa _c`$, which is found in the superconducting doping levels of $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$.Nakamura et al. (1991) There is a seeming correlation between this suppression and superconductivity because a phononic low temperature peak reappears in overdoped, non-superconducting $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ as well as in $`\mathrm{La}_{2\mathrm{x}\mathrm{y}}(\mathrm{RE})_\mathrm{y}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$, where superconductivity is suppressed in favor of static stripe order.Büchner et al. (1994); Baberski et al. (1998); Klauss et al. (2000) Baberski et al. Baberski et al. (1998) qualitatively explained these observations based on the idea that in superconducting compounds fluctuating stripes provide a new scattering channel for phonons. In recent studies by Sun et al., such a scattering channel plays an important role in the data interpretation.Sun et al. (2003a, b, b)

There is clear-cut evidence that in isostructural stripe ordering $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{NiO}_4`$ the phonon thermal conductivity $`\kappa _{\mathrm{ph}}`$ is closely correlated with both the dynamics and the periodicity of stripes: While $`\kappa _{\mathrm{ph}}`$ is almost unaffected in the presence of static and long range ordered stripes, it is strongly suppressed as soon as the stripes become disordered or dynamic.Hess et al. (1999); Cassel et al. (1999); Lee et al. (2002) Apparently, in these compounds the thermal conductivity is indeed a probe for stripe correlations. One might question, however, whether this is true also in the cuprates for two reasons: First, the electron-phonon coupling in the nickelates is much stronger than in the cuprates.Eisaki et al. (1992); Anisimov et al. (1992); Bi and Eklund (1993); Chen et al. (1993); Zaanen and Littlewood (1994); Cheong et al. (1994) Therefore, the effect of stripes on $`\kappa _{\mathrm{ph}}`$ can be expected to be much smaller in the cuprates. Second, from a structural point of view the situation in the nickelates is much simpler than in the cuprates. In the latter a structural instability exists and as a consequence a number of structural phase transitions occur as a function of temperature as well as of Sr- and RE-content.Büchner et al. (1994) Since the structural instability involves soft phonon modes,Birgeneau et al. (1987); Pintschovius (1990); Martínez et al. (1991); Keimer et al. (1993) enhanced scattering of the heat carrying phonons is likely.

In this article we reinvestigate the phonon thermal conductivity $`\kappa _{\mathrm{ph}}`$ of doped $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ and present new experimental results on Eu-doped single crystals. The single crystalline data allow us to investigate the anisotropy of the $`\kappa `$-tensor and provide more precise absolute values of $`\kappa `$. In previous measurements on polycrystalsSera et al. (1997); Baberski et al. (1998) this information was not available. It is, however, necessary in order to judge the strength and therefore the importance of various scattering mechanisms for phonons. The analysis of our data yields compelling arguments that both impurities (dopants) and soft phonons, which are associated with the lattice instability in these compounds, strongly scatter phonons and therefore must not be neglected in the data interpretation. Data on $`\mathrm{La}_2\mathrm{NiO}_4`$ corroborate this conclusion and allow us to qualitatively understand $`\kappa _{\mathrm{ph}}`$ of $`\mathrm{La}_{2\mathrm{x}\mathrm{y}}(\mathrm{RE})_\mathrm{y}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ in a wide doping range. In particular, it is not necessary to include a stripe induced scattering channel.

The structure of this paper is as follows: After a brief description of the experimental details in section II, we review previous results on the thermal conductivity of doped $`\mathrm{La}_2\mathrm{CuO}_4`$ in section III, before we proceed to the presentation of our new experimental results on $`\mathrm{La}_{1.8\mathrm{x}}\mathrm{Eu}_{0.2}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ in section IV. Our main results will be discussed in section V.

## II Experimental

We have prepared single crystals of $`\mathrm{La}_{1.8\mathrm{x}}\mathrm{Eu}_{0.2}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ ($`x=0`$, 0.08, 0.15, 0.2) as well as of $`\mathrm{La}_2\mathrm{NiO}_4`$ utilizing the traveling solvent floating zone technique. $`\kappa `$ of these crystals was measured as a function of temperature $`T`$. We used a standard steady state method on pieces cut along the principal axes with a typical length of 2$`\mathrm{mm}`$ along the measuring direction and of about 0.5$`\mathrm{mm}`$ for the two other directions. The thermal gradient was determined by measuring the temperature difference $`\mathrm{\Delta }T`$ between the junctions of a differential Au/Fe-Chromel thermocouple. The junctions of this thermocouple have been glued onto the sample using GE varnish.foot1  $`\mathrm{\Delta }T`$ varied between 0.5% and 2% of the absolute temperature, which has been stabilized for each data point. Errors due to radiation loss, which could occur at higher temperatures, are avoided in our experimental setup. Stoichiometric oxygen contents in $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$ and $`\mathrm{La}_2\mathrm{NiO}_4`$ were achieved by annealing in high vacuum and $`\mathrm{CO}/\mathrm{CO}_2`$-atmosphere, respectively.

## III Previous results

Prior to discussing our new experimental results, we review previous resultsNakamura et al. (1991); Baberski et al. (1998) on the striking doping dependence of the thermal conductivity of doped $`\mathrm{La}_2\mathrm{CuO}_4`$. At first, we concentrate on the low temperature phononic peak in the thermal conductivity of $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$. As is evident from the lower panel of Fig. 1, which reproduces $`\kappa _c`$ of $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ ($`x=0`$, 0.1, 0.15, 0.2, 0.3) as published by Nakamura et al.,Nakamura et al. (1991) this peak evolves non-monotonically with increasing Sr-content: A well pronounced phononic low-$`T`$ peak in $`\kappa _c`$ is only present at $`x=0`$ and $`x=0.3`$, whereas at intermediate doping levels ($`0.1x0.2`$) a peak structure is hardly identifyable. Note that the material is a superconductor in this doping range whereas it is insulating and metallic at $`x=0`$ and $`x=0.3`$, respectively.

Baberski et al. have pointed out that the suppressed low-$`T`$ peak at intermediate Sr-doping reappears upon RE-doping, provided that this doping induces the so-called LTT-phase (Low Temperature Tetragonal) and thereby suppresses superconductivity.Baberski et al. (1998) This is illustrated in the upper panel of Fig. 1, where $`\kappa `$ of polycrystalline Pr- and Nd- doped $`\mathrm{La}_{2\mathrm{x}\mathrm{y}}(\mathrm{RE})_\mathrm{y}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ at the finite Sr-content $`x=0.12`$ is shown. The Pr-doped compound does not undergo the transition to the LTT-phase. Like $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$, it remains in the so-called LTO-phase (Low Temperature Orthorhombic) and is a superconductor.Schäfer et al. (1994) Its thermal conductivity monotonically decreases with decreasing $`T`$ and hence is very similar to the aforementioned findings by Nakamura et al. for $`\kappa `$ of $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ at $`0.1x0.2`$. The situation is completely different in the non superconducting, Nd-doped compound, which is in the LTT-phase below $`T_{LT}80`$ K: Here, $`\kappa _{\mathrm{ph}}`$ abruptly enhances at $`T_{LT}`$ and exhibits a well pronounced phononic peak around 25 K.

The doping dependence of $`\kappa `$ described above must be attributed to strong changes in the phonon thermal conductivity $`\kappa _{\mathrm{ph}}`$ of this material: Non-phononic contributions to $`\kappa `$ can be excluded in the out-of-plane direction, since the electrical conductivityNakamura and Uchida (1992, 1993) and the magnetic couplingThio et al. (1988) along the $`c`$-axis are too small to allow significant electronic and magnetic thermal conduction. Concomitantly, the doping dependence of the low-$`T`$ peak is evident along all crystal directionsNakamura et al. (1991) (c.f. also the data below).

In crystalline materials, the low-$`T`$ peak in $`\kappa _{\mathrm{ph}}`$ is very sensitive to scattering of phonons. Generally, the height of this peak reduces with increasing phonon scattering rates. One important mechanism in this regard is scattering of phonons by impurities. In alloyed compounds like doped $`\mathrm{La}_2\mathrm{CuO}_4`$ this phonon-impurity scattering is induced by non-uniform ions on one lattice site. Upon alloying, this should lead to a gradual reduction of the phonon peak.Berman and Brock (1965) Even though phonon-defect scattering inevitably must be present in this material, it is not a scattering mechanism which solely dominates $`\kappa _{\mathrm{ph}}`$, because the phonon peak evolves non-monotonically upon Sr-doping and even reappears at additional RE-doping.

Another scattering mechanism for phonons in this material is suggested by the abrupt change of $`\kappa _{\mathrm{ph}}`$ at the structural phase transition, which occurs in the Nd-doped compound. It appears that phonons are scattered stronger in the LTO-phase. This interpretation is confirmed by neutron scattering studies on Nd-doped $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ at $`T_{LT}`$: The acoustic phonon line width abruptly decreases at the transition from the LTO- into the LTT-phase, signaling a proportional decrease of scattering processes.Pintschovius and Reichardt (1998) Anomalous phonon thermal conductivity in the vicinity of structural phase transitions is well known from a number of perovskite oxides as, for example, $`\mathrm{SrTiO}_3`$ and $`\mathrm{KTaO}_3`$. There, a suppression of $`\kappa _{\mathrm{ph}}`$ is caused by enhanced scattering of acoustic phonon modes due to their energetic degeneracy with soft optical phonon modes.Steigmeier (1968); Barret and Holland (1970) Indeed, soft phonon branches do exist in the LTO-phase of doped $`\mathrm{La}_2\mathrm{CuO}_4`$,Birgeneau et al. (1987); Pintschovius (1990); Martínez et al. (1991); Keimer et al. (1993) which could cause a suppression of $`\kappa _{\mathrm{ph}}`$ (in the following this scattering mechanism will briefly be named ’soft-phonon scattering’). The change of $`\kappa _{\mathrm{ph}}`$ at $`T_{LT}`$ then would follow from the discontinuous hardening of the soft phonon branch in the LTT-phaseMartínez et al. (1991); Keimer et al. (1993) and an associated reduced scattering rate of acoustic phonons.foot2

Since all compounds with a suppressed phonon peak are in the LTO-phase at low $`T`$, soft-phonon scattering could be important for understanding the suppression of the peak as well.foot2a  There is, however, not a one to one correlation between the LTO-phase, i.e., the possible presence of soft phonon scattering, and the suppression of the peak. This is most obvious in $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$, because the soft phonon properties only slightly change for $`x0.2`$, whereas the phonon peak is maximum for $`\mathrm{La}_2\mathrm{CuO}_4`$.

Baberski et al. have noticed that both impurity scattering and soft-phonon scattering separately do not allow to understand the suppression of the phonon peak. They therefore suggested, that the suppression of the phononic peak could be correlated with the superconducting properties, since the peak suppression occurs whenever the material is superconducting.Baberski et al. (1998) In their model, they proposed that stripes couple to phonons and thereby cause an unconventional scattering mechanism for phonons depending on the dynamics of stripes: A suppression of the phonon peak then is the consequence of fluctuating stripes, which are present in the superconducting compounds, while $`\kappa _{\mathrm{ph}}`$ exhibits a usual phonon peak when stripes are static or even absent in the non-superconducting cases. There, static stripes are only present in the LTT-phase. Therefore, instead of being caused by an abrupt softening of optical phonon modes, the jump of $`\kappa _{\mathrm{ph}}`$ at $`T_{LT}`$ could just as well originate from a change of the stripe dynamics from static (LTT) to fluctuating (LTO). Even though a consistent interpretation of the data discussed afore is possible within such a model, its validity is questionable because the actual role of impurity- and soft-phonon scattering remains unclear.

## IV New results

In order to elucidate the importance of this conventional scattering mechanisms we now turn to our measurements on single crystals, which provide profound information in this regard.

In Fig. 2 we present the complete data set of $`\kappa `$ of $`\mathrm{La}_{1.8\mathrm{x}}\mathrm{Eu}_{0.2}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ ($`x=0`$, 0.08, 0.15, 0.2) measured along the $`c`$-axis ($`\kappa _c`$) and parallel to the $`ab`$-planes ($`\kappa _{ab}`$) in the temperature-range $`7300`$ K. All compounds are in the LTT-phase at $`T<T_{LT}`$, with a variation of $`T_{LT}`$ between about 110 K and 130 K. For $`T>T_{LT}`$ and $`x0.15`$ the structure is LTO in the investigated temperature range. At $`x=0.2`$ the compound undergoes a further structural phase transition from the LTO- into the so called HTT-phase (High-Temperature-Tetragonal) at $`T_{HT}220`$ K.Klauss et al. (2000)

$`\kappa _c`$ of $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$ \[cf. panel (a) of Fig. 2\] exhibits a low-temperature peak around 20 K with a falling edge that continuously extends to $`T_{LT}130`$ K. A jumplike decrease occurs at $`T_{LT}`$, followed by a constant $`\kappa _c`$ up to room temperature. The inset of Fig. 2 depicts a cooling- and heating curve of $`\kappa _c`$ in the area of the jump. A clear hysteresis being a characteristic feature of first-order phase transitions is evident.

For $`T50`$ K, the thermal conductivity along the $`\mathrm{CuO}_2`$-planes, $`\kappa _{ab}`$, is comparable to $`\kappa _c`$. The peak centered around 20 K is slightly smaller; we find $`\kappa _c^{max}/\kappa _{ab}^{max}1.3`$, which is similar to the findings in non-doped $`\mathrm{La}_2\mathrm{CuO}_4`$ and isostructural $`\mathrm{La}_{5/3}\mathrm{Sr}_{1/3}\mathrm{NiO}_4`$.Nakamura et al. (1991); Hess et al. (1999) For $`T50`$ K the temperature dependencies of $`\kappa _{ab}`$ and $`\kappa _c`$ differ completely: With rising temperature, $`\kappa _{ab}`$ strongly increases and evolves into a broad peak at room temperature. A step-like anomaly is present at $`T_{LT}`$, which is in a similar way hysteretic as the anomaly in $`\kappa _c`$ (not shown).

The panels (b) and (c) of Fig. 2 show $`\kappa _{ab}`$ and $`\kappa _c`$ of $`\mathrm{La}_{1.8\mathrm{x}}\mathrm{Eu}_{0.2}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ for $`x=0.08`$ and $`x=0.15`$, respectively. For both compounds $`\kappa _{ab}`$ and $`\kappa _c`$ exhibit a low temperature peak, where the peak of $`\kappa _{ab}`$ is again slightly smaller than that of $`\kappa _c`$. The temperature dependence of $`\kappa _c`$ is very similar to that of $`\kappa _c`$ of $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$ in the whole temperature range. As in the case at $`x=0`$, $`\kappa _{ab}`$ deviates from the qualitative $`T`$-dependence of $`\kappa _c`$ above $`50`$ K and evolves into an increase with rising temperature. A steplike anomaly at $`T_{LT}`$ exists in each case, but no high temperature peak as in $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$ is present.

A comparable anisotropy as in the previous compounds is also evident at $`x=0.2`$ \[cf. panel (d) of Fig. 2\], but the low temperature peak in $`\kappa _c`$ is almost completely suppressed and $`\kappa _c`$ slightly increases for $`T75`$ K with nearly constant positive slope. The low temperature peak is even absent in $`\kappa _{ab}`$. Here, $`\kappa _{ab}(T)`$ monotonically increases upon heating in the entire temperature range. At $`T_{LT}110`$ K no anomaly is present, neither in $`\kappa _{ab}`$ nor in $`\kappa _c`$. Apparently, the transition at $`T_{HT}`$ also causes no anomaly in the thermal conductivity.

## V Discussion

### V.1 Anisotropy

As mentioned before, $`\kappa _c`$ is purely phononic, since electronic and magnetic contributions can be ruled out for heat transport along the $`c`$-axis. In order to understand the anisotropies of $`\kappa `$ for $`T50`$ K, it is reasonable to distinguish the cases $`x=0`$ and $`x>0`$. It has been shown in Ref. 18 that in insulating, antiferromagnetic $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$ the high temperature peak of $`\kappa _{ab}`$ must be explained by magnetic contributions. Upon Sr-doping, the doped holes destroy the magnetic order which leads to a strong suppression of the high temperature peak (cf. Fig. 2). Simultaneously, the material becomes electrically conducting within the $`ab`$-planes.Hess et al. (2003); Nakamura et al. (1991) Therefore, the much smaller, but apparently still existing anisotropy in the Sr-doped compounds is due to electronic, rather than magnetic, contributions to $`\kappa _{ab}`$. Indeed, an estimation using the Wiedemann-Franz-law yields electronic contributions of the same order of magnitude as the observed anisotropy for $`\kappa _{ab}`$ and negligible electronic contributions for $`\kappa _c`$. The magnetic and electronic contributions to $`\kappa _{ab}`$ become important for $`T50`$ K. Below this temperature, $`\kappa _{ab}`$ can be considered to be primarily phononic with a similar magnitude of $`\kappa _{\mathrm{ph}}`$ as $`\kappa _c`$.foot3  This is nicely confirmed by the very similar temperature dependence of $`\kappa _{ab}`$ and $`\kappa _c`$ in this temperature range, as shown in Fig. 3. It is remarkable that for all compounds the phonon peak of $`\kappa _c`$ is slightly larger than that of $`\kappa _{ab}`$ as long as a peak is clearly resolved in both quantities. A possible explanation could be related to slightly different velocities of the acoustic phonons along the $`ab`$\- and $`c`$-directions.Pintschovius et al. (1991)

### V.2 Phononic peak and jump at $`T_{LT}`$

The development of the low temperature peak of $`\kappa _{\mathrm{ph}}`$ upon doping can be viewed in Fig. 3. Obviously, the peak size is gradually reduced as the Sr-content increases. A similar result is found for the size of the jump $`\mathrm{\Delta }\kappa _c`$ at $`T_{LT}`$ which is shown as a function of doping in the inset of Fig. 3. While the gradual reduction of the peak size can straightforwardly be explained by phonon-impurity scattering, the reduction of the jump at $`T_{LT}`$ requires further comments. First of all, we stress the presence of a jump in insulating $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$. It clearly proofs that it is caused by soft-phonon scattering in the LTO-phase, which therefore has to be regarded as an important scattering channel for phonons indeed. When this scattering channel becomes active in the LTO-phase, the relative importance of phonon-impurity scattering is reduced. The gradual reduction of the jump size $`\mathrm{\Delta }\kappa _c`$ with increasing Sr-content may therefore be attributed to doping induced phonon-impurity scattering as well. However, we should note that the structural differences between the LTT and LTO-phases diminish with increasing Sr-content. Therefore, apart from phonon-impurity scattering, structural reasons play a role in the suppression of the jump at $`T_{LT}`$. This is consistent with the observation by Baberski et al., that the jump disappears, whenever the tilting angle of the $`\mathrm{CuO}_6`$-octahedra becomes lower than a critical value.Baberski et al. (1998)

It is now very instructive to consider $`\kappa _{\mathrm{ph}}`$ of pure $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ (cf. Fig. 1) for comparison. Apparently, doping Sr into $`\mathrm{La}_2\mathrm{CuO}_4`$ suppresses the phonon peak much more effectively than in the Eu-doped counterpart: On the one hand the phonon peak of non-doped $`\mathrm{La}_2\mathrm{CuO}_4`$ is by a factor of about two larger than in non-doped $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$. On the other hand the phonon peaks in $`\mathrm{La}_{1.8\mathrm{x}}\mathrm{Eu}_{0.2}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ are clearly better developed for finite Sr-contents with higher maximum values as in $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$. One has therefore to conclude that a further scattering mechanism exists in $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ which becomes more important upon doping and which is absent or at least much weaker in the Eu-doped compounds. In this case, this mechanism must cause a more effective scattering than the surely present phonon-impurity scattering induced by the Eu-ions. One plausible candidate (besides scattering due to dynamic stripes) for such a mechanism is soft-phonon scattering connected with the structural instability of the LTO-phase since the soft-phonon energies decrease further upon Sr-dopingfoot4  and hence enhance soft-phonon scattering.

In order to judge the relevance of soft-phonon scattering a measure for its strength with respect to phonon-impurity scattering is necessary. Such a measure could be achieved, for example, by comparison with a compound where neither a lattice instability nor doped impurities are present. Though being no cuprate, isostructural $`\mathrm{La}_2\mathrm{NiO}_4`$ is yet a well suited candidate perfectly fulfilling these requirements. The phonon spectra of $`\mathrm{La}_2\mathrm{CuO}_4`$ and $`\mathrm{La}_2\mathrm{NiO}_4`$ are almost identicalPintschovius et al. (1989, 2001) because the atomic masses of Ni and Cu are very similar. At low temperatures ($`T73`$ K) this compound is in the LTT-phase and hence a structural instability does not exist.Rodríguez-Carvajal et al. (1988) Moreover, no phonon-impurity scattering is present in this non-doped compound.

Fig. 4 presents our result for phononicfoot5  $`\kappa _{ab}`$ of electrically insulating $`\mathrm{La}_2\mathrm{NiO}_4`$ in comparison with $`\kappa _{ab}`$ of $`\mathrm{La}_2\mathrm{CuO}_4`$ as measured by Nakamura et al. and with $`\kappa _{ab}`$ of $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$. It is intriguing that the phonon peak of $`\mathrm{La}_2\mathrm{NiO}_4`$ is about one order of magnitude larger than the peak in both $`\mathrm{La}_2\mathrm{CuO}_4`$ and $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$. At $`T_{LT}73`$ K, the structure changes from LTT to LTO, which also in this material causes a jumplike decrease in the thermal conductivity. Compared to $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$ the jump size is larger by a factor of about 10. In the LTO-phase, $`\kappa _{ab}`$ is of the similar size as in $`\mathrm{La}_2\mathrm{CuO}_4`$.

It immediately follows from this observation that in $`\mathrm{La}_2\mathrm{CuO}_4`$ as well as in $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$ the heat carrying phonons are subject of severe scattering at low temperatures. In both cases this scattering obviously suppresses the peak of $`\kappa _{\mathrm{ph}}`$ by about one order of magnitude. Translated to the strength of the two discussed scattering mechanisms this means that referring to a non-impurity-doped LTT-phase both soft-phonon scattering and phonon-impurity scattering independently of each other reduce the peak of $`\kappa _{\mathrm{ph}}`$ by almost the same amount.

This conclusion now allows us to understand the doping dependence of the phonon peak, when a combination of both, phonon-impurity- and soft-phonon scattering is consired: In $`\mathrm{La}_{1.8\mathrm{x}}\mathrm{Eu}_{0.2}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$, where the thermal conductivity is already strongly influenced by the Eu-ions, Sr-doping simply further increases the phonon-impurity scattering rate. This consequently leads to a gradual reduction of the phonon peak with increasing Sr-content. Soft-phonon scattering is not relevant here.

In $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ with the same Sr-content than a Eu-doped counterpart, the phonon-impurity scattering rate is slightly reduced since no Eu-ions are present. Yet, the resulting phonon thermal conductivity is somewhat smaller. Hence, the effect of a slightly reduced phonon-impurity scattering rate must be overcompensated by soft-phonon scattering. This is indeed reasonable since it is a qualitatively different scattering mechanism whose importance upon Sr-doping grows since thereby the relevant soft modes soften further.Braden et al. (1994)

The reappearance of the phonon peak in $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ at high Sr-concentrations as $`x=0.3`$ can now be understood as a natural consequence of the doping dependence of $`T_{HT}`$. For $`x0.22`$, the LTO-phase disappears in $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ and the structure remains in the HTT-phase for all temperatures.Takagi et al. (1992) Hence, for $`x=0.3`$ soft-phonon scattering associated with the LTO-phase is not active for such overdoped compounds. The result of such decreased phonon scattering is the reappearance of the phonon peak.

It is necessary to mention that the growing density of holes as the Sr-content is increased in principle could be a further source of scattering and therefore contribute to the reduction of $`\kappa _{\mathrm{ph}}`$ and in particular the phonon peak via phonon-electron scattering. However, since the transition to superconductivity in $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ causes no significant anomaly in $`\kappa _c`$ at $`T_c`$ (cf. Fig 1), this scattering mechanism is usually considered to be unimportant in this material and its (RE)-doped relatives.Baberski et al. (1998) Yet if this scattering channel is active in this material no inconsistency to our interpretation arises because in this case the phonon-’impurity’ scattering induced by Sr-ions can be regarded simply as slightly more effective than the phonon-impurity scattering induced by Eu-ions.

### V.3 Stripes as a scattering mechanism for phonons?

The above discussion of $`\kappa _{\mathrm{ph}}`$ of the insulating materials $`\mathrm{La}_2\mathrm{CuO}_4`$, $`\mathrm{La}_{1.8}\mathrm{Eu}_{0.2}\mathrm{CuO}_4`$ and $`\mathrm{La}_2\mathrm{NiO}_4`$ provides unambiguous evidence that phonon-impurity and soft-phonon scattering are important scattering mechanisms for phonons in doped $`\mathrm{La}_2\mathrm{CuO}_4`$. Since the major doping dependencies of the low-$`T`$ peak can be explained based on these mechanisms without any problems, there is no compelling reason to incorporate a stripe induced scattering channel in the data interpretation. It is unlikely though that stripes in doped $`\mathrm{La}_2\mathrm{CuO}_4`$ have no effect at all on $`\kappa _{\mathrm{ph}}`$ because the aforementioned stripe-induced phonon scattering in the nickelatesHess et al. (1999); Cassel et al. (1999) is an unambiguous physical fact. However, it appears extremely difficult to prove the existence and to study the strength of purely stripe induced scattering in doped $`\mathrm{La}_2\mathrm{CuO}_4`$ via $`\kappa _{\mathrm{ph}}`$. This is not only because phonon-impurity and soft-phonon scattering are already dominating $`\kappa _{\mathrm{ph}}`$ in insulating $`\mathrm{La}_2\mathrm{CuO}_4`$; Due to the intimate relation of LTT-phase and stripe dynamics, scattering on fluctuating stripes in the Sr-doped compounds could be viewed as an altered but already existing soft-phonon scattering, i.e., from this point of view these two scattering mechanisms are conceptually indistinguishable.

### V.4 High temperature increase of $`\kappa _{\mathrm{ph}}`$

For completeness, we briefly mention the unusual temperature dependence of $`\kappa _c`$ which is evident at higher temperatures $`T100`$ K. Instead of a decrease as $`T^1`$, which at elevated temperatures is usually expected for thermal conductivity by acoustic phonons, $`\kappa _c`$ is almost temperature independent or even increases with rising $`T`$. As has been shown in Ref. 60, this strong deviation from the usual behavior arises due to thermal conduction by dispersive optical phonons in addition to the usual contribution by acoustic phonons.

## VI Summary

In summary, we have presented experimental results on the thermal conductivity of $`\mathrm{La}_{1.8\mathrm{x}}\mathrm{Eu}_{0.2}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ for a wide doping range of Sr. The analysis of our data suggests that in this material phonons are strongly scattered on the doped impurities, i.e., Sr and Eu, as well as on soft phonons that are present in the LTO-phase of this material. Comparison of our data with the thermal conductivity of $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ and isostructural $`\mathrm{La}_2\mathrm{NiO}_4`$ leads to the conclusion that these scattering mechanisms are most relevant in $`\mathrm{La}_{2\mathrm{x}}\mathrm{Sr}_\mathrm{x}\mathrm{CuO}_4`$ as well. In particular, in contrast to other studiesBaberski et al. (1998); Sun et al. (2003a, b), there remains no direct evidence that the stripe correlations cause a relevant scattering channel in doped $`\mathrm{La}_2\mathrm{CuO}_4`$.

## VII Acknowledgments

C.H. acknowledges support by the DFG through HE3439/3-1. We thank N. Jenkins for proofreading the manuscript.


text:
# Population fragmentation and party dynamics in an evolutionary political game

\[

## Abstract

We examine kinetic symmetry breaking phenomena in an evolutionary political game in which voters, inhabiting a multidimensional ideological space, cast ballots via selection mechanisms subject to the competing social forces of conformity and dissent. Our understanding of the spatiotemporally complex population dynamics is informed by a system of nonlinear replicator equations, discrete deterministic cousin of the original stochastic Seceder Model.

\]

Fragmentation, dispersal & coalescence are crucial processes in a multitude of highly correlated, dynamically complex many-body systems, particularly those far from equilibrium where stochasticity, symmetry, and self-organization can conspire to generate rich pattern formation phenomena . Natural examples abound- the rings of Saturn, its Cassini and lesser divisions, consist of a gravitationally bound myriad of ever-colliding rocky chucks, macroscopically stable within an intricately detailed set of circular bands . Alternatively, consider the coat of the leopard, zebra, or giraffe- each a glorious product of competing reaction-diffusion biochemistry . In these instances, the dynamical forces are physiochemical in nature, but movement can be inspired by equally compelling, though fundamentally different mechanisms, especially in biological, economic and sociological contexts. Visual/verbal cues (i.e., information propagation) within a school of fish  or flock of birds  lead to dynamic domain formation & destruction. Volatile populations of informed, interacting stock market traders can condense, exhibiting herd-like behavior . Rich stochastic dynamics, as well as phase transition phenomena, are evident in various evolutionary minority (e.g., El Farol Bar) , public-goods , and other societal selection games, such as the Seceder Model , which introduces a novel dynamical frustration via competing tendencies of being distinct, yet part of the group. Here, we consider the seceder mechanism within an evolving political populace, uncovering an interwoven set of strange attractors whose reign is governed by kinetic symmetry-breaking within an ideological plane. We find that the conflicting tendencies toward conformity and dissent yield rich spatiotemporal party population dynamics. The seceder Catch-22 is that dissenting parties provide alternatives, but frequently grow in popularity themselves, inevitably suffering rebellion from within. Clearly, these simple models are not intended to capture all important features; the goal, however, is to gain appreciation of robust, universal aspects, as was done, e.g., for Ising systems in the case of ferromagnetic, liquid-gas, and binary alloy critical phenomena , or quadratic maps & Lotka-Volterra coupled ODEs, in the matter of chaotic dynamics .

For ease of presentation, we discuss the illustrative case of 3-party dynamics driven by a triplet-selection mechanism; i.e., voters poll three individuals to inform their political decision. We are lead, initially, to study the fixed points (FPs) associated with the following set of symmetry-broken replicator equations :

$`\dot{R}=R^3+3R(G^2+B^2)+\alpha RGBR`$

$`\dot{G}=G^3+3G(R^2+B^2)+\beta RGBG`$

$`\dot{B}=B^3+3B(R^2+G^2)+\gamma RGBB`$

where $`R,G,\&B1`$ represent the relative populations of the three parties, Red, Green & Black, while probability conservation demands $`\alpha `$+$`\beta `$+$`\gamma `$=6. The rate variables are dictated by the competing tendencies of, resp., homogeneity, distinction, and the democratic dilemma- equal representation of the three parties. More specifically, the purely cubic term results when the triplet selection group has uniform representation of one party; the decision maker switches to (or remains in, as the case may be…) that party. In the event that the polled group has an unequal distribution; i.e., a minority party, outnumbered 2 to 1, the underdog triumphs- thus, the seceder linear-quadratic terms. Finally, the trilinear contributions arise when the opinion group has one individual from each of the three parties. We consider successive levels of symmetry breaking- first, if $`\alpha `$=$`\beta `$, we have effectively, an isosceles geometry, $`R`$ & $`G`$ are equivalent in their relation to $`B`$ and we find the following population FPs: $`(R,G,B)=(\frac{1}{2},\frac{1}{2},0)`$ and $`(\frac{2}{8\alpha },\frac{2}{8\alpha },\frac{4\alpha }{8\alpha }).`$ For $`\alpha <4`$, the 3-party solution is stable, 2-party dynamics unstable; for $`\alpha 4`$ vice versa. Single party domination, e.g., $`(0,0,1),`$ is never stable. Clearly, $`\alpha =2`$ would represent a fully symmetrized situation with the three groups equidistant from each other, an equilateral (EQ) arrangement with evenly distributed populations $`(\frac{1}{3},\frac{1}{3},\frac{1}{3}).`$ By contrast, $`\alpha =3`$ corresponds to the standard Seceder Model in one dimension , with the odd party, $`B,`$ poised at the midpoint of the line segment connecting ideologically opposed groups, $`R`$ & $`G,`$ at the ends; in this case, the stable population FP is the familiar $`(\frac{2}{5},\frac{2}{5},\frac{1}{5}).`$ In fact, any short isosceles (SI) triangle arrangement of the parties, with sides $`RG`$$`>`$$`RB`$$`=`$$`GB,`$ is drawn to a chaotic orbit about that fixed point; see later. Alternatively, the geometric arrangement could be tall isosceles (TI), with $`RG`$$`<`$$`RB`$$`=`$$`GB,`$ so that there is effectively a remote party equidistant from the others; here, connection to the stochastic Seceder Model suggests $`\alpha `$=0; the dynamical

flows are drawn to the attractor $`(\frac{1}{4},\frac{1}{2},\frac{1}{4}).`$ We reveal, in Figure 1, the full phase-diagram of this nonlinear replicator system in the $`\alpha \beta `$-plane, relaxing the isoscelean geometric constraint. Surprisingly, a stable 3-party dynamic persists for quite a range of parameter values, occupying the triangular region bounded by the lines $`\alpha =2,\beta =2,\alpha +\beta =8`$ (i.e., $`\gamma =2`$). The unique, full symmetry EQ FP, $`(\frac{1}{3},\frac{1}{3},\frac{1}{3}),`$ shown as a red diamond at the point $`\alpha =\beta =\gamma =2,`$ is surrounded by a trio of permutation-related TI (SI) FPs, indicated by blue circles (green squares). Oddly enough, 3-party politics exists even for modestly negative values of $`\alpha ,\beta ,`$ or $`\gamma `$; i.e., explicit biasing against a particular party. For sufficiently great forcing, however, the dynamics switches from 3 to 2-party politics; e.g., with $`\gamma <2,`$ all flows are driven to the same $`(R,G,B)=(\frac{1}{2},\frac{1}{2},0)`$ FP, with a 50-50 split in the population between the two surviving parties. Along special lines of symmetry, such as $`\alpha =\beta <2,`$ there is, numerically, a first-order coexistence between competing relevant 2-party FPs; in the mentioned instance, $`(0,\frac{1}{2},\frac{1}{2})`$ and $`(\frac{1}{2},0,\frac{1}{2})`$. All this is in strong contrast to the 3-party region of the phase diagram, where there is a unique FP associated with each choice of $`\alpha \beta \gamma `$ values.

Having investigated the role of symmetry-breaking within the context of the well-mixed, deterministic replicator equations, we now explore its manifestation for a related spatially-extended, stochastic model in which, one might imagine, voters occupy specific positions in a 2D ideological plane, with axes corresponding to a pair of compelling political issues, e.g., taxation and international cooperation, running the gamut from extreme leftist to ultra-right wing conservative positions. We’ll see shortly that for this stochastic Seceder Model, kinetic symmetry breaking (KSB) arises in an entirely intrinsic manner via fluctuations built-up in this far-from-equilibrium system. We make long-time runs, $`t`$=5x105 generations, with substantial population sizes, $`N`$=10$`^3,`$ running simulations in 2d with the entire population initialized in a single party at the origin. As discovered previously , the population rapidly fragments and disperses but within a thousand generations coalesces into three divergent groups, heading radially away from the origin at a slightly sublinear speed  with equally-sized

party populations. In fact, if one carefully monitors party membership over time, there are impressive fluctuations about (as well as away from!) the neighborhood of this dynamical $`(\frac{1}{3},\frac{1}{3},\frac{1}{3})`$ EQ FP. In Figure 2a, we show the entire time series for the three parties- here indicated by the colors red, green and black. We see immediately that for much of $`t<1`$x10$`^5,`$ the average party population is indeed $`N/3`$ but in due course, one of the parties, in this case the Greens @$`t`$$``$78000, suffers a severe drop in membership, suddenly falling to 200; simultaneously, the Red and Black Parties rise, fluctuating about a mean of 400. In Figure 2b, we show an expanded detail of the subsequent, especially active epoch. Apparently, the political system has dynamically migrated to a chaotic stay about the $`(\frac{2}{5},\frac{1}{5},\frac{2}{5})`$ SI FP, where two parties, Red & Black, equally dominate the third, Green. This state of affairs persists until $`t`$$``$91000, when a second abrupt dynamical transition occurs- this time to the $`(\frac{1}{2},\frac{1}{4},\frac{1}{4})`$ TI FP, in which there is a single popular party, here the Reds, with twice the membership of two equal minority parties, Black and Green. Finally, just beyond $`t`$=100000, the system reverts to the original, unique EQ FP. Overlays to Figure 2b reveal relevant geometric aspects, including evolution of various pairwise separations between parties (blue), as well as radial distances from the origin (party colors). We notice, for example, when the Green Party membership plummets, the party actually suffers a stochastic stall (green trace plateaus…), its radial motion ceasing as population redistribution takes place. The SI$``$TI transition is always marked by the simultaneous stall of all three parties (red, green, and black traces all horizontal), while the triangular geometry switches to one short and two (nearly) equal long sides. The return the EQ FP for $`t`$$`>`$100000 is, of course, temporary- see, again, Figure 2a. Note, in particular, the rather long run, $`t`$=2.5-3.9x105, during which the political dynamic is controlled by the $`(\frac{2}{5},\frac{2}{5},\frac{1}{5})`$ SI FP, though even that impressive interlude is punctuated by intermittent, though typically short-lived ($``$ 1-2000 generations) TI episodes, three of which are clearly visible in the figure.

It is informative to study the time evolution of the party populations from the vantage point of a Poincare type section; we show in Figure 3 the fractional populations of the Black and Red Parties for the entirety of the original simulation. For the time interval $`t`$=2.0-4.0x10$`^5,`$ an epoch strongly dominated, recall Figure 2a, by the $`(R,G,B)=(\frac{2}{5},\frac{2}{5},\frac{1}{5})`$ FP, indicated by leftmost green square, a symbol shared by its SI FP permuted cousins. As previously, the TI FPs are shown by blue circles, and the solitary, but evenly distributed, EQ FP is the red diamond. The time interval was chosen specifically to show governing role of a particular SI FP, which gives rise to a compelling crescent-shaped feature in the Poincare trace. Note, especially, the position of complementary blue TI FPs in the wings of the crescent, which faithfully track the relatively less frequent excursions of the system to states in which a single majority party strongly dominates the political dynamics; here, Red or Green Parties with memberships in the neighborhood of 500- in such instances, the Black Party moves up momentarily to 250. It is apparent that the interval $`t`$=2.0-2.4x105 is responsible for the faint black cloud centered about the EQ FP, complementary to the solid crescent. For the given interval, the system spends little time with either the Green or Red Parties in the weak minority, so those sections of the Poincare plot are rather sparse, but in the infinite time limit, or even if one considers the full 1/2 million generations implicit in Figure 2a, the Poincare section (black and gray dots taken together) has an easily perceived elliptical form resulting from the superposition of three symmetrically permuted cusps pinned at each of the green square, SI FPs. In this way, wings overlap, where the crescents share a common, blue circle TI FP. Of course, this ellipse (or ellipsoid, in the full 3d plot with (R,G,B) on the three axes…), does not have a sharp boundary, though the fall-off is rather rapid. The soft asymptotic ellipse, as well as the short-

timescale crescent constitute the signature details of the Seceder Model Poincare section.

We can decipher these matters a bit by studying the party membership PDFs, obtained by isolating intervals of time in which a single FP dominates the dynamics. In Figure 4, the population PDFs associated with the EQ, SI, TI attractors are shown, resp., in red, green, and blue. We construct the red $`(\frac{1}{3},\frac{1}{3},\frac{1}{3})`$ EQ PDF by sampling over some 200000 noncontiguous generations, primarily $`t`$$`<`$80000 and $`t`$$`>`$420000 and find three overlapped, essentially identical, asymmetric, nongaussian, skewed distributions, peaked near $`N/3`$=333, low slung with very substantial widths. By contrast, focussing on the 150000 generation interval commencing at $`t`$=240000, where $`(\frac{2}{5},\frac{2}{5},\frac{1}{5})`$ FP reigns, we uncover the membership PDF associated with the green SI FP. Clearly a different beast altogether, it is bimodal, with 2 superposed, short, but broad components (mean@400), long tails on the low side- still substantial even at 250, and a single sharp, nearly symmetric peak centered at 200, considerably narrower, with rather small probability to rise above 250 or fall below 150. These features, esp. the latter, are evident in the crescent boundaries of Figure 3. Finally, the TI PDF, in blue, also bimodal, is characterized by a very narrow peak for the solitary majority party @500,

and two identical overlapped broader contributions with mean at 250 for the equal pair of minority parties. The nature of the various tails, as well as the mixed population size, go far in explaining why certain dynamical transitions, such as EQ$``$SI or SI$``$TI, are considerably more common than others; e.g., the system rarely tunnels EQ$``$TI. We seek a deeper understanding of the transition probabilities associated with these important global shifts in the population dynamics and have taken a step in this direction, studying the decay modes and lifetime distribution of an initial EQ state- see Figure 5. We discover that the branching ratio is highly biased in favor of the EQ$``$SI decay route; indeed, 97.6%(2.4%) of the time the transition is to the SI(TI) attractor.

Interestingly, the essential aspects of KSB can be analyzed via the replicator equations; in particular, we consider the case of $`d+1`$ symmetrically arrayed parties in $`d`$ dimensional space (i.e, 2D equilateral triangle, 3D tetrahedron, etc.) and examine what happens if a single party, with fractional membership $`z`$, let’s say, kinetically breaks the symmetry, moving closer/farther from the remaining parties, each of whom possess equal membership $`x=(1z)/d.`$ For this characteristic situation, the rate equation for the anomalous party reads: $`\dot{z}=z^3+3dz(1+\delta (d1))x^2z,`$ where the symmetry-breaking parameter $`\delta =0,\frac{1}{2},1,`$ in SI, EQ, TI geometries, with associated fixed point values $`z^{}=\frac{3d}{3+d},\frac{1}{d+1},\frac{1}{2}`$. The $`\delta =\frac{1}{2}`$ result is expected, the total population distributed evenly amongst the $`d+1`$ groups. If, however, the anomalous party distances itself from the rest moving farther away, this subpopulation grows to a dimension-independent percentage of 50%, the remaining parties sharing equally the other half of the population; of course, for the 2D case we’ve focussed on thus far, this corresponds to the TI FP $`(\frac{1}{4},\frac{1}{4},\frac{1}{2}).`$ Although fluctuations can permit higher fractions for this dominant party- our earlier simulations indicated, see Figures 2&4, upwards of 3/5, there are apparently dynamical constraints that strongly limit the maximum membership of this favored party. Should the party in power become too fashionable,

the seceder mechanism inceases the likelihood of desertion from within. Finally, in the situation where the stochasticity draws the errant party closer to the others, a stable SI FP exists only stable below the upper critical dimension $`d_c=3!`$ Consequently, fluctuations in $`d3`$ bring about the successive demise of all but three of the remaining parties, the system suffering cascadal collapse & dimensional reduction to a 2D ideological hyperplane .

Financial support for Tim HH has been provided by NSF DMR-0083204, Condensed Matter Theory.


text:
# Phase field approach for modeling intracellular dynamics

## Abstract

We introduce a phase field approach for diffusion inside and outside a closed cell with damping and with source terms at the interface. The method is compared to exact solutions (where possible) and the more traditional finite element method. It is shown to be very accurate, easy to implement and computationally inexpensive. We apply our method to a recently introduced model for chemotaxis by Rappel et al. \[Biophys. J. 83, 1361 (2002)\].

In dealing with free boundary problems, the so called phase field approach pf  appears as a method of choice. It has successfully been applied to various problems ranging from dendritic solidification kr98 , viscous fingering pf-hs  and crack propagation pf-crack . In the spirit of time-dependent Ginzburg-Landau models, the method avoids the tracking of the interface by introducing an auxiliary field that locates the interface and whose dynamics is coupled to the other physical fields through an appropriate set of partial differential equations. In comparison to the more traditional boundary integral methods, the method is much simpler to implement numerically.

In this Brief Report we introduce a phase field model for intracellular dynamics i.e. diffusion inside and outside a stationary, closed domain with source terms at the interface. We apply the method to a recently introduced model for the response of a Dictyostelium amoeba bill  following stimulation with the chemoattractant cAMP rltl02 . In rltl02 , due to the need to use a finite element method the numerical implementation of the model was limited to two space dimensions and the cells were treated as disks. As we will see below, the phase field method is capable of faithfully capturing no-flux boundary conditions. Thus, it becomes feasible to investigate more realistic cell shapes in three dimensions.

Before we introduce our approach, we like to point out some possible extensions of our methodology. The phase field approach can be easily modified to include problems where the domain boundary is not stationary. For example, force generation on cell membranes, leading to shape changes, can be incorporated within the phase field approach in a straightforward manner. This would require adding an additional equation for the phase field, but does not require the explicit calculation of a boundary. We should stress that attempting to model this type of problem using conventional techniques, with explicit boundary tracking, become quite cumbersome.

Let us first introduce the most salient ingredients of our method. Our purpose is to describe the situation where some chemoattractant diffuses only inside the cell, i.e. its concentration $`c`$ obeys the diffusion equation inside a closed domain:

$$\frac{c}{t}=D\stackrel{}{}^2c$$

and satisfies zero-flux boundary conditions:

$$\widehat{n}\stackrel{}{}c=0$$

As a phase field we take a function that takes the value $`\varphi _{in}=1`$ inside the cell and $`\varphi _{out}=0`$ outside the domain and varies smoothly across the interface. A possible choice in one dimension for a cell between $`x=a`$ and $`x=a`$ is:

$$\varphi (x)=\frac{1}{2}+\frac{1}{2}\mathrm{tanh}((a|x|)/\xi )$$
(1)
The variable $`\xi `$ denotes the interface width. We then define a field $`u`$ that is to obey the equation:

$$\varphi \frac{u}{t}=\stackrel{}{}\left[D\varphi \stackrel{}{}u\right]$$
(2)
Our claim is now that inside the domain the field $`u`$ behaves very similarly to $`c`$. Let us therefore first show, for simplicity in one dimension, that in the thin interface limit one recovers the no-flux boundary condition. Integration of (2) over the interface yields:

$$_{a\xi }^{a+\xi }𝑑x\varphi \frac{u}{t}D\frac{u}{x}|_{x=a\xi }$$

since $`\varphi (a\xi )1`$ and $`\varphi (a+\xi )0`$. From this we deduce:

$$\frac{u}{x}|_{x=a}\xi $$

and thus, in the sharp interface limit $`\xi 0`$, the reflective boundary conditions are recovered. It is also clear that inside the domain where $`\varphi `$ is constant $`u`$ satisfies the diffusion equation.

Another interesting property of Eq. (2) is that the quantity $`\varphi u𝑑\stackrel{}{x}`$, which can be interpreted as the total concentration inside the cell, is conserved under the dynamics. Thus, even if the field $`u`$ may become non-zero outside the cell (and indeed it does), the total concentration inside the cell remains constant. In fact, the asymptotic solution of Eq. (2) (together with phase field Eq. (1)) is $`u(x)=A`$ where the constant $`A`$ is determined through $`2aA=𝑑x\varphi u(t=0)`$. The real asymptotic solution is of course $`u(x)=\stackrel{~}{A}=_a^a𝑑xu(t=0)/(2a)`$. Because $`\varphi =1`$ inside the domain and $`\varphi =0`$ outside of it, the error $``$ in the asymptotic solution depends solely on the amount of concentration initially near the interfaces: $`_{\pm a\xi }^{\pm a+\xi }𝑑x\varphi u(t=0)/2a`$. Since $`\varphi `$ is an antisymmetric function around the interface, the error is minimal if we take in the phase field approach an initial condition $`u(t=0)`$ which is locally symmetric around the boundary.

This raises the following important point. While the solution $`u`$ outside the cell is not of physical interest, it is essential to keep track of it. In practice we solve (2) where the phase field $`\varphi `$ exceeds a small threshold $`\delta `$ (typically $`\delta 10^8`$).

As a first test for our model we now solve Eq. (2) numerically in two space dimensions. For the phase field we take here:

$$\varphi (r)=\frac{1}{2}+\frac{1}{2}\mathrm{tanh}((r_0r)/\xi )$$
(3)
where $`r`$ is the radius in polar coordinates. In the case of a radially symmetric initial condition we can compare the field $`u`$ to the analytic solution $`v`$ which is expressed in terms of the Bessel functions

$$v(r,t)=\underset{n}{}a_nJ_n(\lambda _nr)\text{e}^{\lambda _n^2t}$$

where $`\lambda _n`$ is defined as the smallest number for which $`J_n(\lambda _na)=0`$ and $`a_n`$ is the projection of the initial condition on the set of Bessel functions. We also compare both fields $`u`$ and $`v`$ with the solution $`w`$ of a finite element method obtained with MATLAB’s PDE Toolbox (The Mathworks, Natick, MA).

In Fig. 1a we show the spatial profile of the phase-field together with the exact solution at a given time. It can be seen that the agreement is excellent, with a relative error of less than $`1\%`$ (see insert). In Fig. 1b we show the fields $`u`$,$`v`$ and $`w`$ at three different points. Again the agreement is excellent.

In view of these promising results, let us now consider diffusion in presence of a production term and damping. Including a source term at the interface is relatively easy. We add to the right hand side of (2) a term $`b(\stackrel{}{}\varphi )^2/𝑑\stackrel{}{x}(\stackrel{}{}\varphi )^2`$. The factor $`(\stackrel{}{}\varphi )^2`$ ensures that it only acts at the interface and the denominator is a normalization factor. To compensate for the production we also include a damping term of the form $`\mu \varphi u`$. The equation of motion then becomes:

$$\varphi \frac{u}{t}=\stackrel{}{}\left[D\varphi \stackrel{}{}u\right]\mu \varphi u+b\frac{(\stackrel{}{}\varphi )^2}{𝑑\stackrel{}{x}(\stackrel{}{}\varphi )^2}$$
(4)
We have compared the phase field method with these two supplementary ingredients with the finite element method, again in the case of a two-dimensional circular cell with $`r_0=5`$. As can be seen in Fig. 2 the result is excellent.

We now use our phase field method to solve a biological model for the response of a Dictyostelium amoeba following stimulation with the chemoattractant cAMP rltl02 . In recent experiments the establishment of an asymmetry within a few seconds after a rise of extracellular cAMP was demonstrated. The cAMP however diffuses rapidly around the cell and the applied signal is several orders of magnitude larger than the value required to elicit a response. This strongly suggests the presence of an inhibitory mechanism that suppresses responses at the back.

In rltl02  an abstract model for the initial response of the cell to the chemoattractant was proposed: it was supposed that the membrane can be characterized in terms of three states: quiescent (with density $`\rho _q`$), activated (with density $`\rho _a`$) and inhibited (with density $`\rho _i`$). Initially the entire membrane is in the quiescent state. As the cAMP reaches the front of the cell the membrane becomes activated at rate $`\alpha [\mathrm{cAMP}]`$ and an inhibitor, in rltl02  suggested to be cGMP, is produced (at rate $`\sigma _g\rho _a`$). The inhibitor diffuses toward the back of cell where it turns the membrane from quiescent to inhibited with rate $`\beta _r[\mathrm{cGMP}]`$. Furthermore both activated and inhibited state decay spontaneously to the quiescent state at rates $`\delta `$ and $`\beta _f`$ respectively. The equations for the membrane state variables are thus:

$`{\displaystyle \frac{\rho _q}{t}}`$ $`=`$ $`\alpha c\rho _q+\beta _f\rho _i\beta _rg\rho _q`$ (5)
$`{\displaystyle \frac{\rho _a}{t}}`$ $`=`$ $`\alpha c\rho _q\delta \rho _a`$ (6)
$`{\displaystyle \frac{\rho _i}{t}}`$ $`=`$ $`\beta _f\rho _i+\beta _rg\rho _q+\delta \rho _a`$ (7)
The reactants cGMP and cAMP diffuse respectively inside and outside the cell. At the membrane they satisfy zero-flux boundary conditions. There is an source term for the cGMP field that accounts for the production of cGMP at the interface. Both cAMP and cGMP fields are damped at rates $`\mu _c`$ and $`\mu _g`$. The phase-field method is thus well suited to solve the dynamical equations for the cAMP and cGMP concentrations. For the phase field corresponding to the cAMP (which diffuses outside the cell) we simply take the complement of $`\varphi `$ given by Eq. (3): $`\varphi _c=1\varphi `$. The equations for the membrane variables are solved on all lattice sites where $`(\stackrel{}{}\varphi )^2`$ exceeds a certain threshold, namely $`10^4`$.

We now present a comparison of the phase field approach and the results obtained with a finite element method in rltl02 . A circular cell of diameter 10 $`\mu `$m is placed in a square domain of 30 $`\mu `$m by 30 $`\mu `$m. The diffusion constant of cAMP and cGMP was taken to be identical: $`D_c=D_g=250\mu `$m2/s. The values of the other parameters can be read in the caption of Fig.3.

In order to mimic the asymmetric cAMP stimulus we maintain the cAMP concentration at a value well above threshold at the upper left corner of our domain. As expected from our earlier results the agreement of the cAMP fields, that solely diffuse around the cell, is excellent, see Fig. 3a. We observe a slight discrepancy for the cGMP field (also Fig. 3a), which grows with time. This is related to the larger difference (of around 10%) for the membrane variables, the production of cGMP being proportional to $`\rho _a`$. The origin of the discrepancy might lie in the way the state variables are calculated: on a ring of finite width in the phase field model, and on 40 points on the interface in the finite element method. At any rate, since the model is of a quite abstract nature and since its predictions are only qualitative, a detailed investigation is beyond the scope of this paper.

From a computational perspective, we note that in the phase field method the CPU time required is linear in the number of lattice sites whereas it is quadratic in the number of nodes in the finite element method. One thus expects the phase field method to be much faster. However this naive result must be taken with some caution as in the finite element method the nodes are not distributed uniformly in space. To resolve some particular region in space with high accuracy we are thus obliged to take a small lattice spacing in the phase field method, such that the number of lattice points exceeds the number of nodes. Other parameters that will affect the accuracy of both methods are the time step $`\mathrm{\Delta }t`$ and the integration algorithm. Again a detailed comparison is beyond the scope of the Report. In practice it turns out that the phase field method where the lattice spacing is equal to the minimal distance between two nodes of our finite element method is about one order of magnitude faster.

Finally, we turn to the chemotaxis model in three space dimensions. For the sake of simplicity we take a spherical cell, i.e. with a phase field like in Eq. (3) but where $`r`$ is now the radius in spherical coordinates. We consider again a cell of radius $`r_0=5\mu m`$, in a box of dimensions $`30\times 30\times 30\mu m`$. Now the stimulus is applied by setting the cAMP initially well above threshold at one side of the box, here taken to be $`x=15\mu m`$. As is illustrated in Fig. 4, the phenomenology of the two dimensional case is reproduced here. As the cAMP front progresses toward the back of the cell, only the front of the cell is activated.

In conclusion we have proposed a phase field model for intracellular dynamics. Our method is shown to be very accurate, easy to implement and computationally inexpensive. Another advantage lies in its much greater flexibility with respect to other methods, like the finite element method. We revisited a chemotaxis model and, when considering three dimensional cells, we reach the same conclusions as in rltl02  for two-dimensional cells. Whereas we have restricted ourselves to circular and spherical domains, the extension to other geometrical forms poses no major problems, the only task being to generate a phase field $`\varphi `$. Even better, our approach can easily be extended to deal with non-stationary boundaries. This situation arises in a multitude of biological problems. For example, the Dictyostelium cells change their shape continuously during chemotaxis. Another example are shape transformations observed in vesicles seifert . In these cases, the phase field becomes a dynamic variable, that evolves under the appropriate Ginzburg-Landau type of equation. We are currently working along these lines of research.

This work was supported by the NSF sponsored Center for Theoretical Biological Physics (grants PHY-0216576 and 0225630) and the NSF Biocomplexity program (grant MCB 0083704).


text:
# Cosmology of Nonlinear Oscillations

## Acknowledgements

We would like to thank D.K. Hong and A. Zee for useful comments. This work was supported in part under DOE contract DE-FG06-85ER40224.


text:
# Fast proof of functional equation for 𝜁(𝑠)

(Date: 13 May 2003)

The functional equation for $`\zeta (s)`$ seems certain to remain a most important topic. More than a few dozen proofs for it are already known; perhaps this is a lower bound only. In Récoltes et Semailles Grothendieck admonishes mathematicians to look anew at old concepts in solitude and in absolute, childlike innocence. Following his dictum I found a very short proof of the functional equation, only to realize later -one expects to be forgiven- that it was an even shorter version of Titchmarsh’s ( p. 15) first proof of the functional equation based on the fundamental Fourier transform view of $`\zeta (s)`$ that appears, easily overlooked, on page 14 formula (2.1.5) of his treatise. Namely, with $`\rho (x)`$ standing for the fractional part of $`x`$,

(0.1)
$$\frac{\zeta (s)}{s}=_0^{\mathrm{}}x^{s1}\rho (x)𝑑x,(0<\mathrm{}s<1),$$

which, by the way, is the basis for the Nyman-Beurling approach ,  to the Riemann hypothesis.

Starting from (0.1) we can write

$`(2^s1){\displaystyle \frac{\zeta (s)}{s}}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}x^{s1}(\rho (x)\rho (2x))𝑑x`$
$`=`$ $`{\displaystyle _0^{\mathrm{}}}x^{s1}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n\pi }}(\mathrm{sin}4n\pi x\mathrm{sin}2n\pi x)dx.`$
valid for $`0<\mathrm{}s<1`$. On account of the uniform boundedness of the partial sums we can interchange integral and sum to obtain

$`(2^s1){\displaystyle \frac{\zeta (s)}{s}}`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n\pi }}{\displaystyle _0^{\mathrm{}}}x^{s1}(\mathrm{sin}4n\pi x\mathrm{sin}2n\pi x)𝑑x`$
$`=`$ $`(2^s1)2^s\pi ^{s1}\mathrm{\Gamma }(s)\mathrm{sin}{\displaystyle \frac{\pi s}{2}}\zeta (1s),`$
where the last equation is valid at first when $`1<\mathrm{}s<0`$, but it obviously provides the analytic extension to the whole plane.

Luis Báez-Duarte
Departamento de Matemáticas
Instituto Venezolano de Investigaciones Científicas
Apartado 21827, Caracas 1020-A, Venezuela