text:
# 1 Introduction

## 1 Introduction

The debate over what stabilises complex food webs has been an active one for over a quarter of a century (McCann, 2000). While it is found that the stability of a randomly linked food web model tends to decrease with the proportion of links and the number of species (May, 1972), real food webs display a high degree of stability, in spite of being very complex. Clearly, real food webs have features that are ignored in randomly linked models. For example, the distribution and strength of links in real food webs are far from random. More realistic approaches use link values taken from real food webs (Yodzis, 1981), or assemble a food web by repeated addition of species from a large species pool that contains different types of species such as “plants”, “herbivores”, “carnivores”, and “top predators”(Morton & Law, 1997; Law, 1999). Using Lotka-Volterra type dynamics, assembly models of this type lead to large complex webs. These Lotka-Volterra type models lack nonlinear effects such as predator saturation, which many consider important for food web dynamics. The species pool is usually composed of a large number of species which are generated in an ad hoc manner, by pre-defining the possible types and value ranges of links. A complementary approach to food web stability focuses on more realistic dynamics (that include predator saturation, for instance) of small systems consisting only of a few species. These small systems can be viewed as a small part of a food web. It has been found that the presence of many weak links and the use of nonlinear growth rates, that prevent predators from feeding successfully on many prey at the same time, stabilise the dynamics of these small systems by reducing population fluctuations (McCann et al., 1998). Again, the possible types and values of links are pre-defined in these studies.

In real food webs, the link strengths and the linkage pattern can change with time through invasion and replacement with related species, and through adaptive evolutionary processes (Thompson, 1998). It is therefore desirable to investigate food web models that implement such long-term changes due to modification of existing link strengths and linkage patterns. In the following, we call this long-term change of the food web “evolution”, with the meaning of “development of the food web structure in time”, while the “true” biological evolution is of course happening on the level of individuals. Our approach is different from the species assembly models insofar that the new species are not picked irrespective of the web composition, but are similar to existing species.

Considering the long-term change of food webs adds a new dimension to the complexity-stability debate: even though species die out and are replaced with others during the course of time, ecosystems as a whole persist in time. Except for rare catastrophic events due to external causes, food webs preserve their general structure. This means that typical structural features of food webs, such as the total number of species or the number of trophic layers, do not collapse to zero, but rather fluctuate around some mean value far away from zero. This definition of stability is perhaps more relevant for understanding the complexity of ecosystems, than some of the other definitions found in the literature. Persistence in time of a complex ecosystem of course also requires a certain kind of stability of the population dynamics: while population dynamics does not need to reach fixed points, oscillations must be small enough such that the population sizes of the majority of species remain positive after the introduction of a mutant or invading individual. If this were not so, each change in the web composition would lead to a collapse of the web, and complex ecosystems could never arise.

In this paper, we therefore study evolutionary food web models and investigate under which conditions these models show a persistent complex structure with several trophic layers. The evolutionary model that we use was studied previously with specific types of population dynamics (Caldarelli et al., 1998; Drossel et al., 2001). Here, we explore the effect of changing the functional response used in the model. We chose a representative selection of functional responses found in the literature and modified them such that they are suitable for many-species communities. Our aim here is to show that biologically realistic requirements, in particular the capability of predators to adjust their feeding rates and focus on their best prey, are essential to the construction of large stable complex webs.

## 2 The model

The governing equation for the population $`N_i(t)`$ of species $`i`$ is taken to have the form

$`{\displaystyle \frac{dN_i(t)}{dt}}`$ $`=`$ $`\lambda {\displaystyle \underset{j}{}}N_i(t)g_{ij}(t){\displaystyle \underset{j}{}}N_j(t)g_{ji}(t)`$ (1)
$`d_iN_i(t){\displaystyle \underset{j}{}}\alpha _{ij}N_i(t)N_j(t),`$
which naturally accounts for the four processes which lead to a change in the population number $`N_i(t)`$. The function $`g_{ij}(t)`$ is the functional response, i.e. the rate at which an individual of species $`i`$ feeds on species $`j`$; it depends on the population sizes, and its analytical form will be specified below. The constant $`\lambda `$ is the ecological efficiency at which consumed prey are converted into predator offspring. The first term describes population growth due to food consumption, while the second term describes population decline due to predation. The death rate, $`d_i`$, will be assumed to be equal to 1 for all species. The last term describes direct interference competition between two predators $`i`$ and $`j`$. It is needed to represent intra-specific interference competition, which ensures bounded solutions when the functional response is independent of predator density. This term also allows the incorporation of inter-specific interference competition into such models, which facilitates comparison with the predator dependent functional responses.

On a larger time scale, species undergo changes, either by replacement with invading similar species, or by evolutionary change. In order to implement these changes, each species was characterised by a set of 10 out of 500 possible features (Caldarelli et al., 1998). This representation gives a measure of similarity between species (the number of features they have in common) and allows for “mutations” or “invasions” by randomly replacing one feature of one individual with another. Scores between features were assigned in a random and asymmetric fashion. They are a measure of how useful a feature is for its carrier at feeding on a species carrying another feature. The interactions, or scores, $`a_{ij}`$, between species were obtained by adding the mutual scores between all pairs of features carried by the two species. Positive scores indicate that the first species can feed on the second species and negative scores mean that the first species is consumed by the second. In the latter case the score was set to zero, since this information was already contained in the (positive) score where the two species are interchanged. Therefore all scores were such that $`a_{ij}0`$. The external resources were represented as an additional species of fixed (and large) population size, which does not feed on any species. The evolutionary dynamics consisted of the following steps: starting from some initial species configuration (usually the external resources plus one basal species feeding on it), one individual was picked at random, and one of its features was changed. Under the dynamics (1), this new species either died out, or added to the system, or drove one or several other species to extinction. When the dynamics had reached an equilibrium, the next “mutation” occurred, and the process was repeated.

In the following, we list the different functional responses used in our computer simulations. Lotka-Volterra dynamics are recovered by taking

$$g_{ij}=a_{ij}N_j.$$
(2)
A Holling type II functional response, which implies saturation of consumption rates at high prey abundance, is given by

$$g_{ij}=\frac{a_{ij}N_j}{1+_kb_{ik}N_k},$$
(3)
where the sum in the denominator is taken over all prey $`k`$ of species $`i`$.

More complicated functional responses for multi-species systems can only be found in the recent literature. Arditi and Michalski (1996) suggest the following generalised Beddington form:

$$g_{ij}=\frac{a_{ij}N_j}{1+_kb_{ik}N_k+_lc_{il}N_l},$$
(4)
where the first sum is again taken over all prey $`k`$ of species $`i`$, and the second sum is taken over all those predator species $`l`$ that share a prey with $`i`$. We chose the $`c_{il}`$ such that individuals belonging to the same species competed more strongly with each other than individuals belonging to different species.

Arditi and Michalski (1996) also suggest the following ratio-dependent functional response, which implements the idea that predators share the prey:

$$g_{ij}=\frac{a_{ij}N_j^{r(i)}}{N_i+_{kR(i)}b_{ik}N_k^{r(i)}},$$
(5)
with the self-consistent conditions

$$N_j^{r(i)}=\frac{\beta _{ji}N_i^{C(j)}N_j}{_{kC(j)}\beta _{jk}N_k^{C(j)}},N_k^{C(j)}=\frac{h_{jk}N_j^{r(k)}N_k}{_{lR(k)}h_{lk}N_l^{r(k)}}.$$

Here $`\beta _{ij}`$ is the efficiency of predator $`i`$ at consuming species $`j`$, $`h_{ij}`$ is the relative preference of predator $`i`$ for prey $`j`$, $`R(i)`$ are the prey species for predator $`i`$, $`C(i)`$ are the species predating on prey $`i`$, $`N_j^{r(i)}`$ is the part of species $`j`$ that is currently being accessed as resource by species $`i`$ and $`N_k^{C(j)}`$ is the part of species $`k`$ that is currently acting as consumer of species $`j`$.

Finally, in a recent paper (Drossel et al., 2001), we used the ratio-dependent expression

$$g_{ij}(t)=\frac{a_{ij}f_{ij}(t)N_j(t)}{bN_j(t)+_k\alpha _{ki}a_{kj}f_{kj}(t)N_k(t)},$$
(6)
where $`f_{ij}`$ is the fraction of its effort (or available searching time) that species $`i`$ puts into preying on species $`j`$. These efforts are determined self–consistently from the condition

