text:
# Torsional oscillations in the solar convection zone

## 1 Introduction

An important feature of the solar convection zone is the presence of differential rotation in the form of a decrease in angular velocity from equator to the pole. This has been observed both in the surface layers (e.g. Snodgrass 1984) and deeper in the convection zone, as inferred from helioseismic measurements (e.g. Thompson et al. 1996). Furthermore, the differential rotation on the surface has been observed to vary with time (e.g. Howard & LaBonte 1980; Snodgrass, Howard, & Webster 1985). These so called torsional oscillations, which have periods of about 11 years, manifest themselves in the form of four alternating latitudinal bands of slightly faster and slower than average zonal flows which migrate towards the equator in about 22 years. These oscillations have also been confirmed by the analysis of the helioseismic data from the Solar and Heliospheric Observatory (SOHO) spacecraft for the present solar cycle (Kosovichev & Schou 1997; Schou et al. 1998).

Recent analysis of the helioseismic data, from both the Michelson Doppler Imager (MDI) instrument on board the SOHO spacecraft and the Global Oscillation Network Group (GONG) project, has also produced evidence that this banding signature is not merely a surface feature, but extends into the convection zone, to a depth of at least 8 percent in radius (Howe et al. 2000). These authors present data on departures of the reconstructed rotation rate from its temporal averages - the residuals - as a function of latitude at several target depths, which behave in a manner similar to the migration of sunspots during the solar cycle (the ‘butterfly diagram’). This finding is also supported by the analysis of Antia & Basu (2000), who use different data sets from GONG and independent inversion techniques. The time-base of these observations is only a few years, less than a complete solar cycle.

These torsional oscillations are thought to be produced as a consequence of the nonlinear interactions between the magnetic fields and the solar differential rotation. A zero order ‘mean’ differential rotation is assumed to be maintained by the Reynolds stresses of the turbulence; this can be included as a constant ‘background’ angular velocity, or explicitly parametrized, e.g. by the so called $`\mathrm{\Lambda }`$–effect (e.g. Rüdiger 1989). Attempts have been made to explain the surface oscillations in terms of the effects of the Lorentz force exerted by the large scale magnetic field on the azimuthal velocity field (e.g. Brandenburg & Tuominen 1988), or as a consequence of the ‘quenching’ of the turbulence-dependent quantities by the magnetic field (Kitchatinov et al. 1999).

Here we study this nonlinear coupling in the context of a two dimensional axisymmetric mean field dynamo model, in a spherical shell, in which the only nonlinearity is the action of the azimuthal component of the Lorentz force of the dynamo generated magnetic field on the solar angular velocity. Obtaining torsional oscillations in this way is also of interest in view of the fact that the Lorentz force is of second order in the magnetic field, thus naturally leading to the excitation of hydrodynamical oscillations of about half the period of the magnetic oscillations.

In the next section we introduce our model. Section 3 contains our results and section 4 gives a brief discussion.

## 2 The model

Here we shall assume that the gross features of the large scale solar magnetic field can be described by a mean field dynamo model, with the standard equation

$$\frac{𝐁}{t}=\times (𝐮\times 𝐁+\alpha 𝐁\eta \times 𝐁),$$
(1)
where $`𝐮=v\widehat{\varphi }\frac{1}{2}\eta `$, the term proportional to $`\eta `$ represents the effects of turbulent diamagnetism, and where the velocity field is taken to be of the form

$$v=v_0+v^{},$$
(2)
where $`v_0=\mathrm{\Omega }_0r\mathrm{sin}\theta `$, $`\mathrm{\Omega }_0`$ is a prescribed underlying rotation law and the component $`v^{}`$ satisfies

$$\frac{v^{}}{t}=\frac{(\times 𝐁)\times 𝐁}{\mu _0\rho r\mathrm{sin}\theta }.\widehat{\varphi }+\nu D^2v^{},$$
(3)
where $`D^2`$ is the operator $`\frac{^2}{r^2}+\frac{2}{r}\frac{}{r}+\frac{1}{r^2\mathrm{sin}\theta }(\frac{}{\theta }(\mathrm{sin}\theta \frac{}{\theta })\frac{1}{\mathrm{sin}\theta })`$ and $`\mu _0`$ is the induction constant. The assumption of axisymmetry allows the field $`𝐁`$ to be split into toroidal and poloidal parts, $`𝐁=𝐁_T+𝐁_P=B\widehat{\varphi }+\times A\widehat{\varphi }`$, and results in Eq. (1) yielding two scalar equations for $`A`$ and $`B`$. Nondimensionalizing in terms of the solar radius $`R_{}`$ and time $`R_{}^2/\eta _0`$, where $`\eta _0`$ is the maximum value of $`\eta `$ and letting $`\mathrm{\Omega }=\mathrm{\Omega }^{}\stackrel{~}{\mathrm{\Omega }}`$, $`\alpha =\alpha _0\stackrel{~}{\alpha }`$, $`\eta =\eta _0\stackrel{~}{\eta }`$, $`𝐁=B_0\stackrel{~}{𝐁}`$ and $`v^{}=\mathrm{\Omega }^{}R_{}\stackrel{~}{v}^{}`$, results in a system of equations for $`A,B`$ and $`v^{}`$, with the dynamo parameters $`R_\alpha =\alpha _0R_{}/\eta _0`$, $`R_\omega =\mathrm{\Omega }^{}R_{}^2/\eta _0`$, $`P_r=\nu _0/\eta _0`$, and $`\stackrel{~}{\eta }=\eta /\eta _0`$, where $`\mathrm{\Omega }^{}`$ is the solar surface equatorial angular velocity (see Moss & Brooke 2000 for details). Here $`\nu _0`$ and $`\eta _0`$ are the turbulent magnetic diffusivity and viscosity respectively and $`P_r`$ is the turbulent Prandtl number. The density $`\rho `$ is assumed to be uniform, and stress free boundary conditions ensure angular momentum conservation.

