Luca Placidi
email: luca.placidi@uninettunouniversity.net

Francesco Dell'isola

Emilio Barchiesi

Heuristic Homogenization of Euler and Pantographic Beams

In the present contribution, we address the following problem: is it possible to find a microstructure producing, at the macro-level and under loads of the same order of magnitude, a beam which can be both extensible and flexible? Using an asymptotic expansion and rescaling suitably the involved stiffnesses, we prove that a pantographic microstructure does induce, at the macro-level, the aforementioned desired mechanical behavior. Thus, in an analogous fashion to that of variational asymptotic methods, and following a mathematical approach resembling that used by Piola, we have employed asymptotic expansions of kinematic descriptors directly into the postulated energy functional and a heuristic homogenization procedure is presented and applied to the cases of Euler and pantographic beams.

Introduction

While in the standard finite deformation Euler beam theory the energy functional depends only on the material curvature, i.e., the normalized projection of the second gradient of the placement on the normal vector to the current configuration, the energy functional for the nearly inextensible pantographic beam model depends also on the projection of the second gradient of the placement on the tangent vector to the current configuration. Thus, the full decomposition of the second gradient of the placement is present in the latter model. In order to analyze this fact, a heuristic homogenization

Introduction

Customarily, the theory of nonlinear beams is either postulated by means of a suitable least action principle in the so-called "direct way" or is deduced, by means of a more or less rigorous procedure, starting from a three-dimensional elasticity theory. The first example of direct model can be found in the original paper by Euler [START_REF] Euler | Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes sive solutio problematis isoperimetrici latissimo sensu accepti[END_REF]. Many epigones of Euler used this approach: a comprehensive account for this procedure can be found in Antman [START_REF] Antman | Nonlinear problems of elasticity[END_REF]. On the other hand, by following the procedure described by De Saint-Venant, one can try to identify the constitutive equation of a Euler type (1D) model in terms of the geometrical and mechanical properties, at microlevel, of the considered mechanical systems. This is done, in more modern textbooks, using a more or less standard asymptotic micromacro identification procedure, which generalizes the one used by De Saint-Venant for bodies with cylindrical shape (see, for instance, [START_REF] Placidi | An inverse method to get further analytical solutions for a class of metamaterials aimed to validate numerical integrations[END_REF]. It can be rigorously proven, under a series of well-precised assumptions, that only flexible and inextensible beams can be obtained [START_REF] Murat | Comportement asymptotique des solutions du système de l'élasticité linéarisée anisotrope hétérogène dans des cylindres minces[END_REF][START_REF] Mora | A nonlinear model for inextensible rods as a low energy γ-limit of three-dimensional nonlinear elasticity[END_REF][START_REF] Jamal | Théorie asymptotique des tiges courbes anisotropes[END_REF][START_REF] Pideri | Asymptotics of a non-planar rod in non-linear elasticity[END_REF][START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF].

Long fibers are often modeled as Euler beams. Here, we will define a Euler beam from a continuum point of view for the extensible and for the inextensible cases. A discrete model for the same beam will be also introduced and a heuristic homogenization procedure, see, e.g., dell 'Isola et al. (2016d), applied. A rescaling law will be derived for the extensible and for the inextensible cases.

Continuous Euler Beams

Kinematics

At each point S of C 0 , see Fig. 3.1, is associated a copy of the rigid section R through O such that C 0 and R are orthogonal. B 0 is the reference configuration of a beam. B is the present configuration, which is defined as (i) A vector function χ (S) that gives the present position of q 0 (S). (ii) An orthogonal tensor field R (S) that gives the rotation of R from the reference to the present configuration.

The kinematics is therefore defined by the following fields (Fig. 3.2):

χ (S) , R (S) .

(3.1)

The admissible motion is, e.g., for a cantilever, those kinematic fields such that χ (0) = 0, R (0) = 0.

Action

Physical intuition and definition of the action functional

A = L 0 W χ, R, χ , R -W ext (χ, R) -W ext S S=0,S=L ,
where W is the strain energy that is assumed to depend upon the fundamental kinematical fields and their derivative. W ext is the energy of the distributed forces and W ext S that of the concentrated ones.

Objectivity and Representation of the Invariants

Let us assume that the fundamental kinematical fields in one frame of reference are represented in (3.1). In another frame of reference they are as

χ = U + Qχ, R = QR, (3.2)
where U and Q are the translation and the rotation of the second frame of reference with respect to the first one. The derivative of (3.2) yields

χ = Qχ , R = QR . (3.3)
Let us define the following two fields in the first frame of reference

E = R T R (3.4) e = R T χ -q 0 (S) . (3.5)
In the second frame of reference, they are from (3.2) to (3.5)

Ẽ = R T R = R T Q T QR = R T R = E ẽ = R T χ -q 0 (S) = R T Q T Qχ -q 0 (S) = R T χ -q 0 (S) = e,
which means that they are invariant. In order to represent E and e, we assume a Cartesian frame of reference. The origin of such a frame of reference is q 0 (0)with basis

D 1 = q 0 , D 2 D 3 .
If S is a curvilinear abscissa of straight frame of reference q 0 = S D 1 that yields

q 0 = D 1 = ⎛ ⎝ 1 0 0 ⎞ ⎠ . The 2D assumption is χ ∈ span {D 1 , D 2 } .
We define the displacement vector field

u = χ -q 0 ,
and therefore its derivative

u = χ -D 1 ⇒ u = u 1 u 2 = χ 1 -1 χ 2 ,
so that the derivative of the placement is

χ = χ 1 χ 2 = 1 + u 1 u 2
and its squared modulus,

χ 2 = 1 + u 1 2 + u 2 2 .
A representation of R is given in terms of the rotation angle ϕ,

R = cos ϕD 1 ⊗ D 1 -sin ϕD 1 ⊗ D 2 + sin ϕD 2 ⊗ D 1 + cos ϕD 2 ⊗ D 2 (3.6) = cos ϕ -sin ϕ sin ϕ cos ϕ . (3.7)
Thus, the two invariants are represented as follows:

E = R T R = cos ϕ sin ϕ -sin ϕ cos ϕ -sin ϕ -cos ϕ cos ϕ -sin ϕ ϕ (3.8) = 0 -ϕ ϕ 0 = 0 -κ κ 0 e = ε γ = R T χ -q 0 (S) = (3.9) = cos ϕ sin ϕ -sin ϕ cos ϕ 1 + u 1 u 2 - 1 0 = 1 + u 1 cos ϕ + u 2 sin ϕ -1 -1 + u 1 sin ϕ + u 2 cos ϕ ,
which means in terms of κ, ε, and γ. A representation of the internal energy W that is compatible with the indifference frame principle is given by the function g

W χ, R, χ , R = g (κ, ε, γ) .
Let us give a representation for κ andε whether the beam is assumed to be shear un-deformable, i.e., with γ = 0. Thus, from the second equation of (3.9) we have

1 + u 1 sin ϕ = u 2 cos ϕ, ⇒ tan ϕ = u 2 1 + u 1 , (3.10) that means ϕ = arctan u 2 1 + u 1 = arctan χ 2 χ 1 . (3.11)
Besides from the first equation of (3.9), we have

ε = -1 + 1 + u 1 cos ϕ + u 2 sin ϕ (3.12) Keeping in mind that tan 2 ϕ = sin 2 ϕ cos 2 ϕ = 1 -cos 2 ϕ cos 2 ϕ = 1 cos 2 ϕ -1 that yields 1 + tan 2 ϕ = 1 cos 2 ϕ ,
we have from (3.10)

cos 2 ϕ = 1 1 + tan 2 ϕ = 1 1 + u 2 1+u 1 2 = 1 + u 1 2 1 + u 1 2 + u 2 2 , that yields cos ϕ = 1 + u 1 1 + u 1 2 + u 2 2 (3.13) and sin ϕ = tan ϕ cos ϕ = u 2 1 + u 1 1 + u 1 1 + u 1 2 + u 2 2 = u 2 1 + u 1 2 + u 2 2 . (3.14)
Therefore, from (3.12), (3.13), and (3.14)

1 + ε = 1 + u 1 1 + u 1 1 + u 1 2 + u 2 2 + u 2 u 2 1 + u 1 2 + u 2 2 = (3.15) = 1 + u 1 2 + u 2 2 1 + u 1 2 + u 2 2 = 1 + u 1 2 + u 2 2 = .χ (3.16)
Besides, the derivative of (3.11) is (3.17) where the 90 • rotation matrix is defined as follows:

κ = ϕ = 1 1 + χ 2 χ 1 2 χ 2 χ 1 -χ 2 χ 1 χ 2 1 = χ • * χ χ 2 ,
* = 0 -1 1 0 , so that χ • * χ = χ 1 χ 2 0 -1 1 0 χ 1 χ 2 (3.18) = χ 1 χ 2 -χ 2 χ 1 = -χ 2 χ 1 + χ 2 χ 1 . (3.19) Let us call ê = χ χ , c = χ χ . (3.20)
Thus, the curvature (3.17

) is κ = ϕ = * ê • c. (3.21)

Macroscopic Strain Energy for the General Case

A quadratic form of the strain energy in terms of the two invariants κ, from (3.21), and ε, from (3.16), is

U macro = l 0 1 2 K b κ 2 + 1 2 k e ε 2 dS (3.22) = l 0 ⎡ ⎣ 1 2 K b χ • * χ χ 2 2 + 1 2 K e χ -1 2 ⎤ ⎦ dS. (3.23)
It is worth to be noted that the strain energy is of second gradient type only for the normal component * ê • χ . The tangential component ê • χ do not have any contribution in the strain energy. In pantographic structures, we will see that also this tangential contribution is able to accumulate strain energy.

Macroscopic Strain Energy for the Inextensible Case

For inextensible beams, χ is the unit vector ê,

χ 2 = 1 = χ • χ (3.24) that means χ • χ + χ • χ = 0 =⇒ χ • χ = 0,
and

χ 2 = χ • χ = χ • * χ • χ * χ + χ • χ χ = = χ • * χ • χ * χ = χ • * χ 2 .
(3.25) Thus, from (3.24) and (3.25), the strain energy (3.23) for the inextensible case is

U macro = l 0 ⎡ ⎣ 1 2 K b χ • * χ χ 2 2 ⎤ ⎦ dS = l 0 1 2 K b χ 2 dS. (3.26)

Discrete Henky-Type Beam

Microscopic models for the inextensible Euler beams in the reference configuration are plotted in Fig. 3.3 (bottom). The bars are rigid and of length ε Piola's Ansatz

p i = χ (P i ) , (3.27)
Thus, the position of the other points is

i p i+1 = χ (P i+1 ) = χ (P i ) + εχ (P i ) + 1 2 ε 2 χ (P i ) .
(3.28)

p i-1 = χ (P i ) -εχ (P i ) + 1 2 ε 2 χ (P i ) , (3.29)
so that the cosine of the angle θ i is defined, in the inextensible case (3.24), by

cos θ i = (p i+1 -p i ) • (p i-1 -p i ) p i+1 -p i p i-1 -p i (3.30) = εχ (P i ) + 1 2 ε 2 χ (P i ) • -εχ (P i ) + 1 2 ε 2 χ (P i ) ε 2 , (3.31) = -ε 2 χ 2 + 1 2 ε 3 χ • χ -1 2 ε 3 χ • χ + 1 4 ε 4 χ • χ ε 2 (3.32) = -1 + 1 4 ε 2 χ • χ . (3.33)
Discrete energy in the inextensible case is defined as

U micro = N i=1 k b (1 + cos θ i ) = N i=1 1 4 k b ε 2 χ 2 .
With the homogenization formula,

N i=1 ε (•) -→ l 0 (•) dS, (3.34)
we have

U micro = l 0 1 4 k b ε χ 2 dS. (3.35)
An identification of (3.26) and (3.35) implies

k b ε = 2K b . (3.36)
Thus, in order to have a finite macro-energy, even in the limit ε → 0, we need to impose the scaling law (3.36). This means that, if one wants a finite macro-energy, then the lower the size of the cell, the higher is the rigidity of the rotational spring. Besides, in the limit of ε → 0, k b should be infinite, i.e., k b → ∞.

For extensible beams, the distance between internal hinges is not fixed to be equal to ε. The extensible bar at the right-hand side of P i has the following length:

p i+1 -p i = εχ (P i ) + 1 2 ε 2 χ (P i ) • εχ (P i ) + 1 2 ε 2 χ (P i ) (3.37) = ε 2 χ 2 + 1 2 ε 3 χ • χ + 1 2 ε 3 χ • χ + 1 4 ε 4 χ • χ (3.38) = ε χ 2 + εχ • χ + 1 4 ε 2 χ • χ = ε χ + ε 2 χ • χ 2 χ . (3.39)
The extensible bar on the left-hand side of P i has the following length:

p i-1 -p i = -εχ (P i ) + 1 2 ε 2 χ (P i ) • -εχ (P i ) + 1 2 ε 2 χ (P i ) = ε 2 χ 2 - 1 2 ε 3 χ • χ - 1 2 ε 3 χ • χ + 1 4 ε 4 χ • χ = ε χ 2 -εχ • χ + 1 4 ε 2 χ • χ = ε χ -ε 2 χ • χ 2 χ .
The representation of the cosine of the angle θ i in the extensible case, from (3.31) to (3.33), is

cos θ i = (p i+1 -p i ) • (p i-1 -p i ) p i+1 -p i p i-1 -p i (3.40) = εχ (P i ) + 1 2 ε 2 χ (P i ) • -εχ (P i ) + 1 2 ε 2 χ (P i ) ε χ + ε 2 χ •χ 2 χ ε χ -ε 2 χ •χ 2 χ (3.41) = -ε 2 χ 2 + 1 4 ε 4 χ • χ ε 2 χ 2 -1 4 ε 4 χ •χ χ 2 = -χ 2 + 1 4 ε 2 χ • χ χ 2 -1 4 ε 2 χ •χ χ 2 . (3.42)
The Taylor series expansion of the function f (x) around x = 0

f (x) = -a + bx a -cx = -1 + f (0) x, (3.43) f (x) = b (a -cx) -(-a + bx) (-c) (a -cx) 2 , (3.44) f (0) = ba -ac a 2 = b -c a (3.45) that imply cos θ i = -1 + χ • χ -χ •χ χ 2 4 χ 2 ε 2 . (3.46)
Discrete energy for the extensible case is therefore

U micro = N i=1 k e 2 ( p i+1 -p i -ε) 2 + k b (1 + cos θ i ) ,
that, because of (3.39) and (3.46)

U micro = N i=1 k e 2 ( p i+1 -p i -ε) 2 + k b (1 + cos θ i ) (3.47) = N i=1 k e 2 ε 2 χ -1 2 + k b χ • χ -χ •χ χ 2 4 χ 2 ε 2 . (3.48)
The last addend, because of the definitions (3.20), is rearranged as

χ • χ -χ •χ χ 2 χ 2 = c • c -c • ê 2 (3.49) = c • c -ê c • ê = c • c⊥ = c⊥ • c⊥ = c • * ê 2 , (3.50)
which imply another form of the discrete strain energy

U micro = N i=1 k e 2 ε 2 χ -1 2 + 1 4 k b χ • * χ χ 2 2 ε 2 . (3.51)
With the homogenization formula (3.34), the equation (3.51) yields

U micro = l 0 1 2 εk e χ -1 2 + 1 4 k b ε χ • * χ χ 2 2 dS.
Thus, in order to have a finite macro-energy, we need to impose the following rescaling law:

k b ε = 2K b , εk e = K e .
(3.52) Thus, in order to have a finite macro-energy, even in the limit ε → 0, we need to impose the scaling laws (3.52). This means that, if one wants a finite macro-energy, then the lower the size of the cell, the higher is the rigidity of the tensional and of the rotational springs. Besides, in the limit of ε → 0, k b should be infinite, i.e., k b → ∞.

Pantographic Beams

Introduction

In this section, we discuss the discrete micro-mechanical model which is employed throughout this paper. We begin giving a geometrical description and then we give a mechanical characterization, by choosing a deformation energy. It is a Hencky-type spring model with the geometrical arrangement of a pantographic strip. Once the energy of the micro-model is chosen in its general form, we assume a particular asymptotic behavior for some relevant kinematic quantities, i.e., the elongation of oblique springs, as will be clear in the sequel. We consider the quasi-inextensibility case, i.e., the relative elongation of the oblique springs is small. As a further specialization, the inextensibility case is considered. Finally, after having defined a micro-macro identification, we express the energy of the micro-system in terms of macroscopic kinematic descriptors to prepare the field to the homogenization procedure which will be discussed in details in the next section.

Discrete Micro-model

Geometry

In the spirit of dell'Isola et al. (2016d), [START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: A rigorous proof[END_REF], and [START_REF] Alibert | Extensional elastica in large deformation as -limit of a discrete 1d mechanical system[END_REF], in this section, we introduce a discrete-spring model (also referred to as the micro-model, since it resembles the features of a specific microstructure). The topology and features of the undeformed and deformed discrete-spring system are summarized in Figs. 3.4 and 3.5, respectively. In the undeformed configuration, N + 1 material particles are arranged upon a straight line at positions P i 's, i ∈ [0; N ], with a uniform spacing ε. The basic ith unit cell centered in P i is formed by four springs joined together by a hinge placed at P i . Between two oblique springs, belonging to the same cell and lying on the same diagonal, a rotational spring opposing to their relative rotation is placed. Rotational springs are colored in Fig. 3.4 in blue and red. We denote with p i the position in the deformed configuration corresponding to position P i in the reference one. In order to completely describe the kinematics of the micro-model, we have to introduce other descriptors. At this end, the length of the oblique deformed springs, indicated with l αβ i , is introduced, the indices α and β belonging, respectively, to the sets {1, 2} and {D, S} and referring to the first and second diagonal and left and right, respectively. Referring to Fig. 3.5, we consider the ith node, notwithstanding that the same quantities can be defined for each node. We define α i as the angle between the vectors p ip i-1 and p ip i+1 , respectively. We define as ϑ α i the angle measuring the deviation of two opposite oblique springs from being collinear. In order to illustrate the definition of ϕ αβ i , we consider the case α = 1 and β = D. The quantity ϕ 1D i is the angle between the vector p i+1p i and the upper oblique spring hinged at p i . By means of elementary geometric considerations, we have that

ϑ 1 i = α i + ϕ 1D i -ϕ 1S i ϑ 2 i = α i + ϕ 2S i -ϕ 2D i , i ∈ [0; N ] . (3.53)
In the undeformed configuration, see Fig. 3.4, we have

l αβ i = √ 2 2 ε, α = 1, 2 β = D, S i ∈ [0; N ] ϑ 1 i = ϑ 2 i = 0 p i -p i-1 = ε, i ∈ [0; N ] . (3.54) Considering that ϕ αD i , ϕ αS i ∈ [0, π],
by means of the law of cosines, we get

ϕ 1D i = cos -1 p i+1 -p i 2 + l 1D i 2 -l 2S i+1 2 2l 1D i p i+1 -p i ϕ 2D i = cos -1 p i+1 -p i 2 + l 2D i 2 -l 1S i+1 2 2l 2D i p i+1 -p i ϕ 1S i = cos -1 p i -p i-1 2 + l 1S i 2 -l 2D i-1 2 2l 1S i p i -p i-1 ϕ 2S i = cos -1 p i -p i-1 2 + l 2S i 2 -l 1D i-1 2 2l 2S i p i -p i-1 . (3.55)

Mechanical Model

The micro-model energy, written as a combination of the elastic energy contributions of the springs, is defined as

M = i α,β k e αβ,i 2 l αβ i - √ 2 2 ε 2 + i α k f α,i 2 ϑ α i 2 + + i k m i 2 ( p i+1 -p i -ε) 2 . (3.56) Reminding that ϑ α i = α i + (-1) α ϕ αS i -ϕ αD i
, then (3.56) recasts as

M = i α,β k e αβ,i 2 l αβ i - √ 2 2 ε 2 (3.57) + i α k f α,i 2 α i + (-1) α ϕ αS i -ϕ αD i 2 + + i k m i 2 ( p i+1 -p i -ε) 2 .
In the next subsections, we will specialize this form of the energy by means of assumptions on the properties of the micro-system. In particular, we will discuss in detail the representation of the micro-energy for the quasi-inextensibility assumption that will be made clear next and, subsequently, for the (complete) inextensibility cases.

Toward the Continuum Model

Asymptotic Expansion and Quasi-inextensibility Assumption

We postulate that the following asymptotic expansion holds for l αβ i :

l αβ i = ε lαβ i1 + ε 2 lαβ i2 + o ε 2 , (3.58)
where the constant (with respect to ε) term is not present. We now turn to what we refer to as the quasi-inextensibility case. It consists in fixing the value of the first-order term in (3.58) as

lαβ i1 = √ 2
2 . Moreover, to lighten the notation, we drop the subscript "2" of lαβ i2 , i.e., lαβ i = lαβ i2 . Hence, (3.58) reads as

l αβ i = √ 2 2 ε + ε 2 lαβ i + o ε 2 . (3.59)

Piola's Ansatz

The reference shape of the macro-model is a one-dimensional straight segment S and we introduce on it an abscissa s ∈ [0, B] -where B = N ε is the length of S which labels each position in S. Proceeding as in the pioneering works of Gabrio Piola, an Italian mathematician and physicist who lived in the 1800s (see Dell'Isola et al. 2015 for a historical review), we introduce the so-called kinematical maps, i.e., some fields in the macro-model that uniquely determine p i and lαβ i :

χ : [0, B] → E lαβ : [0, B] → R + , (3.60)
with E the Euclidean space on V ≡ R 2 . We choose χ to be the placement function of the 1D continuum and, hence, it has to be injective. The current shape can be regarded as the image of the (sufficiently smooth) curve χ : [0, B] → E and, unlike the reference shape, it is not parameterized by its arc length and it is not a straight line in general. In order for these fields to uniquely determine the kinematical descriptors of the micro-model (i.e., p i and lαβ i ), we use the Piola's Ansatz and impose

χ (s i ) = p i lαβ (s i ) = lαβ i , ∀i ∈ [0; N ] . (3.61)

Micro-model Energy as a Function of Macro-model Descriptors

In this subsection, we obtain the micro-model energy for the quasi-inextensibility case in terms of the macroscopic kinematical maps. Assuming that χ is at least twice continuously differentiable with respect to the space variable in s i 's, we have

χ (s i+1 ) = χ (s i ) + εχ (s i ) + ε 2 2 χ (s i ) + o ε 2 χ (s i-1 ) = χ (s i ) -εχ (s i ) + ε 2 2 χ (s i ) + o ε 2 . (3.62)
Plugging (3.61) in (3.59) and (3.62), we get the following expressions:

l αβ i = √ 2 2 ε + ε 2 lαβ (s i ) + o ε 2 p i+1 -p i = εχ (s i ) + ε 2 2 χ (s i ) + o ε 2 p i-1 -p i = -εχ (s i ) + ε 2 2 χ (s i ) + o ε 2 . (3.63)
Substituting (3.63) into (3.55) and expanding ϕ αS i -ϕ αD i up to first order with respect to ε, we get

ϕ αS i -ϕ αD i = √ 2 4 χ (s i ) 2 + l(3-α)D (s i-1 ) -l(3-α)S (s i+1 ) χ (s i ) 1 -χ (si) 2 2 ε + + χ (s i ) 2 -1 lαS (s i ) -lαD (s i ) χ (s i ) 1 -χ (si) 2 2 ε + o (ε) . (3.64)
Finally, substituting (3.64) in (3.57) yields the micro-model energy M as a function of the kinematical descriptors χ and lαβ of the macro-model

M = i α,β k e αβ,i ε 4 2 lαβ i 2 + i k m i ε 2 2 χ i -1 2 (3.65) + i α k f α,i ε 2 2 ϑ (s i ) + (-1) α √ 2 4 χ (s i ) 2 + l(3-α)D i (s i-1 ) -li (3-α)S (s i+1 ) χ (s i ) 1 -χ (si) 2 2 + + (-1) α χ (s i ) 2 -1 lαS i (s i ) -li αD (s i ) χ (s i ) 1 -χ (si) 2 2 2
, where α i = εϑ (s i ) has been used and

ϑ = χ ⊥ • χ χ 2 ,
with χ ⊥ the 90 • anticlockwise rotation of χ , is the material curvature, i.e., rate of change with respect to the reference abscissa of the orientation of the tangent χ (s) = ρ (s) [cos ϑ (s) e 1 + sin ϑ (s) e 2 ] to the deformed centerline. We remark that the micro-model energy, when written in terms of macroscopic fields, contains already a contribution from the second gradient of χ(s). Finally, it is worth to be noticed that, for a fixed ε, Eq.(3.65) provides an upper bound for ||χ ||, i.e., χ < √ 2, even if no kinematic restrictions directly affect ||χ ||.

The Case of Inextensible Fibers

We consider now the case of inextensible oblique springs. This translates in considering lαβ i = 0 and it is referred as the inextensibility case. Moreover, for the sake of simplicity, we consider the elastic constants of the rotational springs to satisfy k 1

F,i = k 2 F,i := k F,i , ∀i ∈ [1; N ]. We remark that lαβ i = 0 implies, through a purely geometric argument, that ϕ S1 i+1 = ϕ S2 i+1 = ϕ D1 i = ϕ D2 i := ϕ i .
Once the kinematic restrictions implied by the inextensibility assumption have been presented, we are ready to define the micro-model energy (3.57) as

I = i k f i α α i + (-1) α (ϕ i -ϕ i-1 ) 2 2 (3.66) + i k m i 2 ( p i+1 -p i -ε) 2 .
Proceeding in analogy with the previous construction, we introduce the kinematical map

ϕ : [0, B] → 0, π 2
and, then, we perform the Piola's Ansatz by imposing

ϕ (s i ) = ϕ i , ∀i ∈ [0; N ] . (3.67)
Assuming both χ and ϕ to be at least one time continuously differentiable with respect to the space variable in s i and taking into account the Piola's Ansatz (3.67), we have

p i+1 -p i = εχ (s i ) + o(ε) ϕ i-1 -ϕ i = -εϕ (s i ) + o(ε). (3.68)
Substituting (3.68) into (3.66) yields the micro-model energy for the inextensibility case in terms of the kinematical quantities of the macro-model

I = i k f i ε 2 ϑ 2 (s i ) + ϕ i 2 (s i ) + i k m i ε 2 2 χ i -1 2 . (3.69)
We now impose the so-called internal connection constraint:

√ 2ε cos ϕ (s i ) = χ (s i+1 ) -χ (s i ) , (3.70)
which, up to ε-terms of order higher than one, reads √ 2 cos ϕ = χ .

(3.71)

This constraint ensures that, in the deformed configuration, the upper-left spring of the ith cell is hinge joint with the upper-right spring of the (i -1)th cell, and the lower-left spring of the ith cell is hinge joint with lower-right spring of the (i -1)th cell. Due to this constraint, the maps ϕ and χ are not independent and it is possible to rewrite the expression of the micro-model energy in terms of the placement field χ(s) only. Indeed, deriving (3.71) with respect to the space variable yields

- √ 2ϕ (s i ) sin ϕ (s i ) = χ (s i ) , (3.72)
which, in turn, implies

ϕ (s i ) = - χ (s i ) √ 2 sin ϕ (s i ) .
Reminding ϕ ∈ [0, π] and taking into account (3.71), we get

ϕ i = - χ i √ 2 1 -cos 2 ϕ (s i ) = = - χ i 2 -χ i 2 .
Hence, in the inextensibility case, the micro-model energy (3.69) can be recast, as a function of the macro-model descriptor χ only, as

I = (3.73) i k f i ε 2 ⎡ ⎣ ϑ (s i ) 2 + χ (s i ) 2 -χ (s i ) 2 2 ⎤ ⎦ + i k m i ε 2 2 χ (s i ) -1 2 .
Clearly, since the inextensibility case is just a special case of the quasi-inextensibility case, it is possible to show that this expression can be also obtained in a more direct way from (3.65) by setting lαS

(s i ) = 0 and k 1 F,i = k 2 F,i := k F,i .

Continuum-Limit Macro-model

In this section, by performing the final steps of the heuristic homogenization procedure presented throughout this paper, we derive a 1D continuum model, also referred to as the macro-model, associated to the aforementioned microstructure. Besides, we analyze the quasi-inextensibility and inextensibility cases and we obtain the corresponding macro-model energies in terms of the displacement field χ.

Rescaling of Stiffnesses and Heuristic Homogenization

The preliminary step to perform the homogenization procedure consists of the definition of the quantities K e αβ,i , K f α,i , and K m i . These quantities are scale-invariant, meaning that they do not depend on ε. Their role is to keep track of the asymptotic behavior of the stiffnesses k e α,β,i , k f α,i , and k m i of the micro-model springs. More explicitly, we assume

k e αβ,i (ε) = K e αβ,i ε 3 ; k f α,i (ε) = K f α,i ε ; k m i (ε) = K m i ε . (3.74)
We remark that in this rescaling, as ε approaches zero, the ratio between the stiffness k e αβ,i of the oblique springs and the stiffness k f α,i will approach infinity with a rate of divergence in ε equal to two, i.e.,

k e αβ,i k e α,i
∼ ε 2 . Now, we are ready to perform the homogenization procedure. First, we consider the more general quasi-inextensibility case. For simplicity, let us set

K e 1D,i = K e 1S,i = K e 2D,i = K e 2S,i := K e i ; K f 1,i = K f 2,i := K f i . (3.75)
Let us introduce the kinematical maps

K e : [0, B] → R + ; K f : [0, B] → R + ; K m : [0, B] → R +
such that they satisfy the following Piola's Ansatz:

K e (s i ) = K e i ; K f (s i ) = K f i ; K m (s i ) = K m i . ( 3
+ S K e 2 l2S 2 ds + S K e 2 l2D 2 ds + + S K f 2 ϑ + - √ 2 χ 2 -4 l2D -l2S -χ 2 -1 l1D -l1S χ 2 -χ 2 ⎫ ⎬ ⎭ 2 ds + + S K f 2 ϑ + √ 2 χ 2 + 4 l1D -l1S + χ 2 -1 l2S -l2D χ 2 -χ 2 ⎫ ⎬ ⎭ 2 ds + + S K m 2 χ -1 2 ds, (3.77)
which is the continuum-limit macro-model energy for a 1D pantographic beam under the hypothesis of quasi-inextensible oblique micro-springs. It is worth to remark that, when K m = 0, lαβ = 0 and χ (s) = Cse 1 , with C ∈ R, the beam undergoes a floppy mode, i.e., (3.77) vanishes. Thus, under the above conditions, the configuration χ (s) = Cse 1 is isoenergetic to the undeformed configuration for any C. For a fixed ε, considering k m i = 0 and li αβ = 0 in the micro-model energy (3.65), we have that χ(s i ) = Cs i e 1 is a floppy mode for the micro-model as well. This means that the homogenization procedure that we have carried out has preserved a key feature of the micro-model. Up to now, the expression of the continuum-limit homogenized energy depends both on the kinematical maps χ and l. In the next section, we show that, at equilibrium, it is possible to write the macro-energy in terms of the placement field only.

Macro-model Energy as a Function of the Placement Field

We now equate to zero the first variations of (3.77) with respect to lαβ 's, i.e., we look for stationary points, with respect to lαβ , of (3.77). This is a necessary first-order condition for optimality. In the continuum limit homogenized energy, no spatial derivatives of lαβ appear. Such energy depends only by linear and quadratic contributions in lαβ . Hence, this process yields four algebraic linear equations in lαβ . Solving these equations gives lαβ at equilibrium

l1D = √ 2 2 K f χ • C + ϑ D l2D = √ 2 2 K f χ • C -ϑ D l1S = √ 2 2 K f -χ • C -ϑ D l2S = √ 2 2 K f -χ • C + ϑ D (3.78) with C = χ 2K f χ 2 -1 2 K e χ 2 + 8K f D = χ 4 L2 -χ 2 K e L2 χ 2 -2 -2K f χ 2 .
From (3.78), we can get, in some particular cases, interesting information about the properties of the pantographic beam. First, let us notice that l1D = -l1S and l2D = -l2S . Moreover, we also notice that when χ = ρe 1 , with ρ independent of the abscissa s, then, as χ vanishes, lαβ = 0, i.e., the fibers undergo no elongation. Instead, when χ (s) = ρ (s) e 1 , with ρ depending on s, then l1D = l2D = -l1S = -l2S . This remarkable and counterintuitive feature can be used as a possible benchmark test to validate, as ε approaches zero, a numerical scheme based on the discrete micro-model. Let us consider the case of nonzero bending curvature, i.e., ϑ = 0, when χ • C << ϑ D, which implies that l1D = -l2D = -l1S = l2S . If ϑ > 0 then l1D , l2S > 0 and l2D , l1S < 0 while, if ϑ < 0 then l1D , l2S < 0 and l2D , l1S > 0. We are now ready to express the macro-model energy E (χ) as a function of the placement χ only, by substituting (3.78) in (3.77):

E (χ (•)) = min lαβ (•) E = (3.79) S K e K f ρ 2 -2 ρ 2 K e -4K f -2K e ϑ 2 + ρ 2 2 -ρ 2 ρ 2 K e -4K f + 8K f ρ 2 ds + + S K m 2 (ρ -1) 2 ds = = S K e K f χ 2 -2 χ 4 χ 2 K e -4K f -2K e χ ⊥ • χ 2 ds + + S K e K f 2 -χ 2 χ 2 K e -4K f + 8K f χ • χ 2 ds + + S K m 2 χ -1 2 ds.
We observe that, for 0 < ρ < √ 2 and for any choice of the positive macro-stiffnesses K e , K f , and K m , (3.79) is positive definite. Moreover, not only we can classify this homogenized model as a second gradient theory, but we notice that the full second gradient χ of χ contributes to the strain energy. Indeed, beyond the usual term χ ⊥ • χ related to the Lagrangian curvature, also the term χ • χ , deriving from the presence of the oblique springs, appears. There is a remarkable feature in this model which deserves to be discussed. From (3.79), it is clear that in the limit ||χ || → √ 2 the model exhibits a so-called phase transition: it locally degenerates into the model of an uniformly extensible cable, notwithstanding that √ 2 is an upper bound for ρ. Indeed,

ρ 2 -2 ρ 2 K e -4K f -2K e → 0 ρ 2 2 -ρ 2 ρ 2 K e -4K f + 8K f → +∞,
so that no deformation energy is stored for finite bending curvature and, in order for the energy to be bounded for bounded deformations, ρ must approach zero, meaning that the elongation must be locally uniform. Further developments of this model could consist in contemplating a phase transition to a model that, for finite bending curvature, entails a nonzero stored deformation energy.

Nondimensionalization

In order to handle more easily the model in the numerical implementation and in the interpretation of the corresponding results, we turn to the use of nondimensional quantities. Therefore, we introduce the following nondimensional fields:

s = Bs; χ = Bχ; K e = KK e ; K f = KK f ; K m = K m K m .
In terms of these new quantities, we can recast (3.79) as (3.80) where the symbol denotes differentiation with respect to the dimensionless abscissa s.

K B 1 0 K e K f χ 2 -2 χ 4 χ 2 K e -4K f -2K e χ ⊥ • χ 2 ds + + K B 1 0 K e K f χ • χ 2 2 -χ 2 χ 2 K e -4K f + 8K f ds + +K m B 1 0 K m 2 χ -1 2 ds,

The Inextensibility Case

Let us focus now on the inextensibility case. The homogenization procedure follows the same lines of the previous case. Indeed, keeping in mind (3.75) and (3.76), letting ε → 0 in (3.73) yields the continuum-limit macro-model energy for the inextensibility case

S K f ϑ 2 + ρ 2 2 -ρ 2 + K m 2 (ρ -1) 2 ds = = S K f χ ⊥ • χ 2 χ 4 + χ • χ 2 χ 2 2 -χ 2 + K m 2 χ -1 2 ds. (3.81)
This result is consistent with the quasi-inextensibility case. Indeed, we could have found (3.81) also by letting K e → +∞ in (3.79). Let us remark that, also in this case, the homogenized continuum model, due to the richness of the microstructure, gives rise to a full second gradient theory.

Linearization

An interesting connection can be traced with the existing literature on the formulation of 1D continuum homogenized model for microstructured media and, in particular, for pantographic ones. Indeed, this connection is traced by considering a linearization of the pantographic beam energy in the (complete) inextensibility case. We set χ (s) = s 0 + η ũ, with ũ independent of η, i.e., we linearize with respect to the displacement u = χ (s) -s 0 , and K m = 0. By means of simple algebra manipulations, it is possible to derive the deformation energy in Eq. (with K + = K -) of [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF] (see also [START_REF] Seppecher | Linear elastic trusses leading to continua with exotic mechanical interactions[END_REF]:

S K f u 2 ds. (3.82)
We remark that in the linearized energy (3.82) the transverse displacement and the axial one decouple.

Numerical Simulations of the Continuous Model

Preliminaries

Using the so-obtained 1D continuum model, we show some equilibrium shapes exhibiting highly non-standard features, essentially related to the complete dependence of the homogenized continuum energy density functional on the second gradient of the placement field.

In the sequel, K m = 0 will be considered, which means that the standard quadratic additive elongation/shortening contribution to the deformation energy will be turned off. This is made in order to better highlight some non-standard features of the nearly inextensible pantographic beam model. In this section, we show numerical results for the quasi-inextensible and inextensible pantographic beam model and for the geometrically nonlinear Euler model. We remind that these cases stand for K e < +∞ and K e → +∞, respectively. Two benchmark tests are exploited in order to illustrate peculiar and non-standard features of the pantographic beam model. Convergence of the quasi-inextensible pantographic beam model to the completely inextensible one is shown, by means of a numerical example, as the macro-stiffness K e related to elongation of the oblique springs approaches +∞. This is due to the fact that, as it is clear from Eq. (3.77), if K e → +∞, then lαβ → 0. Of course, the same discussion and simulations can be made for the micro-model and this could be the subject of a further investigation. For the sake of self-consistence, we recall that the deformation energy of the geometrically nonlinear Euler model employed in the following simulations is the following:

S K e 2 χ -1 2 + K b 2 χ • χ χ 2 - χ • χ χ 2 2 ds = = S K e 2 (ρ -1) 2 + K b 2 ϑ 2 ds,
and we notice that, while in the nearly inextensible pantographic beam model both ρ and ϑ can be enforced at the boundary, for the nonlinear Euler model it can be done for ϑ only, as no spatial derivative of ρ appears in the energy.

Semi-circle Test

We consider for both the nearly inextensible pantographic beam model and the geometrically nonlinear Euler beam model the reference domain to be the interval [0, 2π].

We enforce the following boundary conditions for both models

1. χ (0) = 0; 2. χ (2π) = 2e 1; 3. ϑ (0) = - π 2 ; 4. ϑ (2π) = π 2
and, for the nearly inextensible pantographic beam model, we also have the following two additional constraints:

5. ρ (0) = ρ 0 ; 6. ρ (2π) = ρ 0 .
In

Three-Point Test

We consider for both the quasi-inextensible pantographic beam model and the geometrically nonlinear Euler beam model the reference domain to be the interval [0, 2]. We enforce the following boundary conditions for both models:

1. χ (0) = 0; 2. χ (1) • e 2 = u; 3. χ (2) = 0; 4. ϑ (0) = 0; 5. ϑ (2) = 0.

In Fig. 3.8, the deformed shapes for the nearly inextensible pantographic beam model (red, light blue) and for the geometrically nonlinear Euler beam model (blue, green) are shown for different values of u in the legend. Figure 3.9 shows, for different values of the parameter u, the elongation ρ -1 versus the reference abscissa for the nearly inextensible pantographic beam model. The parameter u is increasing from bottom to top. We observe that, as u increases, at some point, there is a concavity change in the elongation plot and, increasing further the parameter u, curves start to intersect. This means that, for some points of the beam, an increase of the prescribed displacement u implies a decrease in the elongation. Figure 3.10 shows the pulling force, i.e., Lagrange multiplier associated to the weak constraint χ (1) • e 2 = u, changed of sign, applied at the midpoint in order to vertically displace it of an amount u. In the nearly inextensible pantographic beam model (blue) negative stiffness property, also Figure 3.12 shows the deformed configuration for the nearly inextensible pantographic beam model, while in Fig. 3.13 the elongation ρ -1 versus the reference abscissa for the nearly inextensible pantographic beam model is shown.

Fig. 3 . 1

 31 Fig. 3.1 Reference configuration B 0 . q 0 (S) is the position of the origin O of the section R

Fig. 3 . 3

 33 Fig. 3.3 Microscopic models for the inextensible Euler beams in (bottom) the reference configuration. Definition (top) of the angle θ i

Fig. 3 . 4

 34 Fig. 3.4 Undeformed spring system resembling the microstructure

  Fig. 3.6 (up), the deformed shapes for the nearly inextensible pantographic beam model and for the geometrically nonlinear Euler beam model (GNEM) are shown for different values of ρ 0 reported in the legend. In Fig. 3.6 (down), the elongation ρ -1 for the nearly inextensible pantographic beam model and for the geometrically nonlinear Euler beam model (GNEM) is shown for different values of ρ 0 reported in the legend. It is remarkable that passing from ρ 0 > 1 to ρ < 1, there is a change

Fig. 3 . 6

 36 Fig. 3.6 Semi-circle test. Deformed shapes for the nearly inextensible pantographic beam model and for the geometrically nonlinear Euler beam model (GNEM). (left) Elongation ρ -1 versus the reference abscissa for the nearly inextensible pantographic beam model and for the geometrically nonlinear Euler beam model (GNEM) (right). Numbers in the legends stand for different dimensionless values of ρ 0

Fig. 3 . 8

 38 Fig. 3.8 Three-point test. Deformed shapes for the nearly inextensible pantographic beam model (red, light blue) and for the geometrically nonlinear Euler beam model (blue, green) for different values of u in the legend

Fig. 3 .

 3 Fig. 3.10 Three-point test. Pulling force (i.e., Lagrange multiplier associated to the weak constraint χ (1)• e 2 = u), changed of sign, applied at the midpoint in order to vertically displace it of an amount u (abscissa). In the nearly inextensible pantographic beam model (blue) elastic softening is observed, while in the geometrically nonlinear beam model (green) elastic softening is not observed

Modified Three-Point Test

We consider for both the quasi-inextensible pantographic beam model and the geometrically nonlinear Euler beam model the reference domain to be the interval [0, 2]. We enforce the three-point test boundary conditions for both models 1. χ (0) = 0; 2. χ (1) • e 2 = u; 3. χ (2) = 0; 4. ϑ (0) = 0; 5. ϑ (2) = 0


domain: spi phys

text:
Paul Chapellier
email: paul.chapellier@onera.fr

Pierre Lavenus

Olivier Le Traon

Bernard Dulmet

High Q 2D-Length Extension Mode Resonators for potential Time-Frequency applications

This paper presents recent advances on two dimensional Length-Extension Mode (2D-LEM) quartz resonators providing high quality (Q) factor on resonances at a few MHz. The resonators have been collectively manufactured using one or two steps quartz Deep Reactive-Ion Etching (DRIE) processes. These resonators combine the intrinsic qualities of quartz in comparison to silicon (i.e. high Q factor, low temperature sensitivity and piezoelectricity) and the advantages of Microelectromechanical systems (MEMS) resonators: small dimensions, low power consumption and collective processes. Samples vibrating at frequencies f of 2.2, 3 and 4.5 MHz have shown promising results with very high Q factor. Q factor as high as 180 000 for fundamental mode vibrating at 2.2 MHz and 89 000 for harmonic mode at 8.9 MHz were measured which lead to quality factor and resonance frequency products (Q.f) figure of merit near 10 12 Hz at the state of the art for 2D-LEM quartz resonators and the higher Q factor measured for DRIE made quartz resonators. Two designs, several dimensions and two processes have been investigated. Two main limiting damping mechanisms were identified and one of them is strongly linked to the technological limits of the etching process.

INTRODUCTION

An oscillator is a self-sutaining system capable of generating a signal with a stable frequency. Oscillators are used as frequency references, frequency synthetiser or clocks in almost every electronic devices. They are mandatory in a wide range of fields:watchs, consumer electronics, network equipment, automotive, or embedded systems for space and military applications (Vig2016) An oscillator is composed at least of a frequency selective element, usually an electromechanical resonator that acts like a filter and a gaining element that maintains oscillations.

A resonator is primaly characterized by its resonance frequency in Hz, its quality (Q) factor that reflects the sharpness of the resonance and its serie resistance Rm in Ω. For oscillator applications, high Q factor is mandatory to ensure low phase noise close to the carrier and low Rm is important to ensure the start of the oscillations and low phase noise far from the carrier. Other parameters such as frequency stability with temperature, with time (short term or ageing) or pressure are also important to ensure a stable output frequency. The multibillion dollar market of oscillators is dominated by quartz based oscillators that benefit from significant advantages of quartz such as low internal viscous losses (resulting in high Q factor), linear conversion between mechanical and electrical energy through piezoelectricityand high stability versus temperature for well-chosen crystal cuts.

However quartz resonators can be quite bulky and their integration is off-chip leading to a cumbersome packaging and they are somewhat incompatible with emerging applications such as mobile, new telecommunication networks or Internet of Things (IoT) where small size, low cost, low power consumption and medium to high performances are required (Nguyen2007). To overcome these issues silicon (Si) or Thin-film piezoelectric-on-silicon (TPoS) Micro-Electro-Mechanical Sytems resonators have been developed since the end of the 90' (Jokic2015) Altought Si MEMS resonators can show very high Q factor up to several millions at few MHz (Xereas2015), thermal frequency stability of silicon is pretty poor (around 30 ppm/°C) and has to be compensated using additionnal material layer or doping. Capacitive transduction leads to high serie resistance (at best few kΩ versus 50 to 100 Ω for standard quartz resonators) and induced non-linearities lead to a degradation of the phase noise close from the carrier (Kaajakari2005). TPoS MEMS resonators using AlN or ZnO on Si bulk show low Q factors,lower than 10 000 usually (vanBeek2012). and thermal stability of the stack also needs compensation. On the other hand Rm is low (from tens of Ω up to 1 kΩ) thanks to piezoelectricity. The approach for the Si and TPoS MEMS is to use established MEMS processes on well known but not the most suited materials for oscillator applications.

The approach developed in this paper is different as the aim is to use quartz-one of the best if not the best material for oscillator application-with collective processes still under development to manufacture higly miniaturized quartz resonators in order to integrate them in high stability, low cost and power consumption oscillators. As the dimensions are shrinking the standard wet chemical etching using Hydrofluoric acid (HF) and Ammonium Fluoride (NH4F) is not suitable because of etching induced faceting [START_REF] Pelle | Etch rates of crystallographic planes in Z -cut quartz-experiments and simulations[END_REF] nor is the chemical-mechanical polishing. This paper deals with the designs of new planar quartz resonators for oscillator applications based on a previous high Q factor bulky resonator originally designed to reach the quantum ground state. The challenging collective processes flow using quartz Deep Reactive Ion Etching (DRIE), the characterizations of the resonators and the comparison with the expected performances are furthermore discussed.

PREVIOUS DESIGN AND MOTIVATIONS

Reaching the quantum ground state (QGS) with a macroscopic device thus creating a bridge between quatum and classic theories is a topic that many teams and laboratories have been interested in within the last two decades. Devices used are typically nano-micromechanical resonators such as nanotube resonator (Khosla2019), aluminum micromechanical disk resonator (Pirkkalainen2015), silica toroidal resonators (Rivière2010) or aluminu nitride microresonator (Connell2010). Advanced cooling techniques are required to reach such low energy level. Extremely high sensitivity are mandatory to observe quantum fluctuations as well as very high resonators quality factors to ensure long coherence time and long state life time. In a previous work of our group (Le Traon et al 2011), a three dimensional Length-Extension Mode (3D-LEM) quartz resonator was designed, fabricated and characterized to reach the ground state. Figure 1a) shows a cross-section of the 3D-LEM concept for a cylindrical structure: a pillar is vibrating in a length-extensional motion whereas a ring around this pillar is vibrating in an out-of-phase motion allows to compensate the lateral displacement induced by Poisson effect in tethers. This compensating effect leads to reduced displacements and energy losses at the anchors areas. Figure 1b) shows a realization of nine 3D-LEM resonators on a 1.5-inches square Z-cut quartz wafer using wet chemical etching and Figure 1c) shows a detail of the micro-pillar. a triangular pillar was selected to reduce as much as possible the faceting induced by wet chemical etching. This 3D-LEM quartz resonator exhibited a very high Q factor of 2 000 000 at 4 MHz reaching a Q.f product of 8x10 12 Hz near the theoritical limit of quartz of 3.2x10 13 Hz in the Akhiezer regime where the predominant losses mechanisms is phonon-phonon interaction (Akhiezer 1939). Thanks to its very high Q factor, this resonator might find a relevant application in another field such as Time and Frequency where very high Q factor quartz resonator are usually used to control oscillators for timing and frequency references. The best resonators used in Ultra Stable Oscillators (USO) are usually bulky (>200 mm 3 ) shear mode Stress Compensated (SC) or Temperature Compensated (AT) cut quartz resonators with very high Q up to 2.6x10 6 at 5 MHz (Besson 1995) to ensure high frequency stability. The 3D-LEM resonator performances are close to that and that is why this resonator could be a potential candidate for Time and Frequency applications. Unfortunately, the 3D-LEM resonator is also quite bulky (>100 mm 3 ) and its 3D design does not allow electrical actuation. Compact planar two dimensional length-extension mode (2D-LEM) resonators using similar working principle have been designed to overcome this issues and should be suited for the realization of quartz micro-resonator based MEMS oscillators (Bourgeteau2014).

NEW RESONATORS DESIGNS AND QUALITY FACTOR

3.1

Mode of interest 2D-LEM resonators of type I and II consist of a beam vibrating in extension-compression. The resonance frequency is around a few MHz with a central beam of few millimeters. More precisely, the resonance of the fundamental mode (n=0) and harmonics for Z-cut quartz are given by the well-known equation at the first order (1):

• !"# = !"# $ % # & '' ( = (2* + 1) .•• !"# (1)
where L is the total length of the vibrating beam in mm, $ is the mode number, % the quartz density and & •• the elastic compliance along Y-axis. As seen from ( 1), the resonance frequency can be adjusted depending on the targeted applications by varying L. This should be done though homothetic transformations of the whole structure to maintain the optimized decoupling between the beam and the mounting area. For the ease of processing, resonators are actuated via gold electrodes only deposited on the upper surface.

Thanks to the quartz piezoelectric coupling an electric field applied along X-axis induces a strain along Y-axis as it can be seen in Figure 2.

Resonators designs

3.2.1

First design, "type I" resonator In addition to the central beam, the type I resonator has two large and short lateral beams that vibrate in an in-plane forced flexure movement that compensates for the shear stress induced through Poisson effect by the vibrating central beam. Figure 3a shows the amplified nodal displacements of type I and II resonators using the Finite Element Method (FEM) simulations with the Multiphysics software OOFELIE and the graphical interface SAMCEF. Thus, the confinement of the stress close to the central beam allows to minimize the vibrating energy in the mounting area and reduces the energy loss towards the mounting which is absolutely necessary to reach the expected high intrinsic quality factor of the quartz resonator. Figure 3b shows the X component of the magnified deformation of the beam in a flexure motion and Figure 3c shows the Y component of the deformation. Typical nominal dimensions and targeted frequencies of the type I resonator are shown in table 1. Table 1 Mean dimensions of type I resonator and targeted frequencies

3.2.2

Second design, "type II" resonator In addition to the central beam, the type II resonator has two thin and long lateral beams as shown at the top of Figure 4. The decoupling between the central beam and the mounting area is partly ensured by a flexure movement similar to the type I resonator. The type II 2D-LEM resonator also benefits from an additional compensation mechanism relying on a length-extension vibration of the lateral beams in opposite phase with respect to the central beam [START_REF] Patent | Structure plane de résonateur mécanique découplé par des vibrations de flexion et d'extension compression[END_REF]. This is illustrated in Figure 4b and Figure 4c which show respectively the magnified deformations along X and Y-axis of the type II 2D-LEM resonator. A drawback of this additional decoupling structure is the addition of unwanted eigen modes close to the natural resonance of the central beam corresponding to the eigen modes of the lateral beams which would result in mode coupling if the quality factor is not high enough. Typical nominal dimensions and targeted frequencies of the type II resonator are shown in table 2.

3.3

Theoritical Q factor The Q factor of a resonator can be defined as the ratio of energy stored within the structure over the energy dissipated per cycle via all the different damping mechanisms. Assuming no coupling between the damping mechanisms, energy losses can be written as the sum of the loss induced by the mechanism considered separately, leading to a quality factor that can be written under vacuum as:

! " = ! " #$%& + ! " '() + ! " *), + ! " -.&/01 (2)
where damping mechanisms are briefly described further. A more complete description can be found in [START_REF] Brand | Resonant MEMS-Fundamentals, Implementation and Application[END_REF] and an application in the case of silicon beams is provifed in (Le Foulgoc.et al. 2006). Qanchor corresponds to the strain energy lost at the anchors in the mounting areas and is estimated using FEM simulations. By design Qanchor is very high (>10 8 ) for a flawless structure. QTED is estimated using FEM simulations (Oofelie Solver) able to take into account the thermoelasticity damping. This loss mechanism is due to the bending of the beams which induces a thermal gradient and a strain gradient leading to thermal transport and damping. In the case of 2D-LEM, the flexure motion is not predominant and QTED is very high (>10 7 ). Akhiezer loss QAKE (Akhiezer 1939) also known as the quantum limit corresponds to the energy loss by phonon-phonon interaction needed to return in a thermodynamics equilibrium due to the lattice deformation induced by the displacement. QAKE is estimated via an analytic expression [START_REF] Ghaffari | Quantum limit of quality factor in silicon micro and nano mechanical resonators[END_REF]) and the Q.f product depends on the material. For quartz Q.f reaches 3.2x10 13 Hz leading to a QAKE around 10 7 for structures vibrating at few MHz. Finally, energy loss arises from the viscosity of the gold electrodes which acts as a damping layer. The related Qvisc can be estimated using an analytical model developed in [START_REF] Le Traon | LGS and GaPO4 piezoelectric crystals: New results[END_REF]) and [START_REF] Bourgeteau | Quartz resonator for MEMS oscillator[END_REF]) is expected to be the main damping contribution in the 2D-LEM case. With a perfect geometry, a quality factor around 500 000 is expected for resonators vibrating at a few MHz making possible to reach a Q.f product between 1.1 and 1.5x10 12 Hz.

