text:
# I Introduction

## I Introduction

Reacting complex flows are ubiquitous in the atmosphere, the oceans, and the interior of our planet: Stratospheric ozone chemistry is perhaps the most well known example of chemical dynamics occurring in a stirred geophysical fluid. But the fate of many other reacting tracers is also influenced by mesoscale eddies in the ocean, or by convection in the Earth mantle, to name just a few cases. Biological population dynamics is also a kind of “chemical reaction”. The interactions of nutrients and plankton species in the sea are the first steps in the ocean food chain, and they are among the most important ingredients for the understanding the CO2 interchange between atmosphere and ocean.

Tracers stirred by fluid motion are known to develop strong inhomogeneities, usually in the form of filamental features, arising from a kind of variance cascade from the forcing scale towards smaller scales. These structures are poorly resolved in global atmospheric and, especially, in oceanic models. They provide however sensitive mixing mechanisms and they are in some sense “catalists”  that enhance the chemical or biological activity occurring in geophysical flows. We mention as examples the increase of ozone depletion produced by the filamental structure of chlorine filaments , or the increase in biological production in small-scale ocean structures .

In this set of Lectures, we will present some theoretical ideas for the study of structures arising in advected reacting tracers. The emphasis will be in the modifications that chemical or biological activity introduces in reactive patterns with respect to the ones appearing in passively advected substances.

Examples from plankton populations dynamics will be used through these Lectures. It should be mentioned, however, that if one looks into the literature, the diversity of observations related to plankton distributions is so great that some authors, rather pessimistically, argue that “… One should expect no general results that describe in all (or most) cases how ’the biology’ modifies a spatial pattern that has arisen solely from the dispersal of species …. We use here the biological experimental observations more as a motivation and illustration, than as a detailed check of a particular theory. Our aim is to point out simple and robust mechanisms that can be useful in classifying the large number of observations into a more reduced set of different scenarios. Later quantitative test would need synoptic measurements of the biological variables and of the hydrodynamic flows. This would be certainly more accurately done in laboratory chemical experiments than in open sea biology. Most of the results presented here are of direct application to chemical reactions of the decaying or of the excitable type. Nevertheless, we mention that comparison of some of our theoretical ideas with simultaneous satellite data of both a reactive tracer (sea temperature in interaction with air temperature) and hydrodynamic flow has been very recently initiated .

In Section II we will briefly review some aspects of plankton dynamics that will be useful in the following. The paradigm of chaotic advection, i.e. the chaotic transport of fluid parcels by smooth large-scale flows, will be the one more frequently used here to discuss flow effects, and it is described in Section III. Sections IV and V describe results obtained for two rather general classes of chemical or biological reactions in chaotic flows: the first-order reaction, or linearly decaying tracer, and reactions of the excitable type. Finally (Sect. VI), some warnings will be given about the difficulties in modeling discrete individuals (such as planktonic organisms) in terms of continuous concentration fields.

## II Plankton, patches, and blooms

### A Background

Plankton is the generic name given to a huge variety of aquatic organisms, both from marine and from freshwater environments, comprising from viruses with less than 200 nm in size, to crustaceans or even fish larvae, with sizes above the centimeter. What they have in common, and distinguishes them from other major group called nekton, is their inability to overcome the major ambient hydrodynamic currents. In this sense, they are strongly influenced by the existent hydrodynamic weather. It is often assumed that they are passively transported by the marine currents. This is probably true for the smallest species, but it should be said that most of these organisms present flagellae and other natatory organs that give to them some mobility, especially in the vertical direction. In addition to interactions with the hydrodynamic environment different plankton species interact among them in a variety of ways, such as predator-prey interaction, parasite-host relationship, competition for resources, etc. The biological interaction between species and with the nutrient substances dissolved in water is formally equivalent to a chemical reaction dynamics.

Any attempt to overview here some major theme in plankton biology and its interaction with hydrodynamics would result very partial and incomplete, and we refer the interested reader to some monographs that turn out to be very readable. We just mention that a major division between plankton species arises between the ones that are able to assimilate inorganic material by photosynthesis, the phytoplankton, and the species that necessarily graze on the above, the zooplankton. In the following we will briefly describe just two phenomena, plankton patchiness and plankton blooms, extremely stimulating from the perspective of cross-disciplinary research, and for which some ideas from Nonlinear Science may be of relevance.

Plankton patchiness  is the term used to label the large inhomogeneity observed in plankton distributions. Earlier observations from land and from ships reported water masses, of sizes around say 10 km, distinctly colored by the presence of different species of planktonic organisms. The search for biological or physical mechanisms singling out a characteristic patch size, despite mixing with the surrounding water stimulated the first theoretical approaches. Modern satellite observations reveal that the range of sizes of plankton structures is much broader. The largest features, thousand of kilometers wide, are associated with major ocean gyres and currents. Upwellings of nutrient-rich water from deeper ocean layers, and river discharges, fertilize the upper layer of the ocean and induce associated plankton structures. In the mesoscale range, fronts and eddies are seen as major players, both in bringing to the surface nutrients from lower layers, and in providing horizontal transport and stirring processes.

A closer look at the inhomogeneities in plankton distributions can be performed by analyzing water samples along a ship trajectory (or transect). Higher spatial resolution can then be attained and the result is an extremely intermittent distribution (Fig. 1, lower panel). At this point, it is pertinent to ask about the relevance of the fact that plankton is a living substance, or if perhaps the same kind of distributions can be found for non-biological tracers. A partial answer to this question is given by Fig. 1, in which distributions of physically controlled tracers (temperature and salinity) are compared with the phytoplankton concentration, indicated by water induced fluorescence, at a resolution of about half a kilometer. Although correlations can be appreciated, evident differences are seen between the distributions of both kinds of tracers. A brief walk on the literature shows that the distributions may look very different form place to place in the ocean, and that their characteristics may also change in time. In general one can say that they will be affected by both the characteristics of the fluid flow that transports them, and by the chemical or biological interactions in which they are involved.

The irregular nature of the observed distributions immediately suggest characterization in terms of scaling exponents. This has traditionally been done in terms of the power spectrum of concentration fluctuations, since a comparison with Kolmogorov-like turbulence theories was then direct. In situations in which the phytoplankton concentration spectrum displayed a large and well defined power-law regime (which is not the case of the data in Fig. 1) it decayed as $`k^\beta `$, with $`k`$ the wavenumber. The scaling exponent $`\beta `$ results to be larger than one. We prefer to discuss the scaling properties of concentration fluctuations in terms of structure functions, since a more complete characterization of intermittency properties can be done. These quantities, and their relation with the power spectrum, are introduced in Sect. IV.

A second kind of inhomogeneity found in plankton distributions is temporal inhomogeneity associated to plankton blooms. Seasonal blooms are episodes of rapid phytoplankton growth that occur in mid-latitude seas at the beginning of spring. This is followed by zooplankton (and also bacterial) growth that consumes the available phytoplankton until concentrations return back to the initial low level. The whole cycle finishes with the summer. Sometimes a smaller version of the bloom repeats in autumn. The mechanism for seasonal blooms seems to be understood, at least qualitatively: the beginning of ocean stratification in spring reduces the depth of the layer in which phytoplankton is mixed by turbulent motion, the mixed layer, so that these photosynthesis-capable organisms spend more time exposed to solar light. Fast growth then occurs. The subsequent growth of predators, and the consumption of nutrients in the upper layer, ends the bloom. Destratification, with the destruction of the seasonal thermocline in autumn, brings again to the surface some nutrients from the deeper layers, so that phytoplankton activity may increase again.

There is another kind of blooms that occur more localized in space, and not with a seasonal periodicity, but occasionally and associated to local warming or local increase of nutrients in the water. They are called sometimes red tides since the large amount of organisms in water may give to it some distinct coloration. Time scales for the initial development of the bloom are in the range of days, and the end usually arrives in a time of the order of weeks or month.

### B Models

Modeling the above spatial and temporal evolutions of planktonic populations is a challenge since many physical and biological processes are involved, some of them poorly understood. A variety of levels of description are available, from detailed models based on individual organism behavior, to mean field approaches that consider the whole sea as an homogeneous soup without spatial features. Another aspect of modeling that should be addressed is the degree of aggregation at which groups of species are represented by the same collective variables. It should be said that these starting points, the election of the adequate levels of description, and the relationship between different degrees of detail, is an important subject by itself, and its complexity is just beginning to be grasped. More on this in Sect. VI.

A class of models simple enough to allow some explicit analysis, but still retaining some of the complexities of real-life trophic interactions is the one known as the NPZ class of models (from Nutrients, Phytoplankton and Zooplankton). We describe it here for illustration, and because it will be of use in some of the results presented later. The processes included within this framework are displayed in Fig. 2, and its mathematical modeling is as follows: in each portion of well mixed fluid, the time evolution of the amount of nutrients ($`N`$), phytoplankton ($`P`$), and zooplankton ($`Z`$) is ruled by

$`{\displaystyle \frac{dN}{dt}}`$ $`=`$ $`F_N(N,P,Z)=s(N_0N)f(N)P+cP+ag(P)Z+eZ`$ (1)
$`{\displaystyle \frac{dP}{dt}}`$ $`=`$ $`F_P(N,P,Z)=f(N)PqPg(P)Z`$ (2)
$`{\displaystyle \frac{dZ}{dt}}`$ $`=`$ $`F_Z(N,P,Z)=kg(P)ZmZ^r.`$ (3)
The nutrient input $`N_0`$ affects the nutrient concentration in the mixed region, and this influence is transferred to phytoplankton by the second equation, and to zooplankton by the third. Finally, zooplankton either dies or is consumed by higher organisms, which leads to the last term in the third equation. The exponent $`r`$ in this last term is usually taken to be either $`r=1`$ or $`2`$, depending on the dominant kind of predation on zooplankton one wants to model.

For the description of the spatial dependence, the framework of Advection-Reaction-Diffusion equations is frequently used. They are partial differential equations in which the organisms and their nutrients are described by continuous distributions subjected to advection by a velocity field $`𝐯(𝐱,t)`$, and to diffusion, in addition to the interactions (3). The general form of the equations is

$$\frac{C_i}{t}+\left(𝐯C_i\right)=F_i(N,P,Z)+D_i^2C_i.$$
(4)
The index $`i`$ takes three values, so that $`F_i`$ is either $`F_N`$, $`F_P`$, or $`F_Z`$, and the variables $`C_i`$ are either $`N`$, $`P`$, or $`Z`$ (and now these symbols denote concentrations per unit volume). If the flow is incompressible ($`𝐯=0`$) then the second term in the left-hand-side can be written as $`𝐯C_i`$. At scales large enough, so that the effects of diffusion can be neglected, Eq. (4) admits a Lagrangian representation that in the incompressible-flow case coincides with the system in Eq. (3), where now $`N`$, $`P`$, and $`Z`$ are either the concentration or the amount of substance contained in a fluid element that moves with the flow velocity. The trajectory $`𝐫(t)`$ of this fluid parcel thus satisfies

$$\frac{d𝐫}{dt}=𝐯(𝐫,t).$$
(5)
It may seem that Eqs. (3) and (5) represent uncoupled dynamical systems. This would be the case if none of the parameters in (3) is spatially dependent. But if one of them varies from point to point, it should be evaluated at the position of the fluid element at time $`t`$, that is, at $`𝐫(t)`$. This couples the chemical and the transport dynamical systems. For example, the parameter $`N_0`$ may represent an inhomogeneous nutrient input $`N_0(𝐱)`$ into the upper ocean mixed layer by a localized upwelling, and thus the fluid element following the trajectory $`𝐫(t)`$ would experience a time-dependent input $`N_0(𝐱=𝐫(t))`$.