These equations were solved using the code and boundary conditions described in Moss & Brooke (2000), over the range $`r_0r1`$, $`0\theta \pi `$, with uniform spacing in both $`r`$ and $`\theta `$. The computational domain is the region $`r_0=0.64r1`$; with the solar convection zone proper being thought to occupy the region $`r>0.7`$, the region $`r_0r0.7`$ can be thought of as an overshoot region/tachocline. In the following simulations we used a mesh resolution of 61 points uniformly distributed radially and 101 points uniformly distributed latitudinally (over $`0\theta \pi `$), but test runs were carried out at higher spatial resolutions.

In this investigation, we took $`\mathrm{\Omega }_0`$ to be given in $`0.64r1`$ by an interpolation on the MDI data obtained from 1996 to 1999 (Howe et al. 2000), depicted in Fig. 1. For $`\alpha `$ we took $`\stackrel{~}{\alpha }=\alpha _r(r)f(\theta )`$, where $`f(\theta )=\mathrm{sin}^2\theta \mathrm{cos}\theta `$ (cf. Rüdiger & Brandenburg 1995) and

$$\alpha _r=1;0.7r0.8$$
(4)
with cubic interpolation to zero at $`r=r_0`$ and $`r=1`$, with the convention that $`\alpha _r>0`$ and $`R_\alpha <0`$. Also, in order to take into account the likely decrease in the turbulent diffusion coefficient $`\eta `$ in the overshoot region, we allowed a simple linear decrease from $`\stackrel{~}{\eta }=1`$ at $`r=0.8`$ to $`\stackrel{~}{\eta }=0.5`$ in $`r<0.7`$.

We monitor the time evolution of the total magnetic energy $`E`$ and the global parity of the magnetic field, defined as $`P=\frac{E^SE^A}{E^S+E^A}`$, where $`S`$ and $`A`$ refer respectively to the parts of the magnetic field that have symmetry or antisymmetry with respect to the equatorial plane (see also Brandenburg et al. 1989). Thus $`P=+1`$ and $`1`$ correspond to symmetric and antisymmetric fields respectively.

## 3 Results

We calibrated our model so that near marginal excitation the cycle period was about 22 years. This determined $`R_\omega =6\times 10^4`$, corresponding to $`\eta _02.5\times 10^{11}`$ cm2s-1, given the known values of $`\mathrm{\Omega }^{}`$ and $`R_{}`$. The first solutions to be excited in the linear theory are odd parity ($`P=1`$) limit cycles, which in this case have marginal dynamo number $`R_\alpha 3.12`$. The even parity ($`P=+1`$) solutions are also excited at similar marginal dynamo numbers of $`R_\alpha 3.16`$. We considered two values of the Prandtl numbers, $`P_r=0.1`$ and $`P_r=1.0`$.

With these parameter values, we found that this model, with underlying zero order angular velocity chosen to be consistent with the recent (MDI) helioseismic data (Fig. 1), is capable of producing butterfly diagrams which are in qualitative agreement with the observations. An example of this is depicted in Fig. 2. The polar feature is rather too strong – we have checked that this can be rectified by modifying slightly the spatial dependence of $`\alpha `$, by for example choosing $`f(\theta )=\mathrm{sin}^4\theta \mathrm{cos}\theta `$.

We also found that this model successfully produced torsional oscillations in the convection zone, similar to those deduced from recent helioseismic data. To compare our model with these results, we have plotted in Figs. 3 and 4 the variations of the rotation rate with latitude and time (‘butterfly diagrams’) to reveal the migrating banded zonal flows. In these models, the basic magnetic field parity is odd ($`P=1`$).

For the sake of comparison with observational results, we have also plotted in Fig. 5 the evolution of the residual rotation rate with time, at radius $`0.84`$ and latitude $`60`$ degrees.