PROCESS FLOW

Standard wet etching using Hydrofluoric acid (HF) and Ammonium Fluoride (NH4F) is not suitable for etching the vibrating part due to its high anisotropy and faceting [START_REF] Pelle | Etch rates of crystallographic planes in Z -cut quartz-experiments and simulations[END_REF]. Deep Reactive Ion Etching of quartz has been developed for the last decade thanks to the development of high density plasma etcher. DRIE-made resonators have been reported for various applications such as Quartz Crystal Microbalance [START_REF] Kutsuwada | Sensors and Actuators A : Physical Fabrication of a true-Gaussian-shaped quartz crystal resonator[END_REF] [START_REF] Abe | One-chip multichannel quartz crystal microbalance (QCM) fabricated by Deep RIE[END_REF], gyroscope [START_REF] Grousset | Quartz-based vibrating MEMS fabricated using a wafer-bonding process with sealed cavities[END_REF] or timing applications [START_REF] Kubena | A Fully Integrated Quartz MEMS VHF TCXO[END_REF] [START_REF] Kamijo | Wafer-level quartz dry etching technology[END_REF]. Vertical etching of quartz over several tens microns is still challenging.

4.1

Process 1: Two DRIE steps process Process 1 is described in Figure 5. Chromium/Gold electrodes are evaporated and patterned, then a bilayer of resins (LOR 10B from MICROCHEM and S1813 from MICROPOSIT) is coated and a 3 µm thick nickel mask is sputtered. Lift-off is performed by soaking the samples into hot remover solution (Remover PG from MICROCHEM), once rinsed using isopropanol (IPA) and deionized water samples are ready to be etched. This process uses two DRIE steps: once the patterns of the resonators are made with the first etching step, the wafer is flipped and cavities are etched until the desired thickness and full liberation of the structures. Residual nickel is removed afterwards with nickel etchant solution (TFB from TRANSENE). This DRIE process uses a combination of a fluorocarbon gas tetrafluoromethane (CF4) and dioxygen: vertical etching of quartz takes place thanks to the presence of a passivating polymer [START_REF] Sankaran | Integrated feature scale modeling of plasma processing of porous and solid SiO2. I. Fluorocarbon etching[END_REF]) (HyunHo 1996) on the sidewalls. Because of the highly directional ion-milling, the passivating layer is barely absent at the bottom of the trenches.