The model is still largely undefined until one specifies the parameters and, particularly, the response functions $`f`$ and $`g`$, which really contain information about the species interaction. A rather common form for the function $`f`$, experimentally established by Monod, is the so-called Michaelis-Menten form

$$f(N)=\mu \frac{N}{N_e+N}$$
(6)
which may be also used for the grazing function ($`g(P)=\alpha P/(P_e+P)`$). This form may be justified by arguments taken from enzyme kinetics . Any function that, like (6), presents a linear increase close to the origin followed by saturation at large $`N`$ receives the name of Hollings type-II response function. Another commonly used Hollings type-II response function is the Ivlev one: $`g(P)=\alpha (1e^{\lambda P})`$. Hollings type-III response functions differ from the above by its behavior at small argument: the linear increase is replaced by a slower one, indicating that predators will not be interested in a too small prey concentration. An example of this type is

$$g(P)=\alpha \frac{P^2}{P_e^2+P^2}.$$
(7)
Large uncertainties exist in the form of the response functions and in the numerical values of the parameters involved in (3) for particular ecosystems. We mention however that quantitative determination can be achieved under laboratory conditions .

Many variations and simplifications of model (3) can be found in the literature. It is rather common to reduce the modeling to the predator-prey competition of zooplankton and phytoplankton. In this case the limitation imposed on the phytoplankton growth by the nutrient concentration is included by replacing the second equation in (3) by one containing a logistic term:

$$F_P(P,Z)=\alpha P\left(1\frac{P}{C}\right)+\mathrm{}$$
(8)
The parameter $`C`$ is called the carrying capacity, and represents the maximum amount of phytoplankton the fluid parcel can support.

### C Some proposed mechanisms for patchiness

Probably, among the first theoretical approaches used to address the formation of plankton patches is the one contained in papers by Skellam, and by Kierstead and Slobodkin (the so-called KISS theory). They considered the free growth (at rate $`\mu `$) of phytoplankton (with concentration $`P(𝐱,t)`$), in a finite region of water suitable for growth, with plankton escaping the suitable zone by diffusion:

$$\frac{P}{t}=\mu P+D^2P.$$
(9)
This is a very particular case of (4). It is assumed that $`P=0`$ outside the region. If the patch is too small, the growth can not overcome the diffusive escape flux, and the patch dies. This can be seen, for example, via the exact one-dimensional solution of (9), with boundary conditions such that phytoplankton disappears outside the region $`[0,L]`$ of suitable conditions for growth (i.e. $`P(x=0,t)=P(x=L,t)=0`$):

$$P(x,t)=\underset{n=1}{\overset{\mathrm{}}{}}A_n\mathrm{sin}\left(\frac{n\pi x}{L}\right)e^{t\left(\mu D(\frac{n\pi }{L})^2\right)}.$$
(10)
The constants $`A_n`$ are fixed by the initial conditions. It is clear that if $`L<L_c=\pi \sqrt{D/\mu }`$ no growth occurs and the patch dies, so that $`L_c`$ can be identified with the minimum patch-size able to support growth.

There is a number of criticisms than can be posed to this approach. The first one is that the model is incomplete, since once growth begins it continues without limit. Nonlinear saturation and interactions with predators would be needed to stop this. The diffusion coefficient $`D`$ is certainly not the originated from the Brownian motion of the organic particles, since this would be irrelevant to processes above, say, the millimeter scale. It is rather a turbulent eddy-diffusion coefficient aimed to represent in an averaged way the effect of dispersion by the turbulent flow, with values that would depend on the observation scale as found for example in . One should ask, thus, how can the plankton disperse in such way, whereas the nutrients remain static within the fixed parcel, despite they live in the same flow. In any case, if one introduces a growth rate for phytoplankton of the order of one day ($`\mu 10^5`$ s-1) and an eddy diffusivity appropriate for the 10-100 km range ($`D100`$ m2/s) , one finds $`L_c10`$ km, which seemed to compare well with the earliest observations of plankton patches.

An interesting development that has some common ideas with the KISS approach, but introduces a new and important element, is described in . The idea is that what will control dispersion of a patch is not only eddy diffusivity, but also the geometric characteristics of the mean flow. It turns out that if an incompressible fluid flow induces dispersion in one direction it necessarily produces concentration in another, to conserve the fluid volume. The situation is especially clear in a purely two-dimensional flow, where the rate of expansion of a small fluid blob in one direction should be exactly the same as the rate of contraction in another. Thus, an initially circular patch will take the form of a filament locally aligned with particular directions in the flow, the unstable manifolds of the flow hyperbolic points. These manifolds have been nicely observed in laboratory experiments , and phytoplankton filaments on the sea surface have been observed from space . Reference  proposes the following model to describe the transverse profile of a phytoplankton filament:

$$\frac{P}{t}\lambda x\frac{P}{x}=\mu (t)P+D\frac{^2P}{x^2}.$$
(11)
There are two differences with the KISS model. One is the possible time-dependence of the growth rate $`\mu (t)`$, that is a simplified linear way of modeling nonlinear interactions and predation on phytoplankton after the first growth stages. The other is the advective term $`\lambda x_xP`$, associated with a velocity field $`v_x(x)=\lambda x`$ that models a local strain that, in the direction perpendicular to the filament, compresses it. A solution to Eq. (11) that is stable and attractive can be found analytically : a Gaussian with time-dependent height but fixed width $`w=\sqrt{D/\lambda }`$. Thus, the lateral scale of the filament is not controlled by some externally imposed size of the suitable water, but by the competition between diffusion and advection. It is somehow surprising that the biological growth rate does not affect the filament width in this model. As we will see in Section V this feature is not shared by other models.

A different mechanism that has been invoked as possible source of plankton patchiness  is the Turing mechanism. It appears generically when there is competition between a self-reproducing prey and a predator that diffuses faster than the prey, and leads to patterns with a characteristic periodicity. Although diffusion coefficients of the organic particles have roughly the same very small value for both predator (zooplankton) and prey (phytoplankton), the former has usually a larger vertical mobility. If there is some vertical shear in the water column, the organisms with larger vertical mobility will experience a more variable velocity field, and thus they will be subjected to larger dispersion. If one models such enhanced dispersion by an effective diffusion coefficient, it would be larger for zooplankton than for phytoplankton, and then the Turing mechanism can be at work. The analysis in  implies that for models of the type (3), the Turing mechanism does not occur in the relevant range of parameters. But this may not be the case for modifications of this class . In addition, related mechanisms, that would lead also to plankton distributions spatially periodic, may appear in the presence of differential motion between predator and prey .

Horizontal stirring is another mechanism that produces strong inhomogeneity. This is the one that will be addressed in the following. According to , the attractors of models of the form (3) are, for most of the relevant parameter range, fixed points representing stable species coexistence. In most cases, at least if type III Hollings functions are used, excitable behavior  occurs close to these equilibrium points. This motivates us to consider the effect of stirring on population models displaying either simple relaxation to fixed points (Section IV) or excitable behavior (Section V). Population oscillations also occur at realistic parameter values , but we will not consider this situation here. Before going to the results for the chemical/biological dynamics, we summarize in the next section some results about chaotic advection, which is the kind of fluid flow considered in the developments presented here.

## III The paradigm of chaotic advection

When dealing with transport processes in complex fluid flows, the concept of turbulence, the unsteady and irregular kind of flow in which motion occurs in a large range of scales, is the paradigm that most often comes to mind. Fluid elements, and transported particles, follow intricate trajectories in such velocity fields. It was recognized some time ago  that a turbulent velocity field is not necessary to produce chaotic fluid-element trajectories. In three spatial dimensions the dynamical system (5) has chaotic trajectories  even for very simple steady and smooth velocity fields $`𝐯(𝐱)`$. In two-dimensional flows, a simple periodic time dependence in $`𝐯(𝐱,t)`$ is enough to induce generically chaotic trajectories in (5). This situation of chaotic trajectories produced by a simple and regular velocity field is called chaotic advection or Lagrangian turbulence. The essential characteristic of chaos is sensibility to initial conditions. This means that fluid elements initially very close would follow diverging trajectories. The exponential rate of separation is given by the flow Lyapunov exponent $`\lambda `$. In an incompressible fluid, contraction must occur in another direction to conserve volume. The clearest situation is in hyperbolic two-dimensional flows, where one can define at each point an expanding direction and a contracting direction, with both the contraction and expansion rate occurring at the same exponential rate $`\lambda `$. As a consequence, portions of fluid initially compact become stretched in long and thin filaments, and are also repeatedly folded to remain inside the system. As a result, chaotic advection produces a “cascade” of inhomogeneities from large scales to smaller and smaller scales, in the same qualitative way as a fully turbulent velocity field will do. Quantitative details are, however, different. It turns out that the patterns produced by advection in smooth velocity fields are more singular than the ones produced by singular velocity fields such as Kolmogorov turbulence.

Experimental realizations of (nonturbulent) chaotic advection can be performed in small containers with high viscosity fluids, so that Reynolds number remains small and the flow laminar. It is now recognized that the regime of chaotic advection also applies to turbulent flow at scales smaller than the Kolmogorov scale ( $`l_K(\nu ^3/ϵ)^{1/4}`$, where $`\nu `$ is the kinematic viscosity and $`ϵ`$ the energy dissipation rate). This is the so-called Batchelor regime. The Lyapunov exponent in that range is of the order of the mean strain $`\lambda \sqrt{ϵ/\nu }`$. For water ($`\nu 10^6`$ m2s-1) at conditions typical in the upper mixed layers of oceans and lakes ($`ϵ10^7`$ W kg-1,), $`l_K`$ is between the millimeter and the centimeter. Thus the Batchelor regime is certainly relevant for plankton processes such as feeding, aggregation, encounter, etc. It is however a range much smaller than the structures shown in Fig. 1 or the ones visible from satellites. A third situation where chaotic advection has been used with success is in the modeling of the largest scales of geophysical flows . Smaller motions would then have to be represented by some eddy diffusion.

In the following, we focus on the situation of chaotic advection. This would allow us to progress in our objective of identifying mechanisms that can contribute to classify observations in robust scenarios. Chaotic advection is a framework simple enough to obtain a number of explicit results that turn out to be rather insensible to model details. Of course, explanation of particular geophysical observations would require additional discussion on whether or not the paradigm of chaotic advection applies to the scales considered.

## IV The decaying transported substance: a model for stable dynamics

Most of the simpler models of chemical dynamics, and models such as (3) in a large range of parameters, evolve in such a way that, at long times, the different species reach an equilibrium at which they coexist at some fixed-point concentration values. When some of the coefficients in Eq. (3) are functions of space they act (when coupled with the transport equation (5)) as time-dependent source or sink forcing terms that disturb such equilibrium. Even in this forced case it may happen that the concentration values at each fluid element tend to relax to a unique time-history, determined by the fluid-element trajectory. When this happens we say that we have a stable chemical dynamics. In a sense, the concentrations try to relax to the local-equilibrium values that would correspond to the values of the parameters found along the fluid-element trajectory. The mathematical way to check this stability is by initializing a particular fluid element with two slightly different sets of concentration values, and monitor how these values evolve under the coupled chemical and transport dynamics (3) and (5). The exponential rate of divergence between these two different initial concentrations, under the same fluid trajectory, defines the so-called chemical Lyapunov exponent, $`\lambda _C`$. If this quantity is negative, it identifies convergence to a common evolution, and this is the case of stable chemical dynamics. The usual ergodic properties would guarantee that the same value of the chemical Lyapunov exponent will be obtained for almost all fluid trajectories and initial concentrations.

