A Avila
email: aavila@ufro.cl

Georges Griso
email: georges.griso@wanadoo.fr

Bernadette Miara
email: miarab@esiee.fr

E Rohan
email: rohan@kme.zcu.cz

A Ávila

Multi-scale modelling of elastic waves Theoretical justification and numerical simulation of band gaps

Keywords: elastic waves, composite materials, strong heterogeneities, homogenization, negative mass density, band gaps

Multi-scale modelling of elastic waves. Theoretical justification and numerical simulation of band gaps

Introduction

After the huge impact due to the development of photonic crystals [START_REF] Pendry | [END_REF]14], the development of phononic crystals has received growing interest in recent years. These artificial crystals, which mimic a crystalline atomic lattice, are structured materials formed of periodic microstructures. In the case of phononic crystals considered by Vasseur and al. [START_REF] Vasseur | Experimental evidence for the existence of absolute acoustic band gaps in two-dimensional periodic composite media[END_REF], the media is a two-dimensional binary solid-solid composite made of elastic arrays of Duralumin cylindrical inclusions embedded in a resin epoxy matrix. For this structure, measured transmission exhibit absolute acoustic band gaps. A band gap is a range of frequency in which elastic or acoustic waves cannot propagate; it is surrounded, above and below, by propagating states.

From a technological point of view the main interest of these composites is to help reduce the noise level, they are also good candidates for the design of elastic or acoustic waveguides or filters.

