D Faux

O Thomas

S Grondel

É Cattan

Ivana Kovacic

Dynamic simulation and optimization of artificial insect-sized flapping wings for a bioinspired kinematics using a two resonant vibration modes combination

Keywords: Resonant motion Micro, Nano electromechanical system (M, NEMS) Modal combination Design of structural vibration Optimization

This paper addresses the design of the elastic structure of artificial wings to optimize their dynamical behaviour to reproduce insect wings kinematics. Our bioinspired kinematics is based on the original concept of using the resonant properties of the wing structure in order to combine the motion of two vibration modes, a flapping and a twisting mode, in a quadrature phase shift. One way of achieving this particular combination is to optimize the geometry and elastic characteristics of the flexible structure such that the two modes are successive in the eigenspectrum and close in frequency. This paper first proposes a semi-analytical model, based on assembled Euler-Bernoulli beams, to understand, compute and optimize the artificial wing dynamic vibrations. Then, using this model, it is shown that it is possible to obtain several artificial wing structures with a flapping and a twisting mode close in frequency. Finally, experimental validations are performed on micromachined insect-sized prototypes to validate the model and the concept.

Introduction

Over the course of their evolution, flying insects have developed and improved their fly mechanisms and demonstrate impressive capacities in terms of agility and maneuverability. Indeed, among flying species, insects show remarkable performances in regard to hovering, backward flight or sudden acceleration and are unrivaled in terms of weight and scale. The key to such abilities mainly relies on their wings motions, and more particularly on their specific trajectories. Hence, it is now well established, thanks to the substantially improvement of the insects flight mechanisms understanding, that the generated lift results of a wide range of unconventional aerodynamic mechanisms [START_REF] Srygley | Unconventional lift generating mechanisms in free flying butterflies[END_REF] as well as specific kinematics induced by the flexible characteristics of the wings [START_REF] Wootton | The functional morphology of the wings of odonata[END_REF]. Moreover, the enhancement of the high-speed camera technologies provided the possibility to describe more precisely this particular wing kinematics [START_REF] Ellington | The aerodynamics of hovering insect flight. III. Kinematics[END_REF][START_REF] Ennos | Inertial and aerodynamic torques on the wings of Diptera in flight[END_REF].

An other central element in the understanding of the insect flight relies on their wings structure made of a thin membrane supported by flexible veins. In fact, this structure flexibility is a key element of the insects flight capacities since this characteristic notably allows the wing deformation in order to adapt its shape dynamically during flight [START_REF] Chen | On the natural frequencies and mode shapes of dragonfly wings[END_REF]. However, the exact role of this particularity on the aerodynamic performances remains controversial. In literature, some studies [START_REF] Wootton | The functional morphology of the wings of odonata[END_REF][START_REF] Ennos | The importance of torsion in the design of insect wings[END_REF][START_REF] Ennos | The inertial cause of wing rotation in Diptera[END_REF] demonstrate that this flexibility enables a camber of the wing during the flapping which generates higher lift forces than with rigid wings [START_REF] Mountcastle | Aerodynamic and Functional Consequences of Wing Compliance[END_REF][START_REF] Young | Details of insect wing design and deformation enhance aerodynamic function and flight efficiency[END_REF].

On the contrary, other simulation works [START_REF] Zhao | Aerodynamic effects of flexibility in flapping wings[END_REF] expose a reduction of aerodynamic performances with higher flexibility for small angles of attack. Furthermore, as discussed by Refs. [START_REF] Bao | Improved micromachining of all SU-8 3D structures for a biologically-inspired flying robot[END_REF][START_REF] Bao | Design and fabrication of insect-inspired composite wings for MAV application using MEMS technology[END_REF], this flexibility property is also contributing for their solidity or their resistance to impacts.

Much discussions are also oriented about the resonant mechanisms developed by some insects. Several studies [START_REF] Michelin | Resonance and propulsion performance of a heaving flexible wing[END_REF][START_REF] Masoud | Resonance of flexible flapping wings at low Reynolds number[END_REF][START_REF] Zeng | Two dimensional, non-contact measurement of the natural frequencies of dragonfly wings using a quadrant position sensor[END_REF][START_REF] Sims | A structural dynamic analysis of a Manduca sexta forewing[END_REF] highlight the use of the wing structure natural frequencies to amplify a small amplitude motion. While on the other hand, some works [START_REF] Chen | On the natural frequencies and mode shapes of dragonfly wings[END_REF][START_REF] Ha | Relationship between wingbeat frequency and resonant frequency of the wing in insects[END_REF] tend to indicate that a wing beat that is different from the wings structure resonant frequencies is used. However, it is also demonstrated that the resonant mechanism relies on other body parts such as the thorax [START_REF] Hrncir | Thoracic vibrations in stingless bees (Melipona seminigra): resonances of the thorax influence vibrations associated with flight but not those associated with sound production[END_REF].

When creating an artificial flapping wing prototype, there is a broad design space ahead, outlined by the number of flying insects, and designers have to crawl through it before finding out the correct combination of parameters meeting their requirement. This approach is even more critical when original technology choices are introduced and no known-how is available paving the way, as it is the case for the resonant wings proposed in our original concept [START_REF] Faux | Two modes resonant combined motion for insect wings kinematics reproduction and lift generation[END_REF]. To alleviate this burden, analytic and numerical tools can be envisioned to assist and hasten these preliminary design tasks. Obviously, most of the numerical tools published until today are dedicated to model accurately insect flight [START_REF] Smith | Simulating moth wing aerodynamics -towards the development of flapping-wing technology[END_REF][START_REF] Ho | Unsteady aerodynamics and flow control for flapping wing flyers[END_REF][START_REF] Singh | Dynamics and Aeroelasticity of Hover-Capable Flapping Wings: Experiments and Analysis[END_REF][START_REF] Hamamoto | Application of fluid-structure interaction analysis to flapping flight of insects with deformable wings[END_REF][START_REF] Liu | Micro air vehicle-motivated computational biomechanics in bio-flights: aerodynamics, flight dynamics and maneuvering stability[END_REF][START_REF] Nakata | A fluid structure interaction model of insect flight with flexible wings[END_REF][START_REF] Dawson | Flapping micro air vehicle: wing fabrication and analysis[END_REF]. Those attempts provide a good understanding of the physics involved in flapping wings but seem too specific and accurate for design purpose. Of course, more practical numerical tools [START_REF] Rakotomamonjy | Modelization and kinematics optimization for a flapping-wing micro air vehicle[END_REF][START_REF] Orlowski | Modeling and simulation of nonlinear dynamics of flapping wing micro air vehicles[END_REF][START_REF] Macfarlane | Quasi-steady and computational aerodynamics applied to hovering drosophila dynamics[END_REF] are existing but they are devoted to computational intensive studies such as flight mechanisms. Few of them have been used for sizing purposes [START_REF] Ho | Unsteady aerodynamics and flow control for flapping wing flyers[END_REF][START_REF] Stanford | Aeroelastic optimization of flapping wing venation: a cellular division approach[END_REF][START_REF] Agrawal | Design of bio-inspired flexible wings for flapping-wing micro-sized air vehicle applications[END_REF][START_REF] Isogai | Optimum aeroelastic design of resonance type flapping wing for micro air vehicles[END_REF][START_REF] Khan | Optimal hovering kinematics of flapping wings for micro air vehicles[END_REF] and even fewer for preliminary design tasks [START_REF] Combes | Fluid-structure interaction simulation for the design of bio-inspired micro air vehicle wings[END_REF][START_REF] Stewart | Parametric representation and shape optimization of flapping micro air vehicle wings[END_REF].

Within this context, this article focuses on the design and the optimization of the elastic structure of the artificial wings of a micromachined insect-sized prototype. Following the concept introduced in Ref. [START_REF] Faux | Two modes resonant combined motion for insect wings kinematics reproduction and lift generation[END_REF], we propose to generate a motion inspired by the particular insects wings kinematics with fully flexible artificial wings. More precisely, the leading idea is to combine the motion of two vibration modes, a flapping mode and a twisting mode, close to their resonance, to artificially create a wing kinematics similar to the one of insects. This situation is obtained by designing the geometry and the elastic structure of the wings such that the two modes are successive in the eigenspectrum and close in frequency. Then, only one actuator is necessary and the particular motion of the wings is solely a result of the elastic properties of the wings, actuated at a single well chosen frequency.

The purpose of this paper is first to propose and describe a semi-analytic model to optimize the flexible wing structure of the wings. Then, we aim at demonstrating that it is possible to obtain two vibration modes -a flapping and a twisting mode -successive in the eigenspectrum and close in frequency. This is achieved by selecting several classes of geometry and optimizing their parameters. The paper reports the result of this theoretical study. Finally, for selected optimized geometries, some micromachined insect-sized prototypes have been built. For each ones, the paper reports their vibration properties, that have been experimentally measured, in order to validate the proposed model and the concept.

Vibration mode combination concept

Generally, flapping micro and nano air vehicles feature artificial wings composed of an articulated rigid leading edge attached to the thorax and a flexible membrane [START_REF] Wood | The first takeoff of a biologically inspired at-scale robotic insect[END_REF][START_REF] Raney | Mechanization and control concepts for biologically inspired micro aerial vehicles[END_REF]. To produce an appropriate wing slope and lift off [START_REF] Ramananarivo | Rather than resonance, flapping wing flyers may play on aerodynamics to improve performance[END_REF] the most common mechanism is to impose a large flapping motion on the rigid leading edge, by using a mechanical amplifier [START_REF] Wood | The first takeoff of a biologically inspired at-scale robotic insect[END_REF][START_REF] Wood | Design, fabrication, and analysis, of a 3DOF, 3cm flapping-wing MAV[END_REF][START_REF] De Croon | Design, aerodynamics and autonomy of the DelFly[END_REF], and to exploit the passive twisting of the wings. In this case, the wing is equivalent to a single degree of freedom system in twisting, with the flapping displacement of the leading edge prescribed by the actuator.

Conversely, the approach followed here is based on an actuation principle without any mechanism to transmit and transform the actuator motion. Following the concept introduced in Ref. [START_REF] Faux | Two modes resonant combined motion for insect wings kinematics reproduction and lift generation[END_REF], we propose to consider fully flexible artificial wings and to actuate them in a particular resonant motion, similar to the kinematics of real insect wings. In particular, we consider our wings as a deformable vibrating structure featuring two particular vibration modes, with shapes of the form of a flapping motion for the first mode and a twisting motion for the second one. This structure must be designed in such a way that the eigen frequencies of these two modes appear successive in the eigenspectrum and are close enough to enable their combination. Then, by selecting a particular driving frequency between the two resonances of the modes, it is possible to actuate the two modes simultaneously with a quadrature phase shift between the two motions. Basically, the main idea is to obtain a flapping mode shape which is providing the leading edge displacement and a twisting mode shape which is impacting on the trailing edge one. As shown in Fig. 1 with the motions amplitudes (b) and the associated wing tip tracking (c), by reaching the quadrature (phase shift equal to 𝜋∕2) between the two motions, it is theoretically possible to obtain an asynchronous kinematics close to those observed on insects. Indeed, the wing slope is maximal in the middle of a half flapping cycle and the slope reversal occurs at its end.