The simplest case of stable chemical dynamics is the linearly decaying passive scalar, or first order reaction. It consists in the evolution of a single concentration $`C(𝐱,t)`$ under a linear decay at rate $`b`$:

$$\frac{C}{t}+𝐯(𝐱,t)C=S(𝐱)bC+D^2C,$$
(12)
A space-dependent forcing consisting in a source $`S(𝐱)`$ of the substance is included, to maintain a non-trivial concentration field at long times.

It is easy to check (by using the Lagrangian representation valid for $`D0`$) that $`\lambda _C=b`$. The simple dynamics (12) can be considered either as an approximation to more complex chemical or biological evolutions with stable chemical dynamics, or as a description of simple specific processes such as spontaneous decomposition of unstable radicals, evolution of a radioactively or photochemically decaying substance, or relaxation of sea-surface temperature towards atmospheric values . It also describes the concentration of the $`C`$ reactive in the binary reaction $`B+CD`$ when the $`B`$ concentration is maintained constant.

The first systematic study of model (12) seems to be in . The author analyzed, by a spectral element method, the power spectrum in a number of threedimensional turbulent regimes. In the Batchelor regime a power-law is found with an exponent depending on the quotient between the decay rate $`b`$ and the mean strain rate, which is equivalent to a flow Lyapunov exponent $`\lambda `$. Here we will obtain an essentially equivalent result within the framework of chaotic advection. In addition, intermittency corrections and the differences between open and closed flows will be discussed. Our results are first obtained in terms of the behavior of the concentration differences between spatially close points ($`\delta C=C(𝐱+\delta 𝐱,t)C(𝐱,t)`$) and later in terms of the structure functions, that are the moments of $`\delta C`$ averaged in space.

Exact results in  for the decaying scalar in a closed random Kraichnan flow in the Batchelor regime confirm the generality of the mechanism found: chaotic advection and decaying chemistry produce singular patterns that can be characterized by scaling exponents depending on the chemical decay rate and on the flow Lyapunov exponent (more precisely, on the finite-time Lyapunov exponent distribution). We will present here also results for nonlinear plankton models, showing that the results of the linearly decaying chemistry can be translated to the case of nonlinear stable chemistry at the smaller scales simply by replacing the decay rate $`b`$ by the absolute value of the chemical Lyapunov exponent $`|\lambda _C|`$. Detailed justification of this can be found in  for multiple species in the framework of chaotic advection, or in  for a single concentration in the framework of the random Kraichnan flow. It is important to mention that in the present framework the relevant time scale affecting the distributions structure is identified to be $`\lambda _C`$, and not other time scales frequently used in the literature such as the phytoplankton linear growth rate.

### A The smooth-filamental transition for the decaying passive scalar in a closed flow

We consider (12) with a two-dimensional, incompressible, smooth, and non-turbulent velocity field. Chaotic advection is obtained generically if a simple time-dependence is included in $`𝐯(𝐱,t)`$. We assume that diffusion is weak and transport is dominated by advection. Thus one expects that the distribution on scales larger than a certain diffusive scale is not affected by diffusion, and set $`D=0`$. In this limit the above problem can be described in a Lagrangian picture by the pair of dynamical systems (5) and

$$\frac{dC}{dt}=S[𝐱=𝐫(t)]bC,$$
(13)
where the solution of (5) gives the trajectory of a fluid parcel, $`𝐫(t)`$, while (13) describes the Lagrangian chemical dynamics in this fluid element: $`C(t)C(𝐱=𝐫(t),t)`$.

To obtain the value of the chemical field at a selected point $`𝐱`$ at time $`\overline{t}`$ one needs to know the previous history of this fluid element, that is the trajectory $`𝐫(t)`$ ($`0t\overline{t}`$) with the final condition $`𝐫(\overline{t})=𝐱`$. This can be obtained by the integration of (5) backwards in time. Once $`𝐫(t)`$ has been obtained, the solution of (13) is

$$C(𝐱,\overline{t})=C[𝐫(0),0]e^{b\overline{t}}+_0^{\overline{t}}S[𝐫(t)]e^{b(\overline{t}t)}𝑑t.$$
(14)
Then one can obtain the difference at time $`\overline{t}`$ of the values of the chemical field at two different points $`𝐱`$ and $`𝐱+\delta 𝐱`$ separated by a small distance $`\delta 𝐱`$ in terms of the difference $`\delta C[𝐫(t),t;\delta 𝐫(t)]C[𝐫(t)+\delta 𝐫(t),t]C[𝐫(t),t]`$ for $`0t\overline{t}`$, namely:

$$\delta C(𝐫,\overline{t};\delta 𝐱)=\delta C[𝐫(0),0;\delta 𝐫(0)]e^{b\overline{t}}+_0^{\overline{t}}\delta S[𝐫(t);\delta 𝐫(t)]e^{b(\overline{t}t)}𝑑t$$
(15)
where $`\delta 𝐫(t)`$ ($`0t\overline{t}`$) is the time-dependent distance between the two trajectories that end at $`𝐱`$ and $`𝐱+\delta 𝐱`$ at time $`\overline{t}`$, and $`\delta S`$ is the difference of the source term at points $`𝐫(t)`$ and $`𝐫(t)+\delta 𝐫(t)`$. In the following we omit the first term in the right-hand-side of (15) since it always disappears at long times.

If the trajectory is chaotic, we have in the backwards dynamics, for $`t<0`$ and large, $`|\delta 𝐫(t)|\left|\delta 𝐫(0)\right|e^{\lambda t}`$, where $`\lambda `$ is the positive value of the Lyapunov exponent along the trajectory. In (15), $`\delta 𝐫`$ is obtained backwards starting from $`\delta 𝐱`$ at $`t=\overline{t}`$. In this case:

$$\delta 𝐫(t)\delta 𝐱e^{\lambda (\overline{t}t)},t<\overline{t}.$$
(16)
We note that there is a particular direction for the orientation of the initial displacement $`\delta 𝐱`$ along which the trajectories converge instead of diverging (this is the contracting direction in the backwards flow which is the expanding one in the forward dynamics). For $`\delta 𝐱`$ oriented along it, $`\lambda `$ in (16) should be replaced by $`\lambda `$. For hyperbolic dynamical systems, the value of the Lyapunov exponent is the same for almost all the initial or final conditions $`𝐱`$ . We will show later, however, that the deviations that could occur in sets of zero measure may have some observable consequences.

Taking the limit $`\delta 𝐱0`$ and with substitution of (16) in Eq. (15), one obtains

$$\delta C(𝐫,\overline{t};\delta 𝐱)\delta 𝐱_0^{\overline{t}}S[𝐫(t)]e^{(\lambda b)(\overline{t}t)}𝑑t.$$
(17)
Rewriting $`\delta 𝐱=\overline{𝐧}|\delta 𝐱|`$ so that $`\overline{𝐧}`$ is a unit vector, one finds the directional derivative along the direction of $`\overline{𝐧}`$ as

$$\overline{𝐧}C(𝐫,\overline{t})\overline{𝐧}_0^{\overline{t}}S[𝐫(t)]e^{(\lambda b)(\overline{t}t)}𝑑t.$$
(18)
If $`\lambda <b`$ this derivative remains finite in the $`\overline{t}\mathrm{}`$ limit and the asymptotic field $`C_{\mathrm{}}(𝐱)C(𝐱,\overline{t}\mathrm{})`$ is smooth (differentiable). Otherwise the derivatives of $`C`$ diverge as $`e^{(\lambda b)\overline{t}}`$ leading to a nowhere-differentiable irregular asymptotic field. The exception occurs when $`\delta 𝐱`$ points along the expanding direction of the forward flow: $`\lambda `$ should be replaced by $`\lambda `$ and the directional derivative is always finite. The irregular field, thus, has a filamental character. It should be noted that the limiting distribution $`C_{\mathrm{}}`$ is not a steady field, but one following the time dependence of the flow. For time-periodic flows $`𝐯(𝐱,t)`$, $`C_{\mathrm{}}`$ will also be time periodic. Its singular characteristics however do not change in time.

The characterization of the singular asymptotic field has been performed in . Defining $`\delta C_{\mathrm{}}\delta C(𝐱,\overline{t}\mathrm{})`$ and taking into account that the exponential separation (16) can only persist until it reaches the characteristic scale of the flow field, one finds  in the $`|\delta 𝐱|0`$ limit the scaling

$$\delta C_{\mathrm{}}(𝐫;\delta 𝐱)|\delta 𝐱|^\alpha ,$$
(19)
with a Hölder exponent $`\alpha `$ given by

$$\alpha =\mathrm{min}\{\frac{b}{\lambda },1\}.$$
(20)
This general relation expresses the local Hölder exponent in terms of the local infinite-time Lyapunov exponent and the chemical decay rate. As stated before, for hyperbolic systems this Lyapunov exponent has the value $`\lambda _0`$ everywhere but in a set of zero measure. By decreasing $`b`$ or increasing $`\lambda `$ one finds a transition from a smooth $`\alpha =1`$ distribution to a rough $`\alpha <1`$ one with filamental structure. This is the so-called smooth-filamental transition. For uniform $`\lambda =\lambda _0`$, the relationship between the Hölder exponent and the power spectrum decay exponent ($`\beta `$ in the power spectrum scaling $`S(k)k^\beta `$) is $`\beta =1+2\alpha `$. Thus the result (20) implies an exponent $`\beta `$ larger than one, as it is usually observed. The Batchelor result  $`S(k)k^1`$, valid for the nondecaying scalar in smooth velocity fields is recovered for $`b=0`$.

In Fig. 3 we present snapshots of the asymptotic field $`C_{\mathrm{}}`$ evolving according to (5) and (13). For the flow we take a simple time-periodic velocity field defined in the unit square with periodic boundary conditions by

$`v_x(x,y,t)`$ $`=`$ $`{\displaystyle \frac{2U}{T}}\mathrm{\Theta }\left({\displaystyle \frac{T}{2}}tmodT\right)\mathrm{cos}(2\pi y)`$ (21)
$`v_y(x,y,t)`$ $`=`$ $`{\displaystyle \frac{2U}{T}}\mathrm{\Theta }\left(tmodT{\displaystyle \frac{T}{2}}\right)\mathrm{cos}(2\pi x)`$ (22)
where $`\mathrm{\Theta }(x)`$ is the Heavyside step function. In our simulations $`U=1.2`$, which produces a flow with a single connected chaotic region in the advection dynamics. The value of the numerically obtained Lyapunov exponent is $`\lambda _02.67/T`$. We use the source term $`S(x,y)=1+\mathrm{sin}(2\pi x)\mathrm{sin}(2\pi y)`$. Backward trajectories with initial coordinates on a square grid were calculated and used to obtain the chemical field at each point by integrating (13) forward in time. The smooth and the filamental behavior are found by changing the decay rate $`b`$, in complete agreement with the theory.

