Claudia Comi
email: claudia.comi@polimi.it

Jean-Jacques Marigo
email: jean-jacques.marigo@polytechnique.edu

J.-J Marigo

Homogenization approach and Bloch-Floquet theory for band-gap prediction in 2D locally resonant metamaterials

This paper provides a detailed comparison of the two-scale homogenization method and of the Bloch-Floquet theory for the determination of band gaps in locally resonant metamaterials. A medium composed by a stiff matrix with soft inclusions with 2D periodicity is considered and the equivalent mass density of the homogenized medium is explicitly obtained both for in-plane and out-of-plane wave propagation through two-scale asymptotic expansion. The intervals of frequency where the effective mass is negative identify the band gaps of the material. The Bloch-Floquet problem is then considered and, through an asymptotic analysis, its is shown that it leads to the same prediction of the band gaps. The results are confirmed by some examples and the limits of the asymptotic approach are explicitly given and numerically verified.

Introduction

Locally resonant metamaterials (LRM) are periodic or quasi-periodic media constituted by different materials, which can exhibit a peculiar dynamic behavior with spectral gaps at frequencies well below those characteristic of Braggs diffraction phenomena. The existence of band-gaps, i.e. frequency ranges where the propagation of real elastic waves is forbidden, is due to the presence of local resonances and has recently attracted a lot of interest in view of a number of potential applications, e.g. for the innovative design of vibration insulators [START_REF] D'alessandro | Modelling and experimental verification of a single phase three-dimensional lightweight locally resonant elastic metamaterial with complete low frequency bandgap[END_REF], [START_REF] Xu | Tailoring visbration suppression bands with hierarchical metamaterials[END_REF], of seismic shields [START_REF] Colombi | Forests as a natural seismic metamaterial: Rayleigh wave bandgaps induced by local resonances[END_REF], of impact absorbers [START_REF] Comi | Metamaterials for crashworthiness of small cars[END_REF] and of energy harvesters [START_REF] Lv | Vibration energy harvesting using a phononic crystal with point defect states[END_REF]. Various examples of LRM have been reported in the literature, they can be constituted by a periodic structure, made of a single material, with local resonators within the unit cell [START_REF] D'alessandro | Modelling and experimental verification of a single phase three-dimensional lightweight locally resonant elastic metamaterial with complete low frequency bandgap[END_REF], by three components (e.g. hard inclusions with a soft coating embedded in a hard matrix [START_REF] Liu | Locally Resonant Sonic Materials[END_REF][START_REF] Krushynska | Towards optimal design of locally resonant acoustic metamaterials[END_REF] or cellular materials filled by a foam with hard inclusions [START_REF] Comi | Wave propagation in cellular locally resonant metamaterials[END_REF]) or by two materials, typically soft inclusions in hard matrix, as studied e.g. in [START_REF] Wang | Two-dimensional locally resonant phononic crystals with binary structures[END_REF][START_REF] Hsu | Lamb waves in binary locally resonant phononic plates with two-dimensional lattices[END_REF]. In this latter case of binary composition, which is also considered in the present work, the metamaterials are also referred to as high-contrast composite materials [START_REF] Auriault | Dynamique des composites élastiques périodiques[END_REF].

The study of wave propagation in LRM, aimed to determine the dispersive properties and the presence of band-gaps, can be tackled in different ways. The two-scale homogenization approach was first proposed in [START_REF] Auriault | Dynamique des composites élastiques périodiques[END_REF] for high-contrast binary elastic composite materials in the long wavelength regime, and then developed in different directions. The extension of the the approach to three-component composite materials was proposed in [START_REF] Auriault | Long wavelength inner-resonance cut-off frequencies in elastic composite materials[END_REF], while in [START_REF] Marigo | Effective Dynamic Properties of a Row of Elastic Inclusions: The Case of Scalar Shear Waves[END_REF][START_REF] Pham | Two scale homogenization of a row of locally resonant inclusions -the case of anti-plane shear waves[END_REF] the approach for a row of locally resonant inclusions was developed. In all the above referenced works, the equivalent homogenized medium is characterized by an equivalent mass density that depends on the frequency and that can become negative in the vicinity of the resonant frequencies of the inclusions, thus identifying the stop bands. To overcome the classical hypothesis of wavelength much greater then the characteristic lattice length, the high-frequency asymptotic procedure for homogenization was proposed in [START_REF] Craster | High-frequency homogenization for periodic media[END_REF] and then used in [START_REF] Antonakakis | Homogenisation for elastic photonic crystals and dynamic anisotropy[END_REF] to analyze wave propagation through elastic media that contain periodic arrangements of traction free, or clamped, inclusions.

A different approach, widely employed in various contexts, considers the wave propagation in an infinite periodic media and uses the Bloch-Floquet principles to calculate the dispersion surfaces and to determine the pass and stop bands [START_REF] Brillouin | Wave propagation in periodic structures[END_REF]. In this framework, a popular strategy is the direct numerical simulation, using finite elements, of a unit cell of the periodic material with Bloch-Floquet boundary conditions [START_REF] Åberg | The usage of standard finite element codes for computation of dispersion relations in materials with periodic microstructure[END_REF], [START_REF] Phani | Wave propagation in two-dimensional periodic lattices[END_REF].

The main advantage of homogenization methods with respect to direct numerical solutions of the Bloch-Floquet problem is that they can lead to analytical predictions of the stop and pass bands and thus are more suitable in the design of the metamaterial itself. They allow to obtain a priori information on the dependence of the band gaps on the mechanical parameters of the constituents and on the geometry and thus open the way to optimization of the metamaterial for specific applications. On the other hand, homogenization methods are restricted by several hypothesis and the actual range of application is not fully defined.

To the author knowledge, there are very few studies devoted to the comparison of the two approaches. A significant exception is the work of Allaire and coworkers [START_REF] Allaire | A comparison between two-scale asymptotic expansions and Bloch wave expansions for the homogenization of periodic structures[END_REF] which, however, is focused on the discussion of the equivalent forth order stiffness tensor and not on the band gap prediction.

The aim of the present paper is to establish the link between the two-scale homogenization method and of the Bloch-Floquet theory for the determination of band gaps in 2D LRM. To this purpose in Section 3 we revise the homogenization approach for in-plane and out-of-plane wave propagation in a binary periodic material with circular soft inclusions in a stiff matrix and we obtain the effective mass density in closed form: the intervals of frequency for which the effective mass density becomes negative identify the band gaps. In Section 4, the Bloch-Floquet problem is considered and its properties are presented. The asymptotic analysis of that problem, for given wave number, provides an estimate of the band gaps which turns out to coincide with that obtained through homogenization. The results obtained are illustrated in Section 5 by some examples. The dispersion plots are numerically computed by finite element analysis on a single cell, imposing the Bloch-Floquet boundary conditions with wave vectors belonging to the boundary of the first Brillouin zone and the band gaps thus obtained are compared with the intervals defined with the homogenization approach. The limit of validity of these latter is also set and verified in the numerical examples.

Framework

We consider the wave propagation problem in a two-constituents elastic material with 2D periodicity. The main assumptions are the following ones:

1. The body Ω is heterogeneous, made of two materials: the matrix m and the cylindrical inclusions (also referred to as fibers, f ) periodically distributed in the plane x 1x 2 , see Figure 1a; 2. The two dimensional Bravais lattice of the material, shown in Figure 1b is characterized by the primitive vectors a 1 and a 2 (which can be non-orthogonal and have different moduli). The unit cell Y of the lattice is defined as (see Figure 1c):

Y = {x ∈ R 2 : x = 2 i=1 α i a i , -1/2 ≤ α i < 1/2 ∀i ∈ {1, 2}};
3. The reciprocal lattice is also 2D, its two basis vectors are denoted by A i , i = 1, 2 and defined from the a i 's by the relations

A i • a j = 2πδ ij , for i, j = 1, 2
where δ ij denotes the Kronecker symbol. The Wigner-Seitz cell of the reciprocal lattice defines the First Brillouin Zone, as shown in Figure 1d. 4. The size a = a 1 ∧ a 2 of the unit cell Y is small with respect to the considered wavelengths in the matrix material L m : a = L m with a dimensionless small parameter; 5. The fiber Y f is constituted by a linear elastic isotropic material with low stiffness. The Lamé's constants of this material are 2 λ f and 2 µ f , its volume mass density is ρ f . The radius of the inclusion is R = R; 6. The remaining part Y m of the cell is made of a stiff isotropic material with elastic moduli λ m , µ m and mass density ρ m ; 7. The characteristic length of the body in the out of plane direction x 3 is very large with respect to the size a.

x 1 The last assumption allows to treat the body as unbounded in the x 3 direction and hence to decouple the problem of wave propagation into two separate problems, one for longitudinal and shear in-plane propagation, and the other for shear waves with out-of-plane polarization. This decoupling is valid for both the homogenization and the Bloch-Floquet approaches, hence in the following we will discuss first the in-plane problem and then we will give the results for the out-of-plane problem.

x 2 x 3 ⌦ S x 1 x 2 a 1 a 2 Y ✏ m Y ✏ f R ✏ Y ✏ unit cell First Brillouin Zone (a) (b) (c) (d) A 1 A 2
Notation. Vectors and tensors are represented by bold face letters. An upper index T denotes transposed quantities. Summation over repeated indices is never implicitly assumed. The complex numbers (or complex fields) are denoted with sans serif letters, like i for √ -1; the complex conjugate of a complex number is indicated by a superposed bar; the modulus of the complex number u is denoted by |u|. The angular frequency ω is often called just frequency in the text for the sake of brevity.

3 Homogenization approach 3.1 In-plane wave propagation

Main assumptions

We consider a plane motion at given angular frequency ω. Specifically, the displacement field has a vanishing component in the direction 3 and components in the other two directions depending only on (x 1 , x 2 ) (plane strain assumption). That motion is governed by the Helmholtz equation

divσ + ρ ω 2 u = 0 in Ω (1)
where the stress field σ is related to the small strain deformation field by the elastic constitutive relation