From a mathematical point of view, the homogenization approach (which consists in replacing a composite by a limit homogeneous material) is relevant for the modelling of such periodic structures. Let us note that the property for a periodic structure to always present band gaps (Floquet-Bloch's theory), disappears in the case of an homogeneous material. However, the asymptotic analysis shows that in the case of photonic crystals the limit "homogeneous permeability" is negative for certain wavelengths and hence band gaps appear [START_REF] Bouchitté | Homogenization near resonances and artificial magnetism from dielectrics[END_REF]. In the case of phononic crystals our study aims at justifying the existence of band gaps for certain wavelengths; this result is a consequence of the non positivity of the limit "homogeneous mass density". The homogenization method we use to obtain this result relies on the unfolding method [START_REF] Cioranescu | Periodic unfolding and homogenization[END_REF] that combines the dilatation technique with ideas from finite element approximations.

From a numerical point of view some computational works have been developed to optimize the shape of the inclusions [START_REF] Jensen | On maximal eigenfrequency separation in twomaterial structures: The 1D and 2D scalar cases[END_REF][START_REF] Sigmund | Systematic design of phononic band-gap materials and structures by topology optimization[END_REF] with a classical approach (different from ours which is based on the micro-macro study). For the time being, we present in this paper the numerical simulation of the band gaps, with an emphasis on weak band gaps (propagation in certain directions only) and strong band gaps (no propagation in any directions). Further studies on the sensitivity analysis of these forbidden bands are under preparation [START_REF] Rohan | Microstructural sensitivity of acoustic band gaps[END_REF][START_REF] Miara | Shape optimization of phononic materials[END_REF].

To present the problem under study let us start with the description of the geometry of the composite whose reference configuration Ω of the elastic body is supposed to be stress-free. The bounded domain Ω ⊂ R 3 with micro structures of size ε > 0, is split into a domain Ω ε 1 occupied by the matrix made of material 1, and a domain Ω ε 2 (with Lipschitzcontinuous boundary denoted by ∂Ω) which contains periodically distributed inclusions made of material 2, hence Ω = Ω ε 1 ∪ Ω ε 2 with Ω ε 1 ∩ Ω ε 2 = ∅. We note that the whole domain Ω is independent of ε, whereas the domains occupied by the matrix and the inclusions are both ε-dependent. As for the study of three-dimensional periodic structure, let us introduce the reference cell Y = [0, 1[ 3 with its elementary inclusion Y 2 , Y 2 ⊂ Y, Y 1 = Y \Y 2 . Therefore, material 2 occupies the domain Ω ε 2 obtained by ε-periodicity and material 1 occupies the remaining domain, Ω ε 1 = Ω\Ω ε 2 :

Ω ε 2 = k∈K ε ε(Y 2 + k), K ε = {k ∈ Z 3 , ε(Y 2 + k) ⊂ Ω}.
Let us denote by u ε (ω) the static elastic displacement field that the body undergoes at a fixed wavelength ω (the displacement is indexed by ε since, obviously, it depends on the microstructure size ε). Our paper deals with the convergence of the sequence {u ε (ω)} ε when ε goes to zero. In section 2 we recall the propagation equations for elastic waves solved by u ε (ω) for positive values of ε. In section 3 we introduce the unfolding operator and give its essential properties so that we are in a position, in section 4, to establish the main convergence theorem which gives the propagation equations solved by the limit displacement field u(ω). In section 5 we discuss the possibility of a "negative" mass density and its consequence for the existence of forbidden propagation bands. All the proofs of the existence and convergence theorems are given in section 6. Finally, in section 7, numerical simulations illustrate the influence of the change of some parameters of the micro structures (such as shape of the inclusions, average mass density of the composite, fill-in coefficient).

Propagation of elastic waves

In this section we recall the equations of the elastic waves propagation in the composite material described previously and next we give the equilibrium equations in the static case with fixed wavelength. Let T > 0; under the action of applied forces 1 3 which is solution to the evolution problem:

F = (F m ) : Ω × (0, T ) → R 3 the body undergoes an elastic displacement field U ε = (U ε m ) : Ω × (0, T ) → R
         r ε (x) ∂ 2 ∂t 2 U ε m (x, t) - ∂ ∂x n (c ε mnkl (x)e kl (U ε (x, t))) = F m in Ω × (0, T), U ε (x, t) = 0 on ∂Ω × (0, T), U ε (x, 0) = U 0 (x), ∂ ∂t U ε (x, 0) = U 1 (x) in Ω, with initial conditions U 0 : Ω → R 3 , U 1 : Ω → R 3
and where

e kl (V (x, t)) = 1 2 (∂ k V l (x, t) + ∂ l V k (x, t)),
is the linearized deformation tensor. The mass density r ε : Ω -→ R and the elasticity tensor of the structure c ε = (c ε mnkl ), with c ε mnkl : Ω -→ R possess the classical properties of any elastic body: -There exists two positive constants ρ ε -, ρ ε + such that:

ρ ε -≤ r ε (x) ≤ ρ ε + for all x ∈ Ω.
-The elasticity tensor is symmetric and coercive, i.e., :

c ε klmn = c ε mnkl = c ε nmlk ,
and there exists α ε > 0, β ε > 0 such that, for any symmetric matrix (X mn ), we have

α ε X mn X mn ≤ c ε mnkl (x)X mn X kl ≤ β ε X mn X mn for all x ∈ Ω.
For any fixed ε > 0, and standard assumptions on the regularity of the data

r ε ∈ L 2 (Ω), c ε mnkl ∈ L 2 (Ω), F ∈ L 2 (Ω × (0, T )), U 0 ∈ H 1 0 (Ω), U 1 ∈ L 2 (Ω),
the associated variational problem has a unique weak solution

U ε ∈ C(0, T ; H 1 0 (Ω)), ∂ ∂t U ε ∈ C(0, T ; L 2 (Ω)).
Let us now consider an incident wave period ω fixed and independent of ε, and the periodic solution U ε (x, t) = u ε (x, ω)e iωt associated to periodic applied forces F (x, t) = f (x)e iωt and compatible initial conditions; in the sequel we denote by u ε (x), instead of u ε (x, ω), the amplitude of the elastic wave.

Hence, for any fixed ε > 0, the elastic field u ε : Ω → C 3 is given by the stationary problem,

ω 2 r ε (x)u ε m (x) + ∂ n (c ε mnkl (x)e kl (u ε (x))) = -f m in Ω, u ε (x) = 0 on ∂Ω.
Because of the linearity of the problem, we consider in the sequel only real-valued displacement field u ε : Ω → R 3 . According to Fredholm alternative, for each fixed value of ω different from the resonance values (square root of the eigenvalues of the elasticity problem) which depend upon ε, the variational problem:

ω 2 Ω r ε u ε • Φ - Ω c ε mnkl e kl (u ε )e mn (Φ) = - Ω f • Φ ∀Φ ∈ H 1 0 (Ω). (1)
(where

u • v = u k v k ) has a unique solution u ε ∈ H 1 0 (Ω).
The rest of the paper aims, first of all, at showing that the sequence {u ε } ε of solutions to (1) converges (in a certain sense that will be made more precise later) to the solution of an homogeneous problem, and next at interpreting the theoretical and numerical properties of this limit solution.

3 Unfolding operator and heterogeneities

Definitions and basic properties of the unfolding operator

In conjunction with the elementary cell Y , there exists, for all z ∈ R 3 , a unique decomposition z = [z] + {z} with an integer part [z] and a remaining part {z} such that {z} = z -[z] ∈ Y and equivalently, we have the unique decomposition z = ε[ z ε ] + ε{ z ε }. We introduce the unfolding operator T ε related to the study of periodic structures [START_REF] Cioranescu | Periodic unfolding and homogenization[END_REF][START_REF] Damlamian | An elementary introduction to Periodic Unfolding[END_REF]. For all v ∈ L 2 (Ω) extended by 0 outside Ω :

T ε : v ∈ L 2 (Ω) -→ T ε (v)(x, y) = v(ε[ x ε ] + εy), x ∈ Ω, y ∈ Y.
The main properties of this operator are given below (complete proofs can be found in [START_REF] Cioranescu | Periodic unfolding and homogenization[END_REF]).

• If {v ε } ε is uniformly bounded in L 2 (Ω) then, there exists v ∈ L 2 (Ω × Y ) such that, up to subsequence still denoted with the same indices, we have the convergence:

T ε (v ε ) v weakly in L 2 (Ω × Y ).
• If the sequence {v ε } ε is uniformly bounded in L 2 (Ω) and the sequence {ε∇v ε } ε is uniformly bounded in L 2 (Ω; R 3 ) then, there exists a limit field v ∈ L 2 (Ω; H 1 per (Y )) such that, up to subsequences still denoted with the same indices, we have the convergence:

T ε (v ε ) v weakly in L 2 (Ω × Y ), T ε (ε∇ x v ε ) ∇ y v weakly in L 2 (Ω × Y ; R 3 ),
where

∇ x v = ∂v(x, y) ∂x i 1≤i≤3 and ∇ y v = ∂v(x, y) ∂y i 1≤i≤3 .
• If {v ε } ε is uniformly bounded in H 1 (Ω) then, there exists a limit field v ∈ H 1 (Ω) and a corrector v ∈ L 2 (Ω, H 1 per (Y )) such that, up to subsequences still denoted with the same indices, we have the convergence:

   v ε v weakly in H 1 (Ω), T ε (v ε ) v weakly in L 2 (Ω × Y ), T ε (∇ x v ε ) ∇ x v + ∇ y v weakly in L 2 (Ω × Y ; R 3 ).
Periodic problems have also been studied by the two-scale method [START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Nguetseng | A general convergence result for a functioanl related to the theory of homogenization[END_REF]

Heterogeneities

We are now in a position to state the dependence of the materials characteristics in terms of ε. When these characteristics are not scaled by ε the limit homogeneous model does not exhibit band gaps. Since we are interested by the modelling of the bad gap structure we assume in the sequel that there exists new functions r 1 , r 2 , c 1 , c 2 independent of the size ε of the micro structures such that:

       (In the matrix) r ε (x) = r 1 ({ x ε }), c ε mnkl (x) = c 1,mnkl ({ x ε }) x ∈ Ω ε 1 , (In the inclusions) r ε (x) = r 2 ({ x ε }), c ε mnkl (x) = ε 2 c 2,mnkl ({ x ε }) x ∈ Ω ε 2 .
in other words, functions r 1 , r 2 , c 1 , c 2 are defined in the elementary cell Y by the relations,

   T ε (r ε )(x, y) = r 1 (y), T ε (c ε mnkl )(x, y) = c 1,mnkl (y), x ∈ Ω ε 1 , y ∈ Y 1 , T ε (r ε )(x, y) = r 2 (y), T ε (c ε mnkl )(x, y) = ε 2 c 2,mnkl (y), x ∈ Ω ε 2 , y ∈ Y 2 .
The scaling ε 2 , which appears in front of the elasticity tensor, is the expression of the strong heterogeneity that exists between the elastic properties of the matrix and that of the inclusions ones. Different kinds of scalings are possible, however this scaling is the only one that gives rise to a limit model with significant physical meanings (in this direction, see [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] and [START_REF] Arbogast | Derivation of the double porosity model of single phase flow via homogenization theory[END_REF] where an example for which this type of assumption is used to solve a double porosity problem). This scaling and the continuity and coercivity constants introduced in Section 3.1 now read: There exists positive constants ρ -, ρ + , α, β independent of ε such that:

ρ -≤ r 1 (y) ≤ ρ + for all y ∈ Y 1 , ρ -≤ r 2 (y) ≤ ρ + for all y ∈ Y 2 ,
and such that, for all symmetric matrix (X mn ), we have:

αX mn X mn ≤ c 1,mnkl (y)X mn X kl ≤ βX mn X mn for all y ∈ Y 1 ,
and the same property holds for the elasticity tensor c 2 , αX mn X mn ≤ c 2,mnkl (y)X mn X kl ≤ βX mn X mn for all y ∈ Y 2 .

The limit problem

We recall that the solution u ε to (1) exists for all frequencies ω different from resonance values (which depend upon ε). Therefore, before showing that when the size ε of the micro structures in the composite goes to zero the sequence {u ε } ε converges to the solution of limit problem, we show that there exist an admissible set of frequencies W (independent of ε) such that {u ε } ε exists for all frequencies in W for ε small enough. However this is not possible for all structures, hence we have to restrict our study to a class of structures which allow limit wave propagation.

Assumption on the data. We assume that the density r 2 , the elastic characteristics c 2 and the geometry of the elementary inclusion Y 2 and the density r 1 assure the existence of a non empty open set of frequencies W ⊂ R + (this will be made more precise in Step 7 of Section 6), that allows us to prove the existence of a solution to (1) for all frequency ω ∈ W and next its convergence to the solution of a limit problem.

Existence Theorem 1. For all ω ∈ W , there exists a positive value ε 0 (ω) such that, for all ε ∈]0, ε 0 (ω)] problem (1) has a unique solution bounded in L 2 (Ω).

The proof, for the sake of clarity, is postponed to Step 7 of Section 6.

Convergence Theorem 2.

For all values of the incident wave ω ∈ W there exists two limit vector fields

u 1 ∈ H 1 0 (Ω), u 2 ∈ L 2 (Ω; H 1 0 (Y 2 )) such that: (i) The sequence T ε (u ε ) strongly converges to u = u 1 + u 2 in L 2 (Ω × Y ). (ii) The limit u 1 ∈ H 1 0 (Ω)
is the unique solution to the variational problem:

ω 2 Ω A * (ω)u 1 (x) • Φ(x)dx - Ω c * mnkl e kl (u 1 (x))e mn (Φ(x))dx = - Ω B * (ω)f (x) • Φ(x)dx ∀Φ ∈ H 1 0 (Ω), (2)
where the homogenized elasticity tensor c * , the homogenized generalized mass density matrix A * (ω) and the homogenized matrix B * (ω) are given in ( 6) and ( 9) below in section 6.

(iii) The limit displacement field

u 2 ∈ L 2 (Ω; H 1 0 (Y 2 )
) is the unique solution to the variational problem:

ω 2 Y 2 r 2 (y)u 2 (x, y) • Ψ(y)dy - Y 2 c 2,mnkl (y)e kl,y (u 2 (x, y))e mn,y (Ψ(y))dy = -ω 2 u 1 (x) • Y 2 r 2 (y)Ψ(y)dy -f (x) • Y 2 Ψ(y)dy ∀Ψ ∈ H 1 0 (Y 2 ). (3)
The existence of a unique solution u 1 to problem (2) is proved in Step 4, Lemma 3 of Section 6, the existence of a unique solution u 2 to problem (3) is proved in Step 3 of Section 6.

Problem (2) has the same form as problem [START_REF] Allaire | Homogenization and two-scale convergence[END_REF]. However it is important to notice that the mass density scalar r ε has been replaced by the homogenized tensor A * (ω) which depends upon the wavelength considered ω, and that the elasticity tensor c * = (c * mnkl ) is now homogeneous and independent of ω.

Negative mass density and band-gaps

In view of the limit model given by the solution of problem ( 2), the propagation of waves in the homogenized structure that occupies the domain Ω depends on the structure of the "mass density matrix" A * (ω). Therefore we begin by the examination of its properties in order to prove the existence of the elastic band-gaps, and next we introduce the notion of the so-called weak or strong band-gaps specially suitable for applications.

Properties of the homogenized "mass density matrix" A * (ω)

Let us first give the expression of A * (that will be justified in Step 3 of section 6),

A * (ω) = j∈J A * ,j (ω) + r * I, r * = Y 1 r 1 (y)dy + Y 2 r 2 (y)dy.
The elements of each matrix A * ,j = (A * ,j pq ) are given by:

A * ,j pq (ω) = -ω 2 ω 2 -λ j Y 2 r 2 (y)ϕ j p (y) Y 2 r 2 (y)ϕ j q (y),
where {ϕ j , λ j } j≥1 are the eigenelements associated to the elasticity operator (r 2 , c 2 )

posed in the domain Y 2 , see equation ( 7) below, and

J = {j ≥ 1, Y 2 r 2 (y)ϕ j (y)dy = 0}.
Hence, matrix A * (ω) is symmetric, its eigenvalues are real; however when this matrix is not positive definite (i.e., when it has at least one negative eigenvalue), the limit problem may have evanescent solutions which means that there is no wave propagation in certain directions corresponding to any linear combination of eigenvectors associated to the negative eigenvalues.

The sign of eigenvalues of matrix A * (ω) will now be investigated in each elementary interval ]

√ λ j , √ λ j+1 [. The smallest eigenvalue µ ∇ (ω) is given by the infimum of the function z → A * (ω)z • z with z = (z p ) ∈ R 3 and ||z|| = 1, µ ∇ (ω) = inf ||z||=1 j∈J -ω 2 ω 2 -λ j 3 p=1 z p Y 2 r 2 (y)ϕ j p (y) 2 + r * .
For any j ∈ J, the function ω → µ ∇ (ω) is strictly increasing in the interval ] √ λ j , √ λ j+1 [. We can distinguish the following two cases: either λ j is of order of multiplicity equal to one, or it is of order of multiplicity greater than one (this situation happens, for example, when the inclusions present geometrical symmetries).

• In the first case, matrix A * ,j (ω) is of rank one, hence the infimum of A * ,j (ω)z

• z vanishes. Since in each interval ] √ λ j , √ λ j+1 [ all quantities A * ,k (ω)z • z are bounded for all k = j, we infer that the smallest eigenvalue increases from -∞ to a finite value µ j ∇ = µ ∇ ( √ λ j+1 ). If this value µ j ∇ is non negative, there exists a value ω j ∇ ∈] √ λ j , √ λ j+1 [ such that µ ∇ (ω j ∇ ) = 0, hence ] √ λ j , ω j
∇ [ is a band gap in the sense that it may happen that some waves are not progressive. By contrast, in the interval ]ω j ∇ ,

√ λ j+1 [ all waves propagate. If µ j ∇ is strictly negative then the whole interval ] √ λ j , √ λ j+1
[ is a forbidden band. • In the second case, λ j is of order of multiplicity M > 1 and (ϕ j,m ) m=1,M are the associated eigenvectors. Matrix A * ,j (ω) = (A * ,j pq (ω)) is given by:

A * ,j pq (ω) = -ω 2 ω 2 -λ j M m=1 Y 2 r 2 (y)ϕ j,m p (y) Y 2
r 2 (y)ϕ j,m q (y).

If matrix A * ,j (ω) is of full rank, (contrary to the previous case), µ ∇ (ω) increases from -∞ to +∞, and there always exists a value

ω j ∇ ∈] √ λ j , √ λ j+1 [ such that µ ∇ (ω j ∇ ) = 0.
In other words there always exists a band gap ] √ λ j , ω j ∇ [. In the case where matrix A * ,j (ω) is not of full rank the behavior of its eigenvalue is the same as in the previous case, i.e., bounded at ω =

√ λ j+1 . To sum up, in each interval ] √ λ j , √ λ j+1 [, the smallest eigenvalue of A * (ω) takes negative values in each intervals ] √ λ j , ω j ∇ [, with possibly ω j ∇ = √ λ j+1
; in this last case, the band gap extends to the whole interval ]

√ λ j , √ λ j+1 [.
By the same way, we establish that the largest eigenvalue µ ∆ (ω) given by the supremum of the same function z → A * (ω)z • z, z ∈ R 3 , ||z|| = 1, behaves in the following way:

• Either, it increases from a finite value µ j ∆ to +∞, • Or, it increases from -∞ to +∞. Hence, as before, in each interval ] √ λ j , √ λ j+1 [, it may exist a value ω j ∆ such that the largest eigenvalue of A * (ω) takes negative values in each interval ] √ λ j , ω j ∆ [ with possibly ω j ∆ = √ λ j ; in the latter case, the largest eigenvalue is positive in the whole domain.

Strong and weak band gaps

With the notations introduced in the previous section, we call strong band gap the interval ] √ λ j , ω j ∆ [ and weak band gap the interval ]ω j ∆ , ω j ∇ [, with possibly

ω j ∆ = √ λ j or ω j ∇ = √ λ j+1
. This means that, in a strong band gap, matrix A * is negative definite and there is no propagation in any direction; in a weak band gap, matrix A * is neither positive nor negative, there is propagation in at least one direction corresponding to its positive eigenvalue. Finally, in the interval ]ω j ∇ , [. Let us insist on the introduction of these definitions which are justified by their importance on the applications (such as for example noise suppression or reduction in one or all directions).

√ λ j+1 [ matrix A * is

Proof of the existence and convergence Theorems.

The convergence relies partly on Bouchitté and Feldbacq's results established in the Helmoltz diffusion case [START_REF] Bouchitté | Homogenization near resonances and artificial magnetism from dielectrics[END_REF]. We generalize their approach to the framework of linearized elasticity. The proof of the convergence Theorem is broken into 7 steps. In Step 1 we assume that the sequence of solutions {u ε } ε is uniformly bounded in L 2 (Ω). This yields, in Step 2, to the existence of two limit fields u 1 , u 2 , and of a corrector u which are coupled solutions to the limit problem. In Step 3 we solve this limit problem, so that the limit field u 1 is solution to a wave propagation equation and we identify u 2 and u as solution to variational problems. In Step 4, we establish the existence of a unique solution for the limit problem (2). In Step 5 we show the strong convergence of the sequence {T ε u ε } ε , and finally in Step 6 we show by contradiction, that the a priori bound is satisfied, this conclude the proof. Finally, in Step 7, we prove the existence of a unique solution for problem [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] for small enough values of ε.

Step 1. We begin by an a priori assumption. Let us suppose that the sequence {u ε } ε is bounded uniformly in ε, i.e., there exists a constant C > 0 such that for all ε > 0,

||u ε || L 2 (Ω) ≤ C.
Hence, by taking in the stationary problem (1), the test function Φ = u ε , we get:

Ω ε 1 c ε mnkl e kl (u ε )e mn (u ε ) + Ω ε 2 c ε mnkl e kl (u ε )e mn (u ε ) = ω 2 Ω r ε |u ε | 2 + Ω f • u ε ,
and the scaling ( given in Section 3.2 ) and coercivity condition (given in Section 2) on the elasticity tensor c ε yields:

α Ω ε 1 e kl (u ε )e kl (u ε ) + ε 2 Ω ε 2 e kl (u ε )e kl (u ε ) ≤ ω 2 ρ 2 Ω |u ε | 2 + Ω f • u ε .
Therefore, the a priori assumption yield the following majoration:

||e kl (u ε )|| L 2 (Ω ε 1 ) + ||εe kl (u ε )|| L 2 (Ω ε 2 ) ≤ C, 1 ≤ k, l ≤ 3, (4)
where C is a constant independent of ε.

Step 2. Convergence of the unfolded sequences

For all open set O ⊂ R 3 let us introduce the elasticity semi-norm (equivalent to the H 1 0 (O) norm):

|v| E,O = i,j ||e ij (v)|| L 2 (O) .
Since we have the inclusion Y 2 ⊂ Y 1 and the boundary of Y 2 is Lipschitz-continuous, there exists a linear and continuous extension operator P :

H 1 (Y 1 ) -→ H 1 (Y ).

Lemma 1. Extension of a bounded vector field

Let v ∈ H 1 (Ω ε 1 ), there exists an extension v ∈ H 1 (Ω) that satisfies the bound:

|v| E,Ω ≤ C|v| E,Ω ε 1 ,
with C independent of ε.

Proof. First we consider a displacement field v ∈ H 1 (Y 1 ). There exists a rigid displacement r such that:

||v -r|| H 1 (Y 1 ) ≤ C|v| E,Y 1 .
Hence we can define the extension v by:

v = (v -r) + r in Y 1 , P(v -r) + r in Y 2 ,
and obviously we get

|v| E,Y ≤ C|v| E,Y 1 .
Next we consider a displacement field v ∈ H 1 (Ω ε 1 ). With the same extension operator we can define the extension v ∈ H 1 (Ω) and get the majoration |v| E,Ω ≤ C|v| E,Ω ε 1 with C independent of ε.

Corollary of Lemma 1. Decomposition of the displacement field u ε

There exists two displacement fields u ε,1 and u ε,2 such that the solution u ε of problem (1) can be decomposed as:

u ε = u ε,1 + u ε,2 . (i) The displacement field u ε,1 ∈ H 1 0 (Ω) coincides with u ε in Ω ε 1 and satisfies the bound: |u ε,1 | E,Ω ≤ C|u ε | E,Ω ε 1 with a constant C independent of ε. (ii) The displacement field u ε,2 ∈ H 1 0 (Ω) vanishes in Ω ε 1 . (iii) Moreover, we have the bounds |u ε,1 | E,Ω ≤ C, |u ε,2 | E,Ω ≤ C ε , which yield ||u ε,1 || H 1 (Ω) ≤ C, ||u ε,2 || L 2 (Ω) ≤ C. (5)
Proof. We denote by v ε the restriction of u ε to Ω ε 1 . By Lemma 1, let u ε,1 be the extension of v ε to Ω, u ε,1 = v ε , and we get the majoration of step (i). Since u ε,1 = u ε in Ω ε 1 we get step (ii).

Step (iii) is a consequence of bounds (4).

The homogenization method presented in section 3 and majorations (4) and ( 5) yield the following convergence result.

Lemma 2.

(i) There exists two limit vector fields u 1 , u 2 and a corrector u,

u 1 ∈ H 1 0 (Ω), u 2 ∈ L 2 (Ω; H 1 0 (Y 2 )), u ∈ L 2 (Ω; H 1 per (Y )), Y u = 0
such that, up to subsequences still denoted with the same indices, we have the convergences:

                   T ε (u 1,ε ) u 1 weakly in L 2 (Ω × Y ),
T ε (e kl (u 1,ε )) e kl,x (u 1 ) + e kl,y (u) weakly in L 2 (Ω × Y ),

T ε (u ε,2 ) u 2 weakly in L 2 (Ω × Y 2 ), εT ε (e kl (u ε,2 )) e kl,y (u 2 ) weakly in L 2 (Ω × Y 2 ).
(ii) The three fields (u 1 , u 2 , u) solve the following three coupled variational problems:

                                             Y 1 c 1,mnkl (
y) e kl,x (u 1 (x)) + e kl,y (u(x, y)) e mn (Φ(y))dy = 0 for all Φ ∈ H 1 per (Y 1 ),

ω 2 Y 2 r 2 (y)(u 1 (x) + u 2 (x, y)) • Ψ(y)dy - Y 2 c 2,mnkl (y)
e kl,y (u 2 (x, y))e mn,y (Ψ(y))dy Proof. Weak convergences (i) are obtained by using the bounds (4), ( 5) and the properties of the unfolding operator. The limit problems of part (ii) are obtained with appropriate test-functions in problem (1). More precisely:

= -f (x) • Y 2 Ψ(y)dy for all Ψ ∈ H 1 0 (Y 2 ), ω 2 Ω u 1 (x) • Φ(x)( Y r(y)dy)dx + ω 2 Ω Φ(x) • (
-For the first problem we choose test-functions of the form:

εw(x)Φ({ x ε }), x ∈ Ω, with w ∈ D(Ω), Φ ∈ H 1 per (Y ).
-For the second problem we choose test-functions of the form:

w(x)Ψ({ x ε }) x ∈ Ω ε 2 , 0 x ∈ Ω ε 1 . with w ∈ D(Ω), Ψ ∈ H 1 0 (Y 2 ).
-For the last problem we choose test-functions of the form: Φ ∈ H 1 0 (Ω).

Step 3. In this step we solve, successively, each limit variational problem.

• First we solve the problem posed in Y 1 to compute the corrector u. The corrector u is expressed as the linear combination u(x, y) = e mn,x (u 1 (x))z mn (y) where the basis functions z mn ∈ H 1 per (Y 1 ) are solutions to the variational problems (by symmetry, there exists only 6 different problems in three-dimensional elasticity and 3 different problems in two-dimensional elasticity):

Y 1 c 1,ijkl e kl,y (z mn ) + δ kl mn e ij,y (Φ))dy = 0 for all Φ ∈ H 1 per (Y 1 ),
hence, we get the homogeneous (independent of x) tensor c * :

c * ijkl = Y 1 c 1,ijmn e mn,y (z kl ) + δ kl mn dy, (6)
and δ kl mn is the Kronecker symbol δ kl mn = 0 for m = k or n = l, and δ mn mn = 1. It is easy to show that tensor c * has the same properties of symmetry and coercivity of the initial one c ε . Let us remark that c * is independent of ω and only depends upon the value of the elasticity tensor c 1 within the matrix and the shape of the matrix, more precisely the shape of the elementary inclusion Y 1 , the same result would have been obtained with a perforated domain, Y 2 being the hole. It is important also to note that the corrector u is not determined in Y 2 .

• Next we solve the problem posed in Y 2 , this allows the computation of u 2 :

ω 2 Y 2 r 2 (y)u 2 • Ψ - Y 2 c 2,mnkl (y)e kl,y (u 2 )e mn,y (Ψ) = -ω 2 u 1 • Y 2 r 2 (y)Ψ -f • Y 2 Ψ for all Ψ ∈ H 1 0 (Y 2 ).
Let us examine the spectral properties of the previous problem in u 2 . We note {ϕ j , λ j } j≥1 the eigenelements associated to the elasticity operator. The positivity of r 2 and the coercivity of tensor c 2 imply that the eigenvalues {λ j } j≥1 are real and positive, let us range them in ascending order, 0

≤ •• ≤ λ j ≤ λ j+1 ≤ ••            Y 2 c 2,mnkl (y)e kl,y (ϕ j )e mn,y (Ψ) = λ j Y 2 r 2 (y)ϕ j • Ψ for all Ψ ∈ H 1 0 (Y 2 )
, without summation on j, and with the orthogonality condition

Y 2 r 2 (y)ϕ p • ϕ q = δ q p , . (7)
We decompose u 2 in the basis of the eigenvectors {φ j } j≥1 . Hence for all frequencies different from the resonance ω 2 = λ j , j ∈ J with

J = {j ≥ 1, Y 2 r 2 (y)ϕ j (y)dy = 0}, (8)
the displacement field u 2 can be explicitly given in terms of u 1 by the series:

u 2 (x, y) = j≥1 -ω 2 u 1 (x) • Y 2 r 2 (y)ϕ j (y) ω 2 -λ j ϕ j (y) - f (x) • Y 2 ϕ j (y)
ω 2 -λ j ϕ j (y).

• Finally we solve the problem posed in Ω to get the limit elastic field u 1 :

ω 2 Ω u 1 • Φ Y r(y) + ω 2 Ω Φ • Y 2 r 2 (y)u 2 (x, y) - Ω c * mnkl e kl (u 1 )e ij (Φ) = - Ω f • Φ for all Φ ∈ H 1 0 (Ω).
We replace u 2 by the expression obtained previously, this yields

Y 2 r 2 (y)u 2 (x, y) • Φ(x) = -ω 2 A(ω)u 1 (x) • Φ(x) -B(ω)f (x) • Φ(x),
where matrices A = (A qp ) and B = (B qp ) (of order 3) are given by:

                       A qp (ω) = j∈J Y 2 r 2 (y)ϕ j p (y) Y 2 r 2 (y)ϕ j q (y) ω 2 -λ j , B qp (ω) = j∈J Y 2 ϕ j p (y) Y 2 r 2 (y)ϕ j q (y) ω 2 -λ j ,
and a straightforward computation leads to the limit problem (2) with

A * (ω) = -ω 2 A(ω) + r * I, r * = Y 1 r 1 (y)dy + Y 2 r 2 (y)dy, B * (ω) = -ω 2 B(ω) + I. (9)
Step 4. As stated in section 5.2, for certain types of data (geometry of Y 2 and values of c ε and r ε ) it may happen that there is no zone of wave propagation, this is the case when for all j ∈ J, ω j ∇ = √ λ j+1 . Therefore to establish the existence of the solution to the limit problem, we assume that there exists a set, still denoted J, such that ω j ∇ < √ λ j+1 for all j ∈ J. With this assumption we have the following existence lemma. Lemma 3. There exists a non empty, open set W ⊂ R + such that for all frequencies ω ∈ W the limit problem (2) has a unique solution.

Proof. Matrix A * (ω) is positive definite for all ω ∈]ω j ∇ , √ λ j+1 [, j ∈ J.
In each interval of this type there exists a countable set of resonance frequencies {ω j k } k∈K j such that for all ω = ω j k the following problem has a unique vanishing solution:

ω 2 Ω A * (ω)v(x) • Φ(x)dx - Ω c * mnkl e kl (v(x))e mn (Φ(x))dx = 0 ∀Φ ∈ H 1 0 (Ω). Let W = {ω ∈]ω j ∇ , √ λ j+1 [, ω = ω j k , k ∈ K j , j ∈ J}.
Hence, by Fredholm alternative, problem (2) has a unique solution for all ω ∈ W . An easy way to assure the existence of such a set W is to increase the mass density r 1 which implies the existence of a non empty domain of frequencies where A * is definite positive.

Step 5. In this section we establish the strong convergence of the sequence

{T ε (u ε )} ε in the space L 2 (Ω × Y ).
According to Corollary of Lemma 1 and Lemma 2, up to subsequences still denoted with the same indices, we have the convergence:

u ε,1 -→ u 1 strongly in L 2 (Ω), T ε (u ε,2 ) u 2 weakly in L 2 (Ω × Y ).
Using the decomposition u ε = u ε,1 + u ε,2 we rewrite problem (1) in the following way:

ω 2 Ω ε 2 r ε u ε,2 • Φ -ε 2 Ω ε 2 c ε mnkl e kl (u ε,2 )e mn (Φ) = ε 2 Ω ε 2 c ε mnkl e kl (u ε,1 )e mn (Φ) - Ω ε 2 (f + ω 2 r ε u ε,1 ) • Φ ∀Φ ∈ H 1 0 (Ω ε 2 ).
Since the unit cells are disjointed, the previous system reduces to a problem posed in

Y 2 : For almost all x ∈ Ω find T ε (u ε,2 )(x, y) ∈ H 1 0 (Y 2 ) such that: ω 2 Y 2 r 2 T ε (u ε,2 )(x, .) • Φ - Y 2 c 2,mnkl e kl,y (T ε (u ε,2 )(x, .))e mn,y (Φ) = ε Y 2 c 2,mnkl e kl (T ε (u ε,1 )(x, .
))e mn,y (Φ) -

Y 2 T ε (f + ω 2 r 2 u ε,1 )(x, .) • Φ ∀Φ ∈ H 1 0 (Y 2 ).
We make use of the eigen basis (ϕ j ) j≥1 introduced in (7) to express the solution T ε (u ε,2 ) as:

T ε (u ε,2 )(x, y) = j≥1 ε Y 2 c 2,mnkl e kl (T ε (u ε,1 )(x, s
))e mn,y (ϕ j (s))ds

ω 2 -λ j ϕ j (y) - j≥1 Y 2 T ε (f + ω 2 r 2 u ε,1 )(x, s) • ϕ j (s)ds ω 2 -λ j ϕ j (y).
And by linearity of the unfolding operator T ε we get s))e mn,y (ϕ j (s))ds

T ε (u ε,2 ) -u 2 = j≥1 ε Y 2 c 2,mnkl e kl (T ε (u ε,1 )(x,
ω 2 -λ j ϕ j (y) - j≥1 Y 2 (T ε (f ) -f )(x, s) • ϕ j (s)ds ω 2 -λ j ϕ j (y) - j≥1 Y 2 ω 2 r 2 (T ε (u ε,1 ) -u 1 )(x, s) • ϕ j (s)ds ω 2 -λ j ϕ j (y).
Hence, from the convergence

T ε (f ) -→ f strongly in L 2 (Ω × Y ), T ε (u ε,1 ) -→ u 1 strongly in L 2 (Ω × Y ),
we get the strong convergence,

T ε (u ε,2 ) -→ u 2 strongly in L 2 (Ω × Y ), T ε (u ε ) -→ u 1 + u 2 strongly in L 2 (Ω × Y ), which implies Ω |u ε (x)| 2 dx -→ Ω×Y |u 1 (x) + u 2 (x, y)| 2 dxdy.
Step 6. We are now in a position to justify the a priori boundedness assumption of Step 1. Let us assume, by contradiction, that ||u ε || L 2 (Ω) -→ ∞ and let us con-

sider the displacement field ũε = u ε ||u ε || L 2 (Ω)
. This field is solution to a problem similar to the initial one (1) but where the right-hand side f has been replaced

by

f ε = f ||u ε || L 2 (Ω) , || f ε || L 2 (Ω) -→ 0. The sequence {|| ũε || L 2 (Ω) } ε is uniformly
bounded, the convergence Theorem 2 can be applied for all admissible frequencies, ω ∈ W , to show that the sequence {T ε ( ũε )} ε strongly converges to a vanishing displacement field, this states the contradiction.

Step 7. Finally we prove, for ε small enough, the Existence Theorem 1 for the initial problem [START_REF] Allaire | Homogenization and two-scale convergence[END_REF]. Let us proceed by contradiction. We consider a vanishing applied force f = 0 and a sequence of strictly positive numbers {ε p } p∈N * converging to zero and such that, for all p ∈ N * , there exists a displacement field u εp ∈ H 1 0 (Ω), solution to problem (1) which satisfies ||u εp || L 2 (Ω) = 1, thus:

ω 2 Ω r εp u εp • Φ - Ω c εp mnkl e kl (u εp )e mn (Φ) = 0 ∀Φ ∈ H 1 0 (Ω).
From the strong convergence of Step 5 we get

Ω |u εp | 2 -→ Ω×Y |u 1 + u 2 | 2 .
By assumption (ii), for all ω ∈ W , we have u 1 as the unique (vanishing) solution to the limit problem:

ω 2 Ω A * (ω)u 1 (x)•Φ(x)dx- Ω c * mnkl e kl (u 1 (x))e mn (Φ(x))dx = 0 ∀Φ ∈ H 1 0 (Ω)
and next the computation gives u 2 = 0 which implies a contradiction.

Numerical illustration in the two-dimensional case

In this section we provide some numerical simulations to illustrate the acoustic band gaps determined by the eigenvalues of A * (ω) and the effect of changing the parameters of the model, in particular we investigate the influence of the average material mass density and of the geometry of the microstructure.

The theoretical results obtained in the previous sections were given in the framework of three-dimensional elasticity, of course they apply as well in the twodimensional case of in-plane vibrations. Therefore, in order to reduce the computational effort (since in this case we have the explicit formulae of both eigenvalues of A * at hand, i.e., without any further computing) we restrict our work to the twodimensional case which exhibits the most important characteristics of the band gaps structure.

The numerical identification of the band gaps is done through the following steps:

1. Computation by the finite element method of the approximate value ( λj , φj ), j ∈ N of the eigenelements (λ j , ϕ j ) of the elasticity problem [START_REF] Jensen | On maximal eigenfrequency separation in twomaterial structures: The 1D and 2D scalar cases[END_REF] posed over Y 2 .

2. Determination of the reduced index set J by eliminating the eigenvalues that do not contribute to the expression of A * , thanks to the introduction of a threshold τ ,

J = {j ≥ 1, | Y 2
r 2 (y) φj (y)dy| > τ }.

3. In each interval ] λj , λj+1 [, j ∈ J, and for selected frequencies ω,

• Compute the entries of matrix à * (ω) by replacing the infinite sum over the index set J, see ( 8) by the finite sum over J, à * (ω) = j∈ J à * ,j (ω) + r * I.

• Compute explicitly the largest and the smallest eigenvalues denoted by μ∆ (ω) and μ∇ (ω) of matrix à * (ω).

• Localize numerically the frequencies denoted by ωj ∇ (respectively ωj ∆ ), for which the smallest eigenvalue (respectively the largest eigenvalue) of A * (ω) vanishes. Hence the strong and weak band gaps and the wave propagation zone can easily be identified.

The numerical examples presented below have been obtained by using an in-house software based on the MATLAB computational tools. For analysis of the eigenvalue elasticity problem [START_REF] Jensen | On maximal eigenfrequency separation in twomaterial structures: The 1D and 2D scalar cases[END_REF] defined in domain Y 2 we computed the approximation of the displacement eigenfunctions with linear finite elements on triangular meshes.

Numerical simulation of strong and weak band gaps

In Figure 1 we display the successive resonance frequencies λj , j ∈ J, for elliptic inclusions. Each frequency band ] λj , λj+1 [ is decomposed into one or two zones with no wave propagation (the strong and weak zones) followed by a wave propagation zone. In a weak band, the largest eigenvalue μ∆ of à * is positive and the other one μ∇ is negative, there is propagation only in the direction of the eigenvector ψ ∆ associated to μ∆ . This direction may change when the frequency ω varies. With the same elliptic inclusions as in Figure 1 we display, in Figure 2, the variation (with respect to ω) of the orientation angle of the eigenvector ψ ∆ . For this example, the numerical experiment shows that, within the entire weak band gap, there is no change of the direction of ψ ∆ ; hence the direction of propagation remains the same. It is worth noting that, due to this property, it makes sense for applications to use the whole weak gap interval as the propagation zone for suppressing vibrations in the direction orhogonal to ψ ∆ .

However, for more complicated geometries of the inclusions [START_REF] Miara | Shape optimization of phononic materials[END_REF], there is a change in the direction of ψ ∆ . In such cases the weak band gaps behave as the strong ones.

Influence of some microstructural parameters in the band gaps distribution

We illustrate how the acoustic bands depend on some selected features characterizing the microstructure; in particular we study the effects of changing

• the averaged material density given by r * = (1 -|Y 2 |)r 1 + |Y 2 |r 2 , when homogeneous materials with mass density r 1 and r 2 are considered respectively in Y 1 Y 2 ;

• the shape of the inclusions, i.e. the shape of Y 2 ; The value of q is set, from top to bottom, to 25%, 70%, 100% and 150%.

• the volume fraction |Y 2 |/|Y |, i.e. the (relative) size of the inclusions when keeping their shape fixed.

Averaged material density

The influence of the material density in the inclusion is rather complex. Obviously a change in density r 2 of the material in Y 2 re-scales the distribution of all resonant frequencies λj , j ∈ J and influences the magnitude of the frequencydependent part of tensors A * ,j (ω), as well as its isometric part r * I. By contrast, a change in r 1 , i.e. in the density of the matrix component, is easy to foresee. It results in a modification of the average density without any impact on the distribution of the resonant frequencies. Nevertheless, such a modification leads to a shift in the bounds of both the weak and strong gaps,whereby the quality of the gaps may change also, for example a strong band becoming a weak one or a weak one becoming a full propagation zone. This effect is captured in Figure 3, where the influence of changing the averaged density is tested on an elliptical geometry of Y 2 with r * = (1 -|Y 2 |)qr 1 + |Y 2 |r 2 the value of q is set to 25%, 70%, 100% and 150%. It can well be observed that the lighter the matrix is, the larger the band gaps are; more precisely, a smaller density r 1 results in an increase of the band gap widths.

Shape of the inclusion.

We perform the computation of the band gaps for different shapes of inclusion Y 2 for both symmetric geometries (circles, squares) and non symmetric geometries (ellipses, rectangles); in the latter cases the weak band gaps are distinguishable, in contrast with the case of symmetric inclusions where only strong band gaps (and of course propagating zones) appear.

In Figure 4, the distribution is displayed of the predicted band gaps for the first frequency band. For symmetric micro structures, i.e. those with more than two axes of the symmetry,we obtain µ j ∇ ≡ µ j ∆ , which recalls the analogy with the case of diffusion [START_REF] Bouchitté | Homogenization near resonances and artificial magnetism from dielectrics[END_REF].

For highly elongated ellipses the strong gap disappears, see also Figure 4 for the rectangular domains. All the examples displayed here were obtained with a square unit cell Y , however measurements show that other types of lattices, such as hexagonal ones, can be more appropriate to enlarge the band gaps. It would be of interest to investigate such configurations.

Volume fraction of the inclusion.

The effect of the "fill-in" coefficient |Y 2 |/|Y | of the inclusion can be studied analytically, so that having solved the eigenvalue problem in a domain Y 2 for a fixed volume fraction, the gap distribution can be predicted for other micro structures with re-scaled size of the inclusion. A change in the ratio |Y 2 |/|Y | influences the averaged density, but also the magnitudes of the resonance frequencies. When this ratio is increasing, the eigenfrequencies are decreasing and, thereby, the band gaps are "shifted" to lower bands, as illustrated in Figure 5.

As a brief conclusion for this simulation part, we remark that the numerical simulations presented in this section show the sensitivity of the distribution of the wave propagation bands with respect to some physical parameters. The next challenging question to address would be to optimize the design of the micro structures according to some figures of merit (position, reduction, enlargement, shift of the band gaps). The first step of this analysis which is the microstructural sensitivity approach is actually under study [START_REF] Miara | Shape optimization of phononic materials[END_REF].

Y 2 r 2 1 c 1 ,

 2211 (y)u 2 (x, y)dy)dx -Ω×Y mnkl (y) e kl,x (u 1 (x)) + e kl,y (u(x, y)) e mn (Φ(x))dx = -Ω f (x) • Φ(x)dx for all Φ ∈ H 1 0 (Ω).

Y 1 c 1 ,

 11 mnkl (y) e kl,x (u 1 ) + e kl,y (u) e mn (Φ) = 0 for all Φ ∈ H 1 per (Y 1 ).

Fig. 1 :

 1 Fig. 1: Band structure for elliptic inclusions. The resonance frequencies √ λ j are displayed in red. The largest μ∆ (solid) and smallest μ∇ (dashed) eigenvalues of matrix A * delineate the wave propagation zones. The bands of unlimited wave propagation are displayed in green, the strong band gaps are displayed in yellow and the weak band gaps are in white. The eigenvalues with a (almost) vanishing contribution in the set J are represented by a circle.

Fig. 2 :

 2 Fig. 2: Band structure for elliptic inclusions. Orientation angle of the eigenvector associated with the largest eigenvalue µ ∆ (ω).

Fig. 3 :

 3 Fig. 3: Influence of the average density variation r * = (1 -|Y 2 |)qr 1 + |Y 2 |r 2 on an elliptic inclusion Y 2 . The value of q is set, from top to bottom, to 25%, 70%, 100% and 150%.

Fig. 4 :

 4 Fig. 4: Influence of the shape of the inclusion Y 2 . Left: Band gaps distribution. Right: The 2nd resonance eigenmode for the corresponding shapes is illustrated in terms of the von Mises stress generated by the eigenfunction ϕ 2 . Note that the 2nd mode ϕ 2 corresponds to the eigenfrequency √ λ 2 representing the lower bound of the first strong band gap, when it exists.

Fig. 5 :

 5 Fig. 5: Influence of the volume fraction |Y 2 |/|Y |. The microstructures #1, #2, #3, #4 correspond to volume fraction equal to 21, 30, 40 and 50 %.

  positive definite, there is propagation in all directions. Thus, in each interval ] √ λ j+1 [, four situations may happen: either the whole interval is formed by only one weak band gap, or by a weak band gap ]

	√ a propagation zone ]ω j λ j , ∇ ,	√	λ	√	λ j , ω j ∇ [ followed by

j+1 [, or by a strong band gap ] √ λ j , ω j ∆ [ followed by a weak band gap ]ω j ∆ , √ λ j+1 [, or for the last configuration, by a strong band gap ] √ λ j , ω j ∆ [ followed by a weak band gap ]ω j ∆ , ω j ∇ [ and next by a propagation zone ]ω j ∇ , √ λ j+1

Latin exponents and indices take their values in the set {1,

2, 3}. Einstein convention for repeated exponents and indices is used. Bold face letters represent vectors or vector spaces.

Acknowlegment

This work has been supported by the European project "Smart Systems" HPRN-CT-2002-00284 and by the Chilean projects Fondecyt 1020298 et 7040209 .


Emmanuelle Mlayah

Lacaze

H Ayeb

J Grand

H Sellame

S Truong

J Aubard

N Felidj

A Mlayah

E Lacaze

Gold Nanoparticles in a Cholesteric Liquid Crystal Matrix: Self-organization and Localized Surface Plasmon Properties

Gold nanoparticles in a cholesteric liquid crystal matrix:

self-organization and localized surface plasmon properties

Introduction

During these last years, a large number of hybrid systems combining soft matter films and nanoparticles (NPs) have been described, in which soft matter, for example gel-like materials, polymer brushes, is often used to accurately control the inter-NPs interaction [1][2][3] . Many outcomes are in particular expected from hybrid liquid crystals (LCs)/nanoparticles (NPs) systems 4,5 . The LC and NPs properties are expected to be modified and even new properties may come out from the LC/NPs interaction, core-shell NPs embedded in nematic liquid crystals may exhibit negative refraction 6 . Several works [7][8][9] reported that ferroelectric and ferromagnetic NPs may drive the phase transitions and the optical properties of LCs under applied electric or magnetic fields [7][8] . As pointed out in these studies 7,8 , aggregation of the NPs within the LC matrix should be avoided in order to control the physical properties of the LC. Gold NPs have also been used to modify the electro-optical response or memory effects of LCs 5 . To avoid aggregation of gold NPs, usually grafting of appropriate molecules is performed on NPs 4,5 . One drawback of these gold NPs with grafted molecules is the size-dependent electron-NP surface scattering which leads to a strong surface plasmon damping and a resonance broadening, since these NPs are rarely of diameter larger than 3nm.

One main interest in hybrid LC/NPs systems lies in the possibility of controlling the metal NPs. In particular, the localized surface plasmon resonance (LSPR) of metal NPs is very sensitive to the optical index of the surrounding medium. Thus, owing to the anisotropy of the optical index of CLCs, the so-called LC birefringence, fine and reversible tuning of the LSPR wavelength and intensity by an applied electric field 10 or by light polarization 11 can be foreseen. However, a well-defined orientation of the LC molecules around the NPs is mandatory, though being difficult to realize for large NPs without any particular grafting of molecules. To achieve this proper LC orientation, high magnetic fields have been used 12 . Alignment layers, mainly rubbed polymers [13][14][15] , have been used to force the orientation of the LC molecules located in the vicinity of the NPs. Nevertheless, the question of the local molecular disorder around the NPs remains open 14- 15 . Despite the local orientational disorder of the LC molecules, a shift of LSPR associated with the variation of the incident light polarization has been evidenced in smectic LCs [14][15] . These results open the route for smart systems based on complex LC textures. Spherical gold NPs have been widely used in association with LCs due to their well-controlled synthesis. The main drawback of spherical particles of few tens of nanometers is the rather small shift of the LSPR associated with typical birefringence of liquid crystals. It is thus interesting to consider other nano-objects. Nano-materials of non spherical shape are particularly promising 11 , but still difficult to synthesize. It has been shown that NPs agglomeration may enhance the sensitivity of the LSPR to changes in the optical index 16,17 . On the other hand, anisotropic agglomeration allows for controlling the LSPR when light polarization is monitored with respect to the agglomeration orientation 18 . However, if suitable organic molecules are not grafted on NPs, when embedded in thermotropic nematic LCs, uncontrolled agglomeration of NPs takes place and eventually leads to a phase segregation between the NPs and LC 4 .

In this paper we report on the study of the optical properties of a hybrid LC/NP nanomaterial consisting of gold nanospheres embedded in a Cholesteric Liquid Crystal (CLC) matrix. We take advantage of the presence of a structural modulation close to the CLC/air interface, to control the NPs optical properties. We address the following experimental issues: (i) What is the impact of the NPs distribution induced by the LC, on their surface plasmon properties and (ii) could the LC be used for tailoring these properties? To carry out this study we have combined optical imaging, Scanning Electronic Microscopy (SEM), extinction spectroscopy and numerical simulations.

Experimental method 2.1 Preparation and optical imaging

Polyvinyl alcohol (PVA) was spin coated onto clean glass slides and subsequently rubbed unidirectionally to allow for planar LC alignment. The CLC was prepared by mixing 2% in weight of a chiral agent S-2011 (also known as MLC-6247 from Merck), with the nematic liquid crystal E7 (from Merck) in order to set the helix pitch of the CLC film at p= 5.5µm. The E7 liquid crystal is in the nematic phase for temperature ranging from Tmin =-20°C to Tmax=58°C with a high optical anisotropy ∆n = ne -no =0.23 (ne=1.75, parallel to optical axis and no=1.52, perpendicular to optical axis, where ne and no are respectively the extraordinary and ordinary optical indices). The prepared CLC was spin-coated on the unidirectionally rubbed PVA layer to obtain a CLC film with an average thickness of 12µm. The CLC film thus yields a planar anchoring onto the substrate and a homeotropic arrangement (perpendicular CLC director) at the CLC/air interface. The cholesteric texture was observed using Optical Microscopy (Leica DMRX) and crossed polarizer and analyzer. Periodic dark and bright textural stripes can be seen on a large scale (Fig 1).

This CLC texture has been interpreted earlier by Baudry et al. 19 using a phenomenological model. According to this model, the CLC director, which corresponds to the average molecular orientation, rotates, in the bulk of the film, in a plane parallel to the substrate; the cholesteric helix is thus perpendicular to the substrate (Fig 1 . b). While approaching the CLC/air interface, the CLC helix rotates in order to let the director match the homeotropic arrangement at the surface. This rotation occurs within an elastic boundary layer of thickness p/2, as shown in the schematic of Fig 1 .b 19,20 . In this layer, the orientation of the molecules alternatively changes from perpendicular to parallel with respect to the plane of the substrate, thus forming stripes of alternating planar and homeotropic molecular orientations. These stripes are visible under optical microscopy (fig. 1a) and are associated with optical index varying between n o and n e when the light polarization is parallel to the stripes. Their orientation is determined by the planar anchoring and the CLC film thickness 19,20 . Λ is the distance between two dark or bright lines. One expects Λ ~ 2p with p the CLC pitch 19,20 . From perpendicular to the stripes, along the white dashed line. In the lower part of the sample, the director is in a planar orientation and the cholesteric helix is perpendicular to the substrate. When approaching the CLC/air interface, the helix rotates and the director matches the homeotropic anchoring at the LC/air interface. This leads to the presence of a boundary layer of thickness p/2 in which the CLC helix is parallel to the substrate 19,20 . 20 nm spherical NPs coated with citrate ions and prepared in water solution (25.10 14 NPs/l) were synthesized following Reference 21 . PolyVinylPyrrolidone (PVP at 0.2% in wt) was added in order to avoid NPs aggregation in water 22 . To improve the miscibility of the NPs in the liquid crystal, 50% of acetonitril was added to the NPs/0.2%PVP solution. This NPs concentration corresponds to the sample referred to as S1. Another sample, S2, with a NP concentration 10 times larger than S1, was also studied. 20 µl of NPs solution was dropped onto the liquid crystal heated up at 55°C. This temperature was chosen to favor the NP mobility and to have the LC in its cholesteric phase. In all cases, the surface modulation of the CLC texture (Fig 1) is preserved after incorporation of the NPs. The modulation period is increased to about 15 µm compared to the initial value Λ= 12.4µm with no agglomeration between NPs visible by Optical Microscopy.

Extinction spectroscopy and SEM

The optical properties of the hybrid CLC/NPs system have been investigated by micro-spectrometry. The sample is illuminated using linearly polarized white light. In order to avoid the twisted structure of the CLC film that is likely to change the polarization of the incoming light, we performed the measurements by illuminating the CLC surface. The transmitted light is collected with a 50X objective (0.55 numerical aperture) and dispersed with a LOT Oriel MS260i spectro-photometer.

Along with a 2 mm pinhole, this objective allows to probe a 100x100 µm 2 area. In order to extract the extinction spectrum of the gold NPs, a CLC film without gold NPs was used as reference. Both the CLC/NPs sample and the CLC reference were carefully aligned with respect to the incident polarization. All spectra were collected at room temperature in the 400-1000 nm spectral range.

We have also performed SEM measurements on a field emission, Zeiss Supra 40, Scanning Electron Microscope (SEM).

Since a conductive substrate is required, we have prepared a sample similar to S2 but with the CLC film deposited on top of a conductive MoS 2 substrate. Planar anchoring 23 and formation of modulation stripes were confirmed by optical microscopy.

After deposition of the NPs the sample was heated at 80°C, in isotropic phase, in order to let the CLC slowly penetrate the bulk layered structure of MoS 2 until it is completely drained. The agglomeration of the NPs, directly deposited on top of the substrate, could then be observed by SEM. It is now a well-documented fact that the electromagnetic interaction between metal NPs leads to a hybridization of the surface plasmons which resonance wavelength increases with decreasing separation between the NPs (i.e. with increasing interaction) 18,[24][25][26][27] . Such interaction is responsible for red shift, broadening and asymmetry of the LSPR as observed in Fig is indeed known that a localization close to the air/ interface can correspond to a stable localization for NPs embedded in cholesteric materials 28,29 . The presence of the NPs close to the CLC/air interface modifies the LC anchoring at air, which explains the observed variation of the modulation period Λ upon deposition of the NPs 19,20 .

Results and discussion

1 Non interacting and weakly interacting nanoparticles

In order to determine the environment of the NPs, we have simulated the optical extinction spectra using the discrete dipole approximation (DDA) as implemented in ddscat software 30 . We consider an isolated spherical gold NP, 20 nm in diameter coated with a shell of PVP 31 . This core-shell NP is surrounded by a medium with optical index n o = 1.52 or n e = 1.75. We used the dielectric properties of bulk gold 32 corrected for the size dependent plasmon damping 33 .

The The initial concentration of NPs in the colloidal solution has been increased by a factor 10 and a drop was deposited on the CLC film (sample S2). The SEM measurements (Figure 3) show that NPs clusters (from few tens to few hundreds of nanometers in size, with an average length of the order of 10 NPs) are formed as well as NPs chains. In these aggregates, the separation between NPs is well defined; the average separation was determined 0.6 ± 0.5 nm from the SEM.

Optical spectra were acquired from 10 different areas of sample S2 in order to probe the sample homogeneity. Unlike the case of sample 1, we found a good homogeneity attested by small spectral changes, suggesting that the solvent evaporation becomes more homogeneous when the NPs concentration is increased. Figure 4 shows typical extinction spectra measured with incident polarization either parallel or perpendicular to the modulation stripes. The spectra from sample S2 exhibit two LSPRs (Fig. 4a). One is located at 536 ±2 nm, and was already observed in sample S1 (Fig. 2a). In sample S2, there is an Since in sample S2 the NPs density is 10 times higher than in sample S1, interactions between the NPs are more likely to occur and it is well-known that such interactions may lead to red-shifted LSPRs 18,[24][25][26][27]34 . In order to relate the LSPR band at 706 nm (Fig 4a) to an average separation between NPs we have performed DDA simulations. The fact that this LSPR comes out for both polarizations suggests that interactions between NPs take place in these both directions. Therefore, we considered two types of NPs arrangements: (i) NPs chains with randomly distributed orientations and (ii) NPs clusters. These two types of arrangements were indeed observed by SEM.

According to previous theoretical studies 35,36 Figure 5a shows the calculated spectra for nanoparticle chains with 1nm gap between nanoparticles. The incident polarization is parallel to the chain and the optical index of the surrounding medium is ne eff . As already reported in several studies 18,[24][25][26][27] , because of the electromagnetic interaction between the nanoparticles, the collective electronic oscillations along the chain give rise to a longitudinal surface plasmon resonance. The resonance wavelength is red shifted with increasing length of the chain as shown in Fig 5a . The red shift becomes very pronounced with decreasing separation between the NPs (see inset in Fig. 5a). The interesting point is that, for a given separation between the NPs, the red shift comes to saturation with increasing chain length [37][38][39] . This saturation effect occurs because the surface plasmons have a finite coherence length that decreases with increasing separation between neighboring NPs.

For polarization perpendicular to the stripes, the NPs embedded in the modulated texture should experience the effect of a surrounding medium with no eff as well as the NP located at the LC/air interface. In the former case, the LC director rotates around an axis perpendicular to the stripe and is therefore always perpendicular to the polarization (see figure 1). Figure 5b shows that, for 21 NPs chains and gaps between NPs of 1nm, when the optical index is changed from ne eff to no eff , the LSPR wavelength drops down to 670nm. This value is somewhat smaller but close to the experimental value of 706 nm suggesting a separation between the NPs smaller than 1nm and a PVP shell compressed in the gap region between NPs, which could be indeed the case since PVP consists in quite flexible organic chains. This result is in agreement with the SEM images, which indicate a separation of 0.6 nm between NPs. We have thus calculated the LSPR wavelength at saturation (31 NPs) for a 0.6nm separation between NPs and an optical index no eff . We found a value λ = 753nm which is higher than the measured one. This indicates that the average chain length is smaller than the considered chain length (31 NPs). As mentioned above, the NPs could also form aggregates. The interaction between the NPs within an aggregate may lead to a red-shifted LSPR as well. Therefore we simulated the extinction spectra of aggregates consisting of 21 NPs with separations between first neighboring NPs ranging from 1 to 3.5 nm (Figure 6). The NPs inside each cluster are distributed as shown in Figure 6c. Both red-shifted and single particle like LSPR come out (for both incident polarizations) because the cluster contains NPs chains, that are both parallel and perpendicular to the incident polarization, and quasi-isolated NPs (at the cluster periphery) with resonance close to the 536 nm isolated NP LSPR. As expected, the resonance wavelength of the red-shifted LSPR decreases with increasing separation between the NPs (Figure 6). It is around 640 nm for 1nm separation between first neighboring NPs and for no eff =1.38. This value is smaller than the saturation wavelength 670nm calculated for the chain consisting of 21 NPs separated by 1 nm, because the 21 NPs cluster has a smaller size than the equivalent chain.

wavelength is very close to the one of a 5-NPs chain (fig. 5a and6b). This means that to account for the measured LSPR wavelength (706 nm ± 8 nm), 0.6 nm separation between NPs is required as shown by SEM results, but also cluster lengths in average smaller than 31 NPs, still in agreement with SEM results (Figure 3).

Let us now comment on the dependence of the red-shifted LSPR lineshape on the polarization orientation reported in Fig. F o r ( i ) , w h e n t h e p o l a r i z a t i o n i s c h a n g e d f r o m p e rpendicular to parallel to the modulation stripes, the optical index increases within the modulation texture and could lead to a red-shift of the 706nm LSPR, as evidenced by the simulations (fig. 5b, fig. 6). However it appears that the separation between NPs is around 0.6 nm whereas the LC molecular length is of the order of 2 nm. Thus LC molecules cannot be parallel to the coupling direction due to the small inter-particle space. In other words, the optical environment in the gap between NPs should be associated with an optical index close to no eff corresponding to a LC optical axis tilted with respect to the coupling direction. Therefore, we do not expect a significant influence of the birefringence in the case of strongly coupled NPs.

For (ii), the modulated texture of the LC could drive the formation of NPs chains and/or clusters orientated along the stripes. NPs chains aligned along the stripes should lead to red-shifted LSPR associated with increase of extinction when polarization is changed from perpendicular to parallel to the stripes, in agreement with the experimental observations in figure 4. This is also expected for elongated clusters orientated along the stripes. Therefore, we may conclude that in our system, part of the NPs chains and/or clusters is oriented by the modulation stripes associated with the LC texture. According to the SEM, chains and anisotropic clusters are indeed present. However, no preferential orientation has been evidenced (Fig. 3). On MoS 2 , the orientation of the NPs chains and clusters may have been lost due to the spreading/removal of the CLC required for SEM visualization.

These results finally demonstrate that, for high concentration, the NPs penetrate the modulated texture. It has been reported for small NPs (diameter 3nm), that demixing of NPs in CLC is favored, thus leading to localization of the NPs at the air/CLC interface 28,29 . Such demixing may also take place in our system as observed at low NPs concentration (sample S1).

However this demixing may be competing with cluster formation induced by Van der Waals interactions between the NPs.

Large NPs clusters, observed by optical microscopy, have been shown to coexist with NPs self-organization at the CLC/air interface for small NPs 27 . However, in our system, the cluster size remains limited, of the order of 10NPs in length, in agreement with the fact that they remain invisible by optical microscopy. It is known that in nematic LCs 4,5 , contrary to smectics 15 , aggregation between NPs occurs easily. Due to the large CLC pitch (compared to the NPs diameter), the environment of the NPs is close to a nematic environment, but no segregation occurs. The size of the clusters remains limited with a preferred orientation along the modulation stripes. This phenomenon may be associated with the attraction between LC defects and colloids in nematics 40 which is known to increase with the size of the colloids 41 . A critical size of the NPs clusters may be defined, below which aggregation between NPs is favored, and above which attraction occurs between clusters of NPs and the highly deformed area located between bulk and modulation. This would lead to the observed stabilization of the hybrid system with an anisotropy of the clusters, preferably oriented along the modulation stripes. In order to confirm and better explore this phenomenon, further studies in which cholesteric pitch would be varied are now highly desirable.

Conclusion

We have investigated a hybrid LC/NPs nano-material that makes use of the LC texture modulation, characteristic of CLC films, for the control of NPs agglomeration and localisation. We found that, owing to the penetration of the NPs within the CLC texture at high concentrations, strong interactions between the NPs take place. This was evidenced by the observation of a strongly red-shifted LSPR that exhibited sensitivity to the polarization of the incident light. The interaction between the NPs within simple configurations of chains and regular aggregates has been investigated using DDA simulations. The latter allowed for addressing the main physical parameters that determine the red-shifted LSPR wavelength (chain length, aggregate size, inter-particle separation, optical index of surrounding medium). Our results show that in the hybrid LC/NPs nanomaterial the control of the optical response of the metal NPs could be achieved not only by monitoring the LC birefringence, but also by driving the formation of NPs aggregates using the LC texture. Our results point out the potential use of the molecular organization in LCs for the engineering of new smart nanomaterials.

  Fig 1a we find Λ = 12.4 ± 0.3 µm which is indeed slightly larger than 2p = 11 µm. Since the modulation orientation can be controlled by the CLC film thickness, our aim is to take advantage of the CLC modulated texture (Fig 1a) in order to monitor the NPs LSPR.

Figure 1 .

 1 Figure 1. Upper panel: Texture of the CLC film deposited on a glass plate coated with a rubbed PVA layer. The texture shows periodic dark and bright lines (modulation).

( a )

 a Zoom on the modulation texture. The dashed lines mark the CLC modulation period Λ = 12.4 µm. (b) Schematic of the cholesteric director rotation in side view. According to Baudry's model, the figure shows a cross-section in a plane

Figure 2 .

 2 Figure 2. Extinction spectra (b and d) measured on sample S1 for incident light polarizations parallel (red) and perpendicular (black) to the modulation stripes. Spectra (b) and (d) were recorded on two different points of the sample. Their first derivative are plotted in (a) and (c), respectively. In (d) are also shown the spectra calculated for 20 nm gold NP coated with 2nm PVP and surrounded by a medium with optical index n o = 1.52 (blue) or n e = 1.75 (green).

Figures

  Figures2b and 2dshow the extinction spectra of the hybrid CLC/20nm gold NPs film recorded on two different typical areas of sample S1 and for two perpendicular incident polarizations oriented with respect to the modulation stripes. From the inflexion points of the first derivative spectra shown in Fig2a and 2c we found the LSPRs at 544 ± 1nm and 534 ± 1nm, respectively (for both polarizations). The LSPRs are red-shifted with respect to that of NPs in water which has been found at 525 ± 1nm. It is well-known that increasing the optical index of the surrounding medium (from water to CLC) leads to a redshift of the LSPR, in agreement with a CLC environment around the NPs. However, here the LSPR red-shift depends on the

  Fig 2c,d and Fig 2a,b, respectively). Indeed, weakly interacting NPs with randomly fluctuating inter-particles separations give rise to a spectral distribution of red-shifted surface plasmon resonances. The incident light polarization has been oriented with respect to the modulation stripes (see Fig.1) in order to monitor the evolution of the LSPR with respect to the change of the CLC optical index from n e to n o . The extinction spectra recorded with incident light polarization parallel and perpendicular to the CLC modulation are very similar (Fig 2): the LSPRs are located at nearly the same wavelength and have the same intensity for both polarizations (Fig2b and 2d). This indicates that the NPs did not penetrate the modulated texture and should be located close to the CLC/air interface in a homeotropic environment. It

