Alexandru-Razvan Luzi

Jean-Marc Biannic

Jean Mignot

Dimitri Peaucelle

Structured adaptive attitude control applied on a Myriade simulation benchmark

This paper considers the problem of reaction-wheel attitude control inside the mission mode of the Myriade satellites. A structured adaptive algorithm, allowing to extend the operating domain of a static proportionalderivative controller is presented and conditions for designing a stabilizing, continuous-time adaptive law are given. In view of implementation, a discrete-time adaptive algorithm is derived and tested on a benchmark of the DEMETER satellite, which was part of the Myriade program. Simulation results show that the structured adaptation and the use of the σ-modification allow the adaptive closed-loop to follow attitude step references of up to 20 degrees without saturating the reaction wheels. This allows the adaptive law to cover the whole mission mode and replace the currently-implemented switched-based control strategy, thus potentially simplifying the verification and validation. process.

I. Introduction

The Centre National dEtudes Spatiales (CNES) Myriade series is a concept of modular microsatellites weighing less than 200 kg that afford quick, low-cost access to space. The aim of Myriade is to bring space research within easier reach of the scientific community, particularly for scientists working on large-scale international missions and facing long lead times between the planning and implementation stages. DEMETER (detection of electromagnetic emissions transmitted from earthquake regions), [START_REF] Buisson | The DEMETER program: a pathfinder to a high performance microsatellite line[END_REF][START_REF] Cussac | The Demeter microsatellite and ground segment[END_REF] the first microsatellite of the Myriade series, was launched in 2004 and was taken off service in 2011. Its complete model of DEMETER is now open for the scientific community [START_REF] Pittet | DEMETER: a benchmark for robust analysis and control of the attitude of flexible microsatellites[END_REF] and will be used in this paper as a benchmark. Two other missions, Parasol and Picard are still orbiting and the two next one Taranis and Microscope will be launched respectively in 2015 and 2016. The evolutions of the Myriade series was intend to enhance the payload capability. The flight domain has been redefined in order to cope with the probable missions over this period: Earth observation and monitoring, sun and space observation (from Low Earth Orbit (LEO)), defense program. These led to specify a range of 500 -800 km altitude, Sun Synchronous Orbit (SSO) orbits for all Local Time at Ascending Node (LTAN). The specified in orbit lifetime is 5 years, with extended qualification to 7 years when possible.

The on-orbit life of the Myriade satellites is decomposed, as for all satellites in general, into phases, allowing the convergence from the separation state to the mission state, the orbit control, and the safe state in case of fault. For each phase or mode, a specific set of attitude pointing requirements is to be met by using a specific set of actuators and sensors. These requirements define the associated Attitude Control mode. Inside a mode, several attitude control and estimation algorithms can be designed, in order to meet increasing performance levels requirements. Over the past thirty years, the number of Attitude and Orbit Control System (AOCS) modes has been reduced, so that generic architectures now typically consist of three modes (safehold, mission and orbit control mode). In the same time the software complexity has increased to meet ever more demanding requirements. Thus, the switching between different AOCS modes requires exhaustive verification and validation in order to demonstrate the robustness of the system with stringent kinematics conditions and switching of attitude control sensors.

In this paper, the problem of attitude control inside the mission mode of the Myriade satellites is considered. Due to the limitations of the reaction wheels used for actuation, the use of a single LTI controller over the whole flight domain covered in the mission mode would not be adequate: the controller would either demand huge actuator solicitation at large pointing errors or generate a very slow closed-loop response. Therefore, the desired closed loop behaviour can only be obtained with a varying controller : when the pointing error is large the controller has to ensure a response that is not over-demanding for the actuators, while a rapid and precise response is required when the pointing error is smaller. The existing controller of the Myriade satellites satisfies this demand by switching between two control laws at large pointing errors a speed tracking loop (momentum bias) is used and at small angles the controller switches to an attitude proportional-derivative law. A satisfying response time and a limited control activity (in particular a limitation of the reaction wheel rotation speed and a null torque for transient regimes, once the travelling speed is reached) are thus ensured. A discontinuity in the control torque can however appear if this travelling speed has not been reached at the switching point. The switch also induces abrupt changes in the closed-loop plant and makes the stability analysis difficult. Based on these remarks, the replacement of the switched-based strategy by a direct adaptive control law is proposed and validated by simulations.

This work was induced by previous studies dealing with the application of robust control techniques to space activities. The suitability of such new techniques has been proven for very various space applications, including the use and validation of robust techniques for the attitude control of Demeter with his flexible appendixes. [START_REF] Pittet | LMI based multi-objective H∞ control of flexible microsatellites[END_REF][START_REF] Peaucelle | Robust analysis of DEMETER benchmark via quadratic separation[END_REF][START_REF] Peaucelle | Simple Adaptive Control without Passivity Assumptions and Experiments on Satellite Attitude Control DEMETER Benchmark[END_REF] This paper considers the direct adaptive control scheme, [START_REF] Fradkov | Adaptive Stabilization of a Linear Dynamic Plant[END_REF][START_REF] Kaufman | Direct adaptive control algorithms[END_REF] consisting of directly tuning the controller gains, based solely on the measured outputs of the system. With respect to more complex estimation/gain scheduling scheme, it has the advantage of simplicity and for this reason it is also referred to as simple adaptive control. [START_REF] Kaufman | Direct adaptive control algorithms[END_REF] As this scheme is often based on passivity properties, it is equally called passivity-based adaptive control. [START_REF] Peaucelle | Passification-Based Adaptive Control of Linear Systems: Robustness Issues[END_REF][START_REF] Fradkov | Passification-based adaptive control with implicit reference model[END_REF] The considered adaptive output feedback law writes, in the general cas, u(t) = K(t)e(t), where e(t) = y(t)-y r (t) is the error with respect to a reference signal y r and where K(t) is the adaptive gain. In this paper, gain and output structuring is considered, such that u(t) = ī i=1 K i (t)e i (t), e i (t) = y i (t)-y ri (t), where each gain component K i (t) is adapted independently, according to a differential equation. This equation contains a gradient-type term -G i ee T i , which drives the gains to stabilizing values, according to passivity properties. The use of different G i matrices is inspired by Ref. 11 and is an extension of the G-passivity property introduced in Ref. 12. The choice of e i and constraints that can be imposed on G i offer important degrees of freedom in designing the adaptive law and prove to be particularly useful for satellite attitude control. The second term of the adaptation equation is a barrier term -Φ(K i -F 0i ), that keeps the gains inside a bounded domain. A term of the type -σ i (K i -F 0i ) (usually called σ-modification [START_REF] Kaufman | Direct adaptive control algorithms[END_REF] ) is also considered. This term is generally used [START_REF] Barkana | Gain conditions and convergence of simple adaptive control[END_REF] for persistent noise compensation as well as for driving the adaptive gains to the "nominal" values F 0i , when the tracking errors become weak.

Structuring the adaptation and the use of the σ-modification allows, for the considered application, obtaining a controller with an enlarged operating domain. In this way a single adaptive algorithm can be applied over the whole mission mode of the Myriade satellites.

This paper is organized as follows. First, a section describes the Myriade attitude control problem and the DEME-TER benchmark. Section 3 presents the method allowing to design stabilizing continuous-time structured adaptive laws, and presents a possibility for deriving a discrete-time law, suitable for implementation on on-board computers. The discrete adaptive law is applied in Section 4 on the DEMETER benchmark and the simulation results are given. Finally, some concluding remarks are driven. Notations R m×p is the set of real m by p matrices. A T is the transpose of the matrix A. Tr(A) is the trace of A. 1 and 0 are respectively the identity and the zero matrices of appropriate dimensions. For symmetric matrices, A > (≥)B means that A -B is positive (semi) definite. For a non symmetric matrix