σ = E ε in Ω (2) with ε = 1 2 (grad u + grad T u ) in Ω (3)
In ( 1)-( 2), E and ρ denote the (periodic) repartition of the elastic stiffness tensor and the density in the body. In view of the isotropy of the two constituents materials, E only depends on their Lamé's constants λ and µ and (2) in index notation (i,j = 1, 2) reads

σ ij = λ (ε 11 + ε 22 )δ ij + 2µ ε ij in Ω. (4)
Specifically, one has

λ = λ m in Y m 2 λ f in Y f , µ = µ m in Y m 2 µ f in Y f , ρ = ρ m in Y m ρ f in Y f (5)
In such a plane strain context, the other components of the stress field are given by

σ 13 = σ 23 = 0, σ 33 = λ (ε 11 + ε 22 ).
Therefore, throughout the present section, all the indices will range in {1, 2} and all the vectors and tensors will be considered as two-dimensional. The cross-section of the three dimensional domain Ω is denoted by S.

Asymptotic expansions and some properties of the first order terms

The hypotheses listed above allows to use the two-scale asymptotic method [START_REF] Allaire | Homogenization and 2-scale convergence[END_REF][START_REF] Sanchez-Palencia | Introduction aux méthodes asymptotiques et à l'homogénéisation[END_REF] to construct the homogenized (or effective) equation of motion of the body. Denoting by x = (x 1 , x 2 ) the (in-plane) macroscopic coordinate, let us introduce the microscopic one y = x/ which lives in the rescaled cell Y = Y / . The classical two-scale asymptotic method consists in expanding the field u in the following form:

u (x) = u 0 (x, x/ ) + u 1 (x, x/ ) + 2 u 2 (x, x/ ) + • • • (6)
where the fields u i = u i (x, y) are defined on S × R 2 and are Y -periodic with respect to y. In the same manner, the stress field σ is expanded in

σ (x) = σ 0 (x, x/ ) + σ 1 (x, x/ ) + 2 σ 2 (x, x/ ) + • • • (7)
where the fields σ i = σ i (x, y) are defined on S × R 2 and Y -periodic with respect to y. Note that one assumes in (7) that the stress expansion starts at order 0 while it should start a priori at order -1 by virtue of (2), ( 3) and [START_REF] Auriault | Long wavelength inner-resonance cut-off frequencies in elastic composite materials[END_REF]. In fact, it is possible to prove that σ -1 vanishes and assuming that property allows us to simplify the presentation. Since the terms of the expansions depend on the microscopic and macroscopic coordinates, we must distinguish the partial derivatives with respect to x or y. Specifically, for a vector function u of (x, y), one denotes by grad x u and ε x (u) the gradient of u and its symmetric part with respect to x, and by grad y u and ε y (u) the gradient of u and its symmetric part with respect to y. In the same manner, for a tensor function σ of (x, y), one denotes by div x σ and div y σ the divergence of σ with respect to x and y, respectively. In the case where the function depends on one variable only and when there exists no risk of confusion, the gradient and the divergence are denoted classically, without the index.

In the rescaled cell Y , one denotes by Y m the part occupied by the hard material and by Y f the part occupied by the soft material of radius R = R / , see Figure 2. 2 The rescaled cell of the periodic material with a compliant circular inclusion inside a stiff matrix From the expansion (6), one gets

Y Y f R Y m n Fig.
ε(u ) = -1 ε y (u 0 ) + ε y (u 1 ) + ε x (u 0 ) + ε y (u 2 ) + ε x (u 1 ) + • • • (8)
and from the expansion [START_REF] Brillouin | Wave propagation in periodic structures[END_REF], one gets divσ = -1 div y σ 0 + (div y σ 1 + div x σ 0 ) + (div y σ 2 + div x σ 1 ) + • • •

The above expansions of the displacement, (6), of the stress, [START_REF] Brillouin | Wave propagation in periodic structures[END_REF], and of their derivatives, (8) and [START_REF] Comi | Metamaterials for crashworthiness of small cars[END_REF], are substituted into the governing equations (1)-(3), then the terms of the same order in are considered separately. Inserting [START_REF] Brillouin | Wave propagation in periodic structures[END_REF] and ( 8) into (2)-(3) gives at order -1

ε y (u 0 ) = 0 in Y m
which, combined with the periodicity conditions, implies that u 0 depends only on x in S × Y m , say

u 0 (x, y) = U 0 (x) in S × Y m . (10)
Equation ( 1) gives, at order -1 in :

div y σ 0 = 0 in S × Y
whereas the expansion of (2)-(3) gives at order 0

σ 0 = E m ε x (U 0 ) + ε y (u 1 ) in S × Y m 0 in S × Y f ( 11
)
where E m is the isotropic elastic stiffness tensor of the matrix. The term ε 0 ≡ ε x (U 0 ) can be interpreted as a constant eigenstrain within the cell of independent components ε 0 ij (namely ε 0 11 , ε 0 22 and ε 0 12 = ε 0 21 in the plain strain conditions here considered). By virtue of (10) and using the periodic conditions that u 1 and σ 0 n must satisfy on the boundary of Y (n denoting the unit outer normal), by linearity, one deduces that u 1 can be written in Y m as

u 1 (x, y) = 2 i,j=1 ε 0 ij (x)χ ij (y) + U 1 (x) in S × Y m , (12)
where χ ij , for i, j = 1, 2, are the displacement fields solutions of the so-called elementary cell problems:

     div y (E m ε y (χ ij )) = 0 in Y m E m (ε y (χ ij ) + e i e j ) n = 0 on ∂Y f χ ij periodic, (E m ε y (χ ij ))n anti-periodic on ∂Y (13)
where e i is the unit vector in the direction i and denotes the symmetric tensorial product (2e i e j = e i ⊗ e j + e j ⊗ e i ). The solution of each elementary cell problem is defined up to a constant which can be fixed by imposing that Ym χ ij (y)dy = 0.

The elementary cell problems are classical in homogenization theory, note however that here problems (13) correspond to the case where the soft material is replaced by a cavity.

By virtue of ( 11) and [START_REF] Comi | Wave propagation in cellular locally resonant metamaterials[END_REF], the stress field σ 0 can read as

σ 0 (x, y) =      2 i,j=1 E m (ε y (χ ij )(y) + e i e j )ε 0 ij (x) in S × Y m 0 in S × Y f . (14)
Taking the average value of σ 0 on Y and using (13) one obtains

σ 0 (x) := 1 |Y | Y σ 0 (x, y)dy = E * ε 0 (x) (15)
where E * denotes the fourth order tensor whose components are given by

E * ijkl = 1 |Y | Ym E m ε y (χ ij ) + e i e j • ε y (χ kl ) + e k e l dy. (16)
In [START_REF] Liu | Locally Resonant Sonic Materials[END_REF] the dot denotes the tensor scalar product:

A • B = 2 r,s=1 A rs B rs .
The tensor E * represents the in-plane effective stiffness tensor of the matrix with a cavity inside. That tensor is symmetric and definite positive. It corresponds to the stiffness of an elastic anisotropic material, which can be isotropic (in the plane) if the cell is sufficiently symmetric.

In-plane eigenvibrations of the inclusion fixed at its boundary

Before to continue the asymptotic analysis of the equation of motion, let us characterize the in-plane eigenmodes of vibration and the associated eigenfrequencies of an inclusion with Dirichlet boundary conditions. The eigenmode η * = 0 and the associated eigenfrequency ω * are solution of

(λ f + µ f )grad y div y η * + µ f ∆ y η * + ρ f ω 2 * η * = 0 in Y f η * = 0 on ∂Y f , (17)
where

∆ y ≡ ∂ 2 ∂y 2 1 + ∂ 2 ∂y 2 2
is the two dimensional Laplace operator. From the classical spectral theorem for self-adjoint compact operators, it follows that: (i) there exists a countable family of eigenfrequencies; (ii) the eigenmodes associated to an eigenfrequency constitute a linear space of finite dimension (called the eigenspace); (iii) the family of all the eigenmodes constitutes an orthogonal basis (in the sense of the L 2 (Y f ) inner product) of the space of smooth in-plane displacement fields defined in Y f and vanishing on ∂Y f .

In the present case of a circular inclusion Y f , that eigenvalue problem can be solved in closed form. We only give the main results without the details of the calculations. In view of the geometry, polar coordinates r and ϑ are adopted:

y = r cos ϑe 1 + r sin ϑe 2 , r ∈ [0 , R], ϑ ∈ [0 , 2π].
The eigenmodes can be expressed as

η * (y) = f * (r) cos(pϑ + θ 0 )e r + g * (r) sin(pϑ + θ 0 )e ϑ ( 18
)
with p ∈ N, θ 0 an arbitrary constant, e r and e ϑ the unit vectors in the radial and tangential directions.

The two functions of r are given by

     f * (r) = a * J p-1 (k s * r) + J p+1 (k s * r) + b * J p-1 (k l * r) -J p+1 (k l * r) g * (r) = a * J p-1 (k s * r) -J p+1 (k s * r) + b * J p-1 (k l * r) + J p+1 (k l * r) (19)
where J p denotes the Bessel function of first kind and order p, a * and b * are constants. In [START_REF] Pham | Two scale homogenization of a row of locally resonant inclusions -the case of anti-plane shear waves[END_REF], k l * and k s * are the rescaled longitudinal and shear wave numbers for the inclusion, corresponding to the eigenfrequency ω * and defined as

k l * = ω * ρ f λ f + 2µ f = ω * c l , k s * = ω * ρ f µ f = ω * c s (20)
where c l and c s are the scaled velocities of propagation of longitudinal and shear waves. To fulfill the boundary condition η * = 0 on ∂Y f , the pair (a * , b * ) = (0, 0) must satisfy

J p-1 (k s * R)a * = -J p-1 (k l * R)b * , J p+1 (k s * R)a * = J p+1 (k l * R)b * , this is possible if ω * is a solution of den p ≡ J p-1 (k s * R)J p+1 (k l * R) + J p-1 (k l * R)J p+1 (k s * R) = 0. ( 21
)
This latter equation gives the eigenfrequencies with mode number p; they constitute a countable family of positive real numbers ω p n , n ∈ N * . Moreover, for p > 0, the eigenspaces are (at least) of dimension 2. As far as the properties of orthogonality of the eigenmodes with the translations is concerned, one must distinguish the case p = 1 from the other ones.

