text:
# Scaling laws of creep rupture of fiber bundles

## I introduction

Under high steady stresses fiber composites may undergo time dependent deformation resulting in failure called creep rupture which limits their lifetime, and hence, has a high impact on the applicability of these materials in construction elements. Both natural fiber composites like wood  and various types of fiber reinforced composites  show creep rupture phenomena, which have attracted continuous theoretical and experimental interest over the past years. The underlying microscopic failure mechanism of creep rupture is very complex depending on several characteristics of the specific types of materials, and is far from being well understood. Theoretical studies encounter various challenges: on the one hand, applications of fiber composites require the development of analytical and numerical models which are able to predict the damage histories of loaded composites in terms of the specific parameters of constituents. On the other hand, creep rupture, similarly to other rupture phenomena, presents a very interesting problem for statistical physics, it is still an open problem to embed creep rupture into the general framework of statistical physics and to understand the analogy between rupture phenomena and phase transitions.

Creep failure tests are usually performed under uniaxial tensile loading when the specimen is subjected either to a constant load $`\sigma _\mathrm{o}`$ or to an increasing load (ramp-loading) and the time evolution of the damage process is monitored by recording the strain $`\epsilon `$ of the specimen and the acoustic signals emitted by microscopic failure events . In the present paper we study the creep rupture of fiber composites in the framework of fiber bundle models  focusing on general characteristics of creep rupture. Our main goal is to reveal the universal aspects of a creeping system which do not depend on specific material properties. We will derive scaling laws which emerge when approaching the failure stress both from below and above, furthermore, the finite size scaling of the time to failure is studied. Beyond the theoretical understanding, exploring universal features of creep helps to evaluate experimental data and can serve as a guide to extract the relevant information invariant under the variation of specific material properties.

We consider only the case of global load sharing (GLS) for the redistribution of load following fiber failure . For several types of materials GLS provides an adequate approach, and it has the advantage that many of the important GLS results can be obtained in closed analytic forms. During creep at the critical point and above it, stress localization occurs which can be captured by localized load sharing in fiber bundle models  but this is beyond the scope of our investigations. The present study is focused on analytic modeling in the framework of GLS. We also develop efficient simulation techniques for testing the analytic results, furthermore, to study finite size system. Since similar studies with local load sharing require large scale computer simulations and a large amount of data processing it will be presented in a forthcoming publication.

## II models

In order to work out a theoretical description of creep failure we introduce two novel fiber bundle models improving the classical fiber bundle models  which have proven very successful in the study of fracture of disordered materials . The two models differ in the microscopic mechanisms assumed to be responsible for the macroscopic creep behavior: $`(i)`$ in the first approach the fibers themselves are viscoelastic and they break when their deformation exceeds a stochastically distributed threshold value (strain controlled breaking), $`(ii)`$ in the second model the fibers are linearly elastic until they break stochastically in a stress controlled way, however, after breaking their relaxation is not instantaneous but the sliding of the broken fiber with respect to the matrix material, and the yielding or creeping of the matrix introduces an intrinsic time scale for the relaxation.

The mathematical formalism of the first model (viscoelastic bundle) is more simple, hence, the steps of analysis will be presented in details for this case. For the second model the corresponding results will be summarized.

### A Viscoelastic fiber bundle

#### 1 Analytic model

Our model consists of $`N`$ parallel fibers having viscoelastic constitutive behavior. For simplicity, the pure viscoelastic behavior of fibers is modeled by a Kelvin-Voigt element which consists of a spring and a dashpot in parallel (see Fig. 1) and results in the constitutive equation

$`\sigma _\mathrm{o}=\beta \dot{\epsilon }+E\epsilon ,`$ (1)
where $`\beta `$ denotes the damping coefficient, and $`E`$ the Young modulus of fibers, respectively. Eq. (1) provides the time dependent deformation $`\epsilon (t)`$ of a fiber at a fixed external load $`\sigma _\mathrm{o}`$

$`\epsilon (t)={\displaystyle \frac{\sigma _o}{E}}\left[1e^{Et/\beta }\right]+\epsilon _\mathrm{o}e^{Et/\beta },`$ (2)
where $`\epsilon _o`$ denotes the initial strain at $`t=0`$. It can be seen that $`\epsilon (t)`$ converges to $`\sigma _\mathrm{o}/E`$ for $`t\mathrm{}`$, which implies that the asymptotic strain fulfills Hook’s law. To incorporate breaking in the model we introduce a strain controlled failure criterion for fibers: a fiber fails during the time evolution of the system if its strain exceeds a damage threshold $`\epsilon _d`$, which is an independent identically distributed random variable of fibers with probability density $`p(\epsilon _d)`$ and cumulative distribution $`P(\epsilon _d)={\displaystyle _0^{\epsilon _d}}p(x)𝑑x`$. Due to the validity of Hook’s law for the asymptotic strain values, the formulation of the failure criterion in terms of strain instead of stress implies that under a certain steady load the same amount of damage occurs as in the case of stress controlled failure, however, the breaking of fibers is not instantaneous but distributed over time. When a fiber fails its load has to be redistributed to the intact fibers. Assuming global load sharing, the time evolution of the system under a steady external load $`\sigma _o`$ is finally described by the equation

$`{\displaystyle \frac{\sigma _\mathrm{o}}{1P(\epsilon )}}=\beta \dot{\epsilon }+E\epsilon ,`$ (3)
where the viscoelastic behavior of fibers is coupled to the failure of fibers in a global load sharing framework .

For the behavior of the solutions of Eq. (3) two distinct regimes can be distinguished depending on the value of the external load $`\sigma _\mathrm{o}`$: When $`\sigma _\mathrm{o}`$ is below a critical value $`\sigma _\mathrm{c}`$ Eq. (3) has a stationary solution $`\epsilon _s`$, which can be obtained by setting $`\dot{\epsilon }=0`$ in Eq. (3)

$`\sigma _\mathrm{o}=E\epsilon _s[1P(\epsilon _s)].`$ (4)
It means that until this equation can be solved for $`\epsilon _s`$ at a given external load $`\sigma _\mathrm{o}`$, the solution $`\epsilon (t)`$ of Eq. (3) converges to $`\epsilon _s`$ when $`t\mathrm{}`$, and no macroscopic failure occurs. However, when $`\sigma _\mathrm{o}`$ exceeds the critical value $`\sigma _c`$ no stationary solution exists, furthermore, $`\dot{\epsilon }`$ remains always positive, which implies that for $`\sigma >\sigma _c`$ the strain of the system $`\epsilon (t)`$ monotonically increases until the system fails globally at a time $`t_f`$ .

In the regime $`\sigma _\mathrm{o}\sigma _\mathrm{c}`$ Eq. (4) also provides the asymptotic constitutive behavior of the fiber bundle which can be measured by controlling the external load $`\sigma _\mathrm{o}`$ and letting the system relax to $`\epsilon _s`$. It follows from the above argument that the critical value of the load $`\sigma _c`$ is the static fracture strength of the bundle which can be determined from Eq. (4) as $`\sigma _c=E\epsilon _c[1P(\epsilon _c)]`$, where $`\epsilon _c`$ is the solution of the equation $`d\sigma _\mathrm{o}/d\epsilon _s|_{\epsilon _\mathrm{c}}=0`$ . Since $`\sigma _\mathrm{o}(\epsilon _\mathrm{s})`$ has a maximum of the value $`\sigma _\mathrm{c}`$ at $`\epsilon _\mathrm{c}`$, in the vicinity of $`\epsilon _\mathrm{c}`$ it can be approximated as

$`\sigma _\mathrm{o}\sigma _\mathrm{c}A(\epsilon _\mathrm{c}\epsilon _s)^2,`$ (5)
where the multiplication factor $`A`$ depends on the probability distribution $`P`$. A complete description of the system can be obtained by solving the differential equation Eq. (3). After separation of variables the integral arises

$`t=\beta {\displaystyle 𝑑\epsilon \frac{1P(\epsilon )}{\sigma _\mathrm{o}E\epsilon \left[1P(\epsilon )\right]}}+C,`$ (6)
where the integration constant $`C`$ is determined by the initial condition $`\epsilon (t=0)=0`$.

The creep rupture of the viscoelastic bundle can be interpreted so that for $`\sigma _\mathrm{o}\sigma _\mathrm{c}`$ the system suffers only a partial failure which implies an infinite lifetime $`t_f=\mathrm{}`$ and the emergence of a macroscopic stationary state, while above the critical load $`\sigma _\mathrm{o}>\sigma _\mathrm{c}`$ global failure occurs at a finite time $`t_f`$, which can be determined by evaluating the integral Eq. (6) over the whole domain of definition of $`P(\epsilon )`$.

Below the critical point $`\sigma _\mathrm{o}\sigma _\mathrm{c}`$ the bundle relaxes to the stationary deformation $`\epsilon _s`$ through a decreasing breaking activity. To find the characteristic time scale of this relaxation process the behavior of $`\epsilon (t)`$ has to be analyzed in the vicinity of $`\epsilon _s`$. It is useful to introduce a new variable $`\delta `$ as $`\delta (t)=\epsilon _s\epsilon (t)`$. The governing differential equation of $`\delta `$ can be obtained from Eq. (3) by expanding it around $`\epsilon _s`$

$`{\displaystyle \frac{d\delta }{dt}}={\displaystyle \frac{E}{\beta }}\left[1{\displaystyle \frac{\epsilon _sp(\epsilon _s)}{1P(\epsilon _s)}}\right]\delta .`$ (7)
The solution of Eq. (7) has the form $`\delta \mathrm{exp}\left[t/\tau \right]`$, where $`\tau `$ is the characteristic time scale of the relaxation process

$`\tau ={\displaystyle \frac{\beta }{E}}{\displaystyle \frac{1}{\left[1{\displaystyle \frac{\epsilon _sp(\epsilon _s)}{1P(\epsilon _s)}}\right]}}.`$ (8)
It is a very important question how the relaxation time $`\tau `$ changes when the external driving approaches the critical point $`\sigma _\mathrm{c}`$ from below. Based on Eq. (5) it can be simply shown that

$`\tau \left(\sigma _c\sigma _\mathrm{o}\right)^{1/2},\text{for}\sigma _o<\sigma _\mathrm{c},`$ (9)
which means that approaching the critical point from below the relaxation time of the system diverges according to a universal power law with an exponent $`1/2`$ independent on the form of disorder distribution. Note that a similar power law divergence of the number of successive relaxation steps was found in Refs.  for a dry fiber bundle subjected to a constant load.

Above the critical point the behavior of the lifetime of the system can be analyzed analogously when $`\sigma _\mathrm{o}`$ goes to $`\sigma _\mathrm{c}`$ from above. When $`\sigma _\mathrm{o}`$ is in the vicinity of $`\sigma _\mathrm{c}`$, i.e. $`\sigma _\mathrm{o}=\sigma _\mathrm{c}+\mathrm{\Delta }\sigma _\mathrm{o}`$, where $`\mathrm{\Delta }\sigma _\mathrm{o}<<\sigma _\mathrm{c}`$, it can be expected that the curve of $`\epsilon (t)`$ falls very close to $`\epsilon _\mathrm{c}`$ for a very long time and the breaking of the system occurs suddenly. Hence, the total time to failure, i.e. the integral in Eq. (6), is dominated by the region close to $`\epsilon _\mathrm{c}`$ when $`\mathrm{\Delta }\sigma _\mathrm{o}`$ is small. Making use of the power series expansion Eq. (5) the integral in Eq. (6) can be rewritten as