This concept, introduced and detailed in Ref. [START_REF] Faux | Two modes resonant combined motion for insect wings kinematics reproduction and lift generation[END_REF], has then two remarkable properties: the capacity to reproduce an insect wing kinematics thanks to a combination of flapping and twisting structural deformations and the possibility to obtain large amplitude since the actuation on a quadrature frequency is close to the resonances of the modes. However, all this is possible (i) if two modes in the spectrum have the shape of a flapping mode and a twisting mode; (ii) if those modes appear successively in the spectrum and (iii) if their resonances are close enough so that the quadrature frequencies are also close enough to the two resonances peaks for the two modes to be actuated with a suitable amplitude. The remaining of the paper addresses those three points.

Wing structure modelling

The developed prototypes (a design example is shown in Figs. 1(a) and 2) are composed of several parts. The wings are based on a polymeric flexible skeleton, which has the shape of assembled beams, on which a thin membrane is deposited. The skeleton is linked to a thick thorax, also made in the same polymeric material. A single electromagnetic actuator, composed of a magnet glued on the central part of the skeleton and a coil supported by the thorax, creates the necessary force to drive the wings in vibration.

As explained in the previous section, our concept is based on the actuation of two vibration modes of the flexible wings. We choose to consider in our modelling only the flexible structure of the wings, clamped to the thorax, that is assumed infinitely rigid since its thickness is ten times the one of the beams of the skeleton. Moreover, we also choose to neglect the wing membranes in our modelling. Their structural influence, in term of added mass and stiffness, can be neglected since they are very thin (400 nm for the membrane, 80 μm for the beams of the skeleton, giving a 200 ratio between thicknesses) and not stressed at rest. However, their aerodynamic influence has to be carefully addressed. When the wings are vibrating, their membranes stir the fluid, the effect of which is mainly to add mass and aerodynamic damping to the vibration. Those effects are experimentally addressed in section 3.1. It is shown that the added mass effect of the fluid lowers the resonance frequencies of the structure whereas the aerodynamic damping lowers their amplitude and sharpness. However, the main resonant response of the system does not change fundamentally, the two resonance peaks staying close to each other. Those aerodynamic effects are theoretically difficult to take into account in a simple model. Since our goal is the design of the flexible structure of the wings to bring close two particular vibration modes, we simply choose to neglect in our model the influences of the membranes, both structural and aerodynamic, on the dynamics of the wing structures. It will be shown in section 4.3 that this modelling can be successfully used as a pre-design tool.

Our design procedure is thus mainly based on a modal analysis of the wing structure, to identify its vibration modes. The wing structure is modeled by an assembly of Euler-Bernoulli beams, with degrees of freedom restricted to bending, in the flapping -out of plane -direction (z in Fig. 1(a)), and twisting. We opt for a semi-analytic modelling, for which the exact solution for the deformation of each straight beam is combined to suitable boundary and assembly conditions. In comparison with a finite elements modelling with a commercial code, this kind of method allows to easily restrict the number of degrees of freedom to exploit, and so to limit the number of calculated non-contributory modes. In this way, only the modes inducing a deformation in the out of plane -flapping -direction (direction z in Fig. 1(a) perpendicular to the wing plane) are considered, the others (with deformed shapes in the in plane directions) being useless, especially since the chosen actuation solution cannot allow their excitation (the actuator imposes a force in z direction). This restriction is mandatory when an optimization loop is under concern, as described in section 4.4. With the results of the modal analysis, it is then possible to compute the frequency response function (FRFs) of the structure to identify the driving frequencies that leads to a quadrature phase shift between the two modes under interest. To obtain FRFs similar to the one obtained experimentally at air pressure, we simply add a linear damping in the model, with a quality factor of the order of 10, as observed in Fig. 3(b) (refer to section 3.1 for more details).

It must also be noted that the structure of the wings of our prototypes are similar to those of real insects: a skeleton with a membrane. Moreover, as shown in Refs. [START_REF] Bao | Improved micromachining of all SU-8 3D structures for a biologically-inspired flying robot[END_REF][START_REF] Bontemps | Design and performance of an insect-inspired nano air vehicle[END_REF][START_REF] Dargent | Micromachining of an SU-8 flapping-wing flying micro-electro-mechanical system[END_REF], the SU-8 polymer material selected for our wings has mechanical properties very close to those of natural insects. Despite of that, since we rely on the mode combination mechanism to create the suitable wing kinematics and produce lift, and also because our design procedure is based on parametric optimization, our wing structures are much simpler that those of real insects, with less beam members and a geometry that could be considered different from that of insects, although there is a very important variety of insect wing shape. According to this point, we thought that our wings will get closer and closer to some insect wings.

Aerodynamic effects on the vibration

To experimentally estimate the aerodynamic effects on the vibratory response of the wings, we measured the frequency response of one point of the wing structure for several air pressures, by using a vacuum chamber. Precisely, a prototype is placed in a vacuum chamber with a glass window, through which the laser beam of a laser vibrometer (Polytec OFV-5000) enables to measure the velocity signal of a selected point of the wing structure of the prototype. The latter is driven in vibration by a white noise signal sent to the electromagnetic actuator. Several measures are performed for selected values of the air pressure, from 6 μbar to the atmospheric pressure (1 bar). Fig. 3(a) shows some spectrum of the velocity signal (obtained by a simple fast Fourier transform), on a frequency band around the two modes of interest, for seven values of the air pressure. One can observe that the two resonances are clearly visible. Moreover, there is an attenuation of the two resonances and a frequency shift to the low frequencies, as the air pressure increases. With a simple fit of the low frequency resonance peak with a Lorentzian function, its resonance frequency f 1 and quality factor Q are estimated. Fig. 3(b andc) show the evolution of f 1 and Q as a function of the air pressure. This quantitatively assess the effect of the aerodynamic on the vibratory response of the prototype: an increase of the damping and a frequency shift of the resonances toward the low frequencies.

Quantitatively, the obtained results on the increase of the quality factor and the decrease of the resonance frequencies are related to a given prototype at a given driving amplitude (the aerodynamic force is at first order a quadratic function of the amplitude, which leads to a linear increase of the quality factor as a function of the amplitude). More qualitatively, (i) the two peaks are still visible at atmospheric pressure and (ii) they show almost the same frequency shift, which leads to a frequency difference almost independent of the atmospheric pressure. We choose to consider those two qualitative results valid for any of the proposed prototype architecture.

Continuum problem

We first consider a simple "L" shaped wing structure to explain our modelling strategy. It is modeled with an assembly of two linear Euler-Bernoulli beams with deformations in bending and twisting, as shown in Fig. 4 with beams 1 and 2 representing the leading edge and the chord respectively.

Based on classical mechanical vibrations theory [START_REF] Géradin | Mechanical Vibrations: Theory and Applications to Structural Dynamics[END_REF], we consider a straight beam i with a transverse displacement along z axis w i (x i , t) and a section twisting rotation along the beam axis 𝜙 i (x i , t) (with x i the coordinate along the beam axis t the time).

They satisfy the following differential equations:

