text:
# Anomalous scaling and universality in hydrodynamic systems with power-law forcing

## 1 Introduction

The understanding of the small scale statistics of turbulent systems is a problem of considerable interest . With turbulent systems, we mean both the dynamics of velocity fields in high-Reynolds number flows, and the advection of scalar or vector fields by turbulent flows as, e.g., temperature or magnetic fields. In the last few years, much progress has been reached both in experimental  and numerical  investigations of turbulent systems. We have now plenty of observations showing that two properties generally holds. First, turbulent systems are intermittent , as quantified by the anomalous scaling of moments of field increments. Second, the anomalous scaling exponents display universality with respect to the boundary conditions and to the large scale forcing mechanisms . Some subtle points arise when isotropy is broken at large scale by the injection mechanisms . In that case, for universality to hold it is required that anisotropic contributions are subleading with respect to the isotropic one. However, evidences for universality of both isotropic and anisotropic scaling exponents have been found .

Nonetheless, the mechanisms responsible for anomalous scaling and universality in turbulent systems are still not fully comprehended with the remarkable exception of linear problems such as the passive transport of a field, for which intermittency and universality have been systematically understood (see Ref.  for an exhaustive review). In this context, in particular for the class of Kraichnan models , closed equations for the correlation functions can be derived. These are linear partial differential equations, whose homogeneous solutions (zero modes) generally exhibit anomalous scaling. On the other hand, the inhomogeneous solutions, constrained by the external forcing, possess dimensional (non-anomalous) scaling. Universality results then from the decoupling between the zero modes scaling and the forcing. Remarkably, the zero modes can be interpreted as statistically preserved structures, i.e. functions that do not change in time once averaged over the velocity field realisations and particle trajectories. This allowed for successfully testing the entire picture also for the advection by realistic velocity fields , and in the context of shell models for passive transport .

On the theoretical side, it is tempting to export the concepts issuing from linear turbulent systems to nonlinear ones, i.e. to see whether the same mechanisms for anomalous scaling and universality are at work also in nonlinear hydrodynamic systems such as Navier-Stokes turbulence. A possible test case to probe universality with respect to the forcing mechanisms is to study turbulent fluctuations stirred at all scales by a self-similar power-law random forcing. Indeed, the non-analytical properties of the forcing may alter the energy exchange between the turbulent field fluctuations. For instance, in some cases, the energy injection mechanism may prevail over the energy cascade process, modifying the inertial range physics.

This problem was pioneered by means of Renormalisation Group (RG) methods , in the late ’70s, for the $`d`$-dimensional Navier-Stokes equations. In the RG calculations, for $`d=3`$, the forcing spectrum is chosen as $`E_f(k)k^{3y}`$ with $`y`$ playing the role of a small parameter in perturbative expansions. Unfortunately the interesting physical case, obtained for $`y=4`$ and corresponding to the Kolmogorov spectrum for the velocity field, lies in a range where convergence of the RG expansion is not granted . Notwithstanding extensions of the RG formalism to $`yO(1)`$ values have been tempted by different approaches , the problem is still open. Recent numerical simulations tried to shed light on this issue  but, because of the limited resolution, their results are not conclusive.

The aim of this paper is to investigate the general issue of the small scale statistical properties of linear and nonlinear hydro-dynamical systems, in the presence of stirring acting at all scales with a power-law spectrum. A systematic study of the scaling behaviour at varying the forcing spectrum allows indeed for understanding the interplay of dimensional and anomalous scaling in turbulent fields.

In Sec. 2, as an instance of linear problems, we consider passive scalars stirred at all scales by a power-law forcing, and advected by a $`2d`$ turbulent velocity field. Our results confirm the universality scenario originating from the zero modes picture which predicts two distinct regimes. Forcing-dominated regime: the scaling of low order structure functions is non-anomalous, with exponents dimensionally related to the forcing spectrum; for the higher order moments, scaling is anomalous and dominated by the zero modes. Forcing-subleading regime: the dimensional scaling related to the balance with the forcing is subleading, at any order, with respect to the anomalous one, similarly to the case of a standard large scale injection. Hence, anomalous scaling is observed for any order statistics. A subtle technical point revealed by the passive scalar case is the existence of many important power-law terms that contribute to the scaling properties. As clarified in Sec. 2, for both regimes, to disentangle the authentic scaling behaviours, it is necessary to take into account the leading as well as the subleading terms. Concerning nonlinear systems, in Sec. 3, we review the numerical study of the $`3d`$ Navier-Stokes equations done in Ref. , where due to natural limitation in the resolution only semi-quantitative results were obtained. Then, to overcome this difficulty, we consider shell models for turbulence. These are nonlinear models that maintain most of the richness of the original problem but allow for reaching higher Reynolds numbers. Shell model results coherently fit with the passive scalar and NS ones. Conclusions follow in Sec. 4.

## 2 Linear dynamics: passive scalar transport

The evolution of a passive scalar field $`\theta `$ advected by an incompressible flow $`𝒗`$ is governed by the advection-diffusion equation :

$$_t\theta +𝒗\mathbf{}\theta =\kappa \mathrm{\Delta }\theta +f,$$
(1)
where $`\kappa `$ is the molecular diffusivity. The standard phenomenology is as follows. Scalar fluctuations are injected at large scale by the source term $`f`$; then, by a cascade mechanism induced by the advection term, they reach small scales where dissipation takes place balancing the input, and holding the system in a statistically stationary state.

Since the problem is passive, it can be studied also with a synthetic velocity field mimicking some of the features of realistic turbulent flows. Much attention has been recently devoted to the Kraichnan model , where $`𝒗`$ is a Gaussian, incompressible, homogeneous, isotropic and $`\delta `$-correlated in time field of zero mean. The only reminiscence of real turbulent flows is the existence of a scaling behaviour that, due to the assumed self-similarity, is fixed by a unique exponent, i.e. $`(\delta _r𝒗\widehat{𝒓})^2r^\xi `$ with $`0<\xi <2`$ (where $`\delta _r𝒗=𝒗(𝒙+𝒓)𝒗(𝒙)`$). For the sake of analytical control, the forcing $`f`$ is also taken as a Gaussian, homogeneous, isotropic field of zero mean and with correlation function: $`f(𝒓,t)f(\mathrm{𝟎},0)=\delta (t)(r)`$, where $`(r)`$ decays rapidly for $`rL_f`$, identifying $`L_f`$ as the large scale of the problem.

Thanks to the above simplifying assumptions, a closed equation for the generic $`p`$-point correlator, $`C_p(𝒓,t)=\theta (𝒓_1,t)\theta (𝒓_2,t)\mathrm{}\theta (𝒓_p,t)`$, can be derived. It formally reads

$$_tC_p=_pC_p+C_{p2},$$
(2)
where $``$ is the forcing spatial correlation defined above, and $`_p`$ is a linear differential operator coming from the diffusion and the advection terms (see Ref.  for details). The stationary solution of (2) satisfies the equation $`_pC_p=C_{p2}`$. As usual, a linear equation is solved by the superposition of the solution of the associated homogeneous equation, $`_p𝒵_p=0`$, and the solution of the inhomogeneous one, which will be denoted as $`C_p^{}`$. In the inertial range of scales both $`C_p^{}`$ and $`𝒵_p`$ display a scaling behaviour, i.e.

$$C_p^{}(\lambda \underset{¯}{𝒓})\lambda ^{\zeta _p^{dim}}C_p^{}(\underset{¯}{𝒓})\text{and}𝒵_p(\lambda \underset{¯}{𝒓})\lambda ^{\zeta _p}𝒵_p(\underset{¯}{𝒓}),\text{with}\underset{¯}{𝒓}=(𝒓_1,\mathrm{},𝒓_p),$$
(3)
where the scaling exponents $`\zeta _p^{dim}=p(2\xi )/2`$ are fixed by dimensionally matching the inertial operator with the forcing, while the $`\zeta _p`$, not constrained by any dimensional requirements, are independent of the forcing. One is usually interested in the $`p`$-point irreducible component of the correlation function, $`C_p`$, that is the structure function $`S_p(r)=(\delta _r\theta )^p`$. Still in the framework of Kraichnan model, the zero modes scaling exponents contributing to $`S_p`$ have been calculated perturbatively in . It has been found that $`\zeta _p<\zeta _p^{dim}`$ for $`p>2`$, i.e. the leading contribution to the structure function scaling is anomalous.