$`t_f\beta {\displaystyle 𝑑\epsilon \frac{1P(\epsilon )}{\mathrm{\Delta }\sigma _\mathrm{o}A(\epsilon _\mathrm{c}\epsilon )^2}},`$ (10)
which has to be evaluated over a small $`\epsilon `$ interval in the vicinity of $`\epsilon _\mathrm{c}`$. After performing the integration it follows

$`t_f(\sigma _\mathrm{o}\sigma _\mathrm{c})^{1/2},\text{for}\sigma _\mathrm{o}>\sigma _\mathrm{c}.`$ (11)
Thus, $`t_f`$ has a power law divergence at $`\sigma _\mathrm{c}`$ with a universal exponent $`\frac{1}{2}`$ independent of the specific form of the disorder distribution $`P(\epsilon )`$, similarly to $`\tau `$ below the critical point.

#### 2 Finite size effect

In the above analytic treatment the size of the system, i.e. the number of fibers in the bundle, is infinite. However, it can be expected that the lifetime of a finite bundle has a non-trivial size scaling even in the case of global load sharing. During the creep rupture process the fibers of a finite bundle break one-by-one in the increasing order of their breaking thresholds. Let $`\epsilon _j,j=1,\mathrm{},N`$ denote the breaking thresholds assigned to the fibers in a realization of the bundle. Since fibers break one-by-one, the actual load on the intact fibers after the failure of $`i`$ fibers is $`\sigma _i=\sigma _oN/(Ni)`$ where $`i=0,\mathrm{},N1`$, and the time $`\mathrm{\Delta }t(\epsilon _i,\epsilon _{i+1})`$ between the breaking of the $`i`$th and $`i+1`$th fibers (in the ordered series) reads as

$`\mathrm{\Delta }t(\epsilon _i,\epsilon _{i+1})={\displaystyle \frac{\beta }{E}}\left[\mathrm{ln}\left(\epsilon _{i+1}{\displaystyle \frac{\sigma _i}{E}}\right)\mathrm{ln}\left(\epsilon _i{\displaystyle \frac{\sigma _i}{E}}\right)\right].`$ (12)
The lifetime $`t_f`$ of a sample of $`N`$ fibers takes the form

$`t_f(N)={\displaystyle \underset{i=0}{\overset{N1}{}}}\mathrm{\Delta }t(\epsilon _i,\epsilon _{i+1}).`$ (13)
In order to determine how $`t_f`$ depends on $`N`$, Eq. (13) has to be averaged over many realizations of the disorder distribution, which can be performed analytically. The details of the analytic calculations are summarized in the Appendix. Finally, the average lifetime $`t_f(N)`$ of a bundle of $`N`$ fibers can be cast in the form

$`t_f(N)t_f(\mathrm{})+{\displaystyle \frac{\beta \sigma _o}{N}}{\displaystyle \frac{d\epsilon E\epsilon P(\epsilon )\left[1P(\epsilon )\right]}{\left[\sigma _oE\epsilon (1P(\epsilon ))\right]^3}}.`$ (14)
Eq. (14) shows that for finite bundles the average lifetime $`t_f(N)`$ converges to the lifetime of the infinite bundle $`t_f(\mathrm{})`$ as $`1/N`$ with increasing number of fibers $`N`$. It is interesting to note that in the case of global load sharing the average strength of the bundle $`\sigma _c`$ does not have a size dependence.

#### 3 Simulation technique

Most of the analytic results of the previous sections, except for the finite size scaling of $`t_f`$, were obtained for infinite bundles. Computer simulations of the creep rupture of finite bundles are needed to justify the validity of analytic predictions for finite systems, and to be able to model the rupture process of realistic finite systems.

In the framework of GLS an efficient simulation technique can be worked out for the failure process. Based on the arguments of the previous subsection, the GLS simulation of the creep process of a bundle of $`N`$ fibers proceeds as follows: $`(i)`$ random breaking thresholds $`\epsilon _i,i=1,\mathrm{},N`$ are drawn from a probability distribution $`p`$, then the thresholds are put into increasing order. $`(ii)`$ The time $`\mathrm{\Delta }t(\epsilon _i,\epsilon _{i+1})`$ between the breaking of the $`i`$th and $`i+1`$th fibers is calculated according to Eq. (12). $`(iii)`$ The time elapsed till the breaking of the $`i`$th fiber is obtained as $`t(\epsilon _i)=_{j=0}^{i1}\mathrm{\Delta }t(\epsilon _i,\epsilon _{i+1})`$, from which the deformation as a function of time $`\epsilon (t)`$ can be determined by inversion. The lifetime $`t_f`$ of the finite bundle can be obtained making use of Eq. (13).

For simulations we considered a uniform distribution of failure thresholds between $`0`$ and a maximum value $`\epsilon _m`$ with the probability density function $`p(\epsilon )=1/\epsilon _m`$, and distribution function $`P(\epsilon )=\epsilon /\epsilon _m`$. In this case the stationary solution, the critical load and the corresponding critical strain can be obtained as $`\sigma _\mathrm{o}=E\epsilon [1\epsilon /\epsilon _m]`$, $`\sigma _c=E\epsilon _m/4`$, $`\epsilon _c=\epsilon _m/2`$, respectively. $`\epsilon _m=1`$ was set in all the simulations. The deformation-time diagram $`\epsilon (t)`$ obtained by simulations is presented in Fig. 2 for several different values of $`\sigma _o`$ below and above $`\sigma _c`$. The two regimes depending on the value of the external load $`\sigma _o`$ can be clearly distinguished.

To test the validity of the universal power law behavior of $`t_f`$ as a function of the distance from the critical load given by Eq. (11), simulations were performed with various disorder distributions, i.e. besides the uniform distribution the Weibull distribution of the form $`P(\epsilon )=1\mathrm{exp}\left[\left(\epsilon /\lambda \right)^\rho \right]`$ was employed. The value of the characteristic strain $`\lambda `$ was set to one, and the shape of the distribution was controlled by varying the value of $`\rho `$. The results are presented in Fig. 3, where an excellent agreement of the simulations and the analytic results can be observed. Fig. 3 supports that the exponent of $`t_f`$ as a function of $`\sigma _o\sigma _c`$ is universal, it does not depend on the specific form of the disorder distribution.

To study the finite size scaling of the time to failure $`t_f`$ a uniform distribution was used for the failure thresholds. The value of the external load was fixed above $`\sigma _c`$ and the number of fibers $`N`$ was varied from $`5\times 10^2`$ to $`10^7`$. Averages were calculated over $`10^4`$ samples for each system size $`N`$. The results obtained by simulations are presented in Fig. 4, where an excellent agreement of simulations and analytic results can be observed for four orders of magnitude in the system size $`N`$.

### B Slowly relaxing fibers

#### 1 Analytic model

Another important microscopic mechanism which can lead to macroscopic creep is the slow relaxation following fiber failure. In this case, the components of the solid are linearly elastic until they break, however, after breaking they undergo a slow relaxation process, which can be caused, for instance, by the sliding of broken fibers with respect to the matrix material or by the creeeping matrix. To take into account this effect, our approach is based on the model introduced in Refs. , where the response of a viscoelastic-plastic matrix reinforced with elastic and also viscoelastic fibers have been studied.

The model consists of $`N`$ parallel fibers, which break in a stress controlled way, i.e. subjecting a bundle to a constant external load fibers break during the time evolution of the system when the local load on them exceeds a stochastically distributed breaking threshold $`\sigma _i,i=1,\mathrm{},N`$. Intact fibers are assumed to be linearly elastic i.e. $`\sigma =E_\mathrm{f}\epsilon _\mathrm{f}`$ holds until they break, and hence, for the deformation rate it applies

$`\dot{\epsilon }_\mathrm{f}={\displaystyle \frac{\dot{\sigma }}{E_\mathrm{f}}}.`$ (15)
Here $`\epsilon _f`$ denotes the strain and $`E_f`$ is the Young modulus of intact fibers, respectively. The main assumption of the model is that when a fiber breaks its load does not drop to zero instantaneously, instead it undergoes a slow relaxation process introducing a time scale into the system. In order to capture this effect, the broken fibers with the surrounding matrix material are modeled by Maxwell elements as illustrated in Fig. 5, i.e. they are conceived as a serial coupling of a spring and a dashpot which results in a non-linear response

$`\dot{\epsilon }_\mathrm{b}={\displaystyle \frac{\dot{\sigma }_\mathrm{b}}{E_\mathrm{b}}}+B\sigma _\mathrm{b}^\mathrm{m},`$ (16)
where $`\sigma _b`$ and $`\epsilon _\mathrm{b}`$ denote the time dependent load and deformation of a broken fiber, respectively. The relaxation of the broken fiber is characterized by three parameters $`E_b,B,`$ and $`m`$, where $`E_b`$ is the effective stiffness of a broken fiber, and the exponent $`m`$ characterizes the strength of non-linearity of the element. We study the behavior of the system for the region $`m1`$.

Assuming global load sharing for the load redistribution, the constitutive equation describing the macroscopic elastic behavior of the composite reads as

$`\sigma _\mathrm{o}=\sigma (t)\left[1P(\sigma (t))\right]+\sigma _\mathrm{b}(t)P(\sigma (t)).`$ (17)
Eq. (17) takes into account that broken fibers carry also a certain amount of load $`\sigma _b(t)`$, furthermore, $`P(\sigma (t))`$ and $`1P(\sigma (t))`$ denote the fraction of broken and intact fibers at time $`t`$, respectively . It can be seen from Eq. (17) that under a constant external load $`\sigma _0`$, the load of intact fibers $`\sigma `$ will also be time dependent due to the slow relaxation of the broken ones.

Due to the boundary condition illustrated in Fig. 5, the two time derivatives have to be always equal

$`\dot{\epsilon }_\mathrm{f}=\dot{\epsilon }_\mathrm{b}.`$ (18)
The differential equation governing the time evolution of the system can be obtained by expressing $`\sigma _\mathrm{b}`$ in terms of $`\sigma `$ from Eq. (17) and substituting it into Eq. (16) and finally into Eq. (18)

$`\dot{\sigma }\left\{{\displaystyle \frac{1}{E_\mathrm{f}}}{\displaystyle \frac{1}{E_\mathrm{b}}}\left[1{\displaystyle \frac{1}{P(\sigma )}}+{\displaystyle \frac{p(\sigma )}{P(\sigma )^2}}(\sigma \sigma _\mathrm{o})\right]\right\}=`$ (19)
$`B\left[{\displaystyle \frac{\sigma _\mathrm{o}\sigma \left[1P(\sigma )\right]}{P(\sigma )}}\right]^m.`$ (20)
In order to determine the initial condition for the integration of Eq. (19) the breaking process of fibers has to be analyzed. Subjecting the undamaged specimen to an external stress $`\sigma _o`$ all the fibers attain this stress value immediatelly due to the linear elastic response. Hence the time evolution of the system can be obtained by integrating Eq. (19) with the initial condition $`\sigma (t=0)=\sigma _o`$. Since intact fibers are linearly elastic, the deformation-time history $`\epsilon (t)`$ of the model can be deduced as $`\epsilon (t)=\sigma (t)/E_f`$, which has an initial jump to $`\epsilon _o=\sigma _o/E_f`$. It follows that those fibers which have breaking thresholds $`\sigma _i`$ smaller than the externally imposed $`\sigma _o`$ immediately break.