{A} S = A + A T . ||A|| D = Tr(A T DA) is the weighted Frobenius norm of A with weight D > 0. Φ D,α (R m×p → R m×p ) is a function defined by Φ D,α (K) = ψ α (||K|| 2 D ) • DK, where ψ α ([0 αβ[→ [0 + ∞[
) is a scalar function having the expression of a dead-zone around the origin ψ(0 ≤ k ≤ α) = 0 and that of an increasing barrier going to infinity as k tends to αβ. An example of such a function is

ψ(α ≤ k < αβ) = (k -α)/(αβ -k).

II. Myriade -Attitude control

This section briefly presents the AOCS control loop of the Myriade satellites, with an emphasis on the currently implemented switching-based controller, represented in Fig. 1. Models of the DEMETER satellite (see [START_REF] Pittet | DEMETER: a benchmark for robust analysis and control of the attitude of flexible microsatellites[END_REF] for a full full description of the DEMETER benchmark), which was part of the Myriade program, are used.

J J f J T f 1 δθ η = 0 0 0 -C S δθ η + 0 0 0 -K S δθ η + 1 0 T C (1)
where δθ ∈ R 3 is the vector of angular deviations from a reference δθ = θ -θ r around the axes x, y, z of the satellite; η ∈ R 8 is the state vector associated to the flexible modes (two states per flexible mode). J and J f are respectively the axes inertia and the coupling matrix between δθ and the flexible modes, composing the generalized mass matrix.

Off-diagonal terms of J are in first approximation negligible making the three axes dynamics almost decoupled. The matrices

C s = diag k=1...4 2ζ k ω k 0 0 2ζ k ω k , K s = diag k=1...4 ω 2 k 0 0 ω 2 k
are diagonal generalized matrices of damping and stiffness respectively. [START_REF] Manceaux-Cumer | Minimal LFT form of a spacecraft built up from two bodies[END_REF] T c is the control torque applied to the satellite.

The synthesis models are given by SISO representations of each axis (the inter-axes coupling is thus neglected at the control design stage). The inertias are chosen equal to the nominal value and one flexible mode is considered, defined by ω 1 = 0.4 • 2π rad/s and ζ 1 = 2.7 • 10 -3 (the next flexible modes have less influence and are neglected for the control design). For illustration, the expression of the associated transfer function for the x axis is given:

H sat,x (s) = 0.03933s 2 + 0.0005437s + 0.2485 s 4 + 0.01706s 3 + 7.797s 2 (2)
Once the design of adaptive controllers done on these simplified satellite models, the control law is tested of the 3-axes models with all four flexible modes defined by : (3) Notice that the natural frequency and damping of the first flexible mode is different in simulation from the one used at the design stage. This contributes to testing robustness.

II.B. Actuator and sensor model

The controlled torque for each axis is delivered by accelerating a reaction wheel. Its model is given in Fig. 2

Two saturation blocks can be observed in this model. The one on the rotation speed is critical: if the wheel reaches the maximal speed then the generated torque is zero. The satellite becomes non-controllable, which is critical, as it is open-loop unstable (double integrator). The wheel model used for the controller synthesis is given by H RW (s). The control law is then tested in simulation with the complete model represented in Fig. 2.

The attitude measure is given by a star tracker, modelled in simulation by a 0.45 seconds delay and an additive noise. For the synthesis model, this measure is considered perfect. The delay corresponding to the reaction wheel and the star tracker could be taken into consideration for synthesis by the means of Pade approximations. This is not done here, the synthesis model assuming zero delays. Robustness with respect to these neglected dynamics is then tested by simulation, where all experiments are done including the delays.

II.C. Existing flight software

The attitude control is achieved in a decentralized way, each axis being controlled independently. For each axis, the flight software is composed of an angular velocity estimator and a switching control law followed by a stabilizing linear filter.

The satellite angular velocities are estimated from the angular positions via a high pass filter defined by the transfer function H estim (s) = s 1+0.5s . For each axis, the reaction wheel control system consists of a nonlinear law for large pointing errors which switches to a proportional-derivative controller when the pointing error becomes small:

T ai = -k 0 (δω i + ω d sign(δθ i )), if |δθ i | > θ L (5) T ai = -(F 0θi δθ i + F 0ωi δω i ), if |δθ i | ≤ θ L (6)
where δθ i is the pointing error around axis i ∈ {x, y, z}.

The threshold θ L and speed reference ω d allow to optimize a trade-off between convergence time and control activity. They are chosen the same for all axes. For DEMETER, the numerical values are: ω d = 0.015 deg/s and θ L = 0.3 deg and the scalar parameter k 0 = 1. The gains F 0θi and F 0ωi are chosen so that the produced torque is null, when the reference speed has been reached, and continuous at the switching point. This constraint implies:

F 0θi θ L = F 0ωi ω d (7)
and the numerical values are: F 0θi = F 0θ = 0.1 and F 0ωi = F 0ω = 2 for all three axes i ∈ {x, y, z}. This design generally guarantees the non-saturation of the reaction wheel, that of its speed in particular.

In order to improve stability and to guarantee a good precision level at low pointing errors, a linear filter has been added after the switching structure. [START_REF] Pittet | LMI based multi-objective H∞ control of flexible microsatellites[END_REF] For the x satellite axis, it is given by the transfer function:

H filter,x (s) = 3.039s 2 + 1.457s + 0.09635 0.3333s 4 + 1.371s 3 + 1.263s 2 + 0.4489s (8)
These filters (one per axis) are designed taking into account flexible modes, to satisfy an H 2 /H ∞ trade-off of the linear model close to the equilibrium (|δθ| ≤ θ L ). For small pointing errors (θ L = 0.3 deg), the decentralized control T pi = -H filter,i (s) (F 0θ δθ i + F 0ω δω i ) can therefore be considered as optimal.

The control algorithm runs on the on-board computer at a frequency of 4Hz and discrete representations of the stabilizing filters are obtained from the corresponding continuous-time transfer functions by using the bilinear transformation.

II.D. Towards adaptive control

Even though the existing switching structure answers most control objectives (fair system response time, no actuator saturation, good robustness properties) it has however several drawbacks: the first one regards the possible discontinuity in the computed torque at the switching point. This can happen for instance when switching from the coarse to the fine pointing mode, if the reference speed has not been reached. Also, if this law needs to be modified, several parameters need to be changed at the same time, in order to verify the continuity condition (7). For instance, if the gain F 0θ needs to be changed, then either θ L or ω d have to be modified. Redesign can therefore prove to be delicate.

These are the main reasons for which adaptive control is investigated in this paper. One of the advantages of using such a strategy, is to have a control law covering a larger functioning domain. This would allow replacing the switching-based strategy and potentially simplify the validation process: only one control law would have to be validated instead of two, as in the case of the switching strategy. Adaptive algorithms also have the advantage of gains that vary continuously in time, thus avoiding the switching related problems. Finally changing the parameters of an adaptive algorithm that is unique over the whole operating domain will impact the global behaviour of the system. This would allow, for instance, the improvement of the dynamical performances.

Based on these considerations, an adaptive law based on the proportional-derivative controller is proposed. For each axis i ∈ {x, y, z} this law reads as:

T pi (t) = -H filter,i (s) (K θi (t)δθ i (t) + K ωi (t)δω i (t))
where the gains K i (t) = K θi (t) K ωi (t) are adapted around the values F 0 = F 0θ F 0ω = 0.1 2 known to be optimal for small pointing errors.

When designing the adaptation rules for the gains K θi (t) and K ωi (t), a key issues is to into account the constraints regarding avoidance of reaction wheel saturations, that of its angular speed in particular. Standard approaches to cope with such actuators limitations are based on anti-windup design techniques. [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF][START_REF] Biannic | Optimization and implementation of dynamic anti-windup compensators in aircraft control systems with multiple saturations[END_REF][START_REF] Zaccarian | Modern Anti-windup Synthesis: Control Augmentation for Actuator Saturation[END_REF] The central idea behind these techniques -initially developed in the early 1960's and extensively revisited in the recent literature -consists of synthesizing a controller augmentation that has no effect when saturation does not occur and that otherwise attempts to provide satisfactory performances despite possibly large or fast varying control input signals. The anti-windup device can then be viewed as an additional loop which is designed a posteriori to modify the nominal control laws when saturations appear in the closed-loop system. The adopted strategy is somehow related but indirect. It is based on the following facts.

• At small pointing error the gains F 0 = F 0θ F 0ω = 0.1 2 are chosen for fast compensation of perturbations. If applied at large pointing error (|δθ i | large) the first term gives a value that goes beyond torque saturation (and would hence, also generate the wheel speed saturation). In term of adaptation rules, this indicates that for large pointing error K θi (t) should be decreased such that the possibility of saturating the reaction wheels is reduced. When |δθ i | becomes small, K θi (t) should return the nominal value F 0θ . Note that in order to deal with large pointing errors, K θi (t) should be allowed to take values as small as possible. In other words, the domain inside which the gain is allowed to evolve should be as large as possible.

• With respect to the gains K ωi (t), the physics indicate that globally the speed of reaction wheels follows the angular speed of the satellite. The errors on the angular speed of the satellite should therefore be kept as small as possible to avoid the speed of reactions wheels to saturate. This imposes to choose δω i = ω i -ω ri with ω ri = 0 and have K ωi (t) increase when |δω i | is large. When |δω i | becomes small, K ωi (t) should return to the value F 0ω such that the nominal behaviour is kept.

III. Structured adaptive control

This section gives the theoretical results necessary for designing adaptive laws satisfying the constraints of the considered attitude problem. Conditions for obtaining a stabilizing continuous-time adaptive law are given. The continuous law is then sampled in order to obtain a discrete algorithm suitable for implementation on the satellite computer.

The following systems are considered:

   ẋ = Ax + Bu y = Cx = C T 1 ... C T  T x = y T 1 ... y T  T (9)
where x ∈ R n is the state, u ∈ R m the control input vector and y ∈ R p the measured output, structured in sub-vectors y j ∈ R pj , with y j = C j x, j = 1 . . . .  is the number of sub-vectors in which the output is divided and it can be chosen by the designer. The existence of a stabilizing static output feedback law is assumed:

u = F 0 e =  i=j
F 0j e j (10) where F 0 = F 01 ... F 0 is the stabilizing static output feedback and each F 0j is the gain component corresponding to e j . y r represents the reference, e = y -y r the tracking error and vectors y r and e are structured in a similar manner to the output: e j = y j -y rj , j = 1 . . . . F 0 guarantees stability of the closed loop state matrix:

A(F 0 ) = A + BF 0 C.

III.A. Structured adaptive control

When considering the attitude control problem for one axis models of the satellite, the static law u(t) = F 0 e(t) corresponds to the proportional-derivative controller (6) e = [ δθ δω ] T , which is optimal when the system is closed to the equilibrium point, but has a limited operating domain (|δθ| ≤ 0.3deg). Starting from this static law, u(t) = F 0 e(t), we aim at obtaining an adaptive algorithm, u(t) = K(t)e(t), which could cover an extended operating domain.

The considered adaptive control is of the following type:

   Kj (t) = -G j e(t)e j (t) T Γ j -Φ Dj ,αj (K j (t) -F 0j )Γ j u(t) = K(t)e(t) =  j=1 K j (t)e j (t) , j = 1 . . .  (11)
It is structured as follows. The control action is decomposed in a sum of K j (t)e j (t) adaptive feedbacks with the aim of tuning separately each of the related adaptations laws. This allows to include engineering knowledge about the process to stabilize. It can for example allow to have fast adaptation with respect to errors measured via e 1 and slower adaptation with respect to the other errors. More than just the speed of adaptation (which is tuned via the positive definite gains Γ j ), each adaptation law can be given a 'direction' by the choice of matrices G j . Imposing G 1 e = g 1 e 1 where g 1 is a positive scalar will force K 1 to decrease when e 1 e T 1 is large (go for lower gain when the error is large). Imposing G 1 e = -g 1 e 1 will on the contrary force K 1 to grow.

As in the case of attitude control, for one-axis models, e ≡ [ δθ δω ] T , with e 1 ≡ δθ and e 2 ≡ δω,  = 2, the proposed structure proves to be valuable for this application. Indeed, the requirements regarding the gains evolution, presented in the previous section, can be taken into account as each gain component is adapted separately. By imposing G θ = g θ 0 , g θ > 0 and G ω = 0 g ω , g ω < 0, K 1 (t) ≡ K θ (t) will always decrease when δθ is non-zero and K 2 (t) ≡ K ω (t) will always grow when δω is non-zero. Combined with a modification which will be detailed in Section III.B, this will allow to have the gains evolve as required.

The second term in the adaptation law involving the Φ Dj ,αj function works as a barrier preventing the gains to escape some neighbourhood of the a priori given values of static gain F 0 . More precisely, the following property is proved in. [START_REF] Peaucelle | Robust adaptive L 2 -gain control of polytopic MIMO LTI systems -LMI results[END_REF] Lemma 1. If e(t) is bounded for all t ≥ 0 then ||K j (t) -F 0j || 2 Dj < α j β for all t ≥ 0.

Besides guaranteeing that the gains will always stay bounded, the barrier function plays an important role in proving the stability of the adaptive algorithm. [START_REF] Peaucelle | Simple Adaptive Control without Passivity Assumptions and Experiments on Satellite Attitude Control DEMETER Benchmark[END_REF] Indeed, the stability property is guaranteed by an appropriate choice for the parameters (G j , D j , α j ) j=1... , where G j defines the first term used in the adaptation equation (11) and D j and α j define the barrier function Φ Dj ,αj . The following results give a strategy for finding such values, based on passivity properties of the closed-loop. The proofs are not reproduced here in order to give sufficient space for the application. Some elements for the proofs can be found in. [START_REF] Peaucelle | Simple Adaptive Control without Passivity Assumptions and Experiments on Satellite Attitude Control DEMETER Benchmark[END_REF][START_REF] Peaucelle | Robust Simple Adaptive Control with Relaxed Passivity and PID control of a Helicopter Benchmark[END_REF] See also [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF] for the positive-real lemma on which following proposition is based:

Proposition 1. If ẋ = A(F 0 )
x is asymptotically stable then there exist (P, D j , G j ), j = 1 . . . , solution to the following LMI problem: P > 0,

      {P A(F 0 )} S P B -C T G T 1 • • • P B -C T G T  B T P -G 1 C -D 1 • • • 0 . . . . . . . . . . . . B T P -G C 0 • • • -D        < 0. (12)
Matrices (G j , D j ) then define weighted combinations of inputs w and outputs y = Cx such that the following system is strictly passive:

           ẋ = A(F 0 )x + Bw z =    G 1 . . . G     Cx + 1 2    D 1 . . . D     w (13)
Theorem 1. Let β > 1 be a given scalar. If (F 0 , G j , D j ) j=1... is a feasible solution of the LMI problem (12), then there exist (Q > 0, R j , T j , F j , α j ) j=1... solution to the following LMI problem:

R j QB -C T G T j B T Q -G j C D j ≥ 0, i = 1 . . . ī (14)
T j (F j -F 0j ) T (F j -F 0j ) D -1 j ≥ 0, Tr(T j ) ≤ α j , j = 1 . . .  (15)
{QA(F 0 )} S +  j=1 α j βC T j C j + R j + C T j (F j -F 0j ) T G j C S < 0. ( 16
)
The solution of the LMIs is such that u = F e, where F = F 1 ... F ī , is a stabilizing static output feedback verifying ||F j -F 0j || 2 Dj ≤ α j . Moreover, whatever positive definite Γ j > 0, the parameters (Γ j , G j , α j , D j , β) j=1... define a stabilizing adaptive controller given by (11).

Based on these results, the synthesis procedure works as follows: starting from a known stabilizing static gain F 0 , Proposition 1 combined to any additional linear constraints, allows to design by simple LMI solving some (G j , D j ) j=1... matrices. These matrices being chosen, Theorem 1 outputs the set of scalars (α j ) j=1... which are the last coefficients to be precisely chosen for guaranteeing stability of the adaptive control. The role of the scalars α j is to parametrize the distance from the static gain F 0 that can take the adaptive gain K(t) (see Lemma 1): the larger the α j , the more adaptation is allowed, without braking stability properties of the closed loop.

Since the sets {K j : ||K j -F 0j || 2 Dj ≤ α j β} not only depend of the values α j , but also on the values of D j , a combined search of (D j , α j ) is needed if one wants to maximize the sets of possible excursions of the gains K j around F 0j . The following 2-step heuristic procedure is suggested when performing the design of the adaptive law:

1. Minimize  j=1 w j Tr(D j ) subject to (12) and other possible linear constraints on the (G j , D j ) matrices, were the w j > 0 are positive weights chosen depending on what parameters are expected to wary much or less.

Maximize 

j=1 w j α j subject to (14-16).

Finally, [START_REF] Peaucelle | Robust adaptive L 2 -gain control of polytopic MIMO LTI systems -LMI results[END_REF] suggests taking β closed to 1 in order to limit conservatism of the LMIs and proposes the value β = 1.1. This value is used throughout the results presented here.

The LMIs involved in the design procedure can be solved in polynomial time with efficient semi-definite programming algorithms such as SeDuMi. [START_REF] Sturm | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF] The results presented in the following sections have been coded with Matlab c by using the YALMIP interface. [START_REF] Löfberg | YALMIP : A Toolbox for Modeling and Optimization in MATLAB[END_REF]

III.B. The σ -modification

The σ-modification has been introduced since the first results on adaptive control in order to drive back the control gains to nominal values once the system has stabilized. It also plays a role in compensating persistent noise. For the considered structured adaptive laws, the σ-modification reads as follows:

Kj (t) = -G j ee T j Γ j -σ j (K j (t) -F 0j )Γ j -Φ Dj ,αj (K j (t) -F 0j )Γ j ( 17
)
where the σ j are positive scalars. When neglecting the effects of the other terms, the σ-term drives the control gains to F 0j with first order linear dynamics. The larger is σ j the faster are these dynamics. There exist few formal results for the stability of adaptive control under σ-modification. Convergence of the state to an equilibrium point is hard to prove in general; only results stating the converge to a bounded set are provided in. [START_REF] Kaufman | Direct adaptive control algorithms[END_REF] However, in our considered case, since the static gain F 0 -towards which the σ-term pushes the adaptive gain K(t)is stabilizing, one can expect this modified law to guarantee asymptotic stability, at least locally. More precisely, let us consider the properties of the adaptive law depending on the values of σ j (K j -F 0j ). When this term is negligible stability is then proved by the LMI conditions of Theorem 1. When on the contrary σ j (K j -F 0j ) is dominant then K converges to F 0 , which is a stabilizing static gain. The stability properties of the closed loop are therefore expected to be preserved.

The σ-modification makes the controller evolve between two values: the nominal gain F 0 and K(t) -adaptive gain for large errors. This type of law is therefore useful for the satellite attitude control problem: for weak errors, the nominal behaviour is kept, with the controller varying little around F 0 . When the errors are important, the system drifts from the nominal functioning point and K(t) follows a law which can be approximated by (11), limiting the actuator solicitations. When the errors re-become weak, K returns to F 0 and the nominal behaviour is recovered.

III.C. Discrete implementation of the adaptive law

In view of implementation, the algorithm ( 17) is sampled, based on the Euler transformation (see [START_REF] Ben Yamin | Robust Discrete-Time Simple Adaptive Tracking[END_REF][START_REF] Barkana | On Stability and Gain Convergence in Discrete Simple Adaptive Control[END_REF] ). The barrier function is replaced by a projector whose effect is illustrated, for a two-dimensional case, on Figure 3 The following adaptive algorithm is thus proposed:

     Kj (k) = K j (k -1) -G j e(k)e j (k) T Γ j T s -σ j (K j (k -1) -F 0 )Γ j T s K j (k) = proj Dj ,αj ( Kj (k), F 0 ) u(k) =  j=1 K j (k)e j (k) (18)
where T s is the sampling period.

Remark: If the gain K j is scalar, the projector can be expressed as a saturation function:

proj Dj ,αj ( Kj , F 0j )| Kj ∈R = sat Dj ,αj ( Kj , F 0j ) =          K , if F 0 - αj β Dj ≤ K ≤ F 0 + αj β Dj F 0 - αj β Dj , if K < F 0 - αj β Dj F 0 + αj β Dj , if K > F 0 + αj β Dj (19)
IV. Adaptive control for DEMETER

The adaptive control strategy exposed in Section III is now applied for the satellite. The general methodology is to replace the six static gains (K θi , K ωi ) i∈{x,y,z} with adaptive gains. Preliminary results for only one axis control, with unstructured adaptive control and without σ-modification have been published in. [START_REF] Peaucelle | Simple Adaptive Control without Passivity Assumptions and Experiments on Satellite Attitude Control DEMETER Benchmark[END_REF] Here we consider the 3-axes case and improve the results thanks to the structured adaptive control and do not deteriorate the optimality for low pointing errors with an appropriate usage of the σ-modification.

IV.A. Design of the adaptive control law

The synthesis model of each axis is given by:

1 H estim (s) H sat,i (s)H RW (s)H filter,i (s) , i ∈ {x, y, z} (20)
The state space matrices (A i , B i , C i ) are then obtained and the 2-step design procedure described in Section III.A is applied three times, once for each axis.

For the first step, solving the LMIs ( 12), the constraints G ωi = 0 g ωi , g ωi < 0 and G θi = g θi 0 , g θi > 0 are added such that the gains evolve according to the specifications described in Section II.D.

The continuous-time adaptation law for each scalar gain reads:

Kθi (t) = -g θi δθ 2 i + σ θi (K θi (t) -F 0θ ) + Φ D θi ,α θi (K θi (t) -F 0θ ) Γ θi , i ∈ {x, y, z} Kωi (t) = -g ωi δω 2 i + σ ωi (K ωi (t) -F 0ω ) + Φ D ωi ,α ωi (K ωi (t) -F 0ω ) Γ ωi (21)
This gives the following discrete adaptive algorithm:

Kθi (k) = K θi (k -1) -g θi • δθ 2 i (k) + σ θi (K θi (k -1) -F 0θ ) Γ θi T s (22) Kωi (k) = K ωi (k -1) -g ωi • δω 2 i (k) + σ ωi (K ωi (k -1) -F 0ω ) Γ ωi T s (23)
K θi (k) = sat D θi ,α θi ( Kθi (k), F 0θ ) (24)
K ωi (k) = sat D ωi ,α ωi ( Kωi (k), F 0ω ) (25) T ai (k) = -(K θi (k) • δθ i (k) + K ωi (k) • δω i (k)) (26) T pi (k) = H filter,i (z)T ai (k) (27)
for i ∈ {x, y, z}. The sampling period T s = 0.25s in the case of the Myriade AOCS and H filter,i (z) are the discrete transfer functions obtained from the continuous-time filters H filter,i (s) by applying the bilinear transformation. The control structure of each axis thus contains a proportional-derivative adaptive controller (22-26), with gains adapted around the nominal values F 0θ = 0.1, F 0ω = 2. This adaptive structure is followed by the invariant filter (22). The values of the different parameters, obtained by solving the LMIs (12-16), are summarized in Table 1. The weights w θ = 100, w ω = 1 have been chosen in order to give more importance to large variation domains on the K θi gains. The domains inside which the gains for each axis are allowed to evolve are given in Table 2. There are almost 100% admissible variations of the adapted K θi gains and about 20% admissible variations of the adapted K ωi gains. Another important choice is that of σ θi and σ ωi . The following reasoning illustrates, how this choice is made for the parameters σ θi : suppose a large attitude reference is applied, for instance, on the x axis. After this, the large value of the δθ x term in the adaptation equation (22) causes the gain K θx to decrease and reach the border of the variation domain. At this point its value is constant: K θx = K θx,b = 0.0071, as the saturation function (24) stops the gain from decreasing further. K θx will then leave the domain limit and evolve towards the value F 0θ = 0.1 when the σ-term balances the error term. This occurs when the attitude error reaches a threshold value, δθ thresh,x , verifying:

-g θx • δθ 2 thresh,x = σ θx (K θx,b -F 0θ ) (28)
As g θx , K θx,b and F 0θ are known, the parameter σ θx allows defining the point at which the gain returns to its nominal value. In the following simulations, two sets of numerical values are used: {σ θx = 4.4, σ θy = 14.4, σ θz = 19} and {σ θx = 1.1, σ θy = 3.6, σ θz = 4.8} corresponding to thresholds of 5 and 2.5 degrees respectively for each axis.

The same reasoning applies for choosing the σ ωi values. The choices σ ωx = 5.66 • 10 -4 , σ ωy = 1.65 • 10 -4 and σ ωz = 6.76 • 10 -4 are made, which correspond to a speed threshold of 0.03 deg/s for each axis.

The last parameters to be chosen are the Γ θi and Γ ωi adaptation gains. The simulations are done with the following values Γ θx = 0.15 , Γ θy = 0.05 , Γ θz = 0.03 Γ ωx = 9.7 , Γ ωy = 43.9 , Γ ωz = 10.8

that have been chosen in order to have reasonable values of the gains variation, K θi (k) -K θi (k -1) and K ωi (k) -K ωi (k -1), over the whole functioning domain.

IV.B. Simulation results

Simulation results corresponding to the satellite in closed loop configuration with the adaptive law (22-27) are given next. A 3-axes simulator of the DEMETER AOCS, involving the flexible satellite, actuators and sensors models is considered. The satellite model with the inertias having nominal values and flexible modes defined by ( 3) is first considered. Responses obtained with the adaptive law with two sets of values for the σ θ parameters are compared with results given by the application of the existing control strategy, based on the switching law (5-6), described in Section II.C, Figures 4 and5 show the satellite response to step attitude references of amplitude of 20, 10 and 15 degrees for the x, y and z axis respectively. The first step reference covers the whole domain of the Demeter mission mode. Figure 4 shows the signals corresponding to the x axis: measured satellite angular position and speed, the corresponding adaptive gains as well as the actuator response (angular reaction wheel speed). The same signals corresponding to the y and z axes are plotted in Figure 5.

These results show that the proportional-derivative adaptive law allows tracking of attitude references of large values without saturating the actuators: even if the adaptive law generates larger a control action than the switch-based law, the angular speeds attained by the reaction wheels do not surpass 40 % of the maximal allowed value. The adaptation of the gains thus allows to extend significantly the functioning domain which could be covered by the static proportional-derivative law. The switch-based law could therefore be replaced by the proposed adaptive algorithm: this would also allow a gain in response time, as it can be seen on the attitude plots of Figures 4 and5.

When analysing the plots of the adaptive gains, three phases can be observed: first, after the reference is applied the gains evolve rapidly and reach the barrier of their allowed variation domains. Follows a phase when they remain constant, at values on the barrier (for the x axis, this occurs for K θx between t = 0 and t ≈ 500, and for K ωx between t = 0 and t ≈ 250). Finally, the gains return to their nominal values once that the attitude error and the angular speed respectively are lower than the thresholds defined by the values of the σ parameters.

The influence of these parameters is easy to observe. On the x axis for instance, for σ θx = 4.4, K θx starts evolving towards its nominal value when δθ x = 5deg (t ≈ 500s), while if σ θx = 1.1, K θx returns to its nominal value when δθ x = 2.5deg (t ≈ 700s), as predicted by (28). When K θx starts increasing an acceleration phase can be observed on the attitude plots. The moment when this acceleration appears thus depends on the value of σ θx and this parameter influences in this way the response time. For σ θx = 4.4 the acceleration phase occurs earlier (t ≈ 500s) than for σ θx = 1.1 (t ≈ 700s) and thus the response is faster when σ θx is larger. At the same time, one can also note that the actuator response is also of larger amplitude for larger values of σ θx . This analysis shows that the values of σ influence the closed loop behaviour of the satellite. σ θi in particular can thus be used as design parameters.

Robustness tests

Results of robustness tests of the adaptive law are now presented. First consider the robustness with respect to parametric uncertainties. The previous results considered a simulation model containing four flexible modes, where the satellite inertia J was equal to its nominal value. The case where the satellite inertia J has different values is now considered.

The uncertainty on this parameter was estimated at 20 % of the nominal value and several cases were tested for J having values between J s = 0.8 • J nom and J b = 1.2 • J nom , where J nom is the nominal value of the inertia matrix. The adaptive law is able to handle these uncertainties and the output responses are very close to the ones obtained with the nominal inertia value. As they are very similar to the previous figures, they are not represented here. The actuator responses however are significantly different and are represented in Figure 6 for the adaptive law with {σ θx = 4.4, σ θy = 14.4, σ θz = 19}. Two responses are plotted, for the value of the inertia at the extremity of the uncertainty interval: J = J s = 0.8 • J nom and J = J b = 1.2 • J nom . With respect to the nominal case a smaller control action can be observed for J = J s . On the contrary, when J = J b , the actuator response is of larger amplitude, remaining however far from the saturation value of 293 rad/s. The influence of measurement noise on the adaptive closed-loop system is now analysed. The following simulations consider the case where the satellite is in steady-state and the attitude reference is equal to zero. In the presence of noise the tracking errors for the satellite's attitude and speed are never identically zero. As the error terms used in the adaptation equations (22-23) contain the squares of these errors, the gains are expected to drift in a direction defined by the sign of g θi and g ωi respectively. This drift is balanced by the σ-terms which act as a low-pass filters and push the gains back to their nominal values. Even if the drift is balanced, the equilibrium point K θi = F 0θ , K ωi = F 0ω is lost and an offset with respect to F 0 , along with transmission of the noise, is expected to appear on the adaptive gains.

This behaviour is illustrated in Figure 9, where the adaptive gains for the z axis are plotted, similar responses being obtained for the other axes. The gains stay very close to their nominal values, {0.1, 2} and as the measurement noise is augmented by the speed estimator (high-pass filter), K ω suffers more from the disturbance. The offset on this gain is however of less than 0.1% with respect to its nominal value. The computed torque T az obtained in steady-state with the switch-based law and with the adaptive law respectively is shown in Figure 8. Note that for the switch-based law only the proportional-derivative branch (6) having constant gains is active, as the pointing error is much smaller than the threshold θ L = 0.3deg. The influence of noise can be seen on both control signals and it can be noticed that the noise levels are the same in both cases. Thus, even though the adaptive gains are influenced by the measurement noise, the impact on the computed control signal is negligible with respect to the case where the gains are constant. The optimal behaviour guaranteed at low angles by the static controller is thus preserved with the adaptive law. The new adaptive algorithm has a discrete-time expression and uses only the measured outputs of the system, not needing any additional excitation. The required computational effort is limited,making the adaptive algorithm suitable for implementation on on-board satellite computers. The simulation results confirmed the interest of using an adaptive control law: step references of up to 20 degrees can be tracked without engendering the saturation of the reaction-wheel, while model uncertainties as well as persisting noise can be coped with. Finally, the adaptive algorithm disposes of several design parameters and a connection between the values of σ θi and the response time for each axis was shown.

Based on these preliminary positive results, in-flight tests on the Picard satellite are being considered.

Figure 1 .

 1 Figure 1. Myriade AOCS loop

ω 1 ω

 1 2 ω 3 ω 4 = 0.2 0.3 0.45 0.6 • 2π rad/s ζ 1 ζ 2 ζ 3 ζ 4 = 30 50 5 25 • 10 -4 .

  3 :

Figure 2 .

 2 Figure 2. Reaction wheel model

  : if the computed gain Kj is such that || Kj -F 0j || 2Dj ≥ α j β, then Kj is projected on the domain border ||K j -F 0j || 2 Dj = α j β, along the line connecting Kj to F 0j .

Figure 3 .

 3 Figure 3. Projector function

Figure 4 .Figure 5 .

 45 Figure 4. Tracking of attitude references -X axis

Figure 6 .

 6 Figure 6. Reaction wheels response for different values of the inertia

Figure 7 .

 7 Figure 7. Influence of measurement noise on the adaptive gains -Z axis

Figure 8 .

 8 Figure 8. Influence of measurement noise on computed torque -adaptive and switch-based law -Z axis

Figure 9 .

 9 Figure 9. Influence of measurement noise on computed torque -adaptive and switch-based law -Z axis

Table 1 .

 1 LMI solving results

	g θi	g ωi	D θi	D ωi	α θi	α ωi
	x axis 53.52 -941.44 1135.46 9683.27	8.9	1831
	y axis 175.71 -208.04 3615.9	30954	28.56 3370.7
	z axis 231.04 -843.44 4754.61 40796.6 36.17 4338.59

Table 2 .

 2 LMI solving results

	x axis	y axis	z axis
	K θi [ 0.0071 0.1929 ] [ 0.0072 0.1927 ] [ 0.0072 0.1927 ]
	K ωi [ 1.5439 2.4561 ] [ 1.6580 2.3927 ] [ 1.6539 2.3461 ]

of 16 American Institute of Aeronautics and Astronautics

of 16 American Institute of Aeronautics and Astronautics

of 16 American Institute of Aeronautics and Astronautics


Nicolas Jeannin

Maura Outeiral

Laurent Castanet

Carlos Pereira

Danielle Vanhoenacker-Janvier

Carlo Riva

Carlo Capsoni

Lorenzo Luini

Mario Cossu

Antonio Martellucci

Atmospheric Channel Simulator for the Simulation of Propagation Impairments for Ka Band Data Downlink

Keywords: Earth observation data downlink, ka band, attenuation ground station, rain, numerical weather forecast models. I

This paper describes a methodology to model the propagation conditions for Earth observation data downlink operating at Ka band. It relies on the use of numerical weather forecast models to perform local high resolution reanalysis of the meteorological conditions on which the propagation effects can be computed. From the meteorological simulations spanning long durations, time series representative of attenuation between an orbiting satellite and a ground station are extracted, knowing orbital and RF characteristics of the system.

INTRODUCTION

The current X band data downlink of Earth Observation (EO) satellites is a bottleneck for the always increasing data acquisition rate, driven mainly by the increase in sensors resolution and swath. The lack of bandwidth dedicated to those kinds of links at X band motivates the use of a portion of Ka band for the generation of observation satellite currently in development both in US and in Europe where a much larger bandwidth is available [START_REF] Cossu | VCM performance for Ka-band downlink with multiple ground stations[END_REF][START_REF] Toptsidis | Link adaptation for Ka band low Earth orbit Earth Observation systems: A realistic performance assessment[END_REF][START_REF] Arapoglou | A tool for synthesizing rain attenuation time series in LEO Earth Observation satellite downlinks at Ka band‖[END_REF].

One of the technical difficulties raised by the use of this frequency bands is induced by the higher propagation impairments induced by the troposphere. An accurate assessment of the impairments affecting those links is essential for a proper design of the physical layer and for the performance evaluation of the system. So far there exists only a limited collection of data representative of this configuration [START_REF] Kostulski | Variable slant-path Ka-band propagation measurements on the Australian LEO microsatellite ‗FedSat[END_REF][START_REF] Chen | ROCSAT-1 Ka-band experimental communication payload operations[END_REF]. Though, on the contrary, the propagation channel at Ka band has been characterized for a relatively long time for geostationary satellites, models and measurements collected in this framework can hardly be employed to investigate the performance of high-frequency EO transmission systems due to the different geometry of the Earth-space link. Especially an accurate representation of the propagation channel at low elevation angle and of the coupling between the evolution of the propagation conditions around the receiving ground station and the evolution of the link seems of prime significance. This requires a depiction of the meteorological conditions in space and in time around the ground station.

All the atmospheric parameters needed to describe the propagation conditions can not be quantified on a generic basis with remote sensing data. Thus, to characterize the channel state for this kind of links, a solution relying on the use of non hydrostatic meso-scale weather forecast models has been developed. The main principles and preliminary results from the simulator are described in this paper.

II. DESCRIPTION OF THE SIMULATOR

The aim of the simulator developed is to perform local reanalysis of past weather conditions using WRF (Weather Research and Forecasting) model and then to turn them into propagation parameters for a LEO to ground link [START_REF] Garcia | Use of WRF model to characterize propagation effects in the troposphere[END_REF][START_REF] Skamarock | Coauthors 2008: A description of the Advanced Research WRF version 3[END_REF]. The main algorithmic steps performed by the simulators are described on Figure 1 The first step is to define the period of simulation and position of the ground station. WRF model is initialized using given boundary conditions from ERA-interim reanalysis for the area around the ground station and the considered period. Three nested domains with a finer resolution are used by the model as illustrated on Figure 2. The final resolution of the output is of 2 km in a square of 180 km around the position of the ground station. To get a temporal depiction of the channel state compatible with the available data storage and the time step of the model used for microphysic computation, data are generated every 5 min. Generated data include among other profiles of temperature, pressure, specific humidity, wind and hydrometeor content on model levels. Temperature and pressure profile are needed to compute oxygen attenuation while the water vapour attenuation requires also the water vapour content. Among the hydrometeors related variables considered in the microphysic scheme used in WRF (single moment 6 classes), rain specific content and cloud liquid water content are used to account respectively for rain attenuation and cloud attenuation. The way to convert the various meteorological fields at a given time into specific attenuation can be found in [START_REF] Garcia | Use of WRF model to characterize propagation effects in the troposphere[END_REF]. Unlike the attenuation computation algorithm from numerical weather prediction models proposed in [START_REF] Hodges | An attenuation time series model for propagation forecasting[END_REF], the specific attenuation is directly computed using the hydrometeor content given by the model, using a lookup table containing the specific attenuation function of the specific rain water content and of the temperature, assuming Mie scattering.

A methodology to compute the scintillation characteristics from the output of the models is also described in [START_REF] Pereira | Simulation of tropospheric scintillation on LEO satellite link based on space-time channel modeling‖[END_REF]. The specific attenuation fields are then turned into 2D fields representing the azimuth and elevation for each pointing direction of the antenna. Knowing the orbitography of the satellite of interest and the position of the ground station, the time series of elevation and azimuth can be computed. The attenuation time series are extracted by a tri-linear interpolation in the time, azimuth, elevation domain of the generated attenuation fields as illustrated on Figure 3.

III. PRELIMINARY RESULTS AND VALIDATION METHODOLOGY

The simulator described in the previous section has been used in order to generate time series of total attenuation for long duration on locations that could be foreseen to host a ground station of future Earth observation satellites (mainly in high latitude areas to maximize the contact time with quasi polar satellites). An example is shown on Figure 4 for a ground station located in Fairbanks (64°50'N, 147°42'W). It has to be kept in mind that the visibility time for satellites at low altitude represents not more than 10% of the time (maximum 16 min, each epoch for a LEO satellite at 800 km of altitude). Thus the uncertainty associated to conditional statistics on the elevation angle for instance is high due to this short duration. To increase the duration of the time-series in the analysis, the inclusion of time series generated with slightly different orbital parameters (as the longitude ascending node) but keeping the same azimuth elevation distribution can be interesting. It could allow increasing the duration of the analysis for a duration spanning the one of the meteorological simulations.

To get an idea of the accuracy of the methodology, statistics generated with the simulator can be compared to existing models and experimental results. An example of the statistics obtained for various elevation angles from one year of simulations are compared on Figure 5 with the complementary cumulative distribution function given by ITU-R Rec P.618-10. For most elevation, angles the statistics extracted from the simulation are comparable to ITU-R model (it has to be noticed that the year 2012 was very close to the long term average in term of rainfall amount in Fairbanks).

The ability of the simulator to reproduce results from past propagation experiment is also under investigation. It has been parameterized to reproduce the measurements made in Spino d'Adda (Italy) during the ITALSAT campaign. Outputs of the simulator over Spino d'Adda for 1996 for the same azimuth, elevation, polarization and frequency as the experimental link have been compared to the experimental results.

A comparison of the statistics collected during the ITALSAT experiment with the statistics computed from the simulator is illustrated on Figure 6. Figure 6 shows that the total attenuation statistics at 18.7 GHz computed from the simulated time series are comparable to the ones computed from the time series compute from ITALSAT experiment. It tends to imply that the outputs of the simulator are statistically representative of the attenuation in this area.

A scatter plot illustrating the dispersion between the simulated attenuation values to the concurrent experimental values is illustrated on Figure 7. An extremely large dispersion can be observed on the scatter plot with a linear correlation coefficient slightly higher than 0.2. The results of Figure 7 are not excessively surprising as meteorological models fail often to get a good representation of the conditions at a particular time and a particular location.

It is however anticipated that the structure of the precipitations, clouds and water vapor will be representative of the one observed over the location of interest due to the physical modeling of those processes.

IV. CONCLUSION AND WAY FORWARD

This paper has presented a methodology to simulate propagation conditions from the output of high resolution nonhydrostatic meteorological forecast models to supplant the lack of data for the sizing of LEO to ground links.

The outputs from the meteorological models are used to describe the propagation medium, and classical propagation computation methods are used to account for the attenuation losses. The main interest of this methodology is to have a simultaneous description of all the atmospheric parameters at a sufficient resolution to compute the propagation effects directly from them. The description of the meteorological parameters is supposed to obey to local climate specificities due to the inclusion of a high resolution land-cover and terrain model in the model. Furthermore seasonal or inter-annual trends are reproduced using this methodology as intrinsically contained in the large scale input to the model. With this methodology attenuation fields with a resolution around 1 km and a temporal resolution of 5 min can be generated. This methodology could be applicable to various kinds of problems and may possibly act as fictitious experimental data for some problems (gateway margin optimization, ground to UAV links…), but need first to be further validated.

Validation is on-going and expected to give bounds on the current accuracy of such an approach. Encouraging results have been obtained by comparing statistical results from the ITALSAT campaign with simulations made so that it reproduces the experimental configuration. Simulator's outputs have been found to be statistically representative of ITALSAT measurements for one year. Further test on other propagation experiments in various places and for longer durations will be carried out in a near future. A particular emphasis will be put on areas of interest for future ground station for the reception of Earth observation data download link.

The presented simulator shows obvious limitations in terms of ease of usability and of IT infrastructure needed to get the results as for instance the simulation of one year of data requires more than 3To of memory and approximately 400 CPU days combining simulation and post-processing. More flexible and generic approach could be of interest in order to get a channel description at a lesser computational cost. Those models will be developed in the next phase of the study from a statistical analysis of the outputs of the simulator.

Figure 1 :

 1 Figure 1: Flow chart of the simulator

Figure 2 :

 2 Figure 2: Example of simulation domains defined to reach a resolution of 2 km from ERA-interim reanalysis data.

Figure 3 :

 3 Figure 3: Example of azimuth elevation gridded attenuation fields and extraction of the time series for three consecutive visibility periods.

Figure 4 :

 4 Figure 4: Example of simulated attenuation time series at 26GHz between a satellite with a METOP like orbit and a ground station in Fairbanks. Left picture: one year time series. Right picture: zoom on specific visibility period.

Figure 5 :

 5 Figure 5: Comparison of ITU Rec P.618 for Fairbanks for links at 26GHz at various elevations with statistics computed from simulations from 2012

Figure 6 :

 6 Figure 6: Comparison of the empirical distributions computed from ITALSAT data and from the simulator on the same configuration

Figure 7 :

 7 Figure 7: Sample by sample comparison of the data generated by the simulator with measured ITALSAT data

  .

	S e le ctio n o f
	a g ro u n d	C o a rse
	sta tio n a n d o f sim u la tio n d a te	m e tre o ro lo g ica l d a ta (E R A -in te rim )	syste m W R F p re -p ro ce ssin g	B o u n d a ry co n d itio n s fo r W R F
		S ta tic d a ta
		(a ltitu d e , la n d -u se … )		C o m p u ta tio n w ith W R F
				H ig h re so lu tio n
				m e te o ro lo g ica l d a ta
				(kilo m e tric )
				P o st p ro ce ssin g a lg o rith m s	E M ch a ra cte ristics
					o f th e lin k
				P ro p a g a tio n
				p a ra m e te rs fie ld s
				T im e se rie s e xtra ctio n a lg o rith m s	O rb ito g ra p h y
				P ro p a g a tio n tim e
				se rie s fo r L E O to
				g ro u n d lin ks

ACKNOWLEDGMENT

The authors are grateful to ESA for funding partly this activity through the contract ESTEC 4000105326/12/NL/LvH "Channel modelling for design of Earth observation Ka-band data downlink systems".


Matthieu Bellec

Ulrich Kuhl

Gilles Montambaux

Fabrice Mortessagne
email: fabrice.mortessagne@unice.fr

Manipulation of edge states in microwave artificial graphene

Keywords: numbers: 42.70.Qs, 71.20.-b, 73.22.Pr, 73.20.r

Edge states are one important ingredient to understand transport properties of graphene nanoribbons. We study experimentally the existence and the internal structure of edge states under uniaxial strain of the three main edges: zigzag, bearded, and armchair. The experiments are performed on artificial microwave graphene flakes, where the wavefunctions are obtained by direct imaging. We show that uniaxial strain can be used to manipulate the edge states: A single parameter controls their existence and their spatial extension into the ribbon. By combining tightbinding approach and topological arguments, we provide accurate description of our experimental findings. A new type of zero-energy states appearing at the intersection of two edges, namely corner states, are also observed and discussed.

Introduction

Edge states play an essential role in condensed matter physics for both fundamental aspects and electronic transport applications. Recently, due to their immunity against disorder or impurities (absence of backscattering), topologically protected edge states have raised interest in quantum spin Hall systems [START_REF] Kane | Quantum spin Hall effect in graphene[END_REF][START_REF] Bernevig | Quantum spin Hall effect and topological phase transition in HgTe quantum wells[END_REF][START_REF] Bernevig | Quantum spin Hall effect[END_REF] as well as in topological insulators [START_REF] Haldane | Model for a quantum Hall effect without Landau levels: condensed-matter realization of the "parity anomaly[END_REF][START_REF] Hasan | Topological insulators[END_REF]. Edge states have been first predicted in graphene ribbons [START_REF] Nakada | Edge state in graphene ribbons: Nanometer size effect and edge shape dependence[END_REF][START_REF] Brey | Electronic states of graphene nanoribbons[END_REF][START_REF] Kohmoto | Zero modes and edge states of the honeycomb lattice[END_REF] and later observed along "zigzag" edges [START_REF] Kobayashi | Observation of zigzag and armchair edges of graphite using scanning tunneling microscopy and spectroscopy[END_REF][START_REF] Niimi | Scanning tunneling microscopy and spectroscopy of the electronic local density of states of graphite surfaces near monoatomic step edges[END_REF]. Although not strictly speaking topologically protected, edge states in graphene possess a topological origin coined by the Zak phase and remain robust against weak chiral symmetry perturbations [START_REF] Ryu | Topological Origin of Zero-Energy Edge States in Particle-Hole Symmetric Systems[END_REF][START_REF] Delplace | Zak phase and the existence of edge states in graphene[END_REF].

These peculiar features are due to the multicomponent (spinorial) structure of the wavefunction. While for scalar wavefunctions (like for free particles in a box), hard wall boundary conditions impose the vanishing of the wavefunction at the edges, for a two-component wavefunction only one component has to vanish, leaving the possibility of a finite amplitude for the other component. In graphene, the two components are simply the amplitudes of the wavefunction on each of the two atoms of the bipartite honeycomb lattice.

The richness of edge-state physics is not limited to condensed matter. Any finite system characterized by multi-component wavefunctions can constitute a good candidate. Over the past few years, pertinent realizations of artificial graphene have emerged in various contexts such as 2D electron gases in molecular assemblies [START_REF] Gomes | Designer Dirac fermions and topological phases in molecular graphene[END_REF] or in nanopatterned semiconductors [START_REF] Singha | Two-dimensional Mott-Hubbard electrons in an artificial honeycomb lattice[END_REF], ultracold atoms in optical lattices [START_REF] Tarruell | Creating, moving and merging Dirac points with a Fermi gas in a tunable honeycomb lattice[END_REF], polaritons in semiconductor microcavities [START_REF] Jacqmin | Direct observation of Dirac cones and a flatband in a honeycomb lattice for polaritons[END_REF] as well as light and microwaves in photonic crystals [START_REF] Rechtsman | Strain-induced pseudomagnetic field and photonic Landau levels in dielectric structures[END_REF][START_REF] Bittner | Observation of a Dirac point in microwave experiments with a photonic crystal modeling graphene[END_REF][START_REF] Kuhl | Dirac point and edge states in a microwave realization of tight-binding graphene-like structures[END_REF][START_REF] Bellec | Topological transition of Dirac points in a microwave experiment[END_REF] (see [START_REF] Polini | Artificial honeycomb lattices for electrons, atoms and photons[END_REF] for a recent review). The main advantages of these analogue systems lie in their high tunability and the control of their lattice properties. The synthetic honeycomb lattices offer the possibility to investigate phenomena that are hardly reachable in genuine graphene, and more particularly those appearing at the edges. For instance, edge states at bearded terminations (not stable in real graphene) have been first observed in photonic lattices [START_REF] Plotnik | Observation of unconventional edge states in "photonic graphene[END_REF]. As pointed out in [START_REF] Hafezi | Robust optical delay lines with topological protection[END_REF][START_REF] Fang | Realizing effective magnetic field for photons by controlling the phase of dynamic modulation[END_REF], the manipulation and control of edge states may lead to promising photonic applications.

Moreover, graphene band structure, and consequently transport properties, can be engineered via lattice strain [START_REF] Guinea | Energy gaps and a zero-field quantum Hall effect in graphene by strain engineering[END_REF][START_REF] Guinea | Strain engineering in graphene[END_REF][START_REF] Wunsch | Dirac-point engineering and topological phase transitions in honeycomb optical lattices[END_REF][START_REF] Montambaux | Merging of Dirac points in a two-dimensional crystal[END_REF][START_REF] Pereira | Tight-binding approach to uniaxial strain in graphene[END_REF]. While not observable in the realm of genuine graphene, bandgap opening and associated topological phase transition, in uniaxial strain honeycomb structures, as well as pseudo-magnetic field creation via inhomogeneous strain, have been observed in various artificial systems [START_REF] Gomes | Designer Dirac fermions and topological phases in molecular graphene[END_REF][START_REF] Tarruell | Creating, moving and merging Dirac points with a Fermi gas in a tunable honeycomb lattice[END_REF][START_REF] Rechtsman | Strain-induced pseudomagnetic field and photonic Landau levels in dielectric structures[END_REF][START_REF] Bellec | Topological transition of Dirac points in a microwave experiment[END_REF][START_REF] Rechtsman | Topological creation and destruction of edge states in photonic graphene[END_REF].

In this paper, we propose an experimental manipulation of edge state properties by controlling a uniaxial strain. We use a photonic implementation of honeycomb lattice in the microwave regime [START_REF] Bellec | Topological transition of Dirac points in a microwave experiment[END_REF][START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF]. Fig. 1(a) shows a realization of an artificial graphene ribbon exhibiting the three usual types of edges, namely zigzag, bearded and armchair. The sites of the lattice are occupied by dielectric microwave resonators with a cylindrical shape (diameter of 8 mm, height of 5 mm). The resonance frequency of an isolated resonator ν 0 is around 6.65 GHz and corresponds to the on-site energy of atoms The experimental setup and the tight-binding description of the microwave artificial graphene are detailed in [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF].

The paper is organized as follows. In Section 2, we first focus on zigzag and bearded boundary geometries. We show experimentally how uniaxial strain acts as a switch between zigzag and bearded edge states. Based on a tight-binding analysis, a diagram of existence of edge states is theoretically proposed. We recall in Section 3 the topological origin of the three types of considered edge states, namely zigzag, bearded and armchair. A geometrical analysis in the k-space allows to predict the presence of edge states and their evolution under strain. Section 4 is dedicated to armchair geometries. A quantitative experimental and theoretical analysis is done. The existence of a new type of states, appearing at the intersection of two type of edges, namely corner states, is eventually discussed.

Zigzag and bearded edges in honeycomb lattice under uniaxial strain

The lattice presented in Fig. 1(a) exhibits three different edges: Armchair, zigzag and bearded. We will consider ribbons uniaxially strained along one lattice axis (horizontal direction in Figs 2 and6) where the strain changes one of the three nearestneighbor couplings only. The modified coupling is denoted by t ′ and the anisotropy parameter by β = t ′ /t. Armchair edges are along the strain axis and consequently will not support any edge state whatever the anisotropy as will be discussed in Sec. 3. Fig. 1(b) shows a typical DOS measured in an unstrained ribbon, i.e. β = 1. The Dirac frequency ν D is obtained by following the procedure described in Ref. [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF] and defines the frequency origin. The peak observed at the origin corresponds to "zeroenergy" modes in the condensed-matter context, we will call them "zero-modes" in the following. Experimentally, we can extract the intensity distributions by means of reflection measurements (see [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF] for details). Fig. 2 shows the intensities of the wavefunctions associated to zero-modes for different values of β: The zero-modes are all located along edges. In case of the unstrained lattice β = 1, Fig. 2(b), the intensity is clearly distributed along both zigzag and bearded boundaries. Then, the anisotropy parameter β controls the relative weight between the two types of zero-modes. For β = 0.4 [Fig. 2(a)], bearded edges are dominant whereas they are totally absent for β = 2.5 where only the zigzag edge is illuminated [Fig. 2(d)]. For intermediate anisotropy [β = 1.5, Fig. 2(c)], both edge types are excited, with an opposite relative strength compared to the β = 1 case.

To go beyond this qualitative discussion, we propose in the following part a theoretical description by means of a tight-binding model. Although both zigzag and bearded edges are built experimentally on the same ribbon sample, we will consider theoretically two independent semi-infinite lattices. We demonstrated in a previous work the validity of the TB model to describe our artificial graphene realization [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF]. Here, to address the issue of edge states, we will restrict the TB model to first nearest neighbor couplings. Let us first consider a zigzag edge [Fig. 3(a)]. The A-B unit cell, appropriated to describe such an edge, is represented by the dashed box. The edge is built by translating the dimer with vectors m a 2 (m integer). According to the Bloch theorem, an additional phase e imk is acquired after m translations. Here, k ≡ √ 3k a, where a is the lattice spacing and k is the 1D wave vector pertaining to the edge. The 1D Brillouin zone (BZ) is defined by k

∈ -π/ √ 3a, π/ √ 3a corresponding to k ∈ [-π, π].
For a given m, the bulk sites correspond to translations of the dimer with vectors -n a 1 (n positive integer). The presence of the edge breaks the translation symmetry and the Bloch theorem does not apply anymore. We label the site positions with the index n from n = 0 (edge) to n = N → ∞ (bulk) [see Fig. 3(a)].

Taking into account only nearest-neighbor couplings, each sublattice gives a tightbinding recurrence equation which defines the zero-modes [START_REF] Kohmoto | Zero modes and edge states of the honeycomb lattice[END_REF]:

i t i A i = 0 and i t i B i = 0, where i counts the three nearest-neighbors of a given B-site and Asite, respectively, t i the corresponding coupling strength, and A i and B i denoting the amplitudes at the corresponding sites. In Fig. 3(a), the missing column of sites nearest to the edge belongs to the A sublattice: A -1 ≡ 0. Thus, for zigzag edge, the tight-binding recurrence equation

βA n + (1 + e ik )A n+1 = 0 (1)
implies that all A-sites are identically null. For the B-sites, amplitudes of zero-modes fulfill the following condition:

βB n+1 + (1 + e -ik )B n = 0 . (2)
The amplitudes B n decay as

|B n | 2 = 2 β cos k 2 2n ≡ e -n/ξzz , (3)
which defines a localisation length ξ zz (k, β) plotted with a red colorscale in Fig. 3(c). The zigzag edge states exist when

|k| > 2 arccos β 2 . (4)
Conversely, for the bearded edges, sketched in Fig. 3(b), all B-sites are identically null.

For the A sites, the recurrence relation is the same as above (Eq. 1) and the nonvanishing amplitudes A n decay as

|A n | 2 = 2 β cos k 2 -2n ≡ e -n/ξ bd (5)
and

|k| < 2 arccos β 2 . ( 6
)
From Eqs ( 4) and ( 6), we obtain the existence diagram of edge states depicted in Fig. 3(c). The red and green color scales give respectively the spatial extension of zigzag and bearded edge states extracted from Eqs. ( 3) and [START_REF] Hasan | Topological insulators[END_REF]. For β = 1, the diagram shows that zigzag zero-modes occupy 1/3 of the 1D Brillouin zone whereas the proportion is 2/3 for the bearded states [START_REF] Nakada | Edge state in graphene ribbons: Nanometer size effect and edge shape dependence[END_REF][START_REF] Brey | Electronic states of graphene nanoribbons[END_REF][START_REF] Kohmoto | Zero modes and edge states of the honeycomb lattice[END_REF][START_REF] Delplace | Zak phase and the existence of edge states in graphene[END_REF]. Thus, the TB prediction for the ratio between number of bearded and zigzag states is 2. The diagram also shows that the zero-modes are more localized along zigzag edges than along bearded edges. The square dots in Fig. 2 (bottom row) represent normalized measured intensities integrated over vertical lines of resonators. Green and red zones are used as a guide for the eyes and indicate the bulk extension of the bearded and zigzag zeromodes, respectively. The total intensity in each zone is proportional to the number of corresponding zero-modes [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF]. For β = 1, the measured ratio is 2.3 in close agreement with the expected value and the bearded edge states have indeed a larger extension. For β = 0.4, we observe a majority of strongly localized bearded edge states and only a few extended zigzag zero-modes. For β = 1.5, both edge states are equally present, the zigzag ones being slightly more localized. Finally, for β = 2.5, only zigzag zero-modes exist (the intensity along the bearded edge is strictly zero) with a larger extension compared to the previous case. These observations are in good accordance with the features depicted in the diagram of Fig 3(c) and demonstrate that the anisotropy parameter allows for an accurate manipulation of edges states.

In the next section, we propose a topological argument, first presented in [START_REF] Ryu | Topological Origin of Zero-Energy Edge States in Particle-Hole Symmetric Systems[END_REF][START_REF] Delplace | Zak phase and the existence of edge states in graphene[END_REF], to address the existence of zero-modes, not only for bearded and zigzag but also for armchair edges.

Zak phase

A simple geometrical way to describe the existence of edge states is to relate their existence to a topological quantity which is a 1D winding number called the Zak phase [START_REF] Zak | Berry's phase for energy bands in solids[END_REF]. The spectrum of an infinite ribbon of finite width containing M dimers A-B consists in 2M one-dimensional bands ǫ j (k ). Therefore for each wave vector k , there are 2M states. Among these 2M states, two zero-modes may be localized along the edges. We recall below the relation between their existence and the Zak phase [START_REF] Ryu | Topological Origin of Zero-Energy Edge States in Particle-Hole Symmetric Systems[END_REF][START_REF] Delplace | Zak phase and the existence of edge states in graphene[END_REF].

Starting from the infinite unstrained system, in the (A, B) basis, the effective Hamiltonian has the general form

H k = -t 0 f * (k) f (k) 0 , (7)
where the function f (k) describes the coupling between atoms of one sublattice with the three nearest neighbors belonging to the other sublattice (see for instance Ref. [START_REF] Montambaux | A universal Hamiltonian for motion and merging of Dirac points in a two-dimensional crystal[END_REF]).

The wavefunctions have the form

ψ k (r) = 1 √ 2 ±1 e iφ(k) e ik•r , (8)
where φ(k) = arg[f (k)] and ± corresponds to positive and negative energies. The winding of the relative phase φ(k) in the reciprocal space has quite interesting properties. First of all, as seen on Fig. 4, around each Dirac point, the phase rotates by ±2π. The circulation of φ(k) along a surrounding closed path is quantized: the Berry phase is defined by 1/2 dk ∇ k φ(k) = ±π [START_REF] Neto | The electronic properties of graphene[END_REF]. Our purpose here is to stress that the winding of the φ(k) phase in reciprocal space carries an additional information related to the existence of edge states. Let us consider a ribbon geometry and define the directions parallel ( ) and perpendicular (⊥) to the ribbon length. For a ribbon of finite width containing M dimers, the perpendicular wavevector of a bulk state k ⊥ is quantized. For a two-component wavefunction, the quantization condition reads k ⊥ (M + 1)a -φ(k) = κπ with κ = 1, • • • , M [START_REF] Delplace | Zak phase and the existence of edge states in graphene[END_REF]35]. This equation has M or (M -1) roots depending on the winding of the phase under variation of k ⊥ . When only (M -1) solutions exist, the missing solution corresponds to an edge state. More precisely one has the correspondence

Z(k ) = π ←→ edge states 0 ←→ no edge states ,
where the Zak phase Z(k ) is the phase accumulated in the first 1D Brillouin zone along the k ⊥ direction:

Z(k ) = 1 2 BZ dk ⊥ ∂φ(k , k ⊥ ) ∂k ⊥ . ( 9
)
Therefore existence of edge states may be read immediately from a plot of the phase φ(k , k ⊥ ) as plotted in Fig. 4. For a given k , if the rotation of the phase is 2π, there is an edge state, if the total rotation is 0, there is no edge state. The function φ(k) depends on the boundary. For a given type of edge, the ribbon is constructed by the translations of an elementary dimer [dotted box in Figs 3(a), 3(b), and 7(a)], so that the writing of the bulk Hamiltonian has a form which depends on the considered edge [START_REF] Ryu | Topological Origin of Zero-Energy Edge States in Particle-Hole Symmetric Systems[END_REF][START_REF] Delplace | Zak phase and the existence of edge states in graphene[END_REF]. Let us consider first the zigzag edge. The function f (k) reads in this case

f zz (k) = 1 + βe ik•a 1 + e ik•a 2 = 1 + βe i( √ 3 2 k -3 2 k ⊥ )a + e i √ 3k a , (10)
where the two elementary vectors a 1 and a 2 are shown on Fig. 3(a). Similarly, the functions f (k) for the bearded [Fig. 3(b)] and armchair [Fig. 7(a)] cases read, respectively

f bd (k) = β + e ik•a 1 + e ik•a 2 = β + e i( √ 3
2 k + 3 2 k ⊥ )a + e i(- √ 3 2 k + 3 2 k ⊥ )a , (11)
f ac (k) = 1 + βe ik•a 1 + e ik•a 2 = 1 + βe i( 3 2 k - √ 3 2 k ⊥ )a + e i( 3 2 k + √ 3 2 k ⊥ )a . ( 12
)
Fig. 4 shows φ(k) for the different edges and for various deformations characterized by the parameter β. With increasing distortion, the density of edge states increases in the zigzag case, while it decreases in the bearded case as indicated by the red and green zone, respectively. In the armchair case, while there are no edge states in the undistorted lattice [START_REF] Nakada | Edge state in graphene ribbons: Nanometer size effect and edge shape dependence[END_REF][START_REF] Brey | Electronic states of graphene nanoribbons[END_REF][START_REF] Kohmoto | Zero modes and edge states of the honeycomb lattice[END_REF], they appear in the presence of a distortion which is not parallel to the edge. If the distortion is along the edge, the function f (k) reads f ac (k) = β+e ik•a 1 +e ik•a 2 instead of [START_REF] Delplace | Zak phase and the existence of edge states in graphene[END_REF]. The corresponding plot of φ(k , k ⊥ ) is shown on Fig. 5 and we see that Z(k ) = 0 everywhere, as all lines of phase jumps are parallel to k ⊥ , confirming the absence of edge states in this case. A thorough experimental investigation of armchair zero-modes is presented in the following section.

Study of armchair edges in honeycomb lattice under uniaxial strain

In previous work, we addressed the issue of topological phase transition in strained artificial graphene and incidentally observed armchair edge states [START_REF] Bellec | Topological transition of Dirac points in a microwave experiment[END_REF]. Here, we perform a quantitative study on ribbons with only armchair edges where the anisotropy axis is along their widths. Fig. 6 shows the zero-modes for β ranging from 0.4 to 3.5. It is worth noting that no zero-modes appear on armchair edges along the anisotropy axis whatever the value of β. Moreover, for the case β = 1, the ribbon does not support any armchair edge states. These features were already observed with an hexagonal flake in Ref. [START_REF] Bellec | Topological transition of Dirac points in a microwave experiment[END_REF] and are in accordance with the Zak phase analysis developed in Sec. 3. For β = 1, zero-modes along the two oblique edges are clearly present. Two main features have to be stressed: (i ) Edge states live only on one of the two triangular sub-lattices. The excited sub-lattice depends whether β is less or larger than 1, as seen in Fig. 6(a) and Fig. 6(c-f), respectively. The switch between the two sub-lattices is clearly visible in the insets. (ii ) The zero-mode localization length along the edges depends on the anisotropy and, for β > 1, decreases with β.

To understand the observations (i ) and (ii ), we propose a tight-binding analysis comparable to what has been reported in Sec. 2. Let us consider a semi-infinite honeycomb lattice with armchair edge [Fig. 7(a)]. The appropriate A -B unit cell is represented by the dashed box. Now, the armchair edge is built by translating the dimer with vectors m (a 1 + a 2 ) (m integer). Here, k ≡ 3k a, where a is the lattice spacing and k is the 1D wave vector pertaining to the armchair edge. The 1D Brillouin zone is defined by k

∈ -π 3a , π 3a corresponding to k ∈ [-π, π].
The bulk lattice is built by translating the dimer with vectors n a 2 (n positive integer). The tight-binding recurrence equations for zero-modes read [see Fig. 7(a)]:

A n+2 + A n+1 + βA n e ik = 0 , (13)
B n+1 + B n + βB n+2 e -ik = 0 . (14)
According to Eq. ( 13), for β > 1, the amplitude increases with n so that A-sites cannot support any edge states, we have A n ≡ 0, ∀n. Inversely, for β < 1, Eq. ( 14) implies B n ≡ 0. Here, we have provided a simple explanation of the switching mechanism experimentally observed and described above.

For β > 1, with B 0 = 1 as initial condition and B 1 = -e ik /β to ensure B i<0 ≡ 0, the intensities on B-sites are given by

|B n | 2 = |λr n -+ µr n + | 2 ≡ e -n/ξac , (15)
with

r ± = -1 ± 1 -4βe -ik 2βe -ik , (16)
λ = µ * = B 0 r + -B 1 r + -r - . (17)
For |r + | < |r -| < 1, we obtain decaying solutions with The diagram of existence of armchair edge states, defined by condition [START_REF] Bittner | Observation of a Dirac point in microwave experiments with a photonic crystal modeling graphene[END_REF], is plotted in Fig. 7(b). The blue color scale gives the transverse spatial extension of the armchair zero-modes extracted from Eq. ( 15). First, as pointed out in Section 3 and observed in Fig. 6(b) armchair edges along the parallel axis do not support any edge states for β = 1. Then, we see that their existence is not related to the topological transition observed in Ref. [START_REF] Bellec | Topological transition of Dirac points in a microwave experiment[END_REF]: zero-modes appear as soon as β = 1. In the case β = 0.4, A-sites support zero-modes occupying a limited range of k, as given by Eq. ( 18 18). For β = 2.5 (green diamonds) and β = 3.5 (red triangles), the agreement is very good. For β = 0.4 (blue circles) and β = 1.8 (orange squares), the difference between measured and calculated profiles is more important. To explain these discrepancies, it is worth mentioning that for β > β c zero-modes actually appear in the bandgap and are therefore clearly isolated in the DOS. When β < β c , bulk states are very close to the Dirac point and spoil the zero-mode intensity extraction. Moreover, the calculations are performed with semi-infinite lattices whereas experiments use finite ribbons. Consequently, in the experiments, we have to take into account a discrete sampling of the 1D Brillouin zone, the precise k-values depending on the ribbon length. Thus, the discrete sum of such zero-modes having different localization lengths can lead to a larger overall extension compared to the continuous case. Here again, β stands as control parameter allowing notably to switch edge states from one sublattice to the other.

|k| < arccos β 2 (3 -β 2 ) . (18)
To conclude this section, we would like to raise the existence of a new type of zero-modes in finite size systems located at the intersection of two edges, namely corner states. The comparison between Fig. 6(a) and 6(f) clearly underlines that the distribution of the zero-mode along the parallel edge is homogeneous for β = 3.5 (i.e. all edge B-sites are illuminated) whereas the intensity decreases from the top-right corner for β = 0.4. Once more, a tight-binding analysis for zero-modes allows to obtain non-homogeneous edge states. In Fig. 9(a), we build for β < 1 the amplitude of a k = 0 edge state and obtain a decaying solution. Along the oblique edge, the amplitude reads A m = mβ m-1 , m being the A-row index. Since localized edge states exist only for small k-values [dark blue zone in Fig. 7(b)], this k = 0 decaying solution dominates the overall intensity distribution. On the contrary, for β = 3.5, such decaying solutions may also exist for k = 0 but are compensated by uniform solutions associated with larger k-values, thus leading to an homogeneous overall intensity distribution.

One could expect that corner states exist for various type of edge intersections. For example, Fig. 9(b) indicates a corner state generated at the intersection of zigzag edges for β > 2. The hierarchy of its amplitudes follows the rules of a Pascal triangle. The same hierarchy can be observed in Fig. 2(d) along the direction going diagonally from the corner along a zigzag axis to the bulk. Corner states present strong similarities with zero-modes associated to localized impurities in anisotropic graphene [START_REF] Dutreix | Friedel oscillations at the Dirac cone merging point in anisotropic graphene and graphenelike materials[END_REF] and will be the subject of further studies.

Conclusion

We have studied edge states in artificial microwave graphene ribbons under uniaxial strain. By directly imaging the zero-modes, we have shown how the anisotropy parameter allows to create and manipulate edge states. Based on a nearest-neighbor tight-binding analysis and supported by topological arguments (Zak phase), we have drawn diagrams of existence and localization length variation for zigzag, bearded and armchair edge states in remarkable agreement with experiments. Our results show also that higher order nearest-neighbor couplings, inherently present in the experiments [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF], have no significant impact on the edge states. We also discuss the existence of a new type of zero-modes at the intersection of two type of edges, namely the corner states whose study is left for a future work.

Figure 1 .

 1 Figure 1. (a) Picture of an unstrained artificial graphene ribbon with zigzag, bearded and armchair edges. The lattice constant is 15 mm. (b) Corresponding experimental density of states (DOS). The arrow indicates the zero-modes appearing at the Dirac frequency ν D .

Figure 2 .

 2 Figure 2. (top row). Experimental zero-mode intensities for β ranging from 0.4 to 2.5. The anisotropy axis is horizontal. Bearded and zigzag edges are displayed on the same lattice (respectively at the left and the right). (bottom row). Normalized intensities integrated over vertical lines. Green (resp. red) areas indicate qualitatively the bearded (resp. zigzag) zero-mode extension.

Figure 3 .

 3 Figure 3. Schematic representation of a semi-infinite lattices with zigzag [(a), red] and bearded [(b), green] edge along the axis. The anisotropy axis is horizontal, here along the ⊥ axis. The corresponding coupling strength is denoted t ′ . A n and B n are the wavefunction amplitudes in the two sublattices. k is the 1D Bloch wavevector (see text for details). (c) Zigzag -bearded edge states existence diagram for various k and β. The colormap corresponds to the inverse localization length in the zigzag (red, 1/ξ zz ) and bearded (green, 1/ξ bd ) cases, the darker the color, the stronger the states are localized. The dashed lines correspond to the values of β used in the experiments.

Figure 4 .

 4 Figure 4. Stream and density plots of the phase φ(k , k ⊥ ) in reciprocal space, for the three edges (zigzag, bearded and armchai) and for various deformations (see the sketches in Figs 3 and 7). The first 2D (bulk) BZ is indicated by a white contour. The red dots are the positions of the Dirac points. The first 1D BZ associated to the ( ) and (⊥) directions are shown by solid and dashed black lines, respectively. For a given k , the existence of edge states is related to the winding of the phase along the k ⊥ direction in the BZ: there are edge states when the total rotation of the phase is 2π as seen in the shaded areas [Z(k ) = π]. There are no edge states otherwise [Z(k ) = 0].

  armchair

Figure 5 .

 5 Figure 5. Stream and density plots of the phase φ(k , k ⊥ ) in reciprocal space, for the armchair edge and a deformation parallel to the edge. Marker and line legends are the same as in Fig. 4. The Zak phase is 0 everywhere, confirming the absence of edge state in this case.

Figure 6 .

 6 Figure 6. Experimental armchair zero-mode intensities for β ranging from 0.4 to 3.5. (a) β < 1. (b) β = 1, no edge states. (c -f) β > 1. Insets in (a) and (f) emphasize the switch between the two triangular sublattice excitation.

Figure 7 .

 7 Figure 7. (a) Schematic representation of a semi-infinite lattice with armchair edge along the axis. The anisotropy axis is horizontal, the corresponding coupling strength is denoted t ′ . A n and B n are the wavefunction amplitudes in the two sublattices. k is the Bloch wavevector (see text for details). (b) Armchair edge states existence diagram for various k and β. The colormap corresponds to the inverse localization length 1/ξ ac . Light (resp. dark) colors means delocalized (resp. localized) state. The white area means no edge states. For β > 1, only B-sites are excited. For β < 1, only A-sites are excited.

Figure 8 .

 8 Figure 8. Experimental (a) and theoretical (b) armchair zero-mode intensity profiles for β ranging from 0.4 to 3.5. For β < 1 (resp. β > 1), A (resp. B) sites are considered.

  ), almost half of them being well localized (dark blue colors in the diagram). For 1 < β < β c = 2, the edge states supported by the B-sites are still belonging to a finite range of k and are mainly delocalized (light blue colors in the diagram), as observed in the Fig. 6(b) and (c). For β > 2, zero-modes run over the full 1D Brillouin zone and their localization lengths decrease with β. To be more quantitative, we plotted in Fig. 8 the intensities of B n (resp. A n ), for β > 1 (resp. β < 1) and n ranging from 0 to 5. Fig. 8(a) shows the experimental profiles extracted from Fig. 6 by integrating over vertical lines of sites. Fig. 8(b) depicts the intensities calculated from Eqs. (15 -17) with an integration of k over all possible values defined by Eq. (

1 Figure 9 .

 19 Figure 9. (a) Armchair zero-mode amplitude components for β < 1 and k = 0. The wavefunction is non-zero only on A-sites (empty circles). m indicates the A-row index. The anisotropy axis is vertical, thick lines represent the corresponding coupling. (b) Zero-mode amplitude components for β > 2 and k = 0 appearing at the intersection of two zigzag edges.

Acknowledgments

The authors are thankful to Pierre Delplace for enlightening discussions concerning the Zak phase during the international workshop on "Artificial Graphene" held at the Centro Internacional de Ciencias A.C (Cuernavaca, Mexico) in 2013. This work was supported by the DFG via the Forschergruppe 760.