As can be seen, consistent with the observations, in each hemisphere there are alternating latitudinal bands, with the width of approximately 10 degrees, of slightly faster and slower than average zonal flows. These migrate towards the equator in about 22 years, and extend deep into the convection zone. The amplitudes of these oscillations increase with depth below the surface and depend on the parameters of our model, in particular the Prandtl number. For $`P_r=1.0`$, these amplitudes range from about 0.07 nHz at the surface to more than $`0.4`$ nHz towards the bottom of the convection zone, somewhat lower than, but in principle compatible with, the results of Howe et al. 2000. The torsional oscillations present in our model have periods half that of the period of the global magnetic field which is compatible with the observed period of the oscillations at the surface and consistent with the observed behaviour inside the convection zone.

The strictly odd parity models presented have an equatorial feature in the ‘butterfly diagrams’ for the velocity perturbations which does not appear to be present in the current inversions of the observational data. This becomes weaker as Prandtl number increases. We also note in passing that these figures are almost identical for cases where even parity solutions are found, except that this feature is then absent. We further point out that the large scale solar magnetic field is probably of mixed global parity (predominantly odd) – see, e.g., Pulkinnen et al. (1999).

## 4 Discussion

We have studied a solar dynamo model, calibrated to have the correct cycle period, with a mean rotation law given by recent helioseismic observations. This model produces butterfly diagrams in qualitative agreements with those observed.

We have shown that a nonlinear coupling between the magnetic field and the solar differential rotation, where the nonlinearity is due to the action of the azimuthal component of the Lorentz force, is capable of producing torsional oscillations, with period of about 11 years, which penetrate into the convection zone and which migrate towards the equator in about 22 years. The period of these oscillations is about half that of the period of the global magnetic fields. This is in agreement with the observed period of the torsional oscillations at the surface. For the oscillations inside the convection zone, this is a testable prediction, not contradicted by the current helioseismic observations, which so far extend over an interval less than a solar cycle. We note also that solutions with even parity ($`P=+1`$) show slightly larger amplitudes without the equatorial feature.

The current inversions of the helioseismological data seem to suggest that the torsional oscillations largely disappear below about $`R=0.9R_{}`$, in contrast to our model oscillations. However, there are uncertainties in these inversions, specially at the deeper levels. At the same time our dynamical model is oversimplified and substantial improvements can be made. In particular, the predicted amplitudes for torsional oscillations in our model are likely to be affected by our assumption of uniform density in Eq. (3). Nevertheless, we find it interesting that some of the major features in the torsional oscillations can be readily reproduced. A more detailed study of these oscillations, including more realistic density profiles, is in progress.

###### Acknowledgements.

We would like to thank Rachel Howe and Mike Thompson for many useful discussions regarding the helioseismic observations and for providing us with the data displayed in Fig. 1. EC is supported by a PPARC fellowship and RT benefited from PPARC UK Grant No. L39094.


format: text

text:
# Phase separation in a chaotic flow

\[

## Abstract

The phase separation between two immiscible liquids advected by a bidimensional velocity field is investigated numerically by solving the corresponding Cahn-Hilliard equation. We study how the spinodal decomposition process depends on the presence — or absence — of Lagrangian chaos. A fully chaotic flow, in particular, limits the growth of domains and for unequal volume fractions of the liquids, a characteristic exponential distribution of droplet sizes is obtained. The limiting domain size results from a balance between chaotic mixing and spinodal decomposition, measured in terms of Lyapunov exponent and diffusivity constant, respectively.

\]

A system of two immiscible fluids at rest will gradually phase-separate, forming domains whose size grows algebraically with time. Everyday experience, however, shows that by continuously stirring or shaking the fluids the domains or droplets of the phases (say, oil and vinegar) break and coalesce, leading to a dynamic stationary state with domains of finite size.

A first approach consists in modelling this situation by subjecting the binary fluid to a homogeneous shear velocity field . However, experiments , numerical simulations , and more recently analytical approaches  show that in such a situation infinitely long domains aligned with the flow are formed. The effect of the velocity field is to counter the Rayleigh instability, stabilising lamellar and (in certain cases) even cylindrical domains . Domain breakup in those situations takes place only at large Reynolds numbers, and is generally attributed to inertial effects . Studying these inertial effects numerically is difficult, as a realistic description of the feedback of domain shape on the flow is required . The corresponding simulations are therefore limited by finite size effects.

In this paper we investigate a different mechanism by which domains of finite size can be stabilized in a demixing system. In particular, we show that a saturation of the average length scale takes place even in the absence of inertial effects if the flow has Lagrangian chaos, (i.e. if the trajectories of nearby starting points diverge exponentially with time). This is interesting for two reasons: First, with an appropriate time-dependence of the velocity field one can still have Lagrangian chaos in a situation of high viscosity in which inertial effects are negligible — this is how one mixes pastes. Secondly, it is possible in that case to decouple the hydrodynamic problem from the phase separation. This problem of a passive, phase separating scalar field  is of course much simpler, so that simulations using large systems are possible. Our approach therefore extends earlier extensive studies of passive scalar advection by periodically driven chaotic flows .

Our study is also related in spirit to earlier studies of advection by ‘synthetic’ velocity fields tuned to model turbulent flows . Phase separation was studied in this context in Ref. . An essential difference, compared to our work, is that in such turbulent flows the separation between nearby tracer particles appears to increase algebraically, rather than exponentially, with time.