To characterize the macroscopic behavior of the composite the solutions $`\sigma (t)`$ of Eq. (19) have to be analyzed at different values of the external load $`\sigma _o`$. Similarly to the previous model, two different regimes of $`\sigma (t)`$ can be distinguished depending on the value of $`\sigma _o`$: if the external load falls below a critical value $`\sigma _\mathrm{c}`$ a stationary solution $`\sigma _s`$ of the governing equation exists which can be obtained by setting $`\dot{\sigma }=0`$ in Eq. (19)

$`\sigma _o=\sigma _s\left[1P(\sigma _s)\right].`$ (21)
This means that until Eq. (21) can be solved for $`\sigma _s`$ the solution $`\sigma (t)`$ of Eq. (19) converges asymptotically to $`\sigma _s`$ resulting in an infinite lifetime $`t_f`$ of the composite. Note that Eq. (21) provides also the asymptotic constitutive behavior of the model which can be measured by quasistatic loading. If the external load falls above the critical value the deformation rate $`\dot{\epsilon }=\dot{\sigma }/E_f`$ remains always positive resulting in a macroscopic rupture in a finite time $`t_f`$. It follows from Eq. (21) that the critical load $`\sigma _c`$ of creep rupture coincides with the static fracture strength of the composite.

The behavior of the system shows again universal aspects in the vicinity of the critical point. Below the critical point the relaxation of $`\sigma (t)`$ to the stationary solution $`\sigma _s`$ is governed by a differential equation of the form

$`{\displaystyle \frac{d\delta }{dt}}\delta ^m,`$ (22)
where $`\delta `$ denotes the difference $`\delta (t)=\sigma _s\sigma (t)`$. Hence, the characteristic time scale $`\tau `$ of the relaxation process only emerges if $`m=1`$, furthermore, in this case also $`\tau (\sigma _c\sigma _o)^{1/2}`$ holds when approaching the critical point. However, for $`m>1`$ the relaxation process is characterized by $`\delta (t)=at^{1/1m}`$, where $`a0`$ with $`\sigma _0\sigma _c`$.

Similarly to the previous model, it can also be shown that the lifetime $`t_f`$ of the bundle has a power law divergence when the external load approaches the critical point from above

$`t_f\left(\sigma _0\sigma _\mathrm{c}\right)^{(m1/2)},\text{for}\sigma _o>\sigma _\mathrm{c}.`$ (23)
The exponent is universal in the sense that it is independent on the disorder distribution, however, it depends on the stress exponent $`m`$, which characterizes the non-linearity of broken fibers.

#### 2 Simulation technique

Subjecting a finite bundle of $`N`$ fibers to $`\sigma _o`$ external stress those fibers whose failure threshold falls below $`\sigma _o`$ break immediately. The number $`N_o`$ of intially breaking fibers can be estimated from the disorder distribution as $`N_oNP(\sigma _o)`$. In the presence of broken fibers the system slows down and the remaining fibers of the bundle break one-by-one in the increasing order of their breaking thresholds $`\sigma _{N_o+1}<\sigma _{N_o+2}<\sigma _N`$. In order to construct an efficient simulation technique one has to determine the time elapsed between two consecutive breakings during the creep process.

The macroscopic constitutive equation for a system of $`N`$ fibers when $`i`$ fibers have already failed can be written as

$`\sigma _\mathrm{o}=\sigma {\displaystyle \frac{Ni}{N}}+\sigma _\mathrm{b}{\displaystyle \frac{i}{N}}.`$ (24)
Making use of Eqs. (16,18), the differential equation describing the time evolution of the load of intact fibers $`\sigma `$ can be cast in the form

$`\dot{\sigma }\left[{\displaystyle \frac{1}{E_f}}{\displaystyle \frac{1}{E_b}}\left(1{\displaystyle \frac{N}{i}}\right)\right]=B\left({\displaystyle \frac{N}{i}}\right)^mf_i(\sigma )^m,`$ (25)
where $`f_i(x)`$ is introduced for brevity as

$`f_i(x)=\sigma _o{\displaystyle \frac{Ni}{N}}x`$ (26)
The time $`\mathrm{\Delta }t(\sigma _i,\sigma _{i+1})`$ elapsed between the breaking of the $`i`$th and $`i+1`$th fiber can be determined by integrating Eq. (25) from $`\sigma _i`$ to $`\sigma _{i+1}`$, which yields for $`m1`$

$`\mathrm{\Delta }t(\sigma _i,\sigma _{i+1})={\displaystyle \frac{K_i}{(m1)}}\left[f_i(\sigma _{i+1})^{1m}f_i(\sigma _i)^{1m}\right],`$ (27)
and the multiplication factor $`K_i`$ reads as

$`K_i={\displaystyle \frac{N}{Ni}}\left({\displaystyle \frac{i}{N}}\right)^m{\displaystyle \frac{1}{B}}\left[{\displaystyle \frac{1}{E_f}}{\displaystyle \frac{1}{E_b}}\left(1{\displaystyle \frac{N}{i}}\right)\right].`$ (28)
For $`m=1`$ the corresponding equation has the form

$`\mathrm{\Delta }t(\sigma _i,\sigma _{i+1})=K_i\left[\mathrm{ln}f_i(\sigma _{i+1})\mathrm{ln}f_i(\sigma _i)\right].`$ (29)
Then the simulation proceeds as in the case of viscoelastic bundles but in the above formulas the number of broken fibers $`i`$ varies as $`i=N_o,N_o+1,\mathrm{},N1`$, so the time as a function of $`\sigma `$ can be obtained as $`t(\sigma _i)=_{j=N_o+1}^i\mathrm{\Delta }t(\sigma _j,\sigma _{j+1})`$ from which the deformation as a function of time $`\epsilon _i(t)`$ can be determined, since $`\sigma =E_f\epsilon _f`$ always holds. The lifetime $`t_f`$ of the system can be obtained by summing up all the $`\mathrm{\Delta }t`$’s.

For the purpose of explicit calculations a uniform distribution was prescribed for the breaking thresholds $`\sigma _i`$ between 0 and 1. The deformation as a function of time is plotted in Fig. 6 for several different values of the external load below and above the critical load. Similarly to the previous model the two regimes of the creeping system can be clearly distinguished.

To study the behavior of the time to failure as a function of the distance from the critical point, simulations were performed for several different values of the exponent $`m`$. In Fig. 7 the results are presented for $`m=1.5`$ and $`m=2.5`$. The slope of the fitted straight lines agrees very well with the analytic predictions of Eq. (23).

The size scaling of the time to failure $`t_f`$ was analyzed by simulating the creep rupture of bundles of size $`N=510^210^7`$ setting a uniform distribution for the breaking thresholds. We found that $`t_f(N)`$ converges to the lifetime of the infinite system $`t_f(\mathrm{})`$ according to the universal law Eq. (14) independently on the value of the exponent $`m`$. In Fig. 8 the best fit was obtained for both curves with slope $`1.0\pm 0.05`$ for both $`m`$ values.

## III Conclusions

The creep rupture of fibrous materials occurring under a steady external load is microscopically a rather complex phenomenon depending on a diversity of possible material specific mechanisms. Therefore, on the one hand, it is impossible to work out a general theoretical framework which takes into account all the features of the process and has predictive power, and on the other hand, it is very important to reveal universal aspects of the creep process, which do not depend on specific material properties relevant at the micro level.

In the present paper we studied the creep rupture of fibrous materials in the framework of fiber bundle models taking into account two possible microscopic mechanisms of creep: $`(i)`$ in the first approach the fibers themselves are viscoelastic and they break when their deformation exceeds a stochastically distributed threshold value, $`(ii)`$ in the second model the fibers are linearly elastic until they break, however, after breaking their relaxation is not instantaneous but the creeping matrix introduces an intrinsic time scale for the relaxation. The first model can be relevant for natural fiber composites like wood which are composed of viscoelastic fibers , while the second model can provide an adequate description of metal matrix composites reinforced by brittle fibers . Analytical and numerical calculations showed in both models that increasing the external load on a specimen a transition takes place from a partially failed state of infinite lifetime to a state where global failure occurs at a finite time. The critical load turned to be the static fracture strength of the material. It was found that although the details of the two model systems are different, in the vicinity of the critical point they exhibit universal behaviour which can be relevant for experiments.

## IV Acknowledgment

This work was supported by the project SFB381, by the NATO grant PST.CLG.977311 and by OTKA T037212. F. Kun acknowledges financial support of the Bólyai János Foundation of the Hungarian Academy of Sciences and of the Research Contract FKFP 0118/2001.

## V Appendix

Here we provide the derivation of the average lifetime for the general case when the lifetime $`t_f`$ of a bundle with a specific realization of the disorder can be cast in the form

$`t_f={\displaystyle \underset{i=0}{\overset{N1}{}}}\left[G({\displaystyle \frac{i}{N}},x_{i+1})G({\displaystyle \frac{i}{N}},x_i)\right],`$ (30)
i.e. $`t_f`$ is a sum of terms which depend on the number of broken fibers $`i`$ and on a single breaking threshold $`x_i`$ that can be given as strain or stress. The $`x_i`$-s are obtained by choosing $`N`$ breaking thresholds independently from a cumulative probability distribution $`P(x)`$ and putting them into increasing order. This treatment includes both models discussed in the present paper. The expectation value of a function $`f(x_i)`$ can be determined as

$`f(x_i)`$ $`=`$ $`{\displaystyle \frac{N!}{(i1)!(Ni)!}P(x)^{i1}\left[1P(x)\right]^{Ni}}`$ (32)
$`\times p(x)f(x)dx.`$
The probability distribution in Eq. (32) that the value of the $`i`$th largest breaking threshold falls between $`x`$ and $`x+dx`$ has a sharp peak for large $`N`$ values for each $`i`$. The above integration can be carried out by expanding the distribution about its peak. After expanding the result in terms of $`1/N`$ and neglecting higher order terms we arrive at

$`f(x_i)`$ $`=`$ $`f(\overline{x}_i)+{\displaystyle \frac{1}{2N}}P(\overline{x}_i)(1P(\overline{x}_i))`$ (34)
$`\times \left[{\displaystyle \frac{f^{^{\prime \prime }}(\overline{x}_i)}{[P^{^{}}(\overline{x}_i)]^2}}{\displaystyle \frac{f^{^{}}(\overline{x}_i)P^{^{\prime \prime }}(\overline{x}_i)}{[P^{^{}}(\overline{x}_i)]^3}}\right],`$
where $`\overline{x}_i`$ is defined implicitly by $`P(\overline{x}_i)=i/(N+1)`$. Applying Eq. (34) to Eq. (30) the resulting summation can be approximated by integrals replacing $`i/N`$ by the equivalent $`P(\overline{x}_i)(1+1/N)`$. Neglecting corrections higher order in $`1/N`$ after straightforward calculations we arrive at