Figure 3 .

 3 Figure 3. SEM image of sample prepared in the same conditions as sample S2, and deposited on a conducting MoS 2 substrate. The CLC has been removed prior to SEM measurements.

  additional LSPR around 706 ± 8 nm. The indicated uncertainty reflects the LSPR wavelength variation with the probed areas on the sample for 10 measured areas. In contrast with sample S1, sample S2 exhibits changes of the LSPRs with orientation of the incident light polarization. No clear evidence for a red-shift of the LPSRs is found when changing the polarization from parallel to perpendicular with respect to the modulation stripes: the LSPR wavelengths, as determined from the first derivative of the extinction spectra (Fig 4b), are nearly the same for both polarizations. However, one notices a change of the red-shifted LSPR lineshape with polarization. Indeed, for polarization parallel to the modulation stripes, the extinction spectrum systematically exhibits an excess of red components (Fig 4a) as shown by the red-shifted band in the difference spectrum (Fig 4c). Moreover, for nearly half of the probed areas, when the incident polarization is changed from perpendicular to parallel to the modulation stripes, we observe an increase of the LSPRs intensities at 536 and 706 nm. This intensity dependence is systematically observed when both LSPR bands are spectrally well separated as in Fig 4a. This effect is pointed out in the difference spectrum shown in Fig 4c and was completely absent for isolated and weakly interacting NPs (Fig 2).