We consider a two dimensional flow that can be tuned to be regular, mixed, or fully chaotic. Specifically, the incompressible velocity field $`𝒗(x,y,t)`$ is a modified version of the so-called time-dependent Harper map  (related to the ‘partitioned-pipe mixer’, a special case of ‘eggbeater flow’ ). The dynamics takes place on a square of side $`L`$, with periodic boundary conditions. The velocity field is an alternating sequence of shears in the $`x`$ and in the $`y`$ direction with a time period $`T`$,

$`v_x`$ $`={\displaystyle \frac{2\pi \alpha L}{T}}\mathrm{sin}\left({\displaystyle \frac{2\pi y}{L}}\right);v_y=0;`$ $`n<{\displaystyle \frac{t}{T}}<n+{\displaystyle \frac{1}{2}}`$ (1)
$`v_y`$ $`={\displaystyle \frac{2\pi \alpha L}{T}}\mathrm{sin}\left({\displaystyle \frac{2\pi x}{L}}\right);v_x=0;`$ $`n+`$ $`{\displaystyle \frac{1}{2}}<{\displaystyle \frac{t}{T}}<n+1.`$
The parameter $`\alpha `$ controls the chaoticity of the trajectories. If $`\alpha `$ is small, the two semicycles are composed into the smooth, laminar velocity field: $`v_x=\frac{\pi \alpha L}{T}\mathrm{sin}\left(\frac{2\pi y}{L}\right)`$; $`v_y=\frac{\pi \alpha L}{T}\mathrm{sin}\left(\frac{2\pi x}{L}\right)`$. For larger values of $`\alpha `$ the trajectories stretch and fold, and the flow becomes chaotic. In order to visualise this, it is convenient to follow the position of a point at the end of each cycle. This ‘kicked Harper’ map is shown in Fig. 1 for several values of $`\alpha `$. For $`\alpha 0.2`$ the flow is a mixture of laminar and chaotic regions, and becomes fully chaotic around $`\alpha 0.4`$. In the chaotic situation, it is convenient to characterise the flow by the Lyapunov exponent $`\lambda `$, defined by the fact that nearby starting points separate as $`e^{\lambda t}`$. We have computed $`\lambda `$ as in Ref.  and found that the relation $`\lambda 1.96\mathrm{ln}(3.35\alpha )/T`$ is a good approximation throughout the chaotic regime, $`\alpha 0.4`$.

The spinodal decomposition of the two-component fluid is described by the Cahn-Hilliard equation

$$\frac{\varphi (𝒓,t)}{t}+𝒗(𝒓,t)\mathbf{}\varphi (𝒓,t)=\mathrm{\Gamma }^2\left(\frac{\delta F[\varphi ]}{\delta \varphi (𝒓,t)}\right).$$
(2)
Here $`\varphi `$ is a dimensionless concentration field, the concentrations of the species are $`[1\pm \varphi ]/2`$. We work at $`T=0`$, since temperature is irrelevant in this process . The free energy functional is of the Ginzburg-Landau form and reads

$$F[\varphi ]=\mathrm{d}^d𝒙\left[\frac{\xi ^2}{2}(\mathbf{}\varphi )^2\frac{1}{2}\varphi ^2+\frac{1}{4}\varphi ^4\right].$$
(3)
Here, $`\xi `$ is the equilibrium correlation length controlling the width of the interfaces, and $`d`$ is the number of spatial dimensions. We consider two topologically different situations (i) $`\varphi 0`$ : a species is less abundant than the other and forms disconnected droplets, and (ii) $`\varphi =0`$: the two phases are in equal quantity and form a bicontinuous structure. Situation(i) has been studied experimentally .

In a chaotic flow, the passive scalar mixes rapidly, whereas in the case of phase separation this tendency is opposed by surface tension. The competition between these two effects can be quantified through two adimensional parameters, $`D\mathrm{\Gamma }T/\xi ^2`$ (the adimensional transport coefficient of the Cahn-Hilliard equation) and the chaoticity parameter $`\alpha `$, or alternatively the adimensional Lyapounov exponent $`\lambda T`$. A large $`D`$ means that appreciable diffusive transport will take place during each laminar half cycle. A large $`\lambda T`$, on the other hand, means that the mixing process is efficient within a few cycles. Note that $`\lambda `$ can also be interpreted as an average elongation or shear rate experienced by the fluid particles.

Equation (2) is integrated numerically with the velocity field (1), using the implicit spectral method developed and discussed in Ref. . The results are presented with time and length units chosen as the cycle period $`T`$ and the interfacial thickness $`\xi `$, respectively. The system size, lattice parameter and time step are $`L=512\xi `$, $`\mathrm{\Delta }x=\xi `$, and $`\mathrm{\Delta }t=510^4T`$ respectively.