The zero modes dominance scenario explains both the anomalous scaling exponents and their universality. Indeed the forcing does not enter the definition of $`𝒵_p`$, but fixes only the multiplicative constants in front of $`C_p^{}`$ and $`𝒵_p`$, necessary to match the large scale boundary conditions. This implies that, even though scaling exponents are universal, probability density functions (PDF) of scalar increments are not. The above picture based on zero modes is now recognised to apply in all linear hydrodynamic problems. As confirmed by investigations of passive advection by realistic velocity fields , and of shell models for passive transport .

Let us now focus on the main issue of this work: the scaling behaviour in the presence of forcing fluctuations directly injected in the inertial range of scales. We consider the case of a source term $`f(𝒓,t)`$ which is a random Gaussian scalar field, with zero mean, white-in-time and characterised by the spectrum

$$E_f(k)=\pi k|f(𝒌)|^2k^{1+\beta }.$$
(4)
In the following, it is always assumed the presence of ultraviolet and infrared cutoffs for the forcing, i.e. expression (4) holds only in the range $`k_1<|𝒌|<k_2`$ of wave-numbers, where $`k_1`$ and $`k_2`$ are of the order of the inverse of the largest scale in the system and of the dissipative scale, respectively.

Remaining in the framework of the Kraichnan model, Eq. (2) still holds, and its stationary solution is again the superposition of the homogeneous solution, $`𝒵_p`$, which remains unchanged (viz. with the same scaling exponents) and of the inhomogeneous one, $`C_p^{}`$, whose scaling exponents now depend on the slope of the forcing spectrum, $`\beta `$. Indeed by dimensional reasoning, we have $`\zeta _p^{dim}(\beta )=p(2\xi \beta )/2`$.

Two regimes can be identified. If $`\beta <0`$, the scaling of the inhomogeneous solution is always subleading with respect to the anomalous one. The scaling exponents measured from the structure functions are the same as those obtained with the standard large scale forcing. On the other hand, for $`\beta >0`$, there exists a critical order $`p_c`$ such that for $`p<p_c`$: $`\zeta _p^{dim}(\beta )\zeta _p`$, and for $`p>p_c`$: $`\zeta _p\zeta _p^{dim}(\beta )`$. In other words, the scaling behaviour of the low-order structure functions is non-anomalous and dominated by the forcing. The appearing of anomalous scaling for $`p>p_c`$ can be understood due the fact that $`\zeta _p`$ as a function of $`p`$ is concave. Moreover, it is known numerically  and to some degree analytically  that in the Kraichnan model the scaling exponents above a certain order saturate to a finite value, $`\zeta _{\mathrm{}}`$, which makes even more intuitive the existence of a finite $`p_c`$. Some subtle points may arise for positive $`\beta `$ larger than $`1`$. In such a case, the strong UV components of the forcing spectrum may allow a matching with zero modes, disregarded for a large scale forcing , exploding at small scales. Finally, the case $`\beta =0`$ is marginal, because the exponents $`\zeta _p^{dim}`$ coincide with the large scale prediction, up to possible logarithmic corrections.

Checking the validity of the above predictions in numerical simulations of a realistic flow is interesting for two reasons. First, it is a further demonstration, in the Eulerian framework, of the zero modes picture for anomalous scaling and universality in linear problems, which was previously assayed with Lagrangian studies . Second, it offers a controlled testing ground for interpreting some aspects of the nonlinear hydrodynamics which will be considered in the next section.

As an example of realistic velocity field, we consider two-dimensional incompressible Navier-Stokes equations in the inverse cascade regime :

$$_t𝒗+𝒗\mathbf{}𝒗=\mathbf{}P+\eta \mathrm{\Delta }𝒗\alpha 𝒗+𝒇_v.$$
(5)
The terms in the r.h.s. of Eq. (5) have the following meaning: $`P`$ is the pressure field; $`𝒇_v`$ injects kinetic energy at scale $`L_0`$; $`\nu \mathrm{\Delta }𝒗`$ (where $`\nu `$ is the viscosity) dissipates enstrophy at small scale; $`\alpha 𝒗`$ removes energy at the large scales allowing for a statistically stationary state. In the inverse energy cascade regime, the velocity statistics is self-similar (not-intermittent), with $`(\delta _r𝒗\widehat{𝒓})^pr^{p/3}`$, but temporal correlations are non trivial. Moreover, precise numerical measurements of passive scalars with large scale forcing  have shown that: $`\zeta _2=2/3`$ while for $`p>2`$ the exponents are anomalous, and saturation is observed for $`p10`$ with $`\zeta _{\mathrm{}}1.4`$. It is noteworthy that saturation seems to be generic in passive scalars, as found also in experiments . From a physical point of view, it means that the strongest scalar fluctuations are statistically dominated by front-like structures, i.e. huge variations of the field $`\theta `$ at very small scales.

Once a power law forcing (4) is considered, the predicted dimensional scaling is $`\zeta _p^{dim}(\beta )=p(2/3\beta )/2`$. Therefore, the two above described regimes –statistics dominated by zero modes or by the forcing– should appear for $`\beta <0`$ and $`\beta >0`$, respectively.

We performed three sets of direct numerical simulations (DNS) of Eqs. (1),(5): for run (a), we used a standard large scale forcing; for runs (b) and (c), we considered a forcing as in Eq. (4) with $`\beta =0.3`$ and $`\beta =0.3`$, respectively. More details on the DNS are given in the caption of Fig. 1, where we show the snapshots of the scalar field $`\theta `$ for all runs. Already at a first sight, it is possible to observe that runs (a) and (b), where the zero modes dominate over the inhomogeneous solutions at any order, display qualitatively similar features: the field is organised into large scale structures or plateaux characterised by a good mixing, separated by sharp fronts. Differently in run (c), which corresponds to the case of forcing dominated statistics (with $`p_c6`$), the scalar fluctuations develop a wider range of scales, and the most evident large scales plateaux disappear.

In order to make the above observations more quantitative, we studied the scaling of the structure functions $`S_{2p}(r)=(\delta _r\theta )^{2p}`$ (odd orders are zero due to the isotropy of the forcing). For run (a), we found a very good scaling range of about one decade, which allowed us for accurately measuring $`\zeta _p`$ up to $`p1012`$. The quality of the local slopes and the measured values agree with previous investigations , which were performed with larger statistics and higher resolution. For the power-law cases, runs (b) and (c), we could not find evidence of a good scaling as in run (a). This is observed already from the second order structure function $`S_2`$, whose local slope is plotted in Fig. 2 for the three runs.

We understand the poor scaling behaviour as due to the competition between the scaling of the anomalous part and the forcing dominated one. This becomes clear by looking at the scalar flux, $`\mathrm{\Pi }_\theta (r)[\delta _r𝒗\widehat{𝒓}](\delta _r\theta )^2`$. For this quantity under rather general hypothesis, an exact analytical prediction can be derived for a generic forcing, through the Yaglom equation . For a power law forcing $`f`$ as in (4), we obtain:

$$\mathrm{\Pi }_\theta (r)[\delta _r𝒗\widehat{𝒓}](\delta _r\theta )^2c_1r+c_2r^{1\beta },$$
(6)
where we have kept only the first two leading contributions. It is important to notice that constants $`c_1`$ and $`c_2`$, which depend on the space dimension and on the details of the forcing spectrum, turn out to have opposite signs. In (6), in addition to the standard linear term, there is the first leading term induced by the forcing (4). As one can see, the critical value $`\beta =0`$ naturally arises in the flux expression to separate the inertial dominated ($`\beta <0`$) from the forcing dominated ($`\beta >0`$) regime. Indeed, for $`\beta <0`$ the spectral flux is dominated by the low wavenumber components of the forcing spectrum and saturates, in Fourier space, to a constant value as a function of $`k`$. Scalar fluctuations are transferred down-scale via  an intermittent cascade and the statistics is dominated by the anomalous scaling, the forcing contribution even if present at all scales is subleading. For $`\beta >0`$, the scalar spectral flux no longer saturates to a constant: the direct input of energy from the forcing mechanism affects inertial range statistics in a self-similar way, down to the smallest scales where dissipative terms start to be important.

In addition to $`\mathrm{\Pi }_\theta (r)`$, we also studied the third moment of the flux variable $`\mathrm{\Pi }_\theta ^3(r)=[(\delta _r𝒗\widehat{𝒓})(\delta _r\theta )^2]^3`$ whose scaling is anomalous for the large-scale forcing case, i.e. $`\mathrm{\Pi }_\theta ^3(r)r^{2.4}`$ (see also ). Results are summarised in Fig. 3.