∀x i , ∀t { EI i w ′′′′ i + 𝜌S i ẅi = P i GJ i 𝜑 ′′ i -𝜌J i φi = 0 , i = 1, 2 (1)
with • ′ = ∂ • ∕∂x and ̇• = 𝜕 • ∕𝜕t. P i (x i , t) represents the transverse force (no external forcing in twisting is considered), E the material Young's modulus, S the beam cross section area, I and J the cross section moment in bending and twisting respectively, 𝜌 the material density, and G = E∕(2(1 + 𝜈)) the Coulomb modulus with 𝜈 the Poisson ratio.

With the beams oriented with x 1 = OA∕ ∥ OA ∥ and x 2 = -AB∕ ∥ AB ∥, the bending moment M i , the twisting moment N i and the shear force T i are expressed as:

⎧ ⎪ ⎨ ⎪ ⎩ M i = E i I i w ′′ i N i = G i J i 𝜑 ′ i T i = -M ′ i = -E i I i w ′′′ i (2)
Recalling that the section bending rotation for a linear Euler-Bernoulli beam is defined by w ′ , we can write this assembled beams system with the following boundary and continuity conditions:

Beam 1 clamped end (x 1 = 0 in O) ⎧ ⎪ ⎨ ⎪ ⎩ w 1 (0) = 0 w ′ 1 (0) = 0 𝜑 1 (0) = 0 (3) Continuity beam 1/beam 2 (x 1 = L, x 2 = 0 in A) ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ w 1 (L) = w 2 (0) w ′ 1 (L) = -𝜑 2 (0) 𝜑 1 (L) = w ′ 2 (0) M 1 (L) = N 2 (0) N 1 (L) = -M 2 (0) T 1 (L) = -T 2 (0) (4) Beam 2 free end (x 2 = -H in B) ⎧ ⎪ ⎨ ⎪ ⎩ M 2 (-H) = 0 T 2 (-H) = 0 N 2 (-H) = 0 (5)
Eq. ( 3) describes the clamped condition of beam 1 with a null transverse displacement, bending section rotation and twisting rotation. Continuity Eq. ( 4) represent the coupling between the bending of one beam and the twisting of the other one (and inversely) by binding on the connection an equal transverse displacement for the two beams, a twisting rotation of a beam equal to the bending section rotation of the other one (in absolute value), a twisting moment of one equal to the bending moment of the other (in absolute value), and the shear forces equilibrium. Finally, Eq. ( 5) describe the free condition of beam 2 with a null bending moment, twisting moment and shear force.

Modal analysis

A modal analysis of the structure is performed with the aim of determining its eigenmodes (eigenfrequencies and mode shapes). In order to generalize and simplify the development, equations are rendered dimensionless with regard to a characteristic length, chosen as the length L of the leading edge (beam 1), which also defines the wingspan. Dimensionless variables (with an overline) are then written:

x = x L , t = t L 2 √ EI 1 𝜌S 1 , P i = P i L , Q i = Q i L H = H L , w i = w i L (6)
To simplify the reading of the equations, the dimensionless variables will now be written without the overline (except the length ones as H to avoid any misleading).

For the modal analysis, Eq. ( 1) are rewritten dimensionless and in free vibrations (P i = 0 and Q i = 0), as:

∀x i , t { w ′′′′ i + 𝛿 4 i ẅi = 0 𝜑 ′′ i -𝛾 2 φi = 0 , i = 1, 2 (7)
with

𝛿 4 1 = 1, 𝛿 4 2 = S 2 I 1 S 1 I 2 , 𝛾 2 = EI 1 GS 1 L 2 . ( 8
)
Variables 𝛿 i and 𝛾 represent the stiffness ratio in bending and twisting respectively, and the choice of the variable t give 𝛿 4 1 = 1.

It can be noted that 𝛾 is not dependent of the considered beam.

By separating space and time variables, the deformations in bending and twisting are expressed as:

{ w i (x, t) = Φ i (x)q(t) 𝜑 i (x, t) = Ψ i (x)q(t) , i = 1, 2 (9)
with q(t) the modal coordinate. Eq. ( 1) in bending and twisting are then turned into:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ q = -𝜔 2 q Φ ′′′′ i -(𝛿 i 𝛽) 4 Φ i = 0 Ψ ′′ i -(𝛾𝛽 2 ) 2 Ψ i = 0 , i = 1, 2 (10)
with 𝛽 4 = 𝜔 2 and 𝜔 the angular frequency of the associated eigenmode. The solutions of the system are then expressed as:

{ Φ i (x) = A i cos(𝛿 i 𝛽x) + B i sin(𝛿 i 𝛽x) + C i cosh(𝛿 i 𝛽x) + D i sinh(𝛿 i 𝛽x) Ψ i (x) = E i cos(𝛾𝛽 2 x) + F i sin(𝛾𝛽 2 x) , i = 1, 2 (11)
where A i , B i , C i , D i , E i , F i are constants. Based on Eq. ( 2), the boundary and continuity conditions of the structure can then be turned into:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Φ 1 (0) = 0 Φ ′ 1 (0) = 0 Ψ 1 (0) = 0 Φ 1 (1) = Φ 2 (0) Φ ′ 1 (1) = -Ψ 2 (0) Ψ 1 (1) = Φ ′ 2 (0) EI 1 GJ 2 Φ ′′ 1 (1) = Ψ ′ 2 (0) GJ 1 EI 2 Ψ ′ 1 (1) = -Φ ′′ 2 (0) EI 1 EI 2 Φ ′′′ 1 (1) = -Φ ′′′ 2 (0) Φ ′′ 2 (-H) = 0 Φ ′′′ 2 (-H) = 0 Ψ ′ 2 (-H) = 0 (12)
To determine the eigenmodes, Eq. ( 11) and the above conditions are written under the following matrix form:

M(𝛽) V = 0 (13) with V = [A 1 B 1 C 1 D 1 A 2 B 2 C 2 D 2 E 1 F 1 E 2 F 2 ]
T , and the matrix M(𝛽) defined by:

M = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝚽 1 (0) 0 0 0 0 0 0 0 0 𝚽 ′ 1 (0) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 𝚿 1 (0) 0 0 𝚽 1 (1) -𝚽 2 (0) 0 0 0 0 𝚽 ′ 1 (1) 0 0 0 0 0 0 𝚿 2 (0) 0 0 0 0 -𝚽 ′ 2 (0) 𝚿 1 (1) 0 0 EI 1 𝚽 ′′ 1 (1)
0 0 0 0 0 0 -GJ 2 𝚿 ′ 2 (0) 0 0 0 0 EI 2 𝚽 ′′ 2 (0) GJ 1 𝚿 ′ 1 (1) 0 0 EI 1 𝚽 ′′′ 1 (1) EI 2 𝚽 ′′′ 2 (0) 0 0 0 0 0 0 0 0 𝚽 ′′ 2 (-H) 0 0 0 0 0 0 0 0 𝚽 ′′′ 2 (-H) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 𝚿 ′ 2 (-H) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (14)
with each line corresponding to a condition and each column to one of the coefficients

[A i B i C i D i E i F i ].
To enhance the matrix reading, the vectors 𝚽 and 𝚿 correspond to

𝚽 i (x) = [cos(𝛿 i 𝛽x) sin(𝛿 i 𝛽x) cosh(𝛿 i 𝛽x) sinh(𝛿 i 𝛽x)] and 𝚿 i (x) = [cos(𝛾𝛽 2 x) sin(𝛾𝛽 2 x)]).
To determine the structure's eigenmodes, we are looking for the non-trivial solutions of system [START_REF] Michelin | Resonance and propulsion performance of a heaving flexible wing[END_REF], so the ones which provide det(M) = 0. 𝛽 values which nullify the M matrix determinant then correspond to the eigenfrequencies of the structure, with the following formula which results from 𝛽 4 = 𝜔 2 and Eq. ( 6):

f = 𝛽 2 2𝜋L 2 √ EI 1 𝜌S 1 . ( 15
)
The mode shapes associated to those eigen modes, corresponding to the solutions [START_REF] Bao | Improved micromachining of all SU-8 3D structures for a biologically-inspired flying robot[END_REF], can then be determined by calculating the unknown coefficients

[A i B i C i D i E i F i ]
for a given 𝛽 value. To this end, since the determinant det(M) is null, system [START_REF] Michelin | Resonance and propulsion performance of a heaving flexible wing[END_REF] is homogeneous and singular and have an infinity of solutions V. In order to determine the non-null values of the solutions

V = [A i B i C i D i E i F i ],
it is possible to normalize the eigenmodes by setting arbitrarily the value of one of them (for example F 2 = 1), then define a non-homogeneous new system (with a unique solution) to solve. Once the coefficients are determined, we can obtain the mode shapes by introducing them in Eq. [START_REF] Bao | Improved micromachining of all SU-8 3D structures for a biologically-inspired flying robot[END_REF]. It must be noticed that their amplitudes will be relative to the value of the used normalization coefficient. The above described methodology is a special case of spectral methods [START_REF] Lee | Spectral Element Method in Structural Dynamics[END_REF]. Eq. ( 13), solved by the proposed method, can result in a badly conditioned problem which prevents the solving and the obtention of the coefficients. In order to overcome this kind of problems, some other algorithms [START_REF] Wittrick | A general algorithm for computing natural frequencies of elastic structures[END_REF][START_REF] Maurini | Numerical methods for modal analysis of stepped piezoelectrics beams[END_REF] can be used with an indirect determination of the eigen frequencies. In the particular case of the structures described in this article, no such problems have been encountered. The 𝛽 values are simply obtained in two steps. First, several initial values of 𝛽 that nullify det(M) are obtained by computing det(M) in a given 𝛽 range with a coarse sampling. Then, for each initial beta value, a Newton-Raphson algorithm is used to obtained 𝛽 with a given suitable accuracy.

Frequency response functions

In order to compute the deflection shapes of the structure at a given driving frequency, we will now focus on the Frequency Response Functions (FRF) in bending. Based on the classical vibratory mechanics [START_REF] Meirovitch | Fundamentals of Vibrations[END_REF], bending and twisting deformations variables w(x, t) and 𝜙(x, t) verify Eq. ( 1) and can be expressed as a modal superposition as following:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ w i (x i , t) = N ∑ n=1 Φ (n) i (x i )q n (t) 𝜑 i (x i , t) = N ∑ n=1 Ψ (n) i (x i )q n (t) (16)
with i the beam number, n the mode number, N the number of modes taken into account, Φ (n) i (x i ) and Ψ (n) i (x i ) the mode shapes in bending and twisting of the mode n, and q n (t) the modal coordinate. By substituting Eq. ( 16) into Eq. ( 1) and using the following orthogonality properties:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 2 ∑ i=1 ∫ L i 0 [𝜌S i Φ (p) i Φ (q) i + 𝜌J i Ψ (p) i Ψ (q) i ]dx = m p 𝛿 pq , 2 ∑ i=1 ∫ L i 0 [EI i Φ (p) i Φ ′′′′(q) i + GJ i Ψ (p) i Ψ ′′(q) i ]dx = 𝜔 2 p m p 𝛿 pq , ( 17
)
with L i the ith beam length, 𝛿 pq the Kronecker symbol (𝛿 pq = 1 for p = q and 0 for p ≠ q) and m n the modal mass of mode n, the nth modal coordinate verify:

qn + 2𝜉 n 𝜔 n ̇qn + 𝜔 2 n q n = E n m n , E n = 2 ∑ i=1 ∫ L i 0 P i Φ (n)
i dx [START_REF] Hrncir | Thoracic vibrations in stingless bees (Melipona seminigra): resonances of the thorax influence vibrations associated with flight but not those associated with sound production[END_REF] with 𝜉 n the damping ratio of mode n and E n the modal forcing.

By modelling the structure forcing as a local periodic force applied at a point x 0 of the kth beam with an amplitude P 0 such as P j (x, t) = P 0 𝛿(xx 0 )cos(Ωt), we can then obtain from Eq. ( 16) the expression of the FRF at a point x of beam i such as:

H(x, Ω) = ŵi (x, Ω) P 0 = N ∑ n=1 Φ (n) i (x)Φ (n) k (x 0 ) m n ⏟⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏟ A (n) i 1 𝜔 2 n -Ω 2 + 2j𝜉 n Ω𝜔 n ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ H n (Ω) (19)
where ŵ(x, Ω) is the complex amplitude of w(x, t), such that:

w(x, t) = Re[ŵ(x, Ω)e jΩt ]. (20)
The extension of this modelling methodology to the entirety of a prototype skeleton is obtained by considering more beams and nodes. Appendix A gives as an example the M matrix for a four beams node.

Wing geometry design

In this section, we apply the modelling strategy exposed in the previous section to simple wing geometries for which full parametric studies are at hand, to analyse their frequency response and the frequencies of its eigenmodes. Then, a general optimization strategy is introduced for more complex wing geometries with more than two free geometry parameters.

"L" shaped wing

On a first approach, a simple geometry is selected to illustrate the methodology and to perform a parametric analysis. A wing with a "L" shaped geometry is chosen, as shown in Fig. 5, which could look like an insect wing with its leading edge and chord modeled by equivalent beams (respectively beams 4 and 5). This wing model is completed with the flexible thorax links (beams 2 and 3) and the actuation beam (beam 1), based on prototypes tested in section 6. Since our prototypes and their actuation are symmetrical, we can only model half of them. As depicted in Fig. 5, this symmetrical condition is a guided boundary conditions on beam 1 such as:

⎧ ⎪ ⎨ ⎪ ⎩ w ′ 1 (0) = 0 T 1 (0) = 0 N 1 (0) = 0 (21)
These conditions ensure a free transverse displacement with a locked bending rotation and a free twisting rotation on this end of the beam.

We assume the thorax as a completely rigid element since it is at least three times thicker than any other part of the prototype. The boundary conditions used on the thorax links (beams 2 and 3) are then clamped ones. These conditions ensure no transverse displacement and a locked rotation in bending like in twisting on these beams ends. As seen in the previous section, the chord end (beam 5) is set as a free end. The continuity between the different beams is ensured with the continuity conditions depicted in Eq. (A.12).

The material parameters used are those of the SU-8 polymer with a density of 𝜌 = 1227 kg m -3 , a Young modulus of E = 2 GPa and a Poisson ratio of 𝜈 = 0.22 (ratio used to define Coulomb modulus G = E∕(2(1 + 𝜈)))

Modal analysis

We use a length L = 1 cm, characteristic of our prototypes, which defines the leading edge length l 4 (beam 4), an actuating spacing from the thorax of 5 mm (length of beam 1, l 1 = L∕2), thorax links of 1 mm long (length of beams 2 and 3, l 2,3 = L∕10) and a chord length of 5 mm (length of beam 5, l 5 = L∕2). Square sections of 80 𝜇m × 80 𝜇m for the veins (beams 1, 4 and 5) and 40 𝜇m × 40 𝜇m for the links (beams 2 and 3) are chosen since the fabrication process impose to use parallelepipedic sections.

The mode shapes of this wing structure, obtained with the modal analysis depicted in section 3, are shown in Fig. 6. The structure shows at low frequencies two eigenmodes with mode shapes corresponding to the seeked motions: a flapping mode (Fig. 6a), with leading edge and chord ends in phase, and a twisting mode (Fig. 6b), where these ends are in phase opposition. By observing the higher eigen modes (Fig. 6(c)-6(j)), we can note that these are mostly characterized by an ondulation of the leading edge beam. These deformations of the leading edge render them unusable for our application. We can also notice that for the flapping mode, the twisting along the leading edge stays almost constant while for the twisting mode it increase monotonously by moving away from the thorax links. In addition, no mode made of only pure twisting can be observed. This behaviour seems logical since the first twisting modes of a beam are generally higher in frequency than the bending ones. In our case, from Eq. ( 10) and for a unique beam of length L and thickness h, we can deduce the ratio between the bending and the twisting frequencies as:

𝜔 twist 𝜔 bend = 𝛿 2 𝛾 = √ 6 1 + 𝜈 L h . ( 22
)
Then, since L ≫ h, 𝜔 twist ≫ 𝜔 bend . This modal analysis highlights three remarkable property of this kind of structure. Firstly, there are two modes with mode shapes which reproduce the two required motions: a flapping and a twisting. Secondly, those two modes appear successively in the spectrum and thirdly, they are the two fundamental modes (the two first ones in the spectrum). It is then possible to look for a design that brings those two modes with closer eigenfrequencies, to apply our strategy of mode combination. In the following, we explore this idea successively with a parametric analysis and a numerical optimization.

Parametric analysis

Since the previous modal analysis demonstrates that the two first modes are suitable, we can then narrow our study on the frequency range which contains them. The two selected modes are illustrated in Fig. 6(a) and (b). Notice that the wing membrane is only plotted as an interpolation of the beams positions to ease the understanding of the deformed shapes, reminding that its structural effect is neglected in this study. By looking at the FRF between the leading edge end point and the actuator (Fig. 7(c)) as well as its modal components, we can observe that the two modal frequencies are quite far from each other, with a gap above 100 Hz.

The outcome of this frequency gap is highlighted by the phases given in Fig. 7(d). Indeed, by examining the phase shift (curve (4)) between the flapping (curve (1)) and the twisting (curve (2)) modes, it clearly appears that the quadrature frequencies (𝜋∕2 phase shift) between the two modes are practically the same as their eigen frequencies. By looking at the modes amplitudes on these frequencies (Fig. 7(c)), we also notice that this frequency gap implies that, on each quadrature frequency, one of the modes own by far a larger amplitude than the other one. Therefore, an actuation at these frequencies will result in a deformation in which the mode shape of the corresponding mode will be predominant. Thus the first quadrature will only provide a flapping kinematics with practically no twisting while the second one will produce a twisting motion without any flapping. Consequently, this wing configuration do not allow to combine properly the two modes of flapping and twisting as wanted.

This observation then demonstrates that it is mandatory to bring closer in frequencies the two modes in order to produce a combination of their mode shapes in a quadrature phase shift. Thus it is essential to identify the relevant geometrical parameters of the structure regarding this frequency criterion. For this purpose, as shown in Fig. 8(a), we process to a parametric study in order to observe the influence of two geometrical parameters of the structure: the chord length H and its position along the leading edge x. Links positions and beams sections are kept constant and beams 1,2 and 3 lengths stay the same as in the previous example in order to limit the number of parameters. The frequency evolution of each mode regarding these two parameters is given in Fig. 8(b) and(c). From the parametric surface describing the twisting mode frequency evolution (Fig. 8(c)), we can notice that the position x variation along the leading edge practically do not impact the twisting mode frequency while the chord length H variation affects it in a much more relevant way. Indeed, an increasing chord length is followed by a lowering of the twisting mode frequency regardless of the chord position. These trends can be explained by the monotonous evolution of twisting rotation along the leading edge. Thus, the chord position can not impact the twisting mode frequency since the added inertia on this motion does not vary in regard of the chord position. However since a chord length variation modifies this inertia it will then impact the associated mode frequency.

The flapping mode frequency, given in Fig. 8(b), shows on the other hand a high dependency to the two parameters. While the flapping mode frequency is not that much sensitive to the chord length H when it is located near to the thorax links (i.e. small x∕L ratio), this sensitivity arises with its distance from the thorax. Hence, the raise of the distance chord-thorax (x∕L ratio) is followed by a growing impact of the chord length on the flapping mode frequency, which is then lowering relevantly. These trends could be explained by the transverse displacement rise along the leading edge, in such a way that an inertia located on a small displacement area (i.e. close to thorax) will have a small impact on the frequency while an inertia located on a high displacement area (i.e. on the leading edge end) will get a strong influence on the mode frequency.

As the aim is to get the two modes close in frequency, by raising the flapping mode and lowering the twisting mode, these observations imply the following settlement: getting a wing structure with a chord near to thorax, in order to minimize its length impact on the flapping mode, while keeping a sufficient length to lower the twisting mode frequency. This length should though be less than the leading edge one. Indeed, we can notice in Fig. 8(b) that, for a chord length longer than the leading edge one, the flapping mode frequency will fall off even combined with a chord position close to the thorax.

The calculation of a wing structure in agreement with these criteria, shown as dot circles in Fig. 8(b) and (c), confirms this choice since we can still observe the two modes with the seeked motions, as shown in Fig. 9(a) and (b). By studying the FRF observed on a point at the end of the leading edge (given in Fig. 9(c andd)), it appears that the two modes are closer in frequency as predicted (20 Hz gap against more than 130 Hz in Fig. 9(c andd)). Moreover, as an achievement, the amplitudes of each modes are not anymore neglectible at the quadrature frequencies.

To temper this result, by comparing Fig. 9(b) and (c) with Fig. 7(b) and (c), we can see that the twisting mode show less leading edge bending deformation than in the previous configuration. This deformation loss will then be materialized as a loss in amplitude on the actual prototype. Additionally, it can also be noticed that this wing geometry seems not suitable to produce 4) difference (a. u. means arbitrary units). some lift. Indeed, it does not look like insects ones and does not possess any membrane on the high amplitude area which produce the maximum lift force [START_REF] Wang | Dissecting insect flight[END_REF][START_REF] Lehmann | The mechanisms of lift enhancement in insect flight[END_REF] (i.e. on the leading edge and chord ends). So, even if this wing geometry is adapted to generate the wanted kinematics, this parametric study shows the necessity of a further optimization process according to other criteria than only getting the two studied modes close in frequency.

Numerical optimization

In order to improve the wings kinematics, we choose to process to a parametric optimization on a more complex geometry shown in Fig. 2. This choice is motivated according to our objective which is to trig on the geometrical parameters of the wings, as shown in Fig. 2, which are the position (H 1 ) and length (L) of thorax links, the position (H 2 ), the angles (𝜃 1 , 𝜃 2 and 𝜃 3 ) and length (W 1 , W 2 and W 3 ) of the different veins of the structure. This optimization is based on the minimization of a non-metaheuristic cost function settled accordingly to the criteria chosen to obtain an optimal solution. Indeed, although metaheuristic cost functions allow to obtain a global solution by using algorithms based on physical analogy, as genetic algorithms, they stay relatively complex to implement and present a sizeable calculation time [START_REF] Glover | Handbook of Metaheuristics[END_REF][START_REF] Collette | Optimisation Multiobjectif[END_REF]. Consequently, even if it will provide local solutions depending on the starting point, the choice of this cost function type is motivated by the small calculation time required. To this end, the Matlab optimization function "fminsearch" (which is minimizing a Rosenbrock's function) will be used.

In order to define the adapted cost function, we need to specify the criteria to fulfill. As exposed in section 2 and demonstrated in section 4.2, to combine properly the two flapping and twisting modes, it is necessary to ensure their frequency closeness. Thus, our first criterion is to get the two quadrature frequencies as close as possible. The second criterion is to set a working frequency range for the prototype, with a frequency target around 100 Hz. This criterion also ensures that the algorithm will not provide absurd results such as a null frequency. The third criterion is finally to be able to generate the largest motion amplitude while keeping an equivalent contribution of each mode in the combined kinematics. Indeed, a combined motion with a largely predominant mode would be pointless. It is then mandatory to ensure an equivalent excitability (ratio between the wing and the actuator amplitudes) in order to produce a kinematics composed of flapping as well as twisting. The following cost function is then settled to fulfill these criteria:

C = 𝛼|f target -f quad1 | + 𝛽|f quad2 -f quad1 | + 𝛾 | | | | | max |f (A flap ) -max |f (A twisting ) max(max |f (A flap ), max |f (A twisting ) | | | | | (23)
max |f (•) corresponds to the maximal value of (•) on the frequency range f. The different terms of this function are defined as:

-𝛼, 𝛽 and 𝛾: weight coefficients to emphasize or devalue a criterion regarding the others with respect to its influence on the cost function, those criteria are here taken as equal importance such that the coefficient set 𝛼 = 𝛽 = 𝛾 = 1 is chosen.

-f target : the target frequency for the two modes, it is arbitrarily settled to 100 Hz, ensuring that the algorithm will not converge on a physically absurd frequency such as a null one (Fig. 10(a)). -f quad1 and f quad2 : the two quadrature frequencies, they are deduced from Eq. ( 19) by evaluating the phase shift between the two modal coordinates (arg[H twist (f quad )] -arg[H flap (f quad )] = 𝜋∕2), the objective is to get them close from each other (Fig. 10(b)). It must be noticed that trigging directly on the two modes frequencies to get them close could be harmful in regard of this criterion. Indeed, by looking at the two unimodal phases ( 1) and ( 2) and their difference (3) provided in Fig. 10(c), we can observe that if the two modes are too close in frequency their phase shift will never reach the quadrature value 𝜋∕2 which then lead to the impossibility to combine properly the two modes in the seeked way. The choice is thus made to minimize the frequency gap between the two quadrature frequencies which will ensure the modes frequencies closeness with the warranty of at least one existing quadrature frequency.

-Since these optimization criteria rely on one FRF data, it is essential to choose the used structure point wisely. The chosen point (shown in red in Fig. 2) is located on the leading edge extremity in order to ensure that the flapping and twisting amplitudes are well evaluated by avoiding to observe a vibration node which could distort the amplitudes comparison. A flap and A twist are then the flapping and twisting amplitudes of the respectives unimodal FRFs, they are obtained by considering only the given mode in Eq. ( 19). The term |max |f (A flap ) -max |f (A twist )| is employed to ensure an equivalent contribution of each mode in the wing kinematics by minimizing the gap between the flapping mode amplitude at the resonance max |f (A flap ) and the twisting mode at the resonance max |f (A twist ), as illustrated in Fig. 10(c). This criterion is then a tool to avoid the predominance of one mode in regard of the other one on the quadrature frequencies by enforcing equivalent amplitudes.

Finally, the 1∕

| | | max(max |f (A battement ), max |f (A torsion ) | | |
term allows to maximize the two modes amplification by maximizing their amplitudes at the resonance. This last criterion then warrants the largest amplitude possible during the prototype actuation.

The prototype geometry shown in Fig. 18 has been obtained by this optimization method.

Experimental setup

In order to validate these theoretical conclusions and results, we tested several prototypes and used specific experimental protocols. The fabrication of the prototypes and the experimental setup are described here. The prototypes are fabricated with dedicated MEMS technologies and are composed of a 3D skeleton made of multiple layers of SU-8 photoresist, with thicknesses ranging between 40 𝜇m and 150 𝜇m, and wings membranes, made of a 400 nm thin film of Parylene C, which results in a 20 mm-30 mm wingspan. All the details can be found in Refs. [START_REF] Bao | Design and fabrication of insect-inspired composite wings for MAV application using MEMS technology[END_REF][START_REF] Faux | Couplage modal pour la reproduction de la cinématique d'une aile d'insecte et la génération de portance d'un nano-drone bio-inspiré, french (english title: Modes coupling to reproduce insect wing kinematics and generate lift with a bio-inspired nano-air vehicle[END_REF]. The obtained fully polymeric wings structure can be considered as equivalent to natural insect wings [START_REF] Bao | Improved micromachining of all SU-8 3D structures for a biologically-inspired flying robot[END_REF][START_REF] Bontemps | Design and performance of an insect-inspired nano air vehicle[END_REF][START_REF] Dargent | Micromachining of an SU-8 flapping-wing flying micro-electro-mechanical system[END_REF]. The wings are then driven simultaneously by a single electromagnetic actuator placed at the middle of the thorax and actuated by a sine current. This electromagnetic actuator has the advantage of being easy to manufacture and integrate into the prototypes, having a simple geometry, fast response speed, high bandwidth, and more particularly, the possibility of tuning a wide range of frequencies. The total weight of the prototypes vary from 15 mg to 30 mg, which places them in the insects dimensions range.

The experimental study of these prototypes wings kinematics is then divided in two steps. The first, processed at small amplitude, is used to acquire the FRF on several points of the wings in order to reconstruct their mode shapes. The second one has then for purpose to record and study the high flapping amplitude kinematics.

FRF and mode shapes acquisition

To study the kinematics of a prototype and compare it to the theoretical results, we are looking to measure its structure FRFs and especially of its wings. The purpose of this measurement is then to detect the different existing vibration modes and to observe the quadrature motions. In order to limit the prototype handling, a single measurement bench (shown in Fig. 11) on which the prototype is glued on a cantilever beam has been developed.

The FRFs are measured with a scanning laser vibrometer (Polytec PSV400), which has the advantage of being a non contact measurement method. This device, based on laser interferometry, allows to measure the velocity on a point of a vibrating structure. It is then possible to calculate the FRF on this point by comparing the Fourier decomposition of the velocity signal to the actuation force one. In our case, this force signal is replaced by the supply current in the coil, since this one is linearly related to the magnet displacement and so to the actuation force. From those measurements on different points of the structure, it is then possible to observe the deformed shapes of the prototype wings, since these shapes are deduced from amplitudes and phases values on the measured points for a given frequency.

Given that this measurement is based on the interference of the emitted laser beam with the reflected one, this reflected beam must then be precisely directed on the vibrometer sensor. According to this constraint, it is necessary to use a small flapping amplitude for this measurement. Indeed, a high flapping amplitude will result in a wing angle unsuitable to redirect the reflected beam in the sensor, thus the flapping angle must stay small to enable this measurement.

High amplitude kinematics recording

In order to observe the high amplitude motions, unlike the previous setup, the prototype is actuated at a given frequency with a high intensity sinusoïdal current (300 mA) without damaging the coil. As illustrated in Fig. 12, a speedcam Vision Phantom V 7.4 is used to record the wing motion. This kind of camera allows to record 9152 frames per second with a resolution of 400 × 800 pixels. Since the prototypes working frequency range is between 0 and 300 Hz, we can then obtain at least 30 frames per wings motion cycle and so study its kinematics. As shown in Fig. 12, a lighting system is also added with a reflecting panel to ensure an appropriated light exposure as needed for this kind of recording.

Experimental results and discussion

In this section, we describe the results of both theoretical design and experimental analysis for three selected families of prototypes, of increasing complexity. The aim is to validate the model and optimization results exposed in section 4. The first family of "L wing" prototypes has the simplest shape. Then, the "F wing" family is the result of the parametric optimization of the "L wing", by changing the chord length and position. Finally, the "Wind rider" family is the result of a numerical optimization on nine free parameters.

"L wing" prototypes

Since the main objective of this research is to develop a bio-inspired nano air vehicle, the first prototype features "L" shaped wings with a geometry close to those observed on insects as shown in Fig. 13. As discussed in section 4.2, this type of geom- etry provides effectively a flapping and a twisting modes. However, we also show that the frequencies of these two modes do not allow to combine properly their motions with enough amplitude to reproduce the seeked kinematics, although, it is still necessary to confirm experimentally these theoretical conclusions as well as to validate the wing structure modelling and its hypothesis. For this purpose, we then crafted and tested firstly these "L Wing" prototypes. Results are given in Fig. 14.

As expected, we can observe on the experimental FRF (blue line Fig. 14(c)) the existence of two resonance peaks located at 30 Hz and 100 Hz. The deformed shapes at these frequencies are illustrated in Fig. 14(a) and (b) respectively. We can easily notice that these two deformed shapes are well corresponding to a flapping motion for the first one and to a twisting motion for the second one. Indeed, Fig. 14(a) shows that the leading edge and the chord ends are moving in phase while on Fig. 14( these ends are moving in a phase opposition. This observation then validates the existence of the two wanted motions as well as the structure dynamic behaviour predicted by the modelling. Furthermore, by comparing the experimental (blue line Fig. 14(c)) and theoretical (red line Fig. 14(c)) FRFs, it clearly appears that the resonance frequencies predicted by the modelling and the measured ones are relatively close. Indeed, the flapping mode is theoretically located at 28 Hz against 30 Hz on the measurement while the twisting mode is estimated at 105 Hz by the modelling and measured at 100 Hz. These gaps between theoretical and experimental frequencies can be explained by the fact that the used modelling only takes into account the prototype skeleton dynamic behaviour and neglects the aerodynamic effects linked to the wings membranes. However even if this part, which can seem primordial in the prototype behaviour, is neglected the modelling still demonstrates its reliability for the identification of the two modes. These observations then validate this dynamic modelling for prototypes dimensioning.

Other resonance peaks can also be noticed, with smaller amplitudes, between the flapping and twisting modes (dotted ellipses Fig. 14(c)). These two resonances correspond to the prototype antisymmetrical modes. Their existence is linked to the magnet positioning during the prototype assembling. Indeed, if the magnet is perfectly centered these modes cannot be excited neither observed. If we induce a small dissymmetry in its positioning to produce a dissymmetrical forcing, it will imply the emergence of this supplementary modes in the calculation as well as in measurements. This observation then shows that the magnet is lightly misplaced on the tested prototype. However, it still must be noticed that these antisymmetrical modes amplitudes are far smaller than the symmetrical ones and so that their influence on the prototype dynamic behaviour is neglectible.

The important experimental frequency gap (70 Hz) and the significant amplitude fall between the two modes linked to this gap then confirming that this wing design is not suitable to combine properly the two motions with enough amplitude.

"F wing" prototypes

In order to obtain a suitable mode combination, we choose to get the two modes close in frequency by modifying the wing design. The parametric study on the chord length and position described in section 4.3 leads to a wing design based on an inverted "F" wing geometry with high length chord located near the thorax, as shown in Fig. 15. Like in the previous results, we can notice on an experimental FRF of this new prototype (blue line Fig. 16(c)) the existence of two resonance peaks. Again and as expected these two resonance frequencies correspond to a flapping and a twisting modes of which measured deformed shapes are illustrated in Fig. 16(a) and (b) respectively. As previously, we can observe that in the case of Fig. 16(a) leading and trailing edges ends are well in phase while they are in phase opposition in Fig. 16(b).

By comparing the experimental (blue line Fig. 16(c)) and theoretical (red line Fig. 16(c)) we can notice that, like for the previous prototypes, the modelling allows to estimate accurately the frequency range in which the two modes are located. However, contrarily to the previous case, although theoretical and experimental frequencies of flapping mode are still close (respectively 100 Hz and 105 Hz), the inaccuracy on the twisting mode frequency estimation is on the contrary amplified (estimated at 120 Hz and measured at 145 Hz). These results point out the modelling limits. Indeed this modelling is firstly a pre-design tool but is not made for a precise prediction of the measurement. It is called to suggest the most suitable geometries for a proper modal combination, the experimental tests then allowing to measure amplitudes and adjust frequencies. These errors of estimation are then not harmful.

The main result of this experimental FRF is that the two modes are effectively closer in frequency than in the previous wing configuration. Moreover, the analyse of this FRF shows that the amplitude between the two modes does not fall off as roughly as for the "L Wing" prototype. The essential effect is that the amplitude on the quadrature combination frequencies is then non-neglectible (same order as for resonance frequencies).

With a frequency sweep and by looking at the deformed shapes obtained with the laser vibrometer, it is then possible to determine the frequencies which feature a modal combination. A recording of this combined motion, made with an actuation frequency of 130 Hz, is shown frame by frame in Fig. 17. By tracking the chord motion (dotted orange line) and comparing it with its initial position (dotted blue line), we can observe that the wing (colored in red) kinematics is well composed of a flapping motion providing the lateral displacement (i.e. right to left in sequences 1 to 3 and left to right in sequences 4 to 6) and a twisting motion acting on its slope (i.e. the slope reversal occurs around sequence 4). This recording confirms that the concept of combining these two modes in a quadrature phase shift is exploitable for the purpose of reproducing a kinematics close to insects wings one, although a few restrictions must be noticed. Firstly, the flapping amplitude obtained is still far smaller than insects one since the recorded flapping angle is here of 10• against near 60• for insects. Moreover, this wing design is departing from the ones observed in nature and from an aerodynamic perspective this configuration does not seems adapted to lift generation. Indeed this wing geometry does not feature a lot of membrane surface in the high amplitude area, and so with high velocity, which are the leading and trailing edges ends.

These prototypes then validate experimentally the combination concept from a kinematical point of view. These results demonstrate the feasibility of designing a structure which allows to combine two motions in a quadrature phase shift when vibrating. Hence, the keypoint is to improve these results in terms of amplitudes, frequencies and useful wingspan to produce lift.

"Wind rider" prototypes

A new family of prototypes is proposed here in order to improve the vibratory behaviour of the previous "F wing" prototypes. The new geometry features an added vein and a tilted leading edge which enable to extend noticeably the wing surface as illustrated in Fig. 18. This configuration permits to redistribute the membrane surface on the high velocity area while preserving the seeked dynamic behaviour. Inspired by insects wings specificities [START_REF] Ennos | A comparative study of the flight mechanism of Diptera[END_REF][START_REF] Mountcastle | Wing flexibility enhances load-lifting capacity in bumblebees[END_REF][START_REF] Bethoux | Nomenclature, homology, phylogeny, and the question of the insect wing venation pattern[END_REF], a section variation along the leading edge and intermediate veins are also added. It provides to insects a rigidity modulation along the leading edge and enable to maximize its flapping amplitude. For this "F wing" prototypes, we also leave free up to nine parameters, defined in Fig. 2, and we perform a numerical optimization as described in section 4.4. The purpose of the optimization stays unchanged. We aim to get the flapping and twisting modes frequencies close while maximizing their amplitudes and keeping a working frequency range in the order of 100-150 Hz.

A theoretical FRF of this optimized wing is given in Fig. 19(c) (red line). We can verify with the deformed shapes, illustrated in Fig. 19(a) and (b), that the two measured modes are well corresponding to flapping and twisting modes of the wing. By comparing theoretical and experimental FRFs (respectively red and blue lines in Fig. 19(c)), we can observe that once again the modelling allows to estimate properly the working frequency range of the prototype but cannot predict precisely modes frequencies. As well as in the previous case, this error of estimation is larger on the twisting mode frequency (calculated at 150 Hz and measured at 195 Hz) than on the flapping mode frequency (calculated at 130 Hz and measured at 142 Hz). The expansion of the wing surface with this configuration and thus the accentuation of the aerodynamic effects influence on the prototype behaviour can explain this growing discrepancy between calculation and measurement. In addition, although the gap between the two modes frequencies is of the same order as for the "L Wing" case the motion amplitude between the modes does not fall off. This can be explained by the higher frequency range than in "L Wing" case and by the rise of aerodynamic damping provided by the wing membrane expansion. Furthermore, this discrepancy between the error made on the flapping frequency and on the twisting frequency can be explained by the specific motions induced by these vibration modes. Indeed, since the modelling takes only into account the structural dynamic effects and neglects the aerodynamic interactions, a mode with a high aerodynamic impact will be estimated with a bigger inaccuracy. Here the twisting motion may produce more aerodynamic effects than the flapping one by creating a wide range of unsteady phenomenons.

By analyzing the wing extremity tracking (Fig. 20) and by comparing it with its initial position, we can verify the existence of frequencies for which the two modes are well combined in a quadrature phase shift. This then provides a kinematics involving lateral wing displacement (right to left in sequences 1 to 4 and left to right in sequences 5 to 8) and a wing slope reversal (occurring around sequence 5). It is also noticeable that these frequencies are singular and that other frequencies of the working frequency range do not provide this specific kinematics. Indeed, in this case, a higher actuation frequency will result in a kinematics with a predominant twisting motion since this quadrature frequency is located near the twisting mode one.

As illustrated in Fig. 20, this kind of prototype also induce a larger flapping amplitude than with previous ones, in the order of 30 •. This result is obtained thanks to optimization process supplied with a wider range of geometrical parameters than previously and a wing shape more suitable for aerodynamic interactions. We can notice as well on these sequences recorded with the high-speed camera that the wing membrane swells during the slope reversal (i.e. sequences 1 and 5). This swelling implies also that the wing extremity does not stay straight during the motion, explaining the camber observed on its tracking. This prototype then demonstrates the expected kinematics while featuring a high frequency range (around 150 Hz) and providing a large flapping amplitude.

Conclusions

In this paper, we tested an original concept combination of two vibration modes of artificial wings structure in a quadrature phase shift to reproduce an insect-like wings kinematics. For this purpose, we developed a semi-analytic model to understand and compute the dynamic behaviour of the prototypes. We then demonstrated theoretically the existence of natural vibration modes with mode shapes that correspond to suitable deformations to generate flapping and twisting motions. Furthermore, these two modes are always the two first ones in the eigenspectrum of the prototype wing structures, enabling their combination without interferences with other modes. Thanks to a parametric study of several wing geometries, we confirmed the existence of configurations for which it is possible to get these two modes close in eigenfrequency, the existence of particular frequencies for which a modal combination in a quadrature phase shift is possible and for which a large motion amplitude is obtained, thanks to the resonant motion. From these observations, several criteria according to the performances on these combination frequencies have been established to define an optimal wing geometry.

To confirm these theoretical conclusions, experimental measurements of the dynamic behaviour of different prototypes generations have been processed and compared to the simulation results. Using a laser vibrometer, we validated the predicted mode shapes for small amplitudes, demonstrating the existence of experimental flapping and twisting modes. By comparing measured and calculated FRFs, we also demonstrated that the model is estimating with an acceptable accuracy the working frequency range and modes frequencies of the prototypes. This result is even more valuable in regard with the hypothesis made in the model, i.e. a modelling of the wings skeleton with Euler-Bernoulli beams neglecting the membrane aerodynamic effects. Finally, with a high-speed camera, we verified the existence of frequencies for which the two modes are well combined in a quadrature phase shift and produce a wing kinematics close to insects one. In addition, as shown in Ref. [START_REF] Faux | Two modes resonant combined motion for insect wings kinematics reproduction and lift generation[END_REF], it must be pointed out that the "Windrider" family of prototypes is able to generate a lift force equal to the weight of the prototype, 22 mg, thus validating our concept of mode combination and also the modelling strategy that leads to their optimization and design.

Future works must now focus on the increase of the flapping amplitude in order to enable the take-off of the prototype as discussed in Ref. [START_REF] Faux | Two modes resonant combined motion for insect wings kinematics reproduction and lift generation[END_REF]. To achieve this goal, we must investigate not only on the wing structure but also on the actuation and transmission parts of the prototype. The multi-criteria optimization presented must also be completed with aerodynamic parameters. Moreover, aerodynamic studies, numerical as well as experimental, will be perform to give more understanding of the aerodynamics of our prototypes, with regard to their particular wing kinematics based on the vibration mode combination. The same way, since the moment M toti of a beam i is written as the sum of a bending moment M i along the y i axis and a twisting moment N i along the x i axis, we can write:

⎧ ⎪ ⎨ ⎪ ⎩ M tot 1 = N 1 x 1 + M 1 y 1 M tot 2 = N 2 x 2 + M 2 y 2 M tot 3 = N 3 x 3 + M 3 y 3 (A.7)
By projecting these moments on (x 1 , y 1 ) axes, we have then: As used previously, assembled beams can then be defined by elementary matrix based on the following continuity conditions, which are deduced from Eqs. (2), (A.2), (A.5), (A.6), (A.10) and (A.11):

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ M tot 1 = N 1 x 1 + M 1
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Φ 1 (1)
=Φ 2 (0) Φ 1 (1) =Φ 3 (0) Φ ′ 1 (1) =cos(𝛼 2 )Φ ′ 2 (0) -sin(𝛼 2 )Ψ 2 (0) Φ ′ 1 (1) =cos(𝛼 3 )Φ ′ 3 (0) -sin(𝛼 3 )Ψ 3 (0) Ψ 1 [START_REF] Srygley | Unconventional lift generating mechanisms in free flying butterflies[END_REF] =sin(𝛼 2 )Φ ′ 2 (0) + cos(𝛼 2 )Ψ 2 (0) Ψ 1 [START_REF] Srygley | Unconventional lift generating mechanisms in free flying butterflies[END_REF] =sin(𝛼 3 )Φ ′ 3 (0) + cos(𝛼 3 )Ψ 3 (0) EI 1 Φ ′′ 1 (1) = EI 2 cos(𝛼 2 )Φ ′′ 2 (0) + EI 3 cos(𝛼 3 )Φ ′′ 3 (0) +GJ 2 sin(𝛼 2 )Ψ ′ 2 (0) + GJ 3 sin(𝛼 3 )Ψ ′ 3 (0) GJ 1 Ψ ′ 1 (1) = GJ 2 cos(𝛼 2 )Ψ ′ 2 (0) + GJ 3 cos(𝛼 3 )Ψ ′ 3 (0) -EI 2 sin(𝛼 2 )Φ ′′ 2 (0) -EI 3 sin(𝛼 3 )Φ ′′ 3 (0) EI 1 Φ ′′′ 1 (1) = EI 2 Φ ′′′ 2 (0) + EI 3 Φ ′′′ 3 (0) (A.12)

The first two conditions, based on Eq. (A.2), bind the transverse continuity displacement. The four following ones bind the beams rotations continuity from Eqs. (A.5) And (A.6). The two following ensure the momentum equilibrium from Eq. (A.10). The last one, based on Eq. (A.11), finally impose the shear forces equilibrium.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Illustration of quadrature mode combination concept. a) Example of a prototype design with the wing tracked leading edge, trailing edge and chord. Periodic oscillations of leading (red) and trailing (blue) edges in a quadrature phase shift during flapping motion: (b) leading and trailing edges motion magnitude versus time for a flapping cycle (a. u. means arbitrary units), (c) associated wing motion tracking from a wing tip perspective (wing membrane plotted in orange for the upstroke, left to right half flapping cycle, and in green for the downstroke, right to left half-cycle). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3 .

 3 Fig. 3. (a) Spectrum of the velocity signal measured on one point of the wings, for several air pressures (as specified on the graph) (b) Estimated quality factor and (c) resonance frequency of the low frequency mode as a function of the air pressure.

Fig. 4 .

 4 Fig. 4. Scheme of a "L" shaped wing modeled by a two Euler-Bernoulli beams assembly in bending (w) and twisting (𝜙).

Fig. 5 .

 5 Fig. 5. "L" shaped wing with thorax links and actuation.

Fig. 6 .

 6 Fig. 6. Mode shapes in bending (beams deformations) and twisting (color variations) of the 10 first modes of a "L" shaped wing structure (undeformed structure in dotted lines) (a. u. means arbitrary units). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7 .

 7 Fig. 7. Simulation of a "L" wing structure with thorax links and actuation. a) Flapping and (b twisting mode shapes (undeformed structure in dotted lines). FRF taken on the leading edge end: c) modal coordinates amplitudes in (1) flapping and (2) twisting and their (3) superposition amplitude; d) phases in (1) flapping and (2) twisting and their (4) difference (a. u. means arbitrary units).

Fig. 8 .

 8 Fig. 8. Parametric study processed on "L" shaped wing with thorax links and actuation: a) used geometrical parameters diagram; variations of b) flapping mode frequency, c) twisting mode frequency and their d) difference according to chord length H (ordinate) and position x (abscissa); optimal area in dotted circles.

