Yvon Maday
email: maday@ann.jussieu.fr

Noura Morcos
email: morcos@ann.jussieu.fr

AND Toni Sayah
email: tsayah@fs.usj.edu.lb

REDUCED BASIS NUMERICAL HOMOGENIZATION FOR SCALAR ELLIPTIC EQUATIONS WITH RANDOM COEFFICIENTS: APPLICATION TO BLOOD MICRO-CIRCULATION

Keywords: Porous media, stochastic homogenization, reduced basis method, microcirculation, blood flow AMS subject classifications. 74Q05, 74SQ5, 65N15, 35J20

We consider a non periodic homogenization model designed to simulate the blood flow at the level of the micro-vascularised tissues. We solve elliptic partial differential equations with two length-scales on the domain and we use the reduced-basis method to speed up the numerical resolution. Finally, we show numerical results and comparisons in 2D and 3D.

1. Introduction. Our goal in writing this paper is to develop a novel numerical technique for approximating the solution of second order scalar elliptic boundary value problems with highly oscillating coefficients amenable to take into account some randomness. Such mathematical models are involved in many applications, including thermal diffusion in porous media, constitutive laws of composite materials, behavior of biological tissues. The high oscillations are represented by a small parameter denoted by ε ∈ IR + that actually may represent several distinct scales, ε being only a notation for representing the small components that are present in the phenomenon. In what follows we only focus on the two scale problem (one macro scale and one micro scale). Solving directly such problems is generally quite expensive, even impossible. For instance, with plain finite element methods, since the discretization mesh should be small enough in order to represent well the oscillations meaning that the mesh size should be smaller than ε, the computation load would scale at least as O(ε) -3 in 3-dimensions, which is difficult to manage especially when ε → 0.

It should first be noticed that, for most of the practical applications, there is only one problem to be solved, corresponding to a unique value of ε say ε = ε 0 , however, if ε 0 is very tiny, with respect to the macroscopic size of the problem, it may be smart to invent a family of problems depending on ε and discover the asymptotic limit behavior of the solution as ε → 0, we then speak about homogenization theory and homogenized problem.

In such cases, the determination and the mathematical justification of the limit behavior is more or less known in some particular cases related to the dependency of the problem in ε, with a full rigorous theory that is available only in few cases. Recent advances have allowed to incorporate randomness in the structure of the oscillations improving the range of applications of these models. In order to go even further, numerical homogenization techniques have been introduced to cope with some current holes in the theory. This subject is very much actual and active due to the variety of applications that can be considered.

In a nutshell, most of the homogenization techniques, when two scales are separated, allow to replace the resolution of one problem with two scales into the resolution of two classes of coupled problems each dealing with one scale, in particular the fine scale problem depends on the large scale and its solution enters in the definition of the large scale limit problem. Numerical homogenization techniques then require the simulation of a huge number of micro scale problems. This paper builds up on the contribution by Boyaval [12] that proposes and illustrates the use of reduced basis methods for solving these micro scale problems. Reduced basis methods introduced in [START_REF] Fink | On the error behavior of the reduced basis technique for nonlinear finite element approximations[END_REF][START_REF] Noor | Reduced basis technique for nonlinear analysis of structures[END_REF], offer a particularly well suited solution to the challenges of this many-query framework, aiming to achieve efficiently the simulation of this huge amount of similar problems. A fundamental observation and assumption utilized by RBM is that the solutions at the micro scale level constitute a family of low dimension parameterized by the macro scale variable. Therefore, one can expect that every particular solution, at the micro scale level, can be well approximated by a finite and low dimensional vector space. In the context of the RBM, this low dimensional space is chosen as being spanned by (very well) approximated solutions to the same problem for very few, well chosen, macro scale spacial variables. In addition, we incorporate in this paper one important new feature: the stochastic dimension in the definition of the microscopic structure.

In the second part of this paper, we present the particular problem that actually has led us to develop the stochastic two scale homogenization: it is related to the modeling of blood flow at the level of the micro-vascularised tissues. Many example referring to homogenization exist in the literature describing such models of microcirculation that account for the highly oscillating media composed of microstructures. Biological tissues are certainly a place where randomness is illustrated and mathematical models have to account for this feature. The model we propose here enters in this stochastic homogenization frame and leads to methods for constructing the limit homogenized behavior. The stochastic frame and the particularity of the material involving components in which the flow does not exist, lead to some new specific difficulties and contributions (in terms of numerical analysis and limit behavior) that we handle here.

Basics of homogenization theory.

2.1. Abstract homogenization. Let D be an open piecewise C 1 bounded domain of IR d with a boundary Γ = ∂D. We denote by n the exterior unit normal to Γ defined a.e. and we assume that Γ is decomposed into two non overlapping open parts Γ D and Γ N . We are interested in the behavior of a sequence of scalar functions u ǫ that satisfy the following scalar elliptic equation in divergence form with random coefficients

-div(A ǫ (.; ω) ∇u ǫ (.; ω)) = f, in D, u ǫ | Γ D = g D ; A ǫ ∇u ǫ . -→ n | Γ N = g N ; (2.1)
for a sequence of scalars ǫ > 0. Here ω ∈ Ω, with (Ω, F, IP) being a probability space over Ω. We assume that, e.g., f ∈ L 2 (D), g D ∈ H 1/2 (Γ D ) and g N ∈ L 2 (Γ N ) and we look for the asymptotic limit of the sequence u ǫ (.; ω) when ǫ → 0. For every ǫ > 0 and almost surely in ω, we assume that A ǫ (.; ω) ∈ L ∞ (D, M α,β ), where M α,β is the set of uniformly positive definite square matrices with uniformly positive definite inverse, that is, matrices B satisfying, for all x ∈ D,

0 < α|ξ| 2 ≤< B(x)ξ, ξ >, ∀ξ ∈ IR d , 0 < β|ξ| 2 ≤< (B(x)) -1 ξ, ξ >, ∀ξ ∈ IR d . (2.2)
The problem (2.1) is almost surely well posed and for every ǫ > 0, there exists a unique solution u ǫ (.; ω) ∈ H 1 (D) verifying the relation

||u ǫ (., ω)|| H 1 (D) ≤ C(D, α, β) ||f || L 2 (D) + g D H 1/2 (Γ D ) + g N L 2 (Γ N ) .
Classical results about H-convergence (see e.g. [START_REF] Jikov | Homogenization of differential operators and integral functionals[END_REF], [START_REF] Murat | H-convergence[END_REF], [START_REF] Allaire | A multiscale finite element method for numerical homogenization[END_REF], [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF]) state that, under the above hypothesis, there exists a subsequence, still denoted as ε, and a homogenized matrix A * (., ω) ∈ L ∞ (D, M α,β ) such that, almost surely in ω, A ε (., ω) H-converges towards A * (., ω). Let us restrict ourselves momentarily to the case where Γ N = ∅ and g D = 0, since homogenization does not depend on the type of boundary conditions that are imposed (nor on f ). In this case, this means that, for ε → 0, the limit being true in

D ′ (D; IR n ) u ε (., ω) ⇀ u * (., ω), A ε ∇u ε (., ω) ⇀ A * ∇u * (., ω), (2.3)
where u * (., ω) is the solution to

-div(A * (.; ω) ∇u * (.; ω)) = f, in D, u * = 0 on ∂D. (2.4)
It should be noticed that, in many cases, the matrix A * (and thus u * ) is actually deterministic: this remarkable property of the matrix A * results from natural assumptions of stationarity and ergodicity. However, it is possible to build materials for which such a stochastic averaging effect does not hold.

Besides, except in some particular cases (e.g. periodic case, or in the stochastic frame of [START_REF] Papanicolaou | Boundary value problems with rapidly oscillating random coefficients[END_REF]) an explicit formula for the limit A * is not available which leaves a large room for investigation. Actually the result above goes only a little bit further by introducing the following abstract elements defined implicitly as the solutions ŵε i to the problems

-div(A ǫ (.; ω) ∇ ŵε i (.; ω)) = -div(A * (.; ω) e i ) in D, ŵε i (.; ω) = x i on ∂D, (2.5)
where ( e i ) 1≤i≤d denotes the canonical basis of IR d , and (x i ) 1≤i≤d the associated coordinates.

We then have the limit properties almost surely in ω

ŵε i (.; ω) ⇀ x i (2.6)
and

A ǫ (.; ω) ∇ ŵε i (.; ω) ⇀ A * e i . (2.7)
In addition (see e.g. [START_REF] Allaire | A multiscale finite element method for numerical homogenization[END_REF]), these solutions, ŵε i , named "oscillating test functions" allow to get the following corrector result:

∇u ε (.; ω) = d i-1 ∇ ŵε i (.; ω) ∂u * (.; ω) ∂x i + r ε (., ω), (2.8)
where r ε (., ω) converges strongly to zero in L 1 (D). It results that, whenever u * (.; ω) ∈ H 2 (D)

u ǫ (.; ω) = u ⋆ (.; ω) + d i=1 ( ŵε i (.; ω) -x i ) ∂u * (.; ω) ∂x i + R ǫ (., ω), (2.9)
where the remainder R ε (., ω) converges strongly to zero in W 1,1 (D). This powerful theory is somehow weakened for the applications by the fact that A * is hard to construct.

The above analysis has lead Kozlov [START_REF] Kozlov | Averaging of random operators[END_REF] to suggest a change of variable, the so called "harmonic coordinates". This has been further used by [START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF], [START_REF] Allaire | A multiscale finite element method for numerical homogenization[END_REF], [START_REF] Engquist | The heterogeneous multiscale method: A review[END_REF], [START_REF] Gloria | An analytical framework for the numerical homogenization of monotone elliptic operators and quasiconvex energies[END_REF] in various forms in order to propose a numerical integration method for approximating the construction of A * (giving rise to the methods that are known now as MsFEM, HMM..) .

In the following section, we restrict the generality of the frame and assume that there exists a clear separation between the only 2 scales in presence and that the situation is either periodic or a (stochastic) perturbation of a periodic situation.

Two-scale homogenization.

Let us now introduce the two-scale homogenization theory (without randomness as for now) to get explicit expressions for the homogenized problem ( [START_REF] Allaire | Homogenization and two-scale convergence[END_REF], [START_REF] Nguetseng | Homogenization of perforated domains beyond the periodic setting[END_REF]). We assume that the tensors A ǫ are functions of two coupled variables on the set locally defined by a fast microscopic variable ǫ -1 x coupled with the slow macroscopic variable x in D.

A ǫ (x) = A(x, x ǫ ), (2.10)
where