$$f_{ij}(t)=\frac{g_{ij}(t)}{_kg_{ik}(t)}.$$
(7)
This condition is such that no individual can increase its energy intake by putting more effort into a different prey. The sum in (6) is over all species $`k`$ which are predators of $`j`$. The competition strength $`\alpha _{ik}`$ is set equal to one only for $`i=k`$, and is smaller than 1 otherwise. For this model, or for any other model with a functional response which depends on predator density, the last term in Eq.(1) should be dropped.

## 3 Simulation results: food web structure

Let us first present results for the simplest type of population dynamics, which is Lotka-Volterra dynamics, Eq. (2). Starting with one species and the external resources, input at a rate $`R`$, we allowed for enough time for the system to evolve towards its characteristic stationary structure. A typical food web resulting after sufficiently long time is shown in Fig. 1. All species feed on the external resources, and additionally on other species. Only for short transient periods of time does the web have species on the second trophic level, and virtually never on the third. We also performed simulations that started from a complex web with several trophic layers, which was stable under the population dynamics. However, the evolutionary dynamics caused this structure to collapse, and the resulting structure was again similar to that shown in Fig. 1. These results were obtained without inter-specific direct competition, i.e. with $`\alpha _{ij}=c\delta _{ij}`$ where $`\delta _{ij}`$ is the Kronecker delta defined by $`\delta _{ij}=0`$ for $`ij`$ and $`\delta _{ii}=1`$. They are still true when direct inter-specific competition is included, as in a different evolutionary model described by Lässig et al (2001) and also in our model if we set $`\alpha _{ij}=c\rho _{ij}(1+q_{ij})/2`$ in Eq. (1). Here $`\rho _{ij}`$ is equal to one if $`i`$ and $`j`$ share at least one prey and is zero otherwise, and $`q_{ij}`$ is the fraction of discrete features shared by the two species.

Our interpretation of these results is as follows. When a modified species that can feed on more prey than a parent species is introduced into the model, the modified species takes over. This is due to the linear growth rates of Lotka-Volterra systems that allow a species to feed on each of several prey just as easily as it would feed on one prey, if only this prey was present. If allowed to evolve, an artificially constructed Lotka-Volterra food web will therefore tend to a simple web structure where each species feeds on many other species and in particular on the most abundant resource, the external environment. Clearly, this feature of the model is unrealistic as there are evolutionary trade-offs that do not allow species to feed very easily on an arbitrarily large number of species. In particular, species cannot normally feed on external resources and on many other types of prey.

In the light of these arguments, a Holling or Beddington type functional response is more realistic, since it naturally limits the total amount of prey taken by a predator. Indeed, the need for nonlinear growth rates in the context of multispecies communities has often been discussed in the literature (Pimm, 1991; Rosenzweig, 1995; Vandemeer et al., 2002). Functional responses of the form (3) were used by McCann *et al* (1998) in their investigation of the effects of weak links. They have the biologically meaningful feature that predators cannot maintain a high feeding efficiency on many prey at the same time, and they have a stabilising effect on the dynamical equations, at least for subsystems consisting of up to four species (McCann et al., 1998). Additional direct inter-specific competition can be introduced in a natural way by using functional responses of the Beddington form (4). However, when evolutionary dynamics by modification of existing species are added to either of these models, we find again that complex webs with several trophic layers are not stable and that after some time the webs consist merely of basal species, with a structure very similar to that shown in Fig. 1.

It appears that even for the Holling or Beddington type functional response too many species can feed on the same food source, and that the model does not have a mechanism that prevents possible links (i.e. $`a_{ij}>0`$) from being realised. In the models we are considering here, there is always a possible predator-prey interaction between any two species $`i`$ and $`j`$, either $`a_{ij}>0`$ and $`a_{ji}=0`$ or $`a_{ji}>0`$ and $`a_{ij}=0`$. Thus the ecosystems are potentially highly connected with one half of all inter-specific links realised. However, species usually focus on those prey to which they are best adapted. We therefore must introduce a mechanism which would limit the number of species a predator actually feeds on, with potential links becoming active only when the preferred prey is scarce or unavailable, or after a change in the composition of the web. One way to (artificially) realise this feature is to allow only the best predators to feed on a given prey. We did this by replacing the interactions $`a_{ij}`$ with adjusted interactions $`a_{ij}^{}=a_j^{\mathrm{max}}(1(a_j^{\mathrm{max}}a_{ij})/\delta )`$ (Caldarelli et al., 1998), with $`\delta `$ being a small parameter and $`a_j^{\mathrm{max}}`$ being the largest interaction against $`j`$. Negative $`a_{ij}^{}`$ are set to zero. Fig. 2 shows a food web obtained from Lotka-Volterra dynamics with this artificial constraint. Table 1 shows the mean number of species, of links per species, and the mean occupation numbers of the trophic levels for this model, for $`R=1\times 10^4`$, $`\lambda =0.1`$, $`c=1.0`$ and $`\delta =0.2`$. For comparison, the results for the original Lotka-Volterra model with the parameters for $`R=1\times 10^4`$, $`\lambda =0.1`$, and $`c=3.0`$ are also shown. Similar web structures were obtained including direct inter-specific competition or using Holling and Beddington forms with the same type of adjusted interactions.

This shows that the ability of predators to concentrate on the prey that they are best suited to exploit, rather than on all possible prey, is essential for the production of realistic food web structures.

Finally, let us discuss evolutionary food web models with ratio-dependent functional responses. Ratio-dependent functional responses naturally limit the number of actual prey of a predator by dividing the available prey among its predators. One such ratio-dependent form, given by (5), was suggested by Arditi and Michalski (1996). However, it leads to such strong competition that only one predator species can feed on a given prey species and thus cannot give rise to complex food web structures. (One can conclude this directly from the published webs (Arditi & Michalski, 1996) or from an analytical calculation with two predator species and one prey species). In reality, however, different species exploit a food source in different ways and therefore competition between individuals of different species is less than between individuals of the same species.

This feature is implemented in the functional response (6), which was introduced by some of us (Drossel et al., 2001). Furthermore, equations (6) and (7) have several other biologically meaningful features: when prey is very abundant, the consumption rates saturate due to the term containing $`b`$ in the denominator. Additionally, predators divide their effort (or time) among the available prey in order to obtain the maximum possible amount of food. The condition (7) is an evolutionarily stable strategy and allows species to swiftly adjust to a modified situation, thus endowing the food web model with the flexibility found in real webs. Finally, equations (6) and (7) are invariant if identical species are aggregated into a single species. This feature is shared by the other models used in this paper if the last term in equation (1) is either set to zero or present for any pair of species but not if, for instance, only inter-specific competition is included.

Using the functional response (6) in the evolutionary model, we rapidly obtain large and complex food webs with a stable structure, although there is an ongoing turnover of species in the system. The dynamics (1) converge quickly to a fixed point. The food webs generated by this model share many features with real food webs, such as the fraction of top, intermediate and basal species, and the mean number of links per species. This is discussed in detail elsewhere (Drossel et al., 2001; Quince et al., 2004a).

The webs generated by this model often preserve their structure under changes in the species composition. Fig. 3 shows a web before and after a basal species was deleted. One can see that in this case all other species manage to survive, with the strength of several links being modified, illustrating the flexibility and resilience of the model food webs. Of course, the deletion of a randomly chosen species is not always buffered in this perfect way. A useful measure of the ecosystem stability to deletion is the fraction of species that can be removed without causing further extinctions (Pimm, 1979). This statistic depends on the model parameters; for the values used to generate Fig. 3 it is 62%. For the deletions which do cause further extinctions the distribution of event sizes decays roughly exponentially with a characteristic size of just a few species (Quince et al., 2004b).

## 4 Simulation results: link strength distributions

There has been an increasing realisation that food webs have a large proportion of weak links (Paine, 1992; Tavares-Cromar & Williams, 1996; Berlow et al., 1999; Neutel et al., 2002), and that weak links tend to stabilise population dynamics (McCann et al., 1998). For this reason, we evaluated the link strength distribution in our model webs. We found that for all functional responses that generated large ecosystems either the link strength distributions were skewed towards zero or a large fraction of links were zero. This is a highly non-trivial result as in contrast to other work on the topic of weak links, the link strength distributions in our model are an emergent property of the system and not put in by hand. It is a strong indication that weak links are the natural outcome of long-term ecosystem evolution coupled to population dynamics.