Figure 4 .

 4 Figure 4. Left: (a) Extinction spectra from sample S2. The incident polarization is either parallel (red) or perpendicular

  , small NP aggregates, like NPs dimers and trimers, cannot account for the large LSPR red-shift observed in Fig 4a, while keeping within physically reasonable inter-particle separations. This suggests that the NPs agglomerates should contain a large number of NPs. Therefore, we have calculated the optical extinction spectra of NPs chains and aggregates containing up to 21 NPs. The calculations were performed for 20 nm gold NPs without the PVP shell in order to keep within a reasonable computation time. Nevertheless, instead of using n o and n e for the optical indices of the surrounding medium, we have determined effective optical indices no eff = 1.38 and ne eff = 1.53. These valuesgive the same LSPR wavelengths for a 20 nm gold NP and for a 20 nm gold/2 nm PVP core-shell NP surrounded by medium with optical index equal to n o or n e . In that way the effect of the PVP shell on the LSPR wavelength is somehow taken into account.

Figure 5 .

 5 Figure 5. (a) Simulated extinction spectra for NP chains separated by 1 nm gap, in an optical index ne eff . The number of NPs within the chain is increased from N = 2 (dimer, in the lower trace) to N = 20 by increment of 2. The spectra are normalized to the geometric section πa 2 , a being the radius of the sphere of equivalent volume. The inset shows the red-shift of the LSPR wavelength with increasing N and for 3 nm (blue), 2 nm (green), 1.5 nm (orange) and 1 nm (red) separations between the NPs. (b) Simulated extinction spectra for 21-NPs chains for 1.5 nm (blue lines) and 1 nm (red lines) gaps between NPs and embedded either in no eff =1.38 (full lines) or ne eff = 1.53 (dashed lines).