A(x, .) is Y = [0, 1] d -periodic, x ∈ D and A ∈ L ∞ (D, M α,β
). We denote by x the macroscopic variable and by y = x ǫ the microscopic variable. As A(x, x ǫ ) is locally periodic, one possible manner to get explicit expression for the homogenized problem is to perform a formal two-scale analysis with the following Ansatz:

u ǫ (x) = +∞ j=0 ǫ j u j (x, x ǫ ) = u 0 (x, x ǫ ) + ǫu 1 (x, x ǫ ) + ǫ 2 u 2 (x, x ǫ ) + ....... (2.11)
where, for any x ∈ D, the functions u j (x, .) are Y -periodic. Inserting the Ansatz (2.11) into equation (2.1) allows us to write u ǫ as:

u ǫ = u ⋆ + ǫ d i=1 w i ∂ i u ⋆ + r ǫ , (2.12)
where (w i (x, .)) 1≤i≤d are d Y -periodic, with zero average over Y , cell functions, parameterized by their macroscopic position x ∈ D that verify the following n cell problems (that now can be used as it is explicit)

-div y (A(x, y)[ e i + ∇ y w i (x, y)]) = 0, ∀y ∈ Y. (2.13)
The function u 0 = u ⋆ does not depend on the fast variable y = ǫ -1 x any more and verifies the homogenized problem:

   -div(A * (x) ∇u * (x)) = f, ∀x ∈ D, u * | Γ D = g D , A * (x) ∇u * . -→ n | Γ N = g N .
(2.14) the entries (A * ij (x)) 1≤i,j≤d of the homogenized matrix A * can now be explicitly computed with the cell functions w i (x, .),

A * ij (x) = Y A(x, y)( e i + ∇ y w i (x, y)) e j dy, (2.15)
where, provided u ⋆ ∈ W 2,∞ (D), the correction error r ǫ can be estimated to locally scale as ǫ (far enough from the boundary), and to globally scale as √ ǫ (see [START_REF] Allaire | Homogenization and two-scale convergence[END_REF]):

||r ǫ || H 1 Γ D (ω) ≤ C 1 ǫ||u * || W 2,∞ (ω) , ∀w c ⊂ D, ||r ǫ || H 1 Γ D (D) ≤ C 2 √ ǫ||u * || W 2,∞ (D) , (2.16)
with constants C 1 and C 2 depending only on D.

The local periodicity assumption (2.10) allows to completely determine the homogenized problem through explicit two-scale expressions. The derivation of the homogenized equation in the case of locally periodic coefficients serves as a basis in many numerical homogenization strategies. For the algorithm concerning the two-scale homogenization strategy, we refer to [START_REF] Engquist | The heterogeneous multiscale method: A review[END_REF].

For the numerical approximation of these problems, finite element methods can be proposed. The construction of the associated stiffness matrix of the macroscopic problem requires the construction of A * at every point in D or more precisely at every integration points on (say) each tetrahedra composing the macroscopic mesh; in what follows we shall denote this set as I N (D). Each of these evaluations involves the resolution of a micro scale problem (2.13). If a standard approximation method is used at the level of the microscopic problem, this may be a very expensive strategy, even though all of these problems are mutually independent and thus can be solved in a parallel way. This has motivated the use of reduced basis method to solve this series of macroscopic problems (see [START_REF]Reduced-Basis approach for homogenization beyond the periodic setting[END_REF]).

3. The reduced-basis method for non periodic homogenization problems.

Basics for the reduced-basis methods for 2-scale homogenization.

Let us first write the problem of interest in an abstract form so as to explain the basics of the Reduced Basis Method (RBM) more easily.

Let X , be the quotient space H 1 (Y )/IR homeomorphic to the space of all Yperiodic functions with zero average over Y , that belong to the Sobolev space H 1 (Y ).

X is equipped with the norm ||u||

X = ( Y ∇u.∇u) 1 2 .
For any integer i, 1 ≤ i ≤ d, the i-th cell problem (2.13) for the cell functions w i (x, .) rewrites in the following weak form: Find w i (x, .) ∈ X solution for:

a(w i (x, .), v; x) = f i (v; x), ∀v ∈ X . (3.1)
where x ∈ D appears to play the role of a parameter, a is a continuous and coercive bilinear form in X × X parameterized by x ∈ D defined as a(u, v; x) = Y A(x, y)∇u(y).∇v(y)dy, and f i are linear and continuous forms defined as