### B Smooth-filamental transition in a nonlinear plankton dynamics model

As we mentioned before, all the results presented in the previous section are valid for a more complex but stable chemical or biological dynamics, as far as the decay rate $`b`$ is replaced the chemical Lyapunov exponent $`|\lambda _C|`$. In this subsection we illustrate the smooth-filamental transition with a simple model of plankton dynamics immersed in a meandering jet flow. Additional results for nonlinear models can be found in .

As in the previous subsection, diffusion is neglected so that a Lagrangian representation in terms of a chemical and a transport dynamical system is possible. The chemical plankton model is the one used in . It is similar to (3), but with phytoplankton growth ruled by a logistic term (8), and a relaxational dynamics imposed on the carrying capacity $`C`$:

$`{\displaystyle \frac{dC}{dt}}`$ $`=`$ $`\gamma \left(C_0(𝐱)C\right)`$ (23)
$`{\displaystyle \frac{dP}{dt}}`$ $`=`$ $`P\left(1{\displaystyle \frac{P}{C}}\right)PZ`$ (24)
$`{\displaystyle \frac{dZ}{dt}}`$ $`=`$ $`PZ\delta Z^2.`$ (25)
In the absence of flow, the only stable fixed point of model (25) is given by $`C^{}=C_0(𝐱)`$, $`P^{}=C_0\delta /(\delta +C_0)`$, and $`Z^{}=P^{}/\delta `$.

The transport dynamical system is

$`{\displaystyle \frac{dx}{dt}}`$ $`=`$ $`{\displaystyle \frac{\psi }{y}}`$ (26)
$`{\displaystyle \frac{dy}{dt}}`$ $`=`$ $`{\displaystyle \frac{\psi }{x}},`$ (27)
given in terms of the following streamfunction :

$$\psi (x,y)=1\mathrm{tanh}\left(\frac{yB(t)\mathrm{cos}\left[k(xct)\right]}{\left(1+k^2B(t)^2\mathrm{sin}^2\left[k(xct)\right]\right)^{\frac{1}{2}}}\right).$$
(28)
It describes a jet flowing eastwards, with meanders in the North-South direction. These meanders are also advected by the jet at a phase velocity $`c`$. $`B(t)`$ and $`k`$ are the (properly adimensionalized) amplitude and wavenumber of the undulation in the streamfunction. As we are interested in a closed flow, we impose periodic boundary conditions at the ends of the interval $`L_x<x<L_x`$. Particles leaving the region through the right boundary are reinjected from the left. Nutrients are injected in and out continuously from fluid elements as these traverse the different regions of the carrying capacity source

$$C_0(x,y)=1+A\mathrm{sin}(2\pi x/L_x)\mathrm{sin}(2\pi y/L_y).$$
(29)
Chaotic advection appears in this model if $`B`$ is made to vary in time, for example periodically: $`B(t)=B_0+ϵ\mathrm{cos}(\omega t+\theta )`$. In our calculations we use the parameter values $`B_0=1.2`$, $`c=0.12`$, $`k=2\pi /L_x`$, $`L_x=7.5`$, $`L_y=4.0`$, $`\omega =0.4`$, $`ϵ=0.3`$, and $`\theta =\frac{\pi }{2}`$. These values guarantee the existence of ‘large scale chaos’, i.e, the possibility that a test particle crosses the jet passing from North to South or viceversa. This is weaker than the requirement of hyperbolicity, but is enough to illustrate the general aspects of our theory. In the chemical subsystem we use $`A=0.2`$, and $`\delta =2.0`$, and vary the value of $`\gamma `$, which controls the relaxational time-scale. It turns out that the resulting chemical dynamics is stable .

Fig. 4 shows snapshots of the long-time phytoplankton distributions. The smooth-filamental transition is clearly observed. More quantitative discussions can be found in reference .

### C Open flows

A finite region is said to be traversed by an open flow if almost all fluid elements enter the region and leave it forever in a finite time. A prototype is a stream passing around a cylindrical body. If the inflow velocity is high enough vortices form in the wake of the cylinder and make the flow time-dependent in this region, while the flow remains steady in front of the cylinder or in the far downstream region. Advected particles enter the unsteady region, undergo transient chaotic motion , and finally escape and move away downstream on simple orbits. The time spent in the mixing region, however, depends strongly on the initial coordinates, with singularities on a fractal set corresponding to particles trapped forever in the mixing region. This is due to the existence of a non-attracting chaotic saddle (a fractal object of zero measure) formed by an infinite number of bounded hyperbolic orbits in the mixing region. The stable manifold of this chaotic saddle contains orbits coming from the inflow region but never escaping from the mixing zone. Thus, permanent chaotic advection is restricted to this fractal set of zero Lebesgue measure. Points close to the unstable manifold of the chaotic saddle have spent a long time in the mixing region of the flow moving near chaotic orbits with a positive Lyapunov exponent. For points precisely on this unstable manifold, the backwards trajectories (the ones from which the Lyapunov exponent in (20) should be computed) remain in the chaotic saddle, thus leading to $`\lambda _0>0`$. The other trajectories spend in the mixing region just a finite time (both in the forward as in the backwards time direction), so that they can not contribute to the development of singularities in chemical distributions ((18) is singular only in the $`\overline{t}\mathrm{}`$ limit). In fact almost all trajectories are characterized by a long-time Lyapunov exponent equal to zero.

Thus open flows provide a rather clear example of strong space-dependence of Lyapunov exponents. According to Eq. (20), the Hölder exponent may be different from $`1`$ only on the unstable manifold of the chaotic saddle, thus implying that the transition from smooth to filamental structure now only takes place in this fractal set of zero measure. The background chemical field is always smooth, independently of the value of $`b`$.

Patterns of chemically decaying substances in a cylinder wake are presented in . Here we present the case of the stable biological model (25) under the jet flow (27)-(28). It is made open simply by not imposing the periodic boundary conditions of the previous subsection. The chaotic mixing region is not restricted to a finite region, but we restrict the source forcing (29) to $`L_x<x<L_x`$. This is the region shown in Fig. 5, being $`C_0=0`$ outside. We let the fluid particles to enter this region with very small $`C`$, $`P`$ and $`Z`$ concentrations. Most of them remain in the forcing region only for a finite time, so that, as in the cylinder wake case, they do not develop singularities. Only particles remaining in the chaotic saddle, i.e. the set of nonscaping orbits, forever will lead to diverging gradients. Fig. 5 shows a phytoplankton pattern for the same parameter values as before and $`\alpha =0.025`$. Smooth zones coexist with singular features. The behavior is distinct from the closed flow case (Fig. 4) and represents a new scenario for the development of chemical or biological patterns.

Since the irregularities now appear only on a set of measure zero, one could ask if they can have any significant effect on measurable quantities. In order to clarify this, we consider for simplicity the case of a single chemical reactive $`C`$, and instead of the previous characterization of the point-wise strength of the singularities by the Hölder exponent, let us investigate the scaling with distance $`\delta x|\delta 𝐱|`$ of the spatial average of the differences $`\delta C_{\mathrm{}}`$. For simplicity, let us assume that, on the saddle, there is no distribution in the local infinite-time Lyapunov exponents, i.e., that the advection on the chaotic saddle is characterized by a single Lyapunov exponent $`\lambda _0`$. In this case the partial fractal dimension (i.e. the dimension of intersections of the set with a line) of the manifolds of the chaotic saddle is

$$\stackrel{~}{D}=1\frac{\kappa }{\lambda _0}.$$
(30)
Here $`\kappa `$ is the escape rate, that is the rate of the exponential decay ($`e^{\kappa t}`$) of the number of fluid elements spending time longer than $`t`$ in the mixing or in the forcing region. On a one-dimensional transect of unit length the total number of segments of length $`\delta x`$ is $`(\delta x)^1`$ while the number of segments containing parts of the unstable manifold (with partial fractal dimension $`\stackrel{~}{D}`$) is $`\delta x^{\stackrel{~}{D}}`$. Thus, according to (20) the spatial average of $`\delta C_{\mathrm{}}`$ along this line, $`\delta C_{\mathrm{}}(𝐱;\delta x)`$, can be written as

$`\delta C_{\mathrm{}}(𝐱;\delta x)=(\delta x)(\delta x)^{\stackrel{~}{D}}(\delta x)^{b/\lambda _0}+`$ (31)
$`+(\delta x)[(\delta x)^1(\delta x)^{\stackrel{~}{D}}](\delta x)`$ (32)
where the first term pertains to the singular component ($`\alpha =b/\lambda _0`$), while the second one pertains to the smooth component ($`\alpha =1`$). In the limit $`\delta x0`$ the dominating behavior is

$$\delta C_{\mathrm{}}(𝐱;\delta x)\delta x^\zeta ,$$
(33)
with

$$\zeta =\mathrm{min}\{1,1+b/\lambda _0\stackrel{~}{D}\}=\mathrm{min}\{1,\frac{b+\kappa }{\lambda _0}\}$$
(34)
showing that if $`\stackrel{~}{D}<\frac{b}{\lambda _0}`$ (or, $`b+\kappa >\lambda _0`$) the average will be dominated by the smooth component, but if the fractal dimension of the singular set is large enough it contributes to the scaling of $`\delta C_{\mathrm{}}(𝐫;\delta r)`$. We see that moments of $`\delta C`$ may be sensible to fractal inhomogeneity in $`\alpha `$, or intermittency.

### D Anomalous scaling of structure functions

The strongly intermittent structure of singularities in open flows is an extreme example. There are additional inhomogeneities affecting both to the open and to the closed flows: although, in the long-time limit the Lyapunov exponent is the same for almost all trajectories in an ergodic region, deviations can persist on fractal sets of measure zero, and as we saw above such sets can contribute significantly to the global scaling. The origin of these inhomogeneities can be traced back by analyzing the finite-time distribution of Lyapunov exponents. In general, the finite-time stretching rates, or local Lyapunov exponents , have a certain distribution around the most probable value. This distribution approaches the time-asymptotic form :

$$P(\lambda ,t)t^{1/2}e^{G(\lambda )t}$$
(35)
where $`G(\lambda )`$ is a function characteristic to the advection dynamics, with the property that $`G(\lambda _0)=G^{}(\lambda _0)=0`$ and $`G(\lambda )>0`$, and $`\lambda _0`$ is the most probable value of the Lyapunov exponent. At infinitely-long times all the measure becomes concentrated at this single value $`\lambda _0`$, as stated before. There are however fractal sets of zero measure that do not share this unique value. The partial dimension of the set with Lyapunov exponent value $`\lambda `$ is  $`\stackrel{~}{D}(\lambda )=1G(\lambda )/\lambda `$. The coexistence of such different and interwoven fractal sets is the signature of multifractality.

For a robust quantitative characterization of the filamental structures, accessible to measurements and sensible to the intermittency features, we consider now the scaling properties of the structure functions associated with the chemical fields. For a single species $`C`$, the $`q`$th order structure function is defined as

$$S_q(\delta x)=|\delta C_{\mathrm{}}(𝐱;\delta x)|^q$$
(36)
where $`\mathrm{}`$ represents averaging over different locations $`𝐱`$, and $`q`$ is a parameter (we will only consider structure functions of positive order ($`q>0`$)). In the absence of any characteristic length over a certain range of scales the structure functions are expected to exhibit, as $`\delta x0`$, a power-law dependence