When p = 1, the eigenmodes are orthogonal to the translations, i.e. the following equalities hold true:

Y f η * (y) • e 1 dy = Y f η * (y) • e 2 dy = 0. ( 22
)
When p = 1, then an eigenfrequency ω * is solution of the following equation:

J 0 (k s * R)J 2 (k l * R) + J 0 (k l * R)J 2 (k s * R) = 0. ( 23
)
Accordingly, there exists a sequence denoted ω n , n ∈ N * (where the superscript p = 1 is omitted to simplify the notation), of eigenfrequencies solutions of [START_REF] Xu | Tailoring visbration suppression bands with hierarchical metamaterials[END_REF]. Those frequencies ω n can be expressed as

ω n = ζ n c l R (24)
with ζ n real numbers depending on the ratio between the longitudinal and the shear velocities, i.e. on the Poisson's coefficient of the inclusion ν f The two corresponding eigenmodes are given by:

c l c s = 2(1 -ν f ) 1 -2ν f := c(ν f ). (25)
η 1 * (y) = f * (r) cos ϑe r -g * (r) sin ϑe ϑ , η 2 * (y) = f * (r) sin ϑe r + g * (r) cos ϑe ϑ (26)
where the two functions f * and g * read

     f * (r) = a * J 0 (k s
So, the eigenmodes are orthogonal to the translations if and only if J 2 (k l * R)J 2 (k s * R) = 0, condition which can be written as

J 2 (ζ n )J 2 c(ν f )ζ n = 0 (29)
and depends on the Poisson ratio ν f and the number of the eigenfrequency. In general, the two n th eigenmodes are not orthogonal to the translations, except for particular values of the Poisson ratio. Numerical computations show that for n = 1 and n = 2, equation (29) has no solution in the range 0 ≤ ν f < 0.5 (see Figure 4a for n = 1), hence the first (n = 1) and the second (n = 2) eigenmodes are never orthogonal to the translations. On the contrary, as one can see on Figure 4b, both J 2 (ζ 3 ) and J 2 (c ζ 3 ) vanish for the value ν f = ν * ≈ 0.2035 and hence the third eigenmodes are orthogonal to the translations for

ν f = ν * . 0.1 0.2 0.3 0.4 0.2 J 2 ⌫ f J 2 (⇣ 1 ) J 2 (c⇣ 1 ) (a) 0.1 0.2 0.3 0.4 -0.2 0.2 J 2 ⌫ f ⌫ ⇤ J 2 (c⇣ 3 ) J 2 (⇣ 3 ) (b)

The motion problem for the soft inclusion

The expansion of the equation of motion [START_REF] Åberg | The usage of standard finite element codes for computation of dispersion relations in materials with periodic microstructure[END_REF] gives at order 0 in the inclusion

div y σ 1 + ρ f ω 2 u 0 = 0 in S × Y f (30)
where we have used the fact that σ 0 = 0 in Y f . The expansion of the constitutive relation at order 1 in

Y f gives σ 1 = E f ε y (u 0 ) in S × Y f , (31)
where E f is the rescaled stiffness tensor of the inclusion. Inserting (31) into (30) and using the continuity of the displacement at the interface between the soft and the hard material lead to the following local problem governing the motion of the inclusion:

(λ f + µ f )grad y div y u 0 + µ f ∆ y u 0 + ρ f ω 2 u 0 = 0 in S × Y f u 0 = U 0 on S × ∂Y f . ( 32
)
For a given ω and a given U 0 , problem (32) admits a unique solution u 0 provided that ω is not a frequency of resonance of the clamped inclusion. Specifically, in that case, by linearity, u 0 can read as

u 0 (x, y) = 2 i=1 U 0 i (x)η i (y) in S × Y f ( 33
)
where η i , for i = 1 or 2, must satisfy (at given ω) the following PDE:

(λ f + µ f )grad y div y η i + µ f ∆ y η i + ρ f ω 2 (e i + η i ) = 0 in Y f (34)
with the boundary conditions η i = 0 on ∂Y f . (35) Let us first discuss the existence and the uniqueness of the solution of problem (34)-(35) depending on the value of ω. One has the following alternatives:

(i) If ω is not an eigenfrequency, the solution of (34)-( 35) is unique.

(ii) If ω is an eigenfrequency whose eigenspace is orthogonal to the translations, the problem (34)-(35) admits a (non unique) solution defined up to an arbitrary eigenmode. (iii) If ω is an eigenfrequency with some eigenmodes non orthogonal to the translations, the problem (34)-( 35) does not admit a solution. The proofs of those properties are direct consequences of the modal decomposition of (34). This problem can be solved in closed form in the case of a circular inclusion.

In case (i), where ω is not an eigenfrequency, the unique solution is of the form

η 1 (y) = η r (r) cos ϑe r -η ϑ (r) sin ϑe ϑ , η 2 (y) = η r (r) sin ϑe r + η ϑ (r) cos ϑe ϑ . ( 36
)
One deduces from (34) that the two functions of r, η r and η ϑ , are given by

     η r (r) = a J 0 (k s r) + J 2 (k s r) + b J 0 (k l r) -J 2 (k l r) -1 η ϑ (r) = a J 0 (k s r) -J 2 (k s r) + b J 0 (k l r) + J 2 (k l r) -1 (37)
In (37), k l and k s are the rescaled longitudinal and shear wave numbers for the inclusion, corresponding to a given angular frequency ω and defined as

k l = ω c l , k s = ω c s . (38)
By virtue of (35), the constants a and b must be such that η r (R) = η ϑ (R) = 0. One deduces from that linear system that the constants a and b are given by

a = J 2 (k l R) den , b = J 2 (k s R) den (39)
provided that the denominator den which reads

den ≡ J 0 (k s R)J 2 (k l R) + J 0 (k l R)J 2 (k s R) (40)
does not vanish. Accordingly, by virtue of ( 23), a and b are well defined by (39) provided that ω = ω n , n ∈ N * , ω n being the eigenfrequencies of the inclusion associated with the eigenmode number p = 1.

In the case (ii) one should distinguish two situations: case (iia) the eigenfrequency ω * does not satisfy [START_REF] Xu | Tailoring visbration suppression bands with hierarchical metamaterials[END_REF]; case (iib) the eigenfrequency ω * satisfies [START_REF] Xu | Tailoring visbration suppression bands with hierarchical metamaterials[END_REF] and its eigenspace is orthogonal to the translations.

In the case (iia) which corresponds to an eigenmode number p = 1, the solution of (34)-( 35) is defined up to an eigenmode. Accordingly, the solution reads as

η 1 (y) = η r (r) cos ϑe r -η ϑ (r) sin ϑe ϑ + η 1 * (y), η 2 (y) = η r (r) sin ϑe r + η ϑ (r) cos ϑe ϑ + η 2 * (y) (41)
where η r and η ϑ are still given by (37)-(40) whereas η 1 * and η 2 * are two arbitrary eigenmodes. In the case (iib), the eigenmode number p = 1, the solution of (34)-(35) still exists and is defined up to an eigenmode (its general expression is not given here). Such a case holds for instance when ν f = ν * ≈ 0.2035 and ω = ω 3 as discussed in Section 3.1.3. In this particular case where both J 2 (k s R) and J 2 (k l R) vanishes, the solutions η 1 and η 2 of (34)-(35) are given by (36)-(37) with a and b related by

aJ 0 (k s R) + bJ 0 (k s R) = 1.
(42) Finally, when ω is an eigenfrequency which satisfies [START_REF] Xu | Tailoring visbration suppression bands with hierarchical metamaterials[END_REF] and whose eigenspace is not orthogonal to the translation (case (iii)), the problem (34)-(35) has no solution.

The effective in-plane motion problem for the composite

One can only consider the cases (i) and (ii) where the problem (34)-(35) admits a solution, because in the case (iii) where the problem (34)-(35) has no solution, from (32) one should have U 0 = 0.

The expansion of the equation of motion gives at the order 0 in the whole cell

div y σ 1 + div x σ 0 + ρω 2 u 0 = 0 in S × Y. (43)
Integrating this equation over Y leads to the following expressions for the different terms:

Y div y σ 1 dy = ∂Y σ 1 n dS = 0 by virtue of the periodic conditions div x σ 0 = div σ 0 = div E * ε(U 0 ) by virtue of ( 15)

ρu 0 = ρ * (ω)U 0
by virtue of ( 10) and ( 33)

where ρ * (ω) is a 2 × 2 matrix whose entries are given by

ρ * ij (ω) = ρ st δ ij + ρ f |Y | Y f η i j (y)dy. ( 44
)
In (44) ρ st denotes the equivalent static mass density of the material,

ρ st = ρ m |Y m | |Y | + ρ f |Y f | |Y | . ( 45
)
Finally the effective equation of motion reads

div E * ε(U 0 ) + ρ * (ω)ω 2 U 0 = 0 in S. (46)
Thus, ρ * (ω) can be interpreted as the matrix of the (in-plane) effective mass density of the composite material. It depends on the frequency because of the dependence of η i on ω. Let us first note that ρ * (ω) is well defined even when the solution η i of (34)-( 35) is not unique. Indeed, that corresponds to the case (ii) where ω is an eigenfrequency whose eigenspace is orthogonal to the translations and then η i is defined up to an eigenmode η * . Therefore, since Y f η * •e 1 dy = Y f η * •e 2 dy = 0, Y f η i j dy is determined. Let us also note that ρ * (ω) does not depend on the stiffness of the hard material, but only on the volume fraction and the mass density of the materials, and on the geometry and the stiffness of the soft inclusion.

In the case of a circular inclusion, ρ * (ω) is obtained in closed form. Let us first consider the cases (i) and (iia) of the previous subsection 3.1.4. Using (36) one easily checks that Y f η i j (y)dy = 0 for i = j and that ρ * 11 = ρ * 22 . Therefore ρ * (ω) can read as

ρ * (ω) = ρ * (ω)I.
Using (37)-(40) one obtains

ρ * (ω) = ρ st + ρ f |Y f | |Y | 2J 2 (k l R)J 2 (k s R) J 0 (k s R)J 2 (k l R) + J 0 (k l R)J 2 (k s R) for ω = 0. ( 47
)
When ω = 0, one finds directly from (34)-( 35) that η i = 0 in Y f and hence ρ * (0) = ρ st .

Figure 5 shows the variation of the effective mass density (normalized with the static one) with the normalized frequency ω/ω 1 . Note that this normalized plot only depends on ν f , which is fixed to 0.4 in Figure 5. Thus, ρ * (ω) is not defined when ω = ω n = ζ n c l /R with n ∈ N * , tends to -∞ and +∞ when ω tends to ω n from above and from below, respectively, and it is strictly increasing in every interval (ω n , ω n+1 ). Consequently, the effective mass density is not positive in a countable family of intervals I n of the form where θ n are the frequencies corresponding to ρ * = 0.

I n = ζ n c l R , θ n , n ∈ N * , (48
In the particular case (iib) of the previous subsection 3.1.4 when ν f = ν * ≈ 0.2035 and ω = ω 3 , using (36), (37) and ( 42) one has Y f η i j dy = 0 and hence ρ * (ω 3 ) = ρ st . Figure 6 compares the variation of the effective mass density in the cases ν f = ν * = 0.2035 (left) and ν f close (but not equal) to ν * (right). When ν = ν * , even though ω 3 ≈ 2.6ω 1 is an eigenfrequency, since its eigenspace is orthogonal to the translations, ρ * (ω 3 ) is well defined (and equal to ρ st ) and there is no interval near ω 3 where ρ * is negative. Therefore, as it will be discussed in the next subsection, no band gap exists opening at ω 3 . If one considers ν f slightly different from ν * (equal to 0.23 in Figure 6 on the right), since the eigenspace associated with ω 3 is no more orthogonal to the translations, then ρ * (ω 3 ) is no more defined and it exists an interval on the right of ω 3 where ρ * is negative. Accordingly, there is a band gap opening at ω 3 . However the width of the band gap is small (within the line thickness in Figure 6, right). In fact the closer ν to ν * and the smaller the width of the band gap near ω 3 ; at the limit, when ν = ν * the band gap disappears.

The band gaps

Let us consider the effective equation of motion (46) and a solution corresponding to an harmonic wave of the form

U 0 (x, t) = exp i(km • x -ωt)ξ,
where ω and m denote the given angular frequency and direction of propagation whereas ξ and k are the non zero polarization vector and wave number to be determined. Such a wave with wave vector k = km can propagate inside the composite material provided that k and ξ satisfy the following problem:

Q * (m)ξ = ρ * (ω) ω 2 k 2 ξ, (49)
where Q * (m) denotes the effective acoustic tensor associated with the direction of propagation m:

Q * ih (m) = 2 j,l=1 E * ijhl m j m l . (50)
Since Q * (m) is positive definite, the wave number k is a real positive number if and only if the effective mass density ρ * (ω) is positive. Therefore, by virtue of the previous analysis, harmonic waves can only propagate if the frequency ω lies outside the intervals I n , n ∈ N * , defined by (48) which thus identify the band gaps.

Let us return to the real physical quantities instead of the rescaled ones. The real wave velocities in the inclusion are c l = c l , c s = c s and the real diameter of the inclusion is 2R = 2 R. Therefore c l /R = c l /R and the frequencies ω n previously introduced can be interpreted as characteristic frequencies of the real inclusion. Hence the intervals I n turn out to depend only on the mass ratio between the hard and the soft material and on the Poisson's coefficient of the inclusion. Figure 7 shows the dependence of the first band gap on ν f , the relative amplitude increases as ν f decreases. 0.5 1.0 1.5 7 Normalized effective mass vs. normalized frequency for ν f . The intervals of negative effective mass give the first band gap. The mass ratio between matrix and inclusion is 1.41 and corresponds to the material considered in Section 5.

-1 0 1 2 3 0. .1 .2 .3 ⌫ f = .4 ⇢ ⇤ ⇢ st ! ! 1 band gaps Fig.

Out-of-plane wave propagation

The same procedure followed for in-plane polarized waves can be used for waves with out-of-plane polarization.

Asymptotic expansions and some properties of the first order terms

Specifically, one assumes now that the displacement field has only a component in the third direction and that it depends only on (x 1 , x 2 ) (anti-plane strain condition):

u (x) = u (x 1 , x 2 )e 3 .

Out-of-plane eigenvibrations of the inclusion fixed at its boundary

An out-of-plane eigenmode η * = 0 and the associated eigenfrequency ω * are solution of

µ f ∆ y η * + ρ f ω 2 * η * = 0 in Y f η * = 0 on ∂Y f . (59)
As for the in-plane problem: (i) there exists a countable family of eigenfrequencies; (ii) the eigenmodes associated to an eigenfrequency constitute a linear space of finite dimension; (iii) the family of all the eigenmodes constitutes an orthogonal basis (in the sense of L 2 (Y f ) inner product) of the space of smooth out-of-plane displacement fields defined in Y f and vanishing on ∂Y f . In our case of a circular inclusion, that eigenvalue problem can be solved in closed form. For every p ∈ N, there exists a countable family of eigenfrequencies ω * which are solutions of

J p (k s * R) = 0, k s * = ω * c s . ( 60
)
The associated eigenmodes expressed in polar coordinates are

η * (y) = J p (k s * r)(a * cos pϑ + b * sin pϑ) (61)
where a * and b * are arbitrary constants (with a * = 0 when p = 0, and (a * , b * ) = (0, 0) when p = 0). As far as the property of orthogonality of the eigenmodes with respect to the constants (which correspond to the translations in the direction e 3 ) is concerned, one must distinguish the case p = 0 from the other ones.

(i) When p = 0, the eigenmodes are orthogonal to the constants, i.e. the following equality holds true:

Y f η * (y)dy = 0.
(ii) When p = 0, then the eigenmodes are purely radial and not orthogonal to the constants. Indeed, by virtue of ( 59)-( 60), one has

Y f η * (y)dy = a * π R 0 rJ 0 (k s * r)dr = - a * πR k s * J 0 (k s * R) = 0,
the last inequality being due to the fact that J 0 (x) = 0 when J 0 (x) = 0.

Accordingly, when p = 0, there exists a sequence of eigenfrequencies ω n , n ∈ N * , which can be expressed as

ω n = z n c s R ( 62
)
where z n are the roots of J 0 (z) = 0: Let us note that, unlike the in-plane case, these z n 's are independent from ν f .

z 1 = 2.

The out-of-plane motion problem for the soft inclusion

In the inclusion, the homogenized first order displacement u 0 must satisfy the following local problem :

µ f ∆ y u 0 + ρ f ω 2 u 0 = 0 in S × Y f u 0 = U 0 on S × ∂Y f (63)
For a given U 0 and a given ω which does not satisfy (60), problem (63) admits as unique solution u 0 that can be expressed as

u 0 (x, y) = U 0 (x)η(y) in S × Y f ( 64
)
where η is the solution of:

∆η + (k s ) 2 η = 0 in Y f η = 1 on ∂Y f (65)
In the case of a circular inclusion, (65) can be solved in closed form and η reads

η(y) = J 0 (k s r) J 0 (k s R) with r = y < R. ( 66
)
In the case where ω is an eigenfrequency ω * , but with J 0 (k s * R) = 0 (hence associated with an eigenmode number p = 0), the solution of ( 65) is defined up to an eigenmode and reads as

η(y) = J 0 (k s * r) J 0 (k s * R) + J p (k s * r)(a * cos pϑ + b * sin pϑ), (67)
where a * and b * are arbitrary constants.

In the case where ω is an eigenfrequency ω * such that J 0 (k s * R) = 0 (hence associated with a radial eigenmode), the problem (65) has no solution and the problem (63) has a solution only if U 0 = 0.

The effective out-of-plane motion problem for the composite

As explained in Section 3.1.2, the integration over the cell of the equation of motion gives at the order 0 an effective equation of motion (see (46)). Similarly, for out-of-plane polarized waves, that equation becomes div

M * grad U 0 + ρ * * (ω)ω 2 U 0 = 0 in S ( 68
)
where M * denotes the second order (out-of-plane) effective stiffness tensor whose components are given by

M * ij = 1 |Y | Ym µ m grad y χ i + e i • grad y χ j + e j dy, (69)
and ρ * * (ω) denotes the out-of-plane effective mass density given by

ρ * * (ω) = ρ m |Y h | |Y | + ρ f |Y | Y f η(y)dy. (70)
In the case of a circular inclusion, ρ * * (ω) is obtained in closed form. Using (66) and (67), one gets

ρ * * (ω) = ρ st + ρ f |Y f | |Y | J 2 (k s R) J 0 (k s R) . ( 71
)
So, ρ * * (ω) is not defined when ω is equal to some ω n defined by (62), and tends to +∞ (resp. -∞) when ω tends to ω n from below (resp. from above). Accordingly, there exists a band gap Ĩn of the form

Ĩn = z n c s R , ϕ n , n ∈ N * , ( 72
)
where ϕ n are the frequencies corresponding to ρ * * = 0. Those band gaps can be seen in Figure 8.

1 2 3 4 5 -1 0 1 2 3 ⇢ ⇤⇤ ⇢ st ! ! 1
Fig. 8 Normalized effective mass vs. normalized frequency for out-of-plane polarized waves. The intervals of negative effective mass give the band gaps. The mass ratio between matrix and inclusion is 1.41 and corresponds to the material considered in Section 5.

4 Bloch-Floquet approach

Setting of the problem

The fundamental problem, denoted here below by P(k), consists in finding, for a given wave vector k, a non trivial (complex) displacement field u, the associated (complex) stress field s and a frequency ω such that

P (k)                div s + ρ ω 2 u = 0 in Y s = E ε(u) in Y 2ε(u) = grad u + grad u T in Y u(x + a i ) = exp(i k • a i )u(x) ∀x on ∂ - i Y , i = 1, 2 s(x + a i )n = exp(i k • a i )s(x)n ∀x on ∂ - i Y , i = 1, 2
Here above, ∂ - i Y and ∂ + i Y denote the opposite sides of the cell Y in the direction i. Accordingly, one has

∂ + i Y = ∂ - i Y + a i . Moreover, if n is the outer unit normal to ∂ + i Y , then -n is the outer unit normal to ∂ - i Y .
The main goal is to find the spectrum, that is to say, the set described by the frequencies when the wave vector k describes R 2 .

By a change of variable, the problem can be reformulated in order to work with periodic fields. Specifically, introducing û by

u(x) = exp(i k • x)û(x), ∀x ∈ Y ,
the problem consists now in finding a non trivial displacement field û, the associated stress field ŝ and a frequency ω such that

P (k)      div ŝ + i ŝ k + ρ ω 2 û = 0 in Y ŝ = E ε(û) + iû k
in Y û periodic and ŝn anti-periodic on ∂Y (73)

The fundamental problem can be written in a variational form. Indeed, let us introduce, for a given wave vector k, the Rayleigh ratio R k defined on the space H 1 # (Y ) of (complex) fields which are in H 1 (Y ) and Y -periodic:

R k (v) = Y E (ε(v) + iv k) • (ε( v) -i v k)dx Y ρ |v| 2 dx . (74)
The following proposition holds true:

Proposition 1 The displacement fields û solution of P (k) are the stationary points of the Rayleigh ratio R k and the associated frequencies ω are given by the Rayleigh ratio of û,

R k (û)(v) = 0 ∀v ∈ H 1 # (Y ), ω 2 = R k (û), (75)
where R k (û)(v) stands for the directional derivative of R k at û in the direction v.

The proof is classical and we only gives the main lines. First, taking the derivative of R k , the stationary condition reads as

Y E (ε(û) + iû k) • (ε( v) -i v k)dx = ω 2 Y ρ û • vdx, ∀v ∈ H 1 # (Y ). ( 76
)
Using the definition of ŝ, (76) becomes

Y ŝ • (ε( v) -i v k)dx = ω 2 Y ρû • vdx, ∀v ∈ H 1 # (Y ).
Making an integration by parts of the term ŝ • ε( v) above and using the fact that v is arbitrary in H 1 # (Y ), one gets by standard arguments of Calculus of Variations that ŝ satisfies the equation of motion in Y and ŝn is anti-periodic on ∂Y . Hence, û, ŝ and ω are solution of P (k). One proves the converse in the same manner.

Some general properties

Let us first indicate some general properties of the fundamental Bloch-Floquet problem.

1. The spectral theorem. Using the variational formulation of the fundamental problem, one gets the so-called spectral theorem:

Proposition 2 (Spectral theorem) At given k, the solutions of P (k) enjoy the following properties: (a) The eigen-frequencies ω constitute a sequence of real numbers ω n , n ∈ N, which tends to infinity without accumulation points; (b) The eigen-displacement fields associated with an eigen-frequency ω n constitute a subspace U n of H 1 # (Y ) of finite dimension; (c) The eigen-spaces are orthogonal to each other in the sense of the scalar product

(û, v) = Y ρ û • vdx and constitute a countable basis of H 1 # (Y ).
The proof is a variant of the classical one for self-adjoint compact operators, see [START_REF] Dautray | Analyse Mathématique et Calcul Numérique pour les Sciences et les Techniques[END_REF]. It is based on the minimization of the Rayleigh ratio.

2. Periodic dependence of ω on k. It is easy to verify that if û, ŝ and ω are solution of the problem P (k), then they are also solution of the problem P (k + A) for all A such that A = pA 1 + qA 2 with p, q ∈ N.

where A i are the basis vectors of the reciprocal lattice, see Figure 1d. The proof is essentially based on the fact that exp(iA • a 1 ) = exp(iA • a 2 ) = 1. Consequently, one can restrict the study of the spectrum to the wave vectors k lying in the First Brillouin Zone (BZ). 3. Smooth dependence of the solution of P (k) on k. By using classical results on linear operators, one knows that the solution of P (k), and hence the frequency ω, depends smoothly (in fact is a C ∞ function) of the wave vector k. 4. ω = 0 is always solution of P (0). Indeed, taking k = 0, one easily checks that û(x) = Û (a constant), ŝ = 0 and ω = 0 are solution of P (0). The eigen-space U 0 corresponds to the translations (the rotations are forbidden by the periodic conditions). 5. ω = 0 is not solution of P (k) when k ∈ BZ, k = 0. Indeed, by virtue of the positivity of the Rayleigh ratio, in order that ω = 0 be a minimizer of R k , one must have R k (û) = 0 and hence ε(û)

+ iû k = 0.
In terms of the original field u, one gets

ε(u) = exp(ik • x)(ε(û) + iû k) = 0.
Hence, u must be a rigid displacement field of the form u(x) = U + w ∧ x. The boundary conditions u(x

+ a i ) = exp(k • a i )u(x) on ∂ - i Y requires that w = 0 and U = exp ik • a i U, for i = 1, 2. Therefore, one must have exp ik • a i = 1 for i = 1, 2 which is possible only if k = 0. 6.
There is no band gap in the neighborhood of ω = 0. Indeed ω = 0 is solution when k = 0 but not when k = 0. If one considers the frequency which minimizes the Rayleigh ratio R k , say ω 0 (k), since it depends continuously on k and since

ω 0 (k) = 0 if k = 0 > 0 if k = 0 ,
there exists a neighborhood of 0 in which all the frequencies can be reached. 7. Dependence of the spectrum on the size of the cell. Let us change the cell Y by λY (homothetic transformation by the ratio λ > 0) and the wave vector k by k/λ. If û and ω are solution of the problem P (k) associated with the cell Y , then ûλ defined on λY by ûλ (x) = û(x/λ) and ω λ = ω/λ are solution of the problem P (k/λ) associated with the cell λY . Therefore, the spectrum associated with λY is obtained from that associated with Y by the factor 1/λ. In other words, a change of size does not modify the structure of the spectrum but simply changes the scale of the eigen-frequencies.

In particular, the existence and the structure of band gaps remain unchanged by changing the size of the cell. To change the structure of the band gaps, one must change the materials inside the cell.

Rescaling

Since, by virtue of the above Property [START_REF] Brillouin | Wave propagation in periodic structures[END_REF], the frequencies scale with the cell dimension, it is more convenient to consider the rescaled cell Y = Y / already introduced in Section 3.1.2, see Figure 2.

Recalling the definition of the rescaled coordinate y = x/ , at û and ŝ defined on Y one can associate û and ŝ defined on Y by û (y) = û( y), ŝ (y) = ŝ( y).

Introducing these rescaling into P (k) given in (73), the rescaled Bloch-Floquet problem consists in finding û = 0, ŝ and ω such that

P (k)                    -1 div ŝ + i ŝ k + ρ f ω 2 û = 0 in Y f ŝ = E f ε(û ) + iû k in Y f û and ŝ n continue on ∂Y f -1 div ŝ + i ŝ k + ρ m ω 2 û = 0 in Y m ŝ = -1 E m ε(û ) + iû k in Y m û periodic and ŝ n anti-periodic on ∂Y (77)
where we distinguish between the two constituents materials. In the equations above, all the derivatives are taken with respect to y. Note that the frequency has not been rescaled, ω represents the "true" angular frequency of the "true" cell Y .

In the next subsection, we make an asymptotic analysis of this rescaled Bloch-Floquet problem by considering in-plane wave vectors k and we compare the results with the homogenization approach.

4.4 Asymptotic analysis of P (k) at given k.

One considers an in-plane wave vector k = k 1 e 1 + k 2 e 2 and searches for solutions which admit the following expansions in terms of :

     û = û0 + û1 + • • • , û0 = 0, ŝ = ŝ0 + ŝ1 + • • • ω = ω + • • •
where we have omitted the superscript 0 for the leading term of the frequency expansion, to simplify the notation. In the expansion above, the displacements and the stresses are three-dimensional (with in-plane and out-of-plane components) but these components depend only on (y 1 , y 2 ).

-Inserting those expansions into P (k) given in (77), the constitutive relation at order -1 reads ε(û 0 ) = 0 in Y m and hence û0 (y

) = U 0 in Y m , (78)
U 0 being an arbitrary (three-dimensional) vector. -The constitutive equations in Y m and Y f give at order 0

ŝ0 = 0 in Y f E m ε(û 1 ) + iU 0 k in Y m (79)
whereas the equation of motion in Y m at order -1 reads

divŝ 0 = 0 in Y m .
Therefore, using the results of the Section 3, û1 can be written in the following form

û1 (y) = i 2 i,j=1 U 0 i k j χ ij (y) + i 2 i,j=1 U 0 3 k j χ i (y)e 3 + U 1 in Y m ( 80
)
where the vectors χ ij are the solutions of the elementary in-plane cell problems (13), the scalars χ i are the solutions of the elementary out-of-plane cell problems (57) and U 1 is an arbitrary (threedimensional) vector.

-In Y f , using the previous results, the constitutive equation at order 1 and the equation of motion at order 0 lead to the following problem for û0 :

     divŝ 1 + ρ f ω 2 û0 = 0 in Y f ŝ1 = E f ε(û 0 ) in Y f û0 = U 0 on ∂Y f (81)
which is a linear combination of The in-plane and out-of-plane problems (32) and (63). Indeed, by linearity (81) can be decomposed into a plane strain problem where û0 3 = U 0 3 = 0 and into an antiplane strain problem where û0 i = U 0 i = 0 for i = 1, 2. Therefore, one can use all the results from Section 3.1.4 up to Section 3.1.6 for the in plane problem, and from Section 3.2.3 to Section 3.2.4 for the out-of-plane problem.

-The equation of motion in the whole cell at order 0 reads as

divŝ 1 + iŝ 0 k + ρω 2 û0 = 0 in Y.
Integrating that equation over Y and using the periodic boundary condition give the following equation:

i ŝ0 k + ω 2 ρû 0 = 0. ( 82
)
Using ( 79)-( 81), (82) becomes

             for i = 1, 2, 2 j,p,q=1 E * ijpq k j k q -ρ * ip (ω)ω 2 U 0 p = 0 for i = 3, 2 j,q=1 M * jq k j k q -ρ * * (ω)ω 2 U 0 3 = 0 (83)
where E * is the effective in-plane stiffness tensor given by ( 16), ρ * (ω) is the effective in-plane mass density given by (44), M * is the effective out-of-plane stiffness tensor given by (69) and ρ * * (ω) is the effective out-of-plane mass density given by (70). In the case of a circular inclusion, ρ * (ω) = ρ * (ω)I with ρ * (ω) given by (47) whereas ρ * * (ω) is given by (71). So, we have found the same equations as in the homogenization approach. The end of the discussion to find the spectrum is essentially the same as in Section 3. We finally recover that the band gaps are the intervals I n and Ĩn , n ∈ N * , defined in Section 3, see (48) and (72).

Examples

The above results are illustrated in this section with reference to a 2D periodic material, characterized by a square elementary cell of side a constituted by an aluminium matrix with polymeric cylindrical inclusions of radius R = 0.415a, Figure 9a. The reciprocal lattice vectors A 1 and A 2 are equal and orthogonal, and the First Brillouin Zone is also a square. Figure 9b shows the Wigner-Seitz cell of the reciprocal lattice with the shaded triangle representing, account taken of the symmetries, the Irreducible Brillouin Zone (IBZ). In the following numerical analyses only wave vectors belonging to the boundary of the IBZ, along the path Γ -X -M -Γ will be considered as proposed e.g. in [START_REF] Phani | Wave propagation in two-dimensional periodic lattices[END_REF]. The properties of the constituents are given in Table 2.

(a) (b) a 1 a 2 a R m=matrix f =fiber X M A 1 A 2 2⇡/a

In-plane waves

The homogenization approach is based on the hypotheses listed in Section 2. In particular hypothesis (3) can be rewritten as

         a L m = aω ρ m λ m + 2µ m << 1 for longitudinal waves a L m = aω ρ m µ m << 1 for shear waves (84)
meaning that the homogenization approach, for a given material, applies only to waves with low frequencies, i.e. waves with ω << ω max = µm ρm /a. Therefore the prediction of the band gaps given by the homogenization approach is a reasonable estimate only if the corresponding frequencies given in eq. ( 48) are smaller than ω max . For the constituents listed in table 2, and a = 1 mm, f max = 1 2π ω max = 523 kHz while the opening and closing frequencies of the first band gap defined in (48) read, respectively, 1 2π

ζ 1 R λ m + 2µ m ρ f = 268 kHz and θ 1 2π = 291 kHz (85)
thus defining an interval verifying conditions (84).

Figure 10 shows the dispersion plot, numerically obtained by imposing the periodic Bloch-Floquet conditions on a unit cell, discretized by finite elements, and the plot of the effective mass ρ * obtained through homogenization. For the sake of comparison, the latter is shown rotated of π/2, as reduced frequency vs effective mass density diagram. The shaded regions in Figure 10a correspond to band gaps; their relative amplitude, defined as the bandwidth divided by the central gap frequency, is also reported. Note that, in view of the Property 7. of Section 4.2, the dispersion plot in term of reduced frequency ωa/2π, is independent from the cell size a. The shaded regions in Figure 10b correspond to the intervals of negative mass density. One can observe that the first band gap in the dispersion plot (Fig. 10a) exactly corresponds to the frequency range of negative effective mass (Fig. 10b); this is not the case for the second one, as it occurs at frequencies that do not satisfy hypothesis (3) in Section 2. The range of validity of the predictions given by the homogenization approach increases as the contrast between the stiffness of the two constituent materials increases. Let us consider for instance inclusions characterized by elastic moduli of one tenth of those of PVC reported in Table 2. The corresponding dispersion plot is shown in Figure 11a, while Figure 11b reports the variation of the equivalent mass density with the frequency. In this case the first three zones of negative equivalent stiffness correspond to frequencies well below the maximum frequency f max = 523 kHz which is set by the matrix properties, hence they give a good prediction of the first three band gaps numerically obtained. However the zoomed view of the dispersion plot shows that the second band gap, between 110.6 kHz and 121.8 kHz is actually separated into two smaller bandgaps (between 110.6 kHz and 118.5 kHz and between 119.0 kHz and 121.8 kHz) by two flat modes. These latter correspond to local resonances of the inclusion endowed with a displacement field inside the inclusion like a dipole, with almost zero displacements in the matrix as shown in Figure 12. These modes have already been evidenced in [START_REF] Wang | Two-dimensional locally resonant phononic crystals with binary structures[END_REF]. Actually the presence of these modes is also predicted by the present homogenization approach, they correspond to the eigenfrequency ω p n with n = 1 and p = 2 which, for the considered example, is equal to 118.98 kHz (black line in Figure 11b). This frequency turns out to fall inside the second interval of negative effective mass as shown by Figure 13 where the curve den 2 (in blue), whose zero gives ω 2 n , is superposed to the effective mass plot. On the same figure the curve den 1 is plotted by a dashed black line: the zero gives the frequencies ω 1 n = ω 1 , ω 2 , ω 3 , .. which correspond to the band-gap opening.

Out-of-plane waves

Consider now the propagation of out-of-plane waves in the material with aluminium matrix and PVC inclusions. The first two intervals of frequencies leading to a negative effective mass, as defined in (44), are [150.6 kHz -177.5 kHz] and [345.7 kHz -360.5kHz]. Since for the constituents of Table 2, and a = 1 mm, the limit frequency is f max = 1 2π ω max = 523 kHz, both intervals are within the low frequency regime and hence give a good prediction of the first two bandgaps. This is confirmed by the comparison between the dispersion plot, numerically obtained by imposing the Bloch-Floquet conditions (Fig. 14a) and the effective mass plot obtained through homogenization (Fig. 14b).

Reduced frequency !a/2⇡ [mm/s]

Conclusions

In this paper, we have provided a detailed study of band gaps formation in 2D binary, locally resonant metamaterials. Under the hypothesis of soft inclusion in hard matrix, the problem has been tackled by two different approaches: the homogenization approach, based on the scale separation, and the Bloch-Floquet analysis. In the first case the band gaps are identified as the intervals of negative effective mass of the homogenized continuum, while in the second band gaps are identified as the frequency intervals where no real solution of the Bloch-Floquet problem exists.

The main findings, which are valid under the hypothesis of high contrast between the constituents and in the long wave length limit, can be summarized as follows:

-the two approaches are fully equivalent and provide the same band-gaps prediction; -in the case of circular inclusion the intervals of negative effective mass density are computed in closed form both for in-plane and out of plane propagation; -the opening frequency of the band gap only depends on the material and geometric properties of the inclusion; -the closing frequency of the band gap also depends on the mass ratio between inclusion and matrix and increases as this ratio increases, thus leading to a wider band gap;

-the elastic properties of the matrix indirectly intervene through the hypotheses of validity of the approach, but do not intervene directly in the band gap limits; -there exist some locally resonant modes that correspond to flat modes that can fall inside the band gaps.

All above results are illustrated through several examples on specific materials. A very good agreement between the numerical results of the analysis of dispersion with Bloch-Floquet boundary conditions and the analytical results obtained through homogenizations is obtained and the range of validity of the approach is confirmed. It is worthnoting that the two scale homogenization adopted requires that the ratio between the stiffness of the inclusions and the stiffness of the matrix is of the order of 2 , with very small. Actually numerical tests have shown that this approach leads to a very good prediction of the first bandgap, already for ≈ 0.04. The closed form expressions obtained are well suited for material optimization and design.

The present approach cannot be used to predict band gaps formation in composite materials with stiff inclusions in soft matrix, nor an extension seem feasible in this direction. On the contrary, the estension of the homogenization approach for band gaps prediction in ternary LRM (e.g. materials with hard inclusions with soft coating embedded in a stiff matrix) seems possible and is currently under study.

Fig. 1

 1 Fig. 1 The 2D periodic material: (a) geometry, (b) 2D lattice, (c) unit cell with basis vectors a i , (d) reciprocal lattice with basis vectors A i and the First Brillouin Zone

Fig. 3

 3 Fig.3Coefficients ζn of the first three eigenfrequencies of the inclusion associated with the eigenmode number p = 1 vs Poisson's ratio ν f of the inclusion.

Fig. 4

 4 Fig. 4 Graphs of (a) J 2 (ζ 1 ) and J 2 (c ζ 1 ) and (b) J 2 (ζ 3 ) and J 2 (c ζ 3 ) as functions of the Poisson coefficient ν f . At ν f = ν * ≈ 0.2035 condition (29) if fulfilled with n = 3, hence the third eigenmodes are orthogonal to the translations for ν f = ν * .

1 Fig. 5

 15 Fig.5Normalized effective mass vs. normalized frequency for ν f = 0.4. The intervals of negative effective mass give the first four band gaps. The mass ratio between matrix and inclusion is 1.41 and it corresponds to the material considered in Section 5.

23 Fig. 6

 236 Fig. 6 Normalized effective mass vs. normalized frequency for ν f = ν * (left) and ν f = 0.23 (right).

  4048, z 2 = 5.5201, z 3 = 8.6537, z 4 = 11.7915, • • • .

Fig. 9

 9 Fig. 9 2D periodic material, PVC inclusions in aluminium matrix: (a) material lattice and elementary cell, (b) reciprocal lattice and IBZ (shaded area). Constituent λ [MPa] µ [MPa] ρ [Kg/m 3 ] matrix (Al) 57375. 27000. 2500. fiber (PVC) 160. 40. 1500.

⇥10 5 Fig. 10 Fig. 11 2D

 51011 Fig. 10 2D periodic material, PVC inclusions in aluminium matrix, in-plane propagation: (a) dispersion plot with shaded bandgaps, (b) effective mass density vs reduced frequency, shaded areas correspond to negative effective mass.

Fig. 12

 12 Fig. 12 2D periodic material, in-plane propagation, very soft inclusions in Aluminium matrix. Representation of the 6th flat eigenmode: (a) contour of the displacements modulus, (b) direction and amplitude of displacements vectors.

Fig. 13

 13 Fig. 13 Graphs of ρ * /ρst in orange, den 1 in black, dashed and den 2 in blue for ν f = 0.4 and mass ratio between matrix and inclusion equal to 1.41.

4

  E↵ective density ⇢ ⇤⇤ [10 3 kg/m 3 ]

Fig. 14

 14 Fig. 14 2D periodic material, PVC inclusions in aluminium matrix, out-of-plane propagation: (a) dispersion plot with shaded bandgaps, (b) effective mass density vs reduced frequency, shaded areas correspond to negative effective mass.

Table 1 and

 1 Figure 3 report the coefficients ζ n of the first three eigenfrequencies for different Poisson's coefficients.

	ν f	ζ 1	ζ 2	ζ 3
	0	2.07516	3.77206 5.31666
	0.2	1.97804	3.28481 5.14909
	0.4	1.74608	2.2812	3.48963
	0.49 0.715047 1.16922 1.59712

Table 1

 1 Coefficients ζn of the first three eigenfrequencies of the inclusion associated with the eigenmode number p = 1 for several values of ν f .

Table 2

 2 Material properties

Acknowledgements -The first Author wishes to thank the Colleagues of École polytechnique for kind hospitality and the Laboratoire de Mécanique des Solides, for financial support.

Therefore the non vanishing stress components are σ 31 = σ 13 and σ 32 = σ 23 . Hence we will consider that the stress field is the vector field σ (x 1 , x 2 ) = σ 31 (x 1 , x 2 ), σ 32 (x 1 , x 2 ) which is defined on S and related to u by

The equation of motion [START_REF] Åberg | The usage of standard finite element codes for computation of dispersion relations in materials with periodic microstructure[END_REF] in this case becomes

The expansions of u and σ with respect to read

where the fields u i = u i (x, y) and σ i = σ i (x, y) are defined on S × R 2 and are Y -periodic with respect to y. Inserting these expansions into (51)-( 52) and following the same procedure as for the in-plane case, one first gets grad y u 0 = 0 in Y m and hence

Then equation of motion (52) gives at order -1

whereas the constitutive equation ( 51) gives at order 0

One deduces by linearity that u 1 can be written in Y m as

where χ i , for i = 1, 2, are the out-of-plane displacement fields solutions of the elementary cell problems:

Therefore the stress field σ 0 can read as


M Yang

O Couturaud

W Desrat

C Consejo

D Kazazis

R Yakimova

M Syväjärvi

M Goiran

J Béard

P Frings

M Pierre

A Cresti

W Escoffier

B Jouault

Puddle-induced resistance oscillations in the breakdown of the graphene quantum Hall effect

We report on the stability of the quantum Hall plateau in wide Hall bars made from a chemically gated graphene film grown on SiC. The ν = 2 quantized plateau appears from fields B 5 T and persists up to B 80 T. At high current density, in the breakdown regime, the longitudinal resistance oscillates with a 1/B periodicity and an anomalous phase, which we relate to the presence of additional electron reservoirs. The high field experimental data suggest that these reservoirs induce a continuous increase of the carrier density up to the highest available magnetic field, thus enlarging the quantum plateaus. These in-plane inhomogeneities, in the form of high carrier density graphene pockets, modulate the quantum Hall effect breakdown and decrease the breakdown current.

Graphene [1,2] shows a unique half-integer quantum Hall effect (QHE) with conductivity plateaus σ xy = 4(m+1/2)e 2 /h where the factor 4 stands for the spin and valley degeneracies and m = 0, ±1, ±2... [3]. The peculiar dispersion E m ±420 |m|B[T] K of the Landau levels (LLs) in graphene induces large energy gaps at low LL indices and allows to explore exotic transport phenomena even at relatively high temperature [4]. The substrate graphene is deposited on plays a central role in determining the features of the observed QHE. In low-mobility graphene on SiO 2 , the presence of electron-hole puddles at the charge neutrality point (CNP) prevents any divergence of the longitudinal resistance at filling factor ν < 2, whereas the Hall resistivity fluctuates around zero due to charge compensation [5]. On the other hand, in high mobility graphene deposited on boron nitride flakes, a complete degeneracy lifting of the Landau levels can be observed and the corresponding spin-valley textures have been identified [6,7]. In this Letter, we consider graphene deposited on top of a SiC substrate (G/SiC). In this system, the quantum Hall plateau at h/2e 2 is exceptionally robust with respect to magnetic field [8]. Moreover, G/SiC shows metrological quantum Hall quantization, with relative accuracies of the quantized resistance better than 10 -9 [9], even at lower magnetic fields and higher temperatures than GaAs-based quantum Hall resistance standards [10]. In this context, it is important to unveil the role of charge reservoirs in the stabilization of the first quantum Hall plateau [11,12] and to identify the mechanisms governing the breakdown of the QHE [8,13,14].

In the following, we report on the results obtained from graphene samples grown by Graphensic AB company on the Si-face of a 4H-SiC substrate. The details concerning the growth conditions can be found in Ref. 15. The samples are tailored into Hall bar geometry with a width W = 100 µm and a total length of 420 µm. The as-grown carrier density of G/SiC is of the order of 1 × 10 13 cm -2 and is reduced to about 5 × 10 11 cm -2 using a polymer gate (photo-chemical doping) [15]. The mobility at T = 4 K is about 7,000 cm 2 V -1 s -1 . The Hall bars are electrically connected using wedge bonding and placed in a low temperature and high (pulsed) magnetic field setup for magneto-transport measurements.

Figure 1 shows the longitudinal and Hall resistances R xx and R xy of one sample, named S1, measured as a function of the magnetic field at a temperature T = 4.2 K. In this Letter, the longitudinal resistance is normalized to a square. For B ≥ 7 T, an extremely wide quantum Hall plateau at R xy = h/2e 2 is observed up to the highest available magnetic field (78 T). As expected, the quantum Hall plateau coincides with a complete vanishing of R xx . A small mixing between R xx and R xy , due to the anisotropy of the conductivity induced by the SiC steps [16,17], is noticed as a small additional bump in R xy at B 7 T. This mixing can been removed by (anti)symmetrizing the measured resistances recorded in the two opposite directions of the magnetic field:

Rxx(xy) = [R xx(xy) (B) ± R xx(xy) (-B)]/2.
Figure 2(a) shows the symmetrized magnetoresistances of sample S1 at different temperatures. The longitudinal resistance Rxx has a minimum around B min 18 T when T increases above 30 K. We assume that, as in two-dimensional semiconductor systems, this minimum occurs when the Fermi energy is at mid-gap between two adjacent Landau levels. This implies that the filling factor ν equals 2 at B = B min and the carrier density increases from 3.4 × 10 11 cm -2 at B = 0 T up to 8.5 × 10 11

cm -2 at B = B min .
This phenomenon is in qualitative agreement with the model proposed in Ref. [18], based on a charge transfer taking place between graphene, the chemical gate and the interface states between SiC and graphene. In particular, it was shown [8,12] that the carrier concentration in similar systems can increase by 200 % from B = 0 T to B = 16 T. However, at higher magnetic fields, when only the m = 0 LL located at the CNP is populated, the model predicts the electron concentration saturation, which is hardly compatible with our experimental findings. Indeed, by assuming a constant carrier density n = 8.5 × 10 11 cm -2 above B min , we would have ν = 0.45 at B = 78 T. Therefore, the large plateau observed in Fig. 1 questions the charge transfer model and deserves further attention.

At low temperatures, the electron diffusion in the quantum Hall regime is often well described by the variable range hopping (VRH) theory in the soft Coulomb gap regime [START_REF] Shklovskii | Electronic properties of doped semiconductors[END_REF], where the conductivity reads:

σ xx = σ 0 T exp(- T 0 T ). (1)
Here, σ 0 is a parameter, T 0 is the hopping temperature, which is related to the localization length ξ by the relation k B T 0 = Ce 2 /4π r 0 ξ, 0 is the vacuum dielectric constant, r 7 is the relative dielectric constant averaged between the dielectric constants of the resist and SiC, and C 6.2 is a constant [START_REF] Lien | [END_REF]21]. VRH has already been observed in various graphene samples below T = 100 K [14,[22][23][24]. The analysis of sample S1, shown in Fig. 2(b), uses only the data in the range T ≥ 30 K, where measurement errors are negligible, and T ≤ 120 K, where conduction by activation to the extended states can be neglected [START_REF]including both VRH and activation is presented in the supplemental material[END_REF]. The data are satisfactorily fitted by Eq. 1 over two decades. The localization length ξ, the σ 0 coefficient and the hopping temperature T 0 are extracted from the VRH fits, see Fig. 2(c,d). At low magnetic fields B 20 T, ξ is comparable to the magnetic length l B = /eB. At higher B, ξ increases only slightly. In the QHE regime, the localization length is expected to vary according to ξ ∝ ∆E -γ , where ∆E is the mobility gap between the nearest delocalized states and the Fermi energy and γ 2.3 is the quantum percolation exponent [START_REF] Huckestein | [END_REF]. Assuming a constant charge density n and a density of localized states independent of energy, the localization length ξ th B would vary as 1/ν 2.3 ∼ B 2.3 , yielding ξ th B=60T /ξ th B=20T 12. By contrast, the experimental variation of ξ in the 20-60 T range is about ξ B=60T /ξ B=20T 2, see Fig. 2(c). This gives n B=60T 2.2 × n B=20T , thus indicating a carrier density increase even for B > B min .

Figure 3 shows R xx of sample S1 for current I in the range 50-400 µA. At I = 100 µA, R xx > 3 Ω over the whole magnetic field range, yielding a longitudinal voltage (V xx > 300 µV) too large to be compatible with metrological measurements. This defines a critical breakdown current I c 1 A/m, one order of magnitude lower than the highest reported values for graphene on SiC [8]. Moreover, R xx (B) does not have any more a unique minimum associated to ν = 2 but shows oscillations periodic in 1/B. The periodicity is evidenced in the Landau plot reported in Fig. 3(b). Assuming a degeneracy of spin and valley, the period ∆(1/B) of these oscillations corresponds to a constant high carrier density n HD = 4(e/h)/∆(1/B) 1.15 ×10 13 cm -2 .

Similar 1/B-periodic oscillations have been observed in other samples. The periodicity corresponds to concentrations in the range 4 -12 × 10 12 cm -2 . The period and amplitude of the oscillations change each time a thermal cycle takes place but the period does not change when the current is increased. In all of the samples, resist cracks appeared after several weeks of measurements at low tem-peratures and several temperature cycles. We postulate that even before the appearance of cracks, the chemical gating is not homogeneous and microscopic detachments of the resist lead to disconnected graphene puddles with native High carrier Density (HD), whereas the rest of the Hall bar remains covered by the resist and has a Low electron Density (LD). Additionally, graphene detachment from the substrate induced by the SiC steps [27], bilayer patches [28] or residual contamination [29] may also induce inhomogeneous carrier density.

For most samples, the I/V xx characteristic is ohmic, the 1/B oscillations appear independently of the injected current and the transport is dissipative. This is not the case for sample S1, which shows a more remarkable behavior. There, the oscillations appear only above a critical current I c 10 µA in the breakdown regime of the quantum Hall effect. At lower current, S1 keeps good metrological properties, as evidenced by additional precision measurements performed in constant magnetic field. At B = 10 T, T = 1.8 K and a DC current I = 10 µA the deviation of the Hall resistance from the expected value h/2e 2 is (∆R H -h/2e 2 )/(h/2e 2 ) = (-0.4 ± 10) × 10 -6 , while the longitudinal resistance is (10 ± 17) mΩ. This precision and this uncertainty correspond to the limit of the setup. This analysis rules out parallel conduction.

Another key feature of such 1/B-periodic oscillations for all samples is that their Berry phase is zero, as evidenced by the Landau plot of Fig. 3(b) for sample S1.

We can take advantage of these unusual features to unravel the role of the inhomogeneity in the carrier density in the stabilization of the quantum Hall plateau. We call S HD and S LD the total areas corresponding to the HD and LD regions. We label n i 0 and µ i 0 the charge density and the chemical potential in each region i ∈ {HD, LD} in the absence of magnetic field. Here µ i 0 is the chemical potential with the energy reference at the CNP of region i. These parameters are constant: n HD 0 11.5 × 10 12 cm -2 from the 1/B-period of the oscillations and n LD 0 = 10 11 cm -2 from the Hall effect at low field B 0.1 T. We label ∆µ HL = µ HD 0 -µ LD 0 = 330 meV the energy difference between the CNP points in the LD and HD regions.

According to Ref. 18, for each homogeneous region i the electron density is given by n i = -n i g + β(A -µ i ) where n i g is the electron density pinned by the electrochemical gate, A is the difference between the work function of undoped graphene and donor states and β is an effective density of states. In the highly doped region, no polymer is present and the charges come from the substrate. Therefore n HD g = 0 and β = n HD 0 /(A -µ HD 0 ). In the low doped region, the polymer acts as a gate and then n LD g = -n LD 0 + n HD 0 (A -µ LD 0 )/(A -µ HD 0 ). Once these parameters have been determined, we let the magnetic field change. The carrier density in each region evolves according to the dispersion of the graphene LLs characterized by energies E m and a Lorentzian broadening Γ. Since the regions are coupled, electrons redistribute to keep the same chemical potential µ everywhere: µ = µ LD = µ HD -∆µ HL . For each value of the magnetic field, µ is determined by solving numerically n = (1 -α) -n LD g + β(A -µ) + αβ A -µ -∆µ HL where n is the average electron concentration: n = (1 -α)n LD + αn HD and α = S HD /(S LD + S HD ). The only remaining fitting parameters are A, Γ and α. We determine the magnetic field-dependent electron density in the two regions as shown in Fig. 4(a). The charge density in LD region increases without saturation up to high magnetic field whereas it oscillates in HD region. These oscillations are 1/B-periodic and in turn induce a modulation of the filling factor with an opposite phase in the LD region. In the transport properties, this phenomenon translates into 1/B-periodic oscillations of R xx (B) with zero apparent Berry phase each time the Fermi energy approaches the LL extended states.

Each HD puddle, with a conductivity much higher than the conductivity of the LD region, behaves as an equipotential region. The direction of the electric field at the HD/LD interface is normal to the interface. In sufficiently high magnetic fields, the current forms a Hall angle close to 90 o with the electric field and cannot enter into the HD puddles which act as electrically insulating regions. This implies a dramatic reduction of the effective Hall bar width W eff over which the current flows. In the VRH theory [30], the current acts as an effective temperature T eff given by the relation:

T eff = eR H Iξ/(2k B W eff ).
( Whereas the bootstrap electron heating theory is most widely accepted to account for the breakdown of the QHE [13], here the underlying VRH mechanism alone is sufficient to explain the observed R xx (B) dependence. The localization length varies as ξ = aµ -γ , where a is chosen to have ξ B=60T = 12 nm, see Fig. 2(c) and µ is given by the calculation, see Fig. 4(b). Introducing ξ in T eff and T eff in Eq. ( 1), R xx is calculated at various currents and magnetic fields, as shown in Fig. 4(c). The 1/B periodic oscillations are indeed reproduced, with the correct experimental phase and amplitude. The best agreement with the data is found by choosing A = 0.4 eV, α = 30% and Γ = 15 meV. The large α value points towards a largely inhomogeneous sample and is consistent with the small W eff /W value. The determination of a more quantitative relation between α and W eff would require a detailed knowledge of the current distribution. The value A = 0.4 eV is widely accepted in the literature, as well as the relation γ = β 0 r /( 0 r -e 2 dβ), where γ is the density of states at the SiC/graphene interface and d is the graphene thickness. However, the condition γ > 0 implies r > 130 × d[nm]. The value r 7 yields d < 0.05 nm, suggesting i) a graphene layer closer to the donor interface states than usually assumed; ii) an enhancement of the dielectric screening by the π electrons. Using the Random Phase Approximation [3], one gets ∞ RPA 13.