f i (v; x) = - Y A(x, y) e i .∇v(y)dy, 1 ≤ i ≤ d.
For many problems depending on some parameter (e.g. for an optimization problem, a control problem, an inverse problem or, like here for the resolution of the small scale problem) the set of all solutions forms a manifold of very small Kolmogorov dimension. In such cases, a reduced basis method can be implemented provided some further ingredients (albeit less fundamental see [START_REF] Chakir | Une mthode combine d'lments finis deux grilles/bases rduites pour l'approximation des solutions d'une E.D.P. paramtrique[END_REF]) are available. Among these ingredients the existence of a cheap error estimator is helpful. Then the RBM can be implemented into two steps through what is known as offline/ online approach. Let us explain it in our frame: Let us first introduce, for any i, the set of all solutions to problem (2.13) S i = {w i (x, .), x ∈ D} that is a manifold over D.

The first step -that may be expensive but done once and for all and, may be, before the RBM approximation is actually used -consists in choosing, e.g. through a greedy method, conducted with an available error estimator, the values of x n that are to be chosen to enrich at most the approximation quality of X iN for S i (see [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF] and algorithm 3.2 below). For these values, we then approximate the associated solutions w i (x n , .), n = 1, . . . , N with some classical (e.g. finite element method) accurate enough. After this is done, still in the offline stage, some scalar products are computed and bilinear forms working on w i (x n , .), n = 1, . . . , N , so that the evaluation of the RBM stiffness matrix will be rapidly achieved. Most of the times, this will involve the use of the Empirical Interpolation Method (EIM) (see [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF], [START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF]).

In a second step, called online stage, for every new values of the integration point x ∈ D, and any i, 1 ≤ i ≤ d , the Galerkin approximation w iN (x, .) ∈ X i,N = Span{w i (x n , .), n = 1, . . . , N } of w i (x, .) satisfies:

a(w iN (x, .), v; x) = f i (v; x), ∀v ∈ X i,N . (3.2)
The stiffness matrix (a(w iN (x n , .), w iN (x m , .); x)) n,m is then built online in a very few operations (since every thing is prepared for this in the offline stage) and the problem is solved in O(N 3 ) operations, providing a RBM approximation for w i (x, .), possibly with an associated error estimate. We refer to [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF] for a general presentation.

The RBM approach for the parameterized cell problems (2.13) should significantly decrease the expense of computations in terms of CPU time for the homogenization problems even if we take into account the construction time of the offline stage as we shall see latter.

3.2. Approximation of the limit coefficients and a posteriori error analysis. The macroscopic stiffness matrix is thus computed by using, at every point x involved in the macro scale numerical integration I N (D), the micro scale functions w iN (x, y). This is done by the following quantities:

A * ij,N (x) = Y A(x, y)( e i + ∇ y w iN (x, y)). e j dy = s ij,N (x) + Y A(x, y) e i . e j dy, (3.3)
where (s ij,N ) 1≤i,j≤d is defined by:

s ij,N (x) = -f j (w i,N (x, .); x) = Y A(x, y) e j , ∇w i,N (x, y)dy. (3.4)
Let us recall now the two following theorems which enable us to define the a posteriori bounds error of this matrix. Theorem 3.1. [START_REF]Reduced-Basis approach for homogenization beyond the periodic setting[END_REF] For 1 ≤ i ≤ d, we define the a posteriori error estimation ∆ N (w i (x, .)) by: ∆ N (w i (x, .)) = ( a(w i (x, .)w iN (x, .), .; x) X ′ )/α, and the corresponding effectivities η i (x) by

η i (x) = ∆ N (w i (x, .))/ w i (x, .) -w iN (x, .) X .
They satisfy the N-independent inequalities:

1 ≤ η i (x) ≤ β -1 α .
Theorem 3.2. [START_REF]Reduced-Basis approach for homogenization beyond the periodic setting[END_REF] Let w i (x, .) (resp. w j (x, .)) be a solution of the problem (3.1) i (resp. (3.1) j ) and w iN (x, .) (resp. w jN (x, .)) a solution of the problem (3.2) i (resp. (3.2) j ) with i = j, we have

| s ij,N (x) -s ij (x) |≤ ∆ s ij,N (x), with ∆ s ij,N (x)
= a(w i (x, .)w iN (x, .), .; x) X ′ a(w j (x, .)w jN (x, .), .; x) X ′ α .

Remark. The above definitions (e.g. ∆ N (w i (x, .))) involve a norm in X ′ of the residuals (a(w i (x, .)w iN (x, .), .; x) = f i (.; x)a(w iN (x, .), .; x)) that may be expensive, a priori, to evaluate but for which there also exists an offline/ online strategy that we shall not dwell upon and refer to [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF]. Such a strategy allows to perform these evaluations in O(N 2 ) operations.

To conclude on this matter, let us summarize the construction of the reduced basis Algorithm 1:

1. For some given x 1 , we compute the finite element approximation w i,h (x 1 , .) of w i (x 1 , .) using (3.1).

2. We set ℓ = 1, and ξ 1 = w i,h (x 1 , .) w i,h (x 1 , .) X , the first reduced function.

3. Then while ℓ < N :

(a) we compute for every x ∈ I N (D) the reduced basis approximations w i,ℓ (x, .) in X i,ℓ = span{ξ j , 1 ≤ j ≤ ℓ} for the cell problem (3.2).

(b) for

x ℓ+1 = argmax x∈I N (D)
∆ s ii,ℓ (x) w i,ℓ (x, .) X we compute the finite element approximation w i,h (x ℓ+1 , .) of w i (x ℓ+1 , .) using the (micro) finite element approximation of (3.1).

(c) we set ξ ℓ+1 (y) = R ℓ+1 (y) R ℓ+1 (y) X , the (ℓ + 1) th reduced basis, where R ℓ+1

is the remainder of the projection of w i,h (x ℓ+1 , .) on the ℓ-dimensional reduced basis:

R ℓ+1 (.) = w i,h (x ℓ+1 , .) - ℓ k=1 (w i,h (x ℓ+1 , .), ξ k ) X ξ k (.).
(d) we do ℓ = ℓ + 1 and go back to (a).

Numerical results.

In this subsection, we present some numerical results illustrating different points presented previously. We consider a case where d = 3, the domain D =]0; 1[ 3 and the permeability is isotropic and diagonal with coefficients on the diagonal equal to

a ii (x, y) = 4 + cos(2y 1 )x 1 + cos(2y 2 )x 2 + cos(2y 3 )x 3 .
We first check that the complexity of the manifold of all micro scale solutions when x varies is small. For this analysis, we perform a large number of finite element approximations of the cell problem for various values of x. The computations are performed with the software FreeFem++ [START_REF] Hecht | FreeFem++[END_REF], using IP 1 finite elements. The considered mesh contains 48.000 elements (tetrahedra). The following figures show the decreasing of the relative maximal error

A N = max 1≤i≤3,x∈I N (D)
∆ N (w i (x, .)) w i (x, .) X , the relative error of the cell functions

B N = max 1≤i≤3,x∈I N (D) w i (x, .) -w iN (x, .) X w i (x, .) X
and the relative error of s ij ,

C N = max 1≤i,j≤3,x∈I N (D) |s ij (x) -s ij,N (x)| |s ij (x)| ,
with respect to the number of basis functions.

As we can see in those figures (Figures 3.1 and 3.2), the considered relative errors decrease rapidly with respect to the number of basis and the one corresponding to s ij,N converge twice as fast as the others. The table (Fig. 3.3) shows the CPU time (in seconds) given by the Freefem++ code to approximate the FE matrix for the homogenized problem either with a direct FE approach or with a RB method for a porous medium, where NB is the bases number, N Y is the number of the degrees of freedom of the cell, N is the number of the degrees of freedom of D, offline part is the time of the bases construction, RBM is the time of the calculation of A * using the reduced basis method and F EM is the time of the calculation of A * using the finite element method.

We can notice that the RB method is much faster than the FEM which shows the efficiency of the RB method, even if we include the offline part, that corresponds to solving less than N 2 finite element solutions which is far less than twice the RBM work, even in the largest case in table 3.3 (N Y = 1, 000 , N = 9, 261).

About randomness.

In order to account for some randomness we choose an approach inspired by [START_REF] Blanc | Une variante de la theorie de l'homogénéisation stochastique des opérateurs elliptiques[END_REF] where the randomness is incorporated in the description of the microstructure by using random deformations of the periodic coefficients previously defined. We thus introduce a (random) parameter dependent function φ θ (y) = -2θy 3 + 3θy 2 + (1θ)y which is a simple one to one mapping from [0, 1] onto [0, 1] for every θ ∈] -2, 1[. Starting from a matrix Â(y), the definition of the previous coefficients A(x, y) in every cell is obtained, for every 1 ≤ i ≤ d, by replacing each component y i by φ θ (y i ) (we denote by Φ θ (y) = (φ θ (y 1 ), . . . , φ θ (y n ))). The coefficients are then A(x, y; θ) = Â(Φ θ(x,ω) (y)) (3.5) meaning that at any macroscopic point x, the associated cell-matrix is randomly modified through the mapping Φ θ(x,ω) .

For the simulation, this does not modify much the homogenization approach : At every point x ∈ I N (D) we compute the micro scale problem (2.13) with a RBM based on a discrete space X N that is constructed, as explained before, but where the parameter x appears now through θ.

Note that this is a model for taking into account some randomness, it is somehow similar to [START_REF] Blanc | Une variante de la theorie de l'homogénéisation stochastique des opérateurs elliptiques[END_REF] but differs from it in the sense that here we do not construct a global random map over IR d to modify the periodic cell structures. Since we are mainly interested in numerical simulations, this cell by cell modification is sufficient. It is out of the scope of this first paper to make the actual link between these two versions or the theory developed in [START_REF] Bourgeat | On the Stochastic Two-Scale Convergence in the Mean and Applications[END_REF][START_REF] Bourgeat | On the double porosity model of single phase flow in random media[END_REF] and further understand to what extent the discrete version presented here is an approximation of the continuous random stationary diffeomorphism introduced in [START_REF] Blanc | Une variante de la theorie de l'homogénéisation stochastique des opérateurs elliptiques[END_REF]. Note that, in opposition to these cited papers, our limit approximation remains random while the limit problem and solution of [START_REF] Blanc | Une variante de la theorie de l'homogénéisation stochastique des opérateurs elliptiques[END_REF] are deterministic (note however that this is also the case of the approximation proposed in [START_REF] Costaouec | Approximation numérique dune classe de problèmes en homogénéisation stochastique (Numerical approximation of a class of problems in stochastic homogenization)[END_REF] that remains also random). We plan to analyze the link between these two notions in a future paper.

The random macroscopic stiffness matrix is now computed by using, at every point x involved in the macro scale numerical integration I N (D), the micro scale functions w iN (x, Φ θ(x,ω) (y)). This is done by the following quantities:

A * ij,N (x, θ) = Y Â(Φ θ(x,ω) (y))( e i + ∇ y w iN (x, Φ θ(x,ω) (y))). e j dy = s ij,N (x, θ) + Y Â(Φ θ(x,ω) (y)) e i . e j dy. (3.6)
Remark. Note that, as is often the case when changes of variables are involved as a parameter for RB approximation the manifold of solutions is defined as

S i = { w i (x, .; θ) = w i (x, [Φ θ(x,ω) ] -1 (.)), x ∈ D} (3.7)
rather than the original

S i = {w i (x, .), x ∈ D};
since, most often, the Kolmogorov dimension is smaller in the former case.

Remark. Note that the matrix  could also depend on x, the random cell matrix would then be A(x, y; θ) = Â(x, Φ θ(x,ω) (y)); there would then be two parameters for the RB frame: x and θ.

Application.

The blood microcirculation in the capillaries can be modelled by a flow in a porous media which is governed by the Darcy law, established by Henri Darcy (1857) as a result of extensive work on the water flow in a layer of sand filter. Considering that the blood is a Newtonian incompressible fluid, we can build the corresponding problem modelling the blood flow in the capillaries in 2D and 3D, using the Darcy law in a domain decomposed in small cells having the same small size ǫ but different structures in each cell:

v ǫ = -A ǫ ∇p ǫ , divv ǫ = 0. (4.1)
Where v ǫ is the velocity, p ǫ the pressure and A ǫ a symmetric tensor representing the quotient between the permeability and the viscosity of the blood. We suppose that the viscosity has a constant value on the domain. There are many applications of the Darcy law in the biology, for instance Vankan and Van Donkelaar [START_REF] Vankan | Poroelasticity of saturedsolid with an application to bloodperfusion[END_REF] describe the blood circulation in the muscles by five layers (arteries, arterioles, capillaries, venules and veins) using Darcy law to represent horizontal flux in each layer, and Baeur D. & all [START_REF] Bauer | First phase microcirculatory reaction to mechanical skin irritation:a three layer model of a compliant vascular tree[END_REF][START_REF] Bauer | A three-layer continuous model of porous media to describe the first phase of skin irritation[END_REF] have described a vascular network by a three-layer model in order to model the effect of the irritation on the microcirculation. For a flow in a periodic porous medium, A ǫ can be represented by a periodic function A ǫ (x) = A( x ǫ ) ( [START_REF] Bakhvalov | Homogenization: averaging processes in periodic media[END_REF][START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF]). To model a porous media with a slight modification of structures from one cell to the other, the tensor A ǫ can be considered as a function of two variables A(x, x ǫ ) (see section 2.2).

4.1. Description of the model. In this section we aim at simulating the blood flow in a network of capillaries. The model we propose will be represented by a non uniform two scale system with cells of size ε where the part of the cell representing the capillaries is permeable and connects the adjacent cells, while the remaining part is not, and thus is represented by a permeability tensor which almost vanishes on these regions of the domain. Each cell will be defined as the (random) deformation of a reference cell: let ŷ = (ŷ 1 , . . . , ŷd ) and let Λ = [0, 1] d be the reference cell where the permeability is defined as follows:

∀ŷ ∈ [0, 1] d , K(ŷ) = K ij (ŷ) 0≤i,j ≤d , with K ij = 0 if i = j and K ii (ŷ) = 1 σ √ 2π e - 1 2
( ŷ1 -0.5 σ

) 2 + • • • + 1 σ √ 2π e - 1 2
( ŷd -0.5 σ

) 2 + ς,
where σ = 0.07 and 0 < ς << 1. We denote by X = H 1 ♯ ( Λ)/IR.

In order to build a random porous media, we will vary randomly the permeability as explained in subsection 3.4 :

K = K θ in each cell by setting K θ (y) = K • Φ θ (y) (remember that θ = θ(x, ω)) i.e. (K θ ) ii (y) = K ii • Φ θ (y) = 1 σ √ 2π e - 1 2
( φ θ (y 1 ) -0.5 σ

) 2 + • • • + e - 1 2
( φ θ (y d ) -0.5 σ and then, we will solve the corresponding cell problem in order to compute u ǫ solution of (2.1), by applying the reduced basis method: For all i ∈ {1, . . . , d}

) 2 + ς
-div(K θ (y)( e i + ∇ y w i (θ, y))) = 0 ∀y ∈ Y =]0, 1[ d , w i (θ, .) ∈ H 1 ♯ (Y )/IR. (4.2)
In this case, as shown in the figure 4.1, at each quadrature points of D, we choose a random value θ in a compact set K of ] -2, 1[. Then at each cell, the matrix K is different and leads to different size of pore.

From the last expression of K, we can deduce that the coercivity coefficient of the bilinear form a is very small, since it is of the order of α A ≃ ς. From theorem 3.1 we derive a non pertinent result because the effectivities η i is bounded by β -1 /α A . To circumvent this problem, we will take into account the fact that • the permeability constant almost vanishes in some part of the domain of computation • the cell design is inherited from a change of variable this allows to map every solution on the reference domain Λ and introduce, over Λ a new scalar product

[[ ŵ, v]] X = Λ( K ∇ ŵ ∇v)
for all ŵ, v ∈ X and we denote by |||.||| X the associated norm.

Lemma 4.1. Let us associate to every w and v ∈ X , the functions

ŵ = w • Φ -1 θ and v = v • Φ -1
θ , and let us assume that we choose a random value θ in a compact set K of ]-2, 1[, then the bilinear form a(w, v; θ) = Y K θ (y)∇w(y).∇v(y)dy is uniformly continuous and elliptic with respect to the previously defined scalar product in the sense that: there exists two constants C(K) and

C ′ (K) such that a(w, v; θ) ≤ C(K)||| ŵ||| X |||v||| X , a(w, w; θ) ≥ C ′ (K)||| ŵ||| X |||v||| X . (4.3)
Proof. We use the change of variable based on Φ θ and write

a(w, v; θ) = Y K θ ∇w.∇v = Y ( K • Φ θ ).∇w.∇v = Λ( K • Φ θ • Φ -1 θ ). ((∇w) • Φ -1 θ ).((∇v) • Φ -1 θ ).|J θ | = Λ K.((∇w) • Φ -1 θ ).((∇v) • Φ -1 θ )(φ -1 θ ) ′ (ŷ 1 ) . . . (φ -1 θ ) ′ (ŷ d ) = Λ K J θ .∇ Φ -1 (w • Φ -1 θ ).J θ .∇ Φ -1 (v • Φ -1 θ ).(φ -1 θ ) ′ (ŷ 1 ). . . . .(φ -1 θ ) ′ (ŷ d ), with J θ the Jacobian matrix of Φ -1 θ , |J θ | = (φ -1 θy ) ′ (ŷ 1 ). . . . .(φ -1 θy ) ′ (ŷ d ) , ∇ Φ -1 the gradient with respect to Φ -1 θ . We then proceed a(w, v; θ) = Λ KJ θ . ∇ ŵ.J θ . ∇v.(φ -1 θ ) ′ (ŷ 1 ) . . . (φ -1 θ ) ′ (ŷ d ) = Λ K11 d ŵ dŷ 1 . dv dŷ 1 .((φ -1 θ ) ′ (ŷ 1 )) -1 .(φ -1 θ ) ′ (ŷ 2 ). . . . .(φ -1 θ ) ′ (ŷ d ) + . . . + Λ Kdd d ŵ dŷ d dv dŷ d .(φ -1 θ ) ′ (ŷ 1 ). . . . .(φ -1 θ ) ′ (ŷ d-1 ).((φ -1 θ ) ′ (ŷ d )) -1 .
And the proof is complete by introducing the constants C(K) (resp. C ′ (K)) as being the maximum (resp. the minimum) value of the quatities

(φ -1 θ ) ′ (ŷ 1 ).(φ -1 θ ) ′ (ŷ 2 ). . . . .((φ -1 θ ) ′ (ŷ i )) -1 . . . . .(φ -1 θ ) ′ (ŷ d ) over all possible 1 ≤ i ≤ d and θ ∈ K: i.e. C ′ (K) = β d-1 A ζ -1 A and C(K) = β -1 A ζ d-1 A with ζ A = max y1∈[0,1] (φ -1 θy ) ′ (y 1 ) and β A = min y1∈[0,1]
(φ -1 θy ) ′ (y 1 ). In what follows, we shall denote by â( ŵ, v, θ) the last bilinear form defined in the reference cell Λ:

â( ŵ, v, θ) = Λ K11 d ŵ dŷ 1 . dv dŷ 1 .((φ -1 θ ) ′ (ŷ 1 )) -1 .(φ -1 θ ) ′ (ŷ 2 ). . . . .(φ -1 θ ) ′ (ŷ d ) + . . . + Λ Kdd d ŵ dŷ d dv dŷ d .(φ -1 θ ) ′ (ŷ 1 ). . . . .(φ -1 θ ) ′ (ŷ d-1 ).((φ -1 θ ) ′ (ŷ d )) -1 .

4.2.

A posteriori error estimator in the degenerate case. In this subsection, we will establish the a posteriori errors corresponding to the new inner product. Theorem 4.2. With this new inner product, we define the a posteriori estimate ∆i,N (θ) = ( f ia(w iN (θ, .), .; θ) X ′ )/C ′ (K) and the effectivities ηi,N (θ) = ∆i,N (θ, .)/ ŵi (θ, .) -ŵiN (θ, .))|| X .