This process suffers from a major issue: the wafer is cracking due to the thermal expansion of nickel submitted to high temperatures during the DRIE etching. This phenomenon has been observed in other work [START_REF] Boy | Collective fabrication of 20 MHz resonators by deep Reactive Ion Etching on 3" quartz wafers[END_REF]. The need for strong on-milling and the low thermal conductivity of quartz lead to this excessive heating up to 250°C on the quartz wafer.

4.2

Process 2: One DRIE step process with cavities etched with wet chemical solution To reduce the cracking, a process using cavities etched with a wet mixture of Hydrofluoric and Ammonium Fluoride (HF-NH4) has been developed and presented in Figure 6. It allows a reduction of the duration of the DRIE step by at least one half. Chromium/Gold layers are evaporated on both faces as they serve as electrodes and etching mask during the wet chemical etching step. Cavities are patterned and etched using the HF-NH4F mixture. Then, as in process 1, a bilayer of resist is coated, nickel is sputtered, lift-off is performed and resonators are etched using single-step DRIE. Only the back side of the vibrating part is left rough as the release of the structure is made by the DRIE step.

RESULTS

2D-LEM resonators with 3 targeted frequencies have been manufactured using both processes. Figure 7 shows type I resonators on wafer after nickel removal and a focus on a 3.0 MHz resonator. Electrical measurements were performed on the wafer using gold wires bonding under vacuum (10 -2 Pa) with impedance analyzer (Keysight EA990A). To estimate resonators performances, Q factor, resonance frequencies and serie resistances Rm are obtained using the standard Butterworth-Van Dyke (BVD) model and are discussed in this section.