Existence of a stationary state. — We first show that a purely chaotic flow does indeed stop the domain growth. In Fig. 2, we show the evolution of a phase-separated sample with $`\varphi =1/2`$, upon turning on a chaotic velocity field ($`\alpha =0.4`$). The large droplets of the initial configuration are broken into smaller droplets, until a stationary state where droplets successively grow and break is reached. Fig. 3 shows the late stages of coarsening of a system with equal concentrations of phases ($`\varphi =0`$) in the four velocity fields of Fig. 1. For $`\alpha =0.1`$, the velocity field is laminar. We observe in that case structures very similar to those found in a homogeneous shear flow, but which now follow the winding flow lines. In the mixed case, $`\alpha =0.25`$, large-domain structures form in the laminar regions of the flow, and break into very small domains in the chaotic ones. In the fully chaotic situation, $`\alpha 0.4`$, a dynamical stationary state is reached, with small domains continuously breaking and reforming. For $`\alpha =1.0`$, the sinusoidal nature of the underlying velocity field becomes apparent. This snapshot nicely illustrates the typical ‘stretch and fold’ processes characteristic of chaotic advection .

Scaling properties in the chaotic flow. — In the stationary regime, it is clear from Figs. 2 and 3 that there exist a typical length scale $`L^{}`$ which depends on the parameters $`D`$ and $`\lambda `$: this will be confirmed below by a quantitative analysis. As in the pure coarsening case , scaling properties are expected in the regime $`\xi L^{}(D,\lambda )L`$. The length scale $`L^{}`$ may be estimated by the following simple argument. In the absence of flow, the domains grow as $`L(t)/\xi (Dt/T)^{1/3}`$, and this growth is stopped by the chaotic flow which introduces a time scale $`\lambda ^1`$. Hence, we estimate $`L^{}L(t=\lambda ^1)`$ and predict

$$L^{}(D,\lambda )\xi \left(\frac{D}{\lambda T}\right)^{1/3}.$$
(4)
Isolated droplets, $`\varphi 0`$. — Following Ref. , we characterize the assembly of droplets by computing the distribution of droplet surfaces $`f(S)`$, where $`f(S)\mathrm{d}S`$ is the probability that the surface occupied by a droplet is between $`S`$ and $`S+\mathrm{d}S`$. In the scaling regime, we expect this distribution to be of the form

$$f(S)\frac{1}{S^{}}\left(\frac{S}{S^{}}\right),$$
(5)
where $`S^{}`$ is a typical droplet area. In Fig. 4, data obtained for a wide range for the values of $`D`$ and $`\lambda `$ are collapsed by using a reduced variable $`S/S^{}`$, with $`S^{}(D/\lambda )^{0.62}`$. The data collapse is satisfactory, and the result for $`S^{}`$ reasonably close to what would be expected from Eq. (4), i.e. $`S^{}(D/\lambda )^{2/3}`$. Finding an exponent slightly smaller than the one expected theoretically is not surprising, since the typical domain sizes are rather small ($`S^{}50\xi ^2`$), so that the asymptotic value for the domain growth exponent in the absence of flow may not be reached. The rescaled distribution functions exhibit an exponential tail, $`(y)e^y`$ (dotted-dashed line in Fig. 4). Such distributions are very is similar to those found in the experiments of Ref. . For the largest droplets, deviations from the exponential fit are observed, indicating either insufficient statistics or a different scaling behaviour for the extreme values of $`S`$ .

Equal concentrations: $`\varphi =0`$. — In the case of equal concentrations, the domains are ramified and extend throughout the sample, so that the area is not a useful measure of domain size. A characteristic domain size can nevertheless be obtained from the two-point correlation function $`C(𝒓,t)L^2\mathrm{d}^2𝒙\varphi (𝒙,t)\varphi (𝒙\mathbf{+}𝒓,t)`$, which is the Fourier transform of the structure factor measured in light scattering experiments. Performing a time average over many configurations shows that the bicontinuous structure is on average perfectly isotropic, as it is in the absence of flow. One can therefore average $`C(𝒓,t)`$ over orientations to obtain a one variable function, $`C(r)`$. The characteristic domain size $`L^{}`$ can be defined by $`C(L^{})=0.5`$. Fig. 5 displays this domain size for various combinations of $`D`$ and $`\lambda `$, as a function of the ratio $`D/\lambda `$. The data can be fitted by $`L^{}\xi \left(D/\lambda T\right)^{0.27}`$ Again, this is in reasonable agreement with the scaling analysis, Eq. (4). Larger simulations, with smaller values of $`\lambda `$, would be necessary to obtain larger domain sizes and avoid the crossover effects which are well known in spinodal decomposition simulations .

A more detailed information on the domain structure is obtained from the full correlation function $`C(r)`$. Here, one could expect from the scaling hypothesis a behaviour of the form

$$C(r)𝒞\left(\frac{r}{L^{}}\right)$$
(6)
with $`𝒞`$ a universal function. This hypothesis is tested in Fig. 6, where $`C(r)`$ is represented for fixed $`\lambda `$ and various values of $`D`$. At fixed $`\lambda `$, a good collapse of the data obtained for different $`D`$ is achieved by using a rescaled variable $`r/D^{0.27}`$. The inset of Fig. 6, however, shows that the shape of the scaling function slightly depends on $`\lambda `$, so that the universal scaling expressed by Eq. (6) is not valid. We attribute this change of the scaling function with the flow pattern to the fact that even in the chaotic regime the flow cannot be considered as being homogeneous and isotropic, but exhibits an underlying sinusoidal structure. This is in contrast with the case of isolated droplets where the droplet distribution was not affected by this structure, recall Fig. 4.