In the literature there are different definitions of link strength, ranging from the biomass flowing through a link per unit time to the response of one population size to a small change in another population size (Laska & Wootton, 1998). Here, we evaluated link strength distributions based on the following two definitions: (i) the proportion of the prey in the predator’s diet after the population dynamics had equilibrated,

$$l_{ij}^{(1)}=\frac{g_{ij}}{_kg_{ik}},$$
(8)
and (ii) the per-capita interaction strength $`l_{ij}^{(2)}`$. These are sometimes known as the elements of the “community matrix” (Laska & Wootton, 1998), although this phrase is also used to refer to the unnormalised matrix of partial derivatives (May 1973). It quantifies the strength of all direct interactions, both predator-prey interactions and interference competition. It is defined as:

$$l_{ij}^{(2)}=\frac{1}{N_i}\left(\frac{\dot{N}_i}{N_j}\right).$$
(9)
Definition (i) typically leads to a U-shaped link-strength distribution, with a peak at weak links and another one near the maximum link strength 1, indicating that many predators have one main prey. (Note that the second peak does not occur if the link strength is defined as the biomass passing through a link, since in this case the link strength is not normalised to $`_jl_{ij}=1`$). Definition (ii) typically leads to broad distributions with or without a peak at zero.

Fig. 4 shows the distribution of link strengths $`l_{ij}^{(1)}`$ for the two Lotka-Volterra models (without and with adjusted scores), and for the model with the ratio-dependent functional response, averaged over many evolutionary time steps and several different model webs. For small link strength, the distribution for the ratio-dependent functional response resembles a power law with an exponent around $`0.74`$, while the Beddington model with adjusted scores shows an exponent around $`0.9`$, and the Lotka-Volterra model with adjusted scores shows an exponent around $`0.29`$ (if a power law fit should be attempted at all). The Lotka-Volterra model without adjusted scores does not have several trophic layers, and its link strength distribution appears to decay exponentially at small values. For the Beddington functional response without adjusted scores, the decay is much steeper than with adjusted scores, but not exponential. All the above mentioned models are skewed towards small values except for the Lotka-Volterra model with adjusted scores where the approach to the origin is quite flat. For this model however only a small proportion, 3%, of links have $`l_{ij}^{(1)}>0`$. We conclude that for the models capable of generating large ecosystems only a small fraction of links were realised or $`l_{ij}^{(1)}`$ had a large weight at small values. We can summarise this by calculating the fraction of links with $`l_{ij}^{(1)}<0.01`$ including zero, which gives 91% and 97% for the Lotka-Volterra models (without and with adjusted scores) and 97% for the ratio-dependent model.

We can investigate these ideas further by examining the link strength according to the definition (ii); this can be either positive or negative. In the case of the Lotka-Volterra models, the only direct inter-specific interactions are between predators and prey with $`l_{ij}^{(2)}=\lambda a_{ij}`$ and $`l_{ji}^{(2)}=a_{ij}`$ if $`i`$ predates $`j`$. Thus in Fig. 5 we only need to show the positive half of the $`l_{ij}^{(2)}`$ distribution. From this we see that, whereas the Lotka-Volterra model without adjusted scores has a distribution of $`l_{ij}^{(2)}`$ values that is skewed towards zero, the model with adjusted scores has a maximum at an intermediate value of $`l_{ij}^{(2)}`$. However, as was the case for $`l_{ij}^{(1)}`$, we find a much higher proportion of non-zero links in the model without adjusted interactions (51% as compared to 3%).

If we introduce direct inter-specific competition into the Lotka-Volterra models as in Section 3 by setting $`\alpha _{ij}=c\rho _{ij}(1+q_{ij})/2`$ in Eq. (1), where $`\rho _{ij}`$ is equal to one if $`i`$ and $`j`$ share at least one prey and is zero otherwise and where $`q_{ij}`$ is the fraction of discrete features shared by the two species, we find that, even if adjusted interactions are used, only small webs can be evolved. The reason for this seems to be that the resulting direct competitive interactions which are allowed by the adjusted interactions are both strong and non-zero for a large fraction of species pairs in the web. This is illustrated by Fig. 6 where the distributions of the quantity $`\alpha _{ij}`$ when $`c=5.0`$ are plotted for the model with and without adjusted scores. We plot $`\alpha _{ij}`$ because it is the contribution of direct competition to $`l_{ij}^{(2)}`$. Note that because $`q_{ij}`$ is discrete so is $`\alpha _{ij}`$ and that the range of non-zero $`\alpha _{ij}`$ is $`5\alpha _{ij}2.5`$.

The model with the ratio-dependent functional response, Eq. (6), also includes direct inter-specific competition, as a result of the implicit sharing of prey between predators. However this does not lead to the strong competitive interactions seen in Fig. 6, probably because this competition occurs through the predator-prey terms, which as was shown in Fig. 4, are themselves very diffuse. This can be seen in Fig. 7, where the distribution of $`l_{ij}^{(2)}`$ values for all competitive interactions and where the shared prey comprised greater than 1% of both competitors diet, is shown. This distribution, although bimodal, is heavily weighted towards the origin. This then explains why the ratio-dependent model with inter-specific competition is, in contrast to the Lotka-Volterra models discussed immediately above, capable of evolving large complex ecosystems.

This finding of ours is complemented by a recent paper by Kondoh (2003), who investigates randomly linked food web models and the cascade model based on Lotka-Volterra dynamics with adjustable foraging efforts. He finds that incorporating adaptive foragers turns a negative relationship between stability and complexity, defined in terms of species number and potential connectance, into a positive one. This result had already been found, albeit over a more limited range of parameters, by Pelletier (2000). These results are strictly different to ours in that the stability criterion, community persistence defined as fraction of species surviving for some long time in a stochastic environment, is different. However the mechanism that allows the generation of large webs with high potential connectance, the reduction in the number of non-zero realised links and the reduction in strength of those links that are realised due to the foraging dynamics, is very similar to that proposed above. In fact this mechanism is arguably just May’s original hypothesis, that strong links or a large fraction of links will destabilise ecosystems, placed in an adaptive context (May, 1972).

## 5 Discussion

The form that a realistic functional response might take has been the subject of a large number of papers in the literature but few, if any, of these suggestions were implemented in a model of multispecies communities in order to test their effectiveness. In this paper we have used an existing model to investigate the effect that different choices for the functional response have on food web structure. We studied models with a wide range of functional responses: Lotka-Volterra with and without direct inter-specific competition, Holling and Beddington forms, all of these but with a mechanism to limit the number of species a predator actually feeds on, and ratio-dependent functional responses. We found that in the first two cases complex webs could not be built up, but in the last two cases they often could. More specifically, unless an (artificial) mechanism was introduced which restricted predator choice to prey which they were best suited to exploit or alternatively a ratio-dependent functional response was used, then stable webs consisting of more than one trophic level could not be built up.

Our second major conclusion was that, given we chose a functional response which gave rise to complex, stable webs as described above, the link strength distributions were skewed towards zero. For the Beddington functional response with adjusted scores or ratio-dependence, the link strength distribution followed a power-law, but this was much less clear in the case of the Lotka-Volterra model, where only a small fraction of the links were realised i.e. a large fraction of the links were zero. Some insight into the mechanisms involved in these effects was also gained by looking at the strength of a direct inter-specific competition introduced into the Lotka-Volterra equations. Without adjusted scores, interactions were present, but not particularly strong. However with adjusted scores, about three-quarters of the interactions (for the parameter values we used) were zero, but the remaining quarter were strong. In this case we found that, even with adjusted scores, large complex webs could not be grown. We concluded that the existence of strong competitive interactions might also destabilise food webs.

In summary, we found that the type of functional response used in the population dynamics of multispecies communities has to be chosen carefully if a large complex community is to be sustained. If this is achieved, then a large proportion of weak links arises naturally from the evolution of the food webs. These conclusions were arrived at within the context of a class of evolutionary food web models. It would be interesting to investigate if these results could be obtained from other starting points or with different model assumptions.