5.1

Resonance frequency and measured dimensions Mean measured values and standard deviations of resonance frequencies versus the inverse of the total measured central beam length are shown in Figure 8. Data points are aligned with the curve obtained using FEM simulations, meaning that the f-L product is constant and the mode of interest is correctly selected. Lateral dimensions were measured with optical microscopy and the resolution is better than 1 µm for the used magnitudes.

Fig. 8 Mean experimental resonance frequencies versus the inverse of the mean experimental total length of the central beam for process 1 (blue triangles), for process 2 (red suares) and simulated frequency versus the inverse of the total length beam of the simulated structure using FEM (dashed black line).

Experimental resonance frequencies standard deviations were quite similar for both processes and larger for smaller devices 5% at 4.5 MHz, less than 2% at 3.0 MHz and less than 1% for the larger resonators vibrating at 2.2 MHz. The trend curve of the simulated structures exhibits a coefficient of 2.725 MHz.mm for the fundamental mode and is in good agreement with the estimation of equation ( 1). Figures 9 a),b),c) and d) show optical microscopy images focusing on type I and type II resonators with measured dimensions for 3 and 2.2 MHz resonators obtained with process 1. In average the measured dimensions of the central and lateral beams are few microns larger than the ones on the lithography mask, this lateral resolution loss is likely due to the resist layers stack and the thick nickel mask. Also the corners of the beams are not straight but slightly rounded. These process-induced flaws are more significant for the smaller resonators at higher frequencies as the resonators dimensions shrink. Some micro-masking is present on the backside cavities corresponding to the dark dots.