We have the inequality

1 ≤ ηi,N (θ) ≤ C(K) C ′ (K) . (4.4)
Proof. First, we have:

f i (v; θ) -a(w iN (θ, y), v; θ) = a(w i (θ, .) -w iN (θ, .), v; θ) ∀v ∈ X . (4.5)
Then lemma 4.1 and the coercivity of the form a allows us to obtain this inequality:

C ′ (K)||| ŵ||| 2 X ≤ a(w, w; θ) = â( ŵ, ŵ; θ).
By replacing w by w i (θ, .)w iN (θ, .) we obtain:

C ′ (K)||| ŵi (θ, .) -ŵiN (θ, .)||| X ≤ |a(w i (θ, .) -w iN (θ, .), w i (θ, .) -w iN (θ, .); θ)| |||w i (θ, .) -w iN (θ, .)||| X ≤ sup v∈ X |f i (v; θ) -a(w iN (θ, y), v; θ)| |||v||| X ≤ |||f i (.; θ) -a(w iN (θ, y), .; θ)||| X ′ .
The continuity of the linear form a allows us to obtain: For all w ∈ X , v ∈ X , Proof. According to the definition of the matrix s and the linearity of the function f we have:

|s ij (θ) -s ij,N (θ)| = |f j (w i (θ, .) -w iN (θ, .), .)|
= |a(w j (θ, .), w i (θ, .)w iN (θ, .); θ)|.

(4.7)

The bilinear form a is symmetric as the matrix A(θ, y). We have:

|a(w j (θ, .), w i (θ, .)w iN (θ, .); θ)| = |a(w i (θ, .)w iN (θ, .), w j (θ, .); θ)|. (4.8)

On the other hand, w jN (θ, .) ∈ X N and thus |a(w i (θ, .)w iN (θ, .), w jN (θ, .); θ)| = 0, then |a(w i (θ, .)w iN (θ, .), w j (θ, .); θ)| = |a(w i (θ, .)w iN (θ, .), w j (θ, .)w jN (θ, .); θ)|. (4.9) Using (4.7), (4.8) and (4.9), we obtain: The inequalities (4.11) and (4.13) give the result.

|s ij (θ) -s ij,N (θ)| = |a(w i (θ, .) -w iN (θ, .
5. Numerical Results. In this section, we show numerical results in D =]0, 1[ d for d = 2, 3, for θ ∈ K = [-0.4, 0.4] and f = 0. These numerical computations are performed with the Freefem++ software [START_REF] Hecht | FreeFem++[END_REF] with classical P1 Lagrange finite elements.

For the implementation of the RBM, we want to note now that if the permeability K θ (y) would have been affine as a function of θ, then a(w, v, θ) would also be affine in θ, resulting in:

a(w iN , v; θ) = k t k (θ)a k (w iN , v), (5.1)
where t k depends on θ while a k only depends on w iN and v. Then an extremely efficient offline-online computational strategy may be developed (see [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF]). In the offline stage we could have form a k (ξ j , ξ i ), 1 ≤ i, j ≤ N ; in the online stage we would need only to assemble and invert a(ξ j , ξ i ; θ) = k t k (θ).a k (ξ j , ξ i ) resulting in computation that would be independent of N the number of degrees of freedom of the finite element solver. Since N ≤ N large computational savings can be achieved. Unfortunately, in our case K θ (y) is not affine in θ, so we use the empirical interpolation method [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF] that recovers the previous online N -independence even in the presence of non-affine parameter dependence. We do not elaborate on this matter since it is quite classical by now.

First, we show the variation in logarithmic scale of: We conclude from the figures (Fig. 5.2) that, the errors of the solution A N computed with the original product scalar depend on the values of ς (note that the factor 4 between ς = 1e -10 and ς = 1e -14 is exactly visible on 5.2). The figure (Fig. 5.3) shows that the errors with the new inner product values, A N is not affected by the variation of ς. Moreover, the figures (Fig. 5.2) and (Fig. 5.3) illustrate clearly that the errors of the solutions computed with the new product scalar is much better than the other one. Hence the advantage of this modified approach. to approximate the FE matrix for the homogenized problem either with a direct FE approach or with a RB method for a porous medium and for a small ǫ.

A N = max
As we can see, the RB method is again much faster than the FE method.

In the figures 5.6-5.9 and just for illustrations, we show color results computed on a coarse mesh (not so fine then that used if figure 5.4), for the problem with different boundary conditions to illustrate different models of applications and to show the powerful of this method. The following figures show some applications when the boundary conditions at the input and the output are not aligned. and the vein. In this case, we take the same boundary conditions as in 2D. As we see in figures (Fig. 5.10) and (Fig. 5.11), in this case (3D), we have the same conclusions as in the two dimensional case about the new product scalar which gives results independent of ς. This color figure (Fig. 5.12) represents the flow repartition in a porous media in three dimensional As in 2D, the table (Fig. 5.13) shows the CPU time (in seconds) given by the Freefem++ code in 3D with the same notations. As we can see, we have the same remarks concerning the RB method which is faster than the FE method.

Conclusion:

In this paper, we use the reduced basis method to solve and speed up the resolution of the Darcy equation in the porous media using the locally homogenization model. This Darcy problem is used as a model for the blood microcirculation in the capillaries. The model is enriched with random properties leading to a random macroscopic stiffness matrix. The microscopic behavior of the micro vessels leads to the introduction of a mathematical model with very small values of the viscosity in some regions. This yields to a non-pertinent a posteriori estimates making the construction of the reduced basis problematic. We propose a new microscopic scalar product to circumvent this difficulty and obtain pertinent results in 2D and 3D. In a future paper we shall analyze more in depth the stochastic behavior (averages, variances . . . ) of the solutions to this model.

Fig. 3 . 1 .Fig. 3 . 2 .Fig. 3 . 3 .

 313233 Fig. 3.1. A N in logarithmic scale with respect of N

Fig. 4 . 1 .

 41 Fig. 4.1. The domain D with variable x on the left with example of cell structures on the right with variable y.

aTheorem 4 . 3 .

 43 (w, v; θ) ≤ C(K)||| ŵ||| X |||v||| X , hence |f i (.; θ)a(w iN (θ, y), .; θ)| = |a(w i (θ, .)w iN (θ, .), v; θ)| ≤ C(K)||| ŵi (θ, .) -ŵiN (θ, .)||| X |||v||| X which leads to:|||f i (.; θ)a(w iN (θ, y), .; θ)||| X ′ ≤ C(K)||| ŵi (θ, .) -ŵiN (θ, .)||| X (4.6)and ends the proof of (4.4).This theorem states that the effectivities ηi,N (θ) are bounded by a term independent of ς. The next theorem gives a bound on s ij,N . With this new inner product, the theorem 3.2 becomes| s ij (θ)s ij,N (θ) |≤ ∆ij,N (θ)with: ∆ij,N (θ) = ||a(w i (θ, .)w iN (θ, .), .; θ)|| X ′ ||a(w j (θ, .)w jN (θ, .), .; θ)|| X ′ C ′ (K) .

  ), w j (θ, .)w jN (θ, .); θ)|. (4.10)Then |a(w i (θ, .)w iN (θ, .), w j (θ, .)w jN (θ, .); θ)| ≤ |||a(w i (θ, .)w iN (θ, .), .; θ)||| X ′ ||| ŵj (θ, .) -ŵjN (θ, .)||| X(4.11)and |||a(w j (θ, .)w jN (θ, .), .; θ)||| X ′ ≥ |a(w j (θ, .)w jN (θ, .), w j (θ, .)w jN (θ, .); θ)| ||| ŵj (θ, .) -ŵjN (θ, .)||| X ≥ C ′ (K)||| ŵj (θ, .) -ŵjN (θ, .)||| X . (4.12) which leads to||| ŵj (θ, .) -ŵjN (θ, .)||| X ≤ 1 C ′ (K)|||a(w j (θ, .)w jN (θ, .), .; θ)||| X ′ . (4.13)