$`t_f{\displaystyle 𝑑x_2G(P(x),x)}`$ (35)
$`+{\displaystyle \frac{1}{2N}}{\displaystyle 𝑑xP(x)\left(1P(x)\right)_1^2_2G(P(x),x)},`$ (36)
where $`_2G(y,x)G(y,x)/x`$, and $`_1^2_2G(y,x)^3G(y,x)/^2yx`$. Substituting the actual form of $`G(y,x)`$ for a specific model the complete form of the size scaling of lifetime can be obtained. However, it can be seen in the general expression Eq. (35) that the first term provides the lifetime of the infinite bundle and the only size dependence is in the prefactor of the second term. Eq. (35) states that if the lifetime can be written in the form of Eq. (30) the lifetime of finite bundles converges to that of the infinite one as $`1/N`$ with increasing number of fibers $`N`$.


format: text

text:
# Towards a new quantization of Dirac’s monopole

## Abstract

There are several mathematical and physical reasons why Dirac’s quantization must hold. How far one can go without it remains an open problem. The present work outlines a few steps in this direction.

magnetic monopole, nonunitary representations, anyons

In his pioneer work Dir  Dirac showed that the existence of just one magnetic monopole would explain electric charge quantization. Nowadays this is known as Dirac’s quantization rule. In sixties Lipkin et al Lip  showed that in the presence of the magnetic monopole the Jacobi identity for the translational group failed. Later Jackiw Jac  found that the Dirac’s quantization restores the associativity of the finite translations, but the infinitesimal generators remain “nonassociative” and the Jacobi identity fails.

There exist various arguments based on the quantum mechanics, theory of representations, topology and differential geometry in behalf of the Dirac’s rule Wu1 ; Wu2 ; Jac ; GS ; Nsh ; K . All of them are related to the necessary conditions of saving the situation and the problem of finding the sufficient conditions of existing of magnetic point like monopole is still open. It is therefore of interest to construct self consistent nonasociative quantum mechanics with an arbitrary magnetic monopole charge or maybe involving the distinct quantization rules related to the magnetic monopole. The present work outlines a few steps in this direction and the problems we have found.

Let us consider a point particle with electric charge $`e`$ and mass $`m`$ moving in the field of a magnetic monopole of charge $`q`$. The non-relativistic classical Hamiltonian may be written as follows Gol

$$H=\frac{1}{2mr^2}(𝐩𝐫)^2+\frac{1}{2mr^2}(𝐉^2\mu ^2)$$

where

$$𝐉=𝐫\times (𝐩e𝐀)\mu \frac{𝐫}{r},\mu =eq$$
(1)
is the conserved total angular momentum and we set $`\mathrm{}=c=1`$. The last term in the formula (1) usually is interpreted as the contribution of the electromagnetic field Gol ; Wil ; Gol2 ; Lyn , which carries an angular momentum

$$𝐋_{em}=\frac{1}{4\pi }𝐫\times (𝐄\times 𝐁)d^3r=\mu \frac{𝐫}{r}$$

Also this is known as Poincaré magnetic angular momentum Ber . A magnetic field of the Dirac’s monopole is

$$𝐁=q\frac{𝐫}{r^3}$$
(2)
and any choice of the vector potential $`𝐀`$ being compatible with (2) must have singularities, the so-called Dirac’s string.

At the quantum level the operator

$$𝐉=𝐫\times \left(ie𝐀\right)\mu \frac{𝐫}{r}$$
(3)
representing the angular momentum $`𝐉`$ has the same properties as a standard angular momentum and obeys the following commutation relations

$`[H,𝐉^2]=0,[H,J_i]=0,[𝐉^2,J_i]=0`$ (4)
$`[J_i,J_j]=iϵ_{ijk}J_k,`$ (5)
Choosing the vector potential as

$$𝐀=q\frac{1+\mathrm{cos}\theta }{r\mathrm{sin}\theta }\widehat{𝐞}_\phi .$$

we find

$`J_\pm =e^{\pm i\phi }\left(\pm {\displaystyle \frac{}{\theta }}+i\mathrm{cot}\theta {\displaystyle \frac{}{\phi }}\mu {\displaystyle \frac{1+\mathrm{cos}\theta }{\mathrm{sin}\theta }}\right),`$ (6)
$`J_0=i{\displaystyle \frac{}{\phi }}+\mu ,`$ (7)
$`𝐉^2={\displaystyle \frac{1}{\mathrm{sin}\theta }}{\displaystyle \frac{}{\theta }}\left(\mathrm{sin}\theta {\displaystyle \frac{}{\theta }}\right){\displaystyle \frac{1}{\mathrm{sin}^2\theta }}{\displaystyle \frac{^2}{\phi ^2}}`$
$`i{\displaystyle \frac{2\mu }{1\mathrm{cos}\theta }}{\displaystyle \frac{}{\phi }}+\mu ^2{\displaystyle \frac{1+\mathrm{cos}\theta }{1\mathrm{cos}\theta }}+\mu ^2`$ (8)
where $`J_\pm =J_x\pm iJ_y`$ are the raising and the lowering operators for $`J_0=J_z`$ satisfying the standard commutations relations

$$[J_0,J_\pm ]=\pm J_\pm ,[J_+,J_{}]=2J_0.$$

It is known that for the Dirac’s monopole problem the Schrödinger equation written in the spherical coordinates admits the separation of variables and the eigenfunctions of the angular part of the Hamiltonian are Wu2

$$Y_{jn\mu }(\theta ,\phi )=\frac{1}{\sqrt{2\pi }}\text{e}^{i(j\mu n)\phi }P_{jn\mu }(\mathrm{cos}\theta ),$$
(9)
with $`n`$ as a non-negative integer and

$`P_{jn\mu }(u)={\displaystyle \frac{1}{\sqrt{h_{jn\mu }}}}(1u)^{\alpha /2}(1+u)^{\beta /2}P_n^{(\alpha ,\beta )}(u),`$
$`h_{jn\mu }={\displaystyle \frac{2^{2(jn)+1}}{2j+1}}{\displaystyle \frac{\mathrm{\Gamma }(j+\mu +1)\mathrm{\Gamma }(j\mu +1)}{n!\mathrm{\Gamma }(2jn+1)}},`$
$`\alpha =j+\mu n,\beta =j\mu n,`$
here $`\mathrm{\Gamma }`$ is the Gamma function, $`j(j+1)`$ and $`jn`$ are the eigenvalues of $`𝐉^2`$ and $`J_z`$ respectively, and $`P_n^{(\alpha ,\beta )}`$ are the Jacobi polynomials.

Using the rotational invariance of the system and restricting by the unitary representations of the group SO(3) one can obtain the Dirac’s quantization rule $`2\mu =`$ integer. This results in the following: the only way to avoid Dirac’s rule is to consider a nonunitary representation of the rotation group.

Further we will try to describe the charge-monopole system as a free particle with an arbitrary spin, the so-called anyons, and relate this with the nonunitary representations of the rotation group. The idea to describe Dirac’s monopole as a particle with a spin is not new Gol ; Gol2 ; Ply . However, in the known models anyon’s statistics is based on the SO(2,1)-symmetry Ply ; Jan . and it is not clear how anyons may be appeared in a non-relativistic theory and its relation with the group SO(3).

Following Ply  we consider the Lagrangian describing a non-relativistic charged particle in the field of a magnetic monopole

$$L=\frac{m}{2}\dot{𝐫}^2+e𝐀\dot{𝐫}$$

where $`𝐀`$ is the vector potential determining the magnetic field (2). In the limit $`m0`$, the action is given by the integral of the one-form $`\theta =eA_idx^i`$ and we have

$$d\theta =\frac{1}{2𝐒^2}\epsilon _{ijk}S_idS_jdS_k,𝐒=\mu \frac{𝐫}{r}$$
(10)
The two-form $`d\theta `$ is closed and nondegenerate, and is interpreted as a symplectic form corresponding to the symplectic potential $`\theta `$ defined on the unit two-sphere. It can be used to introduce Poisson brackets on $`S^2`$ as follows

$$\{S_i,S_j\}=\epsilon _{ijk}S_k,$$
(11)
and thus the angular momentum of the electromagnetic field $`𝐋_{em}`$ can be considered as a classical spin.

Developing this idea, one can give the alternative description of the charge-monopole system as a free particle of the fixed spin with translational and spin degrees of freedom. The corresponding constrained Hamiltonian is given by

$`H={\displaystyle \frac{1}{2m}}\left(𝐩{\displaystyle \frac{1}{r^2}}(𝐫\times 𝐒)\right)^2+\lambda (𝐒𝐫+\mu r)`$ (12)
where the second term describes the spin-orbit interaction. On the constraint surface $`𝐒𝐫+\mu r0`$ the dynamics generated by the Hamiltonian (12) is exactly the same as the dynamics of the initial charge-monopole system Ply .

The above considerations is based on the conventional approach to the spin and therefore involves implicitly the Dirac’s quantization rule. The case of $`\mu `$ being not necessarily integer or half integer implies making use of the nonunitary representations of the group SO(3).

Let the set of $`\{S_\pm ,S_0\}`$ and $`\{L_\pm ,L_0\}`$ form the algebra $`so(2,1)`$ and $`so(3)`$ respectively

$`[S_+,S_{}]=2S_0,[S_\pm ,S_0]=S_\pm `$ (13)
$`[L_+,L_{}]=2L_0,[L_\pm ,L_0]=L_\pm .`$ (14)
We introduce a direct sum $`so(3)so(2,1)`$ as follows

$$J_0=L_0+S_0,J_+=L_+S_+,J_{}=L_{}+S_{}.$$
(15)
The computation shows that the operator $`𝐉`$ obeys the commutation relations of the algebra $`so(3)`$

$$[J_0,J_\pm ]=\pm J_\pm ,[J_+,J_{}]=2J_0$$

and identifying $`J_i`$’s as

$`J_1=L_1+iS_2,J_2=L_2+iS_1,J_3=J_0`$ (16)
we find

$$[J_i,J_j]=i\epsilon _{ijk}J_k.$$

Notice that the group $`\mathrm{SO}(3)\mathrm{SO}(2,1)`$ has previously appeared in Dirac monopole theory, as a dynamical invariance group Jac1 .

We construct the representation by using the generators of the group $`SO(3)`$ with the standard action on the state $`|l,\stackrel{~}{m}`$ given by

$`𝐋^2|l,\stackrel{~}{m}=l(l+1)|l,\stackrel{~}{m},`$
$`L_0|l,\stackrel{~}{m}=\stackrel{~}{m}|l,\stackrel{~}{m},`$
$`L_+|l,\stackrel{~}{m}=\sqrt{(l\stackrel{~}{m})(l+\stackrel{~}{m}+1)}|l,\stackrel{~}{m}+1,`$
$`L_{}|l,\stackrel{~}{m}=\sqrt{(l+\stackrel{~}{m})(l\stackrel{~}{m}+1)}|l,\stackrel{~}{m}1,`$
and the unitary infinite dimensional representation of the group $`SU(1,1)`$ relating to the generators $`S_i`$’s.