5.2

Experimental Q-factor Experimental Q factor were obtained using standard BVD model to fit the data. Figure 11 shows the mean experimental values of Q factors versus the experimental mean resonance frequencies for both types and both processes.

5.2.1

Results of process 1 Resonators vibrating at 2.2 and 3.0 MHz have been manufactured using process 1. As expected from (1), the f-L product is a constant at the first order for the different sizes of manufactured resonators. Experimental quality factors as high as 180 000 (Rm=5.4 kΩ) for type II resonators vibrating at 2.2 MHz and as high 150 000 (Rm=7.5 kΩ) for type II vibrating at 3.0 MHz were measured. This is however less than half of the expected Q meaning that there is at least another source of damping not included or underestimated in the previous estimation.

5.2.2

Results of process 2 The cracking with the process 2 was reduced which allowed the fabrication of smaller resonators -vibrating at 4.5 MHz -and type II resonators vibrating at 4.5 MHz showed best Q.f so far for fundamental mode reaching 5.3x10 11 Hz (Q=117 000, Rm=9.8 kΩ) and improving the state of the art of 2D-LEM quartz resonators. Yet, resonators showed overall lower quality factors as Q did not exceed 130 000 (Rm=8.5 kΩ) at 3 MHz and still type II resonators showed better performances.