Concerning the flux $`\mathrm{\Pi }_\theta (r)`$, for the large scale case run (a), the scaling is good and the Yaglom relation, $`\mathrm{\Pi }_\theta (r)2ϵ_\theta r`$ (being $`ϵ_\theta `$ the scalar dissipation), holds for about one decade (Fig. 3, left panel). On the other hand, for the power-law forced cases, the scaling is poorer, making the identification of the exponents less trivial. Indeed, a scaling behaviour can be properly identified only by taking into account both the leading and subleading terms in Eq. (6). The third order moment is affected by the same problem. So that, on the basis of relation obtained for $`\mathrm{\Pi }_\theta `$, we tried to fit it with the following scaling ansatz:

$$\mathrm{\Pi }_\theta ^3(r)Ar^{2.4}+Br^{3(1\beta )},$$
(7)
where the exponent $`2.4`$ in the first term is associated to the anomalous contribution as measured from run (a). Notice that for run (b) the exponent induced by the forcing is subleading, i.e. $`3(1\beta )=3.9`$, while for run (c) it becomes leading $`3(1\beta )=2.1`$. From the insets of Fig. 3, it is clear that a well defined scaling behaviour is recovered only after having taken into account the two power-law contributions. In all cases the anomalous scaling exponents have values in agreement with those obtained for a large scale forcing, although the interplay of different power laws, when the the forcing is as (4), makes the identification of the exponents very difficult even at considerably high resolution.

Due to these difficulties in measuring the exponents, we looked directly at the PDF of scalar increments $`P(\delta _r\theta )`$. In Fig. 4, we show plots of normalised PDFs of $`\delta _r\theta `$ for different choices of the scale $`r`$ in the inertial range. For runs (a) and (b), the PDFs do not display any rescaling property at changing the scale $`r`$, neither in the core nor in the tails, suggesting intermittent behaviour in the full statistics. Note also that the PDFs in the two runs are different (only scaling exponents are universal while PDFs are not). Conversely, for run (c), where the forcing is dominant, the PDFs display a fairly good rescaling in the core (which governs the low order statistics), while the tails still do not collapse. This agrees with the existence of a critical order, $`p_c`$, above which anomalous scaling appears even in the forcing-dominated case.

As discussed earlier, all differences associated with the two regimes should disappear in the high order statistics, i.e. for $`p>p_c`$ when anomalous scaling exponents imposed by the zero modes should show up irrespectively of the forcing. In testing this point, the presence of saturation for large $`p`$ comes into help, indeed it entails that, for large excursions, the PDFs must approach the form :

$$P(\delta _r\theta )=r^\zeta _{\mathrm{}}Q\left(\frac{\delta _r\theta }{\theta _{rms}}\right)\frac{1}{\theta _{rms}},$$
(8)
with $`\theta _{rms}=\sqrt{\theta ^2}`$ and $`\zeta _{\mathrm{}}1.4`$. Such a prediction is tested and confirmed by the collapse in the PDF tails shown in Fig. 5, which tells us that $`\zeta _{\mathrm{}}`$ is asymptotically approached independently of the forcing.

The above results provide support to the validity of the zero modes picture beyond the boundaries of the Kraichnan model and in agreement with the Lagrangian investigations .

## 3 Nonlinear systems

Observations from the linear problem strongly indicate that the notions of anomalous scaling and universality are closely linked. When we enter the nonlinear world, we do not have any longer a reference theory. Here, we shall rather formulate some working hypothesis suggested by the results from the linear systems and verify their consistency with two cases studied, namely: $`3d`$ Navier-Stokes turbulence, reviewing the work first presented in , and shell models  for turbulence.

### 3.1 The $`3d`$ Navier-Stokes problem

In the $`3d`$ case, we consider random forcing defined by a two-point correlation function that in Fourier space reads

$$f_i(𝒌,t)f_j(𝒌^{},t^{})k^{1y}P_{ij}(𝒌)\delta (𝒌+𝒌^{})\delta (tt^{}),$$
(9)
where $`P_{ij}(𝒌)`$ is the projector assuring incompressibility and the forcing spectrum behaves as $`E_f(k)k^{3y}`$. Hereafter, we use the definition $`y=4\beta `$, referring to the classical notation of the problem as introduced in . Correspondingly, the critical value $`\beta =0`$, separating the two regimes in the linear case, is now $`y_c=4`$.

The relative importance of the stirring on the small scales can be varied by tuning the slope value $`y`$ from $`y0`$ (meaning strong input at all scales) to $`y\mathrm{}`$ (corresponding to a quasi large scale forcing). Also in NS turbulence, a physical insight can be taken from the behaviour of the energy flux, which is constant through scales up to logarithmic corrections for the subleading forcing case, $`y>y_c`$; while it becomes a scale dependent function for $`y<y_c`$, overcoming the energy cascade mechanism.

The case of strong input at all scales, $`y0`$, was originally investigated in Ref.  by means of a RG approach, leading to the following expression for the energy spectrum $`E(k)k^{12y/3}`$, in the domain $`\eta k^1L_f`$ where $`\eta `$ and $`L_f`$ are the viscous scale and large scale of the system, respectively. It is worth noticing that such a prediction, which results also from dimensional analysis , leads to the Kolmogorov spectrum $`E(k)k^{5/3}`$ for $`y=4`$, i.e. quite far from the perturbative region where the RG calculations are under control. Here we want to study the same problems addressed for the linear case, i.e. whether fluctuations are sensitive to the injection mechanism for any $`y`$ and if one observes anomalous scaling for $`y<4`$, when the forcing is dominant. Since, at least with the present knowledge, RG perturbative methods starting at $`y0`$ cannot control anomalous corrections, only numerical simulations at finite $`y`$ values, can possibly give some answers.

In Ref.  a first numerical investigation of $`3d`$ incompressible Navier-Stokes problem was performed. The authors made a set of DNS, varying the spectrum slope from $`y=3`$ to $`y=8`$ at low Taylor’s Reynolds number $`Re_\lambda =22`$. Without considering the issue of anomalous scaling in the region $`y<y_c=4`$, they mostly concentrated on the cases with $`y4`$, and ended up with the conclusion that for $`y4`$ the properties of the statistics are not universal, but varies with the forcing spectrum slope. It is however difficult to consider these as conclusive results, because of the large error bars affecting the data. Later, in Ref. , another numerical study of the same problem at a much higher resolution (up to $`Re_\lambda =220`$) has given some support, even if mostly based on semi-quantitative results, to the scenario drawn in the previous section in the context of passive transport.

In particular, we report here the behaviour of longitudinal velocity increments PDFs for two among the different runs done in , with $`y=3.5`$ and $`y=6`$. For $`y=6>y_c`$, as shown in Fig. 6 left panel, the usual intermittent behaviour for the PDFs is found: the curves $`P(\delta _rv)`$ at different separations $`r`$ do not rescale one onto the other. On the contrary, for $`y=3.5<y_c`$, the PDFs at various scales are almost indistinguishable (see Fig. 6, right panel): a signature of the absence of intermittent corrections at least for the core of the distribution. It is worth noticing that this result is far from trivial, because it coincides with the RG predictions  in the region where RG calculations are well beyond their range of validity. These findings, even if qualitative, fit rather well with those of the passive scalar case (see Fig. 4).

We may try to push forward this indication in the light of the theory for linear systems. For the $`3d`$ Navier-Stokes dynamics, the stationary equations for multi-point correlators $`𝒞_p(𝒓,t)v_{\alpha _1}(𝒓_1)v_{\alpha _2}(𝒓_2)\mathrm{}v_{\alpha _p}(𝒓_p)`$ can be sketched as,

$$\mathrm{\Gamma }_{p+1}𝒞_{p+1}+\nu D_p𝒞_p+_2𝒞_{p2}=0,$$
(10)
where $`\mathrm{\Gamma }_{p+1}`$ is the integro-differential linear operator coming from the inertial and pressure terms, $`D_p`$ is the differential operator describing dissipative effects and $`_2`$ is the two-point forcing correlator. Unfortunately, since the hierarchy (10) is unclosed, no straight-forward analytical approache can be done. Even though nothing can be rigorously said, one may still be tempted to assume that anomalous scaling is brought by the inertial operator. However, the presence of finite size effects, due to the limited inertial range, in DNS data of $`3d`$ Navier-Stokes simulations  provide only a semi-quantitative support to this scenario. As we shall see in the next section, more quantitative statements can be made in the context of shell models for turbulence, for which resolution constraints are much less severe.

### 3.2 Shell models for turbulence

Shell models of turbulence are dynamical systems mimicking Navier-Stokes nonlinear evolution . Their main advantage relies on the possibility of performing high Reynolds number simulations to measure statistical properties, say intermittent corrections, with high accuracy. The model we use, proposed in , is an improved version of the GOY model  (see also  for a recent review). The evolution equation is :