For the representation bounded above we have Jan

$`𝐒^2|\lambda ,n=\lambda (\lambda 1)|\lambda ,n,`$
$`S_0|\lambda ,n=(\lambda +n)|\lambda ,n,`$
$`S_+|\lambda ,n=\sqrt{(2\lambda +n1)n}|\lambda ,n1,`$
$`S_{}|\lambda ,n=\sqrt{(2\lambda +n)(n+1)}|\lambda ,n+1,`$
where $`n=0,1,2,\mathrm{},\mathrm{}`$ and $`\lambda `$ is an arbitrary parameter, and for the highest-weight state $`|\lambda ,0`$ one obtains

$`S_0|\lambda ,0=\lambda |\lambda ,0,S_+|\lambda ,0=0.`$
Thus $`\lambda `$ is the eigenvalue of $`J_0`$ for the state $`|\lambda ,0`$ and characterizes the representation. Recently the infinite dimensional representations of the group SO(2,1) have been used for the description of a theory of anyons where the spin of the one-particle states is taken to be $`s=1\lambda `$ Jan  (see also Ply1 ; Ply2 ).

Starting with the highest-weight state $`|j,j=|l,l|\lambda ,0`$ corresponding to the eigenvalue $`l\lambda `$ of the operator $`J_0`$ and applying the lowering operator $`J_{}`$ one can construct all states as follows

$`J_{}|j,m=\sqrt{(j+m)(jm+1)}|j,m1,`$ (17)
$`J_+|j,m=\sqrt{(jm)(j+m+1)}|j,m1,`$ (18)
$`m=j,j1,\mathrm{},\mathrm{}`$
It is easy to see that the arised representation is nonunitary infinite dimensional and bounded above by $`mj`$, instead of $`jmj`$ well known for the unitary finite-dimensional representations of the rotation group.

The analysis of the nonunitary representation shows that the allowed values of $`j,\mu `$ are $`\mu ^2j(j+1)`$ and $`j\mu `$ = integer, which is an alternative choice to the Dirac rule Nes .

In our model the charge-monopole system is interpreted as a free anyon with translational and spin degrees of freedom. The close approach has been applied in Ply3 ; Ply4  for description of a fractional spin in $`(2+1)`$ dimensions.

Recently Martínez-y-Romero et al Mar ; Mar1  have used a nonunitary representation for a Dirac particle in a Coulomb-like field and Davis and Ghandour Dav  have showed that nonunitary transformations are not so strange in quantum mechanics. The physical consequences of using nonunitary representations in Dirac’s monopole problem are not clear till now and the work is in progress.

###### Acknowledgements.

We would like to thank R. Jackiw and M. Plyushchay for helpful comments on Dirac monopole theory and related subjects.


text:
# A metallicity-flattening relation for dwarf elliptical galaxies

## 1 Introduction

Elliptical galaxies follow a well-known luminosity-metallicity relation, in the sense that more luminous ellipticals are observed to be globally redder, or have larger metallic absorption line strengths in their spectra than less luminous systems (Faber 1973, Visvanathan & Sandage 1977). Obviously this is a mass-metallicity relation, and its canonical interpretation is that less massive galaxies had more significant outflows of metal-enriched gas at an early evolutionary stage (Faber 1973, Mould 1984). The relation seems to hold for the entire luminosity range of spheroidal galaxies, down to the faintest dwarfs (Caldwell 1983, Brodie & Huchra 1991, Caldwell et al. 1992). However, there is also considerable scatter in the relation, and there have been attempts to find a second parameter which governs metallicity (mass being the primary parameter). A good candidate is the internal velocity dispersion (Terlevich et al. 1984, Efstathiou & Fall 1984), which even more strongly correlates with metallicity than mass (Bender et al. 1993).

The largest scatter in the luminosity-metallicity relation is observed in the magnitude range of bright cluster dwarf ellipticals ($`14M_B18`$). While this scatter could plausibly be due to measurement errors, Rakos et al. (2001) find a weak correlation with age for a sample of Fornax cluster dwarfs observed in Stroemgren narrow-band colours.

In this letter we show that the scatter in metallicity for cluster dEs is largely explained by apparent ellipticity (flattening). At a given luminosity, rounder dwarf ellipticals are more metal-rich. There have already been hints that ellipticity might act as a second parameter in normal elliptical galaxies (Terlevich et al. 1984). However, for dEs the effect is so strong that ellipticity appears to be the primary parameter. Such a metallicity-flattening relation for dEs is not implausible. The outflow of metal-enriched gas in stellar systems of intermediate mass ($`M10^9`$ $`M_{}`$) is preferentially occurring along the minor axis, rendering ellipticity a critical parameter for the metallicity of present day dwarfs (De Young & Heckman 1994). Rounder dEs seem to have suffered less significant outflow.

## 2 Colour versus flattening

We first show the evidence from multi-colour photometry, which is how we stumbled over the flattening effect. In Fig. 1 we have plotted integrated colour versus apparent ellipticity for different samples of spheroidal galaxies. The filled circles are $`UB`$ values for 15 Virgo dEs from our own photometry (Barazza et al., in prep.). There is a very clear trend appearing, in the sense that rounder dEs are redder than more flattened ones. Such a correlation is rather unexpected – see for comparison a sample of normal ellipticals from Peletier et al. (1990, triangles). However, our finding is corroborated by the $`UB`$ data for Virgo and Fornax dEs (crossed circles) from Caldwell (1983) and Caldwell & Bothun (1987), and the Stroemgren $`vzyz`$ data for Fornax dEs (open circles) from Rakos et al. (2001). In contrast, local dwarfs with colour data available (from Mateo 1998, asterisks) show no clear relation. All ellipticities used here are taken either from the references given for colour/metallicity data, or from LEDA.

Is the trend for cluster dEs real? Observational bias can be excluded: the samples were not selected by colour or ellipticity. Nor is it conceivable that the photometric errors in colour somehow depend on the ellipticity of the whole galaxy. Another suspicion is that ellipticity is not an independent parameter: if ellipticity were correlated with absolute magnitude, we simply might have recovered the well-known luminosity-colour relation (e.g. Ferguson 1994). However, the ellipticity of dEs is not significantly related to luminosity, nor to surface brightness (Binggeli & Popescu 1995).

We take it for granted that colour is essentially indicating metallicity, as spheroidal galaxies are believed to be dust-free and old. Of course, metallicity can also be measured more directly. We now show that the colour-ellipticity relation is indeed a metallicity-ellipticity relation.

## 3 The luminosity-metallicity relation

Fig. 2 shows the “universal” luminosity-metallicity relation for essentially the same galaxies as in Fig. 1. Normal Es (small triangles) are represented by data from a more recent study (Kobayashi & Arimoto 1999). \[Fe/H\] values for local dwarfs (asterisks) are again from Mateo (1998). For bright cluster dEs ($`14M_B18`$) we have three sources for \[Fe/H\]: 1) Brodie & Huchra (1991) give spectroscopically determined metallicities for nearly 40 Virgo and Fornax cluster dEs (essentially the Caldwell et al. sample) based on several different line indices. However, the estimated errors in \[Fe/H\] are considerable. Data with 1 $`\sigma `$ errors less than 0.6 dex are plotted as filled circles, those with larger errors (but smaller than 1 dex) as filled squares. 2) Held & Mould (1994) provide spectroscopically measured \[Fe/H\] values, with much smaller estimated errors, for 8 nucleated Fornax dEs (crossed circles). 3) Finally, there are the Rakos et al. (2001) $`vzyz`$ Stroemgren colours for Fornax dEs, appropriately transformed into \[Fe/H\] values (open circles), with a small mean error of 0.2 dex.

Naively fitting a line through all data plotted in Fig. 2, we get an expression for the mean (universal) luminosity-metallicity relation: \[Fe/H\] = $`0.1747M_B3.6405`$, which will be used for a residual analysis below. The local dwarfs follow the mean relation surprisingly well. In contrast, the cluster dEs in the intermediate range between normal Es and extreme dwarfs, while falling in place with the universal relation in the mean, show enormous scatter in their individual metallicities. The natural suspicion is that this is simply due to large errors in \[Fe/H\]. However, we now show that this scatter is systematically related to ellipticity, i.e. it must, at least partially, be real.

## 4 Flattening as second parameter

The residual metallicities with respect to the line shown in Fig. 2 are plotted versus ellipticity in Fig. 3. There is clearly no correlation for normal Es and local dwarfs (upper panel), whereas a relation is evident for cluster dEs (lower panel): dEs with positive residuals are throughout round, while highly flattened dEs have always negative residuals. At a given luminosity, rounder dEs tend to be more metal-rich.

Several points have to be noted. First, the trend is followed by all dE samples used, even those where large errors were claimed (the Brodie & Huchra data). Again it is not conceivable why these errors, were they real, should correlate with ellipticity (we believe the errors are simply overestimated). Second, we have tested that the metallicity residuals do not correlate with other parameters, such as effective radius, i.e. there is no other “second” parameter than ellipticity. Third, the true physical relation behind this effect will involve intrinsic ellipticity. Due to random projection this relation should even be stronger than what we see with apparent ellipticity. The distribution of points in Fig. 3 (lower panel) is in accord with what we expect from projection effects: apparently round galaxies with positive residuals will also be intrinsically round, while some of the apparently round galaxies with negative residuals will intrinsically be more flattened (i.e. would be shifted to the right).

## 5 The metallicity-flattening relation

In Fig. 4 we have plotted metallicity directly versus ellipticity, to get an equivalent to the colour-ellipticity diagram of Fig. 1. As expected, no correlation is seen for giant Es, while dwarf objects evidently follow the suspected relation: rounder galaxies tend to be more metal-rich than flatter ones. Owing to some outliers, the trend is less evident for the local dwarfs (no correlation at all was seen in the residual plot). Conceivably, local (field) dwarfs are systematically different from cluster dwarfs (see Sect. 6).

In general, the flattening-metallicity relation is not as striking as the residual plots in Fig. 3. This is because luminosity now acts as a hidden parameter, adding scatter. In Fig. 5 we have plotted metallicity versus logarithmic ellipticity for only cluster dwarfs, which provides an almost linear relation. A linear fit yields \[Fe/H\]=$`0.72\mathrm{log}(ϵ)1.56`$, or in terms of abundance: $`Z/Z_{}ϵ^{0.72}`$.

Including the absolute magnitude and fitting a plane to the cluster dwarfs in the parameter space defined by \[Fe/H\], $`ϵ`$ and $`M_B`$, we get $`Z/Z_{}ϵ^{0.70}L^{0.20}`$, which clearly shows that ellipticity is in fact the primary parameter, having a stronger effect on the metallicity of dEs than the luminosity.

## 6 Discussion of a possible explanation: outflow of metal-rich gas along the minor axis