Triple existe, mais the trip

5.3

Use of harmonics Using harmonics (n=1, 2, 3…in Eq (1)) to increase the Q.f product seems an adequate alternative. The first harmonic vibrating at three times the frequency of the fundamental is shown in Figure 12 for type I and II resonators. The central beam is divided into three smaller beams which is equivalent to a central beam of length three times smaller, hence the vibration at thrice the fundamental frequency. The decoupling mechanisms are still the same as in the case of the fundamental mode. Nevertheless, it is important to note that the design of the resonator and in particular the electrode network was optimized to excite the fundamental mode and not the harmonics. On one hand, the decoupling is less efficient, as can be seen in the Figure 8, lateral beams of the type I resonator and external beams of type II resonator move slightly, which implies an energy loss in the mounting areas. One the other hand, the electrodes array is not optimized for harmonics as a result a low charge collection and high associated motional resistances. Quality factors as high as 89 000 (Rm=42.7 kΩ) at 8.8 MHz (fundamental mode vibrating at 2.9 MHz) and 40 000 (Rm=101 kΩ) at 6.6 MHz (fundamental mode vibrating at 2.2 MHz) have been measured. The maximum Q.f product reached is 7.8x10 11 Hz which is an increase of 50% compared to the Q.f products obtained with fundamental mode. All the results are presented in the graph of the Figure 13. An optimized design and electrode network are necessary to achieve the full potential of harmonics: increasing the resonance frequency without excessively degrading the quality factor. In fact, as the frequency increases the strain energy stored in the vibrating part is increased whereas the energy contained in the gold layers varies slowly with frequency resulting in higher Qvisc. The dependency of the Qanchor for harmonic modes with the sidewall angle and resonator dimensions will be studied in future work.