We have studied the phase separation in conditions in which the species boundaries are passively advected by a incompressible flow. We have shown that a chaotic flow results in a steady state with domains of finite size resulting from the balance between spinodal decomposition and chaotic advection, Eq. (4). This should be contrasted with the situation observed in turbulent flow, where the flow intensity must exceed a threshold in order to stop domain growth . Such a difference can be traced back to the fact that Lyapounov exponents for passive scalar advection are actually 0 in the latter case. The essential approximation in our work, compared to realistic experimental situations, is the assumption that the flow pattern is not modified by the domain growth. This assumption, however, may not be unrealistic if the two fluids have similar viscosities and if the capillary stresses are small compared to viscous stresses. This is measured by the capillary number $`C_a=\eta \lambda /(\gamma /L^{})`$, where $`\eta `$ is the viscosity and $`\gamma `$ the surface tension. In highly viscous fluids, $`C_a`$ is expected to be large, so that the decoupling is possible. This decoupling also makes it possible to consider analytical treatments.

We acknowledge useful discussions with A. J. Bray, B. Cabane, P. Leboeuf, J. F. Pinton and J. E. Wesfreid.


text:
# Epidemic spreading in scale-free networks

## Abstract

The Internet, as well as many other networks, has a very complex connectivity recently modeled by the class of scale-free networks. This feature, which appears to be very efficient for a communications network, favors at the same time the spreading of computer viruses. We analyze real data from computer virus infections and find the average lifetime and prevalence of viral strains on the Internet. We define a dynamical model for the spreading of infections on scale-free networks, finding the absence of an epidemic threshold and its associated critical behavior. This new epidemiological framework rationalize data of computer viruses and could help in the understanding of other spreading phenomena on communication and social networks.

Many social, biological, and communication systems can be properly described by complex networks whose nodes represent individuals or organizations, and links mimic the interactions among them. Particularly interesting examples are the Internet and the world-wide-web, which have been extensively studied because of their technological and economical relevance. These studies have revealed, among other facts, the scale-free nature of these networks . This results in the power-law distribution $`P(k)k^\gamma `$ for the probability that a node of the network has $`k`$ connections to other nodes, with an exponent $`\gamma `$ that ranges between $`2`$ and $`3`$. The importance of local clustering is indeed the key ingredient in the modeling of these networks with the recent introduction of scale-free (SF) graphs .

In view of the wide occurrence of complex networks in nature it is of great interest to inspect the effect of their features on epidemic and disease spreading , and more in general in the context of the nonequilibrium phase transitions typical of these phenomena . The study of epidemics on these networks finds an immediate practical application in the understanding of computer virus spreading , and could also be relevant to the fields of epidemiology  and pollution control .

In this Letter, we analyze data from real computer virus epidemics, providing a statistical characterization that points out the importance of incorporating the peculiar topology of scale-free networks in the theoretical description of these infections. With this aim, we study by large scale simulations and analytical methods the susceptible-infected-susceptible  model on SF graphs. We find the absence of an epidemic threshold and its associated critical behavior, which implies that SF networks are prone to the spreading and the persistence of infections whatever spreading rate the epidemic agents possess. The absence of the epidemic threshold—a standard element in mathematical epidemiology–radically changes many of the standard conclusions drawn in epidemic modeling. The present results are also relevant in the field of absorbing-state phase transitions and catalytic reactions.

The epidemiological analysis of computer viruses has been the subject of a continuous interest in the computer science community , following mainly approaches borrowed from biological epidemiology . The standard model used in the study of computer virus infections is the susceptible-infected-susceptible (SIS) epidemiological model. Each node of the network represents an individual and each link is a connection along which the infection can spread to other systems. This model relies on a coarse grained description of individuals in the population. Individuals exist only in two discrete states, “healthy” or “infected”. At each time step, each susceptible (healthy) node is infected with rate $`\nu `$ if it is connected to one or more infected nodes. At the same time, infected nodes are cured and become again susceptible with rate $`\delta `$, defining an effective spreading rate $`\lambda =\nu /\delta `$. Without lack of generality, we can set $`\delta =1`$. The updating can be performed both with parallel and sequential dynamics. In models with local connectivity (Euclidean lattices and mean-field models), the most significant result is the general prediction of a nonzero epidemic threshold $`\lambda _c`$ . If the value of $`\lambda `$ is above the threshold, $`\lambda \lambda _c`$, the infection spreads and becomes persistent. Below it, $`\lambda <\lambda _c`$, the infection dies out exponentially fast. The epidemic threshold is actually equivalent to a critical point in a nonequilibrium phase transition. In this case, the critical point separates an active phase with a stationary density of infected nodes from a phase with only healthy nodes and null activity. In particular, it is easy to recognize that the SIS model is a generalization of the contact process (CP) model, that has been extensively studied in the context of absorbing-state phase transitions . Statistical observations of virus incidents in the wild, on the other hand, indicate that all viruses that are able to pervade, spread much slower than exponentially, and saturate to a very low level of persistence, affecting just a tiny fraction of the total number of computers . This fact is in striking contradiction with the theoretical predictions unless in the very unlikely chance that all computer viruses have an effective spreading rate tuned just infinitesimally above the threshold. This points out that the view obtained so far with the modeling of computer virus epidemics is very instructive but not completely adequate to represent the real phenomenon.