$$\frac{d}{dt}u_n=i(k_nu_{n2}^{}u_{n1}+bk_{n1}u_{n1}^{}u_{n+1}+(1+b)k_{n2}u_{n+1}u_{n+2})\nu k_n^2u_n+f_n.$$
(11)
The velocity field fluctuation at the wavenumber $`k_n`$, with $`k_n=2^nk_0`$, is expressed in terms of the complex variable $`u_n`$; $`b`$ is a free parameter and $`\nu `$ indicates the viscosity. The number of shells varies as $`n=\{0,\mathrm{},N_{max}\}`$. Remarkably, for large scale forcing, the structure functions show anomalous scaling:

$$S_p(k_n)|u_n|^pk_n^{\zeta _p}$$
(12)
deviating from the dimensional prediction $`\zeta _p=p/3`$. Moreover, anomalous scaling exponents are found to be universal with respect to the forcing, provided it is large scale (see ).

We aim at investigating the properties of the model stirred by a white-in-time, Gaussian field, with zero mean and spectrum

$$|f_n|^2=f_0k_n^\beta ,$$
(13)
being $`f_0`$ the forcing intensity. A similar injection mechanism was proposed in  where authors tried to compare the power-law forced shell models with fractal-grid induced turbulence . In order to understand the effect of forcing as (13) let us start with the energy flux, which obeys an exact equation. Denoting with $`\dot{E}_N=\frac{d}{dt}_{n=0}^N|u_n|^2/2`$ the time derivative of the energy content up to the scale $`N`$, the energy balance relation can be written as :

$$\dot{E}_N+\nu \underset{n=0}{\overset{N}{}}k_n^2|u_n^2|=k_N\mathrm{\Pi }_N+\underset{n=0}{\overset{N}{}}|f_n|^2.$$
(14)
In the previous expression, $`\mathrm{\Pi }_N=𝒯_N`$, where $`𝒯_N`$ is the energy flux through the shell $`N`$ defined as

$$𝒯_N=\mathrm{}\left\{u_N^{}u_{N+1}^{}u_{N+2}+\frac{(1+b)}{2}u_{N1}^{}u_N^{}u_{N+1}\right\}.$$
(15)
For the forcing (13) and for $`\nu 0`$, we get an exact prediction in the stationary state:

$$\mathrm{\Pi }_n=A_1k_n^1+B_1k_n^{1+\beta },$$
(16)
where the constants $`A_1,B_1`$, bound by the equation of motion to have opposite signs, depend on the forcing details ($`f_0`$, $`\beta `$) and on the integral wavenumber $`k_0`$. At this level, where everything is exact, we see that when $`\beta `$ approaches the critical value ($`\beta =0`$), a large number of shells is needed to distinguish leading terms from subleading ones in the third order moment. In our runs, the choice of a large value for $`N_{max}(=40)`$, allow us to have a good control of such effects.

Similarly to the NS case (10), we can write the unclosed hierarchy for the multi-point correlators of shell variables:

$$\frac{d}{dt}C_p(\underset{¯}{n})=M_{p+1}(\underset{¯}{n},\underset{¯}{n}^{})C_{p+1}(\underset{¯}{n}^{})+F_2(\underset{¯}{n},\underset{¯}{n}^{})C_{p2}(\underset{¯}{n}^{}),$$
(17)
where $`C_p(\underset{¯}{n})=u_{n_1}u_{n_2}\mathrm{}u_{n_p}`$ is the generic correlation function of order $`p`$, $`M_{p+1}(\underset{¯}{n},\underset{¯}{n}^{})`$ is the linear operator coming from the inertial terms, and $`F_2(\underset{¯}{n},\underset{¯}{n}^{})`$ is the two-point forcing operator. The main difference with the case of passive scalar transport, as pointed out in the previous subsection, is that now the system (17) is not closed, i.e. we have more unknowns than equations. Without any ambition of a rigorous approach, we may imagine that the general solution is characterised by two terms. The first, $`C_p^{}`$, given by the dimensional matching between the inertial operator and the forcing operator,

$$M_{p+1}(\underset{¯}{n},\underset{¯}{n}^{})C_{p+1}^{}(\underset{¯}{n}^{})=F_2(\underset{¯}{n},\underset{¯}{n}^{})C_{p2}^{}(\underset{¯}{n}^{}).$$
(18)
The second, $`𝒵_p`$, associated to a “zero mode” of the inertial operator, namely a solution of the homogeneous equation, $`M_{p+1}(\underset{¯}{n},\underset{¯}{n}^{})𝒵_{p+1}(\underset{¯}{n}^{})=0`$. Here, as for the linear problem, the non homogeneous terms are the only contributions depending on the forcing, while the homogeneous ones, dimensionally unconstrained, may show anomalous scaling. In particular, for the structure functions $`S_p(k_n)`$, we have

$$S_p(k_n)k_n^{\zeta _p}+k_n^{\zeta _p^{dim}(\beta )},$$
(19)
where $`\zeta _p^{dim}(\beta )=p/3(1\beta )`$. This implies the following behaviour for the second order moment of the flux (equivalent to the sixth-order structure function):

$$\mathrm{\Pi }_n^2=𝒯_n^2=A_2k_n^{1.8}+B_2k_n^{2/3(\beta 1)},$$
(20)
where the power $`1.8`$ coincides with the anomalous exponent measured for sixth-order structure function with a large scale forcing, and the value $`2/3(\beta 1)`$ comes from the dimensional prediction (18).

In Fig. 7 the first and second moments of the flux variables are shown for both the subleading and dominant forcing with $`\beta =0.15`$ and $`\beta =0.3`$, respectively. It is easy to see that to get a good agreement of the numerical results with the prediction (20), both the leading and the subleading power laws have to be taken into account. For both cases with $`\beta =0.3`$ and $`\beta =0.15`$, at least for the order of the statistics we could reach, the anomalous scaling exponents have values in agreement with those found for a large scale forcing, supporting the universality scenario.

## 4 Conclusions

We have discussed the problem of small scale fluctuations in turbulent systems stirred at all scales by a power-law forcing. The main question in our study concerns anomalous scaling and universality of scaling exponents.

For linear systems – i.e. passive transport–, we find clear evidence for the universality of anomalous fluctuations and our results fit well in the zero mode scenario. In the regime in which forcing is subleading anomalous scaling is recovered in quantitative agreement with the case of large scale forcing. In the forcing dominated regime, the dimensional scaling imposed by the injection mechanism overwhelms the anomalous fluctuations of low order moments. Nevertheless, anomalous scaling shows up again for high enough moments, independently of the forcing spectrum slope, as confirmed by the existence of a unique saturation exponent in all regimes.

The results obtained for nonlinear turbulent systems point in the same direction. Indeed both the semi-quantitative results obtained in the context of $`3d`$ Navier-Stokes turbulence and those more quantitative issuing from the investigation of the shell models are compatible with the linear systems scenario for universality. However, in the presence of power-law forcing acting in the turbulent energy cascade range, as shown both in the linear and nonlinear cases, scaling properties are strongly spoiled by the beating of leading and subleading terms. This effect is particularly strong due to cancellations induced by different signed pre-factors in the power-law terms.

Strong cancellation effects, which apparently are a mere technical question, may lead to misinterpretation in analysing data. We wonder, for example, if the observed multi-fractal behaviour in the one-dimensional Burgers equation stirred by a power-law forcing might be a spurious effect .

We are grateful to J. Bec, A. Celani and U. Frisch for useful discussions. We acknowledge J. Davoudi who participated to the early stage of the study for passive scalar problem. This research was partially supported in part by the EU under the grant No. HPRN-CT 2000-00162 “Non Ideal Turbulence” and by the INFM (Iniziativa di Calcolo Parallelo). A. L. acknowledges A. Pouquet for useful discussions when visiting NCAR (Boulder, Colorado), and CNR for the support received under the grant “Short-term mobility of researchers 2003”.

## References


format: text

text:
# Cascade control and defense in complex networks

\[Phys. Rev. Lett. 93, 098701 (2004)\]

## Abstract

Complex networks with heterogeneous distribution of loads may undergo a global cascade of overload failures when highly loaded nodes or edges are removed due to attacks or failures. Since a small attack or failure has the potential to trigger a global cascade, a fundamental question regards the possible strategies of defense to prevent the cascade from propagating through the entire network. Here we introduce and investigate a costless strategy of defense based on a selective further removal of nodes and edges, right after the initial attack or failure. This intentional removal of network elements is shown to drastically reduce the size of the cascade.