Figure 6 .

 6 Figure 6. Simulated extinction spectra of NPs aggregates embedded either in no eff = 1.38 (a) or ne eff =1.53 (b). The spectra are normalized to the geometric section πa 2 , a being the radius of the sphere of equivalent volume. The NPs size is 20 nm and the separation between first neighboring NPs has been varied from 1 nm to 3.5 nm as indicated by the numbers attached to each spectrum. (c) Schematic of the model aggregate.

4 .

 4 As shown in the difference spectrum in Fig 4c, changing the polarization from perpendicular to parallel to the modulation stripes gives rise to additional red spectral components in the extinction spectrum. This polarization effect can be accounted for by (i) the birefringent character of the modulated texture, (ii) an anisotropic distribution of the NPs agglomerations connected with the orientation of the modulation stripes.

Journal of Materials Chemistry

Acknowledgements

This work was supported by the French National Research Agency (ANR contract N° ANR-07-Nano-032, "NP/CL") and by CALMIP Computing facilities at Paul Sabatier University of Toulouse.

Notes and references


Aravinthan Gopalasingham
email: gopalasingham.aravinthan@alcatel-lucent.com

Laurent Roullet
email: laurent.roullet@alcatel-lucent.com

Nessrine Trabelsi
email: nessrine.trabelsi@inria.fr