5 . 1 .

 51 1≤i≤d,θ∈K ∆N (w i (θ, .)) || ŵi (θ, .)|| X with respect of N , when ς = 10 -14 and 10 -10 . Second, we compare this error values in logarithmic scale, of the a posteriori estimate of the solution with the new scalar product A N and the a posteriori estimate of the solution with the original inner product A N = max 1≤i≤d,θ∈K ∆ N (w i (θ, .)) ||w i (θ, .)|| X . Results for d = 2. The repartition of the permeability is given by the next figure 5.1. We have: Γ D = Γ 1 = {0}×]0.35, 0.65[, with the corresponding data g D = 0, and Γ N = Γ 2 = {1}×]0.35, 0.65[, with he corresponding data g N = -1, Γ 3 is the rest of the boundary, with he corresponding data g N = 0, hence the flow is going from the right to the left.

Fig. 5 . 1 .

 51 Fig. 5.1. The domain D representing the capillary bed, a model of artery on the right and model of vein on the left.

Fig. 5 . 2 .Fig. 5 . 3 .

 5253 Fig. 5.2. A N with respect of N , with ς = 1e -10 (left) and ς = 1e -14 (right)

Fig. 5 . 4 .Fig. 5 . 5 .

 5455 Fig. 5.4. Modelisation of the velocity

Fig. 5

 5 Fig. 5.6.

Fig. 5 .Fig. 5 .

 55 Fig. 5.10. A N with respect of N , with ς = 1e -10 (left) and ς = 1e -14 (right)

Fig. 5 . 12 .Fig. 5 . 13 .

 512513 Fig. 5.12. Modelisation of the velocity in 3DNB N Y N RBM FEM 5 125 216 0.13s 41 s 5 1,000 216 0.13 s 326 s 5 125 9,261 8 s 2,500 s 5 1,000 9,261 8 s 20,000 sFig. 5.13. Comparison between the finite element method and the reduced basis method in 3D.


J Barr

COMBINATORIAL MODELS OF RIGIDITY AND RENORMALIZATION

We first introduce the percolation problems associated with the graph theoretical concepts of (k, l)-sparsity, and make contact with the physical concepts of ordinary and rigidity percolation. We then devise a renormalization transformation for (k, l)-percolation problems, and investigate its domain of validity. In particular, we show that it allows an exact solution of (k, l)-percolation problems on hierarchical graphs, for k ≤ l < 2k. We introduce and solve by renormalization such a model, which has the interesting feature of showing both ordinary percolation and rigidity percolation phase transitions, depending on the values of the parameters.

Introduction

Consider an ensemble of bars connected together at joints, around which they can freely rotate, embedded in a two dimensional space. The ensemble of bars and joints is naturally associated with a graph: joints are the vertices of the graph, and bars its edges. The question of whether or not this ensemble constitutes a rigid body is a fascinating one. It found a convenient answer with Laman's theorem [1], which characterizes rigidity in a purely graph theoretical way 1 . It implies for instance that a graph with N vertices must have at least M = 2N -3 edges to be rigid. Here the number 2 corresponds to the number of degrees of freedom of a joint and the number 3 to the number of degrees of freedom of a rigid body, in two dimensions. The concepts of (k, l)sparsity and (k, l)-tightness, introduced a long time ago in graph theory [2], allow to generalize the problem [3]. For instance, bar-joint rigidity in two dimensions corresponds to (k, l) = (2,3). Some values of (k, l) have a clear physical meaning in terms of other models of rigidity, others do not, see Sec. 2.

1 Throughout this paper, we will be concerned only with generic rigidity, which can be reduced to a combinatorial problem. See for instance [21] for a discussion

1
The percolation problem associated with bar-joint rigidity has drawn the attention of physicist's for a long time [4]. It can be stated as follows: take a large lattice, where each site is a joint and each bond a bar; consider that each bar is effectively present with probability p. When p is increased, the system goes from a "floppy" phase where only small rigid clusters exist, to a "rigid" phase, where there is one macroscopic rigid cluster. In between lies the rigidity percolation transition.

To study this problem, Jacobs, Thorpe, Duxbury and Moukarzel have introduced the so called "pebble game" algorithm, which implements in an efficient way the combinatorial characterization of Laman [START_REF] Jacobs | [END_REF]8], see also [9]. The interpretation of the numerical simulations on regular two dimensionnal lattices have fueled some debate, but the following picture seems favored: the rigidity percolation transition in such 2D bar-joints lattices is second order, in a universality class different from ordinary connectivity percolation [START_REF] Jacobs | [END_REF]13]. We note that rigidity percolation may also be studied on trees or random graphs, in which case it is exactly soluble, and displays a first order transition, at variance with ordinary percolation [11]. To our knowledge, there exists no analytical theory, even approximate, for the critical exponents of 2D bar-joint rigidity percolation on regular lattices.

Except for some special cases, such as ordinary percolation and the bar-joint rigidity percolation described above, percolation problems associated to (k, l) sparsity seem to have received relatively little attention in the physics literature (see Sec.2 for a more detailed discussion). In particular, little is known about the associated percolation transitions and their universality classes. A natural idea to get insight into this problem is renormalization theory. The purpose of this paper is to introduce a renormalization tool to study (k, l)-percolation, generalizing the method of [16]. In Sec. 3, we will determine the conditions under which the renormalization procedure is useful. In particular, we will show that it allows a complete solution for some (k, l)-rigidity models on hierarchical graphs [18,16]. In Sec.4, we present and solve explicitly a model of (3,3) percolation (which has a physical interpretation as bodies and bars percolation in two dimensions). We will show that the renormalization transformation on this model has two trivial fixed points, corresponding to the full and empty graph respectively, and two critical fixed points. One fixed point, corresponding to ordinary percolation has three unstable direction, whereas the other one has one unstable direction, and thus governs the critical behavior for almost all values of the parameters. This model provides an example of a system showing both ordinary and rigidity percolation transistions.

Rigidity and graph theory

We first briefly recall the graph theoretic description of bar-joints rigidity [1]. The basic idea is constraint counting: each joint has two translational degrees of freedom, and each bar, by fixing the distance between two joints, removes one degree of freedom. Since a rigid body in two dimensions has three degrees of freedom (two translations and one rotation), an ensemble with N joints needs at least M = 2N -3 bonds to be rigid. In addition, these bars must be "well distributed" among the joints: if there are exactly 2N -3 bonds, there should be no subensemble with n joints containing strictly more than 2n -3 bonds.

Clearly, the numbers 2 and 3 in this description are adapted to barjoints rigidity theory in two dimensions, but it is natural to generalize the idea to other pairs of numbers (k, l). [20,3]. It is also natural to allow for multiple edges in a graph: such graphs are called multigraphs. Although, to our knowledge, they have not appeared in this form in the physics literature, the following definitions about (k, l)-rigidity are not new [3]; we recall them for self-consistency of the article. Throughout the article, we shall always consider multigraphs, although we shall keep calling them graphs.

Given a graph G = (V, E), where V represents the set of vertices of G and E the set of its edges, we say that

G = (V , E ) is a subgraph of G if V ⊂ V , E ⊂ E,
-(k, l)-sparse if every subgraph G of G with n > 1 vertices contains at most m = kn -l edges; -(k, l)-tight if in addition M = kN -l; -(k, l)-spanning if there is a (k, l)-tight subgraph spanning all vertices of G; -(k, l)-redundant if it is not (k, l)-sparse.
Two vertices A and B are said to be (k, l)-rigidly connected if there is a (k, l)-tight subgraph containing A and B.

When there is no ambiguity, we shall sometimes omit the (k, l) indices, and talk about sparse, tight, and redundant graphs.

Notice that if l ≥ 2k, these definitions make no sense: there is no non trivial (k, l)-sparse graph in this case, since the simplest graph with 2 vertices and 1 edge is not (k, l)-sparse. More generally, the edges in a (k, l)-sparse multigraph have multiplicity at most 2k -l.

These graph-theoretic definitions may be immediately translated in the rigidity vocabulary. A (k, l)-tight graph is a (k, l) isostatically rigid, or minimally rigid, graph: it has just the right number of constraints to be rigid, and removing any constraint destroys rigidity. A (k, l)spanning graph is a rigid graph according to the (k, l) rule.

Redundant constraints play an important role in bar and joint rigidity theory. It is straightforward to generalize the concept to (k, l)rigidity. A graph possesses one or more redundant edge if and only if it is not sparse, ie if it contains a subgraph with n vertices, m edges, and m > kn -l. To define and count the number of (k, l)-redundant edges in a graph G = (V, E), one may proceed as follows:

(1) Start with the empty graph (V, ∅), without edges (2) Add the edges of E one at a time (3) Check if this addition creates a redundant subgraph, with n vertices, m edges, and m = kn -l + 1. (4) If this is the case, then reject the newly added edge and add one to the count of redundant constraints (5) Continue until all edges in E have been added. This procedure actually follows the steps of the algorithms used to study (2, 3)-rigidity [9,[START_REF] Jacobs | [END_REF]8,10], and more generally (k, l)-rigidity [3]. A priori, the number of redundant edges computed in this way could depend on the order chosen to add the edges. It is not the case, so that the number of redundant edges is a well-defined concept. This is also the minimum number of edges that one has to withdraw to make a redundant graph sparse. In general, there is a freedom of choice in the edges to be removed to turn a redundant graph into a sparse one, but their number is fixed.

Beyond bar-joint 2D rigidity, which is a physical realization of (2, 3)rigidity, (k, l)-rigidity encompasses several known situations:

• (1, 1)-rigidity corresponds to ordinary percolation, in any dimension. • (3, 3)-rigidity corresponds to the rigidity of systems composed of rigid bodies and bars in 2 dimensions. • More generally, (k, k)-rigidity is related to rigidity of bodies and bars systems in higher dimension [20]. • (3, 6)-rigidity is related to the rigidity of bond-bending networks in 3 dimensions [21,22]. In this case l = 6 = 2k, which is outside the usual range; the definition of sparsity needs to be slightly adapted. • (2, 0)-rigidity corresponds to the graph theoretic concept of 2orientability, which has been used to study bar-joint rigidity of random networks [23].