The problem of attacks on complex networks rv  has attracted a great deal of attention in recent years AJB:2000 ; CEAH:2000 ; CNSW:2000 ; CEAH:2001 ; ML:2002 . The key factor prompting this research is the observation that in many growing networks some nodes evolve to become much more important than others BA:1999 ; ASBS:2000 . From a global perspective, the important nodes are those whose removal may either cause the network to fragment or severely limit the communication between the other nodes. A prime paradigm in this context is that of the scale-free networks (SFNs), which are characterized by an algebraic distribution of degrees (number of edges per node) BA:1999 . While robust against random removals AJB:2000 ; CEAH:2000 ; PV:2001 , SFNs with realistic scaling exponents rv  are likely to fragment into small clusters if a critical fraction of most connected nodes is removed AJB:2000 ; CNSW:2000 ; CEAH:2001 . This happens because of the significant number of edges removed along with the highly connected nodes CEAH:2001 .

But in networks where the flow of a physical quantity is important, such as power grids and computer networks, the fragmentation of the network may be triggered by the removal of one or very few nodes. Indeed, although the removal of few nodes has little effect on the connectivity of the network, these removals may trigger a cascade of subsequent failures which can in turn switch off or disconnect most of the other nodes of the network, such as in the August 14, 2003 event in the northeastern U.S. power transmission network and in recent blackouts in western Europe. A number of important aspects of cascading failures in complex networks have been discussed in the literature, including disturbances in power transmission systems SCL:2000 , the origin of rare events Watts:2002 , the effect of network growth HK:2002 , cascades triggered by intentional attacks ML:2002 , avalanche size distributions MGP:2002 , and congestion instabilities MPVV:2003  (see also Refs. comm ; AAN:2004 ). In particular, a simple model of cascades of overload failures has been introduced ML:2002 . The removal of an even small fraction of highly loaded nodes (due to either attack or failure) has been shown to trigger global cascades of overload failures in SFNs and other networks with heterogeneous distribution of loads. This calls for an investigation on possible strategies of defense to prevent the cascade from propagating through the entire network.

In this Letter we study a general method to reduce the size of cascades of overload failures triggered by attacks on or failures of a small fraction of nodes (hereafter referred to as the initial attack). One such cascade can be divided in two parts: (I) the initial attack, whereby a fraction of nodes is removed; and (II) the propagation of the cascade, where another fraction of nodes is removed due to the subsequent overload failures. The size of the cascade is measured in terms of the ratio $`G=N^{}/N`$, where $`N`$ and $`N^{}`$ are the initial (before (I)) and final (after (II)) number of nodes in the largest connected component, respectively. The method of defense introduced here consists of costless modifications to the network structure after (I) but before (II). In real networks, (I) and (II) are separated in time but this time interval is usually much shorter than the time scale in which the network evolves. Therefore, it is reasonable to consider that no edge or node can be rewired or added to the system after the initial attack because any of these operations would involve extra costs. Accordingly, we assume that, after the initial attack, the only operations allowed in order to reduce the size of the cascade are the intentional removals of nodes and edges. The expression “intentional removals” (IRs) means removals performed after (I) and before (II). The trick point here is that these removals in general reduce even more the final number of nodes in the largest connected component. We show, however, that the IR of carefully chosen nodes and edges can in fact constitute an efficient strategy of defense. Our main result is that the size of the cascade can be drastically reduced with the IRs of nodes having small load and edges having large excess of load (to be defined below). Even though any removal always increases the immediate damage on the network, the resulting $`G`$ is in this case significantly larger (as compared to the case without defense) because these IRs strongly suppress the propagation of the cascade.

For concreteness, we consider the model of overload failures introduced in Ref. ML:2002 , which is defined as follows. For a given network $`𝒲`$, we assume that at each time step one unit of a physical quantity (hereafter called packet) is sent from node $`i`$ to node $`j`$, for every ordered pair of nodes $`(i,j)`$ belonging to the same connected component of $`𝒲`$. We assume that the packet is transmitted along the shortest paths connecting nodes $`i`$ and $`j`$. If there is more than one shortest path connecting two given nodes, the packet is divided evenly at each branching point. The load $`L_k`$ on a node $`k`$ is the total amount of packets passing through that node per unit of time GKK:2001 . Let $`S_k`$ denote the connected component of node $`k`$ and $`L_k^{(i,j)}`$ denote the contribution of the ordered pair $`(i,j)`$ to the load on $`k`$. The load on node $`k`$ is then GKK:2003

$$L_k=\underset{i,j}{}L_k^{(i,j)},$$
(1)
where the sum is over all pairs of nodes in $`S_k`$. Each node $`k`$ is assigned to have a finite capacity $`C_k`$. The node operates in a free-flow regime if $`L_kC_k`$, otherwise the node is assumed to fail and is removed from the network. Now consider that $`𝒲=𝒲(0)`$ is an initially connected network. The initial load $`L_k=L_k(0)`$ is given by Eq. (1) with $`S_k=𝒲(0)`$. The capacity $`C_k`$ of node $`k`$ is assumed ML:2002  to be proportional to the initial load $`L_k(0)`$:

$$C_k=\lambda L_k(0),k=1,2,\mathrm{}N,$$
(2)
where $`\lambda 1`$ is the tolerance parameter and $`N`$ is the initial number of nodes in the network. The cascade can be regarded as a step-by-step process where nodes can be removed at each time step. We start with $`𝒲=𝒲(0)`$ at time $`0`$. The condition $`\lambda 1`$ guarantees that no node of $`𝒲(0)`$ is overloaded, i.e., $`L_k(0)C_k`$ $`k`$. We assume that an initial attack is performed at time 1, whereby a fraction $`p`$ of nodes is removed from $`𝒲(0)`$ and the resulting network is denoted by $`𝒲(1)`$. These removals lead to a global redistribution of loads among the remaining nodes in the network. The updated load $`L_k(1)`$ on a particular node $`k`$ of $`𝒲(1)`$ may then become larger than the capacity $`C_k`$. All the overloaded nodes are removed simultaneously from $`𝒲(1)`$ and the resulting network is denoted by $`𝒲(2)`$. This leads to a new redistribution of loads and subsequent overloads may occur. The overloaded nodes are removed and the resulting network is denoted by $`𝒲(3)`$, and so on. Let $`𝒲(n)`$ denote the updated network at time $`n`$. This cascading process stops only when, for a certain $`n=n^{}`$, the updated load satisfies $`L_k(n^{})C_k`$ for all the nodes $`k`$ of $`𝒲(n^{})`$. The ratio $`G`$ introduced above as a measure of the size of the cascade is defined in terms of the number $`N^{}`$ of nodes in the largest connected component of $`𝒲(n^{})`$.

Here we focus on global cascades \[i.e., $`1G=O(1)`$\] triggered by initial attacks on a small fraction $`p`$ of most loaded nodes. To be specific, we consider a random model of SFNs NSW:2001 , where the degree $`\kappa _i\kappa _0`$ of each node $`i`$ is chosen at random according to the probability distribution $`P(\kappa )\kappa ^\gamma `$, for a given scaling exponent $`\gamma `$ and constant integer $`\kappa _0`$. The degree $`\kappa _i`$ can be regarded as the number of “half edges” emerging from node $`i`$. A network $`𝒲(0)`$ is then generated by randomly connecting half edges to form edges, prohibiting self- and repeated edges sf . The initial attack and the propagation of the cascade correspond to (I) $`𝒲(0)𝒲(1)`$ and (II) $`𝒲(1)𝒲(2)\mathrm{}𝒲(n^{})`$, respectively. Our method of defense consists of an intermediate step $`𝒲(1)𝒲(1)`$ in between (I) and (II), whereby IRs of nodes and/or edges are performed on $`𝒲(1)`$. The network $`𝒲(1)`$ is redefined to incorporate these removals.

We first consider the IR of nodes. In what follows, all the quantities refer to the initial network $`𝒲(0)`$ (unless explicitly mentioned otherwise) and ‘$`(0)`$’ is omitted for simplicity. We recall that $`L_k^{(i,j)}`$ denotes the contribution of the ordered pair of nodes $`(i,j)`$ to the load on node $`k`$. The total contribution of the (unordered) pair $`i`$ and $`j`$ to the load on the network is $`_kL_k^{(i,j)}+_kL_k^{(j,i)}=2(D_{ij}+1)`$, where $`D_{ij}`$ is the shortest path length between nodes $`i`$ and $`j`$. The factor $`2`$ comes from the fact that one packet is sent from node $`i`$ to node $`j`$ and another is sent from node $`j`$ to node $`i`$. We consider that half of this amount is generated by node $`i`$ while the other half is generated by node $`j`$. Then the total load generated by node $`i`$ is