Chung Shue

Chen Abdelkrim Hebbar

Erick Bizouarn
email: erick.bizouarn@alcatel-lucent.com

Generalized Software Defined Network Platform for Radio Access Networks

Keywords: SDN, RAN, LTE/5G, OpenDaylight, optimization

Software defined networking (SDN) is growing rapidly in telecommunications due to its capability to efficiently manage end-to-end networks by decoupling control plane and data plane. The scalability and flexibility can bring benefits to network management and maintenance. However, it has yet to happen for wireless networks since there is a lack of practical solutions for migrating Radio Access Networks (RAN) to software defined wireless networks. Here, we design a wireless SDN controller based on OpenDayLight (ODL). We show how our framework can support typical use case and Self-Optimized Network algorithm. The design is highly flexible and can be leveraged for future wireless networks such as C-RAN and 5G.

I. INTRODUCTION

Wireless data traffic has been growing very rapidly during the last few years due to the dramatic increase of wireless devices, online contents, and new applications. It is predicted that there will be tremendous increase in the number of connected devices by 2020 and mobile broadband network operators will face the challenge of delivering up to 1GB/user/day [START_REF]Nokia Networks APT700[END_REF]. This huge demand will require significant network capacity increase and more importantly to adapt dynamically to various traffics and user patterns.