Clearly, a percolation problem may be associated to each case, by monitoring the largest (k, l)-rigid cluster in a graph. Understanding the general properties of these percolation processes, such as the order of the transition, the universality class, etc... is a largely open question, although several cases have been studied in details: among those, ordinary percolation is of course the best known by far. Percolation for bond-bending networks has been studied, because of its relevance in modeling network glasses and proteins [12], but we are not aware of studies on critical exponents. In [13], C. Moukarzel studies percolation problems on trees and random graphs for various numbers k (g in the notations of [13]) 2 . It is also known that when some bars and joints form a rigid component, they behave in every aspect like a rigid body. This idea has actually been numerically exploited [14]. Then one may expect that body-bar ((k, l) = (3, 3)) and bar-joint ((k, l) = (2, 3)) rigidity percolation belong to the same universality class [15]. Summarizing these informations, one may identify several conjectures and questions:

(1) (1, 1) percolation (that is ordinary percolation) and (2, 3) percolation on 2D regular lattices belong to two different universality classes: this conjecture seems well documented numerically. (2) One may also conjecture that (2, 3)-percolation (bar-joint) and

(3, 3)-percolation (body-bar) on 2D regular lattices belong to the same universality class. (3) On random graphs, it is known that the threshold for (2, 3)percolation is identical to the threshold for (2, 0)-percolation [23].

One may then wonder if, on random graphs, the threshold for (k, l)-percolation could be independent of l, for 0 ≤ l ≤ 2k -1.

A more general question is: given for instance a regular lattice in d dimensions; what does determine the order of the transition and the universality class of (k, l) percolation? Efficient algorithms, generalizing to (k, l)-sparsity the original pebble-game devoted to (2, 3) sparsity, have been recently introduced [3], so that a numerical investigation of this question is in principle possible. To our knowledge, it has not been undertaken. A natural analytical tool to progress in the studies of these percolation processes is the renormalization group. In section 3 we introduce a renormalization transformation, and we prove that in some sense it is equivalent to study the rigidity properties of the original and renormalized graphs, when the parameter range is restricted to k ≤ l. This allows to solve exactly models on hierarchical graphs for k ≤ l; in Sec. 4, we solve such a model for k = l = 3.

Renormalization rule for a graph

The renormalization rule we are about to describe is a generalization of the rule introduced on [16] for (k, l) = (2, 3). In this case, this rule is intuitive enough so that it does not really require a justification relying on graph theory, see [16] (this rule has been rediscovered independently in [17]). Things are not so obvious for general (k, l), and in particular we have to impose the restriction k ≤ l. In order to be sure to make correct statements, we feel it is useful to give precise graph theoretical definitions, and to proceed rigorously, not relying only on intuition. This certainly makes the article heavier to read, but is necessary to identify for instance the restriction k ≤ l, which is not completely intuitive.

From now on, we assume k ≤ l < 2k. We will denote |F | the cardinal of an ensemble F .

3.1.

Renormalization of a whole graph. Consider a multigraph H = (V, E), where V is the set of vertices and E the set of edges. We singularize two sites A, B ∈ V . We would like to replace the graph H by an "equivalent" (in a certain sense) graph with only two vertices A and B, and a certain number of edges M AB linking A and B. We define this renormalization step as follows: Let us call S AB the ensemble of all (k, l)-sparse subgraphs of H containing A and B. For h = (V h , E h ) in S AB , the number of its vertices and edges is respectively Ordinary percolation (k, l) = (1, 1) corresponds to the ordinary percolation case. Then 2kl = 1, and the number of renormalized edges between A and B is 0 or 1. Sparse graphs in this case are graphs without any loop. One sees that a loopless subgraph h

|V h | and |E h |. We define r AB (H) = max h∈S AB [|E h | -k|V h |] Note that since h is (k, l)-sparse, r AB (H) ≤ -l. Then M AB is given by (3.1) M AB = max (r AB (H) + 2k, 0) We have 0 ≤ M AB ≤ 2k -l.
= (V h , E h ) with |E h | ≥ |V h | -1 is necessary B Figure 1.
Example of renormalization of a simple graph with 4 vertices, k = l = 1 (ordinary percolation). The upper graph is renormalized into one edge; the lower graph is renormalized into the empty graph. Note that both graphs contain one redundant edge. for "percolation" across the graph. On the elementary "Wheatstone bridge" graph of Fig. 1, it is easy to see that this coincides with the Bar-joint rigidity percolation In this case (k, l) = (2, 3). Again 2k -l = 1, and the number of renormalized edges between A and B is 0 or 1. The rule says in this case: "Renormalize the graph by a edge if and only if A and B are rigidly connected (ie belong to the same rigid cluster)". This is intuitive, and this is the rule which has been used in [16]. Examples are given on Fig. 2.

Body-bar rigidity percolation In this case (k, l) = (3, 3). 2k -l = 3, so that the number of renormalized edges between A and B is 0, 1, 2 or 3. A and B have to be seen as rigid bodies, with three degrees of freedom each, and three edges are necessary for a rigid connection between them. Examples are given on Fig. 3.

The following property is useful to effectively compute renormalized graphs.

Proposition 3.1. Let H = (V, E) be sparse, with A, B ∈ V . Let e be an edge such that H = (V, E ∪ {e}) is not sparse. Then r AB (H) = r AB (H ).
As a consequence of this proposition, the renormalizations of H and H contain the same number of edges. In other words, adding redundant edges to a graph H does not modify its renormalization. We will use this as a tool to compute renormalized graphs in Sec. 4.

Proof: Clearly, r AB (H) ≤ r AB (H ), since an element of S AB (H) is also in a natural way an element of S AB (H ).

We call u and v the two vertices connected by the edge e; u, v ∈ V .

Let us now consider g = (V , E ) an element of S AB (H ), such that the maximum is attained in the definition of r AB (H ). Thus |E | -k|V | = r AB (H ). We want to construct a subgraph g = (V , E ) in S AB (H) such that |E | -k|V | ≥ r AB (H ). This will prove r AB (H) ≥ r AB (H ). We may assume that e ∈ E , otherwise it is enough to take g = g . Thus the vertices u and v are in V .

Let G 0 = (V 0 , E 0 ∪ {e}) be a minimal redundant subgraph of H . V 0 contains u and v the vertices connected by e, and G 0 = (V 0 , E 0 ) is a tight subgraph of H. We define V = V ∪ V 0 and g = (V , E ) as the subgraph of H induced by the vertices V ∪ V 0 . g is an element of S AB (H). g does not contain e, but we will show that

|E | -k|V | ≥ r(H ).
Assume first V 0 ⊂ V . In this case V = V . Then there is one edge e 0 in E 0 ∪ {e} that is not in E , otherwise g would not be sparse. e 0 is different from e, thus e 0 is also in

E . We conclude |E | -k|V | ≥ |E | -k|V | = r AB (H ).
We assume now that V 0 \V is not empty (we use the notation A\B to denote the elements in set A that are not in set B). We have the relations:

V = V (V 0 \V ) , E ⊃ (E \{e}) (E 0 \E )
where is a disjoint union. Indeed, e ∈ E and e / ∈ E . Thus

(3.2) |V | = |V | + |V 0 \V )| , |E | ≥ |E | -1 + |E 0 \E | Let G 1 = (V 0 ∩ V , E 0 ∩ E )
the subgraph defined by the intersection of V and V 0 . G 1 is a subgraph of H, so it is sparse; it is also a proper subgraph of G 0 , and it contains vertices u and v, so it is not tight, otherwise it would contradict the minimality of G 0 : indeed, (V 0 ∩ V , E 0 ∩ E ∪ {e}) would be a redundant subgraph of H , smaller than G 0 . We conclude that

(3.3) |E 0 ∩ E | < k|V 0 ∩ V | -l We use now V 0 = (V 0 ∩ V ) (V 0 \V ) and E 0 = (E 0 ∩ E ) (E 0 \E ).
Since G 0 is tight, we have Putting this together with (3.2), we have finally

(3.4) |E 0 ∩ E | + |E 0 \E | -k|V 0 ∩ V | -k|V 0 \V | -l = 0
|E |-k|V | ≥ |E |-k|V |-1+|E 0 \E |-k|V 0 \V | ≥ |E |-k|V | = r AB (H ) 3.2.
Renormalization of a subgraph. We have seen how to replace a graph by a certain number of renormalized edges. Now, we will see why and when it is licit to perform this renormalization step on a subgraph of a bigger graph. Consider a graph G = (V, E), V H ⊂ V and H = (V H , E H ) the subgraph of G induced by the vertices V H . We will say that H satisfies the renormalizability condition if there exist two vertices A and B such that all edges in E linking H to the rest of the big graph G are connected to A or B; see Fig. 4. In the following, we always assume this condition is satisfied.

For any subgraph of g ⊂ G, we can define the renormalized graph R(g), as follows: the part of g which is a subgraph of H is renormalized according to the rule detailed in 3.1; the remaining part of g is left unchanged. Note that if A or B, or both, do not belong to g, then the renormalized part of g contains no edge.

Let us call G = R(G) = ( Ṽ , Ẽ) the graph G where the subgraph H has been renormalized. We would like to show that it is in some sense equivalent to study (k, l)-percolation on G and G. For this purpose, we will prove the following: Proposition 3.2. In the setting described above: i) A and B are (k, l)-rigidly connected in G if and only if they are (k, l)-rigidly connected in G.

ii) The number of redundant edges in G and G are related through the formula n red (G) = n red ( G) + n red (H) where we have written n red (F ) for the number of redundant edges in a graph F . This last formula allows then to replace the problem of computing the number of redundant edges in G by the the problem of computing the number of redundant edges in the renormalized graph R(G) = G. Before proving 3.2 i) and ii), we state and prove a few lemmas.

Lemma 3.3. Let G = (V, E) and R(G) = ( Ṽ , Ẽ) be as above. Let g = (ṽ, ẽ) be a subgraph of R(G), and let us call h its renormalized part; h contains m edges. Then there exists g subgraph of G, such that i) g = R(g) ii) If A or B do not belong to ṽ, then, calling h = (V h , E h ) the intersection of g with H, E h is empty (ie h has no edge). iii) If A and B belong to ṽ, then h = (V h , E h ) is sparse and satisfies

m = |E h | -k|V h | + 2k
Proof: Note first that h has at most two vertices, A and B. If A or B do not belong to g, or if m = 0 then we choose for the vertices of h the same vertices as h, and put no edge.

If m > 0, we take h = (V h , E h ) an element of S AB (H) such that |E h | -k|V h | + 2k = m, and h minimal in the sense that no proper subgraph of h has the same properties. Then the renormalization of h produces exactly m edges between A and B. Indeed, suppose it produces m > m such edges; then there exists h = (V h , E h ) subgraph of h, in S AB (H), such that |E h | -k|V h | + 2k = m . By removing edges, we would construct a proper subgraph of h with the same properties as h.