In order to gain further insight on the spreading properties of viruses in the wild, we have analyzed the prevalence data reported by the Virus Bulletin  from February 1996 to March 2000. We have analyzed in particular the surviving probability of homogeneous groups of viruses, classified according to their infection mechanism . We consider the total number of viruses of a given strain that born and die within our observation window. Hence, we calculate the surviving probability $`P_s(t)`$ of the strain as the ratio of viruses still alive at time $`t`$ after their birth and the total number of observed viruses. Fig. 1 shows that the surviving probability suffers a sharp drop in the first two months of a virus’ life. This is a well-known feature  indicating that statistically only a small percentage of viruses give rise to a significant outbreak in the computer community. Fig. 1, on the other hand, shows for larger times a clean exponential tail, $`P_s(t)\mathrm{exp}(t/\tau )`$, where $`\tau `$ represents the characteristic life-time of the virus strain . The numerical fit of the data yields $`\tau 14`$ months for boot and macro viruses and $`\tau 69`$ months for file viruses. These characteristic times are impressively large if compared with the interval in which anti-virus software is available on the market (usually within days or weeks after the first incident report) and indicate that the viral persistence time scale is more related to the implementation of prophylactic safety measures than to the timely availability of the specific anti-virus. These external factors, however, are not possibly competing on the short time scale of the viruses spread (days or weeks), and again we face the very puzzling question of why viruses seem to have access to persistent low prevalence levels but never grow exponentially.

The key point to understand the puzzling properties exhibited by computer viruses resides in the capacity of many of them to be borne by electronic mail as an apparently innocuous attachment . Having this property in mind, it is easy to realize that the topology of the connections between individuals cannot be correctly represented by an Euclidean lattice, or a mean-field model. In this sense, these connections should instead have essentially the topology of the Internet, through which electronic mail travels. The scale-free connectivity of the Internet implies that each node has a statistically significant probability of having a very large number of connections compared to the average connectivity $`k`$ of the network. That opposes to conventional random networks (local or nonlocal) in which each node has approximately the same number of links $`kk`$ . The fact that all virus strains qualitatively show the same statistical features indicates that very likely all of them spread on networks with connectivity properties analogous to those of the Internet . It is then natural to foresee that scale-free properties should be included in a theory of epidemic spreading of computer viruses.

To address the effects of scale-free connectivity in epidemic spreading we study the SIS model on the SF network. We consider the graph generated by using the algorithm devised in Refs. : We start from a small number $`m_0`$ of disconnected nodes; every time step a new vertex is added, with $`m`$ links that are connected to an old node $`i`$ with probability $`k_i/_jk_j`$. After iterating this scheme a sufficient number of times, we obtain a network composed by $`N`$ nodes with connectivity distribution $`P(k)k^3`$ and average connectivity $`k=2m`$. In this work we take $`m=3`$. We have performed numerical simulations on graphs with number of nodes ranging from $`N=10^3`$ to $`N=8.5\times 10^6`$ and studied the variation in time and the stationary properties of the density of infected nodes $`\rho `$ in surviving infections; i.e. the virus prevalence. Initially we infect half of the nodes in the network, and iterate the rules of the SIS model with parallel updating. After an initial transient regime, the system stabilizes in a steady state with a constant average density of infected nodes. The prevalence is computed averaging over at least $`100`$ different starting configurations, performed on at least $`10`$ different realizations of the random networks.

The first arresting evidence from simulations is the absence of an epidemic threshold, i.e., $`\lambda _c=0`$. In Fig. 2 we show the virus prevalence in the steady state that decays with decreasing $`\lambda `$ as $`\rho \mathrm{exp}(C/\lambda )`$, where $`C`$ is a constant. This implies that for any finite value of $`\lambda `$ the virus can pervade the system with a finite prevalence, in sufficiently large networks. In all networks with bounded connectivity the steady state prevalence is always null below the epidemic threshold; i.e. all infections die out. Further evidence to our results is given by the total absence of scaling of $`\rho `$ with the number of nodes that is, on the contrary, typical of epidemic transitions in the proximity of a finite threshold . This allows us to exclude the presence of any spurious results due to network finite size effects. The present result can be intuitively understood by noticing that for usual lattices, the higher the node’s connectivity, the smaller the epidemic threshold. In a SF network the unbounded fluctuations in connectivity ($`k^2=\mathrm{}`$) plays the role of an infinite connectivity, annulling thus the threshold.