The introduction of LTE HetNet is the well known solution for increasing network capacity through higher spectrum spatial reuse by deploying low power small cells (SC) mixed with macro cells (MC). It is expected that 5G will be deployed as heterogeneous network, consisting of various radio access technologies. Besides, there will be a combination of macro and small cells of various size [START_REF]Paper: 5th Generation of Communication Networks[END_REF]. This would generally bring additional management and operational complexity.

3GPP Self-Organizing Network (SON) introduced in LTE is known for simplifying network operations and optimizing radio access network (RAN) through automated mechanism [START_REF] Chen | Joint optimization of radio resources in small and macro cell networks[END_REF]. For example, SON supports low frequency method (e.g., self-configuration at power up, per day self-organization) and also high frequency method (e.g., distributed interference coordination via X2 interface).

Software defined networking (SDN) is the emerging technology for designing, building and managing networks, in decoupling control plane and data plane to better optimize each other. In fixed networks, telecom operators and vendors have already started to align their network management systems (NMS) with SDN paradigm: SDN provides operational advantages like high flexibility and programmability, enables operators to have a complete control over the network from a centralized vantage point, and facilitates operators to easily deploy various applications, services and fine-tuned network policies. Though there is significant research for aligning SDN with RAN (see, e.g., [START_REF] Bansal | OpenRadio: a programmable wireless dataplane[END_REF]- [START_REF] Yap | OpenRoads: empowering research in mobile networks[END_REF]), currently none of the mainstream SDN controllers provides support for 4G RAN. This is disappointing since SDN controller could be both a potential alternate for 3GPP NMS as well as a nice tool for supporting advanced SON features. SDN would be also beneficial to future 5G since SDN can control 5G system from end to end, from access network to core network.

Inspired by SDN and OpenDayLight (ODL) [7], we propose in this paper a SDN platform for RAN with network abstraction. We define a new south bound (SB) protocol called "Radio Net Flow" using 3GPP Stream Control Transmission Protocol (SCTP). We build a prototype using ODL controller and emulated base stations (BSs). We validate our framework with Enhanced inter-cell interference coordination (eICIC) deployed as a north bound (NB) application.

II. ARCHITECTURE AND IMPLEMENTATION

A. Architecture Overview

Our architecture closely follows the SDN paradigm by providing global network view using a logically centralized controller. It is illustrated by Fig. 1. Since the focus of the architecture is to provide vendor-agnostic North Bound (NB) APIs for supporting SON algorithms and optimization, the controller has to provide abstraction of network Key Performance Indicators (KPI) and configuration parameters. The RAN applications would use the controller's NB APIs to:

• discover existing network elements and topology, • retrieve RAN measurements and configuration parameters, and • re-configure RAN using optimal parameters. Radio Net Flow (RNF) is the SB protocol introduced in the ODL controller as a RAN configuration protocol for (i) collecting radio network measurements and configurations, and (ii) re-configuring parameters of those radio elements. We implemented RNF using SCTP since it is the standard transport layer protocol in mobile network architecture, which can facilitate the integration of our SDN controller with standard RAN [START_REF] Ong | An introduction to the stream control transmission protocol (SCTP)[END_REF]. We developed the corresponding protocol agent called Radio Net Flow Agent (RNFA). RNFA can be integrated into any eNodeB as a communication interface to SDN controller. We introduced three service modules in the Service Abstraction Layer (SAL) of ODL: the RAN Configuration Manager, the RAN Statistics Manager, and the RAN Topology Manager. All these service modules interface with the RNF SB protocol using specific libraries. Through RNF, they can collect data from the RAN and also send configurations to the RAN.