Acknowledgements: We wish to thank Paul Higgs for useful discussions. CQ thanks the EPSRC (UK) for financial support during the initial stages of this work.

## References

Arditi, R., Michalski, J., 1996. Nonlinear food web models and their responses to increased basal productivity. In Food webs: Integration of patterns and dynamics. Polis, G. A., Winemiller, K. O. (eds), pp 122–133. Chapman and Hall, New York. Berlow, E. L., Navarrete, S. A., Briggs, C. J., Power, M. E., Menge, B. A., 1999. Quantifying variation in the strengths of species interactions. Ecology 80, 2206–2224. Caldarelli, G., Higgs, P. G., McKane, A. J., 1998. Modelling coevolution in multispecies communities. J. Theor. Biol. 193, 345–358. Drossel, B., Higgs, P. G., McKane, A. J., 2001. The influence of predator-prey population dynamics on the long-term evolution of food web structure. J. Theor. Biol. 208, 91–107. Kondoh, M. K., 2003. Foraging adaptation and the relationship between food web complexity and stability. Science 299, 1388–1391. Laska, M. S., Wootton, T. J., 1998. Theoretical concepts and empirical approaches to measuring interaction strengths. Ecology 79, 461–476. Lässig, M., Bastolla, U., Manrubia, S. C., Valleriani, A., 2001. Shape of ecological networks. Phys. Rev. Lett. 86, 4418–4421. Law, R., 1999. Theoretical aspects of community assembly. In: Advanced ecological theory: principles and applications. McGlade, J. (ed). pp 143–171. Blackwell, Oxford. May, R. M., 1972. Will a large complex system be stable? Nature 238, 413–414. May, R. M., 1973. Stability and complexity in model ecosystems. Princeton University Press, Princeton. McCann, K., 2000. The diversity-stability debate. Nature 405, 228–233. McCann, K, Hastings, A., Huxel, G. R., 1998. Weak trophic interaction and the balance of nature. Nature 395, 794–798. Morton, R. D., Law, R., 1997. Regional species pools and the assembly of local ecological communities. J. Theor. Biol. 187, 321–331. Neutel, A-M., Heesterbeek, J. A. P., de Ruiter, P. C., 2002. Stability in real food webs: weak links in long loops. Science 296, 1120–1123. Paine, R. T., 1992. Food-web analysis through field measurements of per capita interaction strength. Nature 355, 73–75. Pelletier, J. D., 2000. Are large complex ecosystems more unstable? A theoretical reassessment with predator switching. Math. Biosc. 163, 91–96. Pimm, S. L., 1979. Complexity and stability: another look at MacArthur’s original hypothesis. Oikos 33, 351–357. Pimm, S. L., 1991. The balance of nature. The University of Chicago Press, Chicago. Quince, C., Higgs, P. G., McKane, A. J., 2004a. Topological structure and interaction strengths in model food webs. Submitted for publication to Ecol. Model. Quince, C., Higgs, P. G., McKane, A. J., 2004b. Deleting species from model food webs. Submitted for publication to Oikos. Rosenzweig, M. L., 1995. Species diversity in space and time. Cambridge University Press, Cambridge. Tavares-Cromar, A. F., Williams, D. D., 1996. The importance of temporal resolution in food web analysis: evidence from a detritus-based stream. Ecol. Monogr. 66, 91–113. Thompson, J. N., 1998. Rapid evolution as an ecological process, TREE 13, 329–332. Vandemeer, J. *et al*., 2002. Increased competition may promote species coexistence. Proc. Nat. Acad. Sci. 99, 8731–8736. Yodzis, P., 1981. The stability of real ecosystems. Nature 289, 674–676.


format: text

text:
# Vacuum creation of quarks at the time scale of QGP thermalization and strangeness enhancement in heavy-ion collisions

## I Introduction

Multiple particle production in relativistic heavy-ion collisions is not fully understood yet. One of the most popular microscopic approaches to this phenomenon is formulated within the chromoelectric flux tube (or rather string) model CNN79 . The flux tubes are assumed to be initially stretched between the constituents (quarks and diquarks) of the nucleons of colliding nuclei. As the constituents are flying away, the energy of a color tube increases making it unstable against the production of $`q\overline{q}`$ or $`d\overline{d}`$ pairs from the vacuum. New hadrons, created in the course of the color tube fragmentation, can also interact within hot and dense nuclear matter, and the in-medium cascade develops.

The formation and subsequent break-up of color tubes/strings is a common feature adapted by all microscopic string models (see, e.g., lund ; Kaid99 ; qgsm ; venus ; urqmd  and references therein), which are intended to describe dynamics of heavy-ion collisions at relativistic energies. For a uniform chromoelectric field $`E`$ the probability to create a pair of quarks with mass $`m`$, effective charge $`e`$, and transverse momentum $`p_t`$ per unit time and per unit volume reads CNN79

$$P(p_t)d^2p_t=\frac{|eE|}{4\pi ^3}\mathrm{ln}\left\{1\mathrm{exp}\left[\frac{\pi (m^2+p_t^2)}{|eE|}\right]\right\}d^2p_t.$$
(1)
The integrated probability is given by series expansion

$$P_m=\frac{(eE)^2}{4\pi ^3}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1}{n^2}\mathrm{exp}\left(\frac{\pi m^2n}{|eE|}\right),$$
(2)
which reproduces the classical Schwinger result schw  derived in spinor quantum electrodynamics (QED) for $`e^+e^{}`$ production rate in the constant electric field. Usually, only the leading term in Eq. (2) is taken into account. According to this formula, the ratio of production rates of strange to nonstrange quark pairs, widely known as strangeness suppression factor $`\gamma _s`$, is

$$\gamma _s=\frac{P(s\overline{s})}{P(q\overline{q})}=\mathrm{exp}\left[\frac{\pi (m_s^2m_q^2)}{\kappa }\right],$$
(3)
where $`\kappa =|eE|`$ is the so-called string tension. It appears that microscopic models underestimate the yield of strangeness in ultrarelativistic heavy-ion collisions low\_str . This issue is very important, because the abundant yield of strange particles was predicted raf91  as one of the signals of quark-gluon plasma (QGP) creation. According to Eq. (3), the strangeness production can be enhanced either by taking into account string-string interaction, which leads to fusion of strings and formation of color ropes with larger effective string tension bnk84 ; sor92 ; am93 ; soff99 , or by dropping the quark masses blei00 , e.g., due to chiral symmetry restoration. Also, the effects related to the finiteness of the strings can modify the production rates WW88 ; fbs .

In the present Letter we explore another possibility: The real fields emerging in heavy-ion collisions act at the time interval comparable with the Compton scale. Since the dynamics of particle creation in time-dependent homogeneous fields differs from that of stationary fields grib , it is essential to properly modify the system description at early stage of nuclear collisions. The kinetic equation (KE) is a convenient tool to study the nonequilibrium evolution processes. The source term describing the vacuum pair creation process can be incorporated into the KE either in a phenomenological manner gatoff ; asakawa  on the basis of Schwinger-like formula, or derived in a more sophisticated manner from the microscopic equations of motion smol . Although it is believed that the phenomenological source term correctly reproduces qualitative features of the quantum mean field theory kluger98 , such approximation has yet to be verified for, e.g., time dependent fields or multicomponent systems. Compared to the phenomenological treatment, the approach within the framework of a transport equation smol  contains several new dynamical aspects, such as longitudinal momentum dependence of the distribution functions and non-Markovian character of the time evolution. It takes into account effects of the field switching and statistics, as well. Therefore, the abundances of newly produced particles may considerably deviate from the values obtained for the constant and infinite field. The appearance of the non-Schwinger regime in a creation of electron-positron pairs in periodic laser field has been discussed in gan ; pop ; rob . Another noteworthy feature of the modified source term is the suppression of zero momentum bosons smol1 . This circumstance causes the “fermion dominance” effect at the short time scales skok1  and can lead to the abundant production of heavy fermions at expense of light bosons.

In case of the QGP creation the characteristic time of the field variation is estimated to be of order of few fm/$`c`$ heinz  and, therefore, the assumption of the space-time unvarying field is too crude. Hence, it is necessary to elaborate on the dynamics of particle vacuum production at short time scales, which are compatible with the particle inverse masses. Particularly, the dependence of the production rates on the shape and duration of the field pulse should be investigated.

