N Tronko
email: n.tronko@warwick.ac.uk

S V Nazarenko

S Galtier

Weak turbulence in two-dimensional magnetohydrodynamics

Keywords: numbers: 47.27 E, 47.65-d, 52.30 Cv, 52.35 Bj, 95.30 Qd

A weak wave turbulence theory is developed for two-dimensional (2D) magnetohydrodynamics (MHD). We derive and analyze the kinetic equation describing the three-wave interactions of pseudo-Alfvén waves. Our analysis is greatly helped by the fortunate fact that in 2D the wavekinetic equation is integrable. In contrast with the 3D case, in 2D the wave interactions are nonlocal. Another distinct feature is that strong derivatives of spectra tend to appear in the region of small parallel (i.e. along the uniform magnetic field direction) wavenumbers leading to a breakdown of the weak turbulence description in this region. We develop a qualitative theory beyond weak turbulence describing subsequent evolution and formation of a steady state.

I. INTRODUCTION

Magneto-Hydro-Dynamics (MHD) is of great interest for modeling turbulence in magnetically confined and unconfined plasmas. In astrophysics its applications range from solar wind [START_REF] Marsch | [END_REF], to the Sun [START_REF] Priest | Solar Magnetohydrodynamics[END_REF], to the interstellar medium [START_REF] Ng | [END_REF] and beyond [4]. At the same time, MHD is also relevant to large-scale motion in nuclear fusion devices such as tokamaks [5]. One of the pioneering results of incompressible MHD turbulence has been obtained by Iroshnikov [6] and Kraichnan [7] (thereafter IK) who proposed an extension of the Kolmogorov phenomenology [8], originally derived for Navier-Stokes equations for Hydro-Dynamics (HD). For simplicity the assumptions of homogeneity and isotropy were made by Kolmogorov, and the energy cascade was supposed to be dominated by local (in scale) interactions between eddies of similar size. Then the Kolmogorov phenomenology leads to the well-known one-dimensional kinetic energy spectrum E(k) ∼ k -5/3 , where k is the wave vector. The associated cascade properties for its inviscid invariants differ for 3D and 2D turbulence: while in 3D the energy and the kinetic helicity exhibit direct cascades, in 2D the energy cascades inversely -still with a -5/3 scaling -whereas a direct cascade is found for the enstrophy which leads to spectrum E(k) ∼ k -3 at the small scales.

IK modified the Kolmogorov phenomenology by taking into account the magnetic field.

They also assumed turbulence homogeneity, isotropy and locality of interactions. However, there exist fundamental differences between the Kolmogorov and the IK theories. First of all, in MHD the energy cascade is supposed to be dominated not by the interactions between eddies but between Alfvén wave packets propagating in opposite directions; this modification leads to the energy spectrum E(k) ∼ k -3/2 . Furthermore, unlike hydrodynamics the cascades of the ideal MHD invariants exhibit some similarities in their behaviour in 2D and 3D [START_REF] Biskamp | Magnetohydrodynamic turbulence[END_REF]: in both cases, the cascades have the same direction with a direct cascade for the energy and cross-helicity, and an inverse cascade for the magnetic helicity or the anastrophy.

The differences between MHD and HD turbulence go beyond these classical properties.

In the IK theory the large-scale magnetic field is supposed to play the role of external field which is necessary for the existence of Alfvén waves but its main effect, i.e. anisotropy, is not taken into account. The importance of an external magnetic field has been discussed many times during the last two decades [START_REF] Montgomery | [END_REF][11][12][13][14][15][16][17][18] and the anisotropic behavior has been shown in direct numerical simulations for both 2D [11] and 3D [12].

Despite some similarities -like the cascade directions -the question about the identification of differences between 2D and 3D MHD turbulence still represents an important issue.

In early numerical studies [11] mainly 2D simulations were performed because of the limited numerical resources available and the non-accessibility of 3D calculations at high Reynolds numbers. Nowadays, 3D MHD numerical simulations are commonly achieved [19][20][21][22][23][24][25][26]. At the same time, 2D simulations are still used for the illustration of new numerical techniques [27]. There is also an interest for the investigation of freely decaying MHD turbulence because the Reynolds numbers can be higher in 2D than in 3D [23]. When the magnetic field perturbations are small compare to a uniform background magnetic field, the 2D MHD equations are sometimes used to model turbulence. Such a situation is particularly relevant for solar coronal loops as well as for reduced models for fusion plasmas in slab geometry [5].