To conclude, we addressed the robustness of the ν = 2 QH plateau for graphene on SiC. We observed that this plateau can surprisingly persist at least up to 78 T. Faint and rapid oscillations on top of the longitudinal resistance also appear in the breakdown regime. This effect was demonstrated to be intrinsic to graphene and to originate from non-homogeneous in-plane doping of the device. A charge transfer [31] between the SiC substrate and graphene regions of different carrier concentration is entailed by the magnetic field due to a modification of the band structure and the formation of Landau levels. These additional charge reservoirs modulate the breakdown mechanisms, lower the breakdown current and help to stabilize the quantum plateau up to very high fields.

Part of this work was performed at LNCMI under EMFL proposal TSC06-115.

FIG. 1 :

 1 FIG. 1: (a) Longitudinal and Hall magnetoresistances of sample S1 measured at T = 4 K up to B= 80 T. The inset is a sketch of the Hall bar covered by the resists.

FIG. 2 :

 2 FIG. 2: (a) Longitudinal (solid lines) and transverse (dashed lines) symmetrized resistances for sample S1, at T = 1.6 K up to T = 120 K. (b) solid lines, open symbols: T σxx vs T -1/2 from B= 20 T up to B= 60 T with a step of 10 T. Dashed lines: linear fits with the soft Coulomb gap VRH model. The data of Ref. 14 at B = 19 T are reported for comparison (red dashed line). (c) Open symbols: localization length ξ extracted from the VRH fit. The blue and red lines are guides for the eye, following the magnetic length lB and the prediction ξ ∼ 1/ν 2.3 respectively. (d) Parameters T0 and σ0 extracted from the fit.