The paper is organized as follows. Kinetic equation with the source term, describing the vacuum production of fermions, is derived in the collisionless limit in Sec. II. Section III presents study of the influence of the field pulse characteristics on the time-dependent vacuum production rates for the energies of RHIC and LHC. Comparison with the production rates, obtained for the infinite and constant chromoelectric field by the classical Schwinger-like formula, is also performed. Finally, conclusions are drawn in Sec. IV.

## II Kinetic equation

As was mentioned above, the kinetic equation is a standard tool to describe the non-equilibrium evolution of a many body system. The general form of the KE for the distribution function $`f(𝐩,t)`$ in a strong spatially homogeneous time-dependent field is gatoff

$$\frac{f(𝐩,t)}{t}+e𝐄(t)\frac{f(𝐩,t)}{𝐩}=S(𝐩,t)+C(𝐩,t),$$
(4)
where $`S(𝐩,t)`$ and $`C(𝐩,t)`$ are the source term and the collision integral, respectively. The source term describes the vacuum production of particle-antiparticle pairs in the external field, while the collision term governs their rescattering dynamics, which drives the system toward thermal equilibrium. For the sake of clarity, it is relevant to consider the collisionless approximation, where the collision term in the r.h.s. of Eq. (4) is omitted.

The kinetic description of quarks in framework of the Wigner function formalism oh ; qcd  leads to a very complicated system of partial differential equations hard to solve both numerically and analytically. Therefore, following the approach of grib ; smol , we employ the canonical Bogoliubov transformation method within the QED approximation. To derive the source term for the system of fermions in an external electric field, we start from the Dirac equation

$$(i\gamma ^\mu _\mu e\gamma ^\mu A_\mu m)\psi (x)=0.$$
(5)
Using the simple field configuration with vector potential in the Hamilton gauge $`A^\mu =(0,0,0,A(t))`$ and homogeneous electric field $`𝐄(t)=(0,0,E(t))`$, $`E(t)=\dot{A}(t)`$, one looks for the solutions of the Eq. (5) in the form

$$\psi _{𝐤r}^{(\pm )}(x)=\left[i\gamma ^0_0\gamma 𝐤+e\gamma ^3A(t)+m\right]\chi ^{(\pm )}(𝐤,t)R_r\mathrm{e}^{i𝐤\overline{x}},$$
(6)
where the superscript $`(\pm )`$ denotes eigenstates with the positive and negative frequencies. Herein the spinors $`R_r`$ $`(r=1,2)`$ are eigenvectors of the matrix $`\gamma ^0\gamma ^3`$ satisfying the condition $`R_r^+R_s=2\delta _{rs}.`$ The functions $`\chi ^{(\pm )}(𝐤,t)`$ obey the oscillator-type equation

$$\ddot{\chi }^{(\pm )}(𝐤,t)+\left[\omega ^2(𝐤,t)ie\dot{A}(t)\right]\chi ^{(\pm )}(𝐤,t)=0,$$
(7)
where we define the total energy $`\omega ^2(𝐤,t)=\epsilon _{}^2+p_{}^2(t)`$, the transverse energy $`\epsilon _{}^2=m^2+𝐤_{}^2`$, and the longitudinal momentum $`p_{}(t)=k_{}eA(t)`$. The solutions $`\chi ^{(\pm )}(𝐤,t)`$ of Eq. (7) for positive and negative frequencies are fixed by their asymptotic behavior at $`t_0=t\mathrm{}`$, where $`\dot{A}(t_0)=0`$. In this limit, the field operators $`\psi (x)`$ and $`\overline{\psi }(x)`$ can be decomposed by the complete and orthonormalized set of spinor functions (6) as follows:

$$\psi (x)=\underset{r,𝐤}{}\left[\psi _{𝐤r}^{()}(x)b_{𝐤r}(t_0)+\psi _{𝐤r}^{(+)}(x)d_{𝐤r}^+(t_0)\right].$$
(8)
The operators $`b_{𝐤r}(t_0),b_{𝐤r}^+(t_0)`$ and $`d_{𝐤r}(t_0),d_{𝐤r}^+(t_0)`$ describe the annihilation and creation of particles and antiparticles in the in-state $`|0_{in}`$ and obey the standard anticommutation rules. The time evolution leads to the mixing of states with positive and negative energies and, therefore, non-diagonal terms in the Hamiltonian corresponding to Eq. (5) emerge. The diagonalization of the Hamiltonian, which is equivalent to the transition to quasiparticle representation, is performed by the time-dependent Bogoliubov transformation

$`\begin{array}{ccccc}b_{𝐤r}(t)\hfill & =& \alpha _𝐤(t)b_{𝐤r}(t_0)& +& \beta _𝐤(t)d_{𝐤r}^+(t_0),\hfill \\ d_{𝐤r}(t)\hfill & =& \alpha _𝐤(t)d_{𝐤r}(t_0)& & \beta _𝐤(t)b_{𝐤r}^+(t_0),\hfill \end{array}`$ (11)
with the imposed condition $`|\alpha _𝐤(t)|^2+|\beta _𝐤(t)|^2=1.`$ The new operators $`b_{𝐤r}(t)`$ and $`d_{𝐤r}(t)`$ describe the processes of quasiparticle creation and annihilation. By virtue of the Lagrange multipliers, one can find from the equations of motion (7) that the coefficients in the Bogoliubov transformation (11) are connected via the relations grib

$`\begin{array}{ccc}\dot{\alpha }_𝐤(t)\hfill & =& {\displaystyle \frac{eE(t)\epsilon _{}}{2\omega ^2(𝐤,t)}}\beta _𝐤^{}(t)\mathrm{e}^{2i\theta (𝐤,t_0,t)},\hfill \\ \dot{\beta }_𝐤^{}(t)\hfill & =& {\displaystyle \frac{eE(t)\epsilon _{}}{2\omega ^2(𝐤,t)}}\alpha _𝐤(t)\mathrm{e}^{2i\theta (𝐤,t_0,t)},\hfill \end{array}`$ (14)
where the dynamical phase is defined as

$$\theta (𝐤,t_0,t)=_{t_0}^t𝑑t^{}\omega (𝐤,t^{}).$$
(15)
To absorb the dynamical phase it is convenient to introduce new operators

$$\begin{array}{ccc}B_{𝐤r}(t)\hfill & =& b_{𝐤r}(t)e^{i\theta (𝐤,t_0,t)},\hfill \\ D_{𝐤r}(t)\hfill & =& d_{𝐤r}(t)e^{i\theta (𝐤,t_0,t)},\hfill \end{array}$$
(16)
which obey the anti-commutation relations:

$$\{B_{𝐤r}(t),B_{𝐤^{}r^{}}^+(t)\}=\{D_{𝐤r}(t),D_{𝐤^{}r^{}}^+(t)\}=\delta _{rr^{}}\delta _{\mathrm{𝐤𝐤}^{}}.$$
(17)
These operators satisfy the Heisenberg-type equations of motion

$`\begin{array}{ccc}{\displaystyle \frac{dB_{𝐤r}(t)}{dt}}={\displaystyle \frac{eE(t)\epsilon _{}}{2\omega ^2(𝐤,t)}}D_{𝐤r}^+(t)+i[H(t),B_{𝐤r}(t)],\hfill & & \\ {\displaystyle \frac{dD_{𝐤r}(t)}{dt}}={\displaystyle \frac{eE(t)\epsilon _{}}{2\omega ^2(𝐤,t)}}B_{𝐤r}^+(t)+i[H(t),D_{𝐤r}(t)],\hfill & & \end{array}`$ (20)
where $`H(t)`$ is the Hamiltonian of the system of quasiparticles

$$H(t)=\underset{r,𝐤}{}\omega (𝐤,t)\left[B_{𝐤r}^+(t)B_{𝐤r}(t)D_{𝐤r}(t)D_{𝐤r}^+(t)\right].$$
(21)
The first term in the r.h.s. of Eqs. (20) arises because of the unitary non-equivalence of the transition from the representation (8) to the quasiparticle one.

Next consider the evolution of the distribution function of quasiparticles with the momentum $`𝐤`$ and spin $`r`$ defined as