To build g, we complete h in both cases by the part of g not concerned by the renormalization.

Lemma 3.4. Let g = (V g , E g ) be a (k, l)-redundant subgraph of G, minimal (in the sense that no proper subgraph of g is (k, l)-redundant). Then i) or ii) is true: i) g is a subset of H (that is g is included in the part of G to be renormalized). ii) g = R(g) is a (k, l)-redundant subgraph of G = R(G) (that is the image of g by renormalization is still (k, l)-redundant).
Proof: Suppose i) is not true, and A, B ∈ g. Let us call h = (V h , E h ) the intersection of g and H, and h = R(h). We call m the number of É edges of h. The set of vertices and edges of g are called V g and E g respectively. Then

|E g | = |E g| -m + |E h | and |V g | = |V g| + |V h | -2. Thus (3.5) |E g | -k|V g | + l = |E g| -k|V g| + l + |E h | -k|V h | + 2k
Now, h is (k,l)-sparse since g is a minimal redundant graph. Thus h ∈ S AB (h), and

(3.6) m ≥ |E h | -k|V h | + 2k
Putting together Eqs. (3.5) and (3.6), an using

|E g | -k|V g | + l > 0 we get (3.7) |E g| -k|V g| + l > 0
This proves that ii) is true. If i) is not true and neither A nor B are in g, then h is actually empty: otherwise, g would be a minimal redundant subgraph with two disconnected components, which is impossible. Then the renormalization does not modify g and ii) is true.

If i) is not true and A or B, but not both, are in g, then either h has no edge, and we are done since g is not modified by renormalization, or g is a minimal redundant graph made of two components that share only one vertex. It is easy to show by enumerating vertices and edges that this is impossible as soon as l ≥ k. This ends the proof.

Remark: Notice that for 0 ≤ l < k, this lemma is not true. Fig. 5 shows a counter example for (k, l) = (2, 1). It is not clear to us how it is possible to define a useful renormalization transformation in this case.

Lemma 3.5. Let g = (ṽ, ẽ) be a (k, l)-redundant subgraph of G. Then there exists g, (k, l)-redundant subgraph of G such that g = R(g) (that is: g is the image by renormalization of a (k, l)-redundant subgraph of G).

Proof: If ẽ does not contain any A -B edge, then it is enough to take g = (ṽ, ẽ).

We assume now that ẽ contain m > 0 edges linking A and B. This implies that M AB (H) > 0 (since M AB (H) is the number of A-B edges in G).

Using the same reasoning as in Lemma 3.3, we construct h = (V h , E h ) an element of S AB (H) such that |E h | -k|V h | + 2k = m and such that the renormalization of h produces exactly m edges between A and B. Define now g = (V g , E g ), with V g = ṽ ∪ V h and E g contains E h and ẽ minus the m edges between A and B. Then R(g) = g, and

|E g | -k|V g | + l = |ẽ| -m + |E h | -k|ṽ| + 2k -k|V h | + l = |ẽ| -k|ṽ| + l > 0
Then g is redundant.

Proof of 3.2 i):

Rigid connection in G ⇒ rigid connection in G: Suppose A and B are (k, l)-rigidly connected in G. Then there exists

G 1 = (V 1 , E 1 ) a (k, l)- tight subgraph of G containing A and B. Let us call H 1 = (V H 1 , E H 1 )
the intersection of G 1 and the subgraph to be renormalized H; G1 = R(G 1 ) = ( Ṽ1 , Ẽ1 ), H1 = R(H 1 ) = ( ṼH 1 , ẼH 1 ) the subgraphs of G images of G 1 , H 1 by renormalization; H1 has only two vertices A and B, and a certain number of edges between them.

We have

|V 1 | = | Ṽ1 | + |V H 1 | -2 ; |E 1 | = | Ẽ1 | + |E H 1 | -| ẼH 1 | H 1 is a subgraph of G 1 which is tight, so that H 1 is sparse. The subgraph of G 1 defined by the vertices V 1 \V H 1 ∪ {A, B} is also sparse.
This yields the two inequalities

k|V H 1 | -2k ≤ |E H 1 | ≤ k|V H 1 | -l
Then, by renormalization of H 1 , we have

| ẼH 1 | = |E H 1 | -k|V H 1 | + 2k
From this and the fact that G 1 is tight, we obtain finally | Ṽ1 | = k| Ẽ1 |-l. Now we want to prove that G1 is sparse. Let g = (ṽ, ẽ) be a subgraph of G1 . We have to show that |ṽ| -k|ẽ| + l ≤ 0. By Lemma 3.3, we construct g = (v, e). h = (V h , E h ) is the intersection of g with H 1 , and m is the number of edges in the renormalization of h. Then using point iii) of lemma 3.3 and

|v| = |ṽ| + |V h | -2 ; |e| = |ẽ| -m + |E h | , we conclude |ẽ| -k|ṽ| + l ≤ 0
This proves that G1 is (k, l)-tight, so that A and B are rigidly connected in G.

Rigid connection in G ⇒ rigid connection in G:

We assume now that A and B are rigidly connected in G, and want to show this is true also in G. Let G1 = ( Ṽ1 , Ẽ1 ) be a (k, l)-tight subgraph of G containing A and B. Using Lemma 3.3, we construct G 1 = (V 1 , E 1 ) subgraph of G, and call H 1 = (V H 1 , E H 1 ) the intersection of G 1 with H. First, reasoning as above, we have

|E 1 | = k|V 1 | -l.
We now have to prove that G 1 does not contain any (k, l)-redundant subgraph. Suppose it is not the case, and call g 1 a minimal redundant subgraph of G 1 . Then we have the alternative of Lemma 3.4. i) cannot be true because g 1 would be a subgraph of H 1 , which is sparse by Lemma 3.3. ii) cannot be true either, because we would have constructed a redundant subgraph of G1 , which is tight by hypothesis. We conclude that G 1 is tight, and that A and B are rigidly connected in G 1 .

Proof of 3.2 ii):

We want to prove the formula

(3.8) n red (G) = n red ( G) + n red (H)
To prove this, we follow the algorithm defining n red (G), as described in Sec. 2. Starting from the vertices V and no edge, we add the edges in E one by one, starting with the edges which are in H, the subgraph to be renormalized. If the newly added edge is redundant, we discard it and add one to the count of n red (G). This way, we construct a sequence of sparse graphs

G 0 = (V, ∅), G 1 , . . . , G M H , . . . , G M = (V, E M ) ,
where M H = |E H | and M = |E|. We have n red (G) = M -|E M |. E M contains the edges in E which have not been found redundant in the sequential edge addition process. Consider the associated sequence of renormalized graphs

G0 = R(G 0 ), G1 , . . . , GM H , . . . , GM = R(G M )
This is a sequence of sparse graphs thanks to lemma 3.5, and along this sequence, all edges in Ẽ are added one by one. We will show we can use this sequence to count the number of redundant edges in G.

Adding the edges of H: When adding one by one the edges of H, exactly n red (H) edges are discarded, and are added to the count of n red (G). After addition of all edges of H, GM H contains exactly M AB (H) edges between vertices A and B. Along this sequence of edges additions, ( Gn ) n≤M H contains only a number smaller than 2k -l of edges connecting A and B, and no other edge. Thus, the count of n red ( G) remains 0.

Adding an edge e ∈ E\E H (ie an edge in G, but not in H) : We start with the graphs G n = (V, E n ) and Gn = ( Ṽ , Ẽn ), with n ≥ M H . In this case, trying to add an edge in G corresponds to trying to add the same edge in G. We will show that e makes either both G n and Gn redundant, or none of the two.

Case 1 : (V, E n ∪ {e}) is not redundant: the added edge is accepted, G n+1 = (V, E n ∪ {e}), and the count of n red (G) is not modified. Then ( Ṽ , Ẽn ∪ {e}) is not redundant: otherwise, by lemma 3.5, it would be the image by renormalization of a redundant subgraph of (V, E n ∪ {e}), which is impossible. Then Gn+1 = ( Ṽ , Ẽn ∪ {e})), and neither the count of n red (G) nor the count of n red ( G) are modified.

Case 2 : (V, E n ∪ {e}) is redundant: the added edge is discarded, G n+1 = (V, E n ) = G n ,
and the count of n red (G) is increased by 1. Then ( Ṽ , Ẽn ∪ {e}) is redundant. Indeed, let g by a minimal redundant subgraph of (V, E n ∪ {e}). It cannot be included in H, as it contains the edge e. Then, by lemma 3.4, R(g) is a redundant subgraph of ( Ṽ , Ẽn ∪ {e}). Thus the count of n red ( G) is also increased by 1, and Gn+1 = Gn .

Repeating this until all edges in E\E H have been added proves formula (4.2). Figure 6. A hierarchical graph constructed from a "Wheatstone bridge" elementary cell. From left to right, the level is t = 0, t = 1, t = 2. This is the graph which is used in Section 4, where in addition each edge may have multiplicity up to 3.

Notice that formula (4.2) is wrong for the graph of Fig. 5, with (k, l) = (2, 1).

Hierarchical graphs.

If it is possible to iterate the renormalization transformation we have just defined until the graph becomes trivial, then the problems of percolation and counting the number of redundant edges are exactly solved. It is indeed possible to define some graphs on which this procedure can be carried out completely, thus providing exactly solvable models of (k, l) percolation that go beyond trees and random graphs. These graphs are called "hierarchical graphs" [19], and are defined as follows. We start from two vertices, connected by one edge. The graph is then constructed iteratively; at each step, all edges are replaced by a given elementary cell. From each type of elementary cell, one thus constructs a hierarchical graph. A graph where this replacement procedure has been iterated t times will be called a level t hierarchical graph. An example is given on Fig. 6.

Consider now within a hierarchical graph one elementary cell between vertices A and B. It clearly satisfies the renormalizability condition of Sec.3.2: all edges going outside of the cell are connected to the outer vertices A and B, whereas the inner vertices are only linked within the cell. Thus, the renormalization procedure described in Sec. 3.2 is exact when applied to an elementary cell. Furthermore, if the renormalization procedure is applied to all elementary cells of a level t hierarchical graph, the resulting renormalized graph is again a hierarchical graph, of level t -1. These remarks explain why the renormalization procedure allows to solve exactly (k, l)-rigidity models on hierarchical graphs.

An exactly solved model

In this section, we apply the formalism developed in the previous section to a body-bar rigidity model: that is we study (3,3) percolation on a multigraph G, with N vertices and M edges. The multigraph G is constructed starting from the "Wheatstone bridge" hierarchical lattice, described on Fig. 6.