Fig. 9 .

 9 Fig. 9. Simulation of an "inverted L" wing structure. a) Flapping and (b twisting mode shapes (undeformed structure in dotted lines). FRF taken on the leading edge end: c) modal coordinates amplitudes in (1) flapping and (2) twisting and their (3) superposition amplitude; d) phases in (1) flapping and (2) twisting and their (4) difference (a. u. means arbitrary units).

Fig. 10 .

 10 Fig. 10. Optimization criteria applied on the flapping (1) and twisting modes (2) FRFs and their phase difference (3): (a) proximity with a target frequency f target , (b) closeness of the two quadrature frequencies f quad1 and f quad2 , (c) equivalent amplitudes of the two modes and (d) maximization of the flapping amplitude (a. u. means arbitrary units).

Fig. 11 .

 11 Fig. 11. Experimental setup for FRFs and deformed shapes acquisition.

Fig. 12 .

 12 Fig. 12. Experimental setup for high amplitude wing kinematics recording.

Fig. 13 .

 13 Fig. 13. Scheme and photography of a "L Wing" prototype.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Calculation/measurement comparison for a "L Wing" prototype: a) flapping mode experimental deformed shape, b) twisting mode experimental deformed shape, c) theoretical (red line) and experimental (blue line) FRFs obtained on a point of the leading edge, antisymmetrical modes marked by dotted ellipses (a. u. means arbitrary units). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Calculation/measurement comparison for a "F Wing" prototype: a) flapping mode experimental deformed shape, b) twisting mode experimental deformed shape, c) theoretical (red line) and experimental (blue line) FRFs obtained on a point of the leading edge (a. u. means arbitrary units). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 18 .

 18 Fig. 18. Scheme and photography of a "Wind Rider" prototype.