$$f_r(𝐤,t)=0_{in}|b_{𝐤r}^+(t)b_{𝐤r}(t)|0_{in}=0_{in}|B_{𝐤r}^+(t)B_{𝐤r}(t)|0_{in}.$$
(22)
According to the charge conservation the distribution functions for particles and anti-particles are related as $`f_r(𝐤,t)=\overline{f}_r(𝐤,t)`$. Taking derivative in Eq. (22) with respect to time $`t`$ we have

$$\frac{df_r(𝐤,t)}{dt}=\frac{eE(t)\epsilon _{}}{\omega ^2(𝐤,t)}\mathrm{Re}\{\mathrm{\Phi }_r(𝐤,t)\}.$$
(23)
Here the function $`\mathrm{\Phi }_r(𝐤,t)=0_{in}|D_{𝐤r}(t)B_{𝐤r}(t)|0_{in}`$ describes the vacuum production of pairs in external electric field $`E(t)`$. Applying the equations of motion (20), one finds

$$\frac{d\mathrm{\Phi }_r(𝐤,t)}{dt}=\frac{eE(t)\epsilon _{}}{2\omega ^2(𝐤,t)}\left[2f_r(𝐤,t)1\right]2i\omega (𝐤,t)\mathrm{\Phi }_r(𝐤,t).$$
(24)
The solution of Eq. (24) with the initial condition $`\mathrm{\Phi }_r(𝐤,t_0)=0`$ may be written in the following integral form

$$\mathrm{\Phi }_r(𝐤,t)=\frac{\epsilon _{}}{2}\underset{t_0}{\overset{t}{}}𝑑t^{}\frac{eE(t^{})}{\omega ^2(𝐤,t^{})}\left[2f_r(𝐤,t^{})1\right]\mathrm{exp}[2i\theta (𝐤,t^{},t)].$$
(25)
Inserting this result into the r.h.s. of Eq. (23) we obtain the anticipated kinetic equation

$$\frac{df_r(𝐤,t)}{dt}=\frac{eE(t)\epsilon _{}}{2\omega ^2(𝐤,t)}\underset{t_0}{\overset{t}{}}𝑑t^{}\frac{eE(t^{})\epsilon _{}}{\omega ^2(𝐤,t^{})}\left[12f_r(𝐤,t^{})\right]\mathrm{cos}[2\theta (𝐤,t^{},t)].$$
(26)
Since the distribution function does not depend on spin, the subscript $`r`$ can be dropped: $`f_rf`$. Substitution $`𝐩=𝐤e𝐀(t)`$, where the 3-momentum is decomposed onto the transverse and longitudinal components $`𝐩=𝐩(p_{},p_{}(t))`$, yields to the reduction of the KE (26) to Eq. (4) with the source term

$$S(𝐩,t)=\frac{e^2}{2}E(t)\mathrm{w}(𝐩)\underset{t_0}{\overset{t}{}}𝑑t_1E(t_1)\mathrm{w}(𝐩(t,t_1))\left[12f(𝐩(t,t_1),t_1)\right]\mathrm{cos}\left(2\underset{t_1}{\overset{t}{}}𝑑t_2\omega (t,t_2)\right),$$
(27)
where $`\mathrm{w}(𝐩)=\epsilon _{}/\omega (𝐩)`$ and

$`\omega (t,t_1)`$ $`=`$ $`\sqrt{\epsilon _{}^2+p_{}^2(t,t_1)},`$
$`𝐩(t,t_1)`$ $`=`$ $`𝐩e{\displaystyle _{t_1}^t}𝐄(t_2)𝑑t_2.`$ (28)
The source term (27) demonstrates several interesting features, such as the dependence on particle longitudinal and transverse momentum, the account for spin and statistics, and the non-Markovian character of the time evolution. The memory effects are caused by the time integration over the statistical factor $`(12f)`$ and the non-local cosine function smol1 , while the structure of the coefficient $`\mathrm{w}(𝐩)`$ defines the shape of the momentum distribution of created particles.

In the collisionless approximation the kinetic equation (4) with the source term given by Eq. (27) can be transformed into the system of three ordinary differential equations mam

$`2\dot{f}`$ $`=`$ $`eE\mathrm{w}v_1,`$
$`\dot{v}_1`$ $`=`$ $`eE\mathrm{w}(12f)2\omega v_2,`$ (29)
$`\dot{v}_2`$ $`=`$ $`2\omega v_1,`$
where the dot denotes the full time derivative, and the auxiliary functions $`v_1,v_2`$ defined as

$$v_1=2Re\mathrm{\Phi }_r(𝐤,t),v_2=2Im\mathrm{\Phi }_r(𝐤,t).$$
(30)
If the field strength is of order of the critical value, then it is necessary to take into account the back reaction of produced particles on primary field gatoff ; KES . The newly created particles polarize the vacuum and are accelerated by the external field. Their motion generates an internal field that in its turn modifies the initial background field. For the description of this phenomenon, the background field $`E(t)`$ in (27) should be replaced by the sum $`E(t)=E_{ex}(t)+E_{in}(t)`$, where the generated internal field $`E_{in}`$ can be found from the Maxwell equation smol

$$\dot{E}_{in}=\frac{eg_f}{(2\pi )^3}\frac{d^3p}{\omega }\left[2p_{}f+p_{}v_1\right],$$
(31)
with $`g_f=2N_c`$ being the degeneracy factor and $`N_c`$ the number of color degrees of freedom ($`N_c=3`$). The total current density in the r.h.s. of Eq. (31) is the sum of conductivity and vacuum polarization currents, respectively. The integrand in (31) contains the logarithmic divergence which should be removed somehow by means of a regularization procedure mam . We use here the simple ultraviolet cut-off of momentum integration on the border of a grid. The Eqs. (II) and (31) represent the closed system of equations for the numerical analysis of the back reaction problem.

## III Vacuum creation of quarks with different masses

The derived formalism can be applied to study the vacuum creation of quark-antiquark pairs in heavy ion collisions. Of special interest is the analysis of effects caused by the fast change of (chromo)electrical field at the time scale compatible with the inverse quark masses. To investigate a role of the field switching on/off effects, we approximate the time dependence of the flux-tube field by a short pulse oriented along the collision axis of primordial nuclei

$$E(t)=E_m\mathrm{exp}[(2t/\tau )^n],$$
(32)
where $`E_m`$ is the field magnitude and $`\tau `$ is the effective pulse width. The integer exponent $`n`$ governs the steepness of the pulse: for $`n1`$ (we restrict ourselves to $`n_{max}=6`$), the pulse (32) becomes close to the rectangular one. The case with $`n=2`$ reproduces well the soliton-like pulse

$$\begin{array}{ccc}E(t)\hfill & =& E_m\mathrm{cosh}^2(2t/\tau ),\hfill \\ A(t)\hfill & =& E_m(\tau /2)\mathrm{tanh}(2t/\tau ),\hfill \end{array}$$
(33)
for which the Dirac equation allows exact analytical solution grib ; nar .

It is convenient to define the time dependent strangeness suppression factor $`f_s`$ as the ratio of densities of strange $`(s)`$ to nonstrange $`(q)`$ quarks

$$f_s(t)=\frac{n(m_s,t)}{n(m_q,t)},$$
(34)
where in the case of axial symmetry

$$n(m_f,t)=g_f\frac{d^3p}{(2\pi )^3}f(m_f,p_{},p_{}).$$
(35)
The system (II) is integrated by the Runge-Kutta method with the zero initial conditions $`f(𝐩,t_0)=v_1(𝐩,t_0)=v_2(𝐩,t_0)=0`$. The momentum dependence of the distribution function $`f(p_{},p_{},t)`$ is determined by coarse-graining of the momentum space to a 2-dimensional grid; in each of its node the system of equations (II) is solved. The parameters of the grid depend on the field strength, the typical values are $`\mathrm{\Delta }p0.1m`$ (step of the grid) and $`p_{max}1520m`$ (boundary of the grid), so the total number of the solved equations is about $`10^6`$.

Effect of the field pulse duration. To study the influence of the field pulse duration on the ratio of strange and nonstrange quark pairs we compare results, obtained for the time-dependent ratio (34) with the soliton-like pulse (33), with those yielded for the constant chromoelectric field by Eq. (3). Input parameters, such as quark masses and string tension, are chosen according to amel , namely