Strong statements were made by some authors that 2D simulations can be safely used to model 3D situations because the properties of the 2D and the 3D MHD turbulence are essentially the same [23,[START_REF] Biskamp | Nonlinear Magnetohydrodynamics[END_REF]. One of the motivations of the present paper is to test validity of this claim in a special case when the external magnetic field is strong.

In this paper, we consider 2D MHD in the presence of a strong background magnetic field which implies realization of the weak turbulence regime. One of the main advantages of this regime is the fact that it allows one to derive accurate analytical results for the spectrum. An explicit comparison will be made between the weak turbulence regimes in 2D and in 3D; the latter was analyzed rigorously in [START_REF] Galtier | [END_REF]. The weak wave turbulence approach is widely familiar to the plasma physics community [30][START_REF] Akhiezer | Plasma Electrodynamics II: Non-linear Theory and Fluctuations[END_REF][START_REF] Mcivor | [END_REF][33][START_REF] Zakharov | Kolmogorov spectra in weak turbulence problems[END_REF][START_REF] Zakharov | Kolmogorov spectra of turbulence 1. Wave turbulence[END_REF][START_REF] Kuznetsov | [END_REF][37][START_REF] Tsytovich | Nonlinear Effects in Plasma[END_REF][START_REF] Sagdeev | Nonlinear Plasma Theory[END_REF][START_REF] Vedenov | Theory of Turbulent Plasma[END_REF][START_REF] Galtier | [END_REF][42]. It is a statistical description of a large ensemble of weakly interacting dispersive waves. The formalism leads to wave kinetic equations from which exact power law solutions can be found for the energy spectra. There were several reasons which postponed development of weak turbulence for Alfvén waves. The first one is their semi-dispersive nature. Typically, the wave-kinetic approach cannot be used for nondispersive waves since such wavepackets propagate with the same group velocity even if their wavenumbers are different, the energy exchange between such waves may not be considered small, and may lead to possible energy accumulation over a long time of interaction. The Alfvén waves represent a unique exception to this rule because co-propagating waves do not interact and the nonlinear interaction is present only for counter-propagating waves.

The latter pass through each other in some finite time and no long time cumulative effect occur. That is why the Alfvén waves represent a unique example of semi-dispersive waves for which the wave turbulence theory applies. The second reason which renders the weak turbulence theory for Alfvén waves very subtle is the fact that the domination of threewave interactions -as assumed by IK -may be questionable. While in [43] the three-wave interactions were declared absent, the IK argument has been re-established in [13,[START_REF] Galtier | [END_REF]44].

The weak turbulence theory for 3D incompressible MHD was developed in [START_REF] Galtier | [END_REF] (see also [45,[START_REF] Nazarenko | Wave turbulence[END_REF]) where three-wave kinetic equations were derived with their exact solutions via a systematic asymptotic expansion in powers of small nonlinearities.

The main goal of the present paper is to derive the weak turbulence equation for the 2D MHD, analyze it and make a comparison with the 3D case. The crucial technical step which allows a comprehensive theoretical analysis of the solutions consists in transforming the wave kinetic equation into an integrable form by Fourier transforming it and separating the transverse and the parallel dynamics by using a self-similar "effective time" variable. This article is organized as follows. In section II we derive the weak turbulence kinetic equation through a general perturbative procedure. In sections III and IV we proceed with a detailed investigation of its properties in an anisotropic limit. The main goal of such a study is to verify whether or not the turbulence is local. First of all, we will consider the steady state behavior by looking for Kolmogorov-Zakharov type solutions and check their locality.

Next, we will proceed with investigation of an unsteady spectrum evolution by considering two different cases with a gaussian-shaped source and different kinds of dissipation: a uniform friction and a viscosity. Due to integrability of the weak wave-kinetic equation in the first case it is possible to find an exact solution. In the second case, a qualitative analysis for the steady state is complemented by a numerical simulation of the spectrum evolution. The goal of section V is to develop some qualitative reasoning about the turbulent behavior of our system near the applicability margin of wave-kinetic formalism and beyond. Formation of steady state is also discussed. Finally, we present a summary of our results in section VI.

II. WAVE-KINETIC DESCRIPTION A. Alfvén waves

In 3D incompressible MHD there exist two different kinds of Alfvén waves [15]. The first kind, called Shear-Alfvén waves (SAW), have fluctuations of velocity and magnetic field transverse to the initial background magnetic field B 0 , whilst the other kind, called Pseudo-Alfvén waves (PAW), have fluctuations along B 0 . Both waves propagate along B 0 at the same Alfvén speed.

The weak wave turbulence approach for incompressible MHD applies for a small nonlinearity, ǫ ∼ b ⊥ k ⊥ /B 0 k || ≪ 1, were b ⊥ is the perpendicular magnetic field perturbation and k || and k ⊥ are the wave vector components in the parallel and the perpendicular directions to B 0 . Additionally, a strong anisotropy condition is often used, σ = k ⊥ /k || ≫ 1. In the 3D case, it was shown that in the leading order of weak nonlinearity, ǫ ≪ 1, and strong anisotropy, σ ≫ 1, the SAW interact only among themselves and evolve independently from the PAW. At the same time, the PAW scatter from the SAW without amplification or damping, and they do not interact with each other. Such a behaviour does not rule out a possibility for the PAW to be interacting among themselves in the next order of expansion in 1/σ. However in the 3D case, such a process is sub-dominant to a stronger interaction with the SAW and was not considered yet.

In the 2D case, due to the geometrical restrictions, it is only possible to have the PAW and not SAW. In this paper we will see that three-wave interactions of PAW do occur in the 2D case in the next order of expansion in 1/σ and represent the dominant process in the nonlinear evolution.

B. Interaction representation.

The ideal incompressible MHD system in Elsässer variables z = v + sb, with s = ±1, is given by [START_REF] Biskamp | Nonlinear Magnetohydrodynamics[END_REF]

∂ t -sB 0 • ∇ + z -s • ∇ z s = -∇P * , (1)
∇ • z s = 0, ( 2
)
where v is the fluid velocity, b is the magnetic field fluctuation (in velocity units), B 0 is a uniform background magnetic field (also in velocity units, i.e. the Alfvén speed) and P * is the total (thermal plus magnetic) pressure. In what follows we suppose that the background magnetic field is directed along the x axis, B 0 = B 0 x. In the coordinate notations we have

(∂ t -sB 0 ∂ x ) z s j = -ǫ z -s n ∂ n z s j -∂ j P * , (3)
∂ i z s i = 0. ( 4
)
The nonlinear terms in eq.( 3) include Elsässer variables of the opposite signs only. Therefore, the nonlinear interactions take place only between counter-propagating waves.

The first step in the general procedure of the wave kinetic formalism is to identify the linear modes. Neglecting the nonlinear terms in the right hand side (r.h.s.) of the equation

(3) (which includes the pressure term) and looking for solutions in the form of a wave,

z s j ∼ e i(kxx+kyy)-iω s t , (5)
we get two linear modes,

ω s = -sB 0 k x , s = ±1, (6)
which propagate parallel to the background magnetic field (in both directions) with group velocity v s g = -sB 0 . Let us suppose that our system is periodic in the physical space with period L in both x and y, and let us introduce the Fourier series:

z s j (x, t) = k a s j (k, t)e ik•x , (7)
where wave vector k takes values on a 2D grid, k = (k x , k y ) = (2πm x /L, 2πm y /L) where m x , m y ∈ Z. Then, by applying the divergence operation on both sides of the equation ( 3)

and by using (4), we find the expression for the Fourier coefficients of pressure P * :

P * (k) = -k -2 k 1 ,k 2 (k 2 • a -s (k 1 , t))(k • a s (k 2 , t))δ(k 1 + k 2 -k), (8)
where δ(p) is the Kronecker delta: δ(p) = 1 for p = 0 and zero otherwise. Thus, equation

(3) in Fourier space becomes:

(i∂ t -ω s )a s (k, t) = k 1 ,k 2 (k • a -s (k 1 , t)) a s (k 2 , t) - k k 2 (k • a s (k 2 , t)) δ(k 1 + k 2 -k). (9)
Using the incompressibility condition

a s x = -a s y k y k x , (10)
we reduce [START_REF] Biskamp | Magnetohydrodynamic turbulence[END_REF] to one scalar equation,

(i∂ t -ω s )a s x (k, t) = ǫ k 1 ,k 2 k y (k × k 1 ) z (k • k 2 ) k 1y k 2y k 2 a -s x (k 1 , t) a s x (k 2 , t)dk 1 dk 2 . ( 11
)
Let us now introduce the representation of interaction variables,

c s k = i k k y ǫ a s x (k, t)e iω s t , (12)
which represent slowly varying wave amplitudes. Factor e iω s t here compensates for the fastscale oscillations arising due to the linear dynamics. We have introduced a formal constant small parameter ǫ ≪ 1 for easier counting of the powers of nonlinearity, assuming now that c k ∼ 1. Then the system (9) in the interaction representation variables is:

∂ t c ± k = ǫ k 1 ,k 2 V 12k c ∓ k 1 c ± k 2 e 2ik 1x t δ(k 1 + k 2 -k), (13)
with the interaction coefficient:

V 12k = (k • k 2 ) [k 1 × k 2 ] z k k 1 k 2 . ( 14
)
Note that so far we have not used smallness of ǫ and the eq.( 13) is completely equivalent to the initial system (3)-( 4).

C. Wave-kinetic equation

The standard weak turbulence approach [START_REF] Zakharov | Kolmogorov spectra of turbulence 1. Wave turbulence[END_REF][START_REF] Nazarenko | Wave turbulence[END_REF] exploits the smallness of nonlinearity, randomness of phases and the infinite box limit. In Appendix A, we are applying this approach to the system ( 13) and obtain the following kinetic equation for the wave spectrum

n k , ṅ± k = π V 2 k12 n ∓ k 1 n ± k 2 -n ± k δ (k -k 1 -k 2 ) δ(2k 1x )dk 1 dk 2 , (15)
where the interaction coefficient is as in expression (14). Here, δ(p) means Dirac's delta function. In the following sections we proceed with detailed analysis of this equation.

D. Anisotropic limit

One remarkable property of MHD turbulence, which makes it very different from the HD one, is its strong anisotropy in presence of strong background magnetic field. It was illustrated with direct numerical simulations in both 2D [11] and 3D [12]. The wave-kinetic formalism confirms such an anisotropy through the kinetic equation structure. In fact for Alfvén waves, the resonant three-wave interaction [11] is organised in such a way that one member of each triad must have its wave vector perpendicular to the external magnetic field. At the same time, the two other waves in the same triad must have their parallel wavenumbers equal to each other: k || = k 2|| . Formally, this is seen in both the 2D and the 3D kinetic equations whose r.h.s. contains a delta function δ(2k 1|| ), see equation (15) for the 2D case and the equation ( 26) in [START_REF] Galtier | [END_REF] for the 3D case. Using this delta function and integrating over k 1|| we see that the parallel component of the wavenumber enters into the kinetic equation as an external parameter and the spectrum dynamics is decoupled at each level of k || . In other words, there is no energy transfer in the parallel (to the external field B 0 ) direction in the k-space. The initial spectrum is spreading over the transverse wavenumbers k ⊥ , and not in the k || direction. For large time, such a spectrum becomes very flat, pancake-like.

The two-dimensionalisation of the total energy means that, for large time, the energy spectrum is supported on a volume of wave-numbers such that for most of them k ⊥ ≫ k || .

Thus, let us consider the anisotropic limit for kinetic equation (15), which reads in 2D as k y ≫ k x . Taking into the account the resonant interaction conditions for the parallel wavenumbers, we will have a considerable simplification of the interaction coefficient:

V k12 = k x , (16)
and the kinetic equation will take the following form,

ṅ± (k x , k y ) = πk 2 x n ∓ (0, k 1y )[n ± (k x , k 2y ) -n ± (k x , k y )]δ(k y -k 1y -k 2y )dk 1y dk 2y . ( 17
)
This equation describes three-wave interactions of PAW in 2D in the anisotropic limit. One can immediately see that the energy is conserved separately in the "+" and "-" waves separately at each k x :

∂ t n ∓ (k x , k y , t)dk y = 0. ( 18
)
Factor k 2 x in the r.h.s. of relation ( 17) is very important. In the 3D case, there exists a similar term which corresponds to a sub-leading contribution. We remind that in 3D in the leading order of the perturbation theory, the PAW are scattered on SAW and do not interact directly to each other. In the 2D, there is no SAW and, therefore, the r.h.s. of (17) becomes the leading order contribution.

Further, in leading order of 3D there is no k 2

x factor, and substitution n(k ⊥ , k || , t) = n ⊥ (k ⊥ , t)n || (k || ) leads to an equation for n ⊥ (k ⊥ , t) which does not involve k || . In 2D, one can also obtain an equation which does not involve k x but for this, one has to introduce an "effective time" variable τ = πk 2

x t:

∂ τ n ± (k x , k y ; τ ) = n ∓ (0, k 1y ; 0) n ± (k x , k 2y ; τ ) -n ± (k x , k y ; τ ) ×δ (k y -k 1y -k 2y ) dk 1y dk 2y . (19)
Now we can seek solution in the following form,

n(k x , k y , τ ) = µ(k x )η(k y , τ ), (20)
where µ(k x ) represents the parallel (non-evolving) component of the energy spectrum and η(k y , τ ) is the perpendicular one. Without loss of generality we can assume µ(0) = 1.

Substituting ( 20) into (19) we have the following equation for η,

∂ τ η ± (k y , τ ) = η ∓ (k 1y , 0) η ± (k 2y , τ ) -η ± (k y , τ ) δ (k y -k 1y -k 2y ) dk 1y dk 2y . (21)
Thus, like in 3D, we have an evolution equation for the perpendicular part of the spectrum which does not explicitly depend on the parallel part. However, the qualitative difference with the 2D case is that there is an implicit dependence on k x via the effective time variable τ , which in particular leads to the fact that in the r.h.s. of equation ( 21) one of η's is taken at τ = 0, making this equation linear and, as we will soon see, integrable. Another distinct feature arising from such an implicit dependence on k x via τ is sharpening of the spectrum at small k x leading to the breakdown of the wave-kinetic description. Later in this paper we will study this effect and its consequences.

III. KOLMOGOROV-ZAKHAROV SPECTRA AND THEIR LOCALITY

At the first step of our investigation of the wave-kinetic equation ( 21), we seek exact stationary power law solutions, η(k y ; ∞) ± ∝ k ν ± y . For this, we will use so-called Zakharov transformation,

k ′ 1y = k y k 1y k 2y , k ′ 2y = k 2 y k 2y . ( 22
)
We split the integral into the r.h.s. of (21) in two equal parts and we perform the Zakharov transformation in the integrand of one of these parts. This leads to the following conditions on the exponents,

ν + + ν -= -2 . ( 23
)
However, the resulting power law spectra represent true mathematical solutions if and only if the original (before the transformation) integral in the r.h.s. of (21) converges at these solutions. In Appendix B, we perform such a convergence study and show that this integral never converges and therefore, no Kolmogorov-Zakharov solutions are possible for 2D MHD system, contrary to the 3D case. We remark that the balanced turbulence spectrum with ν + = ν -= -1 has a logarithmic divergence. Thus, one could anticipate that such a marginal nonlocality could be "fixed" by logarithmic corrections. We will see later that this is false.

IV. INTEGRATION OF THE KINETIC EQUATION

A remarkable property of the kinetic equation ( 19) is its simplicity. In this section we show that in some physical situations it can be solved analytically. Let us introduce into the equation ( 19) sources and sinks of waves:

∂ τ n ± (k x , k y , τ ) = n ∓ (0, k 1y , 0) n ± (k x , k 2y , τ ) -n ± (k x , k y , τ ) ×δ (k y -k 1y -k 2y ) dk 1y dk 2y + F (k y , k x ) -σ d n(k x , k y , τ ). (24)
The function F (k x , k y ) may represent forcing or dissipation, depending on the choice of the sign before it, and the constant σ d introduces a uniform friction.

In order to use factorisation (20) and eliminate µ(k x ) in both sides of the forced kinetic equation, we assume the following type of force/dissipation function,

F (k y , k x ) = F x (k x )F y (k y ).
Then the parallel component of (20) must be chosen as µ(k x ) = F x (k x ).

Finally, we obtain the following forced/dissipated kinetic equation for the perpendicular component of the energy spectrum,

∂ τ η ± (k y , τ ) = η ∓ (k 1y , 0) η ± (k 2y , τ ) -η ± (k y , τ ) ×δ (k y -k 1y -k 2y ) dk 1y dk 2y + F y (k y ) -σ d η(k y , τ ) . (25)
A. Pseudo-physical space

A considerable simplification of equation ( 25) can be obtained with performing the inverse Fourier transform on η(k y , τ ):

E ± (y, τ ) = η ± (k y , τ ) e ikyy dk y . (26)
We will call E ± (y, τ ) the pseudo-physical space energy, keeping in mind that what is transformed is the spectrum and not the original wave variable.

We arrive at the following representation of equation ( 25) in the pseudo-physical space,

∂ τ E ± (y, τ ) = E ∓ (y, τ ) E ± (y, 0) -E ± (0, 0) -σ d + F (y). (27)

B. General solutions

Let us consider the balanced turbulence case with E + (y, τ ) = E -(y, τ ). Then, the general solution of equation ( 27) can be written as:

E(y, τ ) = C(y)e (E(y,0)-E(0,0)-σ d )τ - F (y) E(y, 0) -E(0, 0) -σ d , (28)
where the first term represents the general solution for the homogeneous equation and the second term is a particular (time independent) solution of the inhomogeneous equation.

Function C(y) has to be fixed by the initial condition,

C(y) = E(y, 0) + F (y) E(y, 0) -E(0, 0) -σ d . (29)
Now let us consider two particular examples of the forcing and dissipation. In both cases we will assume a Gaussian shape forcing, F (y) = σ f e k 2 f y 2 /2 , where constants σ f and k f represent the forcing strength and its characteristic wave vector respectively (in the k y -space the forcing is also Gaussian, centered at k y = 0 and with width k f ). In the first example, the dissipation will be represented by a uniform friction. Here, we can write analytical solutions of the kinetic equation on the pseudo-physical space. In the second case, we will consider a viscous dissipation. Here, a qualitative analysis of the stationary regime can be done. In order to illustrate the spectrum evolution in that case, a numerical solution will be used.

1. Uniform friction.

For the uniform friction case we have:

E (y, τ ) = C(y)e (E(y,0)-E(0,0)-σ d )τ - σ f e k 2 f y 2 /2 E(y, 0) -E(0, 0) -σ d . ( 30
)
For simplicity, let us use a single-wave initial condition

E(y, 0) = 2A cos (k 0 y) , A = const > 0, (31)
which corresponds to two δ-functions in k y -space: at k y = ±k 0 .

Then, we can find a function C(y) using equation ( 29) and substitute it into our solution, which yields:

E(y, τ ) = 2A cos(k 0 y) + σ f e -k 2 f y 2 /2 2A (cos(k 0 y) -1) -σ d e (2A(cos k 0 y-1)-σ d )τ - σ f e -k 2 f y 2 /2 2A (cos(k 0 y) -1) -σ d . ( 32
)
Let us examine the steady state, which corresponds to the limit t → ∞ and, therefore, τ → ∞. Note, however, that the time for the steady state to form becomes longer as k x gets less, and there are always very small k x where the spectrum is evolving at any large time.

In the limit τ → ∞ the solution is given by the second term in the r.h.s. of [START_REF] Mcivor | [END_REF]. Far from the initial and the forcing scales, at k ≫ k 0 and k ≫ k f , which corresponds to y ≪ 1/k 0 and y ≪ 1/k f , we have cos(

k 0 y) = 1 -(k 0 y) 2 /2 + O((k 0 y) 4 ) and e -k 2 f y 2 /2 = 1 + O((k f y) 2 ).
Thus for this range of scales we have the following expression for the steady state solution in the pseudo-Fourier space,

E (y, ∞) = σ f σ d + λy 2 , where λ = Ak 2 0 . (33)
Performing Fourier transform of E(y, ∞) we get the steady spectrum,

η (k y , ∞) = 1 2π ∞ -∞ E (y, ∞) e -ikyy dy. ( 34
)
For wavenumbers in the inertial range k

0 , k f ≪ k ≪ k d = λ/σ d , expression (33)
becomes effectively a delta function in the integrand of [START_REF] Zakharov | Kolmogorov spectra in weak turbulence problems[END_REF],

σ f σ d + λy 2 ≈ πσ f √ σ d λ δ(y), (35)
and we have

η (k y ; ∞) = 1 2 σ f √ σ d λ = 1 2 σ f σ d Ak 2 0 . ( 36
)
Therefore we can reach the conclusion that in the equilibrium state the energy spectrum of our system in the inertial range is flat. Formally, it is a power law with exponent ν = 0 which is very different from the Kolmogorov-Zakharov exponent ν = -1 found in section III. Recall that the Kolmogorov-Zakharov in the balanced case was found to be marginally nonlocal and the common wisdom would suggest that it could be fixed by a log correction.

As we now see this is false: our exact solution has a completely different exponent and has no log factor.

We also see that our exact solution is nonlocal: it does not just depend on the energy flux but it contains information about both the sources and the sinks as well as about the initial conditions.

2. Viscous friction.

Let us now replace the uniform friction by a viscous dissipation keeping the same one-wave initial condition as before. Equation ( 27) becomes:

∂E(y; τ ) ∂τ = 2A (cos(k 0 y) -1) E(y; τ ) + σ ν ∂ 2 E(y; τ ) ∂y 2 + σ f e - k 2 f y 2 2 , ( 37
)
where σ ν denotes now the viscosity coefficient and we have used the initial conditions [START_REF] Akhiezer | Plasma Electrodynamics II: Non-linear Theory and Fluctuations[END_REF].

To realise this estimation, we need first to get the expression for the steady state solution.

Let us examine the steady state concentrating, like before in the uniform friction case, on the scales less than the forcing and the initial scales, which in terms of the pseudo-physical space variables means y ≪ 1/k 0 and y ≪ 1/k f . Performing the same type of expansion in small y as before, we have

σ ν d 2 E(y; ∞) dy 2 -λy 2 E(y; ∞) + σ f = 0. ( 38
)
By performing the following rescaling

y = y λ σ ν 1 4 , E = √ λσ ν σ f E ( 39
)
we obtain

d 2 E( y, ∞) d y 2 -y 2 E( y, ∞) + 1 = 0. ( 40
)
The homogeneous part here is the equation of the parabolic cylinder. Its solutions are the parabolic cylinder special functions, whose properties and asymptotics can be found e.g.

in [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. Qualitatively, the behaviour is similar to the one we had in the previous (friction dissipation) example: it has a maximum at y = 0 and it decays for y → ∞ (faster than in the previous example). In fig.

(1) we present such a solution in the pseudo-physical space obtained using Matlab in the interval with y ∈ [-6.1, 6.1]. In order to obtain the decaying solution we need to take E(0) = 1.3110288959 with high accuracy (to eliminate contribution of the growing parabolic cylinder function).

Respectively, in the k y -space we again have a flat spectrum in the inertial range,

η (k y ; ∞) = Cσ f σ ν λ 3 1 4 for k 0 , k f ≪ k ≪ k ν = (λ/σ ν ) 1 4 , ( 41
)
where C is an order one constant. Once again we see that the spectrum is nonlocal (i.e.

it is dependent on the details of the forcing and the sink parameters rather than just the energy flux), and its exponent is zero (i.e. it is not a log corrected Kolmogorov-Zakharov spectrum).

In order to illustrate the dynamical evolution of the spectrum in the viscous dissipation case, we perform direct numerical simulations for the kinetic equation (25) in the k y -space.

In this equation we take η + = η -, σ d = 0 and F y (k y ) = F f orce -σ ν k 2 y with σ ν = 10 -6 . The result is presented in Fig. [START_REF] Priest | Solar Magnetohydrodynamics[END_REF]. The black curve on this figure represents the initial spectrum with a Gaussian shape with a large-scale forcing F f orce = 3 * 10 -4 /k realised for k ∈ [START_REF] Ng | [END_REF][START_REF] Biskamp | Magnetohydrodynamic turbulence[END_REF].

We see that the system converges toward a steady state with a flat spectrum in the inertial range.

V. BEYOND WEAK TURBULENCE.

In this section we are providing, at a qualitative level of rigor, a description for the energy spectrum behaviour beyond the weak turbulence regime at the late stage of the evolution.

For the wave-kinetic equation to be valid the nonlinearity has to remain weak, i.e. the nonlinear time scale should be much longer than the linear wave period. This results in the following condition,

t nl t L = b y k y B 0 k x < 1. (42)
Here we have used the estimate for the nonlinear evolution time taking the standard hydrodynamic non-linear time scale. In our dynamical equations, the nonlinearity related to the magnetic field is of the same order. Finally, the applicability condition can be rewritten as a limitation for the parallel wave number:

k x > k * x = b y k y B 0 . (43)
It means that for small values of parallel wave vector component near k * x the kinetic equation becomes invalid. The question about applicability of the wave-kinetic equation near k || = 0 has been frequently discussed in the literature, eg. in the 3D MHD turbulence context [START_REF] Galtier | [END_REF].

In particular it was speculated in [START_REF] Galtier | [END_REF] that a sufficient spectrum smoothness for wavenumbers near k || = 0 must be present for the wave-kinetic equation to be applicable. In the 2D case, the smoothness of spectrum near k x = 0 is asymptotically (in time) broken due to a special structure of the kinetic equation.

Indeed, as we mentioned before, the wave-kinetic equation for the 2D MHD system is formulated in terms of a self-similar "time" variable τ = k 2

x t. Therefore dependence on the parallel wavenumber is still present in the perpendicular part of the energy spectrum η(k y , τ ) in an implicit way, via τ . Such a self-similar dependence on k x is manifested, at each fixed k y , in shrinking of the original k x profile along the k x -axis as time grows (see Fig. 3). The spectrum is narrowing and its derivative is growing near small values of k x , and when it is so steep that a significant variation occurs over the range ∼ k *

x , the kinetic equation breaks down. Time estimate for such a breakdown is t ∼ (k *

x ) -2 .

Therefore, the weak turbulence description will break down at the late evolution stages, and the wave-kinetic equation will no longer work. However, it is possible to amend this description to take into account the strongly nonlinear effects and develop a qualitative theory of subsequent evolution leading to a steady state. Below we will present a qualitative argument which will allow us to obtain such a theory.

First of all we note that the three-wave interaction is never exactly resonant: it involves all the quasi-resonant within a certain small distance from the exact resonant frequencyso called nonlinear resonance broadening Γ ∼ t -1 nl . In the other words, the delta function in kinetic equation δ(2k 1x ) = δ (ω k -ω 1 -ω 2 ) should be replaced by a peaked function f (k 1x ) with a small but finite width Γ. For sufficiently smooth spectra the difference from the delta-function can be ignored, but for sharp and narrow spectra the integrand in the kinetic equation's integral become itself peaked and the delta-function broadening becomes important. Its main effect is that of a filter f (k x ) in the k x variable which acts to smoothen any sharp changes over the range ∆k ∼ k *

x . The energy is no longer conserved separately at each fixed k x .

For large wavenumbers k x ≫ k *

x (where the spectrum remains slowly varying even when it is steep at k x ∼ k *

x ) the kinetic equation could be easily amended by replacing

n ∓ (0, k y , 0) = n ∓ (k 1x , k 1y , τ ) δ(k 1x )dk 1x with n ∓ (k y , τ ) = n ∓ (k 1x , k 1y , τ ) f (k 1x )dk 1x .
For small wavenumbers, k x ∼ k *

x , the effect of the resonance broadening is not reduced to such a simple modification of just one function in the integral. It is clear that the spectrum at k x = 0, which was fixed in the wave-kinetic approximation, will suffer changes caused by the smoothing in the direction determined by spectral slope at small k x : if the gradient is positive (negative) the value will increase (decrease), as illustrated in Fig. 4. The details of the evolution at the small wavenumbers are not important because the combined action of the self-similar shrinking and smoothing will lead to a very rapid wipeout of all the gradients in k x and formation of a steady state with η independent of k x . Correspondingly, the values of η at k x = 0 will adjust themselves to the values at k x = ∞. After this moment, when the rapid dependencies on k x disappear, the kinetic equation in its usual weak turbulence form is valid once again, and it can be used for finding the final steady state spectrum. Since η is now independent of k x , the steady state could be readily obtained from the formal condition η(k y , 0) = η(k y , ∞) which simply means our solution independent of τ ; it has nothing to do with the initial/final values of the spectrum in time or k x . Thus, the evolution can be summarized as follows. At the early stage, t ≪ (k * x ) -2 , the evolution is described by the three-wave kinetic equation. Then at the advanced stage, with characteristic times scales t ∼ (k * x ) -2 , the kinetic equation is broken down by its own evolution. Smoothing of strong gradients in k x occurs, which results in spectrum stabilisation and re-emergence of the kinetic equation description at the large time scales t ≫ (k *

x ) -2 . This kinetic equation describes spectrum evolution within the steady state regime. Let us now consider the properties of such a steady state.

A. Spectrum evolution in the steady-state

Let us now analyse the steady state. Based on what was said in the end of the previous section, we will seek a τ -independent solution of the pseudo-Fourier space equation ( 27): (formally coinciding with the condition E(y, 0) = E(y, ∞) = E(y)). Considering this equation at y = 0 we have E(0) = F (0)/σ d . Also we have E(0) = η(k y )dk y > 0 (because η(k y ) ≥ 0). Solving the quadratic equation, we have:

E 2 (y) -E(y) (E(0) + σ d ) + F(y) = 0 (44)
E(y) = 1 2 (σ d + F (0)/σ d ) ± 1 2 (σ d + F (0)/σ d ) 2 -4 F (y) 1/2 . ( 45
)
To satisfy condition E(0) = F(0)/σ d , we must choose "+" if F (0) > σ 2 d and "-" otherwise. We suppose that the forcing decays at infinity, lim y→∞ F (y) → 0, which is the case eg.

for the Gaussian forcing. This means that we do not force the k y = 0 mode. Then we see

that lim y→∞ E(y) → 0 if F(0) < σ 2 d and lim y→∞ E(y) → σ d + F (0)/σ d if F (0) > σ 2 d .
In the second case we have a spectrum with a delta function at k y = 0. Thus we observe an interesting phenomenon of condensation into k y = 0 mode in the cases when the forcing prevails over the dissipation at the small scales y (corresponding to high k y 's). In the first case E(y) is a monotonously decreasing (to 0) function of y, where as in the second case it is monotonously increasing (asymptoting to constant).

The first case is physically more relevant, because in most cases of interest dissipation dominates over forcing at the small scales. In this case E(y) behaves qualitatively similar as in the two examples considered in section IV. Namely, if we take the same Gaussian forcing as in these two examples, we will have E(y) which has a maximum at y = 0, smooth everywhere (including y = 0) and rapidly decaying to zero for y → ∞. However, there is a big difference from the previous examples in that now the characteristic width of function E(y), and respectively the width of the spectrum in the k y variable, is of the same order as the width of the forcing function. Therefore, there is no inertial range in the final steady state considered here. This is an even stronger case of the nonlocal interaction than in the two examples considered before. Both the forcing and the dissipation parameters enter in the final answer, but not the parameters of the initial condition: the steady state beyond the weak turbulence has already forgotten all the initial data.

VI. SUMMARY.

We have shown that the three-wave interactions for the pseudo-Alfvén waves (PAW) in the 2D MHD system are non-empty and it is possible to obtain a three-wave kinetic equation within the weak turbulence approach. These interactions take place in the second order of the anisotropy parameter.

We found Kolmogorov-Zakharov power law spectra for PAW in 2D MHD system and showed that they are not realisable due to divergence of the collision integrals of the kinetic equation. In the balanced case this divergence is marginal. This is an indirect indication that the 2D PAW turbulence is nonlocal: it is dominated by interaction of waves with very different in size wavelengths. Our full analytical solution of the kinetic equation confirms such a nonlocality. It also dispels the myth that all marginally nonlocal spectra can be "fixed" by a logarithmic correction.

The crucial technique for our analysis is passing to pseudo-physical space via Fouriertransforming the kinetic equation and using a self-similar effective time variable. This has allowed us to dramatically simplify the kinetic equation, solve it analytically in some important cases and to fully analyse in the other important cases. The two main examples we analysed have a Gaussian-shaped forcing of low wavenumbers and a dissipation represented by either uniform friction or by viscosity. The first case is solvable analytically, and the second one is shown to possess a similar behaviour. Namely, the spectrum evolves independently at each k x and it tends to a flat steady-state spectrum in the inertial range (which is not a log-corrected Kolmogorov-Zakharov spectrum).

At each fixed k y , the spectrum develops sharp gradients in k x at small k x , which eventually leads to the breakdown of the weak turbulence description. We present a qualitative argument about what follows after this moment. We argue that the effect of strong turbulence is to smoothen the sharp gradients via the nonlinear resonance broadening effect. This will lead to a steady state with no gradient in k x for which the weak turbulence kinetic equation formally works once again, and we present an analytical solution for such a steady state.

On the practical side, one should derive from our work a warning that the 2D and the 3D MHD systems are dramatically different, and one should be careful when extrapolating the 2D results, eg. numerical ones, onto the 3D case. Indeed, in contrast to 2D, in 3D there is no gradient sharpening at small parallel wave numbers, and the Kolmogorov-Zakharov spectrum is a local and well behaved solution.

VII. ACKNOWLEDGEMENT

Sergey Nazarenko gratefully acknowledges support from the government of Russian Federation via grant No. 12.740.11.1430 for supporting research of teams working under supervision of invited scientists.

• when α -< 0 the expression for the main contribution is made by: At the second singular point, after integration out k 2y using the δ-function in (B1), we have: The second integral is always convergent, and the first one is convergent for α -> -1.

Finally, the convergence region for the first collisional integral (B1) in the space of indices is:

{{(α + + α -< -1) ∩ (α -> 0)} ∪ {(α + < -1) ∩ (α -< -0)}} ∩ (α + > -2) ∩ (α -> -1) .

It is represented by the grey trapezoid in Fig. 5.

To find the convergence zone for the second integral of (B1) (with α -in the exponent of k 1y ) we just take reflection of the convergence zone for the first integral with respect to the line α -= α + . Finally, to get the convergence conditions for both collisional integrals one should take the intersection of the both zones.

As we can see in Fig. 5 such an intersection produces a zero set. There are no power law exponents α -and α + for which both collision integrals would be convergent, and there is a single point which corresponds to marginal (logarithmic) divergence, α -= α + = -1. This point corresponds to the Kolmogorov-Zakharov spectrum in the balanced turbulence case.

FIG. 1 .

 1 FIG. 1. Stationary solution in the pseudo physical space.

FIG. 2 .

 2 FIG.2. Spectrum dynamics for a viscous dissipation case τ = 10000.

FIG. 3 .

 3 FIG. 3. Spectrum narrowing for large time scales.

FIG. 4 .

 4 FIG. 4. Gradients smoothing process. Four iterations for spectrum value stabilisation are presented on this figure. At the first stage, the gradient of spectrum in the vicinity of the k x = 0 (see point 1) is positive, the initial value of spectrum η(k y , 0) increases and reaches point 2, and then, after crossing the maximum, it moves to point 3, which corresponds to negative slope of the spectrum. Then, the initial value decreases and arrives at the position 4. This process will continue until the spectrum stabilises at η(k y , 0) = η(k y , ∞).

∞ 1 |k

 1 1y | α + |k y | α -dk 1y ,and is convergent for α + < -1.

ǫ 0 |k

 0 1y | α + (|k y -k 1y | α --|k y | α -) dk 1y ∼ the convergence condition α + > -2.To get this condition, we have performed the series expansion:|k y -k 1y | α -= |k y | α -(1 + α -k 1y /k y ) + .. . , and we have used the fact that the integral ǫ 0 x ν dx is convergent for ν > -1.To obtain the convergence condition for the last singular point, we integrate out k 1y using the δ-function:ǫ 0 |k y -k 2y | α + (|k 2y | α --|k y | α -) dk 2y ∼ (B3) ǫ 0 |k y | α + |k 2y | α -dk 2y -ǫ 0 |k y | α + +α -dk 2y .

Appendix A: Derivation of the wave-kinetic equation

In section II A we wrote the 2D MHD system in the interaction representation (13), which comprises a starting point for derivation of the wave-kinetic equation. Let us define the wave spectrum as

where the average is taken over the random initial conditions, and L 2 is the area of the periodic box. With this normalisation n k tends to a finite limit ∼ ǫ 2 as L → ∞ provided that the wave density is finite and uniform in the 2D physical space.

The next step consists in making use of the time scales separation. We are introducing an intermediate time scale T which should be much smaller than the typical non-linear time, t nl = 2π/(ǫ 2 ω), and much greater than the linear wave period, t L = 2π/ω. Taking T = 2π/(ǫω) will satisfy these conditions, t L ≪ T ≪ t nl . Then we are looking for solutions at time t = T in the form series in small ǫ:

where we suppose that the lowest order amplitudes c ±,(0) k = c ± k (0) correspond to the linear regime.

For the spectrum we have:

After substituting expansion (A1) into equation ( 13) in the first order we have:

where

For the second order we can write:

Next, we are going to assume that the initial amplitudes c (±,0) k are Gaussian random variables which are statistically independent at each k, and use Wick's rule:

We also remember that because the physical space amplitudes are real functions we have

We have:

where we have used abbreviations

and

Let substitute expressions (A8) and (A9) into the eq. (A2),

where we have used that |∆ T (2k 1x )| 2 = sin 2 (k 1x T )/k 2 1x . Now we take the infinite box limit, L → ∞, and pass to the continuous description in the k-space using the rule 1

where δ k 12 in the integrand means Dirac's delta (recall that it is Kronecker delta in the sum). At the next stage of the wave-kinetic procedure we need to use the weakness of the nonlinearity in our system by taking the limit ǫ → 0, which is equivalent to T → ∞. For the r.h.s. of the eq. (A12), we obtain:

Then after multiplying both parts of the eq.(A12) by 1/T , its l.h.s. becomes:

where we took into account that time T is much less than the nonlinear time at which the spectrum evolves.

After these steps we can finally write down the kinetic equation:

Appendix B: Locality study for Kolmogorov-Zakharov solutions.

In order to explore realisability of Kolmogorov-Zakharov spectra, η(k y ; ∞) ± ∝ k ν ± y , we need to proceed with a convergence study of the collisional integrals: Then two cases are possible:

• when α -> 0, the main contribution is made by: ∞ 1 |k 1y | α + +α -dk 1y , which is convergent for α + + α -< -1,


N Choungmo Fofack

Philippe Nain
email: philippe.nain@inria.fr

Giovanni Neglia
email: giovanni.neglia@inria.fr

Don Towsley
email: towsley@cs.umass.edu

Analysis of TTL-based Cache Networks

Many researchers have been working on the performance analysis of caching in Information-Centric Networks (ICNs) under various replacement policies like Least Recently Used (LRU), FIFO or Random (RND). However, no exact results are provided, and many approximate models do not scale even for the simple network of two caches connected in tandem. In this paper, we introduce a Time-To-Live based policy (TTL), that assigns a timer to each content stored in the cache and redraws the timer each time the content is requested (at each hit/miss). We show that our TTL policy is more general than LRU, FIFO or RND, since it is able to mimic their behavior under an appropriate choice of its parameters. Moreover, the analysis of networks of TTL-based caches appears simpler not only under the Independent Reference Model (IRM, on which many existing results rely) but also with the Renewal Model for requests. In particular, we determine exact formulas for the performance metrics of interest for a linear network and a tree network with one root cache and N leaf caches. For more general networks, we propose an approximate solution with the relative errors smaller than 10 -3 and 10 -2 for exponentially distributed and constant TTLs respectively.

INTRODUCTION

Caches are widely used in networks and distributed systems for improving performance. They are integral components of the Web [START_REF] Che | Hierarchical Web caching systems: modeling, design and experimental results[END_REF], DNS [START_REF] Jung | DNS performance and the effectiveness of caching[END_REF], and Content Distribution Networks (CDNs) [START_REF] Saroiu | An analysis of Internet content delivery systems[END_REF]. More recently there has been a growing emphasis on content networks (i.e. CCNs) where content or data is centric and host-to-content interaction is the common case [START_REF] Jacobson | Networking named content[END_REF]. Many of these systems give rise to hierarchical (tree) cache topologies and to even more general irregular topologies. The design, configuration, and analysis of these cache systems pose significant challenges. An abundant literature exists on the performance (e.g. hit probability, search cost) of a single cache running the LRU replacement policy or, its companion, the Move-to-Front (MTF) policy (see [START_REF] Cabe | On serial files with relocable records[END_REF][START_REF] King | Analysis of demand paging algorithm[END_REF][START_REF] Hendricks | The stationary distribution of an interesting Markov chain[END_REF][START_REF] Burville | On a model for storage and search[END_REF][START_REF] Bitner | Heuristics that monotonically organize data structures[END_REF][START_REF] Flajolet | Birthday paradox, coupon collectors, caching algorithms and self-organizing search[END_REF][START_REF] Dan | An approximate analysis of the LRU and FIFO buffer replacement schemes[END_REF][START_REF] Fill | Limits and rate of convergence for the distribution of search cost under the move-to-front rule[END_REF][START_REF] Jelenkovic | Asymptotic approximation of the move-to-front search cost distribution and least-recently used caching fault probabilities[END_REF] for i.i.d. requests and [START_REF] Coffman | Performance of the Move-to-Front algorithm with Markov-modulated request sequences[END_REF][START_REF] Jelenkovic | Critical sizing of LRU caches with dependent requests[END_REF][START_REF] Jelenkovic | Least-recently used caching with dependent requests[END_REF] for correlated requests). With few exceptions, exact models of even single caches are computationally intractable, resulting in the reliance on approximations [START_REF] Dan | An approximate analysis of the LRU and FIFO buffer replacement schemes[END_REF][START_REF] Jelenkovic | Asymptotic approximation of the move-to-front search cost distribution and least-recently used caching fault probabilities[END_REF]. Networks of caches are significantly more difficult to analyze and no exact solution has been obtained so far for even the simple configuration of two LRU caches in series. A few approximations have been proposed, instead, for a simple two-level LRU cache network [START_REF] Che | Hierarchical Web caching systems: modeling, design and experimental results[END_REF] and a general LRU network [START_REF] Rosensweig | Approximate models for general cache networks[END_REF]. However, their inaccuracies can be significant as reported in [START_REF] Rosensweig | Approximate models for general cache networks[END_REF] where the relative error reaches 16%.

In this paper, we focus on a class of caches referred to as Time-To-Live (TTL) caches. When an uncached data is brought back into the cache due to a cache miss, a local TTL is set. The TTL value can be different for different data, but also for the same data at different caches 1 . All requests to that data before the expiration of the TTL are successful (cache hit); the first request for that data to arrive after the TTL expiration will yield a cache miss. In that case, the cache may forward the request to a higher-level cache, if any, or to the server. When located, the data is routed on the reverse-path and a copy is placed in each cache along the path. We strongly support the idea that a TTL policy can be an interesting alternative to more popular schemes like LRU or FIFO for two reasons. First, a TTL policy is more configurable and in particular can mimic the behavior of other replacement policies for an opportune choice of its parameters (the timeout values, see Section 3.2). Second, networks of TTL caches appear to be simpler to study analytically, while networks of LRU or FIFO caches have defeated until now modeling attempts. As we shall show, our TTL policy is general enough and fully configurable, and it appears to be a suitable unified framework for the performance analysis of heterogeneous caching networks where caches may run under different replacement policies.

We develop in particular a set of building blocks for the performance evaluation of hierarchical TTL cache networks with general topologies and where TTLs are set with every request. These blocks allow one to model exogenous requests at different caches as independent renewal processes and to allow for TTL duration to be described by an arbitrary distribution as long as they are independent of each other. The building blocks consist of:

• a renewal theoretic model of a single content TTL cache when fed by a renewal request stream,

• a renewal process approximation of the superposition of independent renewal processes.

The first block forms the basis for calculating cache metrics such as miss and hit probabilities/rates but also for characterizing the output sequence of requests (the miss process) of a cache, while the second block is used to characterize the resulting process of the superposition of exogenous requests and those resulting from a miss process if any. We apply these blocks to the case when TTL values are either exponentially distributed or constant. Then, we focus primarily on linear TTL networks, two level TTL tree networks and combinations of the two. We derive exact results for some cases but when they are not, the relative errors are extremely small and less than 10 -3 in the case of exponentially distributed TTLs and 10 -2 in the case of constant TTLs. Thus, we believe our approach is promising and capable of accurately modeling a richer class of network topologies.

In the literature, the paper closer to our approach is [START_REF] Jung | Modeling TTL-based Internet Caches[END_REF],

where the authors consider a single TTL-based cache fed by i.i.d. requests to a single data. They obtain the hit rate for a constant TTL via the solution of some renewal equation. Despite the increasing interest in CCNs, previous work has mainly focused on global architecture design. [START_REF] Carofiglio | Modeling data transfer in content-centric networking[END_REF] is probably the first attempt to model data transfer in CCNs; the authors develop approximations to calculate the stationary throughput in a network of LRU caches taking into account the interplay between receiver driven transport and per-chunk caching.

The paper is organized as follows. In Section 2, we introduce notation, the model assumptions, and a key result from [START_REF] Lawrance | Dependency of intervals between events in superposition processes[END_REF] regarding the computation of the marginal inter-arrival distribution for a superposition of arrival streams modeled as renewal processes. Section 3 contains our renewal theoretic model along with its application to some network topologies where it leads to exact results. It also shows how the TTL policy can mimic existing policies. We describe in section 4 our approach to model the combined exogenous and miss streams as a renewal request process; and the resulting approximations for a larger class of linear and tree networks. We also report the accuracy of this approximation. Section 7 shows how the TTL policy behaves under finite size buffer. Conclusions are found in Section 8.

Due to space constraints, all the proofs have been omitted as well as a discussion of the computational complexity of our analytic approach, but they can be found in the companion technical report [START_REF] Choungmo Fofack | Analysis of TTL-based Cache Networks[END_REF] together with more examples and a more detailed validation section.

DEFINITIONS AND ASSUMPTIONS

Throughout this paper we mainly focus on particular instances of TTL cache tree networks with infinite buffers size. the words "node" and "cache" will be used interchangeably. Also, a cache will always be a TTL cache unless otherwise specified.

New requests for a data can be generated at any node of the network according to mutually independent renewal processes; these requests are referred to as exogenous requests or arrivals and the sequence of these exogenous request instants is called the exogenous request process. If upon the arrival of a new request the data is not present in the cache, the request is instantaneously forwarded to the next level of the tree and the process repeats itself until the data is found.

In case the data cannot be found along the path toward the root, the root retrieves it from a server. Once the data is found, either at a cache or at a server, a copy of it is instantaneously transmitted to each cache along the path between the cache where the data was found and the cache that issued the request 2 . A new TTL is set for each new copy of the data and the TTL is redrawn at the cache, if any, where the data was found (by convention, the TTL at the server is infinite). This is in contrast with the model in [START_REF] Jung | Modeling TTL-based Internet Caches[END_REF] where there is no TTL reset upon a cache hit. Resetting the TTL also at each cache hit increases the occupancy and the hit probability specially for popular contents (high request rate). This choice is motivated by the CCN paradigm of moving popular documents as close as possible to the users.

We define the miss process at a cache as the successive instants at which misses occur at this cache, namely, the times at which the data is requested and is not found in the cache.

(s) = R ∞ 0 e -st χ(t)dt (s > 0) denotes the Laplace Trans- form (LT) of χ(t).
The following result, taken from [22, Formula (4.1)], will be repeatedly used in this paper .

Theorem 2.1. The CDF R(t), of the inter-event times of the point process resulting from the superposition of K mutually independent renewal processes, is given by

R(t) = K X k=1 α k P K k=1 α k Rk (t) K Y j=1,j =k αj Z ∞ t Rj (u)du,
with R k (t) and α k > 0, the CDF of the inter-event times and the arrival rate of process k, respectively.

EXACT RESULTS

Single cache

We consider a single TTL cache. Requests arrive at the cache according to a renewal process. Without loss of generality, we assume that the first request arrives at time t = 0 and finds an empty cache. We denote by X a generic interarrival time with CDF F (t) and density f (t). We also denote by T a generic TTL duration, with CDF T (t). Since successive inter-arrival times and successive TTLs form two independent renewal sequences and since a miss triggers a new TTL, miss times are regeneration points of the state of the cache. This implies that inter-miss times form a renewal process, with generic inter-miss time denoted by Y and CDF G(t). The stationary hit probability, hit rate and miss rate denoted by hP , hR and mR, respectively, are given by

hP = P (X ≤ T ) = Z ∞ 0 F (t)dT (t), (1)
hR = λhP , mR = 1/E[Y ] = λ(1 -hP ) (2)
respectively, with λ := 1/E[X] the arrival rate.

Propositions 3.1 and 3.2 (Proofs in [START_REF] Choungmo Fofack | Analysis of TTL-based Cache Networks[END_REF]) respectively give the cache occupancy π and the CDF G(t) for the most general result (with arbitrary CDFs F (t) and T (t)).

Proposition 3.1 (Stationary cache occupancy).

π := λE »Z X 0 (1 -T (t))dt - . (3)
Proof. Let V be the time during which the document is in the cache between two consecutive request arrivals. We have

π = E[V ]/E[X] = λE[V ] by renewal theory. Let us find E[V ].
Define the binary rv U (t) to be one if the document is in the cache at time t and zero otherwise. Without loss of generality consider the interval [0, X] corresponding to the inter-arrival time between the first and the second request.

We have

E[V ] = E »Z X 0 U (t)dt - = EX »Z X 0 E[U (t)|X]dt - = = EX »Z X 0 (1 -T (t))dt - where the last equality follows from E[U (t)|X] = E[U (t)] = P (U (t) = 1) = P (T > t). ⋄ Proposition 3.2.
. The CDF G(t) of inter-miss times is the unique bounded solution of the integral equation

G(t) = Z t 0 G(t -x) T (x)dF (x) + Z t 0 T (x)dF (x). ( 4
)
Proof. Let X1 (resp. Y1, T1) denote the first inter-arrival time (resp. first inter-miss time, first TTL) after t = 0. Since Y1 ≥ X1, the event {Y1 < t} may only occur if X1 < t. Therefore, [START_REF] Choungmo Fofack | Analysis of TTL-based Cache Networks[END_REF] where [START_REF] Che | Hierarchical Web caching systems: modeling, design and experimental results[END_REF] follows from the fact that the event {Y1 < t} is true when T1 < X1 < t. It remains to evaluate the probability P (Y1 < t, X1 < t, X1 ≤ T1) in [START_REF] Choungmo Fofack | Analysis of TTL-based Cache Networks[END_REF]. By conditioning on X1 and T1 we obtain

G(t) = P (Y1 < t, X1 < t, X1 ≤ T1) +P (Y1 < t, T1 < X1 < t) = P (Y1 < t, X1 < t, X1 ≤ T1) + P (T1 < X1 < t) (5) = P (Y1 < t, X1 < t, X1 ≤ T1) + Z t 0 T (x)dF (x)
P (Y1 < t, X1 < t, X1 ≤ T1) = (7) Z t x=0 Z ∞ τ =x G(t -x)dF (x)dT (τ ) = Z t 0 G(t -x)(1 -T (x))dF (x), (8)
where the first equality is due to the fact that the TTL is renewed at each request and then Y1 -X1 conditioned to X1 ≤ T1 has the same distribution of Y1.

Suppose that there are two solutions G1(t) and G2(t) satisfying (4). Then

G1(t) -G2(t) = R t 0 (G1(t) -G2(t))(1 - T (x))dF (x)
. By Laplace transforming both sides of this equality, it appears evident that G * 1 (s) -G * 2 (s) = 0 and then the solution is unique. ⋄

If TTLs are exponentially distributed with rate µ, then

π = λ(1 -F * (µ)) µ (9)
from [START_REF] Cabe | On serial files with relocable records[END_REF]. Given that T (t) = 1e -µt , hP = F * (µ) from [START_REF] Bitner | Heuristics that monotonically organize data structures[END_REF], which in turn implies from (2) that hR = λF * (µ) and mR = λ(1 -F * (µ)). Taking the Laplace transform of both sides of (4) yields

G * (s) = F * (s) -F * (s + µ) 1 -F * (s + µ) . (10)
We can check from ( 10) that mR = -(dG * (s)/ds|s=0) -1 .

If the TTL is a constant, equal to T , (4) becomes

G(t) = Z t∧T 0 G(t -x)dF (x) + (F (t) -F (T ))1(t > T ) (11)
with a ∧ b = min(a, b). In this case the hit probability is F (T ).

Relation with other replacement policies

Consider a single cache with capacity B and M contents, whose request processes are independent Poisson processes with different rates λ k for k = 1, 2, . . . , M . We shall tune the TTL policy (i.e. the timer rate µ k for the k-th content) in order to get the same performance metrics of common replacement policies like LRU, FIFO or RND. Without loss of generality, we consider a LRU cache and we denote by p k the corresponding stationary cache occupancy for content k. This distribution and the corresponding one for FIFO have been calculated in [START_REF] King | Analysis of demand paging algorithm[END_REF]. For the exponentially distributed TTL cache, the stationary occupancy of the k-th content is given by

π k = λ k (1 -F * k (µ k ))/µ k , where F * k (s) = λ k /(λ k +s).
It appears evident that, for the infinite capacity TTL cache we can always select µ k = λ k (1/p k -1) such that π k = p k . However, under storage constraints like finite capacity B, we can still take µ k = λ k (1/p k -1) such that P k p k = B = P k π k as we shall see in Sec.7 for the practical TTL policy. From the equality of the stationary cache occupancies, the equality of hit/miss probabilities and rates follows thanks to the PASTA property of request processes. In this sense, the TTL policy is more general than LRU and FIFO, since it can mimic their behavior3 .

Line of caches

Consider the linear cache network in Fig. 1a composed of N TTL caches labeled 1, . . . , N , without exogenous requests at caches 2, . . . , N . Requests arrive to cache 1 according to a renewal process with generic inter-arrival time X and arrival rate λ. TTLs at all caches are mutually independent and exponentially distributed random variables with rate µn at cache n.

The arrival process at cache n is the miss process at cache n-1 since there are no exogenous arrivals. Moreover, it is easily seen that the miss times at cache n -1 form regeneration points, so that the miss process at this cache, and therefore the arrival process at cache n, is a renewal process. We can then apply recursively the results obtained for a single cache. In particular, if we denote by G * n (s) the LST of the intermiss times at cache n, we may apply formula [START_REF] Flajolet | Birthday paradox, coupon collectors, caching algorithms and self-organizing search[END_REF] where the LST of the inter-arrival times is G * n-1 (s). We obtain

G * n (s) = G * n-1 (s) -G * n-1 (s + µn) 1 -G * n-1 (s + µn) (12)
for n = 1, . . . , N , where G * 0 (s) = F * (s). At the cache n the hit probability hP,n = G * n-1 (µn), the miss rate mR,n = mR,n-1(1 -hP,n) and hit rate hR,n = mR,n-1hP,n are derived by using ( 12)

mR,n = λ n-1 Y i=0 (1 -G * i (µi+1)), ( 13
)
hR,n = λ n-2 Y i=0 (1 -G * i (µi+1)) G * n-1 (µn), (14)
with mR,0 := λ. The occupancy of cache n is given by applying (3) with F (t) = Gn-1(t) and λ = mR,n-1, hence

πn = mR,n-1 " 1 -G * n-1 (µn) µn « . (15)

Simple tree

Consider the tree network in Fig. 1b with one root (labeled N + 1) and N children (leaves) labeled n = 1, . . . , N . The exogenous requests arrive at the n-th node according to a Poisson process with rate λn. While the TTL is exponentially distributed with rate µn at the n-th leave, we assume an arbitrary CDF TN+1(t), with LST T * N+1 (s) for the TTL at the root. Using the results in Section 3.3 to study caches n = 1, . . . N in isolation, and Theorem 2.1 to calculate the CDF of the overall request inter-arrival time at cache N + 1, we can derive all metrics of interest [START_REF] Choungmo Fofack | Analysis of TTL-based Cache Networks[END_REF].

APPROXIMATE RESULTS

The exact results in Section 3 cannot be easily extended to general networks. In fact, in the presence of exogenous requests, the aggregated (overall) arrival process at a cache is not a renewal process. Hence, we cannot apply Proposition 3.2 that allows us to characterize the miss requests process. However, we can still determine all performance metrics (e.g: see Sec. 3.4) at the cache by calculating the CDF of the inter-arrival times (Theorem 2.1). In this section we develop an approximation method to overtake this limitation. The solution derived produces highly accurate approximations under more general networks for all metrics considered in Section 3 (hit/miss probabilities, hit/miss rate, cache occupancy). The quality of the approximation is assessed in [START_REF] Choungmo Fofack | Analysis of TTL-based Cache Networks[END_REF] and it is based on the following statement: Approximation A1: the overall arrival process at each node is a renewal process.

As a consequence of A1 and Prop. 3.2, we approximate the miss process at a node by a renewal process. With a slight abuse of notation we will use the notation of Section 3 to denote the corresponding approximate values calculated under Approximation A1. For example, Hn(t) is used to denote the approximate CDF of the overall inter-arrival times and similarly Gn(t), Λn, mR,n, hP,n and hR,n are used to denote approximate quantities at node n. Regarding the total rate Λn, note that

Λn = λn + X i∈C(n) mR,i (16)
where C(n) is the set of children of node n. As in Section 3, exogenous arrivals at each node n form a renewal process, with CDF Fn(t). Thanks to A1, we invoke Theorem 2.1 to get

Hn(t) = 1 - λn Λn Fn(t) Y i∈C(n) νi Z ∞ t Ḡi(u)du - X i∈C(n) νi Λn Ḡi(t)λn Z ∞ t Fn(u)du × Y j∈C(n) j =i νj Z ∞ t Ḡj (u)du. ( 17
)
An approximation of the CDF of the inter-miss times at cache n is obtained from Proposition 3.2

Gn(t) = Z t 0 Gn(t -x) Tn(x)dHn(x) + Z t 0 Tn(x)dHn(x) (18)
where Tn(t) the CDF of the TTL duration at cache n. Eqs ( 17)-( 18) provide a recursive procedure for calculating, at least numerically, approximations of the CDFs Gn(t) and Hn(t) for each cache n of a general network topology, from which we can derive approximate formulas for all metrics.

However, even for small networks the numerical complexity of this procedure can be very high as it requires calculating integrals over infinite ranges (see Eq. ( 17)) and solving integral equations (see Eq. ( 18)).

In order to show that this complexity can be significantly reduced, we focus on a particular class of tree networks, class N and we give explicit results. The TTL at each node n is exponentially distributed with rate µn. A network belongs to class N if, in addition to A1, the following approximation holds:

Approximation A2: The node n is fed by the superposition of two independent request arrival processes: one (stream 1) is the miss rate of a child of cache n and is a generic renewal process and the other one (stream 2) is a renewal process with CDF of the form

Kn(t) = 1 - Mn X m=0 αn,me -βn,mt (19)
where 0 ≤ Mn < ∞ and {βn,m}m is a set of non negative numbers.

In what follows we assume without loss of generality that stream 1 originates from a cache child labeled n -1 and then we denote the CDF of the inter-miss times in stream 1 as Gn-1(t) and the miss rate as νn-1. From ( 19) the arrival rate of stream 2 is ηn := P Mn m=0 αn,mβn,m, the total arrival rate at node n is Λn = νn-1 + ηn. Approximations A1 and A2 together yield the following procedure for approximating G * n (s) and H * n (s).

-s 2 ηnνn-1 Λn

Mn X m=0 αn,m(1 -G * n-1 (s + βn,m)) (s + βn,m) 2 βn,m and G * n (s) = H * n (s) -H * n (s + µn) 1 -H * n (s + µn) . (21)
Proof.

Hn(t) = 1 - ηnνn-1 Λn Mn X m=0 αn,m × Ḡn-1(t) Z ∞ t e -βn,mu du + e -βn,mt Z ∞ t Ḡn-1(u)du
! from which we deduce [START_REF] Laoutaris | Meta Algorithms for Hierarchical Web Caches[END_REF]. ( 21) is obtained from [START_REF] Flajolet | Birthday paradox, coupon collectors, caching algorithms and self-organizing search[END_REF].

⋄ The hit probability is simply hP,n = H * n (µn), the total arrival rate Λn = νn-1 + ηn and the miss rate νn = Λn(1 -

H * n (µn)) are νn = n X i=1 ηi n Y j=i (1 -H * j (µj)), (22) Λn
= n-1 X i=1 ηi n-1 Y j=i (1 -H * j (µj)) + ηn, (23)
Relations ( 20)-( 21) and ( 22)-( 23) provide a recursive procedure for calculating Λn and H * n (µn) for each n, from which we obtain approximations for the hit probability, hit rate, miss rate and stationary occupancy at node n:

hP,n = H * n (µn), hR,n = ΛnH * n (µn), ( 24
) mR,n = Λn(1 -H * n (µn)), πn = Λn " 1 -H * n (µn) µn « .
The latter result follows from [START_REF] Dobrow | The move-to-front rule for self-organizing lists with Markov dependent requests[END_REF].

In [START_REF] Choungmo Fofack | Analysis of TTL-based Cache Networks[END_REF], we show that the network in Fig. 2 belongs to the class N , when TTLs are exponentially distributed and exogenous request processes are Poisson.

VALIDATION

In this section we investigate the accuracy of the approximation method developed in Section 4. Recall that the method consists in assuming that all internal arrival processes at a node (i.e. processes formed of the miss processes of the 17) to calculate the CDF of the inter-arrival times of the superposed process. The miss process at this node can then be characterized by using Eq. ( 18) and the procedure is repeated at the node's parent.

We focus our validation on the case when the TTLs are exponentially distributed, but we also provide some results for constant TTLs.

Exponential timers

We and m M R,2 = λ1p(0, 0, * ) + λ2(p(0, 0, * ) + p(1, 0, * )), where p(i, * ) = P x 2 ,...,x N ∈{0,1} p(i, x2, . . . , xN ) and p(i, j, * ) := P x 3 ,...,x N ∈{0,1} p(i, j, x3, . . . , xN ) are the stationary probabilities that cache 1 is in state i ∈ {0, 1} and caches [START_REF] Bitner | Heuristics that monotonically organize data structures[END_REF][START_REF] Burville | On a model for storage and search[END_REF] are in state (i, j) ∈ {0, 1} 2 , respectively. Due to space constraints we omit the general expressions for these quantities for a generic cache in the line and those for a line of simple tree networks that can be similarly calculated.

= P x∈E,x i =1 p(x) (the superscript "M" stands for "Markov") .
In the rest of this section we compare our approximate results versus the exact ones that can be obtained studying the Markov process. A comparison of the computational costs of the two approaches is in Section 6. We consider first the line network in Fig. 3: it has four nodes (N = 4), exogenous arrivals and exponentially distributed TTLs at node n with rate λn and µn, respectively. We have calculated Figure 3: Line of four caches the absolute relative errors at cache n for the hit probability (EHP,n), the miss rate (EMR,n) and the occupancy probability (EOP,n). The exact value is calculated through the analysis of the Markov process, e.g. EHP,n := |h M P,n -hP,n|/h M P,n . Fig. 4 shows the CDFs of the relative errors at cache 4 for 1001 different parameter vectors ((λn, µn), n = 1, . . . , 4). The values of the exogeneous arrival rates (resp. TTL rates) have been selected in the interval [0.001, 10] (resp. [0.1, 2]) according to the FAST (Fourier Amplitude Sensitivity Test) method (see [21, Sec. VI-C] and references therein). We can observe that the approximation is very accurate: in 99% of the different parameter settings the relative error is smaller than 2 × 10 -5 . Fig. 5 shows how the error changes for different request loads. In this case we have considered the homogeneous scenario where all the caches have the same TTL and the same exogenous arrival rate, i.e. µn = µ and λn = λ for each n. The error is shown as a function of the normalized load ρ = λ/µ. We can observe that the largest error (about 2×10 -4 ) is obtained when arrival rates and timer rates have comparable values (ρ ≈ 1). In this case the different request processes superposed at a node have similar time scales and 6 then the inter-arrival times of the overall request process are more correlated (see also comments below).

We have also investigated the accuracy of the approximation for the line of simple tree networks (defined in Section 3.4) shown in Fig. 6: where nodes 11, 12, 13 (resp. nodes 21, 22, 23) have identical arrival rates and identical TTLs rates. Since the approximation results are exact for all nodes but node 2 we only report results for that node. The empirical CDFs of EHP,2, EMR,2 and EOP,2 are shown in Fig. 7. Like for Fig. 4 exact results have been obtained by considering the Markov process associated with this line of simple trees network. Different request and TTL rates have been selected according to the FAST method respectively in the intervals [0.001, 10] and [0.1, 2]. We used 4921 samples for each rate. Results are analogous to those for a line of caches. The relative errors can be larger in this scenario, but they are probably negligible for most of the applications (< 3 × 10 -4 in 99% of the cases). We have also considered the homogenous scenario also for this topology, the relative errors have the same order of magnitude (< 10 -3 ).

We have shown that Assumption A1 leads to very accurate results when exogenous arrival processes are Poisson and TTL are exponentially distributed. This let us think that the superposition of the request arrival processes at every cache is very 'close' to a renewal process. In order to justify such statement, we have calculated the first autocorrelation lag (r1) for the actual arrival process at node 2 in Fig. 3 using Eq. (6.4) in [START_REF] Lawrance | Dependency of intervals between events in superposition processes[END_REF]. This autocorrelation lag depends on the Figure 8: Tree network arrival rates λ1 and λ2 and the timer µ1. We have found that for any possible choice of these parameters 0 > r1 > -0.015. Simulation results show that the autocorrelation is even less significant at larger lags. We can then conclude that the inter-arrival times are weakly coupled.

Deterministic Timers

When timers are deterministic we need to rely on the general procedure described in Section 4 and based on Eqs [START_REF] Jung | DNS performance and the effectiveness of caching[END_REF] and [START_REF] King | Analysis of demand paging algorithm[END_REF]. There are two sources of errors in this procedure. Firstly, the aggregate request process at a cache is not a renewal process and it is not correct to apply the renewal equation [START_REF] King | Analysis of demand paging algorithm[END_REF]. Secondly, both the steps ( 17) and ( 18) introduce some numerical errors. Two parameters determine the entity of the numerical error: 1) the time interval (τ ) from which the CDF samples are taken, 2) the time distance between two consecutive samples (∆). Clearly the larger τ and the smaller ∆ the smaller the numerical error, but also the larger the computational cost.

We have implemented a Matlab numerical solver that iteratively determines the CDFs in the network as described above, and then the metrics of interest for a cache network. The integrals appearing in Eqs. ( 17) and ( 18) are approximated as simple sums and for simplicity the same values τ and ∆ have been considered for all the CDFs. These parameters are selected as follows: our solver first obtains an approximated solution for the whole network assuming that all the request processes are Poisson and set the parameter τ to 5 times the largest expected inter-arrival time in the network. The parameter ∆ is set to one thousandth of the minimum of the TTL values and the expected interarrival times of the exogenous request processes.

We present some preliminary results for the tree network of Fig. 8 in the following situation: infinite capacity caches, asymmetric traffic conditions and non-exponential (i.e constant) life times. More specifically, the exogenous request processes are Poisson processes with rates λi (i = 5, 6, 7, 8, 9) and TTL values are Ti (i = 1, 2, . . . 9). In order to evaluate the relative error of the estimated metrics, we have considered as correct values those obtained through a long simulation. For example if our method predicts the value hP,n for the hit probability rate at node n and the 99% confidence interval, calculated by simulation, is [h S P,nǫ, h S P,n + ǫ], the relative error is calculated as |hP,n-h S P,n |/h S P,n . The relative incertitude of the simulation (ǫ/h S P,n ) is at most 0.3 × 10 -4 . For all the performance metrics and all the caches the relative error of our approach is less than 10 -2 .

COMPUTATIONAL COST

In this section we perform a preliminary analysis of the computational cost of our approach.

We first address the case of a class N network, and in particular we consider a line of simple tree networks with N trees and M nodes in total as in Fig. 2. Since the computational cost for all the metrics is roughly the same, we focus here on the hit probability. In order to calculate the hit probability at one of the roots of the simple trees, say it cache n, we need to evaluate the LST H * n (µn) (Eq. ( 24)). This requires a number of operations proportional to the number of children of cache n (Rn) and the evaluation of the LST of the miss rate coming from cache n -1 in µn, i.e G * n-1 (µn) (Eq. ( 20)). In turn G * n-1 (µn) can be calculated evaluating H * n-1 (s) in two points (µn and µn + µn-1) (Eq. ( 21)) and so on recursively. This implies that the cost to calculate the hit probability at cache n is O(αRn + 2 n ) for some constant α. When evaluating the hit probability at other caches the same LSTs needs to be evaluated, but in general at different points, then we have that the total cost is O(

P N n=1 αRn + 2 n ) = O(αM + 2 N
). Then, depending on the topology of the network, the cost can be mainly linear in the number of nodes (for a network with small depth, e.g. when there are a few trees each with a lot nodes) or exponential in the number of nodes (for a network with large depth, e.g. for the linear network in Fig. 1).

It is interesting to compare this cost with alternative approaches. For the line of simple tree networks, all the metrics can be exactly calculated solving a Markov process as we mentioned in Section 5. The size of the state space is 2 M , then the cost of determining the steady-state distribution by solving the linear equation system is O(2 3M ) and this is much larger than the cost of our method O(αM + 2 N ). A different approach is to obtain an approximated steadystate distribution of the Markov process using an iterative method. This approach takes advantage of the fact that most of the transition rates have value zero. In fact a state change is triggered by an exogenous request arrival at a cache that does not have the data or by a timer expiration at a cache with the data, i.e. from a given state we can only reach other M states. Then the number of non-zero rates is equal to M 2 M and each iteration of the method requires O(M 2 M ) operations. The total cost of the iterative method is then O(KM 2 M ), where K is the number of iterations until termination and depends on the spectral gap of the matrix used at each iteration and on the required precision, but in general we can expect O(KM 2 M ) << O(2 3M ). Assuming that this is the case, we can observe that our method, even in the worst case of the linear network, is still more convenient than solving the Markov process, because O(2 M ) < O(KM 2 M ). 9 shows the ratio of the computation times to calculate our approximation (T A ) and to solve the Markov chain (T M ) for a line of N caches (with N = 1, 2, . . . 9). Both the methods have been implemented in Matlab, in particular the function linsolve has been used to determine the steady-state distribution of the Markov chain.

Let us now consider the case of a general tree network with constant TTLs (equal to T ). In this case there is no exact solution to compare our approach with, so we consider simulations as an alternative approach. We perform an asymptotic analysis. A meaningful comparison of the computational costs needs to take also into account the incertitude of the solution: both the simulations and our method can produce a better result if one is willing to afford a higher number of operations. In order to combine these two aspects in our analysis we consider as metric the product precision times number of operations. Intuitively the larger this product the more expensive is to get a given precision. For the simulations the computational cost is at least proportional to the number of events that are generated, let us denote it by nE. The incertitude on the final result can be estimated by the amplitude of the confidence interval, that decreases as 1/

√ nE, then the product precision times number of operations is proportional to √ nE for the simulations. In the case of our approach, the heaviest operation is the solution of the renewal equation. If we adopt the same τ and ∆ for all the integrals, we need to calculate the value of the CDF of the miss rate (G(t)) in nP = τ /∆ points and then we need to calculate nP integrals. The integration interval is at most equal to the TTL duration T (see Eq. 11), then each integral requires a number of operations proportional to n ′ P = T /∆. If the value of τ is selected proportionally to T , then the cost of our method is proportional to n 2 P . A naive implementation of the integral as a sum of the function values leads to an error proportional to the amplitude of the time step and then inversely proportional to n ′ P or nP . In conclusion the product precision times the number of operations is proportional to nP . Then, for a given precision, our method would require a number of points much larger than the number of events to be considered in the corresponding simulation (at least asymptotically). The comparison would then lead to prefer the simulations at least when small incertitude is required (then large nE and nP ). In reality integrals can be calculated in more sophisticated ways, for example if we adopt Romberg's method, with a slightly larger computation cost, we can get a precision proportional to n -2 P . In this case the product precision times number of operations is a constant for our method, that should be preferred.

A PRACTICAL TTL POLICY

While the model considered above allows for an arbitrarily large number of contents, a real cache will have a finite capacity B. In this section, we consider a possible practical implementation (TTL Impl) of the ideal TTL policy (TTL Model) we have studied above. The cache uses timers for each content as described above but it does not discard contents whose timer has expired as long as some space in the buffer is available. If a new content needs to be stored and the cache is full the content to be erased is the content whose timer has expired since a longer time (if any) or the content whose timer is going to expire sooner.

We have compared the performance of the TTL Impl with that of our TTL Model (in Sec. 3.3) for a line of two caches with same capacities B1 = B2 = 20 where the requests for each content k = 1, . . . , M = 200 arrive only at the first cache. The request processes are independent Poisson processes with rates distributed according to a Zipf law with parameter α = 1.2. The timers are i.i.d. exponential random variables with rate µ1 and µ2 respectively at the two caches. Each timer rate has been selected so that the expected buffer occupancy for both the TTL Impl and the TTL Model equals the corresponding cache capacity. In other words, µn at cache n = 1, 2 is chosen such that P M k=1 π n,k = Bn. Results show that the performances of the two policies are very similar. The aggregate hit probabilities at the two caches are respectively hP,1 = 0.5584 and hP,2 = 0.4663 for the TTL Model and ĥP,1 = 0.5618 and ĥP,2 = 0.4124 for the TTL Impl. At the content granularity, the matching is even accurate at cache 1, so that we just report in Fig. 10 the results at cache 2. These preliminary results suggest that our analysis can be useful to study TTL policies when caches are constrained.

CONCLUSION

In this paper, we introduced a novel Time-To-Live replacement policy for buffering routers of ICN and we have developed a set of building blocks for the performance evaluation of TTL cache networks based on simple renewal arguments. We derived exact results and closed form formulas to assess the metrics of some hierarchical caching networks like the linear and the simple tree networks. We also provided a recursive and approximate procedure to study a general network topology. For a large class of tree networks with exponentially distributed and deterministic TTLs, we showed that our approximation is highly accurate with relative errors of 10 -3 and 10 -2 respectively. Our approach scales easily and it is promising since it appears as a simple unifying model for accurately modeling a richer class of networks such as heterogeneous caching networks under different replacement policies. We have also demonstrated that our TTL model can be implemented and used to optimize a multi-content cache network under realistic constraints such as the cache size limitation. Ongoing research is investigating the effect of finite capacity cache networks, how to reduce the computational cost of our approach in the case of constant TTLs and how to take into account correlated request processes by using semi-Markov processes.

Figure 1 :

 1 Figure 1: Network architectures

Proposition 4 . 1 (

 41 Approximation for class N ). Under Approximations A1 and A2, for each node n,

Figure 2 :

 2 Figure 2: Line of simple tree networks

  For a line of caches the hit probability and the miss rate at cache 1 are respectively h M P,1 (1) = p(1, * ) and m M R,1 = λ1p(0, * ), while for cache 2 it holds h M P,2 = λ1p(0, 1, * ) + λ2(p(0, 1, * ) + p(1, 1, * )) λ1(p(0, 0, * ) + p(0, 1, * )) + λ2

4 CDFFigure 4 :

 44 Figure 4: CDF of EHP,4, EMR,4, EOP,4 for network in Fig.3

Figure 5 :

 5 Figure 5: EHP,4, EMR,4, EOP,4 for network in Fig.3 with homogeneous nodes (λn = λ = ρµ = ρµn)

Figure 6 :Figure 7 :

 67 Figure 6: Linear tree network

Figure 9 :

 9 Figure 9: Running time comparison

Figure 10 :

 10 Figure 10: TTL-policy model and its implementation

Table 1 :

 1 Glossary of main notation

	This key hypothesis allow us to decouple the management
	of the different contents and study each of them separately.
	For this reason, in what follows we will simply refer to a
	single content or data chunk (simply called the data). The
	effect of finite buffer is considered in Section 7. From now on

  Let us denote by C(n) the set of children of cache n. The (overall) request process, also called the arrival process, at cache n is the superposition of the miss processes of caches in C(n) and of the exogenous request process at cache n,

	if any. We assume that successive TTLs at each cache are
	i.i.d. random variables, that TTLs at different caches are
	mutually independent, and that all TTLs are independent
	of the exogenous arrivals. If Λn is the arrival rate of requests
	at cache n, hP,n (resp. mP,n = 1 -hP,n) and hR,n = ΛnhP,n (resp. mR,n = Λn(1-hP,n)) denote the stationary hit (resp.
	miss) probability and the stationary hit (resp. miss) rate at
	cache n, respectively. We denote by πn the steady-state
	probability that the data is in cache n and we call it the
	occupancy of cache n. Hence, we just have to calculate hP,n
	and Λn.
	For t ≥ 0 and any non-negative random variable X with Cumulative Distribution Function (CDF) χ(t) = P (X < t),
	χ(t) = 1 -χ(t) is the Complementary Cumulative Distri-bution Function (CCDF), χ ⋆ (s) = E[e -sX ] = R ∞ 0 e -st dχ(t) (s ≥ 0) denotes its Laplace-Stieltjes Transform (LST), and χ

This is then different from the TTLs in DNS system where the TTL for a given data is in general set to a common value determined by the authoritative name server.

We observe that our TTL-based policy could be used also in conjunction with other strategies that do not keep a copy at every cache along a path, like those in[START_REF] Laoutaris | Meta Algorithms for Hierarchical Web Caches[END_REF].

We observe that the inter-miss time distributions are very different if timers are exponential, but if we select different TTL distributions we could match also them.

ACKNOWLEDGEMENTS

This research has been supported by Orange Labs in the framework of the contract "Performance models for CCN architecture". We thank in particular Orange Labs researchers Bruno Kauffmann, Luca Muscariello and Alain Simonian for the stimulating discussions and for the valuable feedback.