$$L_i^g=\underset{j}{}(D_{ij}+1)=(\overline{D}_i+1)(N1),$$
(3)
where $`\overline{D}_i`$ is the average shortest path length from node $`i`$ to all the others. The average load on each node is $`\overline{L}=_iL_i/N=_iL_i^g/N=(\overline{D}+1)(N1)`$, where $`\overline{D}=_i\overline{D}_i/N`$ is the average shortest path length between any two nodes. Nodes whose load $`L_i`$ is much larger than $`L_i^g`$ contribute much more to handling than to generating load. These are the most important nodes for the network to operate. The removal of one such node may cause overloads on a number of other nodes. More precisely, if a node $`i`$ is removed from the network, the total load on the remaining nodes increases by at least the amount $`L_i2L_i^g`$, unless the removal of node $`i`$ divides the network into more than one connected component. Along with the observation that nodes having large $`L_i`$ tend to have small $`L_i^g`$ corr , this provides a rationale for attack strategies based on the removal of highly loaded nodes. On the other hand, nodes whose load $`L_i`$ is smaller than $`L_i^g`$ generate more load than they handle. This observation is our starting point to argue that the size of the cascade can be drastically reduced with the IR of a certain fraction $`f`$ of nodes according to any of the following strategies: (1) nodes with smallest $`\mathrm{\Delta }_iL_iL_i^g`$ are removed first; (2) nodes with smallest closeness centrality $`\overline{D}_i^1`$ are removed first; (3) nodes with smallest load $`L_i`$ are removed first; (4) nodes with smallest degree $`\kappa _i`$ are removed first. The IRs are performed on $`𝒲(1)`$, but for simplicity we assume that $`\mathrm{\Delta }_i`$, $`\overline{D}_i`$, $`L_i`$, and $`\kappa _i`$ correspond to the initial values computed on $`𝒲(0)`$. In random SFNs, all these quantities are strongly correlated corr : the quantities $`L_i^g`$ and $`L_i`$ are negatively correlated, while the quantities $`\mathrm{\Delta }_i`$, $`\overline{D}_i^1`$, $`L_i`$, and $`\kappa _i`$ are positively correlated. Therefore we only need to justify one of the strategies of defense (1)-(4). Consider then strategy (1).

Strategy (1) consists in removing a fraction $`f`$ of nodes with most negative $`\mathrm{\Delta }_i`$. In order to reduce the size of a cascade triggered by an initial attack, these IRs have to satisfy two conditions. The first condition is to reduce the load on the remaining nodes in the network. Because $`\mathrm{\Delta }_i`$ is negatively correlated with $`\overline{D}_i`$, the removal of nodes with $`\mathrm{\Delta }_i<0`$ tends to reduce the average shortest path length $`\overline{D}`$ between the nodes that remain connected to the largest component. The average load $`\overline{L}`$ on a node of this component is proportional to both $`\overline{D}+1`$ and the number of other nodes in the same component, and is therefore reduced with strategy (1). This argument is based on the initial network $`𝒲(0)`$ but, for large random networks, the same is expected to hold true on $`𝒲(1)`$ as well. The second condition is that the fragmentation caused by the IRs must be smaller than that otherwise caused by the cascade itself, i.e. $`\stackrel{~}{G}>G^o`$, where $`\stackrel{~}{G}`$ and $`G^o`$ denote the fraction of nodes remaining in the largest connected component right after the IRs and after the cascade without defense, respectively. Because $`\mathrm{\Delta }_i`$ is positively correlated with the degree $`\kappa _i`$, the nodes removed according to strategy (1) tend to be the least connected nodes in the network. It is well known that most of the (unremoved) nodes remain in a single connected component when any fraction of least connected nodes is removed AJB:2000 ; CEAH:2000 . More specifically, $`\stackrel{~}{G}`$ is expected to decrease linearly as $`\stackrel{~}{G}1f`$ with the fraction $`f`$ of nodes removed according to strategy (1). (We assume that the fraction $`p`$ of nodes removed in the initial attack is small enough so that we can neglect the contribution of the initial attack to $`\stackrel{~}{G}`$.) Therefore, the second condition is satisfied even for relatively large $`f`$ insofar as $`f<1G^o`$. The ratio $`G`$ for the cascade with defense is $`G^{}=\stackrel{~}{G}\mathrm{\Delta }\stackrel{~}{G}`$, where $`\mathrm{\Delta }\stackrel{~}{G}`$ is due to the propagation of the cascade. Gathering all these, because $`\stackrel{~}{G}`$ decreases slowly with increasing $`f`$, there should be a certain $`f<1G^o`$ for which $`\stackrel{~}{G}`$ is large (as compared to $`G^o`$) and the average load on the remaining nodes is sufficiently reduced so that the propagation of the cascade is strongly suppressed and $`\mathrm{\Delta }\stackrel{~}{G}`$ is small. Therefore strategy (1) is expected to be effective in reducing the size of global cascades. Because of the correlations mentioned above, the same is expected for the strategies (2)-(4).

On the other hand, the network is expected to be sensitive to the removal of nodes with large $`\mathrm{\Delta }_i`$ and, due to the correlations corr , to the removal of nodes with large $`\overline{D}_i^1`$, $`L_i`$, and $`\kappa _i`$ as well. All these quantities can be regarded as measures of centrality. Therefore our first result could be stated as follows: while the removal of the most central nodes of $`𝒲(0)`$ can trigger global cascades, the removal of the least central nodes of $`𝒲(1)`$ can drastically reduce the size of these cascades.

Now we present numerical verification of our result concerning the IR of nodes. We consider random SFNs with scaling exponent $`\gamma `$ and initial attacks on a small fraction $`p1`$ of most loaded nodes. Strong evidence for our result is presented in Fig. 1(a), where we show the ratio $`G`$ as a function of the tolerance parameter $`\lambda `$ for $`\gamma =3.0`$ and $`p=0.001`$. Without defense, this initial attack on only $`0.1\%`$ of the nodes triggers global cascades even for relatively large values of the tolerance parameter $`\lambda `$ \[Fig. 1(a), stars\]. However, the ratio $`G`$ is shown to be significantly larger when a suitable fraction of nodes is intentionally removed according to any of the strategies of defense (1)-(4) \[Fig. 1(a), open symbols\]. For example, for $`\lambda =1.5`$, we have $`G0.6`$ with defense and $`G0.06`$ without it. A similar improvement is observed for other values of the scaling exponent $`\gamma `$. As a function of the fraction $`f`$ of IRs, the ratio $`G`$ displays a well-defined maximum, as shown in Fig. 1(b) for $`\lambda =1.5`$. When $`f`$ is large, the propagation of the cascade is strongly suppressed and nearly all the damage is caused by the IRs, i.e., $`G`$ is approximately $`1f`$. This explains the linear behavior of $`G`$ for $`f>0.4`$ in Fig. 1(b). When $`f`$ is small, most of the damage is caused by the cascade itself. The maximum of $`G`$ lies in a region of intermediate $`f`$ where the propagation of the cascade is significantly suppressed and the damage caused by the IRs is relatively small. The results presented in Fig. 1(a) correspond to this maximum. The almost perfect agreement between the different strategies of defense in Figs. 1(a) and 1(b) is due to the strong correlations between loads, path lengths, and degrees in random SFNs.

We now turn to the IR of edges. We argue that the size of the cascade can be drastically reduced with the IR of edges not necessarily connected to the nodes removed in the IR of nodes. In analogy to the load on nodes, we define the load $`L_{ij}^e`$ on an edge between nodes $`i`$ and $`j`$ as the total amount of packets passing through that edge per unit of time. The load $`L_i`$ on node $`i`$ can be expressed in terms of the load on edges as