We define the percolation problem as follows: each edge in the hierarchical lattice has multiplicity 0, 1, 2 or 3 with probability p 0 , p 1 , p 2 and p 3 = 1 -p 0 -p 1 -p 2 (an edge with multiplicity 0 is absent). Notice that if edges are either absent or have multiplicity 3 (that is p 1 = p 2 = 0), a subgraph is (3,3) rigid if and only if it is connected. In other words, (3, 3) rigidity with triple edges is equivalent to (1, 1) rigidity, that is ordinary percolation, with simple edges. For the special values of the parameters p 1 = p 2 = 0, this model thus contains ordinary percolation; since (3,3)

) = ϕ[(p (t) 0 , p (t) 1 , p (t) 2 , p (t)
3 )] From an analysis of the renormalization transformation, it is not difficult to obtain the explicit expression of ϕ, but it is very tedious. The details are given in the Appendix. Assuming that we start from a large graph with t 1, we iterate this renormalization transformation. To understand the phase diagram of the model, one needs now to study the renormalization flow induced by ϕ; of particular interest are the fixed points of ϕ. Note that the four dimensional space of parameters is actually easily reduced to three dimensions, since p 0 +p 1 +p 2 +p 3 = 1.

By inspection of the complicated expression for ϕ, three fixed points are easily found: i) p * flop = (p 0 = 1, p 1 = 0, p 2 = 0, p 3 = 0). This is the empty graph fixed point, corresponding to the floppy phase. ii) p * rig = (p 0 = 0, p 1 = 0, p 2 = 0, p 3 = 1). This is the full graph fixed point, corresponding to the rigid phase. iii) p * perc = (p 0 = 0.5, p 1 = 0, p 2 = 0, p 3 = 0.5). This is a critical fixed point, corresponding to ordinary percolation (because in this case, edges are either absent, or have multiplicity 3).

Using a Newton-Raphson scheme and scanning the whole parameter space, we have found another fixed point: iv) p * crit = (p 0 0.326, p 1 0.174, p 2 0.174, p 3 0.326). This fixed point may be found also by noting that the surface p 1 = p 2 , p 0 = p 3 is stable by ϕ. Using the normalization condition p 0 +p 1 +p 2 +p 3 = 1, looking for a fixed point on this surface is then a one dimensional problem. We have not found any other fixed point in the domain

p i ≥ 0, p 0 + p 1 + p 2 + p 3 = 1.
The trivial fixed points p * flop and p * rig are stable. The fixed point p * perc has three unstable directions. The fixed point p * crit has only one unstable direction. The renormalization flow is then as follows:

• The three dimensional parameter space is divided by a critical hypersurface containing p * perc and p * crit . On one side of the surface, the renormalization flows approaches the empty ("floppy", non percolating) fixed point; on the other side, it approaches the full fixed point ("rigid", percolating).

• On the critical hypersurface, the flow is attracted by p * crit , the "rigidity percolation" critical fixed point. We conclude that except for very special choices of parameters p 1 = p 2 = 0, the large scale critical properties of this model are described by the "rigidity percolation" critical fixed point. 4.2. Fixed points analysis. Let us call N red t (p) the number of redundant constraints in a level t hierarchical lattice with a multiplicity distribution p = (p 0 , p 1 , p 2 , p 3 ) for the edges. We introduce N t , B t the number of sites and edges in a level t hierarchical lattice. A simple computation yields B t = 5 t ; N t = 5 t + 3 2 We write r(p) = N red t (p)/N t the number of redundant constraints per site. We actually assume here that we work in the large t limit, and that a Law of Large Numbers holds for the random variable r(p). Following [16], we obtain the following expression

r(p) = +∞ k=0 g(ϕ (k) (p)) 5 k
where g(p) is a regular function associated with the mean number of redundant constraints in one unit cell with multiplicity distribution p, and ϕ (k) is the kth iterate of the function ϕ. We use the convention ϕ (0) =identity. The function r(p) has singularities at unstable fixed points of ϕ, and the exponent α associated to this singularity is related to λ 1 , the largest eigenvalue of the linearization of ϕ at the fixed point, by the equation α = 2 -ln 5/ ln λ 1 . We use here the standard notation for the exponent associated with the free energy singularity. Linearizing ϕ close to p * crit yields the eigenvalues (1.801, 0.568, 0.502), and the exponent α crit -0.74.

Linearizing ϕ close to p * perc yields a scalar matrix, and the eigenvalues (1.625, 1.625, 1.625). The corresponding exponent is α perc -1.32. In principle it is possible to compute the correlation length exponent ν through the hyperscaling relation dν = 2 -α. It is not clear which number should be used as dimension d for a hierarchical lattice. Using d = 2 yields ν crit 1.37, ν perc 1.66. The true exponents are ν 2D perc = 4/3 for 2D ordinary percolation and ν 2D crit = 1.21 ± 0.06 for rigidity percolation (see [7] for this numerical estimate). As expected, the critical exponents of the hierarchical model do not compare very well with the true 2D ones, but they are not completely off. In particular, the hierarchical model recovers ν crit < ν perc .

Finally, we note that the cross-over between standard and rigidity percolation is governed by the two scaling fields p 1 and p 2 , each one being associated with the same exponent ln(1.625) 0.486.

Conclusion

We have introduced and rigorously justified a renormalization transformation adapted to the study of (k, l)-rigidity, for k ≤ l < 2k, which generalizes the well known procedure for ordinary percolation, and the procedure for bar-joint rigidity used in [16]. This method allows to solve exactly (k, l)-percolation problems on hierarchical graphs. We have provided such an example, which has the interesting feature of showing both "ordinary percolation" and "rigidity percolation" behavior for different values of the parameters.

Rather than solving exactly problems on hierarchical graphs, such renormalization transformations might be used to provide approximate solutions for problems on more realistic 2D or 3D lattices. This work would then provide an approximate analytical tool in the general study of (k, l)-percolation problems and their universality classes, a field which is still widely open.

This work is supported by the ANR-09-JCJC-009401 INTERLOP project.

Appendix

We give here the details of the computations yielding the renormalization function ϕ. The problem is simple: apply the rule of Sec. 3.1 to an elementary "Wheatstone bridge" cell. Since there are 5 edges, and each edge may be absent, single, double or triple, there are 4 5 = 1024 different configurations of which we have to study the renormalization according to the rule of Sec. 3.1. The large number of configurations to enumerate is the only difficulty. Prop. 3.1 is useful to simplify these computations, as it allows to remove all redundant edges from the beginning.

In the following tables, we classify the configurations according to the multiplicity of their edges. For instance, an edge multiplicity 32210 means that one edge has multiplicity 3, two have multiplicity 2, one has multiplicity 1 and one is absent. Clearly, depending on how these edges are distributed on the "Wheatstone bridge", the number of renormalized edges may be different. As a consequence, there may be several lines with the same edge multiplicity distribution and different number of renormalized bonds; the first example of this is in the table is for multiplicity distribution 33320.The second column of the tables contains the number of edges of the renormalized cell, and the third column is the combinatorial factor corresponding to the number of configurations with the given edge multiplicity yielding the given number of renormalized edges. Fig. 7 gives some examples to better explain how the table is constructed.

We have collected the results in 6 tables, according to the number of edges with multiplicity 3.

  and the edges in E connect vertices in V . Definition 2.1. For k ≥ 1 and 0 ≤ l ≤ 2k -1, a multigraph G with N vertices and M edges is said to be

Figure 2 .

 2 Figure 2. Example of renormalization of a simple graph with 4 vertices, k = 2 and l = 3 (bar-joint rigidity percolation). The upper graph is renormalized into one edge, as A and B are rigidly connected; the lower graph is renormalized into the empty graph, as A and B are not rigidly connected.

Figure 3 .

 3 Figure 3. Example of renormalization of a simple graph with 4 vertices, k = 3 and l = 3 (body-bar rigidity percolation). The upper graph is renormalized into three edges, as the bodies A and B are rigidly connected; the lower graph is renormalized into two edges only. Note that the upper graph contains one redundant edge, for instance the dashed one.

Figure 4 .

 4 Figure 4. On the left is schematized a graph G; the subgraph H satisfies the renormalizability condition. On the right, H has been renormalized into one edge, the rest of the graph is not modified.

Figure 5 .

 5 Figure 5. For (k, l) = (2, 1), the subgraph H (induced by vertices (A, B, X, Y )) does not contain any redundant edge, and is renormalized into a subgraph with a single edge between A and B. The renormalized subgraph is sparse. However, the subgraph induced by the vertices (A, X, Y, Z, T ) is redundant.

5 1 4Figure 7 .

 17 Figure 7. (A): examples with multiplicity distribution 33210. Depending on the place of the multiple edges in the cell, the number of renormalized edges may be 3, 2, or 1. Note that the number of redundant constraints are then respectively 0, 1 and 2. Simple enumeration yields the corresponding number of configurations: 36, 12, and 12. (B): examples with multiplicity distribution 22110. The number of renormalized edges may be 1, or 0. Enumeration yields respectively 6 and 24 for the correposnding number of configurations.

  rigidity percolation is supposed to belong to a different universality class, it is interesting to study the whole phase diagram of this model.4.1. Renormalization flow.Let us start from a level t hierarchical lattice, with t large. The renormalization rule of Sec. 3 applied to each unit cell transforms the level t lattice to a level t-1 lattice, and induces a transformation of the parameters p i :

	(p	(t-1) 0	, p (t-1) 1	, p (t-1) 2	(t-1) , p 3

  Collecting the information from these tables, we obtain the expression for the renormalization function ϕ. Calling the renormalized probabilities p 0 , p 1 , p 2 , p 3 , we have

		No edge with multiplicity 3
	edges multiplicity in H nbr of renormalized edges nbr configurations
	22222	3	1
	22221	3	5
	22220	2	5
	22211	2	10
	22210	1	20
	22200	1	6
	22200	0	4
	22111	1	10
	22110	1	6
	22110	0	24
	22100	1	6
	22100	0	24
	22000	1	2
	22000	0	8
	21111	0	5
	21110	0	20
	21100	0	30
	21000	0	20
	20000	0	5
	11111	0	1
	11110	0	5
	11100	0	10
	11000	0	10
	10000 edges multiplicity in H nbr of renormalized edges nbr configurations 0 5 3 edges with multiplicity 3 00000 0 1
	33322	3	10
	33321	3	20
	33320	3	16
	33320	2	4
	33311	3	8
	33311	2	2
	33310	3	16
	33310	1	4
	33300	3	8
	33300	0	2

As a tree alone is never rigid, the leafs are usually connected to a busbar to study rigidity for trees.