Finally, we analyze the spreading of infections starting from a localized virus source. We observe that the spreading growth in time has an algebraic form, that is in agreement with real data that never found an exponential increase of a virus in the wild. Noteworthy, by applying the definition of surviving probability $`P_s(t)`$ used to analyze real data, we recover in our model the same exponential behavior in time (see Fig. 3a). The characteristic life-time depends on the spreading rate and the network sizes, allowing us to relate the average lifetime of a viral strain with an effective spreading rate and the Internet size. At the same time, the divergence of lifetimes for larger networks points out the possible increasing of the viruses lifetime during the eventual expansion of the Internet.

We can also approach the system analytically by writing the mean-field equation governing the time evolution of $`\rho (t)`$. In order to take into account connectivity fluctuations, we consider the relative density $`\rho _k(t)`$ of infected nodes with given connectivity $`k`$; i.e the probability that a node with $`k`$ links is infected. The dynamical mean-field reaction rate equations can be written as

$$_t\rho _k(t)=\rho _k(t)+\lambda k(1\rho _k(t))\mathrm{\Theta }(\lambda ).$$
(1)
The creation term considers the probability that a node with $`k`$ links is healthy $`(1\rho _k(t))`$ and gets the infection via a connected node. The probability of this event is proportional to the infection rate, the number of connections, and the probability $`\mathrm{\Theta }(\lambda )`$ that any given link points to an infected node. By imposing stationarity ($`_t\rho _k(t)=0`$) we find the stationary densities

$$\rho _k=\frac{k\lambda \mathrm{\Theta }(\lambda )}{1+k\lambda \mathrm{\Theta }(\lambda )},$$
(2)
denoting that the higher the node connectivity, the higher the probability to be infected. This inhomogeneity must be taken into account in the self-consistent calculation of $`\mathrm{\Theta }(\lambda )`$. Indeed, the probability that a link points to a node with $`s`$ links is proportional to $`sP(s)`$. In other words, a randomly chosen link is more likely to be connected to a node with high connectivity, yielding

$$\mathrm{\Theta }(\lambda )=\underset{k}{}\frac{kP(k)\rho _k}{_ssP(s)}.$$
(3)
Since $`\rho _k`$ is on its turn function of $`\mathrm{\Theta }(\lambda )`$, we obtain a consistency equation that allows to to find $`\mathrm{\Theta }(\lambda )`$. Finally we can evaluate the behavior of $`\rho `$ by solving the second consistency relation

$$\rho =\underset{k}{}P(k)\rho _k,$$
(4)
that expresses the average density of infected nodes in the system. In the SF model considered here, we have a connectivity distribution $`P(k)=2m^2/k^3`$, where $`k`$ as approximated as a continuous variable . In this case, integration of Eq.(3) allows to write $`\mathrm{\Theta }(\lambda )`$ as

$$\mathrm{\Theta }(\lambda )=\frac{e^{1/m\lambda }}{\lambda m}(1e^{1/m\lambda })^1,$$
(5)
from which, using Eq.(4), we find at lowest order in $`\lambda `$:

$$\rho =2e^{1/m\lambda }+h.o..$$
(6)
This very intuitive calculation recovers the numerical findings and confirms the surprising absence of any epidemic threshold or critical point in the model; i.e $`\lambda _c=0`$. Finally, as a further check of our analytical results, we have numerically computed in our model the relative densities $`\rho _k`$ , recovering the predicted dependence upon $`k`$ of Eq.(2) (see Fig. 3b). It is also worth remarking that the present framework can be generalized to networks with $`2<\gamma 3`$, recovering qualitatively the same results .

The emerging picture for epidemic spreading in complex networks emphasizes the role of topology in epidemic modeling. In particular, the absence of epidemic threshold and critical behavior in a wide range of scale-free network provide an unexpected result that changes radically many standard conclusions on epidemic spreading. This indicates that infections can proliferate on these scale-free networks whatever spreading rates they may have. These very bad news are, however, balanced by the exponentially small prevalence for a wide range of spreading rates ($`\lambda <<1`$). This point appears to be particularly relevant in the case of technological networks such as the Internet and the world-wide-web that show scale-free connectivity with exponents $`\gamma 2.5`$ . For instance, the present picture fits perfectly with the observation from real data of computer virus spreading, and could solve the long standing problem of the generalized low prevalence of computer viruses without assuming any global tuning of the spreading rates. The peculiar properties of scale-free networks also open the path to many other questions concerning the effect of immunity and other modifications of epidemic models. As well, the critical properties of many nonequilibrium systems could be affected by the topology of scale-free networks. Given the wide context in which scale-free networks appears, the results obtained here and the proposed investigations could have intriguing implications in biology and social systems.

This work has been partially supported by the European Network Contract No. ERBFMRXCT980183. RP-S also acknowledges support from the grant CICYT PB97-0693. We thank S. Franz, M.-C. Miguel, R. V. Solé, M. Vergassola, S. Zapperi and R. Zecchina for helpful comments and discussions.