$`m_s`$ $`=`$ $`350\mathrm{MeV},`$
$`m_{u,d}`$ $`=`$ $`230\mathrm{MeV},`$
$`\kappa `$ $`=`$ $`|eE|=0.9\mathrm{GeV}/\mathrm{fm}.`$ (36)
Since the critical field is defined as $`eE_{cr}=m^2`$, the chosen value of the field is under-critical one for heavy strange quarks and over-critical one for the light nonstrange quarks. To solve the Maxwell equation (31) the value of the effective charge $`e`$ should be determined. This can be done on basis of the estimates for the initial energy density $`\epsilon `$, which varies from 50 GeV/fm3 at RHIC to 520 GeV/fm3 at LHC cooper . Assuming for simplicity that all initial energy is deposited in the field sector, one obtains $`E_{ex}3.16`$ GeV/fm for RHIC and $`E_{ex}10`$ GeV/fm for LHC. The corresponding values of the effective charge are $`e0.22`$ and $`e0.07`$ for RHIC and LHC, respectively.

The results of calculations of the asymptotic value $`f_s(t\tau )`$ of the suppression factor $`f_s`$ are depicted in Fig. 1 for the field pulse (33). Trivially, the Schwinger-like estimate for the constant and infinite field (3) gives the constant value $`f_s=0.3`$. The time-dependent case is more complex: for pulses shorter than the particle inverse masses the creation of $`s\overline{s}`$ quark pairs is significantly enhanced. This is a direct consequence of the uncertainty relation for the energy and time. It is worth mentioning that already at $`\tau 3`$ fm/$`c`$ the exact result is very close to that given by the Schwinger formula (3). However, for the characteristic time of QGP formation $`\tau _01`$ fm/$`c`$ Bj83 , the production probability of strange quarks in the soliton-like field yield (33) is at least 1.5 times larger than that in the stationary case.

Influence of the field pulse shape. Figure 2 shows the time evolution of the factor $`f_s(t)`$ during action of field pulses (32) with $`n=6`$ and (33) with the same width $`\tau =1`$ fm/$`c`$. We see, that production of strange quarks increases with rise of the power $`n`$ in exponent in Eq.(32). Note also, that the intermediate effective value of the time dependent ratio $`f_s`$ is considerably larger than its final one: it varies from about 0.6 to 0.39 for the soliton-like pulse and from 0.75 to 0.52 for the rectangular one. The evolution of densities of the created particles is presented in Fig. 3. The density of both light and heavy quark pairs in the pulsing field initially increases. As the field saturates and starts to decrease, the process of particle absorption by the field dominates the particle production one. The longer the field pulses, the stronger absorption. Therefore, the way of the field oscillation can significantly change the value of the suppression factor. The similar result concerning the role of pulse shape of periodic laser field on the electron-positron vacuum production rate was obtained in pop  within the framework of the approximate imaginary time method.

In Ref. soff99  two ways of increasing the strangeness production within the framework of the Schwinger mechanism have been discussed, namely (i) either the field strength (string tension) is increasing, or, equivalently, (ii) quark masses are dropping due to chiral symmetry restoration. Assuming a significant reduction of quark masses from (III) to the current quark masses $`m_s=0.23`$ GeV, $`m_q=0.01`$ GeV soff99 , we get from Eq. (3) the enhancement of strangeness production in a constant infinite field from $`\gamma _s0.3`$ to $`\gamma _s0.4`$, which is equivalent to increase of the string tension $`\kappa `$ from 0.9 GeV/fm to 1.22 GeV/fm.

Let us perform similar fitting procedure for the rectangular pulse given by Eq. (32). The corresponding dynamical picture is displayed in Fig. 4. It appears that the scenarios considered above are not fully equivalent in dynamical sense, because the production of light particles at very short times is apparently suppressed. Therefore, the reduction of quark masses due to the chiral symmetry restoration is even more effective for the enhancement of strangeness production than the mechanism of string overlap, which leads to the formation of color rope and increase of the string tension.

The solution of Eqs.(II) and (31) shows that the back reactions play minor role in the pair production processes with the values of parameters given by the set (III), because the current of secondaries is quite weak. The increase of the string tension $`\kappa `$ from 0.9 GeV/fm to 1.22 GeV/fm or the reduction of quark masses to $`m_s=0.23`$ GeV and $`m_q=0.01`$ GeV does not significantly modify this scenario.

## IV Conclusions

Kinetic equation with the source term, describing the vacuum production of fermions in a strong, time-dependent field, is used for the two-component system investigation in the collisionless approximation. It is shown that dynamics of vacuum particle creation at time scale compatible with the inverse mass of particle depends essentially on the specific characteristics of field configuration (a form and a duration of the field pulse). In particular, the Schwinger-like regime of particle creation might not be realized at the typical time scale of QGP formation, $`\tau _01`$ fm/$`c`$. As a consequence, production of heavy strange quarks becomes more abundant. The role of string fusion and reduction of quark masses due to the chiral symmetry restoration which can alter the strangeness production is studied. Within the proposed dynamical scenario the mass reduction mechanism appears to be more effective for the enhancement of strangeness yield than the formation of color rope. Since the current of the produced secondaries is weak, the contribution of the back reactions to the production process is small. Our study shows that the time evolution picture of the gluon field should be incorporated consistently into the color flux tube model for quantitative description of parton production in ultrarelativistic heavy ion collisions.

Acknowledgments. Fruitful discussions with A. Faessler and C. Fuchs are greatfully acknowledged. This work was supported in part by the Russian Federations State Committee for Higher Education under grant E02-3.3-210, Russian Fund of Basic Research (RFBR) under grant 03-02-16877, and Bundesministerium für Bildung und Forschung (BMBF) under contract 06TÜ986.


text:
# Non-Poisson distribution of the time distances between two consecutive clusters of earthquakes

, , , , , ,

## Abstract

With the help of the Diffusion Entropy technique we show the non-Poisson statistics of the distances between consecutive Omori’s swarms of earthquakes. We give an analytical proof of the numerical results of an earlier paper \[Mega et al., Phys. Rev. Lett. 90 (2003) 188501\].

The model usually adopted to describe the time distribution of earthquakes is the Generalized Poisson (GP) model . The GP model assumes that the earthquakes are grouped into temporal clusters of events and these clusters are uncorrelated, and, therefore, completely unpredictable: The clusters are supposed to be distributed at random in time and therefore the time intervals between one cluster and the next one follow a Poisson distribution. The intra-cluster earthquakes are in fact correlated as it is expressed by the Omori’s law , an empirical law stating that the main shock, i.e. the highest magnitude earthquake of the cluster, occurring at time $`t_0`$ is followed by a swarm of triggered earthquakes (after shocks) whose number (or frequency) $`n(t)`$ decays in time as a power law, $`n(t)(tt_0)^p`$, with the exponent $`p`$ being very close to $`1`$. If we denote with the symbol $`\tau `$ the time intervals between one earthquake and the next, then right after a main shock, a short value of $`\tau `$ is followed with a large probability by another short value. For the same reason, far from a main shock and prior to the next one, a long value of $`\tau `$ is followed by another long value of $`\tau `$. This implies that the correlation function $`(\tau _i\tau )(\tau _j\tau )`$ is not zero for $`ij`$ and that it survives for all pairs of seismic events in between two consecutive unpredictable shocks. Omori’s law also implies  that the distribution of $`\tau `$, is a power law $`\psi (\tau )\tau ^p`$. This espression is valid in the time regime inside a swarm, and it is then truncated by a sharp cutoff caused by the arrival of the next swarm.

The catalog we have studied covers the period 1976-2002 in the region of Southern California spanning $`20^0`$ N -$`45^0`$ N latitude and 1000 W 1250 W longitude . This region is crossed by the most seismogenetic part of the San Andrea fault, which accommodates by displacement the primarily strike-slip motion between the North America and the Pacific plates, producing velocities up to $`47`$ mm/yr . The total number of recorded earthquakes in the catalog is $`383687`$.