5.4

Frequency temperature characteristics A 2D-LEM type I resonator manufactured with process 2 is mounted on a TO5 header, packaged in a copper case under vacuum. Packaged resonator is put in climatic chamber to perform frequency temperature characteristics, temperature was measured by a platinum sensor close to the packaged resonator. Figure .14 shows the normalized frequency temperature characteristics of a 2D-LEM type I resonator manufactured with process 2 for fundamental mode vibrating at 2.2MHz for temperatures from -40 to 30°C. The turn over temperature point of the FEM simulations was 10°C whereas the experimental turn over point was -9.5°C. However the trend of the curve is the same in experimental and simulated data. The downshifts of the turn over temperature point by almost 20°C for the fundamental mode and 23°C for the first harmonic are not fully understood yet. An offset of one degree along X axis in the cut angle can almost cause a 10°C variation of the turn over temperature point, as shown in Figure 16 where FEM simulations of normalized frequency temperature characteristics for angle values of ϴ = -2, 0 and +2 ° aound X-axis along with experimental data, turn over point is +30°C for ϴ = -2°, and -15°C for ϴ = +2° (in the FEM simuations the structure is rotated by an angle ϴ around the fixed axis coordinates which is equivalent to a rotation of in the crystalline axis of -ϴ) but it seems unlikely due to the high angle value (≈1.5°) required to fit the data. Rotations around the two other cristalline axes did not lead to significant shift of the turn over point. Another hypothesis implies an effect of the surface roughness that can shift the turn over point as it had been observed in our team for flexural quartz resonators: turn over point was shifted to lower temperatures with increasing surface roughness. Further investigations are necessary to confirm the precise value of the resonators crystalline orientation and to measure frequency temperature characteristics of resonators manufactured with process 1 without the pyramids on the backside cavities