Fig. 19 .

 19 Fig. 19. Calculation/measurement comparison for a "Wind Rider" prototype: a) flapping mode experimental deformed shape, b) twisting mode experimental deformed shape, c) theoretical (red line) and experimental (blue line) FRFs obtained on a point of the leading edge (a. u. means arbitrary units). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 20 .

 20 Fig. 20. Frame by frame high-speed camera recording for a "Wind Rider" prototype with an actuation frequency of 190 Hz. Wing membrane colored in red, wing extremity initial position in dotted blue line, wing extremity current position in dotted orange line, slope reversal occurring around sequences 1 and 5. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

𝜑 1 =

 1 and the rotation continuity between beams 1 and 3 such as:{ 𝝎 1∕0 = 𝝎 3∕0 ⇒ { 𝜑 3 cos(𝛼 3 ) -𝜃 3 sin(𝛼 3 )𝜃 1 = 𝜑 3 sin(𝛼 3 ) + 𝜃 3 cos(𝛼 3 ) (A.6)

y 1 M 1 (A. 8 ) 3 (N 1 ={ T 1 = T 2 + T 3 (

 1183113 tot 2 = (N 2 cos(𝛼 2 ) -M 2 sin(𝛼 2 ))x 1 +(N 2 sin(𝛼 2 ) + M 2 cos(𝛼 2 ))y 1 M tot 3 = (N 3 cos(𝛼 3 ) -M 3 sin(𝛼 3 ))x 1 +(N 3 sin(𝛼 3 ) + M 3 cos(𝛼 3 ))yThe momentum equilibrium of the junction between the three beams is so expressed such as:{ M tot 1 = M tot 2 + M tot (N 2 cos(𝛼 2 ) -M 2 sin(𝛼 2 )) +(N 3 cos(𝛼 3 ) -M 3 sin(𝛼 3 )) M 1 = (N 2 sin(𝛼 2 ) + M 2 cos(𝛼 2 )) +(N 3 sin(𝛼 3 ) + M 3 cos(𝛼 3 )) (A.10)Finally, the shear forces equilibrium of this three beams junction is written such as:

Acknowledgments

This work was supported by ANR-ASTRID CLEAR-Flight (ANR-13-ASTR-0012), RENATECH program, Direction Generale de l'Armement and Haut-de-France Region. Thierry Tison and the LAMIH Laboratory are warmly thanked for the help in the vacuum chamber experiments. Mathieu Colin is also thanked for his careful reading of the equations.

Appendix A. Structure assembling

In order to extend the modelling methodology to the entirety of the prototype skeleton, we give here an example of elementary beams assembly to model the whole structure. The transverse displacement continuity w between the three beams is then written:

Since the rotation vector 𝝎 i/0 of a beam i can be expressed as the sum of a bending rotation 𝜃 i along the y i axis and a twisting rotation 𝜙 i along the x i axis, we can write:

From the projections (A.1), we have then:

The rotation continuity between beams 1 and 2 is then such as:


Angel M Arevalo-Lopez
email: angel.arevalo-lopez@univ-lille1.fr

Yoshiteru Hosaka

Haichuan Guo

Fabio Denis Romero

Takashi Saito

J Paul Attfield

Yuichi Shimakawa

Spin orders in the charge disproportionated phases of A-site layer ordered triple perovskite LaCa 2 Fe 3 O 9

The coupling between spins and charge disproportionation states has been investigated in the LaCa 2 Fe 3 O 9 oxide with neutron powder diffraction. This A-site layer ordered triple perovskite LaCa 2 Fe 3 O 9 undergoes charge disproportionation on cooling and shows two different charge ordering patterns. At 230 K, Fe 3.67+ disproportionates into a 2:1 ratio of Fe 3+ :Fe 5+ which order in a layered manner along the <010> direction of the pseudocubic unit cell. At lower temperatures (T < 170 K), the charge ordering pattern changes to a layered arrangement along the <111> direction. Neutron powder diffraction data shows that in the intermediate temperature range (170 K < T < 230 K) the spins order into a cycloidal structure on the ac plane for the Fe 3+ cations while the Fe 5+ cations remain paramagnetic. For the

lowest temperature range (2 K < T < 190 K), the spin structure follows the charge ordering and evolves to a <111> layered magnetic structure.

Introduction.

Charge ordering is often observed in transition-metal-oxide systems with electronic instabilities of mixed valent cations. The Verwey transition of Fe 3 O 4 with Fe 2+ /Fe 3+ order is a famous example and the phenomenon is evidenced by a drastic change in transport properties. 1,2,3 Similar charge ordering behaviours are observed in a few oxides containing unusually high valent cations due to charge disproportionation (CD) of Fe 4+ . The perovskite CaFeO 3 shows a metal-insulator transition at 290 K due to CD of Fe 4+ to Fe 3+ and Fe 5+ which order in a rock-salt type arrangement. The CD transition in CaFeO 3 changes the uniform Fe 4+ O 6 octahedra to alternating large Fe 3+ O 6 and small Fe 5+ O 6 octahedra, thus the charge order is coupled with the lattice degree of freedom. 4,5 Other examples of CD in Fe-containing perovskite-structure oxides include rock-salt ordered Fe 3+ /Fe 5+ in the A-site ordered CaCu 3 Fe 4 O 12 and checkerboard ordered Fe 3+ /Fe 5+ in the B-site layer ordered Ca 2 FeMnO 6 . 6,7 The CD in high-valent Fe has been rationalized as the localization of ligand holes which are produced by strong hybridization of Fe d with O p orbitals .

8,9

An important aspect of CD behavior concerns coupling of the charge and spin degrees of freedom. In CaFeO 3 a magnetic transition occurs at 210 K, below the CD transition temperature, and a helical magnetic structure of Fe 3+ /Fe 5+ spins is stabilized. 5 In the A-site ordered quadruple perovskite CaCu 3 Fe 4 O 12 , on the other hand, the magnetic transition temperature is the same as the CD transition temperature, and the Fe 3+ and Fe 5+ spins couple ferromagnetically below the CD transition temperature giving rise to ferrimagnetism when coupled with the Cu ions. 10,11 In the B-site layered double perovskite Ca 2 FeMnO 6 , three-dimensional noncollinear magnetic order develops below 95 K, although the two-dimensional CD of Fe 4+ occurs at 200 K. 7 As can be observed, the coupling between the charge and spins in the CD transition varies significantly between oxides containing high-valent Fe ions.

Very recently we found that a new A-site layered triple perovskite LaCa 2 Fe 3 O 9 shows CD from high-valent Fe 3.67+ to a 2:1 ratio of Fe 3+ and Fe 5+ . 12 Interestingly, the compound shows two successive structural transitions due to different CD ordering patterns. The first transition occurs at 230 K, where the charge disproportionated 2:1 Fe 3+ and Fe 5+ are ordered in a layered manner, along the <010> layer stacking direction of the pseudocubic unit cell.

Below the second transition at 170 K, the charge disproportionated 2:1 Fe 3+ and Fe 5+ are ordered along the <111> direction of the pseudocubic unit cell. Magnetic susceptibility data collected as a function of temperature shows that both transitions occur at the same temperatures as the magnetic anomalies. 12 The simple perovskite La 1/3 Ca 2/3 FeO 3 with identical chemical composition to that of the present LaCa 2 Fe 3 O 9 but with random arrangement of the A-site La and Ca ions presents only one structural transition at 217 K to a 2:1 Fe 3+ and Fe 5+ charge ordering pattern along the <111> direction of the pseudocupic perovskite unit cell, similar to the one found at the lowest temperatures in the A-site layer ordered analogue. 13 The magnetic structures of the CD states have not been determined.

In this study, we analyzed the magnetic structures of this newly discovered layered compound which shows unusual CD behaviors. Using neutron powder diffraction data, we have revealed that the two charge ordering patterns result in different magnetic structures.

The magnetic structures of LaCa 2 Fe 3 O 9 and La 1/3 Ca 2/3 FeO 3 are compared to observe the effects of the cation ordering upon the spin orders.

Experimental.

The synthesis of LaCa 2 Fe 3 O 9 and La 1/3 Ca 2/3 FeO 3 was carried out as reported previously. 12,13 Powder neutron diffraction patterns from both samples were measured using the time-of-flight neutron diffractometer WISH at the ISIS spallation source. 14 Data were collected every 50 K between 2 and 300 K and every 5 K in the 150 K < T < 250 K temperature range on warming. All the data were normalised using the MantidPlot program. 15 Crystal and magnetic structures were refined using Fullprof suite software package. 16 Symmetry analysis was performed with BasIrreps within the same software.

Results and discussion.

Figure 1a shows the neutron diffraction patterns (Bank 2) as a function of temperature from 250 to 150 K collected using the WISH diffractometer on warming. The evolution of two different phases due to the successive transitions at 230 and 170 K reported in our previous study can clearly be seen and the patterns are coloured in red and blue respectively. 12 Neutron diffraction data collected at 300 K was fitted well using the structure reported previously. 12 It consists of a a = a + b, b = 6c, c = -a + b Pnma supercell model (with |a| = |b| = |c| = a p being the lattice parameter of the simple perovskite structure ≈ 3.8 Å), necessary to account for the 2:1 layered A-site ordering between Ca and La and the octahedral tilting. The refined occupancies for the La and Ca are 0.67(1)/0.33 and 0.18(1)/0.82 for the different A sites, in agreement with our SXRD results, thus confirming the A-site layered structure. The structure is reported in Supplementary Table 1. The absence of vacancies at any of the oxygen sites is also confirmed by the neutron refinement.

The structural model accounts for all the observed reflections in data collected from LaCa 2 Fe 3 O 9 at temperatures between 240 K ≤ T ≤ 300 K. Figure 2 top shows the long d-spacing bank of the diffraction data at 240 K; other banks and R-factors are shown in Supplementary Information (Figure S1). Although additional superstructure diffraction peaks appear below 230 K, the nuclear diffraction peaks originating from the Pnma supercell fit well the essential diffraction maxima and show a cell-volume decrease at 230 K and an increase at 190 K (Figure 1b). These temperature dependent changes in the crystal structure are consistent with those observed in SXRD data in our previous study. 12 At least five additional diffraction maxima were clearly observed when the sample was cooled below 230 K as shown in Figure 1a. These features are not present in SXRD data collected from LaCa 2 Fe 3 O 9 in the temperature range 170 K < T < 230 K, suggesting the origin of these peaks is likely magnetic. 12 Mo ssbauer spectroscopy data collected in this temperature range show the presence of Fe 3+ and Fe 5+ in a 2:1 ratio and that the Fe 3+ spins are magnetically ordered while the Fe 5+ spins remain paramagnetic. Therefore, an additional magnetic model was included in the fit to the observed NPD in between 190 K < T < 230 K.

All the additional magnetic reflections can be indexed with propagation vector k 1 = [1/3 1/3 1/3] on a cubic a p × a p × a p cell, or with k 1 = [1/3 0 0] on the Pnma supercell that allows for the layered ordering of La 3+ and Ca 2+ cations. As mentioned earlier, our Mössbauer study shows that the Fe 5+ spins are not long range ordered (idle) in this temperature range; thus, the Fe 3+ spins are ordered below the CD transition temperature, whereas the Fe 5+ spins remain paramagnetic. We performed symmetry analysis, [ 17 ] Irreducible Representations are reported in Supplementary Table 2. The best refined magnetic structure (R mag = 2.9 % for Bank 2) consists of a spin cycloid with 120° rotation along the (100) direction and a magnitude of 1.96(2) µ B per Fe 3+ as shown on Figure 3a. The Fe 3+ magnetic layers are AFM coupled along the b axis and the spins are confined to the ac plane.

Spin density wave-like solutions are also plausible with similar R factors as the cycloid model and cannot be distinguished with powder data, however there are no obvious reasons for the moments to take different values and Mössbauer results support only one magnetic moment value for the Fe 3+ , therefore the spin cycloid was preferred. Fe 5+ spins were constrained to zero values since refined moments where in the limit of detection (≈ 0.4 µ B ) and no significant improvement of the R factors was observed, this is also in accordance to Mössbauer data. After refinement of the magnetic structure, it was analysed with FINDSYM [18], the resulting magnetic space group is Pmc2 1 . It is interesting that both the charge ordering pattern and the magnetic structure show two-dimensional features reflecting the Asite layered crystal structure.

Neutron powder diffraction data collected below T N2 = 190 K, where the second structural transition due to the charge ordering pattern change from the <010> to <111> direction occurs, show additional reflections which increase in intensity down to the lowest temperature measured (1.5 K). These peaks were also assumed to be magnetic since our previous Mössbauer results indicated that both Fe 3+ and Fe 5+ spins are magnetically ordered at these temperatures and no additional reflections were observed in SXRD data. These reflections can be indexed with a propagation vector k 2 = [1/6 1/6 1/6] on the a p × a p × a p cubic cell. To account for the charge ordering pattern, the refinement was performed with an 4) and the best refined low temperature magnetic structure model is shown in Fig. 3b (R mag = 5.9 % for Bank 2). The spins of the charge disproportionated Fe 3+ and Fe 5+ ions, which adopt a 2:1 layered arrangement along the <111> direction of the pseudocubic unit cell, are also ordered in a collinear way. The refined moments at 1.5 K are 3.0(1) and 2.4(1) µ B for Fe 3+ and Fe 5+ respectively. The spins of Fe 3+

and Fe 5+ within each layer are ferromagnetically (FM) ordered and this triple-layered spin structure is AFM stacked along the same <111> direction. Because the spin direction of every triple layer is opposite, no net magnetization is given with this magnetic structure, which is consistent with the observed magnetic susceptibility at low temperatures. 12 The magnetic space group resulted in C2/m within a -a -2b , -a , 2c monoclinic cell.

Temperature dependencies of the magnetic moments are shown on Figure 3c. An important point is that the Fe 3+ /Fe 5+ magnetic structure drastically changes when the CD pattern changes.

As we reported previously, the A-site disordered analogue La 1/3 Ca 2/3 FeO 3 shows the same charge ordering pattern along the <111> direction below 217 K. 13 As shown in Fig. S2b, neutron powder diffraction collected at 1.5 K of this disordered analogue is similar to that observed from LaCa 2 Fe 3 O 9 and the data are well accounted for using an identical magnetic structure (R mag = 4.04% for Bank 2). The refined moments at 1.5 K are 3.

µ B for Fe 3+ and Fe 5+ respectively. Therefore, the ground state CD patterns and magnetic structures are the same for the A-site layer ordered LaCa 2 Fe 3 O 9 and the A-site disordered La 1/3 Ca 2/3 FeO 3 . The results also imply that the ground state magnetic structures are primarily determined by the arrangement of Fe 3+ and Fe 5+ ions and their magnetic interactions, and are less influenced by the A-site order/disorder.

The unusual magnetic structure observed at intermediate temperatures, where only Fe 3+ spins are ordered while the Fe 5+ spins remain paramagnetic, is only stabilized in the Asite layer ordered compound. Taking into account the charge ordering pattern at the intermediate temperatures, which implies a 2:1 order of Fe 3+ :Fe 5+ along the <010> direction; a simple analysis with Goodenough-Kanamori-Andersson rules (GKA) explains the unusual spin cycloid magnetic structure. Considering the layered charge ordering pattern, the strongest interaction along the b axis is AFM between Fe 3+ -Fe 3+ (half-filled e g orbitals) and

frustrates the FM coupling to the Fe 5+ (empty e g orbitals) due to the magnetic periodicity along this direction. A similar magnetic frustration has been recently observed in the Ca 0.5 Bi 0.5 Fe 3.5+ O 3 compound. This material shows CD of Fe 3.5+ to a 2:1 ratio of Fe 3+ : Fe 4.5+ and the magnetic order shows AFM coupling between Fe 3+ spins and idle Fe 4.5+ spins. 19 In conclusion, we have determined the intermediate and ground magnetic states in LaCa 2 Fe 3 O 9 which are both coupled to their respective CD pattern. The intermediate magnetic state shows an unusual spin cycloid with AFM layers of 120º rotating spins of Fe 3+

and idle Fe 5+ spins. The ground state on the contrary shows triple FM layers stacked AFM along the <111> direction. These magnetic structures are explained with the magnetic interactions according to the GKA rule. The disordered La 1/3 Ca 2/3 Fe 3 O 9 presents the same magnetic ground state. The A-site order is responsible for the intermediate phase since this is not observed in the disordered analogue; however, the ground state is very robust as the same magnetic structure is observed in both ordered and disordered compounds. Both Fe 3+ and Fe 5+ order on a triple layered antiferromagnetic structure with a propagation vector k 2 = (0, 0, ½) on a √2a p x √2a p x √3a p hexagonal cell shown with lines. Dashed lines show the 3√2a p x 3√2a p x 6a p supercell needed to account for the CD 2:1 Fe 3+ :Fe 5+ , layered A-site cation ordering and octahedral tilting. c) Magnetic moment as function of temperature obtained from Rietveld refinement on LaCa 2 Fe 3 O 9 and the disordered La 1/3 Ca 2/3 FeO 3 .

  equivalent description in a a = a -b, b = b -c and c = a + b + c, α = β = 90°, γ = 120° hexagonal cell with a propagation vector k 2 = [0 0 ½]. Symmetry analysis was performed (Supplementary Table

Figure 1

 1 Figure 1. a) Thermal evolution of the NPD patterns between 150 and 250 K on warming.

Figure 2 .

 2 Figure 2. LaCa 2 Fe 3 O 9 Rietveld refinement plots of neutron powder-diffraction data collected

Figure 3 .

 3 Figure 3. a) Magnetic structure of LaCa 2 Fe 3 O 9 in the intermediate phase, 170 K < T < 240 K.Fe 5+ spins remain idle while Fe 3+ spins order with a propagation vector k 1 = (1/3, 0, 0) on a

Acknowledgments.

This work was partly supported by

Grants-in-Aid for Scientific Research (Grants 16H00888, and 16H02266) and by a grant for the Integrated Research Consortium on Chemical Sciences from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. The work was also supported by Japan Society for the Promotion of Science (JSPS) Core-to-Core Program (A) Advanced Research Networks and Japan Science and Technology Agency (JST), CREST, and by grants from the EPSRC. Support was also provided by the Royal Society and ICR's Short Term Exchange program. We also thank STFC for the provision of ISIS beamtime, and Dr. Dmitry Khalyavin and Dr. Pascal Manuel for assistance with data collection.