$$L_i=\frac{1}{2}\underset{ji}{}L_{ij}^e+(N_i1),$$
(4)
where the sum is over all the edges directly connected to node $`i`$ and $`N_i`$ is the number of nodes in the corresponding connected component $`S_i`$. Now consider the network $`𝒲(1)`$ right after the initial attack. From Eq. (4), for the load $`L_i(1)`$ on a node $`i`$ to be larger than the capacity $`C_i=\lambda L_i(0)`$, the load on at least one of the edges of node $`i`$ must exceed the initial load by a factor larger than $`\lambda 1`$, i.e., $`L_{ij}^e(1)>\lambda L_{ij}^e(0)`$ for some edge connected to $`i`$. We say that this edge has large excess of load. The IR of one such edge is expected to reduce the load on node $`i`$ (except for exceptional cases where the load on $`i`$ remains unchanged). From this follows that the IR of edges having large excess of load can be used to reduce the size of the cascade. However, for these removals to be effective we need to take into account that the removal of an edge may cause overloads elsewhere. As a simple strategy that takes this into account, we propose to remove edges from $`𝒲(1)`$ according to an auxiliary cascade of overload failures on edges. This is not to be taken as a real cascade but instead as a criterion to remove edges. Specifically, we consider the following: (i) each edge has an imaginary capacity $`C_{ij}^e=\xi L_{ij}^e(0)`$, where $`\xi 1`$ is a tunable parameter; and (ii) the resulting cascade on edges consists of a step-by-step process where all the edges satisfying $`L_{ij}^e(1)>C_{ij}^e`$ are removed simultaneously from the network and $`𝒲(1)`$ is redefined at each step to incorporate these removals. Note that the cascade on edges, and hence the IR of edges, takes place after the initial attack (I) and before the propagation of the cascade on nodes (II). Because the condition $`L_{ij}^e(1)>\lambda L_{ij}^e(0)`$ is necessary but not sufficient for node $`i`$ to be overloaded, we expect this strategy to be effective for suitable $`\xi >\lambda `$. Without making effort to optimize the result, in our computations we set $`\xi =2\lambda `$. The size of the cascades is significantly reduced with this strategy, as shown in Fig. 1(a) (solid circles) for random SFNs with scaling exponent $`\gamma =3.0`$. For example, for $`\lambda =1.5`$, the IR of edges increases $`G`$ by a factor of $`6`$. The effectiveness of this strategy for intermediate values of $`\xi `$ expresses a trade-off between the suppression of the cascade on nodes (enhanced for smaller $`\xi `$) and the fragmentation of the network due to the IR of edges (reduced for larger $`\xi `$).

In summary, we have introduced a general method to reduce the size of cascades of overload failures triggered by attacks on or failures of highly loaded nodes. The method is based on the IR of nodes and edges before the propagation of the cascade. We have shown that the size of the cascade can be drastically reduced with the IR of nodes having small load and/or edges having large excess of load. Here we have focused on SFNs and we emphasize that similar results are expected for other networks with heterogeneous distribution of loads. In the model considered above the packets are transmitted along the shortest paths between nodes, but a heterogeneous distribution of load is also expected for many models with local routing. Even in a model where the packets are routed randomly, packets tend to visit highly connected nodes much more often than other nodes NR:2003 . If, on one hand, this heterogeneity in the load distribution makes the system vulnerable to cascading failures, on the other hand, it makes the removal of poorly loaded nodes and excessively loaded edges an effective method to reduce the size of the cascades.

The author thanks Eduardo Guéron for illuminating discussions.


text:
# Neutron Matter : A Superfluid Gas

## 1 INTRODUCTION

Neutron matter is believed to be an unusual superfluid gas with positive pressure at all densities. It has been studied in the context of neutron stars for many years. Recent reviews can be found in references . However, the subject is rapidly progressing; the large number of authors in this paper represents the breadth of recent interest.

The interesting properties of neutron matter follow from the fact that the $`nn`$ scattering length ($`a=18`$ fm) has a large magnitude, so that the dimensionless parameter $`|ak_F|1`$ at densities $`>10^4\rho _0`$, where $`k_F`$ is the Fermi momentum given by $`k_F^3=3\pi ^2\rho `$. We use $`\rho _0=0.16`$ fm-3 to denote the density of nuclear matter in large nuclei. On the other hand, the effective range of the $`nn`$ interaction ($`r_{nn}2.8`$ fm) is not much smaller than the interparticle distance $`r_0=(3\rho /4\pi )^{1/3}`$, in neutron matter at densities of interest. Hence the $`nn`$ interaction can not be easily approximated by a zero range $`\delta `$-function interaction. Realistic models of the $`nn`$ interaction have a repulsive core and strong tensor component; thus simple perturbation theory expansions can not be easily used to describe neutron matter. In addition, due to the strong $`nn`$ attraction in $`{}_{}{}^{1}S_{0}^{}`$ channel, the low-density neutron matter is believed to be an $`S`$wave superfluid, although at higher densities it could be a $`P`$wave superfluid .

Recently it has become possible to study cold dilute Fermi gases in atom traps . In these gases the interatomic interaction can be approximated by a zero range potential characterized by the scattering length $`a`$. The experiments have the ability to tune values of $`ak_F`$ by using Feshbach resonances. These studies essentially probe superfluid gases like neutron matter at very low densities ($`r_02.8`$ fm).

Density-dependent effective interactions are commonly used in nuclear physics to develop mean-field approximations based on energy-density functionals. In absence of experimental data, results of many-body calculations of the equation of state of pure neutron matter  are often used to determine the effective interactions .

In section 2 we briefly describe the progress in quantum Monte Carlo calculations of the properties of dilute Fermi gases and neutron matter. Section 3 is devoted to quantum Monte Carlo calculations of neutrons in a potential well aimed to approximate neutron-rich oxygen isotopes. We discuss effective interactions in neutron matter in section 4, and conclude in section 5.

## 2 UNIFORM SUPERFLUID GASES AND NEUTRON MATTER

Recently it has become possible to calculate the energies and gaps of simple superfluid Fermi gases ab initio from the interparticle interaction, using quantum Monte Carlo methods . The total energy is calculated for $`N`$ particles in a periodic box of volume $`N/\rho `$. The results for dilute (effective range $`r_0`$) Fermi gas with $`ak_F=\mathrm{}`$ are shown in Fig.1. The unit of energy:

$$E_{FG}=\frac{3}{10}\frac{\mathrm{}^2}{m}k_F^2,$$
(1)
corresponds to the energy per particle of noninteracting Fermi gas. Dimensional arguments show that the energy and gap of dilute superfluid gases scale with $`E_{FG}`$ when $`a\mathrm{}`$. The pairing gap $`\mathrm{\Delta }`$ is obtained from the difference between results for odd and even $`N`$ (Fig.1).

These studies of dilute superfluid Fermi gases have been extended to various values of $`ak_F`$ in the domain $`\mathrm{}<(ak_F)^1<\mathrm{}`$. Cubic periodic boxes have shell structure with shell closures at 2, 14, 38, … for spin-$`\frac{1}{2}`$ fermions. In dilute gases the pairing removes shell effects when the interaction is attractive and strong enough so that $`(ak_F)^1>`$0.5. The gap energy defined as:

$$\mathrm{\Delta }(N)=E(N)\frac{1}{2}[E(N1)+E(N+1)],oddN,$$
(2)
does not depend significantly on $`N`$ for $`(ak_F)^1`$1. At negative values of $`ak_F`$, the calculated values of $`\mathrm{\Delta }`$ are in between the estimates obtained from BCS and Gorkov equations . However, at large positive values of $`(ak_F)^1`$, the gap energy is below the Gorkov estimate and approaches $`|E_2|/2`$, where $`E_2`$ is the energy of the two-body ground state. Recall that when the interaction is attractive and $`a>0`$, there must be one or more two-body bound states.

As $`(ak_F)^1`$ varies from $`\mathrm{}`$ to $`+\mathrm{}`$ the character of dilute Fermi gases changes from BCS superfluids with weak pairing to diatomic Bose gases with Bose-Einstein condensation. This transition has been discussed in the literature .

The pairing in nuclei and neutron matter is much weaker than in the dilute Fermi gases with $`(ak_F)^1>`$0.3. However, the energy of neutron matter appears to be close to the energy $`(0.44E_{FG})`$ predicted by quantum Monte Carlo calculations (Fig.1) for large $`|ak_F|`$, over a wide density range. Most calculations of neutron matter with realistic interactions, carried out since 1970 , give $`E(\rho )0.5E_{FG}(\rho )`$ in the density range $`0.01`$ to $`0.5\rho _0`$. The pairing gaps in neutron matter at $`\rho >0.1\rho _0`$ are sensitive to the details of the $`nn`$ interaction. For example, the gaps predicted from only the scattering length and effective range are much larger than those with realistic models with repulsive cores .

Quantum Monte Carlo calculations of neutron matter with realistic interactions are much more difficult than those for dilute gases. In order to treat the effects of the spin, tensor and spin-orbit components in the $`nn`$ interaction exactly, one has to sum over all the $`2^N`$ spin states of the $`N`$ interacting neutrons . Such Green’s function Monte Carlo (GFMC) calculations are limited to systems with no more than 16 neutrons at present. Methods to sample the spin space using auxiliary fields are also being developed. These auxiliary field diffusion Monte Carlo (AFDMC) calculations  are fairly accurate for pure neutron systems, and can be carried out for much larger values of $`N`$.