$$S_q(\delta x)\delta x^{\zeta _q}$$
(37)
characterized by the set of scaling exponents $`\zeta _q`$. The scaling exponent of the power spectrum is given by $`\beta =1+\zeta _2`$.

If the Hölder exponent of the field has the same value everywhere, given by (20) with $`\lambda =\lambda _0`$, the scaling exponents of the resulting mono-affine field are simply

$$\zeta _q=q\alpha _0=q\frac{b}{\lambda _0}.$$
(38)
(we have assumed $`b<\lambda _0`$). In general, the fractal sets of partial dimensions $`\stackrel{~}{D}(\lambda )`$ should be taken into account in an average such as (32) but for a continuous range of values of $`\lambda `$. The result for the scaling exponents in (37) is :

$$\zeta _q=\underset{\lambda }{\mathrm{min}}\{q,1+\frac{qb}{\lambda }\stackrel{~}{D}(\lambda )\}=\underset{\lambda }{\mathrm{min}}\{q,\frac{qb+G(\lambda )}{\lambda }\}$$
(39)
Equation (34) is a particular case of (39) for $`q=1`$ and in the approximation of considering a single value of $`\lambda `$ on the chaotic saddle.

Multifractality thus originates anomalous scaling (nonlinear dependence of $`\zeta _q`$ on $`q`$) of the structure functions. Multifractality and anomalous scaling have been observed in real data from plankton and temperature distributions (see for example . Here we present the scaling exponents for structure functions of the linearly decaying scalar (13) under the flow (22). Fig. 6 shows the numerically calculated scaling exponents (i.e. obtained by direct application of Eq. (36)), and the family of lines corresponding to (39) for different values of $`\lambda `$. Eq. (39) predicts that the actual values of $`\zeta _q`$ are given by the lower envelope of the set of lines, a fact confirmed by the numerical results. The values of $`G(\lambda )`$ can be found in .

Also shown in the same figure are the mono-fractal approximation ($`\zeta _q=qb/\lambda _0`$), which appears to be accurate for small $`q`$, and the function

$$\zeta _q=\sqrt{\left(\frac{\lambda _0}{\mathrm{\Delta }}\right)^2+\frac{2qb}{\mathrm{\Delta }}}\frac{\lambda _0}{\mathrm{\Delta }}.$$
(40)
that results from a parabolic $`G(\lambda )`$:

$$G(\lambda )=\frac{(\lambda \lambda _0)^2}{2\mathrm{\Delta }}.$$
(41)
This can be thought as the first term in a Taylor expansion around $`\lambda _0`$. It gives a good approximation to obtain the small-$`q`$ scaling exponents. Expression (40) becomes exact for the Kraichnan flow .

## V Excitable dynamics

Excitable dynamics  refers to the class of dynamical systems in which one can identify activator and inhibitor variables with the following properties: The activator displays some kind of autocatalitic growth behavior, but the presence of the inhibitor controls it so that the dynamical system has a stable fixed point as unique global attractor. The essence of the excitability phenomenon is the presence of a threshold, such that if there is a perturbation above it the system variables reach the stable fixed point only after a large excursion in phase space. This behavior usually appears when the activator has a temporal response much faster than the inhibitor, which then takes some time before stopping the growth of the activator.

Many chemical systems behave in this way, including the famous Belousov-Zhabotinsky reaction in adequate concentration ranges, or the electrochemical reactions occurring in nerves and muscles . Truscott and Brindley  identified phytoplankton as the fast activator and zooplankton as the slow inhibitor in models of the type (3), and demonstrated excitable behavior when Hollings-III grazing functions are used. This opens the possibility of understanding some plankton dynamics features in terms of well established scenarios found for excitable chemical reactions. In particular  explained features of red tides in terms of the essentials of excitable dynamics. In addition, transport processes coupled to excitable dynamics lead to a rich variety of pattern forming phenomena. The most widely studied is the appearance of excitable waves, of linear, circular, or spiral shape, in excitable media with diffusion.

In the next subsection we summarize the results of  concerning the generation of plankton patchiness (via plankton blooms) in excitable media without transport processes, and in subsection V B we will report some recent results for excitable dynamics under chaotic advection.

### A Localized blooms from excitable dynamics

Excitability is a threshold phenomenon: whenever it is crossed, the system undergoes a large excursion in state space, after which it returns back to the unexcited state. In an extended system, different regions may be in different phases of the excitation cycle. If transport processes are not efficient enough to couple different regions, excitation will occur only in the places where the excitation threshold has been reached, whereas the rest of the system will remain essentially at the equilibrium concentrations. It was realized in  that this simple mechanism leads to localized blooms in plankton models. For adequate parameter values, the temporal evolution of the bloom has the characteristics of red tides.

We illustrate this phenomenon with the model presented in . It is of the form (3) but reduced to two species ($`P`$-$`Z`$) with phytoplankton logistic growth and Hollings-III grazing. In convenient adimensional units it reads

$`F_P(P,Z)`$ $`=`$ $`\alpha P(1P){\displaystyle \frac{P^2}{P^2+P_e^2}}Z`$ (42)
$`F_Z(P,Z)`$ $`=`$ $`\gamma {\displaystyle \frac{P^2}{P^2+P_e^2}}ZmZ.`$ (43)
It presents excitable behavior in a range of parameters. In particular we show in Fig. 7 four stages of the evolution for $`\alpha =0.43`$, $`P_0=0.053`$, $`\gamma =0.005`$, and $`m=0.34`$. The system is distributed in the one-dimensional line, but each point evolves with the same dynamics (43) independently of the others. Its evolution thus depends only on its own initial conditions.

The first panel in Fig. 7 shows this initial condition: Phytoplankton is everywhere at the equilibrium fixed-point concentration, but zooplankton presents values somehow smaller than equilibrium in the central zone. In the region where this depletion is larger than a threshold value, fast phytoplankton growth occurs, followed by growth in zooplankton and subsequent recovery of the equilibrium values. The bloom is localized since there is neither advection nor diffusion to communicate excitation to neighboring places. Related behavior occurs if the initial perturbation consists in a local phytoplankton increase, or in a local increase of the phytoplankton growth rate $`\alpha `$.

### B Excitable dynamics under chaotic stirring

Except in the most quiescent situations, fluid flow will disturb the localized blooms of the previous section. Chaotic advection will deform the static excited patches into long and thin filaments subjected to stretching and folding. Reference  studied this process including also a small diffusion. The chemical dynamics considered was the FitzHugh-Nagumo system, a model originally developed in the context on nerve excitation, but it is clear that the different scenarios found would apply to any excitable chemical or biological dynamics in a chaotic flow.

Three qualitatively different kinds of behavior appeared. In the case in which the stirring time-scale (measured for example by the flow Lyapunov exponent) is much slower than the time scale for the excitation-deexcitation cycle, a localized perturbation propagates essentially as a deformed version of the circular waves that would propagate diffusively in a quiescent medium. In the other extreme case, when the stirring time-scale is faster than the time scale for growth of the activator, stretching of the filaments and the associated lateral contraction occur too fast. The largest gradients associated to filament narrowing increase the diffusive flux out of them, and they can not be filled-up with new excitation fast enough. As a result the stretched filaments become increasingly diluted and concentrations may fall below the excitation threshold. At this point the excitation will disappear. This means that a too vigorous stirring will eliminate a small localized excitation. Perhaps the most surprising behavior appears in a range of intermediate stirring speeds, such as the dilution effect stops the growth of the slow inhibitor, but not of the fast activator: The activator concentration inside the filament remains saturated, the lateral thinning becomes stopped at a finite width by the combined effect of excitation and diffusion, and the filament length increases continuously as a consequence of chaotic stretching. After some time, in a closed flow after repeated folding, the filament will fill-up the whole system. The effect of chaotic advection has changed the character of an initial excitation from localized to global. After the whole system is excited, deexcitation proceeds as in a homogeneous well-mixed medium.

The process is illustrated in Fig. 8, which is obtained  by integration of advection-reaction-diffusion equations (4), with the reaction terms as in the previous subsection (Eq. (43)). Advection and diffusion for a nutrient field is also included (with $`F_N=0`$). The coupling of the nutrient with phytoplankton is made via a dependence of the growth rate $`\alpha `$ on the nutrient concentration ($`\alpha =\alpha _0(1+N/N_0)`$). For the flow (closed and incompressible) a set of randomly seeded eddies has been used. The inset in the second panel of Fig. 8 shows a chlorophyll filament observed from the satellite sensor SeaWiFS 42 days after a fertilization experiment in the Antartic. Reference  interprets this filament as the result of stretching and folding of the initial excitation introduced in the fertilization experiment. This would support the interpretation of the observed dynamics in terms of the concepts of excitability and chaotic advection. The phenomenon reported here, i.e. the occurrence of a global excitation by the effect of chaotic advection on a localized patch, gives a warning about the possibility of responses of unexpected extension arising from very localized perturbations.

The essential object in this kind of phenomena is the stretched filament. In , simplified equations describing the transverse filament profile are discussed. They are essentially of the form (11) adapted to the multicomponent case, and with the simple linear chemistry term $`\mu (t)P`$ replaced by the excitable chemistry. These filament equations can be used to obtain analytical understanding of the transitions between the different regimes, and to get quantitative results. For example, the width of the excited filament is given by $`w=c\sqrt{D\alpha }/\lambda `$, where $`c`$ is a numerical constant depending on the particular excitable model used, $`D`$ is the diffusion coefficient, $`\alpha `$ the phytoplankton linear growth rate, and $`\lambda `$ the exponential contraction given by the strain, that can be identified with the flow Lyapunov exponent.

## VI Some warnings from individual modeling

In the present Lectures we have restricted our considerations to continuous differential equation models describing continuous concentration variables. The use of maps  is an alternative to the use of continuous differential equations. But there is an aspect of real chemical or biological interacting systems that is missed in both approaches: the discrete nature of the interacting entities (chemical molecules or the individual organisms). On general grounds, one would expect just a more noisy dynamics in models in which individuals are resolved. There are cases however in which individual quantization produces more profound effects. We summarize here the results of two recent references  that highlight this fact in very simple examples.

Reference  discusses the following chemical scheme:

$`A+C`$ $``$ $`2C`$ (44)
$`C`$ $``$ $`0`$ (45)
The first autocatalitic equation can be interpreted as a simple process in which the $`C`$ particle eats ‘food’ $`A`$ to self-replicate (the amount of $`A`$ is kept constant), and the second represents the death of $`C`$. The continuous description in terms of a reaction-diffusion equation for the concentration $`C`$ would be

$$\frac{C}{t}=(\lambda A\mu )C+D^2C$$
(46)
$`\lambda `$ and $`\mu `$ are the rates of the first and second reaction in (45) respectively. The continuous model predicts extinction of $`C`$ if the effective growth rate $`\lambda A\mu `$ is negative.

Alternatively one can model the process by a large number of random walkers, each representing a molecule of $`A`$ or $`C`$, that undergo the transformations (45) following some Monte Carlo dynamics. It happens  that if the effective growth rate is not too negative, the amount of $`C`$ increases instead of decreasing. The consideration of the discrete nature of the individuals has changed death to life. What happens is that the $`C`$ individuals die in most of space, but they survive if close enough to ‘food’ particles. The effect is a diverging clustering of $`C`$ particles around $`A`$ random walkers. This strong inhomogeneity is missed in the continuous model, that gives the wrong prediction. The root of the problem is in the use of mean multiplicative growth rates in (46), whereas in fact they vary randomly in space and time, reflecting the positions of the $`A`$ particles. It is known that mild randomness in the multiplicative growth factors in linear equations such as (46) give rise to strongly intermittent distributions. Intermittency reflects here in the very inhomogeneous clustering of the surviving particles.

A related phenomenon is observed in . There, the chemical scheme is

$`C`$ $``$ $`2C`$ (47)
$`C`$ $``$ $`0`$ (48)
which leads to essentially the same continuous description (46), with $`\lambda A`$ replaced by $`\lambda `$. Here again growth is observed in discrete simulations when the continuous model indicates extinction. Here there are no ‘food’ particles to cluster around, but clustering occurs around the initially seeded $`C`$ particles. There are position correlations between a parent particle and its descendants (the ones appearing via the first autocatalitic reaction in (48)), that are completely neglected in the continuous model.

The authors of  check that the reproductive correlations persist in the presence of chaotic advection: the clusters become filaments, but the anomalous survival remains. The authors are able to describe the correct behavior via continuous equations, but not for the average concentration, but for binary correlations. Whether or not a given discrete process can be modeled in a continuum model, and up to which level of statistics one should retain randomness and correlations, is a major question in reacting particle models and in general ecosystem dynamics.

## VII Conclusions

We have presented a number of processes, inspired by concepts in Nonlinear Dynamics such as chaos and excitability, that can be useful to understand generic behaviors in chemical or biological systems in fluid flows.

For the case of the linearly decaying tracer, a rather complete description is available. The scaling exponents of the set of structure functions, which display anomalous (multifractal) behavior, can be expressed in terms of the decay-time constant and the probability distribution of the Lagrangian finite-time Lyapunov exponents of the flow. The decay-time constant makes the tracer power spectrum steeper than the Batchelor law obtained for the passive tracer. The differences between open and closed flows are important, and the relevance of these results for nonlinear plankton models has been be pointed out.

The study of excitable-type reactions is motivated by the presence of excitability in some plankton dynamics models. It contains the essentials of the observed characteristics of plankton blooms. Important changes in the response of excitable models are found when they are placed in a chaotic fluid flow. This may be relevant for the outcome of ocean fertilization experiments.

Finally, some warnings have been given about the difficulties in modeling discrete individuals (such as planktonic organisms) in terms of continuous concentration fields. Aggregation behavior, reproductive correlations, and other phenomena turn out to be difficult to introduce in standard continuous models, whereas they are quite naturally present in individual-based descriptions.

## ACKNOWLEDGMENTS

We acknowledge Susana Agustí for providing us the data shown in Fig.1, Adrian P. Martin for motivating our interest in some of the topics presented here, and Carlos Duarte, Veronique Garçon, Peter H. Haynes, Oreste Piro, Tamás Tél, and Angelo Vulpiani for illuminating discussions. We thank also the organizers of the 2001 International Summer School on Atmospheric and Oceanic Sciences (ISSAOS 2001) for the excellent opportunity to discuss the subject of Chaos in Geophysical Flows. Financial support from MCyT (Spain) projects IMAGEN REN2001-0802-C02-01/MAR and CONOCE BFM2000-1108 is greatly acknowledged. C.L. acknowledges support from MECD (Spain).


format: text

text:
# 1 Introduction

## 1 Introduction

One of the major challenges in QCD is the description of high energy scattering phenomena. In the high center-of-mass energy $`\sqrt{s}`$ and in the perturbative domain when $`\alpha _s1`$ the scattering amplitude is obtained by the summation of diagrams leading in $`\mathrm{log}s`$ . At this level, when $`\alpha _s`$ is frozen, the dependence of the resulting cross section on the energy is governed by the power law $`x^{\omega _{\mathrm{I}\mathrm{P}}}`$ where $`x`$ is the Bjorken variable. The critical exponent $`\omega _{\mathrm{I}\mathrm{P}}=4\mathrm{ln}2\overline{\alpha }_s`$ ($`\overline{\alpha }_s=\alpha _sN_c/\pi `$) is provided by the minimum of the eigenvalue function $`\chi (\gamma )`$ of the BFKL evolution kernel.

The conceptual problem in this approach is the fact that at sufficiently high center-of-mass energies the BFKL amplitude violates the Froissart unitarity bound. This means that the validity of this approach is strongly limited and it has to be modified at very small $`x`$ in order to guarantee the unitarity of the resulting cross section.

The solution to this problem can be provided by including the recombination effects, which are likely to occur at very small values of $`x`$ . By decreasing the value of $`x`$ at fixed gluon virtuality $`k_t^2`$, the density of partons becomes so large that they start to overlap. In this case the gluon splitting process must be supplemented by a competing gluon recombination. In terms of evolution the master equations become nonlinear with an additional quadratic term which suppress the growth of the amplitude with energy and restore the unitarity222Recently  it has been pointed out that the situation can be actually more complicated in a sense that the Balitsky-Kovchegov equation  could lead to the local saturation but not to the unitarisation due to the fact that the target radius in impact parameter space could grow as fast as a power with energy..

There have been extensive studies on this problem, see , -, which result in the nonlinear evolution. One of the important outcomes of these studies is the existence of the saturation scale $`Q_s(x)`$ which is a characteristic scale at which the parton recombination effects become important. In particular case of the Balitsky-Kovchegov equation  the existence of such scale has yet another important impact on the picture of the BFKL evolution. The diffusion into the infra-red, which is the charactersitic property of BFKL evolution is strongly limited due to the existence of the saturation scale . In fact, in the regime when the gluon transverse momenta $`k<Q_s(x)`$, the solution to the nonlinear equation  becomes a function of only one combined variable $`k/Q_s(x)`$ -. In the regime of high momenta, $`k>Q_s(x)`$ the parton density is small and the evolution is governed by a normal linear equation.

In this talk we present the analytical and numerical analysis of the Balitsky - Kovchegov  equation which is a nonlinear evolution equation in the leading $`\mathrm{log}s`$ limit. We illustrate the emergence of the saturation scale and scaling and show that it leads to the suppression of the infra-red diffusion. We also consider the case with additional NLL effects such as kinematical constraint and running coupling.

The results presented in this talk have been obtained in the collaboration with K. Golec - Biernat and L. Motyka. For the details of the calculation the reader is referred to .

## 2 Nonlinear evolution equation

The Balitsky-Kovchegov equation  has been derived as an evolution equation for the dipole-nucleus amplitude in the dipole picture by a summation of multiple Pomeron exchanges in the leading $`\mathrm{log}s`$ level and in the large $`N_c`$ limit.

The resulting evolution equation reads

$$\begin{array}{c}\frac{N(𝐫,𝐛,Y)}{Y}=\overline{\alpha }_s(KN)(𝐫,𝐛,Y)\hfill \\ \hfill \overline{\alpha }_s\frac{d^2𝐫^{}}{2\pi }\frac{r^2}{r^2(𝐫+𝐫^{})^2}N(𝐫+𝐫^{},𝐛+\frac{𝐫^{}}{2},Y)N(𝐫^{},𝐛+\frac{𝐫+𝐫^{}}{2},Y),\end{array}$$
(1)
where $`\overline{\alpha }_s=N_c\alpha _s/\pi `$, and the linear term is determined by the BFKL kernel

$$\begin{array}{c}(KN)(𝐫,𝐛,Y)=\hfill \\ \hfill =\frac{d^2𝐫^{}}{\pi r^2}\left\{\frac{r^2}{(𝐫+𝐫^{})^2}N(𝐫+𝐫^{},𝐛+\frac{𝐫^{}}{2},Y)\frac{r^2}{r^2+(𝐫+𝐫^{})^2}N(𝐫,𝐛,Y)\right\}.\end{array}$$
(2)
The function $`N(𝐫,𝐛,Y)`$ is the dipole-nucleus amplitude for the scattering of the dipole with transverse size $`𝐫`$ at impact parameter $`𝐛`$ and at rapidity $`Y`$.

In the linear approximation, when each dipole scatters only once off the nucleus, the BFKL equation in the dipole picture is obtained. The non-linear term in (1) takes into account multiple scatterings and is essentially determined by the triple pomeron vertex  in the large $`N_c`$ limit. Eq. (1) unitarizes the BFKL pomeron in the sense that at $`x0`$ and $`Q^2`$ fixed,

$$F_2Q^2\mathrm{ln}(1/x).$$
(3)
Thus the power-like rise with energy for the BFKL pomeron is tamed .

For the subsequent analysis we shall assume the approximation of the big nucleus, i.e. $`𝐫𝐛`$ which allows us to factorize the impact parameter dependence in Eq.(1). We also consider the spherical symmetric solutions in $`𝐫`$ and transform the equation (1) into the momentum space by performing the Fourier transform

$$\varphi (k,Y)=\frac{d^2𝐫}{2\pi }\mathrm{exp}(i𝐤𝐫)\frac{N(r,Y)}{r^2}=_0^{\mathrm{}}\frac{dr}{r}J_0(kr)N(r,Y),$$
(4)
where $`J_0`$ is the Bessel function. In this case the following equation is obtained

$`{\displaystyle \frac{\varphi (k,Y)}{Y}}=\overline{\alpha }_s(K^{}\varphi )(k,Y)\overline{\alpha }_s\varphi ^{\mathrm{\hspace{0.17em}2}}(k,Y),`$ (5)
and the action of the BFKL kernel is given by

$$(K^{}\varphi )(k,Y)=_0^{\mathrm{}}\frac{dk^2}{k^2}\left\{\frac{k^2\varphi (k^{},Y)k^2\varphi (k,Y)}{|k^2k^2|}+\frac{k^2\varphi (k,Y)}{\sqrt{4k^4+k^4}}\right\},$$
(6)
where now $`k`$ and $`k^{}`$ are the transverse momenta of the exchanged gluons in the BFKL ladder.

## 3 Saturation scale and geometric scaling

In order to study the diffusion and scaling properties of the solution we shall consider the function $`k\varphi (k,Y)`$. In the case of the linear BFKL equation this function is a Gaussian concentrated around some initial scale $`k_0`$ and with width increasing with rapidity, leading to a diffusion. In Fig.1 we illustrate this distribution for the case of the solution to the linear BFKL and in the nonlinear Balitsky-Kovchegov equation for different choices of rapidity. As an initial condition we have chosen a delta function $`\delta (kk_0)`$. In the case of the linear BFKL evolution the solution is always peaked at $`k=k_0`$ and exhibits the well known diffusion pattern, in which the momentum distribution has an increasing width with growing rapidity. In the nonlinear case however the solution behavies quite differently. The peak of the distribution $`k_{\mathrm{max}}`$ moves towards higher $`k`$ values as the rapidity increases and the solution becomes washed out from the $`k<k_0`$ regime. Only in the region $`kk_{\mathrm{max}}`$ it coincides with the linear evolution.

The impact of unitarization of the BFKL pomeron on the infra-red behaviour can be also visualised by studying the properties of the following normalised distribution

$$\mathrm{\Psi }(k,Y)=\frac{k\varphi (k,Y)}{k_{\mathrm{max}}(Y)\varphi (k_{\mathrm{max}}(Y),Y)},$$
(7)
and by performing the projection of this function onto the ($`\mathrm{ln}k/k_0,Y`$) plane, Fig. 2.

Again, for small $`Y`$, when the non-linearity in the BK equation is negligible, the re-normalized solutions (7) of the BFKL and the BK equations coincide. With increasing Y, when the non-linear effects become important, the difference between them in the region of small $`k`$ becomes fully visible.

Note that in certain region of ($`\mathrm{ln}k/k_0,Y`$) in Fig. 2 the contours become parallel straight lines. This means that $`\mathrm{\Psi }(k,Y)`$ in this region is a function of the combination

$$\xi =\mathrm{ln}(k/k_0)\lambda Y=\mathrm{ln}\left(\frac{k}{k_0\mathrm{exp}(\lambda Y)}\right)\lambda >0,$$
(8)
instead of $`k`$ and $`Y`$ separately. This is referred as the geometric scaling property. We note, that this scaling holds in the regime where $`k<k_{\mathrm{max}}`$ and becomes violated when $`k>k_{\mathrm{max}}`$. This suggests that we can identify the saturation scale $`Q_s(Y)`$ with $`k_{\mathrm{max}}`$

$$Q_s(Y)k_{\mathrm{max}}(Y)=Q_0\mathrm{exp}(\lambda Y),$$
(9)
with the exponential dependence on rapidity governed by the value of the scaling parameter $`\lambda `$.

The solution to Eq. (5) has the same scaling property as the function $`\psi (k,Y)`$ namely

$$\varphi (k,Y)=\varphi (k/Q_s(Y)),$$
(10)
provided that $`\varphi (k_{max},Y)=\mathrm{const}`$. We have checked that this condition is satisfied.

We have checked that the scaling coefficient $`\lambda `$ (9) is a universal quantity and it does not depend on the type of the input distribution whereas the normalisation $`Q_0`$ is determined by the initial conditions.

The transition between the scaling and linear regime can be perhaps best seen in Fig. 3 where the function $`k/Q_s(Y)\varphi (k,Y)`$ is plotted as a function of the scaling variable $`k/Q_s(Y)`$ for different values of the rapidity $`Y`$. The scaling behaviour is represented by a common line to the left of the point $`k/Q_s(Y)=1`$. The slow departure from scaling for $`k>Q_s(Y)`$ is clearly visible. In this sense the line $`k=Q_s(Y)`$ is only an approximation characterizing the position of the transition region in the $`(k,Y)`$-plane. However, the choice based on the position in $`k`$ of the maximum of $`k\varphi (k,Y)`$ as a function of $`Y`$ is the most natural one.

From the numerical solution we have extracted the value of the scaling parameter and found $`\lambda 2.05\overline{\alpha }_s`$ which is in agreement with the estimates of  for the Balitsky-Kovchegov equation and also with the previous work .

The geometric scaling is also the property of the Golec-Biernat and Wüsthoff saturation model  which succesfully described the data on inclusive and diffractive proton structure function. In this model it is assumed that the whole energy dependence of the dipole cross section $`\sigma _d(r,Y)=\sigma _d(r^2Q_s^2(Y))`$ is driven through the saturation scale $`Q_s(Y)`$ (or saturation radius $`1/Q_s(Y)`$). It then results in an approximate scaling property of the data for the total photon-proton cross section

$$\sigma _{\gamma ^{}p}(x,Q^2)=\sigma _{\gamma ^{}p}(Q^2/Q_s^2(Y)).$$
(11)
Such scaling law was found in the small-$`x`$ DIS data .

## 4 Analysis beyond LL level

The Balitsky-Kovchegov equation has been derived in the leading $`\mathrm{log}s`$ level, therefore it is important to study the impact of the subleading corrections. Although it has been argued  that the unitarity corrections might become important before next-to-leading ones, the study of the latter ones is crucial due to the large numerical value of the subleading series. After the evaluation of the NLL contribution to the BFKL amplitude has been performed  it turned out that the correction is very large and need to be resummed  in order to stablize the final result.

We have tested two important types of the NLL corrections: the inclusion of running of the coupling $`\overline{\alpha }_s(k^2)`$ and the so called kinematical constraint.

In the case of linear BFKL the running of the coupling poses serious problems due to the existence of the Landau pole and thus it is necessary to regularise $`\overline{\alpha }_s(k^2)`$ at small scales. This results in a strong dependence of the solution on the cut-off (or freezing ) parameter $`k_0`$.

The intercept of the BFKL Pomeron turns out to be dominated by the values of $`\overline{\alpha }_s(k_0^2)`$ and instead of the typical diffusion pattern one has a factorized behaviour of the solution $`k\varphi (k,Y)\mathrm{exp}\{\lambda Y\}\frac{1}{k}[\mathrm{ln}\left(k^2/k_0^2\right)]^\nu `$ at large rapidities. The distribution of the gluon momenta is dominanted by the virtualities in the infrared regime $`k_0`$.

In the case of the Balitsky-Kovchegov equation with the running coupling one has still problem of regularisation around the Landau pole, however the solution is itself much more stable with respect to the details of the phenomenological regularisation. It turns out, that since the saturation effects are very strong in the regime of small values of $`k`$, they tend to decrease the rapid rise of the amplitude for the values of running coupling evaluated at scales around the cutoff scale $`k_0`$. The rapidity - dependent saturation scale, $`Q_s(Y)`$ shifts the momentum distribution out of the infrared regime into the perturbative domain. This phenomenon can be best visualized by means of similar contour plots as before, Fig.4. In the linear BFKL case one observes that the solution initially undergoes the diffusion pattern, slightly asymmetric due to corrections coming from the running coupling, and then suddenly it drops into the small scales regime $`kk_0`$ where it exhibits factorized (in $`Y`$ and $`\mathrm{ln}k/k_0`$) behaviour. On the contrary, the Balitsky -Kovchegov equation which initially coincides with BFKL, also undergoes some form of transition but then the distribution moves away from the infrared regime due to the generation of the saturation scale.

We have estimated the saturation scale in the case of Balitsky-Kovchegov equation with running coupling. By taking the ansatz that the local exponent $`\lambda (Y)=d\mathrm{ln}(Q_s(Y)/\mathrm{\Lambda })/dY`$ takes the form $`\lambda (Y)=2\overline{\alpha }_s(Q_s^2(Y))`$ where $`\mathrm{\Lambda }=\mathrm{\Lambda }_{QCD}`$ we have found that the saturation scale is given by the solution to the differential equation

$$\frac{d\mathrm{ln}(Q_s(Y)/\mathrm{\Lambda })}{dY}=\frac{12}{b_0\mathrm{ln}(Q_s(Y)/\mathrm{\Lambda })},$$
(12)
with the initial condition $`Q_s(Y_0)=Q_0`$ and $`Y_0`$ chosen in the region where scaling sets in. The solution takes the form

$$Q_s(Y)=\mathrm{\Lambda }\mathrm{exp}\left(\sqrt{\frac{24}{b_0}(YY_0)+L_0^2}\right),Y>Y_0,$$
(13)
where $`L_0=\mathrm{ln}(Q_0/\mathrm{\Lambda })`$. It follows that the local exponent $`\lambda (Y)`$ decreases with increasing rapidity, and $`\lambda (Y)1/\sqrt{Y}`$ for very large $`Y`$.

Such dependence is indeed seen in the numerical analysis.

In Fig. 5 we illustrate scaling in the case with running coupling by showing the function $`(k/Q_s(Y))\varphi (k,Y)`$ plotted versus $`k/Q_s(Y)`$ for different values of rapidity. $`Q_s(Y)`$ is given by formula (13) with the initial condition $`Q_s(Y=0)=2\text{GeV}`$. The overlapping curves at low values of the scaling variable clearly indicate that for $`k<Q_s(Y)`$ scaling is satisfied to a very good accuracy, thus justifying our ansatz (13) for the saturation scale. We have also tested the impact of the so called kinematical constraint onto the solution of Balitsky-Kovchegov equation and found that though qualitative features are the same (with generation of the saturation scale and scaling) the absolute numerical value of the solution is strongly decreased and also the saturation scale is shifted towards smaller value of gluon momenta.

To summarize, the diffusion into infrared in the case of the Balitsky-Kovchegov equation is strongly suppressed due to the emergence of the saturation scale which can be approximately identified with the maximum of the momentum distribution $`k\varphi (k,Y)`$. The solution to the nonlinear equation has the property of the geometric scaling in the regime where $`k<Q_s(Y)`$ whereas in the case when $`k>Q_s(Y)`$ the solution enters the linear regime. In the case of running coupling the nonlinear equation turns out to be more stable as compared with the linear BFKL evolution. The sensitivity to the treatment of the infra-red region is much smaller than in the linear case due to the appearance of the saturation scale. In contrast to the BFKL equation with the running coupling, the large $`Y`$ asymptotics of the dipole scattering amplitude arising from the BK equation is governed by gluon momenta in the perturbative domain.

## Acknowledgments

This research was supported in part by the EU Framework TMR programme, contract FMRX-CT98-0194 and by the Polish Committee for Scientific Research grants Nos. KBN 2P03B 120 19, 2P03B 051 19, 5P03B 144 20.


text:
# Fluctuation of the Correlation Dimension and the Inverse Participation Number at the Anderson Transition

## Abstract

The distribution of the correlation dimension in a power law band random matrix model having critical, i.e. multifractal, eigenstates is numerically investigated. It is shown that their probability distribution function has a fixed point as the system size is varied exactly at a value obtained from the scaling properties of the typical value of the inverse participation number. Therefore the state-to-state fluctuation of the correlation dimension is tightly linked to the scaling properties of the joint probability distribution of the eigenstates.

It is already well established MJ ; rev ; MEM  that the eigenstates at the localization–delocalization transition driven by disorder are multifractals characterized by a set of generalized dimensions, $`D_q`$. This is reflected both in the statistical properties of each eigenstate and in the joint probability distribution of the eigenvector components. In the first case the multifractal analysis is performed on each eigenstate using the usual box–counting algorithm, while in the second case the scaling properties of the distribution function are analyzed.

The generalized dimensions are associated to the scaling of the $`q`$th moment of the eigenstates, $`\psi _\alpha (𝐫)`$. These moments are also termed as generalized inverse participation numbers MJ ; rev

$$I_q=d^dr|\psi _\alpha (𝐫)|^{2q}_\alpha L^{D_q(q1)},$$
(1)
where $`L`$ is the linear dimension of the system. The inverse of $`I_2`$, for example, roughly equals the number of nonzero wave function components, therefore it is a good and widely accepted measure of the extension of the states. For delocalized states $`I_21/N`$, while for localized ones $`I_2O(1)`$. At criticality it scales with $`L`$ with exponent $`D_2`$ (1). This scaling may be found for each state $`\psi _\alpha (𝐫)`$ but then the corresponding exponents fluctuate.

Another property of multifractality is related to the probability overlap of eigenstates with energy separation substantially exceeding the mean level separation. In particular $`D_2`$ describes these density correlations: hence the name of correlation dimension. It has been shown that for multifractal eigenstates these correlations decay slowly FM  no matter how sparse these states are. The role of $`D_2`$ is even more important as it has been related to the level compressibility at criticality CLK ; VK  that is a remarkable relation that has been corroborated at least for the case of weak multifractality.

The multifractal properties of the wave functions at criticality also shows up in the anomalous spreading of a wave packet, where the value of $`D_2`$ plays an essential role Bodo .

The importance of the generalized dimensions and especially of $`D_2`$ therefore is obviously hard to overestimate. This is the reason why it came as a surprise when in P-S  the authors demonstrated that the correlation dimension, $`D_2`$, has a relatively broad distribution even with increasing systems size of a three dimensional Anderson model. The authors concluded that this broadness remains in the thermodynamic limit and therefore there is no unique value of $`D_2`$ in this limit. Such a statement would imply serious objections against one–parameter scaling theory which is corroborated by many other techniques rev . Based on numerical simulations on a random matrix model having eigenstates that exhibit multifractality, however, it has been shown in E-M ; EC  that the fluctuations of $`I_2`$ can indeed be characterized by a universal probability distribution function (PDF), $`\mathrm{\Pi }(I_2)`$ once $`I_2`$ is normalized with the typical value, $`\stackrel{~}{I}_2=I_2/I_2^{typ}`$, where $`I_2^{typ}=\mathrm{exp}\mathrm{ln}I_2`$:

$$\stackrel{~}{\mathrm{\Pi }}(\mathrm{ln}\stackrel{~}{I}_2)=\mathrm{\Pi }_N(\mathrm{ln}I_2/I_2^{typ}).$$
(2)
So the PDFs are simply horizontally shifted by

$$\mathrm{ln}I_2^{typ}=\mathrm{ln}I_2\stackrel{~}{D}_2\mathrm{ln}N$$
(3)
to put on top of each other and obtain the universal curve of $`\stackrel{~}{\mathrm{\Pi }}(\mathrm{ln}\stackrel{~}{I}_2)`$ showing that there is indeed a unique value of $`\stackrel{~}{D}_2`$ to characterize the scaling of $`I_2`$. In the present work we wish to reconcile the apparent contradiction existing about the fluctuations of the correlation dimension $`D_2`$ and to show how to correct the misinterpretation of P-S  and obtain results that are compatible with E-M ; EC .

To this end we performed numerical simulations on the power law band random matrix ensemble that has been recently intensively investigated E-M ; EC ; EC2 ; prbm . This model is more suitable for a quantitative analysis since it allows to obtain a wider range of linear extension ($`N=L`$) of the system as compared to the usual Anderson model defined over an $`N=L^3`$ lattice of linear size $`L`$. Moreover, the model is capable to show many of the features of the traditional Anderson transition rev ; VK ; prbm . We will show that indeed the PDF of the correlation dimension of the states is a broad function which eventually has a fixed point at a value very close to $`\stackrel{~}{D}_2`$ obtained from the scaling of $`I_2^{typ}`$. Moreover, from the average value, $`D_2`$, one can extract the same answer as well.

The elements in these ensembles of $`N\times N`$ complex hermitian matrices are Gaussian distributed random variables with zero mean, $`H_{i,j}=0`$, and a variance

$$(H_{i,j})^2=\left[1+(|ij|/b)^2\right]^1\times \{\begin{array}{cc}1/2,\hfill & ij\hfill \\ 1/\beta ,\hfill & i=j\hfill \end{array}$$
(4)
where $`\beta =1`$ or $`2`$ depending on the global symmetry of the system. In the presence of time reversal symmetry, $`\beta =1`$, and the matrices are real symmetric, in its absence, $`\beta =2`$. The parameter $`b`$ is an effective bandwidth that serves as a continuous control parameter over a whole line of criticality VK . In this work we set $`b=1`$ as this proves to give the best convergence towards a scale invariant distribution function E-M . However, at the end we will give the $`b`$-dependence of the results as well.

The states have been obtained using standard diagonalization routines for the cases of $`N=128,\mathrm{},2048`$. The properties of the states in the vicinity of the middle of the band have been collected. We chose 25% of the states around $`E=0`$, and verified that by narrowing this energy window the results do not alter but obviously the statistics becomes worse.

Each individual eigenstate may show scaling with different exponents $`D_2`$. On the other hand the distribution of the IPN, $`I_2N^{D_2}`$, of the eigenstates shows universal scaling properties rev , thus a relation between these properties is expected. The distribution function of $`D_2`$ can be connected to that of the $`I_2`$ (2) as

$$𝒫_N(D_2)dD_2=\mathrm{\Pi }_N(I_2)dI_2$$
(5)
which yields

$$𝒫_N(D_2)=\mathrm{\Pi }_N(N^{D_2})N^{D_2}\mathrm{ln}(N).$$
(6)
We expect the distribution of $`D_2`$ to tend to a universal distribution as well as that of $`I_2`$. Evers and Mirlin E-M  have proved that in particular $`I_2`$ is a self–averaging quantity and numerically confirmed that for the case of the present model. The same has been shown for the case of the Anderson model, too EC ; MEM . Therefore as $`N\mathrm{}`$ we expect that the distribution

$$𝒫_N(D_2)\delta (D_2\stackrel{~}{D}_2).$$
(7)
This is in direct contradiction to the conclusions of P-S .

Furthermore using Eq. (6) for $`𝒫(D_2)`$ and Eq. (3) we can readily derive expressions of the mean value and the variance of $`D_2`$ in relation to the mean value and variance value of $`\mathrm{ln}I_2`$ as

$$D_2\stackrel{~}{D}_2=\frac{\text{const}}{\mathrm{ln}N}0$$
(8)
and

$$\text{var}(D_2)=\frac{\text{var}(\mathrm{ln}I_2)}{(\mathrm{ln}N)^2}0,$$
(9)
where the convergence to zero for the variance is ensured by EC ; EC2

$$\text{var}(\mathrm{ln}I_2)=\sigma _{\mathrm{}}aN^y,$$
(10)
with the irrelevant scaling $`y=\stackrel{~}{D}_2/(2\beta )>0`$ and also $`a>0`$ and $`\sigma _{\mathrm{}}>0`$. These logarithmic convergences, Eqs. (8) and (9), seem to be too difficult to detect for the case of the traditional Anderson model P-S , even though, Ref. EC ; EC2  showed the scaling (10) very clearly.

In Fig. 1 first we show that indeed the scaling of the PDF $`\mathrm{\Pi }_N(\mathrm{ln}I_2)`$ holds and obtain the value of $`\stackrel{~}{D}_2`$ from the typical value, $`\mathrm{ln}I_2`$ for both $`\beta =1`$ and $`2`$. Their values are listed in Table 1. As a next step we turn to the PDF of the correlation dimension, $`D_2`$. The correlation dimension of the individual eigenstates have been obtained using the standard box-counting algorithm MJ . In Figs. 2 and 3 we show how $`𝒫(D_2)`$ evolves by doubling $`N`$. These figures give the same information as Fig. 3 of P-S  for the case of the power law band random matrix ensemble. By increasing $`N`$ the distribution gets narrower, moreover, there is a fixed point right in the interval of $`0.7<D_2<0.8`$ for $`\beta =1`$ and at $`0.85<D_2<0.9`$ for $`\beta =2`$. In the insets of both figures we plot this region enlarged in a semi-log plot. The fixed point clearly appears at $`D_2^{}0.75`$ for $`\beta =1`$ and $`D_2^{}0.87`$ for $`\beta =2`$. The first value is precisely the one obtained via the scaling of the typical value $`\mathrm{ln}I_2`$ in E-M  and is close to our numerical results in Fig. 1.

We argue that due to the self–averaging property of $`I_2`$ E-M  $`𝒫(D_2)`$ should tend to a $`\delta `$–function so for finite $`N,N^{}`$ there should be a fixed point in the distribution function,

$$𝒫_N(D_2^{})=𝒫_N^{}(D_2^{}),$$
(11)
provided $`D_2^{}`$ falls within the set of $`D_2`$’s of the eigenstates. This ensures that $`D_2^{}=\stackrel{~}{D}_2`$ is the unique value of the correlation dimension at $`N\mathrm{}`$. Moreover it is also the average value, $`D_2`$ in this limit (8).

The convergence of $`D_2`$ towards $`\stackrel{~}{D}_2`$ is given in Fig. 4. The value of $`D_2`$ extracted from this fit together with the other estimates of $`D_2`$ are listed in Table 1.

The error bars characterize the goodness of the linear fits as depicted in Figs. 1 and 4. The different approaches give slightly different values, some of them coincide within a few percent difference. We attribute the discrepancies to possible finite size effects and to the fact that the number of samples has not been large enough. For conserved time reversal symmetry, $`\beta =1`$, we observe a coincidence between our $`D_2`$ and $`D_2^{}`$ values, together with the one obtained by E-M  over a larger ensemble of samples and a wider range of matrix size. All the values from different estimates seem to be within a range of 1–8% difference. Therefore we may conclude that despite the relatively large state-to-state fluctuations of the correlation dimension depicted in a wide PDF of $`𝒫(D_2)`$ in Figs. 2 and 3 one can still define a unique value of $`D_2`$ that is characteristic to the critical point.

Let us apply a similar argument to the results presented in Ref. P-S . Inspecting Fig. 3 of that work closely, contrary to what is claimed there, the distribution does not show a universal, broad form but rather gets more and more peaked as $`N`$ increases. Even the authors admit a very large drop of $`𝒫(D_2)`$ in the interval of $`1.5<D_2<2.0`$ that apparently is system size dependent. One may find a similar fixed point as in our Figures 1 and 2 at around $`D_2^{}1.55`$. This value is consistent with most of the presently known estimates of the correlation dimension calculated using different methods and listed in P-S . It is, however, somewhat larger then more recent and accurate estimates EC ; MEM . The possible reason for this discrepancy has been discussed in detail in MEM .

Finally we present the behavior of the correlation dimension as a function of the parameter $`b`$. This is presented in Fig. (5) which shows that in fact $`D_2`$ is a unique function of $`b\beta `$. In this figure we plotted the slope of $`D_2\mathrm{ln}N`$ vs $`\mathrm{ln}N`$ for the various values of $`b`$ and $`\beta `$. The analytical estimates of E-M

$$D_2=\{\begin{array}{cc}1(\pi b\beta )^1,\hfill & b1\hfill \\ 2b\beta ,\hfill & b1.\hfill \end{array}$$
(12)
According to the figure we see that our numerical simulations are in full accordance with the analytical estimates. The discrepancies here also are attributed to the considerably smaller set of samples and smaller range of matrix size $`N`$ as compared to E-M .

In conclusion we have performed numerical simulations of the power law band random matrices in order to obtain the PDF of the correlation dimension. We found that this function has a fixed point, $`D_2^{}`$, very close to the value of $`\stackrel{~}{D}_2`$ obtained from the scaling of the typical IPN, $`I_2^{typ}`$. Besides a direct relation between the average correlation dimension, $`D_2`$, enabled us to obtain a third estimate of the correlation dimension. All of these are in good agreement with each other that corroborates the expectation, that there is a unique value of the correlation dimension characterizing the states at the Anderson transition. For finite size simulations, the state-to-state fluctuations of both the IPN and the correlation dimension are large, nevertheless, both the PDF of $`I_2`$ and that of $`D_2`$ provide a unique value of $`\stackrel{~}{D}_2`$ in the limit $`N\mathrm{}`$. This way we have reconciled the apparent contradiction between the state-to-state fluctuation of the correlation dimension and the universal scaling of the IPN. Similar results may be achieved for other generalized dimensions, e.g. the information dimension, $`D_1`$ vi2 .

The author thanks financial support from the Alexander von Humboldt Foundation and from the OTKA Nos. T032116 and T034832.