FIG. 3 :

 3 FIG. 3: (a) Longitudinal resistance Rxx for currents higher than the breakdown current Ic 10 µA. 1/B-periodic oscillations are visible. (b) Landau plot for the minima and maxima of these oscillations. The period corresponds to a carrier concentration nHD ≈ 1.15 × 10 13 cm -2 . The Rxxmaxima plotted as a function of their index N go to N = 1/2 when BN → ∞, which is usually a signature of a zero Berry Phase. (c) The observed oscillation damping by the temperature does not agree with the usual Shubnikov-de Haas theory. The oscillation period may change from one cool-down (solid line) to another (dashed lines).

  ) Experimentally, T eff is determined by matching the conductivities σ xx (I) = σ xx (T eff ). At B < 20 T, where the 1/B-periodic oscillations are small, T eff appears to be roughly linear versus I, in agreement with Eq. 2. The relation σ xx (I) = σ xx (T eff ) yields T eff 0.05-0.1 ×(I[µA]) [K] which corresponds to an effective Hall bar width W eff 15 µm considerably smaller than the nominal Hall bar width (100 µm). A similar discrepancy between W eff and W was observed in Ref. 14, 21, again suggesting a very inhomogeneous current distribution in the samples.

FIG. 4 :

 4 FIG. 4: (a) Charge density in LD and HD regions as a function of the magnetic field. The proportion of highly doped regions α is fixed to 30% of the total surface (thick blue and red lines: nHD and nLD respectively), A = 0.4 eV and Γ = 15 meV. The thin dashed lines represent the LL carrier density in LD region. The inset is a sketch of the inhomogneous Hall bar embedding several highly doped regions, a deflected current line is sketched by open arrows. (b) Chemical potential µ as a function of B for α= 30% (red line) and α = 0 (black line). The thin blue dashed lines represent the positions of the Landau Levels in the HD puddles. (c) Longitudinal resistance Rxx as a function of B, at various currents, calculated within the model described in the text. The curves from bottom to top correspond to currents of 100, 250 and 400 µA respectively. The inset is a Landau plot (blue circles: maxima indexed with integers, red circles: minima indexed with half-integers).