The mass-metallicity relation of spheroidal galaxies, evident in Fig. 2, is most likely the result of early chemodynamical evolution. The central idea is gas loss by a galactic wind, which is believed to be particularly metal-enhanced (Vader 1986, 1987). Certainly, the mass (i.e. potential depth) of a galaxy must be a key parameter for the loss of gas (hence metals), and this explains the mass-metallicity relation. However, it has been pointed out by De Young & Heckman (1994) that the shape (i.e. flattening) of a galaxy may play a crucial role as well. These authors show that galaxies of intermediate mass ($`M10^9M_{}`$) would not lose almost all of their ISM in a central starburst, but would preferentially undergo a blowout event, having gas flowing out in the direction of their minor axis. The important point is that the strength of such an outflow would critically depend on the intrinsic flattening of the ISM distribution, i.e. the galaxy. Our cluster dwarfs showing a metallicity-flattening relation have masses precisely in the critical regime discussed by De Young & Heckman (1994): their mass seems to be too low to be completely protected from gas loss (as in giant ellipticals), but also too high to experience a total blowaway (as in dwarf spheroidals). In the intermediate mass range of bright cluster dEs the ellipticity could conceivably be the dominant parameter for gas loss, and hence metallicity.

Of course, reality must be more complicated. For one thing dwarf galaxies are believed to have large amounts of dark matter, so the gas escaping the stellar body should remain bound to the galaxy and be re-accreted on a rather short time scale (Mac Low & Ferrara 1999, Ferrara & Tolstoy 2000). In this context, it is interesting to note that we have found a possible difference between cluster dwarfs and (local) field dwarfs. While the evidence for dark matter in field dwarfs abounds, the situation with cluster galaxies in general is much less clear. X-ray and lensing studies of clusters of galaxies suggest that the dark matter in clusters is not bound to individual galaxies, hence cluster dwarfs might be much less dark than field dwarfs. An additional difference is that cluster dwarfs are subject to ram-pressure stripping by the ICM, which might help to remove the gas flown out (e.g. Murakami & Babul (1999).

## 7 Conclusion

We have presented evidence for a metallicity-flattening relation for dwarf elliptical galaxies based on colour and metallicity data available from the literature. At a given total magnitude, rounder dEs are more metal-rich. In the narrow magnitude range of bright cluster dwarfs ($`14M_B18`$), metallicity is more strongly correlated with ellipticity than luminosity (mass), i.e. ellipticity seems to be the primary parameter for the enrichment history of these galaxies. Possibly this holds only for cluster dwarfs; the effect is not significant for local dwarfs. A possible explanation is provided by the scenario of De Young & Heckman (1994), where the outflow of gas, and hence the regulation of metallicity, depends on the intrinsic shape of a galaxy in the intermediate mass range around $`10^9M_{}`$.

It would be highly desirable to strengthen the evidence with further observations, especially spectroscopically well determined metallicities for many cluster dwarfs. If confirmed, the effect is likely to be of importance for our understanding of the chemodynamical evolution of galaxies.

###### Acknowledgements.

We are grateful to the Swiss National Science Foundation for financial support. We have made use of the LEDA database (http://leda.univ-lyon1.fr).


text:
# Volatility in atmospheric temperature variability

## 1. Introduction

Indications of weather persistence over months and seasons are known . This long term persistence in the atmospheric temperature variability is analysed and quantified by detrended fluctuation analysis (DFA) and wavelet transform techniques with a fluctuation exponent of $`\alpha 0.65`$, which is independent of the location of the site . Also, indications of persistence are confirmed through power spectral analysis . Presence of such a universal persistence law indicates that the processes governing the atmospheric dynamics at different climatological zones are based on similar principles . Here, we are not interested in the temperature fluctuations around the seasonal trend, but in the magnitude of temperature changes between successive days, i.e the temperature volatility. Volatility is the concept often used in econophysics to indicate the fluctuations in price changes . The market is said to be more volatile if the fluctuations in price changes are high . In recent years this concept has been applied to cardiac system and found that cardiac volatility series is long range correlated .

Here we study the volatility of atmospheric temperature data obtained from 10 randomly chosen meteorological stations in Europe, North America, Asia, Russia and Australia, from various climatological zones. Correlations in the volatility series give information about persistence in the changes. If, for example, the change in the temperature between two successive days is small there is a high tendency that the change remains similar for the next consecutive days. Here, we study long-term temperature records (typically 100 years). We use detrended fluctuation analysis (DFA) up to the order of five , which systematically eliminates higher order trends and reveals the correlations present in the highly non-stationary data. Our analysis shows that (i) the persistence, characterised by the correlation $`C(s)`$ of the volatility series separated by $`s`$ days, follows a power law, $`C(s)s^\gamma `$, with roughly the same exponent $`\gamma 0.8`$ for all stations considered, and that (ii) the range of this universal persistence law seems exceed one decade.

## 2. Methodology of Scaling analysis

We have performed the scaling analysis of volatility series for the records of the maximum daily temperature $`T_i`$ of the following weather stations: Prague (218 yr), Melbourne (136 yr), Luling (90 yr), Seoul (86 yr), Kasan (96 yr), Vancouver (93 yr), Tashkent (97 yr), New York city (116 yr), Brookings (99 yr) and St. Petersburg (111 yr). The numbers within (.) are the length of the records. From the daily maximum temperature series, we construct the volatility series $`V_i`$=$`|T_{i+1}T_i|`$. Then we remove the climatological annual cycle  from $`V_i`$ to obtain $`\mathrm{\Delta }V_i=V_i<V_i>`$.

Qualitatively, persistence is clearly seen (patches of high and low volatility) in the plot of $`V_i`$ as shown in Fig. 1a for one year in Prague. Figure 1b shows the volatility series $`\mathrm{\Delta }V_i`$ around the mean. Figures 1e and 1f show the plot of $`V_i`$ and $`\mathrm{\Delta }V_i`$ respectively, obtained from the phase randomised surrogate data (preserving the distribution)  of the temperature increment series $`(T_{i+1}T_i)`$. Quantitatively, persistence in $`\mathrm{\Delta }V_i`$ can be characterised by the (auto)correlation function,

$$C(s)\mathrm{\Delta }V_i\mathrm{\Delta }V_{i+s}=\frac{1}{Ns}\underset{i=1}{\overset{Ns}{}}\mathrm{\Delta }V_i\mathrm{\Delta }V_{i+s},$$
$`(1)`$
where $`N`$ is the number of days in the record. A direct calculation of $`C(s)`$ is hindered by the level of noise present in finite temperature series, and by possible nonstationarities in the data. Following , we do not calculate $`C(s)`$ directly, rather we study the fluctuations in the volatility profile $`Y_n=_{i=1}^n\mathrm{\Delta }V_i`$. We divide the profile into non overlapping time windows of length $`s`$ and determine the squared fluctuations of the profile (as specified below) in each segment. The mean square fluctuations, averaged over all segments of length $`s`$, are related to the correlation function $`C(s)`$ (see below). For this study we employ a hierarchy of methods that differ in the way the fluctuations are measured and possible nonstationarities are eliminated (see e.g., for detailed description of the methods):

(i) In the simple fluctuation analysis (FA), we calculate the difference of the profile at both ends of each segment. The square of this difference represents the square of the fluctuations in each segment.

(ii) In the first order detrended fluctuation analysis, we determine in each segment the best linear fit of the profile. The variance of the profile from this straight line represents the square of the fluctuations in each segment.

(iii) In general, in the $`n`$-th order DFA we determine in each segment the best $`n`$-th order polynomial fit of the profile. The variance of the profile from these best $`n`$-th order polynomials represents the square of the fluctuations in each segment.

The fluctuation function $`F(s)`$ is the root mean square of the fluctuations in all segments. For the relevant case of long-term power-law correlations, $`C(s)s^\gamma `$, with $`0<\gamma <1`$, the fluctuation function increases with $`s`$ according to a power law ,

$$F(s)s^\alpha ,\alpha =1\frac{\gamma }{2}.$$
$`(2)`$
For uncorrelated as well as short range correlated data, we have $`\alpha =\frac{1}{2}`$. For long range correlated data we have $`\alpha >\frac{1}{2}`$.

By definition, FA does not eliminate trends similar to the Hurst method and the conventional power spectral methods . In contrast, DFAn eliminates trends of order $`n`$ in the profile and $`n1`$ in the original time series. Thus, from the comparison of fluctuation functions $`F(s)`$ obtained from different methods one can learn about long term correlations and types of trends, which cannot be achieved by the conventional techniques.

## 3. Scaling analysis of Temperature volatility series

We begin the analysis with the volatility series $`\mathrm{\Delta }V_i`$ for Prague which is the longest series (218 yr) in this study. Figure 1c shows the profile, and Fig. 1d shows the fluctuation functions obtained from FA ($`\mathrm{}`$) and DFA1-5 (o). In the log-log plot, all curves are approximately straight lines with a slope of $`\alpha `$ close to 0.6. This result suggests that there exists a long term persistence expressed by the power law decay of correlation exponent $`\gamma 0.8`$. There is slight upward bend in the $`F(s)`$ curve obtained from FA. This shows that there is a trend in the volatility series, as observed in  for the temperature data. Hence, the effect of the city growth in Prague does not lead only to an increase in the mean temperature of the city, but also to an increase in the temperature volatility. This trend is not removed in FA while DFA1 and subsequently higher orders of DFA have removed it. Figure 1g shows the profile of the volatility series obtained from the surrogate data of the temperature increment series. When compared with Fig. 1c, profile shows much less fluctuations. Figure 1h shows the fluctuation functions $`F(s)`$ obtained from FA ($`\mathrm{}`$) and DFA1-5 (o) for the profile shown in Fig. 1g. In the log-log plot, the fluctuation functions $`F(s)`$ are straight lines with a slope of $`\alpha \frac{1}{2}`$. This shows that the original data have nonlinear features .

Figures 2(a-d) show the results of volatility analysis for four different continental sites: Tashkent from Asia (97 yr), New York city (116 yr), Melbourne (136 yr) and Luling from Texas (90 yr). The top two curves in each panel are the fluctuation curves obtained from FA ($`\mathrm{}`$) and DFA 2 ($``$) for volatility series of temperature data while the third and fourth curves are those obtained from the same analyses (FA and DFA2) of volatility series derived from the surrogate data. For the sake of clarity, we have divided the fluctuation function $`F(s)`$ by $`s^{0.5}`$. The straight line shown in each panel has a slope 0.1. The scaled fluctuation functions $`F(s)/s^{0.5}`$ obtained from the original series $`\mathrm{\Delta }V_i`$ (the first and second curves in all panels shown in Fig. 2) are almost parallel to the line with slope 0.1. The fluctuation functions $`F(s)/s^{0.5}`$ obtained from the surrogate data (third and fourth curves in all panels shown in Fig. 2) are all parallel to x-axis.

We obtained similar results with exponents between 0.58 and 0.63 for all 10 sites analysed. Thus, there exist a quite general persistence law in the volatility series obtained from temperature data. The existence of such a general persistence law in the temperature volatility series indicates that the change in the temperature in different climatic zones may be governed by the same basic principles, leading to similar fluctuations in volatility series in different places.

## 4 Testing the volatility scaling performance of simulated temperature records

We use this result to test the state-of-the art global climate models. In an earlier work, Govindan et al., have shown that temperature data simulated by GCMs violate the observed scaling behavior . Here we consider the scaling analysis of the temperature volatility time series. We concentrate on four models, for which data for all three scenarios viz. Control run (CR), (ii) greenhouse gas forcing only (GHGF) and (iii) greenhouse gas plus aerosol forcing (GHGPS), are available for the same simulation period. The GCMs are: (i) CSIRO-Mk2 (Melbourne), (ii) ECHAM4/OPYC3 (Hamburg), (iii) CGCM1(Victoria, Canada) and (iv) CCSR/NIES (Tokyo) (see  for details).

In CR, the CO2 content is fixed. In GHGF scenario, one mainly considers the effect of greenhouse gases. The amount of greenhouse gas forcing is taken from the historic data until 1990 and then increased at a rate of 1% per year. In GHGPS scenario, the effect of aerosols (mainly sulphates) in the atmosphere is taken into account which can mitigate and partially offset the greenhouse warming. Although this scenario represents an important step towards comprehensive climate simulation, the precise role of aerosols in the mechanism of climate modelling is still unclear.

The temperature data simulated by these four models for three different scenarios are available from the IPCC Data Distribution Center . We extracted data for nine representative sites around the globe (Prague, Melbourne, Seoul, Vancouver, Kasan, Luling, Tashkent, New York and St. Petersburg). For each model and each scenario, we selected the temperature records of the four grid points closest to each site, and bilinearly interpolated the data to the location of the site.

We generate volatility series from the monthly temperature series by the same procedure (explained earlier) for daily temperature data. We removed the seasonal periodic trends in $`T_i`$ (before constructing the increment series) and also in the volatility series $`V_i`$. Since the GCMs simulate temperature data in monthly time scales, for comparison, we averaged the daily temperature (observed) data (of the nine climatological stations) to monthly data. We present our results of scaling analysis of volatility series obtained from temperature data simulated by models, for three different scenarios and also for the observed data in the form of histograms e.g. Fig. 3.

The scaling exponents $`\alpha `$ obtained for the observed data are shown in the topmost panel (for comparison they are plotted repeatedly for three times corresponding to three different scenarios). We can see immediately that all sites occupy the bin corresponding to the $`\alpha `$ value in the range of 0.58 to 0.62 (see Fig. 3 top panels) indicating the universal behavior as found in the daily data. The grey color boxes represent the $`\alpha `$ values of the volatility series derived from surrogate series. For the case of monthly data the $`\alpha `$ values obtained from the volatility series of surrogate data exhibit correlated behavior for some of the sites. However, they are well below the values of the original data. When we compare the results obtained from GCMs, none of the models show a unique behavior as found in the observed data. The $`\alpha `$ values are distributed widely for all the models. The comparison of the exponents of volatility series obtained from original data and that of surrogate data shows their difference is close to zero and even zero for some of the sites. For instance, in CSIRO model, for each of the three different scenarios the exponent values of the volatility series of the original data and that derived from surrogate series, have same value for St. Petersburg(9). Likewise, in ECHAM4, similar conclusion can be drawn for Kasan(5) as well as Luling(6). However, the distribution of $`\alpha `$ values of volatility series obtained from surrogate data, towards higher values is clearly seen for the well established scenarios like GHGF and GHGPS.

It follows from our analysis that there is a universal persistence in the volatility series obtained from temperature increment series with an exponent of $`0.60\pm 0.03`$. When we use this result to test the scaling performance of the virtual climate records simulated by GCMs, we find that (i) models data display wide range of exponent values, (ii) surrogate analysis suggests that models data lack nonlinearity, especially for the well established scenarios, for some of the sites considered here.

This work has been supported by the Deutsche Forschungsgemeinschaft and the Israel Science Foundation. We wish to thank Y. Ashkenazy for useful discussions.


text:
# On the Oosterhoff Classification of the Unusual, Metal-Rich Globular Clusters NGC 6388 and NGC 6441

## 1. Introduction

NGC 6388 and NGC 6441 are unusual in several respects, including the following. (i) In spite of being metal-rich ($`[\mathrm{Fe}/\mathrm{H}]0.55`$ dex), they contain prominent extensions to their red HBs—including both RRL and extended blue tails (Rich et al. 1997). (ii) As pointed out by Sweigart & Catelan (1998), the normally “horizontal” part of these clusters’ HBs is actually strongly tilted, with a $`\mathrm{\Delta }V0.5`$ mag from the lower part of the red HB to the tip of the blue HB—which cannot be accounted for by canonical stellar evolution models. (iii) The mean periods of the RRL in both clusters are extremely long for their metallicity—longer, in fact, even than in metal-poor, OoII globular clusters, thus apparently breaking down the traditional Oosterhoff classification scheme and posing yet another serious challenge to the models (Pritzl et al. 2000, 2001, 2002).

## 2. On the Oosterhoff Classification of NGC 6388 and NGC 6441

Such peculiarities notwithstanding, it has recently been suggested (Clement et al. 2001), on the basis of the RRL’s position on the $`\mathrm{log}PA_V`$ plane—i.e., reportedly coincident with the “line” occupied by “normal” RRL stars in OoII globular clusters (Clement & Rowe 2000)—that these globulars should be classified as OoII, with the RRL evolved away from a position on the blue ZAHB. This suggestion raises some fundamental questions, including the following:

(i) Do blue HB stars spend enough time within the instability strip as they evolve redward to the asymptotic giant branch (AGB) to produce the observed number of RRL in OoII globular clusters?

(ii) Does the predicted period-$`T_{\mathrm{eff}}`$ relation for a given OoII cluster depend on the stellar distribution along the blue HB?

(iii) Is the period-$`T_{\mathrm{eff}}`$ relation independent of metallicity?

To answer these questions, we have computed new evolutionary tracks and HB simulations for OoII globular clusters. Our main results (see Pritzl et al. 2002 for a more thorough discussion) can be summarized as follows:

(i) According to the models, blue HB stars do not spend enough time within the instability strip as they evolve redward to the AGB to produce the observed number of RRL in OoII globulars. This confirms the well-known arguments of Renzini & Fusi Pecci (1988) and Rood & Crocker (1989), which unfortunately have been almost entirely neglected in the literature since the mid-90’s.

(ii) The predicted period-$`T_{\mathrm{eff}}`$ relation for a given cluster depends somewhat on the stellar distribution along the blue HB, leading to the presence of intrinsic scatter in the period-$`T_{\mathrm{eff}}`$ plane and to longer mean period shifts for OoII clusters with bluer HBs. This has to be taken into account, when attempting to classify individual stars into an Oosterhoff type. On the other hand, only stars within a relatively narrow range of ZAHB colors, $`(BV)_0>0`$, have a significant chance of becoming redward-evolving RRL.

(iii) The predicted period-$`T_{\mathrm{eff}}`$ relation for redward-evolving RRL is a weak function of metallicity for OoII metallicities, but not for the metallicities of NGC 6388 and NGC 6441.

### Acknowledgments.

This work has been supported by the National Science Foundation under grants AST 9528080 and AST 9986943. Support for M.C. was provided by Proyecto de Inicio DIPUC 2002-04E. A.V.S. acknowledges support from NASA Astrophysics Theory Program proposal NRA-99-01-ATP-039.

## References

Clement, C., et al. 2001, AJ, 122, 2587

Clement, C., & Rowe, J. 2000, AJ, 120, 2579

Pritzl, B., Smith, H. A., Catelan, M., & Sweigart, A. V. 2000, ApJ, 530, L41

Pritzl, B., Smith, H. A., Catelan, M., & Sweigart, A. V. 2001, AJ, 122, 2600

Pritzl, B., Smith, H. A., Catelan, M., & Sweigart, A. V. 2002, AJ, 124, 949

Renzini, A., & Fusi Pecci, F. 1988, ARA&A, 26, 199

Rich, R. M., et al. 1997, ApJ, 484, L25

Rood, R. T., & Crocker, D. A. 1989, in IAU Colloq. 111, The Use of Pulsating Stars in Fundamental Problems of Astronomy, ed. E. G. Schmidt (Cambridge: Cambridge University Press), 103

Sweigart, A. V., & Catelan, M. 1998, ApJ, 501, L63


text:
# The glueball among the light scalar mesons

## 1 Introduction

The existence of glueballs is among the fundamental predictions of QCD, but the experimental evidence is still in debate. In QCD calculations on the lattice in quenched approximation the lightest glueball is found to have quantum numbers $`J^{PC}=0^{++}`$ and a mass in the region 1400-1800 MeV (for review, see ). The effect of this approximation is still being investigated. In an alternative approach based on QCD sum rules a gluonic state of lower mass around 1000 MeV is required as well .

The search for the lightest glueball should therefore concentrate on the mass region up to about 1800 MeV in the scalar sector. This search has to proceed in parallel with the identification of the low mass scalar $`q\overline{q}`$ nonet(s). The glueball candidate should fulfill some general properties, it should be produced in particular in a gluon rich environment and its decay (for the unmixed glueball) should be “flavour-blind”. Then the scalar states to be identified as members of the $`q\overline{q}`$ nonet or glueball should be found from the list provided by the Particle Data Group
I=0: $`f_0`$(400-1200), ($`\sigma (600)`$?), $`f_0(980)`$, $`f_0(1370)`$, $`f_0(1500)`$, $`f_0(1710)`$
I=$`\frac{1}{2}`$: ($`\kappa (900)`$?), $`K_0^{}(1430)`$, $`K^{}(1950)`$
I=1: $`a_0(980)`$, $`a_0(1450)`$

There are different scenarios for interpretation which include Scenario A:
A starting point is the lattice result, then a glueball with suitable mass is $`f_0(1500)`$ . As this state does not obey all wanted properties one includes the nearby $`f_0(1370)`$ and $`f_0(1710)`$ and all three scalars mix with the glueball and two $`q\overline{q}`$ states (for review of this popular approach, see ). In this case the low mass states like $`a_0(980)`$ and $`f_0(980)`$, possibly also $`\kappa `$ and $`\sigma `$, are considered as multiquark states which should not be included in the $`q\overline{q}`$ spectroscopy.
Scenario B:
One tries to identify the $`q\overline{q}`$ nonet first including the $`a_0`$ and $`f_0`$. In our approach  the $`f_0(980)`$ and $`f_0(1500)`$ are the isoscalars of the $`q\overline{q}`$ nonet with large mixing, very similar to the $`\eta ,\eta ^{}`$ in the pseudoscalar sector. The nonet is completed with $`a_0(980)`$ and $`K_0^{}(1430)`$. The $`\sigma `$ and $`\kappa `$ are not considered as genuine resonances. The remaining states called $`f_0(4001200)`$ and $`f_0(1370)`$ (“red dragon” ) correspond to a single broad state which we suggest to be the lightest glueball. The same scalar states are chosen on the basis of a theoretical model in , except for the $`a_0(980)`$, the existence of a glueball is not accepted though . A similar scheme with a broad state as glueball is considered in , albeit at a higher mass; in this scheme the $`f_0(1370)`$ is taken as $`q\overline{q}`$. A glueball around 1000 MeV which mixes with a nonstrange isoscalar is also obtained in the QCD sum rule approach .

In order to distinguish these possibilities it is of primary importance to verify the existence of the states in question and to determine their constituent structure. As relevant criterion for a resonant state to be included in spectroscopy we require the amplitude to describe a full circle in the complex plane (“Argand diagram”), possibly distorted by a smoothly varying background amplitude. This requires the verification of the appropriate phase motion of the amplitude through an analysis of angular distributions.

## 2 Is there a scalar $`K\pi `$ resonance $`\kappa (900)`$?

Important new information to this discussion has been provided by the FOCUS collaboration on the semileptonic decay $`D^+\mu ^+\nu K^{}\pi ^+`$ . From the forward backward asymmetry in the $`K\pi `$ decay they conclude on the presence of an S-wave component in the mass region studied from 0.8 to 1.0 GeV with constant phase $`\varphi _S=45^{}`$ interfering with a $`K^{}`$ Breit-Wigner resonance. Because of the Watson theorem this phase should equal the elastic $`K\pi `$ scattering phase. This has been determined through the study of $`K^{}pK^{}\pi ^+n`$ by isolating the one-pion-exchange contribution . Indeed, the two results closely agree, the elastic phase shift slowly rises from $`35^{}`$ to $`50^{}`$ in the mass range considered. An S-wave resonance $`\kappa (900)`$ with width around 400 MeV would lead to a phase variation of $`50^{}`$ in this mass region apparently not observed by FOCUS. The elastic $`K\pi `$ amplitude in  has been parametrized by a superposition of the $`K_0^{}(1430)`$ and a background as $`S=T(K_0^{})e^{2i\delta _B}+T_B`$ where the background phase rises slowly and reaches $`\delta _B50^{}`$ at the $`K_0^{}`$. This background we do not consider as evidence for an additional resonance to be included in spectroscopy (see also ).

A related process is the hadronic decay $`D^+K^{}\pi ^+\pi ^{}`$ studied by the E791 Collaboration. The Dalitz plot shows again the presence of the $`K^{}(890)`$ and $`K_0^{}(1430)`$. If the S wave is fitted with an energy independent background amplitude and phase factor multiplying $`K^{}`$ no satisfactory description is obtained. With an additional $`\kappa `$ resonance at mass 798 MeV and width 410 MeV a better fit is obtained. Such a result appears to contradict the above FOCUS result with almost constant phase in the region 800-1000 MeV. We expect that a better fit can be obtained as well if the elastic $`K\pi `$ phase shifts are used, parametrized by a background phase and phase factor for the $`K_0^{}`$ both energy dependent as in .

As an additional check of the fit to the Dalitz plot we suggest comparing not only to the mass spectra with fine binning but also to the higher moments of the decay angular distribution, in particular the first moment $`Y_1^0\mathrm{cos}\theta `$ where cos$`\theta `$, for given $`m^2(K^{}\pi _1^+)`$, is related to $`m^2(K^{}\pi _2^+)`$ (see also ). These moments should reflect the angular distribution of the $`K^{}\pi _1^+`$ channel but with a smooth backround from the resonances in the $`K^{}\pi _2^+`$ channel. Inspection of the Dalitz plot suggests again a strong variation of the asymmetry $`\mathrm{cos}\theta `$ over the $`K^{}(890)`$ region. For the time being we see no compelling evidence for an additional $`K\pi `$ resonance below 1 GeV.

## 3 Is there a scalar $`\pi \pi `$ resonance $`\sigma (600)`$ ?

The elastic $`\pi \pi `$ phase shifts are by now rather well known and a unique solution is established up to $``$ 1400 MeV , similar to the old results . If the rapid phase variation due to $`f_0(980)`$ is removed one finds a slowly moving phase passing $`90^{}`$ at around 1000 MeV and another more narrow structure presumably related to $`f_0(1500)`$. The behaviour of this phase shift may be parametrized by a Breit-Wigner resonance of 1 GeV mass with large width of at least 500 MeV ; if additional background is included the resonance position may shift to higher values around 1400 MeV  or lower values around 600 MeV . The question may be asked whether in other reactions, with different background, a resonance around 600 MeV appears with a corresponding phase variation. Several such proposals have been put forward and we consider some of them in the following.

1. $`J/\psi \omega \pi \pi `$
There is a peak around 500 MeV in the $`\pi \pi `$ mass spectra which may be a signal from a $`\sigma `$ Breit-Wigner resonance . Then the interference term Re$`(SD^{})`$ between the (almost real) D wave which is dominated by $`f_2(1270)`$ and the resonant S wave would change sign at the mass of the $`\sigma `$ and so the angular distribution $`d\sigma /d\mathrm{\Omega }|S|^2+10(3\mathrm{cos}^2\vartheta 1)`$ Re$`(SD^{})+𝒪(|D|^2)`$ would vary accordingly with a sign change of the $`\mathrm{cos}^2\vartheta `$ term (from + to –). The data  do not show any sign change below 750 MeV and therefore there is no indication for a Breit Wigner resonance at 500 MeV.

2. Central $`\pi \pi `$ production in $`ppp(\pi \pi )p`$
The mass distribution of the pion pair in this double Pomeron dominated process peaks shortly above threshold $`400`$ MeV  and has been related to the $`\sigma (600)`$ as well. However, again, there is no related phase variation of the S wave amplitude which should become visible from the $`SP`$ or $`SD`$ interference terms. In this process we understand the origin of the peak. We propose this process to be dominated at low masses by Pomeron Pomeron $`\pi \pi `$ through one pion exchange very much like $`\gamma \gamma \pi \pi `$, the latter process is discussed in . Indeed there is a close similarity between these two processes: the $`I=0`$ component peaks below 400 MeV and the D-wave already near 500 MeV with 1/3 of intensity. As the pion pole is near the physical region the $`\pi \pi `$ angular distribution is very steep, steeper than in more typical interactions mediated by vector ($`\rho `$) exchange. Therefore one estimates that the D wave becomes important not at $`m_{f_2}`$ but already at $`m_{f_2}\times (m_\pi /m_\rho )0.3`$ GeV. This mechanism also explains the low mass peak of the S wave without associated phase variation. The production of a broad state at 1000 MeV as in elastic $`\pi \pi `$ scattering is possible in addition either by rescattering or by direct formation as in $`\gamma \gamma \pi \pi `$ .

3. Decay $`D^+\pi ^{}\pi ^+\pi ^+`$
The $`\pi ^+\pi ^{}`$ mass spectrum presented by the E791 Collaboration  shows three prominent peaks, one just above $`\pi \pi `$ threshold, one related to $`\rho `$ and one to $`f_0(980)`$. Only fits including a light $`\sigma `$ particle have been found successful according to their analysis. In analogy to the decay $`D^+K^{}\pi ^+\pi ^+`$ discussed above we would expect the low mass region to be governed by the elastic $`\pi \pi `$ scattering phase without additional resonance contributions. This should apply strictly to the corresponding semileptonic decay $`D^+\pi ^{}\pi ^+\mu ^+\nu `$. As discussed above for the $`\kappa `$ a resonant $`\sigma (600)`$ should yield a characteristic interference pattern in the projected $`\mathrm{cos}\theta `$ moment which should be analysed. For the moment we consider the question of the $`\sigma (600)`$ in this process as open.

In conclusion, the first two processes show peaks at low mass but definitely no resonant phase motion, for the latter process this question is not definitely answered. We do not discuss here other peaks related to $`\sigma \pi \pi `$ ($`\psi ^{}J/\psi \pi \pi `$, $`Y^{},Y^{\prime \prime }Y\pi \pi `$, $`\tau \nu _\tau 3\pi `$, $`f_0(1370/1500)4\pi `$) which are all lacking a phase analysis. For the moment we suppose the $`\pi \pi `$ phase shifts in all these processes could behave as in elastic $`\pi \pi `$ scattering with a broad state around 1 GeV but no narrower state below 1 GeV.

## 4 Has $`f_0(1500)`$ a strong glueball component?

The main argument  against a strong glueball component in $`f_0(1500)`$ is the observed negative relative phase between the amplitudes

$$T(\pi \pi f_0K\overline{K})=T(\pi \pi f_0\eta \eta ).$$
(1)
These amplitudes have been reconstructed from the measured $`|S|,|D|`$, their relative phase and the absolute phase of the D wave resonances. A similar behaviour of amplitudes (although with different overall phase) is also found in the fits by . This constraint strongly restricts the possible admixture of a glueball component which would contribute with the same sign to all pseudoscalar particle pairs. In particular, among the six models listed in Table 4 of  which describe the $`f_0(1500)`$ as superposition of $`u\overline{u}+d\overline{d},s\overline{s}`$ and glueball only two  yield a negative sign between the amplitudes (1) with a glueball amplitude of 0.22 and 0.01 resp. or a probability of less than 4% for the glueball component.

## 5 The candidate scalar glueball

Having constructed the low mass $`q\overline{q}`$ nonet with $`f_0(980),f_0(1500),a_0(980)`$ and $`K^{}(1430)`$  and not accepting $`\sigma `$ and $`\kappa `$ as genuine resonances below 1 GeV the states left in in our list at low mass are $`f_0(4001200)`$ and $`f_0(1370)`$. We do not discuss here $`f_0(1370)`$, a full circle in the Argand diagram has not been established in our view , also there are problems with an inconsistency of branching ratios . These two objects we consider as representing one single broad resonance as suggested in particular by elastic $`\pi \pi `$ scattering and this is our glueball candidate $`gb(1000)`$.

This state fulfils most standard requirements on glueballs:
1. central production in $`pp`$ scattering;
2. decay of radially excited states $`\psi ^{},Y^{},Y^{\prime \prime }`$ into the respective ground state and $`\pi \pi `$, this can only proceed through gluonic intermediate states;
3. Low energy $`p\overline{p}3\pi `$;
4. $`J/\psi \gamma \pi \pi `$: no prominent signal is observed here up to now which is the only unfavourable aspect of our hypothesis;
5. Flavour properties: The decays of the glueball component $`f_0(1370)`$ favours a glueball over a nonstrange $`q\overline{q}`$ assignment (see, for example, );
6. Small coupling to $`\gamma \gamma `$, see our discussion ; 7. t-channel analysis of elastic $`\pi \pi `$ scattering: The isoscalar exchange amplitude below 1 GeV cannot be saturated by $`q\overline{q}`$ resonances alone .

## 6 Conclusions

We found no definitive evidence for the low mass $`\kappa `$ and $`\sigma `$ Breit-Wigner resonances below 1 GeV, in particular, the semileptonic $`D`$ decays speak against the $`\kappa `$; similar clarity on the $`\sigma `$ could be obtained from semileptonic decays into $`\pi \pi `$. The $`f_0(1500)`$ can only have a small glueball component. The broad state $`gb(1000)`$ is a realistic glueball candidate. Our analysis is mainly phenomenological and driven by simplicity, this does not exclude more complex scenarios, such as a mixing between $`f_0(980)`$ and $`gb(1000)`$. Further experimental results would be helpful, for example the comparative study of quark and gluon jets in their respective fragmentation regions .