In computations using periodic boxes the interparticle interaction beyond $`L/2`$, where $`L`$ is the box length, has to be truncated or approximated. GFMC calculations  use truncated interactions and estimate the contribution of $`v(r>L/2)`$ by approximate methods. It is small for $`\rho \rho _0/2`$, but increases rapidly with density. The studies of pairing gaps in neutron matter have therefore focused on $`\rho \rho _0/4`$. The present, preliminary, results are listed in Table I.

The effect of shell closure at $`N=14`$ is obvious in the results shown in Table I. The energy of the seventh pair, $`E(14)E(12)3.2`$ MeV, is much smaller than that of the eighth pair, $`E(16)E(14)26.8`$ MeV. Calculations with smaller statistical (Monte Carlo) error are necessary to study the dependence of $`\mathrm{\Delta }(N)`$ on $`N`$, in addition to studies of $`\mathrm{\Delta }(N>15)`$ with the AFDMC method. Nevertheless the $`\mathrm{\Delta }(\rho _0/4)`$ appears to be of order $`2.8(5)`$ MeV. This value is in excellent agreement with that obtained using realistic bare $`nn`$ interactions and the BCS equation, 2.6(1) MeV. This agreement is surprising since most estimates  suggested that medium polarization corrections to the bare interaction would suppress the gap by a factor of two or more even at a density of $`\rho _0/4`$.

The BCS equation predicts a maximum gap of 3 MeV at $`0.1\rho _0`$, and the disappearance of the $`{}_{}{}^{1}S_{0}^{}`$ pairing at $`0.7\rho _0`$ . It would be most interesting to verify these predictions with essentially exact ab initio calculations.

## 3 Neutron Drops : Models of Neutron-Rich Nuclei

Neutron drops bound in external wells were studied previously with various objectives . We are presently exploring the possibility to model neutron-rich oxygen isotopes by neutron drops in a well meant to represent the effect of interactions with the protons. The model Hamiltonian for the isotope with $`N`$ neutrons has the form:

$$H=\underset{i=1,N}{}\left(\frac{\mathrm{}^2}{2m}_i^2+V(r_i)\right)+\underset{i<jN}{}v_{ij}+\underset{i<j<kN}{}V_{ijk},$$
(3)
where $`v_{ij}`$ and $`V_{ijk}`$ are realistic bare two-, and three-neutron interactions. The Argonne $`v_{18}`$ and Illinois-2 interactions are used in the present calculations and the single-neutron potential $`V(r)`$ is chosen so that the density distribution of the eight-neutron ground state is close to that of protons in 16O. It is given by:

$$V(r)=\frac{35.5}{1+\mathrm{exp}[(\mathrm{r}3)/1.1]},$$
(4)
in MeV and fm units. This well has only two bound states; a 0s and a 0p at $``$ -12.1 and -4.8 MeV respectively. It can thus bind only 8 noninteracting neutrons. However, bound states of many more interacting neutrons exist in this well due to the $`nn`$ attraction.

The calculated energies of the low-energy bound states of up to 14 neutrons are compared with the experimental energies of oxygen isotopes in Fig.2. A constant of -19 MeV has been added to the calculated energies to match the energy of the 8-neutron drop with that of 16O. The results indicate that neutron-rich ($`N>8`$) oxygen isotopes can be well approximated by neutron drops in a Woods-Saxon well. However, the approximation becomes less accurate for neutron-poor ($`N<8`$) isotopes. Presumably the proton distribution in oxygen isotopes does not change much with $`N`$ when $`N>8`$, but when $`N`$ becomes smaller than 8, the protons become less bound and their distribution changes significantly.

In the traditional shell model, neutron-rich oxygen isotopes are approximated by $`N8`$ neutrons interacting via an effective interaction in the sd-shell. In contrast the present method uses bare interactions among all the neutrons. The single-particle potential $`V(r)`$ is meant to represent only the $`np`$ interaction effects. It can obviously be improved by adding a spin-orbit term to the $`V(r)`$.

In 17O the interactions between neutrons give only $`2/3`$ of the observed $`d_{3/2}d_{5/2}`$ spin-orbit splitting, and the three-neutron interaction contributes approximately half of the calculated splitting. The remaining third must come from the interactions with protons. The variation of the pairing gap with $`N`$ is being studied in more detail.

## 4 The Effective Interaction in Neutron Matter

Density-dependent effective interactions are often used to describe nuclei using the mean-field approximation. They are mostly used in the Hartree, Hartree-Fock or Hartree-Fock-Bogolubov equations. In principle, parts of such effective interactions which depend only on the density of matter, can be chosen such as to reproduce the $`E(\rho )`$ and the pairing gap $`\mathrm{\Delta }(\rho )`$ of uniform neutron matter.

The simplest Skyrme effective interaction between neutrons is of the form:

$$v_{eff}(𝐑,𝐫)=b(\rho (𝐑))\delta (𝐫).$$
(5)
Here $`𝐑`$ and $`𝐫`$ are the center of mass and relative separation between the two interacting neutrons. Commonly used Skyrme interactions have momentum dependent and spin-orbit terms which are omitted here for brevity. Their dependence on the proton density can also be neglected in the present discussion. In the Hartree-Fock approximation, the energy of neutron matter per neutron given by this $`v_{eff}`$ is:

$$E(\rho )=E_{FG}(\rho )+\frac{1}{4}b(\rho )\rho .$$
(6)
The strength of the Skyrme interaction $`b(\rho )`$ can thus be obtained from the $`E(\rho )`$ of neutron matter.

Most realistic calculations of neutron matter give $`E(\rho )0.5E_{FG}(\rho )`$ in the density range $`\rho =0.01`$ to 0.1 fm-3. This result is independent of the interaction model and the calculation details . It simply follows from the fact that in this density range the $`|a_{nn}k_F|1`$. This implies that:

$$b(\rho )\frac{2}{\rho }E_{FG}(\rho )\alpha \rho ^{1/3},$$
(7)
since $`E_{FG}`$ is proportional to $`\rho ^{2/3}`$. This low-density divergence of the effective $`nn`$ interaction is also a consequence of the large magnitude of the $`nn`$ scattering length.

In some models  the neutron-proton effective interaction has a more severe $`(1/\rho )`$ divergence due to the existence of the $`np`$ bound state, the deuteron. However, this divergence is mostly ignored because matter containing neutrons and protons is unstable at low densities; it clusters into nuclei. In contrast the $`\rho ^{1/3}`$ divergence of the effective $`nn`$ interaction could have consequences on the structure of neutron-rich nuclei with a neutron halo and in the neutron gas in a neutron-star crust .

It is known that pairing is very important in determining the structure of neutron-rich nuclei. For example, nuclei with an even number of neutrons dominate the neutron-drip line. It is not possible to use delta-function effective interactions of finite strength to calculate pairing gaps without restricting the Hilbert space. In unrestricted space divergences occur. In stable nuclei reasonable values of the gaps may be obtained using the Skyrme effective interaction and truncating the Hilbert space to a single open shell . However, in nuclei near the neutron drip line, shell structure becomes less dominant, and the problem becomes more difficult. Alternative approaches of using a local density approximation for the pairing fields  are also being pursued.

## 5 Conclusions

Advances in quantum Monte Carlo Methods in the new millennium have made possible essentially exact, ab initio calculations of the equations of state of low-density neutron matter. We are also developing methods to calculate the pairing gaps in low-density neutron matter, and in dilute Fermi gases from the bare interactions. All realistic models of nuclear forces predict similar properties of low-density neutron matter, thus these ab initio predictions are expected to have little model dependence. Ways to define effective interactions in neutron matter from the results of these ab initio calculations are being refined. They can presumably be used to study the structure of neutron-rich nuclei. In addition, studies of the bound states of neutrons in external potential wells meant to represent the effects of the neutron-proton interaction may also provide new insights in the structure of neutron-rich nuclei.

## 6 ACKNOWLEDGMENTS

The calculations were made possible by grants of time on the parallel computers of Argonne and Los Alamos National Laboratory and the National Energy Research Scientific Computing Center. The work of SYC, JM, VRP, and DGR is supported by the U. S. National Science Foundation via grant PHY 00-98353, that of JC by the U. S. Department of Energy under contract W-7405-ENG-36, and that of SCP and RBW by the U. S. Department of Energy, Nuclear Physics Division, under contract W-31-109-ENG-38.