Since the SAL imposes transient data storage, we introduced RAN Inventory that is an external database to maintain long term network state history. We used No-SQL database for performance considerations. During operations, RAN service modules update the RAN Inventory. On the other hand, NB applications access the RAN Inventory for information retrieval and settings via NB REST (HTTP RESTful compliant) APIs. The RAN Inventory is illustrated by Fig. 2. It consists of the following information for NB application to deploy optimization:

• Network Topology: A complete record of network topology contains all the macro and small cells connected to the controller, together with their unique IDs and neighboring cell information. The RAN Topology Manager generates such a graph in the logical table format. The record remains static until there is a change in the network topology (e.g., BS addition/removal/failure). parameters periodically to the base stations.

B. SDN Use Cases for RAN

Our SDN platform is prototyped to support RAN optimization. As shown in Table I, there are three potential use cases of SDN for RAN: radio optimization, network optimization, and load balancing/virtualization. Regarding RAN, we can categorize it into two major classes in terms of coordination requirement: Centralized RAN (C-RAN) and Distributed RAN (D-RAN). The SDN platform can provide dynamic control over fronthaul transport network and also mobile backhaul network to allocate available capacity, control latency, and meet QoS requirements. Besides, the network abstraction using ODL with multiple SB protocols can be used to adjust transport network resources (flow control technique based on OpenFlow [START_REF] Mckeown | OpenFlow: Enabling innovation in campus networks[END_REF] to control bandwidth and QoS assignment) according to the network requirements, e.g., to tackle temporary traffic peak and provide dynamic capacity.

Since the controller has the knowledge of RAN elements, it can trigger handover of mobile users in particular celledge users to lightly loaded neighboring BSs for more radio resource and throughput enhancement. This is a typical load balancing scenario that can be easily achieved by using SDN controller with dynamic resource allocation in fast recovery/response during failure and re-configuration among BBUs. In the following, we describe two use cases of our design.

1) SDN for eICIC: eICIC introduced in 3GPP-LTE is to conduct a time repartition of the physical resource blocks (PRBs) such that during some time slots called Almost Blank Subframe (ABS), macro cells will mute and do not transmit any data traffic. This allows to mitigate strong inter-cell interference that may severely degrade cell-edge user performance. Besides, in HetNets, due to the large difference in the transmit power between macro and small cells, user-cell association should be conducted in a more intelligent manner but not simply be decided due to the strongest received signal. In order to favor user association to small cell against macro cell, Cell Individual Offset (CIO) can be employed. Load balancing between Creation, configuration Balancing/ base stations [START_REF] Gudipati | SoftRAN: software defined radio access network[END_REF] and load balancing Virtualization between vBBUs [START_REF] Yang | OpenRAN: a software-defined RAN architecture via virtualization[END_REF] Our eICIC optimizer is based on a 2-tier model composed of macro and small cells that perform dynamic user-cell association and transmission scheduling optimization, which is to determine optimal CIO and ABS settings for each BS. We use a potential game best-response approach [START_REF] Singh | Distributed downlink resource allocation in cellular networks through spatial adaptive play[END_REF] for network utility maximization. During each optimization iteration at randomly chosen BS, we choose the setting that maximizes a localized potential function, which involves the BS and its neighboring cells. To facilitate the above optimization, each UE reports to its serving BS their experienced channel quality, signal level and interference power. These are sent via the SDN controller to the logically centralized eICIC optimizer.

Through iterations, the global utility of the network improves gradually and converges to an optimal setup. Note that here we can plug different optimization methods (e.g., heuristic or strict global optimization) as applications to the SDN controller, depending on system criteria and cost.

Note that this 2-tier model permits the separation of control decisions between the SDN controller and the local operation of eNBs, so as to facilitate flexible and dynamic network management and optimization. The parameters that would influence the neighboring eNBs need some central coordination in the controller, shown in Fig. 3. Here, CIO re-configuration can trigger handovers and maintain a load balanced network. Determining ABS ratio can coordinate the use of PRBs and mitigate inter-cell interference.

2) Dynamic Point Selection: LTE contains several coordinated multipoint transmission techniques (CoMP) in performing dynamic coordination of transmission and reception among BSs. It aims to enhance system overall performance, resource utilization efficiency, and user service quality (e.g., cell-edge user throughput). However, CoMP requires close coordination between geographically separated eNBs. The main techniques of CoMP are coordinated scheduling/beamforming (CS/CB), dynamic point selection (DPS), and joint transmission (JT). For example, DPS can occur given that there is coordination among BSs so UE can select the best cell in a very fast time scale or even on a subframe basis. This would require a communication and instruction of from which cell the network should direct its data to and also transmission synchronization. SDN is the choice for DPS such that one can maintain a fast switching and dynamic backhauling between the eNBs and benefit from the diversity gain when choosing the best serving point dynamically, e.g., according to channel conditions.

In our proposed architecture, CS/CB and DPS can be implemented as NB applications where they get the required input parameters such as neighbor cell list (NCL), radio resource

III. VALIDATION OF THE FRAMEWORK

To justify the proposed SDN and system design, we tested the first use case: the eICIC optimizer over LTE HetNet delivering optimized CIO and ABS values. We use a Matlab based LTE-compliant simulator developed by TU-Wien to emulate the cellular networks. It generates macro and small cells, mobile users, radio propagation and pathloss, and 3GPP-LTE setup. The network is iterated per Transmission Time Interval (TTI), which generates time varying link quality and traffics, and each eNB runs real-time resource scheduling.

At some point, the MATLAB simulator (shown in Fig. 3) connects with the eNBs emulators via UDP messages to exchange measurements and information. The IP addresses and ports of the emulated eNBs are known by the simulator. An initial connection is done to send the information of network topology (including a list of all cells in the cluster, their neighboring eNBs, the maximum transmission power of each cell). During optimization iteration, the Matlab simulator creates UDP socket and sends any updated messages for each eNB: the CIO values, the ABS patterns, the attachment of mobiles users, and the channel conditions (by the CQI from each UE to its serving cell and neighboring cells). All these messages are then sent via SCTP from the eNBs to the ODL to be stored at the RAN Inventory. The eICIC optimizer therefore sends requests using REST APIs to retrieve the information needed to conduct optimization iterations. The parameters of the tested LTE HetNet for simulation are given in Table II.

The optimization delivers the best CIO values for each cell. Before optimization, the CIO value at each eNB is by default 0. After eICIC optimization, the CIO value at a small cell can be increased to favor attaching to small cell. Meanwhile, the ABS ratio of macro cell is increased to reduce its interference to neighboring small cells in some PRBs. This is to offload UEs to small cells and improve radio resource utilization efficiency with interference management. Via SDN controller, at each TTI, new CIO and ABS values are sent to the different eNBs to update the user attachment and inter-cell interference coordination and tune the system dynamically.

IV. DISCUSSIONS & CONCLUSION

In the performance validation of the framework, we used a 16-core server for running the SDN controller, RAN Inventory and eICIC optimizer. We also measure the latency between the control node and BSs. Fig. 4 shows the result against the number of BSs activated. For 30 BSs, it takes about 2.7 ms for completely updating the database for each cycle of updates from all the BSs to reach the RAN Inventory and about 1.5 ms for sending back the re-configuration parameters from the SDN controller to reach all the BSs after eICIC optimization. In the validation scenario, we have time budget of 50 ms for applying optimal parameters to the BSs, so our SDN design can satisfy the time requirement for eICIC optimization.

As shown in Fig. 4, the latency increases steadily when the number of BSs increases. This indicates that the proposed architecture is scalable. Note that when adding a new BS to the controller, the latency is increased by some nanoseconds, which is mainly due to the time taken in updating the centralized database and querying the database for re-configuration parameters to the BSs. Indeed, it is possible to improve the performance for stringent time budget optimization algorithms by adding more hardware resources at the server for running the controller and database or applying techniques such as selecting proper schema and indexes for the applications and using high-performance storage. This could be a future work.

In the proposed design, we have deployed three different entities in the control node: the controller, the database, and the application. The motivation of using centralized database as an inventory rather than using XML-based yang data model in ODL is to generate historical data as an input for data analytic tools and applications, which can facilitate big data analytics and anticipatory network management. The implementation of REST APIs for NB applications is to enable the deployment of RAN applications in remote server which can be developed in various programming languages. The SB protocol plug-in RNF is implemented using SCTP since SCTP has unique feature of creating multiple and parallel message streams between two end points which makes SCTP robust in comparison to TCP and UDP. Though our SDN platform is primarily developed for supporting LTE SON methods, it can be further extended to support different use cases in Cloud RAN and future 5G.

Fig. 1 :

 1 Fig. 1: SDN platform for RAN architecture.

Fig. 2 :

 2 Fig. 2: RAN Inventory in MongoDB for eICIC optimization.

Fig. 3 :

 3 Fig. 3: Prototype for performance validation. availability, backhauling information (data traffic and also CSI/control traffic) and transmission power setting from the RAN Inventory via REST APIs. Since in ODL, the NB applications can use multiple SB protocols for producing the global abstraction of LTE networks, CS/CB and DPS can easily use our RNF protocol (built on SCTP) for abstracting RAN and use OpenFlow protocol for backhaul network configuration to achieve CoMP and enhanced QoS.

Fig. 4 :

 4 Fig. 4: Latency in the SDN Architecture.

•

  User Attachment and Interference Information: The RAN Statistic Manager is responsible for populating this information in the RAN Inventory. It is updated via RNF protocol once a user performs handover or enters/exits the network. Interference information refers to the interference situation experienced by each user, indicating the user ID, its BS, interfering BSs, and the level of interference. It is sent to the controller periodically.

• Configuration Parameters: The RAN Configuration Manager handles configuration parameters in the RAN Inventory. NB applications store optimized parameters for each BS. The RAN Configuration Manager sends these

TABLE I :

 I SDN platform for RAN use cases

		D-RAN	C-RAN
	Radio	eICIC/ICIC	CoMP
	Optimization	(More generally dynamic/ (e.g., Coordinated scheduling,
		fast SON)	dynamic point selection [11])
	Network	Optimization of mobile	Fronthaul network
	Optimization	backhaul network [12]	optimization [13]
	Load

TABLE II :

 II eICIC simulation parameters Radio 10MHz LTE Release 8 SISO, 3GPP TS36.942 -recommended pathloss and shadow fading models Topology of eNBs Hexagonal, a macro (MC) has 3 sectors & 6 SCs Antenna (MC/SC): Kathrein/omnidirectional Max Power (MC/SC): 40W/1W Mobile users 25 UEs per MC, 10 UEs per SC, total 135 UEs, No mobility, Full buffer model Optimization values CIO choices =[0, 10, -10] ABS choices = [0%, 20%, 40%]

ACKNOWLEDGEMENT

The research leading to these results has received funding from European Unions 7th Framework Programme under grant no. 318115, Connectivity management for eneRgy Optimised Wireless Dense networks (CROWD), http://www.ict-crowd.eu.