DISCUSSION

In this section, the discrepancy between the measured and theoretical quality factor are discussed. As mentioned earlier, resonators made with process 2 showed qualitatively lower Q factor than those obtained with process 1. Moreover, experimental quality factors obtained are lower than those predicted.

6.1

Process 1 and process 2 Resonators manufactured with process 1 showed better Q factors than the ones made with process 2. Figure 17a) shows the backside of a resonator manufactured with process 1, apart from some micro-masking induced by contaminents, the surface is smooth with Ra=4 nm measured with a mechanical profilometer (KLA Tencor D-100). The reduced quality factors of the process 2 are likely due to the presence of pyramids on the back side as illustrated in Figure 17b) induced by the wet chemical etching using HF-NH4F. As there are no pyramids on the top face, the structure obtained with process 2 is less symmetrical than in the case of process 1 and is also way more rougher with Ra=1.9 µm, which may amplify the loss energy through the anchors: this lowers Qanchor hence the overall decrease of the experimental Q factor.

6.2

Experimental Q factors The simulated structure for the calculation of Qanchor is not exactly the one obtained with process 1. Indeed, in the simulations the sidewalls were perfectly vertical whereas DRIE cannot ensure a perfect verticality as illustrated in Figure 8. FEM simulations give the dependency of the Qanchor versus sidewall angle in Figure .8. It appears that energy trapping in the vibrating part becomes rapidly less efficient with an, even slightly, non-vertical sidewall angle. To optimize Qanchor, resonator design must be made so that the mounting areas match vibration nodes where the strain energy is almost null. As soon as the sidewall angle is non-zero, the structure is no longer symmetrical and vibration nodes are slightly displaced whereas the mounting areas are in the same place that for the ideal structures. For an angle around 3° which is likely to be reached with the current process 1, as seen in the detail of Figure 18, Qanchor falls down around 400 000 for the 3 MHz resonator as shown in Figure .18 and 500 000 for the 2.2 MHz and explains the experimental values. One can note that the larger resonators with lower resonance frequencies present better performances than the smaller onesvibrating at higher resonance frequencies -and that was also observed in the work of Kawashima (Kawashima1990). Type II resonators showed better performances than equivalent previous standard type I design, which seems to indicate either a more efficient design with the out-of-phase extension-compression mode of the lateral beams or a less sensitive dependence of Qanchor with sidewall angle.

Conclusions

In this paper, new 2D-LEM quartz resonators with potential Time&Frequency applications have been demonstrated with first prototypes showing encouraging electro-mechanical characteristics providing a product near 10 12 Hz and a Q factor at the state of art for DRIE made quartz resonators. The two main limiting damping mechanisms have been identified: the energy losses arise from the viscous dissipation in the gold layer forming the electrodes and the mounting areas. First characterizations of harmonic modes were performed and showed promising results. Future work will be focused on increasing the resonance frequency and maintaining or improving the quality factor thanks to better control of the sidewall angle.

Fig. 1

 1 Fig.1 a) Concept of the 3D-LEM: a ring around the pillar with an out-of-phase length-extension vibration, b) 1.5-inches square quartz wafer including nine 3D-LEM made using wet etching, c) detail of the micro-pillar.

Fig. 2

 2 Fig.2 Electrode scheme for the actuation of the extension mode of a Quartz Z-cut beam oriented along the Y-axis.

Fig. 3 a

 3 Fig.3 a) Amplified nodal displacements of the fundamental mode of the central beam of a type I resonator during its length-extension mode vibration b) X component of the type I resonator showing the flexure motion of the two lateral beams, c) Y component of the type I resonator showing the extensioncompression displacement of the central beam.Table1Mean dimensions of type I resonator and targeted frequencies

Fig. 4

 4 Fig.4 a) Amplified nodal displacements of the fundamental mode of the central beam of a type II resonator during its length-extension mode vibration b) X component of the type II resonator showing the flexure motion of the two lateral beams, d) Y component of the type II resonator showing the extensioncompression displacement of the central beam and the predominant out-of-phase length extension of lateral beams.

Fig. 5

 5 Fig.5 Process flow 1 using two DRIE steps

Fig. 6

 6 Fig.6 Process flow 2 using wet etching to create cavities and one DRIE step on top face. DRIE parameters are the same as in process 1

Fig. 7

 7 Fig.7 Left, top line type I resonators vibrating at 2.2 MHz, bottom line at 3 MHz. Right, focus on a 3 MHz resonator. Resonators were manufactured using process1.

Fig. 9

 9 Fig.9 a) Optical image of a type I resonator vibrating at 3 MHz magnification factor (MF) = 100, b) focus on the central beam MF = 200, c) Optical image of a type II resonator vibrating at 2.2 MHz MF =100, d) focus on the central beam MF = 200. Resonators were manufactured using process 1.

Figures

  Figures 10 a), b), c), d) ,e) and f) show optical microscopy images focusing on type I and type II resonators with measured dimensions for 4.5, 3 and 2.2 MHz resonators obtained with process 2. Black dots on the top face correspond to some residual attack from the wet chemical echting of the cavities and the darker areas correspond to the pyramids on the backside cavities.

Fig. 10

 10 Fig.10 a) Optical image of a type II resonator vibrating at 4.5 MHz magnification factor (MF) = 50, b) focus on the central beam MF = 200, c) Optical image of a type I resonator vibrating at 3 MHz MF =50, d) focus on the central beam MF = 200, e) Optical image of a type II resonator vibrating at 2.2 MHz magnification factor (MF) = 50, f) focus on the central beam MF = 200. Resonators were manufactured using process 2.

Fig. 11

 11 Fig.11 Mean experimental values of Q factors versus mean experimental resonance frequencies for both types and both processes.

Fig. 12

 12 Fig.12 Amplified nodal displacements of the vibrating part obtained using FEM simulations of the first harmonic mode (n=1). Left: type I resonator, right type II resonator

Fig. 13

 13 Fig.13 Figure of merit Q.f of 2D-LEM resonators investigated in this work for both fundamental and harmonic modes and 2D-LEM quartz resonators from literature (Dinger1981) (Kawashima1990).

Fig. 14

 14 Fig.14 Frequency temperature characteristics of the fundamental mode of type I resonator vibrating at 2.2 MHzA shift is again observed for the first harmonic mode vibrating at 6.6 MHz as shown in Figure15: the FEM simulations turn over point was 123°C whereas experimental turn over point was 99°C.

Fig. 15

 15 Fig.15Frequency temperature characteristics of the first harmonic mode of type I resonator vibrating at 6.6 MHz (fundamental mode vibrating at 2.2 MHz).

Fig. 16

 16 Fig.16 Frequency temperature characteristics of the first harmonic mode of type I resonator vibrating at 2.2 MHz and FEM simulations with differents angles around de X-axis. Indicated values of ϴ correspond to the values entered in the FEM simulation software.

Fig. 17

 17 Fig. 17 a) Optical image of the back side of a resonator using the double steps DRIE process (process 1), b) optical image of the back side of a resonator using wet etching (process 2) and presence of pyramids with several tens of microns diameters MF = 200.

Fig. 18

 18 Fig.18 Simulations using FEM of the Qanchor vs sidewall angle for a 3 MHz type I resonator. In detail, scanning electron microscopy (SEM) of a 90 µm deep cross section beam

Table 2 Mean dimensions of type II resonator and targeted frequencies

 2

Acknowledgment

This work is partially supported by the Centre National d'Etudes Spatiales CNES.