Herein we disprove the GP model, providing evidence for the non-Poisson statistics of inter-cluster times, by applying to the mentioned catalog the Diffusion Entropy (DE) technique . Here we discuss with analytical arguments some issues that in an earlier paper  we have examined by means of a numerical treatment. As in the GP model, we assume that each cluster starts with an unpredictable triggering event (it may or may not be the main-shock). The distance between the $`i`$-th and the $`i+1`$-th cluster is therefore the time distance between such events, which we indicate as $`\tau _i^{[m]}`$, obeying the non-correlation property, $`\tau _i^{[m]}\tau _j^{[m]}=\tau ^{[m]}^2`$ if $`ij`$. The superscript $`m`$ stands for main-shock, but we actually need not to make the assumption that the triggering event is a large earthquake. As we shall see, the DE measures statistical properties of events with no need of identifying them.

Let us recall the definition of the DE functional $`S(t)`$ as the Shannon entropy of $`p(y,t)`$, the probability distribution to observe a fixed number of seismic events $`y`$ in a given time interval . This observation is equivalent to observing the spreading of a number of walkers making one step forward at each time where an event is met. Hence the term “Diffusion Entropy”. Different trajectories are chosen with the usual method of observing different time windows in the sequence, which is herein assumed to be stationary. Let us assume that $`p(y,t)`$ follows the scaling law with the form

$$p(y,t)=\frac{1}{t^\delta }F\left(\frac{y}{t^\delta }\right),$$
(1)
where $`\delta `$ is a positive exponent and $`F(x)`$ is a positive and integrable function of x. As a consequence of this assumption, after a straightforward algebra, we find that

$$S(t)=A+\delta \mathrm{ln}(t).$$
(2)
This means that the entropy of the diffusion process is a linear function of $`\mathrm{ln}(t)`$ and a measure of the slope is equivalent to the determination of the scaling parameter $`\delta `$.

We now show that, in the earthquake series under investigation, the principal source of entropy increase is given by the occurrence of the cluster-initiating seismic events. Let us indicate with $`\varphi (\tau ^{[m]})`$ the probability density function (pdf) of times between clusters and with $`h(x)`$ the pdf of the number of earthquakes in a cluster. The function $`h(x)`$, usually referred to as the Pareto’s law of earthquakes, is known to decay as $`h(x)1/x^{\alpha +1}`$, where $`\alpha `$ is a positive number. In the literature an exponent for clusters size distribution $`\alpha `$ ranging from $`\alpha =1.25`$  to $`3`$  is reported. On the other hand, we assume for $`\varphi (\tau ^{[m]})`$, which is exponential in the GP model, a form $`\varphi (\tau ^{[m]})(1/\tau ^{[m]})^\mu `$. We also assume that $`2\mu <3`$. As it will become clear in the next paragraphs, beyond the upper limit it is impossible for the DE to distinguish between an inverse-power law and an exponential; for $`\mu <2`$, on the other hand, the signal cannot be stationary . The connection between $`x`$ and the early defined variable $`y`$, over which the DE is calculated, is

$$y(t)=\underset{i=1}{\overset{z(t)}{}}x_i,$$
(3)
where the sum is carried over different clusters, $`x_i`$ is the number of shocks in the $`i`$-th cluster and $`z(t)`$ is the number of clusters in the same time window of length $`t`$ considered for $`y`$.

Let us call $`\widehat{\varphi }(s)`$ the Laplace transform of $`\varphi (\tau ^{[m]})`$ and $`\widehat{h}(k)`$ the Fourier transform of $`h(x)`$. Suppose now that the time duration of a cluster is negligible with respect to the mean time distance $`\tau ^{[m]}`$. In this case we can use directly a Continuous Time Random Walk (CTRW) formalism to calculate the probability $`p(y,t)`$ that a random walker, moving of a quantity $`x`$ at the time at which there is a cluster of size $`x`$ (and resting otherwise), is at position $`y`$ after a time $`t`$. For the Fourier-Laplace transform of $`p(y,t)`$ we will have, from the theory of CTRW ,

$$\widehat{p}(k,s)=\frac{1\widehat{\varphi }(s)}{s}\frac{1}{1\widehat{\varphi }(s)\widehat{h}(k)}.$$
(4)
To obtain the asymptotical behavior we write, if $`2<\mu <3`$,

$$\widehat{\varphi }(s)1\tau ^{[m]}s+cs^\gamma ,$$
(5)
where $`\gamma \mu 1`$. Eq.(4) becomes

$$\widehat{p}(k,s)\frac{\tau ^{[m]}}{1\widehat{h}(k)+s\tau ^{[m]}\widehat{h}(k)c\widehat{h}(k)s^\gamma }.$$
(6)
Considering that

$$\widehat{h}(k)1+ikx+bk^\alpha ,$$
(7)
if $`1<\alpha <2`$, we see that Eq. (6) leads to a ballistic scaling in the laboratory reference frame. The DE is insensitive to drifts, so we assume that Eq. (1) is fulfilled in a “detrended” moving reference frame, namely, where the position of the walker $`y`$ fulfills the condition $`y(t)=0`$ at each time $`t`$. To obtain $`\widehat{p}(k,s)`$ in such a moving reference frame we have to perform the following substitution

$$ss+ik\frac{x}{\tau ^{[m]}}.$$
(8)
After that, we apply the diffusive limit $`ks`$ obtaining

$$\widehat{p}(k,s)\frac{\tau ^{[m]}}{s\tau ^{[m]}bk^\alpha c\left(ix/\tau ^{[m]}\right)^\gamma k^\gamma },$$
(9)
thus proving that the most anomalous, namely, the smallest between the two exponents, either $`\alpha +1`$ or $`\mu `$, determines the asymptotic scaling according to the prescription $`\delta =1/\alpha `$ or $`\delta =1/(\mu 1)`$, respectively.

In Fig. 1 we report the results of the DE method. In full circles we plot the entropy $`S(t)`$ as a function of time when all the seismic events of the catalog are considered (independently of their magnitude $`M`$). A fit in the linear region gives a value of the scaling parameter $`\delta =0.94`$. We next consider (open symbols in Fig. 1) only the earthquakes with magnitude larger than a fixed value $`\overline{M}=2,3,4`$. We see that, regardless of the value of the threshold $`\overline{M}`$ adopted, the function $`S(t)`$ is characterized by the same long-time behavior with the same slope. This indicates that we are observing a property of the time location of the main earthquakes. Moreover, if the value of $`\delta `$ were due to $`h(x)`$ we should observe at most $`\delta =1/\alpha =1/1.25=0.80`$. This leads us to conclude that the asymptotic form (and scaling) of $`p(y,t)`$ is determined by the probability $`p(z,t)`$ of finding $`z`$ unpredictable events in a time window of duration $`t`$. This is in turn determined by a $`\varphi (\tau ^{[m]})`$ decaying as an inverse power law with exponent $`\mu =2.06`$. Using the theory of Ref.  we also determine the form of $`p(z,t)`$, which is well approximated by an asymmetric Lévy distribution with Lévy index $`\mu 1`$.

In conclusion, in this paper we have studied the statistical properties of earthquakes time distribution. Inter-cluster distances obey an inverse power law prescription $`\varphi (\tau ^{[m]})(\tau ^{[m]})^\mu `$ with $`\mu =2.06\pm 0.01`$ thus ruling out the GP model. The method proposed is based on the fact that the asymptotic properties of diffusion process generated by the seismic events are scarsely sensitive to the memory stemming from the Omori’s law. They are, on the contrary, sensitive to the anomalous statistics generated by the non-Poisson nature of the time distance between two consecutive large earthquakes. This non-Poisson behavior reflects, in our opinion, the cooperative behavior of the geological processes triggering the main shock, and consequently, to some extent, some sort of predictability. The emergenge of this possibility will be investigated in the future.

As a final remark, we recall that throughout our analysis we made the reasonable hypothesis that the series is stationary. On a formal ground, a scaling value $`\delta 1`$ is compatible with a nonstationary process. These nonstationary contributions may have different orgins. The most trivial is the lack of statistics: we may have for instance, a small number of Omori’s swarms, but very extended in time. For this reason we repeated the analysis using only portions of the series, obtaining the same results. Another, more interesting, source of non-stationarity could be a strong correlation among the series of the triggering processes, namely the condition $`\tau _i^{[m]}\tau _j^{[m]}\tau ^{[m]}^2`$. It is interesting to notice that also this condition might reflect a form of predictability to assess through a properly tailored form of statistical analysis of time series. We leave this as a subject of further investigation.

PG thankfully acknowledges the Army Research Office for financial support through Grant DAAD19-02-0037.