Sophie Tarbouriech

Isabelle Queinnec
email: queinnec@laas.fr

Jean-Marc M Biannic
email: jean-marc.biannic@onera.fr

Christophe Prieur
email: christophe.prieur@gipsa-lab.fr

S Tarbouriech -I

Queinnec

Pilot-Induced-Oscillations alleviation through anti-windup based approach *

The chapter is dedicated to the optimization of a well-known structure of compensators: the anti-windup scheme. This approach belongs to the saturation allowance control class which aims to exploit at the most the actuators capabilities. The objective of this chapter consists of adapting and developing the anti-windup compensator design to some particular classes of nonlinear actuators presenting both magnitude and rate saturations. It is illustrated on the lateral flying case for a civil aircraft in presence of aggressive maneuvering of the pilot. A complete methodology is then proposed comparing several approaches including given anti-PIO filters.

Introduction

Control engineers, where possible, like to work under the assumption of linearity. The mathematics associated with the field of linear systems is well developed and underpins much of the control theory which is applied in industry. Even nonlinear techniques often attempt to generalize linear concepts, and frequently nonlinear systems are linearized to obtain linear models which locally yield good engineering approximations [START_REF] Khalil | Nonlinear systems[END_REF]. The problem with the assumption of linearity is that it is sometimes unrealistic and can lead to erroneous results. Actually, the increasing re-quirements in terms of operational reliability and performance ask to work beyond the linear behavior of the system. Hence, actuators saturations (both magnitude and rate saturations) represent a common nonlinear phenomenon in almost all physical applications, especially in space and aeronautical fields. There are many examples of saturation problems but perhaps the most notorious are those associated with so-called pilot-induced-oscillations (PIO's) in aeronautics (see, for example, [START_REF] Klyde | Pilot-induced oscillation analysis and prediction with actuator rate limiting[END_REF], [START_REF] Duda | Prediction of pilot-in-the-loop osciallations due to rate saturation[END_REF], [START_REF] Klyde | Use of active inceptor cueing to mitigate pilot-vehicle system loss of control[END_REF]). These saturation-induced events have led to the crash of several aircrafts (the SAAB Grippen and the Boeing V22 Osprey are notorious examples [START_REF] Anon | Why the gripen crashed[END_REF]) and several near-misses with others. Actually, recall that a pilot-induced oscillation is a sustained or uncontrollable, undesired oscillation resulting from the action of the pilot to control the aircraft. A common nonlinearity leading to PIO is control surface rate limiting. Then this phenomenon can introduce a delayed response and then the action of the pilot implies that the airplane response is essentially opposite of the command wished by the pilot (see, in particular the recent chapter [START_REF] D'andrea | Challenges Actuated Wingsuit for Controlled, Self-Propelled Flight, volume The Impact of Control Technology[END_REF]) Thus the presence of saturation can lead to performance degradation from the mild to the severe and can also lead to loss of stability [START_REF] Hu | Control systems with actuator stauration: Analysis and design[END_REF], [START_REF]Actuator saturation control[END_REF]. Although this is not always critical, it is clear that some way of predicting the effects of saturation is required and, moreover, that some method of limiting the degradation that occurs is warranted. This reflects the need for the development of new and more complex control techniques in order to meet the new demands.

In the aeronautical literature, there exist some methods mainly based on the addition of filters or estimators designed to predict and reduce the risk of PIO (see, for example the OLOP criterion using describing functions studied in [START_REF] Gilbreath | Prediction of PIO due to actuator rate limiting using the open-loop onset point (OLOP) criterion[END_REF] or the use of a detector based on short time Fourier transform and autoregressive model [START_REF] Liu | Pilot-induced-oscillation detection and mitigation[END_REF]. In this chapter, we choose another route by considering the use of anti-windup technique,with the objective to provide constructive conditions (that is associating theoretical conditions to optimization routines in order to exhibit effective numerical solutions). More specifically, the approach proposed in this chapter is based on the optimization of a well-known structure of compensators: the anti-windup scheme (see e.g., [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF][START_REF] Zaccarian | Modern anti-windup synthesis[END_REF] for an introduction of this notion). This approach belongs to the saturation allowance control class which aims to exploit at the most the actuators capabilities. The basic concept consists of introducing an extra layer to the existing linear controller, accounting for the nonlinearities in order to mitigate the windup phenomenon created by the saturation [START_REF] Hippe | Windup in control. Its effects and their prevention[END_REF]. This strategy, also called anti-windup design, allows the designer to keep the existing linear controller (already validated) and to introduce a compensator which is active only when the nonlinearity arises. In this framework, numerous works have emerged in the context of both magnitude and rate saturation constraints [START_REF] Rundquist | Phase compensation of rate-limiters in unstable aircraft[END_REF], [START_REF] Brieger | Anti-windup compensation of rate saturation in an experimental aircraft[END_REF], [START_REF] Galeani | A magnitude and rate saturation model and its use in the solution of a static anti-windup problem[END_REF] [13], [START_REF] Galeani | A tutorial on modern anti-windup design[END_REF], [START_REF] Tarbouriech | Anti-windup design: an overview of some recent advances and open problems[END_REF]. Such an approach appears to be really attractive as the anti-windup loop may work with existing control laws (a priori designed by the engineers to answer to defined requirements). Indeed, it represents an interesting technique for the controller designers who can use familiar and intuitive techniques for them and then, simply add an extra layer, which will consider the nonlinear behavior in a second step. If originally, results on antiwindup design consisted on ad-hoc methods intended to work with PID controllers [START_REF] Fertik | Direct digital control algorithm with anti-windup feature[END_REF], [START_REF] Åström | Integrator windup and how to avoid it[END_REF], modern anti-windup methods have emerged during the last decade (see, for example, [START_REF] Teel | Anti-windup for exponentially unstable linear systems[END_REF], [START_REF] Grimm | Antiwindup for stable linear systems with input saturation: an LMI based synthesis[END_REF]). Then, the design of such an additional compensator is generally carried out through a static optimization problem of the controller parameters. Thanks to the development of semi-definite programming and convex optimization [START_REF] Boyd | Convex optimization[END_REF], the anti-windup controller design problem can be formulated as the optimization of a multi-objective criterium (corresponding to closed-loop stability and performance specifications) subject to matrix inequalities constraints associated to the dynamical system. Many different techniques exist in control theory to synthesize such anti-windup controllers, among which static and dynamic linear anti-windup augmentation (see [START_REF] Zaccarian | Modern anti-windup synthesis[END_REF]) based on a generalized sector condition representing the saturation [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]. Other anti-windup augmentations are possible as nonlinear synthesis, in particular for control systems equipped with other nonlinearities than magnitude and rate saturations. One can consult [START_REF] Andrieu | Global asymptotic stabilization of systems satisfying two different sector conditions[END_REF] for control systems presenting two different sector conditions, and [START_REF] Tarbouriech | Stability analysis and stabilization of systems with input backlash[END_REF] for a control system with a memory-based input.

The objective of this chapter is to adapt and develop the anti-windup compensator design to a class of systems presenting both magnitude and rate saturations. The techniques proposed first include a modeling of the nonlinear actuator involved to further derive analysis and design conditions. It is illustrated on the lateral flying case for a civil aircraft in presence of aggressive maneuvering of the pilot. A complete methodology is then proposed comparing several approaches including given anti-PIO filters (borrowed mainly from [START_REF] Biannic | Analyse et ajustement de lois de commande en présence de saturations implantation de filtres anti-PIO générés par synthèse anti-windup[END_REF] and [START_REF] Puyou | Application of robust antiwindup design to the longitudinal aircraft control to cover actuator loss[END_REF]).

The chapter is organized as follows. In Section 2, a complete model of plant, actuator and controller involved to address the stability and performance optimization problem is described. Then, the multi-objective problem to be solved in order to design anti-windup loops is stated. Section 3 pertains to the anti-windup design conditions in two cases depending on the signal used as the input of the anti-windup controller. Then, in order to alleviate the PIO risk for a civil aircraft in presence of aggressive maneuvering of the pilot, Section 4 depicts how the previous techniques are very interesting in comparison with classical anti-PIO filters to guarantee stability and performance of the closed-loop system. Several simulations illustrate the benefits provided by the anti-windup compensators, in terms of simple and systematic methods without needing a tuning parameters step. Finally, some concluding remarks end the chapter.

Model description and problem formulation

Anti-windup strategies represent an appropriate framework to mitigate the undesired saturation effects [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF], [START_REF] Zaccarian | Modern anti-windup synthesis[END_REF]. Thus, the general principle of the anti-windup scheme can be depicted in Figure 1, where the (unconstrained) signal produced by the controller is compared to that which is actually fed into the plant (the constrained signal). This difference is then used to adjust the control strategy by preserving stability and performance. The kind of anti-windup controller used in the chapter is specified later and is strongly depending on the considered plant, actuator and controller. Let us first describe the complete model.

Plant model

Unlike most systems in the literature, the outputs of the controller are not affected in a same way by the nonlinear elements. Then, the vector u ∈ ℜ m building the m inputs of the plant is decomposed into two subvectors: the first one, denoted u s ∈ ℜ m s , corresponds to m s saturated inputs, whereas the second one, denoted u ns ∈ ℜ m-m s , corresponds to the linear inputs (unsaturated inputs). The plant model can be defined by: sysP :

       ẋp = A p x p +
where x p ∈ ℜ n p and y p ∈ ℜ p are the state and the measured output of the plant. w ∈ ℜ q generally represents an exogenous perturbation but may also be used to represent a reference signal (or both). Furthermore, z ∈ ℜ l represents the regulated output, which is used to evaluate the performance of the system with respect to the perturbation w via some pertinent optimization criteria.

Controller model

Differently from the classical anti-windup loops, in which the output of the antiwindup controller is injected to the dynamics of the controller and/or the output of the controller, we consider here that the output of the anti-windup controller modifies only partially the dynamics of the controller and/or the output of the controller. Then, with this in mind, the dynamical controller is described as follows:

sysC :

       ẋc = A c x c + B c u c + B cw w + B ca v x y cs = C s c x c + D s c u c + D s cw w + D ca v y y cns = C ns c x c + D ns c u c + D ns cw w (2)
where x c ∈ ℜ n c and u c ∈ ℜ p are the state and the input of the controller. The output of the controller is decomposed into two signals: y cs ∈ ℜ m s , which will be interconnected to u s through a saturated actuator, and y cns ∈ ℜ m-m s , which will be interconnected with the linear (unsaturated) input u ns . Moreover, v x and v y are the additional inputs that will be connected to the anti-windup controller. B ca and D ca are matrices of dimensions n c ×n cr and m s ×m r , and allow to specify what are the n cr states and m r outputs modified by the anti-windup action.

Actuator model

There is an actuator block between the output of the controller y c and the input of the plant u, which is decomposed into two blocks: the first one corresponding to the nonlinear (saturated) part and the second one corresponding to the linear (unsaturated) part. The nonlinear actuator part involves n dz nested saturations, including the case of rate and magnitude saturations, as depicted in Figure 2(a). Such nonlinearities will will be tackled via the use of dead-zone, denoted φ i (.), i = 1...n dz . The dynamical model of the actuator is based on the scheme 2(b) as follows:

sysACT :

       ẋa = v + φ 1 (v) v = T 0 y cs + T 0 φ 0 (y cs ) -T 0 x a u s = x a (3)
with φ 0 (y cs ) = sat u 0 (y cs )y cs and φ 1 (v) = sat u 1 (v)v, where sat u 0 (.) and sat u 1 (.) are classical saturation functions and u 0 and u 1 are the levels of saturation in magnitude and in rate, respectively. The elements of the diagonal matrix T 0 ∈ ℜ m s ×m s classically take values large enough in order to avoid affecting the linear dynamics of the closed-loop system.

Interconnections

The interconnections considered in the chapter can be described as follows:

• linear link between the output of the plant and the input of the controller: u c = y p ;

• the first part of the output of the controller (y cs ) is linked to the corresponding inputs of the plant (u s ) through the actuator model ( 3); • the second part of the output of the controller is directly connected to the corresponding inputs of the plant: u ns = y cns ; • v x and v y are built from the anti-windup compensator (and will be specified later).

Anti-windup compensator

In the DLAW (Direct Linear Anti-Windup) strategy, the anti-windup controller uses as input the difference between the signals issued either from the input and the output of the whole actuator or from the input and the output of the nonlinear elements included in the actuator. Following this, we pursue two strategies to design the antiwindup loops.

• The first strategy is reported in [START_REF] Biannic | Analyse et ajustement de lois de commande en présence de saturations implantation de filtres anti-PIO générés par synthèse anti-windup[END_REF] and considers the difference between the input and the output of the actuator defined by e = u sy cs ∈ ℜ m s . Additionally, one assumes that the anti-windup controller only acts on the dynamics of the controller, which corresponds to v y = 0, or equivalently, m r = 0. The anti-windup controller of order n aw , with v x ∈ ℜ n cr , reads:

AW e : ẋaw = A aw x aw + B e aw (u s -y cs ) v x = C aw x aw + D e aw (u s -y cs ) (4)
• The second strategy considers that the input of the anti-windup controller are the dead-zones associated to each saturation. Hence, the anti-windup controller of order n aw reads:

AW φ :      ẋaw = A aw x aw + B 0 aw φ 0 (y c ) + B 1 aw φ 1 (v) v x v y = C aw x aw + D 0 aw φ 0 (y c ) + D 1 aw φ 1 (v) (5)
where v x and v y are of dimensions n cr and m r , respectively.

Remark 1. The interest of the second anti-windup structure resides in the simplicity of the design conditions. Indeed, the design of a static anti-windup gain (only matrices D 0 aw and D 1 aw are used) is issued from a fully linear problem. In the case of the design of a dynamical anti-windup controller, for a priori given matrices A aw and C aw , the determination of input and transmission matrices is also obtained by solving a linear problem. In the case where n aw = n p + n m s + n c , the resolution of a linear problem can also be considered [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]. At the opposite, the first strategy is more adapted to provide analysis conditions but does not allow to simultaneously compute the matrices of the anti-windup and the matrix of the Lyapunov function through a linear optimization problem, even in the static anti-windup case.

Remark 2. The anti-windup model ( 5) imposes the assumption that the input and output signals of each saturation block is available. To overcome this assumption, alternative strategies can be investigated. For example, the anti-windup may use the difference between the nonlinear actuator and a linear fictitious one (with the same dynamics but without saturation blocks) to explicitly take into account the dynamics of the actuator (present in the rate limiter) [START_REF] Queinnec | Anti-windup design for aircraft control[END_REF]. Another option would be to build an observer to evaluate the internal state of the actuator [START_REF] Tarbouriech | Anti-windup design with rate and magnitude actuator and sensor saturations[END_REF]. In these cases, conditions can be derived in a simpler way than that ones issued from the strategy with ( 4), but they remain more complex than those due to the strategy with [START_REF] Boyd | Convex optimization[END_REF].

Standard formulation

In [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF], a standard formulation of the anti-windup design has been proposed for different kinds of actuators. In the current case, by considering an augmented state of dimensions n = n p + m s + n c + n aw including the state of the plant, the state of the actuator, the state of the controller and the state of the anti-windup controller, the following standard model of the complete closed-loop system can be defined by:

             ẋ = A x + B 0 φ 0 (y c ) + B 1 φ 1 (v) + B 2 w y c = C 0 x + D 00 φ 0 (y c ) + D 01 φ 1 (v) + D 0w w v = C 1 x + D 10 φ 0 (y c ) + D 11 φ 1 (v) + D 1w w z = C 2 x + D 20 φ 0 (y c ) + D 21 φ 1 (v) + D 2w w (6)
Then, depending of the anti-windup scheme under consideration, the matrices of the anti-windup controller are encapsulated into the matrices of system [START_REF] Brieger | Anti-windup compensation of rate saturation in an experimental aircraft[END_REF].

The design procedure of the anti-windup controller consists in optimizing some quantities as the size of the region of stability of the closed-loop system or the guaranteed level of performance. Several optimization problems are then of interest. In particular, the idea by adding the anti-windup loop is to maximize the basin of attraction of the origin for the closed-loop system and/or to minimize the L 2 gain between w and z or to maximize the set of perturbation w, which can be rejected. Then, throughout the chapter, the signal of perturbation is supposed to be bounded in energy as follows:

w 2 2 = ∞ 0 w (t)w(t)dt ≤ δ -1 ; 0 ≤ δ -1 < ∞ (7)
The problem we intend to address in the chapter can be summarized below.

Problem 1. Determine an anti-windup controller and a region E , as large as possible, such that

• Internal stability. The closed-loop system (6) with w = 0 is asymptotically stable for any initial conditions belonging to E (which is a region of asymptotic stability (RAS)); • Performance. The L 2 gain between w and z is finite and equal to γ > 0.

The convex optimization problems associated to Problem 1 are specified in Sections 3.2 and 3.3.

Main Anti-Windup Design Conditions

Solution to standard anti-windup design

The following proposition provides conditions of local stability and L 2 performance for the closed-loop system [START_REF] Brieger | Anti-windup compensation of rate saturation in an experimental aircraft[END_REF]. The result regards existence conditions to solve Problem 1.

Proposition 1. If there exist a symmetric positive definite matrix Q ∈ ℜ n×n , two matrices Z 0 and Z 1 ∈ ℜ m×n , two positive diagonal matrices S 0 and S 1 ∈ ℜ m×m and a positive scalar γ such that the following conditions are verified:

      QA + A Q B 0 S 0 -QC 0 -Z 0 B 1 S 1 -QC 1 -Z 1 B 2 QC 2 -2S 0 -D 00 S 0 -S 0 D 00 -D 01 S 1 -S 0 D 10 -D 0w S 0 D 20 -2S 1 -D 11 S 1 -S 1 D 11 -D 1w S 1 D 21 -I D 2w -γI       < 0 (8) Q Z 0(i) δ u 2 0(i) ≥ 0, i = 1, ..., m (9)
Q Z 1(i) δ u 2 1(i) ≥ 0, i = 1, ..., m (10)
then, 1. when w = 0, the set

E (Q -1 , δ ) = {x ∈ ℜ n ; x Q -1 x ≤ δ -1
} is RAS for the closedloop system (6); 2. when w = 0, satisfying [START_REF] D'andrea | Challenges Actuated Wingsuit for Controlled, Self-Propelled Flight, volume The Impact of Control Technology[END_REF], and for x(0) = 0,

• the trajectories of the closed-loop system remain bounded in the set E (Q -1 , δ ); • the L 2 gain is finite and one gets:

T 0 z(t) z(t) dt ≤ γ T 0 w(t) w(t) dt, ∀T ≥ 0 (11)
In an analysis purpose (the anti-windup controller being given), conditions of Proposition 1 are linear and can be directly used to solve adequate optimization problems. Moreover, in the design context, conditions of Proposition 1 are non convex, matrices A aw , B aw , C aw and D aw being hidden in matrices A , B i , C i , D i j , i, j = 0, 1. Depending on the problem studied, conditions linear in the decision variables can be obtained, more or less directly, by modifying a bit the original conditions or still by considering iterative procedures (including D-K iteration process) allowing to search for Lyapunov matrix and anti-windup matrices. These situations are detailed in the next subsections.

Remark 3. In the sequel, one considers a set X 0 , defined by some directions in the plant state space v i ∈ ℜ n p , i = 1, ..., q, to provide a desired shape of the region E (Q -1 , δ ) to be maximized when solving Problem 1. Then, considering vi = v i 0 ∈ ℜ n , i = 1, ..., q and β a scaling factor such that β X 0 ⊂ E (Q -1 , δ ), an additional condition to those of Proposition 1 have to be considered in the algorithms which follow:

δ 1 β 2 δ v i δ vi Q > 0, i = 1, ..., q (12)

Algorithms for AW e case

From (1), ( 2), ( 3) and ( 4), matrices of system (6) read:

A = A 0 0 0 + B νB A aw C νA + B νB B e aw C νC + B νD C aw C νA + B νD D e aw C νC B 0 = B φ 0 0 ; B 1 = B φ 1 0 ; B 2 = B 2 0 + B νB B e aw C νW + B νD D e aw D 0w C 0 = C 0 0 ; C 1 = C 1 0 ; C 2 = C 2 0 D 00 = 0 ; D 01 = 0 ; D 10 = D 1 ; D 11 = 0 D 20 = 0 ; D 21 = 0 ; D 0w = D 0w ; D 1w = D 1w ; D 2w = D 2w (13)
with

A=       A p + B ns pu ∆ -1 D ns c C p B s pu + B ns pu ∆ -1 D ns c D s pu B ns pu ∆ -1 C ns c T 0 D s c (C p +D ns pu ∆ -1 D ns c C p ) T 0 (D s c D s pu -I ms +D s c D ns pu ∆ -1 D ns c D s pu ) T 0 (C s c +D s c D ns pu ∆ -1 C ns c ) B c C p +B c D ns pu ∆ -1 D ns c C p B c D s pu +B c D ns pu ∆ -1 D ns c D s pu A c + B c D ns pu ∆ -1 C ns c       B 2 =     B pw + B ns pu ∆ -1 (D ns c D pw + D ns cw ) T 0 (D s c D pw + D s cw + D s c D ns pu ∆ -1 (D ns c D pw + D ns cw )) B cw + B c D pw + B c D ns pu ∆ -1 (D ns c D pw + D ns cw )     B φ 0 =   0 T 0 0   ; B φ 1 =   0 I m s 0   ; D 1 = T 0 C 0 = D s c (I p + D ns pu ∆ -1 D ns c )C p D s c (I p + D ns pu ∆ -1 D ns c )D s pu C s c + D s c D ns pu ∆ -1 C ns c C 1 = T 0 D s c (I p + D ns pu ∆ -1 D ns c )C p T 0 D s c (I p + D ns pu ∆ -1 D ns c )D s pu -T 0 T 0 (C s c + D s c D ns pu ∆ -1 C ns c ) C 2 = C z + D ns zu ∆ -1 D ns c C p D s zu + D ns zu ∆ -1 D ns c D s pu D ns zu ∆ -1 C ns c D 0w = D s cw + D s c D pw + D s c D ns pu ∆ -1 (D ns c D pw + D ns cw ) D 1w = T 0 (D s cw + D s c D pw + D s c D ns pu ∆ -1 (D ns c D pw + D ns cw )) D 2w = D zw + D ns zu ∆ -1 (D ns c D pw + D ns cw ) B νB =     0 0 0 I n aw     ; B νD =     0 0 B ca 0     ; C νA = 0 0 0 I n aw C νC = -D s c (I p + D ns pu ∆ -1 D ns c )C p I m -D s c (I p + D ns pu ∆ -1 D ns c )D s pu -C s c -D s c D ns pu ∆ -1 C ns c 0 and ∆ = I m-ms -D ns c D ns pu .
In the analysis context, conditions of Proposition 1 using the AW e structure are linear in the decision variables and can be directly used. On the other hand, in the design context, conditions of Proposition 1 are nonlinear due to, in particular, the products between the Lyapunov matrix Q and the matrices of the anti-windup controller. Then, to address the design and solve Problem 1, some iterative procedure can be applied by considering at the first step a given static (n aw = 0) or dynamic (n aw = 0) anti-windup controller.

The following algorithms can be used.

Algorithm 3.1 Analysis of a given AW e anti-windup controller 1. Give matrices A aw , B e aw , C aw and D e aw . 2. Choose directions to be optimized v i ∈ ℜ n p , i = 1, ..., q and a known perturbation with bound δ .

Solve min

Q,S 0 ,S 1 ,Z 0 ,Z 1 ,γ,µ γ + µ subject to LMI (8), ( 9), ( 10) and [START_REF] Gilbreath | Prediction of PIO due to actuator rate limiting using the open-loop onset point (OLOP) criterion[END_REF] where γ is the L 2 gain between w and z and µ = 1/β 2 .

Algorithm 3.2 Design of a AW e anti-windup controller 1. Select an initial guess for matrices A aw , B e aw , C aw and D e aw . of appropriate dimensions in order to build the desired anti-windup loop. A static anti-windup AW e may also be used by considering n aw = 0. 2. Choose directions to be optimized v i ∈ ℜ n p , i = 1, ..., q and a known perturbation with bound δ . 8), ( 9), ( 10) and [START_REF] Gilbreath | Prediction of PIO due to actuator rate limiting using the open-loop onset point (OLOP) criterion[END_REF] where γ is the L 2 gain between w and z and µ = 1/β 2 . 4. If the solution obtained is satisfactory (some accuracy has to be fixed) or no more improved from the previous steps then STOP. Otherwise, go to the next iteration (the idea is to finish by an analysis step). 5. Synthesis step -Pick the solution Q obtained at Step 3 and solve min A aw ,B e aw ,C aw ,D e aw ,S 0 ,S 1 ,Z 0 ,Z 1 ,γ γ subject to LMI (8), ( 9), ( 10) and (12)

Analysis step

-Solve min Q,S 0 ,S 1 ,Z 0 ,Z 1 ,γ,µ γ + µ subject to LMI (
6. Go to Step 3.

Remark 4. The selection of an initial guess of anti-windup in the Algorithm 3.2 must take care of the dimension of each elements but must also verify that A aw is Hurwitz. Actually, it is not possible to initialize the problem with null matrices of appropriate dimensions (for a given order of the anti-windup scheme n aw ) as the condition on A in the first block of inequality [START_REF] Duda | Prediction of pilot-in-the-loop osciallations due to rate saturation[END_REF] imposes that both the closed-loop linear dynamics of the system and the anti-windup dynamics are asymptotically stable. An option may then be to select any stable dynamical matrix A aw with matrices B e aw , C aw and D e aw equal to null matrices of appropriate dimensions. This initial anti-windup scheme is ineffective but allows to solve the analysis steep and obtain a matrix Q to be used in the synthesis step.

Algorithms for AW φ case

As for the previous case, from (1), ( 2), ( 3) and ( 5), matrices of system (6) are defined by:

A = A B v C aw 0 A aw ; B 0 = B φ 0 + B v D 0 aw B 0 aw ; B 1 = B φ 1 + B v D 1 aw B 1 aw C 0 = C 0 C v0 C aw ; C 1 = C 1 C v1 C aw ; C 2 = C 2 0 D 00 = C v0 D 0 aw ; D 01 = C v0 D 1 aw ; D 10 = D 1 +C v1 D 0 aw ; D 11 = C v1 D 1 aw B 2 = B 2 0 ; D 20 = 0 ; D 21 = 0 ; (14) Matrices A, B 2 , B φ 0 , B φ 1 , D 1 , C 0 , C 1 , C 2 , D 0w
, D 1w and D 2w remain unchanged. Matrices defining the interconnection between the anti-windup loop and the system are:

B v =    0 T 0 D ca 0 I m r B ca I n cr 0    ; C v0 = D ca 0 I m r ; C v1 = T 0 D ca 0 I m r
As in the previous case, the analysis problem (Algorithm 3.3) is linear and the synthesis problem of the anti-windup is nonlinear, including products between decision variables, and in particular between the Lyapunov matrix Q and the anti-windup elements. As for the AW e strategy, a D-K iteration procedure may then be considered for the synthesis problem (Algorithm 3.4). However, differently for the AW e strategy, the synthesis optimization problem may be partialy linearized and, for given matrices A aw and C aw , the design of matrices B i aw and D i aw , i = 0, 1 can be handled via a linear optimization problem (Algorithm 3.5). Algorithm 3.3 Analysis of a given AW φ anti-windup controller 1. Select matrices A aw , B 0 aw , B 1 aw , C aw , D 0 aw and D 1 aw . 2. Choose directions to be optimized v i ∈ ℜ n p , i = 1, ..., q and a known perturbation with the bound δ . 8), ( 9), ( 10) and [START_REF] Gilbreath | Prediction of PIO due to actuator rate limiting using the open-loop onset point (OLOP) criterion[END_REF] where γ is the L 2 gain between w and z and µ = 1/β 2 .

Solve min

Q,S 0 ,S 1 ,Z 0 ,Z 1 ,γ,µ γ + µ subject to LMI (
Algorithm 3.4 Design of a AW φ anti-windup controller 1. Select matrices A aw , B 0 aw , B 1 aw C aw , D 0 aw and D 1 aw of appropriate dimensions in order to build the desired anti-windup loop. A static anti-windup AW φ may also be used by considering n aw = 0. 2. Choose directions to be optimized v i ∈ ℜ n p , i = 1, ..., q and a known perturbation with bound δ . 3. Analysis step -Solve min Q,S 0 ,S 1 ,Z 0 ,Z 1 ,γ,µ γ + µ subject to LMI ( 8), ( 9), ( 10) and [START_REF] Gilbreath | Prediction of PIO due to actuator rate limiting using the open-loop onset point (OLOP) criterion[END_REF] where γ is the L 2 gain between w and z and µ = 1/β 2 . 4. If the solution obtained is satisfactory (some accuracy has to be fixed) or no more improved from the previous steps then STOP. Otherwise, go to the next iteration (the idea is to finish by an analysis step). 5. Synthesis step -Pick the solution Q obtained at Step 3 and solve

min S 0 ,S 1 ,Z 0 ,Z 1 ,B 0 aw ,B 1 aw ,D 0 aw ,D 1 aw ,γ
γ subject to LMI ( 8), ( 9), ( 10) and ( 12)

6. Go to Step 3.

Algorithm 3.5 Design of a AW φ anti-windup controller with fixed dynamics 1. Give matrices A aw and C aw .

2. Choose directions to be optimized v i ∈ ℜ n p , i = 1, ..., q and a known perturbation with the bound δ .

Solve min

Q,S 0 ,S 1 ,Z 0 ,Z 1 , B0 aw , B1 aw , D0 aw , D1 aw ,γ,µ γ + µ
subject to LMI (8), ( 9), ( 10) and [START_REF] Gilbreath | Prediction of PIO due to actuator rate limiting using the open-loop onset point (OLOP) criterion[END_REF] where γ is the L 2 gain between w and z and µ

= 1/β 2 . 4. Compute B 0 aw = B0 aw S -1 0 , B 1 aw = B1 aw S -1 1 , D 0 aw = D0 aw S -1 0 and D 1 aw = D1 aw S -1 1 .
Remark 5. In Algorithm 3.5, condition [START_REF] Duda | Prediction of pilot-in-the-loop osciallations due to rate saturation[END_REF] is not directly applied. The products between B i aw and D i aw with the matrices S i are replaced by the change of variables Bi aw and Di aw , i = 0, 1, which allows to linearize the problem.

Remark 6. An interesting case is the static anti-windup one, for which matrices A aw and C aw are null matrices of appropriate dimensions. It implies that B i aw , i = 0, 1 are also null matrices of appropriate dimensions and only matrices D i aw , i = 0, 1 are computed in the linear Algorithm 3.5.

Remark 7.

Matrices A aw and C aw to be used in Algorithm 3.5 may be selected as the solution of a full-order (n aw = n p + n c + m s ) anti-windup compensator design where the actuator is just a saturation in magnitude (see, for example, the conditions provided in [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]), i.e. via a linear optimization problem. Eventually, an order-reduction step may also be considered in order to select matrices A aw and C aw (see Example 8.5 in [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]). Other procedures developed in [START_REF] Zaccarian | Modern anti-windup synthesis[END_REF] could be used.

Anti-windup and its use for PIO alleviation

The design and analysis algorithms of Section 3 are now applied and compared in the realistic context of lateral maneuvers of a civil transport aircraft. A specific attention will be devoted to aggressive pilot's demands in conjunction with actuator loss.

In the following, pilot's activity is modeled as a static gain K pil . For this application, a normal activity would correspond to K pil = 1. But, in stressful situations, notably in case of actuators loss, a more aggressive pilot's behavior is generally observed, resulting in much higher gains. Here, the gain is set to K pil = 3.

Problem setup and objectives

The two anti-windup structures AW e and AW φ above described are compared in this applicative part of the chapter. Both nonlinear closed-loop Simulink implementations are sketched in Figures 3 and4, respectively. For each design strategy, the state-space models sysP and sysC are readily obtained from the Simulink diagrams of Figures 3 and4 with the help of the Matlab linmod function. The plant corresponds to the "yellow box" depicting the aircraft system while the global controller (including pilot actions) is obtained by extraction of the 3 blue boxes. A standard balanced reduction technique is finally applied to obtain reasonably sized models. The obtained reduced orders, respectively n p = 8 and n c = 20, are compatible with the proposed algorithms. The aircraft system involves 2 inputs (m = 2) with only the ailerons deflection actuator which saturates (m s = 1), and 5 outputs (p = 5), among which the performance output which is set as the roll angle (l = 1). The disturbance of the system is due to the saturation of the input, i.e., B pw = B s pu (q = 1). In the AW e strategy, the anti-windup input is a single scalar signal (m s = 1) which only captures the difference between the input of the nonlinear ailerons deflection actuator and its output. In the AW φ strategy, two signals (one for the magnitude limitation and one for the rate limitation) are used by the anti-windup device. Their generation is detailed in the Simulink implementation of Figure 5.

Whatever the considered approach, both anti-windup controllers act similarly on the internal dynamics of the nominal lateral controller of the aircraft through two scalar signals v p and v b which respectively affect roll and sideslip angles dynamics (v x = [v p v b ] and v y = 0, n cr = 2, m r = 0). Remark yet that the second strategy offers more flexibility with the possibility of a direct anti-windup action at the controller output. However, no significant improvement has been observed with this additional feature which has thus not been further considered in this application.

The main objective of this application is design and evaluate anti-windup systems to improve the aircraft response to roll angle solicitations while limiting oscillations despite actuator loss [START_REF] Puyou | Application of robust antiwindup design to the longitudinal aircraft control to cover actuator loss[END_REF]. During such maneuvers, a significant control activity is observed on the ailerons. This is why the effects of saturations are mod-Fig. 5: Detailed view of the magnitude and rate limitations system eled and into account for these actuators in both diagrams of Figures 3 and4 while no saturation is introduced on the rudders. The effects of saturations become even more penalizing in case of a partial loss of control capability. Assume indeed that the aircraft is controlled by a pair of ailerons on each wing but only one is operational. In that case, the activity of the remaining actuators is doubled as well as the risk of magnitude and rate saturations. Then, the magnitude and rate limits in the following will be halved. We will consider L m = 10deg (instead of 20 in normal conditions) and L r = 20deg/s (instead of 40).

In the sequel, five different anti-windup compensators are implemented and compared:

• A standard anti-PIO filter used in the industry. It is an "open-loop" solution which does not exploit the information relative to the saturation of the signal (see [START_REF] Biannic | Analyse et ajustement de lois de commande en présence de saturations implantation de filtres anti-PIO générés par synthèse anti-windup[END_REF]). This may be considered as the basic solution from the industry. It corresponds to the block REFERENCE in Figures 3 and4; • A dynamic H ∞ AW e anti-windup built by using a structured H ∞ design method [START_REF] Puyou | Application of robust antiwindup design to the longitudinal aircraft control to cover actuator loss[END_REF]. The advantage of such a strategy is that it circumvents some limitations of LMI-based strategies (limitation on the size problem when manipulating LMIs, conservatism of sufficient conditions) but to the detriment of the easiness of construction for engineers not always specialists of advanced control theories; • A dynamic AW e anti-windup designed with the Algorithm 3.2 initialized with the H ∞ AW e anti-windup above-described; • A dynamic AW φ anti-windup designed with the Algorithm 3.5 using matrices A aw and C aw borrowed to the H ∞ AW e anti-windup; • A static AW φ anti-windup designed with the Algorithm 3.5. This strategy is an alternative to the standard anti-PIO filter as it is very easy to implement (no additional dynamical system to introduce in the controller block).

Analysis and design of anti-windup AW e

First, an analysis step (in terms of stability performance) of the anti-windup controller proposed in [START_REF] Biannic | Analyse et ajustement de lois de commande en présence de saturations implantation de filtres anti-PIO générés par synthèse anti-windup[END_REF] (denoted H ∞ AW e anti-windup) is carried out. This antiwindup compensator has been built by using a structured H ∞ design method [START_REF] Puyou | Application of robust antiwindup design to the longitudinal aircraft control to cover actuator loss[END_REF]. Such an anti-windup has provided very good numerical results but no proof of its stability in closed-loop was a priori ensured. By using Algorithm 3.1, one can verify that the conditions are feasible. The optimization problem is then solved by considering the bound on perturbation δ = 0.1 and v 1 = [C p (4, :) 0] as the direction to be optimized over the set E (Q -1 , δ ). Algorithm 3.1 gives the optimal solution:

Analysis H ∞ AW e : γ = 1.7167 ; β = 0.6254

It is interesting to do the same analysis for the closed-loop system without antiwindup. The feasibility is also obtained and the solution is: Analysis without anti-widup: γ = 1.8929 ; β = 0.7851

The solution with the H ∞ AW e anti-windup described through γ and β as performance indicators does not appear as much better than the one without anti-windup: indeed γ is actually decreased in the case with anti-windup but β is slightly degraded. However, simulations presented in [START_REF] Biannic | Analyse et ajustement de lois de commande en présence de saturations implantation de filtres anti-PIO générés par synthèse anti-windup[END_REF] exhibited that the time responses with the H ∞ AW e strategy were very close to the desired behavior (that is without saturation), differently from the case without anti-windup resulting in large overshoot and degraded time evolution. The meaning of this is that the considered criterion of optimization, which does not explicitly include the time response performance, does not exactly fit to the analysis or design of the anti-windup loop. Nevertheless, considering criteria on time response performance is a difficult task and the optimization criterion used here gives a reasonable trade-off between stability guarantee, performance and time response.

In a second step, Algorithm 3.2 is used by considering the previous anti-windup controller at the initialization step (Step 1). After one iteration (after the conditions become unfeasible for numerical reasons), one gets a new anti-windup controller AW e such that:

Design AW e : γ = 1.6887 ; β = 0.6978 Now, we compare the results obtained in response to a step demand of 40 deg on the roll angle, by using the scheme given in Figure 3.

In Figure 6, the time evolutions obtained with the H ∞ AW e anti-windup (used in the first step of analysis) and the AW e anti-windup above designed are compared. The case without saturation is also plotted (denoted "reference" in the figure). The time evolution of inputs δ pc is plotted in both cases and without saturation in Figure 7.

Remark 8. Note that the level of performance obtained with the H ∞ AW e anti-windup cannot be much improved with our strategy as this "initial" anti-windup had been cleverly designed with the structured H ∞ approach. The iterative procedure could be initialized with any other anti-windup controller of order n aw , for example (as suggested before) solution of a design in the case of magnitude saturation only.

The initial choice has however a significant influence on the iterative process and the attainable solution. Many tests have shown that the results obtained are often not much convincing (generally the solution extended from the initial H ∞ AW e antiwindup shows better performance indexes) and then it seems that this anti-windup structure is not much adapted to a design from scratch.

Design of anti-windup AW φ

In this section, we consider the second strategy, whose main advantage is that it is based on a systematic method without tuning parameter. Figure 4 shows how the AW φ anti-windup is implemented. As commented in Section 3, in the context of the AW φ strategy, the design of matrices B i aw and D i aw is cast by solving a linear optimization problem. Then, for the same conditions on δ and X 0 as in the previous case, and by considering matrices A aw and B aw of the H ∞ AW e anti-windup, Algorithm 3.5 gives matrices B i aw and D i aw , Similarly, the time evolutions of δ pc in these cases are depicted in Figure 9.

One can observe that the level of performance of the H ∞ AW e anti-windup is slightly degraded in the case of the design of AW φ , but it remains acceptable. Now, we design a static AW φ anti-windup (only matrices D i aw , i = 0, 1, have to be designed). The main advantage is that we do not need to initialize the algorithm as matrices A aw and C aw do not exist (n aw = 0). Algorithm 3.5 provides the following optimal solution: Static AW φ design : γ = 1.7846 ; β = 0.7772 Figures 10 and 11 illustrate the time evolution of the closed-loop system to a roll solicitation of 40 deg. The responses are compared by considering the case without saturation, a standard anti-PIO strategy (see [START_REF] Biannic | Analyse et ajustement de lois de commande en présence de saturations implantation de filtres anti-PIO générés par synthèse anti-windup[END_REF]) and the static AW φ anti-windup strategy. It is important to underline that a simple static anti-windup strategy allows to obtain better performance than the standard anti-PIO case, which adds dynamics in the system.

Conclusion

In this chapter, an anti-windup analysis and design strategy has been proposed for systems involving both magnitude and rate saturations, and taking into consideration that such saturations elements only affect some of the inputs. Such a situation has much practical interest for many systems issued in particular from aerospace domain. It is illustrated here on a lateral flying model of a civil aircraft in presence of aggressive maneuvering of the pilot. Actually magnitude and rate saturations of the ailerons deflection actuator may lead to an undesirable behavior which is often called Pilot-Induced-Oscillation (PIO). For this class of nonlinear control systems, anti-windup compensators have been adapted through adequate convex optimization schemes. A comparison with given dynamic anti-PIO filters already used for this class of systems has also been provided. This work lets many questions open, such as the design of other anti-windup schemes. Other classes of fruitful anti-windup compensators may include the parameter-varying approach [START_REF] Roos | On-ground aircraft control design using a parameter-varying anti-windup approach[END_REF] or reset controllers [START_REF] Tarbouriech | Anti-windup strategy for reset control systems[END_REF].

Fig. 1 :

 1 Fig. 1: Principle of anti-windup

Fig. 2 :

 2 Fig. 2: (a) Actuator with rate and magnitude saturations. (b) Model used to represent such an actuator (scalar case)

Fig. 3 :Fig. 4 :

 34 Fig. 3: Nonlinear closed-loop Simulink implementation of Anti-Windup AWe for lateral aircraft simulations

Fig. 6 :

 6 Fig. 6: Roll solicitation of 40 deg: comparison of the performance outputs for the cases with H ∞ AW e anti-windup and designed anti-windup AW e

Fig. 7 :

 7 Fig. 7: Roll solicitation of 40 deg: comparison of the inputs for the cases with H ∞ AW e anti-windup and designed anti-windup AW e

Fig. 8 :

 8 Fig. 8: Roll solicitation of 40 deg: comparison of the performance outputs for the cases with H ∞ AW e anti-windup and designed anti-windup AW φ

Fig. 9 :Fig. 10 :Fig. 11 :

 91011 Fig. 9: Roll solicitation of 40 deg: comparison of the inputs for the cases with H ∞ AW e anti-windup and designed anti-windup AW φ

* This work was supported by COCKPIT project, convention ONERA F/20 334/DA PPUJ.


Vinh Quang

Claudia-Lavinia Dang

Ignat

Quality Assessment of Wikipedia Articles: A Deep Learning Approach

Wikipedia is indeed a very important knowledge sharing platform. However, since its start in 2001, the quality of Wikipedia is questioned because its content is created potentially by everyone who can access to the Internet. Currently, the quality of Wikipedia articles is assessed by human judgement. The method is not scalable up to huge size and fast changing speed of Wikipedia today. An automatic quality classifier for Wikipedia articles is required to support users to choose high quality articles for reading and to notify authors for improving their products. While other existing approaches are based on manually predefined specific feature set, we present our approach of using deep learning to automatically represent Wikipedia articles for quality classification.

INTRODUCTION

Today, Wikipedia is the largest 1 and probably the most important knowledge repository in the world. At the time of writing, there are more than five millions articles in English Wikipedia particularly, and around 40 millions article in Wikipedia all languages 2 . On average, ten edits per second are performed on Wikipedia 3 .

The question related to Wikipedia quality was raised since its start in 2001 [START_REF] Denning | Wikipedia risks[END_REF]. Some studies [START_REF] Holman Rector | Comparison of Wikipedia and other encyclopedias for accuracy, breadth, and depth in historical articles[END_REF] claimed that the quality of Wikipedia is not comparable to other traditional encyclopedias.

In order to improve the overall quality of Wikipedia, Wikipedia articles need to be classified based on their quality such that readers can be guided to high quality writing while authors can be notified about low quality texts that require improvements. Several quality classes were defined, from FA as the highest quality class to Stub as the lowest quality class as shown in Table I. Currently, the quality labels are assigned to articles by human judgement 4 . However, due to the huge size of Wikipedia and fast speed of edits, humans cannot manually review them. Therefore, an automatic quality classifier is required.

Existing approaches rely on feature engineering requiring a manual specification of a fea-Class Description FA Professional, outstanding, and thorough; a definitive source for encyclopedic information. GA Useful to nearly all readers, with no obvious problems; approaching (but not equalling) the quality of a professional encyclopedia. B Readers are not left wanting, although the content may not be complete enough to satisfy a serious student or researcher.

C

Useful to a casual reader, but would not provide a complete picture for even a moderately detailed study.

Start

Provides some meaningful content, but most readers will need more.

Stub

Provides very little meaningful content; may be little more than a dictionary definition. Readers probably see insufficiently developed features of the topic and may not see how the features of the topic are significant.

Table I. Description of quality classes of Wikipedia articles 5 . ture set. However, there is no guarantee that the feature set is complete and usually selection of features depends on researchers' expertise.

RELATED WORKS

Existing approaches rely on specific feature set to determine quality of Wikipedia articles.

In other words, each approach proposed a set of features that are believed for correlation with articles quality.

One of the first efforts belongs to the work of [START_REF] Blumenstock | Size matters: word count as a measure of quality on wikipedia[END_REF] where authors used length of articles as an indicator for their quality. [START_REF] Dalip | Automatic quality assessment of content created collaboratively by web communities: a case study of wikipedia[END_REF]] analyzed the effect of the feature set comprising text, review and network on the quality of Wikipedia articles. [START_REF] Anderka | Predicting quality flaws in user-generated content: the case of wikipedia[END_REF] Several studies defined features based on information about authors rather than about articles themselves. [START_REF] Adler | Assigning trust to wikipedia content[END_REF]] used authors reputation as feature set. [START_REF] Suzuki | Quality assessment of Wikipedia articles using h-index[END_REF]] applied the concept of h-index in academic publication for determining quality of However, the definition of the feature set used by existing approaches for automatic classification of Wikipedia articles quality is based on the expertise, experience and knowledge of each research team. There is no gold standard theory to find the best feature set. When a new feature is defined, the only validation method is empirical testing. For instance, [Halfaker and Taraborelli 2015] used the division of number of images by article length as a feature. But there is no evidence that adding some derived features such as division of number of images by square of content length, or by square root of content length would change the performance of the classifier.

A DEEP-LEARNING BASED APPROACH

While traditional machine learning techniques such as k-NN or random forest [START_REF] Kubat | An Introduction to Machine Learning[END_REF]] require manual feature engineering, deep learning techniques can be used to learn features automatically from the dataset rather than defining them before-hand [START_REF] Lecun | Deep learning[END_REF]].

We proposed an approach that uses Doc2Vec [ [START_REF] Le | Distributed representations of sentences and documents[END_REF] for learning features from textual documents and Deep Neural Networks [START_REF] Goodfellow | Deep Learning[END_REF] for classifying Wikipedia articles [START_REF] Dang | Quality assessment of wikipedia articles without feature engineering[END_REF].

Method

We used Doc2Vec [START_REF] Le | Distributed representations of sentences and documents[END_REF] and Deep Neural Networks [START_REF] Goodfellow | Deep Learning[END_REF]. Our proposed approach includes two steps. In the first step we applied Doc2Vec on Wikipedia articles. The idea of Doc2Vec is to convert a variable-length textual document into a fixed-length numerical vector. Instead of applying Doc2Vec on textual documents as the original method presented in [START_REF] Le | Distributed representations of sentences and documents[END_REF], we applied Doc2Vec on raw content of Wikipedia articles. This raw content contains all the necessary information to determine the quality of an article. An example of a raw Wikipedia content is provided in Fig. 1.

In the second step we applied Deep Neural Networks (DNN) on output vectors of Doc2Vec. DNN is defined as an artificial neural networks [START_REF] Mcculloch | A logical calculus of the ideas immanent in nervous activity[END_REF] with multiple hidden layers [START_REF] Goodfellow | Deep Learning[END_REF], as visualized in Fig. 2 7 . DNN receives input values from its input layer, passes the computation through hidden layers and generates predicting values at output layer as follows: the k th layer computes an output vector h k Fig. 2. An example of deep neural networks with three hidden layers using the output h k-1 of the previous layer, starting with the input layer x = h 0 as in Equation 1, where b k is the offset vector and W k is the matrix of weights.

h k = f (b k + W k h k-1 ) (1)
The function f in Equation 1 is called activation function, which decides how each neuron calculates and transfers the signal to the neurons in subsequent layer. We used rectifier as our activation function, defined as follows:

rectif ier(x) = max(0, x) (2)
The values of b and W are set during training process, using gradient-descent technique [START_REF] Bengio | Practical recommendations for gradient-based training of deep architectures[END_REF]].

DNN is an emerging research field, and there are still a lot of open theoretical questions, especially in choosing hyper-parameters for the network [START_REF] Goodfellow | Deep Learning[END_REF]]. In our experiments, we came up with a neural network with four hidden layers as the result of random search technique [START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF].

Results

We validated our approach on the public dataset including 30 000 Wikipedia articles provided by Wikimedia Foundation Research8 . Following the train/test division of existing studies [START_REF] Warncke-Wang | The success and failure of quality improvement projects in peer production communities[END_REF][START_REF] Halfaker | Artificial intelligence service gives Wikipedians 'x-ray specs' to see through bad edits[END_REF], we used 80% of the dataset for training and 20% for testing.

For the classification of Wikipedia articles according to all six quality classes, our method achieved the accuracy of 55% compared to the accuracy of 58% of [START_REF] Warncke-Wang | The success and failure of quality improvement projects in peer production communities[END_REF]] and 60% of [START_REF] Halfaker | Artificial intelligence service gives Wikipedians 'x-ray specs' to see through bad edits[END_REF]. We believe further optimisations can be achieved using the DNN technique.

• 3 Fig. 1 .

 31 Fig. 1. An example of raw Wikipedia content

  Wang et al., with an adaption for French Wikipedia. Finally, Wikimedia Foundation development team used a feature set which includes 24 features for English Wikipedia, and 25 features for French Wikipedia. To our knowledge, Wikimedia ORES service is the state-of-the-art approach in classifying the quality of Wikipedia articles.

built a classifier to detect vandalism on Wikipedia based on cleanup tags. [Warncke-Wang et al. 2013] analyzed a set of 17 features and eventually came up with 11 features and applied random forest technique for the classification. Results are discussed in more details in [Warncke-Wang et al. 2015]. Based on the work of [Warncke-Wang et al. 2015; Warncke-Wang et al. 2013], Wikimedia Foundation 6 built an online API to predict the quality class of Wikipedia articles called ORES(Objective Revision Evaluation Service) [Halfaker and Taraborelli 2015]. The service uses an extended version of the feature set presented by Warncke-

https://wikimediafoundation.org SIGWEB Newsletter Autumn 2016

The image is from[START_REF] Nielsen | Neural networks and deep learning[END_REF]] with the permission of usage.SIGWEB Newsletter Autumn 2016

http://datasets.wikimedia.org/public-datasets/enwiki/ SIGWEB Newsletter Autumn 2016

SIGWEB Newsletter Autumn 2016

ACKNOWLEDGMENTS

The authors would like to thank Dr. Aaron Halfaker, Wikimedia Research Foundation and Morten-Warncke Wang, University of Minnesota for providing the dataset and valuable discussions.

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).

• 5

For the binary classification, our method achieved very good results. The method classifies FA vs Start with the accuracy of 99%, FA-GA vs all with the accuracy of 86% and FA-GA vs C-Start of 90%, much higher than existing approaches. For instance, [START_REF] Xu | Measuring article quality in Wikipedia: Lexical clue model[END_REF] classifies FA vs Start with the accuracy of 84%, [START_REF] Lex | Measuring the quality of web content using factual information[END_REF]] classifies FA-GA vs all with the accuracy of 84%, and [START_REF] Wu | Classifying wikipedia articles using network motif counts and ratios[END_REF]] classifies FA-GA vs C-Start with the accuracy of 66%.

CONCLUSIONS

Feature engineering is a norm in machine learning for a long time [START_REF] Kubat | An Introduction to Machine Learning[END_REF]] with a lot of hands-on tips which are accumulated through a lot of works [START_REF] Zheng | Mastering Feature Engineering: Principles and Techniques for Data Scientists[END_REF]]. The recent research on deep learning brings a new potential to machine learning field to build a system which can learn to represent the data by itself. Our work presented a new way to understand the problem of Wikipedia articles classification, which is an important topic in both research and practice.


Thomas Lacroix

Mansour Karami

Avery E Broderick

Joseph Silk

Céline Boehm

A unique probe of dark matter in the core of M87 with the Event Horizon Telescope

Keywords: numbers: 95.35.+d, 96.50.S-, 95.55.Br

We demonstrate the unprecedented capabilities of the Event Horizon Telescope (EHT) to image the innermost dark matter profile in the vicinity of the supermassive black hole at the center of the M87 radio galaxy. We present the first model of the synchrotron emission induced by dark matter annihilations from a spiky profile in the close vicinity of a supermassive black hole, accounting for strong gravitational lensing effects. Our results show that the EHT should readily resolve dark matter spikes if present. Moreover, the photon ring surrounding the silhouette of the black hole is clearly visible in the spike emission, which introduces observable smallscale structure into the signal. We find that the dark matter-induced emission provides an adequate fit to the existing EHT data, implying that in addition to the jet, a dark matter spike may account for a sizable portion of the millimeter emission from the innermost (sub-parsec) region of M87. Regardless, our results show that the EHT can probe very weakly annihilating dark matter. Current EHT observations already constrain very small cross-sections, typically down to a few 10 -31 cm 3 s -1 for a 10 GeV candidate, close to characteristic values for p-wave-suppressed annihilation. Future EHT observations will further improve constraints on the DM scenario.

I. INTRODUCTION

The dark matter (DM) density profile at the centers of galaxies is critical to indirect searches but remains poorly constrained. In objects such as M87, the DM profile may be significantly enhanced on sub-parsec scales by the central supermassive black hole (SMBH), although there is no direct evidence for such a sharp density increase, referred to as a DM spike.

DM spikes however leave a very distinctive signature in the spectral energy distribution (SED) of a galaxy if DM can annihilate [1]. More specifically, in the case of M87, where a spike should plausibly have formed and survived galaxy dynamics, an anomalous contribution to the SED is expected when the DM annihilation cross-section is larger than ∼ 10 -30 cm 3 s -1 for light (10 to 100 GeV) DM particles, and above ∼ 3 × 10 -26 cm 3 s -1 for candidates as heavy as 100 TeV.

Here we show that it is possible to probe even fainter DMinduced radiation in M87 by using the spatial morphology of the DM-induced synchrotron emission near the central black hole (BH). Due to a lack of angular resolution in existing observational facilities, such a study of the DM-induced signal * lacroix@iap.fr in the inner part of M87 has not been performed yet-nor has it been done in similar objects. However, this is now possible with the advent of the Event Horizon Telescope (EHT).

Heavy DM particles annihilating into Standard Model particles near the central BH are expected to produce synchrotron emission in the frequency range that is currently being probed by the EHT. Here, we show that the synchrotron halo induced by DM is bright enough to be resolved by the EHT, if there is a DM spike. Moreover, this additional radiation enhances the photon ring around the BH shadow, thus making it a prominent feature in the EHT data and a new probe of the DM properties.

In Sec. II, we provide a description of the EHT, then we present our DM model in Sec. III and the EHT data we used in Sec. IV. Our results can be found in Sec. V, and we conclude in Sec. VI.

II. THE EVENT HORIZON TELESCOPE A. General features

The EHT is a global network of millimeter and submillimeter facilities that employs very long baseline interferometry (VLBI) to create an effective Earth-scale high angular resolution telescope [START_REF] Doeleman | 10th European VLBI Network Symposium and EVN Users Meeting: VLBI and the New Generation of Radio Arrays[END_REF][START_REF] Lu | [END_REF]. The purpose of this array is to test general relativity and shed light on physical processes tak-ing place in the vicinity of SMBHs at the centers of galaxies. To date, EHT data for M87 has been reported for a 3-station array comprised of the Submillimeter Telescope (SMT) in Arizona, the Combined Array for Research in Millimeterwave Astronomy (CARMA) in California, and a network of three facilities in Hawaii: the James Clerk Maxwell Telescope (JCMT), the Submillimeter Array (SMA), and the Caltech Submillimeter Observatory (CSO). This configuration has already achieved an impressive angular resolution of order 40 µas at 230 GHz. Presently, the Atacama Large Millimeter/submillimeter Array (ALMA) in Chile, the Large Millimeter Telescope (LMT) in Mexico, the Institut de Radioastronomie Millimétrique (IRAM) 30m, the Plateau de Bure interferometer, and the South Pole Telescope (SPT) 1 will be added to the EHT. Longer term, the Greenland Telescope will join the array. Combined, the EHT will directly access angular scales as small as 26 µas at 230 GHz and 17 µas at 345 GHz.

The angular scale of the Schwarzschild radius of the SMBH at the center of M87 is about 8 µas. 2 The strong gravitational lensing near the black hole (BH) magnifies this by a factor of as much as 2.5, only weakly dependent of BH spin, making M87 a primary target for the EHT [START_REF] Lu | [END_REF].

B. Imaging the shadow of a black hole

Characteristic in all images of optically thin emission surrounding BHs is a "shadow"-a dark central region surrounded by a brightened ring, the so-called "photon ring". This is a direct consequence of the strong gravitational lensing near the photon orbit, and is directly related to the projected image of the photon orbit at infinity. The shadow interior is the locus of null geodesics that intersect the horizon, and thus does not contain emission from behind the BH. The bright ring is in general sharply defined due to the instability of the photon orbit and a consequence of the pileup of higher-order images of the surrounding emission. For a Schwarzschild BH the radius of this shadow is r shadow = 3 √ 3/2R S ≈ 2.6R S [START_REF] Bardeen | Black Holes[END_REF][START_REF] Bozza | [END_REF]. For Kerr BHs, this radius ranges from 2.25R S to 2.6R S , deviating substantially only at large values of the dimensionless spin parameter and viewed from near the equitorial plane, i.e., a 0.9 and θ 60 • [9].

The generic appearance of the shadow, its weak dependence on BH spin, and fundamentally general relativistic origin make it a prime feature in EHT science. Also imprinted on EHT images will be the high-energy astrophysics of the nearhorizon region; the physics of BH accretion and relativistic jet formation. Horizon-scale features have already been observed in early EHT observations of Sagittarius A* (Sgr A*) [10][11][12] and M87 [13,14], demonstrating that such structure exists. Here, we assess the observability of the shadow of the SMBH at the center of M87 in the electromagnetic signal from DM annihilation, and the limits that may be placed on DM properties given such a scenario.

III. PROBING A DARK MATTER SPIKE AT THE CENTER OF M87 WITH THE EVENT HORIZON TELESCOPE

In the context of the observational opportunities offered by the EHT, we now discuss the potential of this instrument in terms of DM searches. In particular, we study the observability of a DM spike at the center of M87, since such a sharply peaked morphological feature is expected to yield strong annihilation signals. At the frequencies of interest for the EHT, typically a few hundred GHz, the main DM signature comes from synchrotron radiation. Therefore, in order to assess the ability of the EHT to probe the inner part of the DM profile of M87, we need to compute the synchrotron emission of electrons and positrons produced in DM annihilations in the inner region. We make the following assumptions:

• The presence of a SMBH at the center is likely to lead to fairly strong magnetic fields, typically around 10-10 2 G [15]. As a result of such strong magnetic fields, synchrotron radiation and advection towards the central BH are the dominant physical processes by which DM-induced electrons and positrons lose or gain energy [16,17], whereas inverse Compton scattering and bremsstrahlung are negligible. Additionally, the timescales associated with synchrotron radiation and advection are much shorter than that of spatial diffusion, so we disregard the latter in the following. Note that even larger magnetic fields-up to ∼ 10 3 G-would arise if the equipartition scenario proposed in Ref. [18] for the center of the Milky Way were realized in M87. Therefore, to account for the uncertainty on the central magnetic field strength, we consider values in the range 10-10 3 G.

• We also disregard two processes that can reduce the synchrotron intensity. On the one hand, the synchrotron self-Compton effect, which would lead to additional energy losses for electrons and positrons, is only relevant for magnetic fields significantly smaller than 0.1 G [16].

For larger magnetic fields, such as the ones we consider here, synchrotron self-Compton losses are negligible with respect to synchrotron losses. On the other hand, synchrotron self-absorption is only relevant below ∼ 10 GHz [16,17], so it can be neglected for the EHT frequency of 230 GHz.

• We assume the presence of a spike in the DM profile, formed through adiabatic growth of a SMBH [19] at the center of a DM halo with power-law density profile.

The existence of such a strong enhancement of the DM density-corresponding to ρ(r) ∝ r -γ sp with typically γ sp = 7/3-is debated. It can actually be weakened by various dynamical processes such as mergers [20] (although the resolution of numerical simulations used to reach that conclusion is not sufficient to set strong constraints on the profile on sub-pc scales), if the BH did not grow exactly at the center of the DM halo [21,22], or if the BH growth cannot be considered adiabatic [22]. Moreover, dynamical heating in the central stellar core would also soften a spike. However, an adiabatic spike is more likely to have survived in a dynamically young galaxy such as M87, in which stellar heating is essentially negligible, as discussed in Refs. [1,23]. Furthermore, other dynamical processes can have the opposite effect of making the survival of a spike more likely; for instance enhanced accretion of DM to counteract the depopulation of chaotic orbits in triaxial halos [24]. Therefore in the following, we assume that an adiabatic spike effectively formed at early times at the center of M87 and has survived until today.

A. Electron propagation in the presence of advection

To derive the intensity of DM-induced synchrotron radiation, we first need to compute the electron and positron spectra from the DM annihilation rate. This is done by solving the propagation equation of DM-induced electrons and positrons which, in the presence of synchrotron radiation and advection, and assuming a steady state reads (see e.g. Refs. [16,17]):

v ∂ f i ∂r - 1 3r 2 ∂ ∂r r 2 v p ∂ f i ∂p + 1 p 2 ∂ ∂p p 2 ṗ f i = Q i , (1)
where f i (r, p) is the distribution function of electrons and positrons in momentum space, at radius r and momentum p, for annihilation channel i. The first, second and third terms correspond to the advection current, the energy gain of electrons due to the adiabatic compression, and the loss term due to radiative losses, respectively. v(r) = -c (r/R S ) -1/2 is the radial infall velocity of electrons and positrons onto the BH in the accretion flow, where R S is the Schwarzschild radius. The minus sign in the expression of the inflow velocity accounts for the direction of the flow, oriented towards the BH. The source function Q i (r, p) is the DM annihilation rate in momentum space for channel i, related to the annihilation rate q i (r, E) in energy space via

Q i (r, p) = c 4πp 2 q i (r, E) (2)
in the ultra-relativistic (UR) regime where E = pc. The UR approximation can be safely used for electrons and positrons for the energy range relevant for this study. The usual annihilation rate in energy space reads:

q i (r, E) = σv i η ρ(r) m DM 2 dN e,i dE (E), (3)
where η = 2 for the case of self-annihilating DM that we consider here. The injection spectrum dN e,i /dE is taken from Ref. [25] and the associated website.

Since radiative losses are dominated by synchrotron losses, the total radiative loss term ṗ = dp/dt reduces to the synchrotron loss term [START_REF] Longair | High Energy Astrophysics[END_REF]

ṗ(r, p) = ṗsyn (r, p) = - 2σ T B 2 p 2 3µ 0 (m e c) 2 . ( 4
)
We assume the intensity of the magnetic field to be homogeneous, i.e. B ≡ B 0 , over the accretion region, which has a size r acc corresponding to the sphere of influence of the BH [17], so typically ∼ 60 pc (as discussed in Ref. [1]), which is also roughly the size of the spike. The resolution of the propagation equation, Eq. ( 1), in the presence of synchrotron losses and advection, in the UR regime, and with the method of characteristics, yields the electron and positron spectrum in terms of the DM annihilation rate [16]

f i (r, p) = 1 c r R S r acc r Q i (R inj , p inj ) R inj R S 5 2 p inj p 4 dR inj , (5)
where the injection momentum p inj ≡ p inj (R inj ; r, p) for a homogeneous magnetic field is given in Appendix B. From there, the electron and positron energy spectrum is given by:

ψ i (r, E) = 4πp 2 c f i (r, p). (6)
We then convolve ψ i with the synchrotron power P syn (ν, E, r) to obtain the synchrotron emissivity:

j syn,i (ν, r) = 2 m DM m e P syn (ν, E, r)ψ i (r, E) dE. (7)

B. Relevance of advection

Advection shapes the inner part of the intensity profile by displacing electrons and positrons towards the BH, thereby accelerating them. This effect is in competition with synchrotron losses. Therefore, depending on the magnetic field, electrons either lose their energy in place through synchrotron radiation, or are first advected towards the center. The dependence on the magnetic field of the size of the region where electrons are affected by advection is obtained by comparing the synchrotron loss term given in Eq. ( 4) with the momentum gain rate due to adiabatic compression, ṗad = -

1 3r 2 ∂ ∂r r 2 v(r) p. ( 8
)
The shape of the emissivity profile is thus governed by advection for:

r   3R 1 2 S µ 0 m 2 e c 4 4σ T B 2 0 E syn   2 3 , (9)
where E syn = 4πm 3 e c 4 ν/(3eB 0 ) 1/2 is the peak synchrotron energy at the frequency ν of interest [START_REF] Longair | High Energy Astrophysics[END_REF]. In terms of the angular distance from the center θ ≡ r/d M87 -where d M87 ≈ 16 Mpc is the distance between Earth and the center of M87this condition reads, for ν = 230 GHz,

θ 7 B 0 10 G -1 µas. (10)
It is therefore essential to include the advection process for B 0 10 G, since in that regime it has a strong impact on the synchrotron intensity in the region of interest for the EHT. For B 0 10 G, advection is negligible since it would only dominate for radii smaller than the Schwarzschild radius of the BH.

C. Dark synchrotron intensity in a curved spacetime

The DM-induced synchrotron intensity for a flat spacetime is computed by integrating the emissivity over the line of sight. However, the actual spacetime accounting for the presence of the SMBH at the center is characterized by the Schwarzschild or Kerr metric respectively if the BH is static or rotating. In these realistic cases, the correct spatial morphology of the synchrotron intensity I ν is obtained using a raytracing technique that accounts for the gravitational lensing effect due to the BH. In our case, this is achieved via the raytracing and radiative transfer scheme described in Refs. [START_REF] Broderick | [END_REF][28][29].

IV. DATA

Here we provide a brief overview of the EHT data used to constrain our models. Details regarding individual observing runs, data calibration, and the uncertainty estimates can be found in Refs. [13,14], to which we direct interested readers.

The EHT, like all interferometers, directly constructs visibilities by cross-correlating observations at pairs of stations. These are directly proportional to the Fourier transform of the image at a spatial frequency to the ratio of the projected baseline distance between the two stations to the observation wavelength. Throughout the night the rotation of the Earth results in a rotation of the projected baseline, changing both its orientation and length, thereby generating a moderate variation in the spatial frequencies probed by any particular pair of sites.

Importantly, like any interferometer the EHT acts as a "high-pass" filter, sensitive primarily to structures on angular scales that lie between λ/u short and λ/u long , where u short and u long are the shortest and longest baselines in the array, respectively. Within the context of the published EHT data on M87 these angular scales are 75 µas to 450 µas, though the high signal-to-noise ratio of the data extends these by roughly a factor of two. Therefore, features that extend over more than a milliarcsecond are effectively invisible to the EHT.

Upon measuring a sufficient number of these visibilities an image can be produced via inverting the Fourier transform. In practice, this is performed via a number of sophis-ticated image-inversion techniques that impose additional requirements on the final image, e.g., positivity, smoothness, etc [30]. However, because the published EHT observations only sparsely sample the spatial-frequency plane (often called the "u-v" plane) we compare directly with the measured visibilities.

By construction the visibilities are complex valued, and therefore described by an amplitude and phase. However, in practice the amplitudes are known much better than the phases as a result of the typically large, and highly variable, atmospheric phase delays. This does not mean that phase information is completely unavailable; "closure phases" constructed from triplets of sites, equal to the sum of the phases over the closed triangle of baselines, are insensitive to site-specific phase errors. Therefore, the two data sets we employ consist of visibility amplitudes and closure phases obtained with the Hawaii-SMT-CARMA array, reported in Refs. [13,14]. In all cases observations were taken at 230 GHz.

Altogether, 160 visibility amplitudes were constructed from data taken on April 5, 6, and 7, 2009, reported in Ref. [13], and March 21, 2012, reported in Ref. [14]. All nights showed visbility amplitudes consistent with a single source structure. Of these 54 visibilities were reported on the CARMA-SMT baseline, 83 on the Hawaii-CARMA baseline, and 23 on the Hawaii-SMT baseline. Note that in all of the reported measurements a 5% systematic calibration error has been added to the uncertainties in quadrature.

On the Hawaii-CARMA-SMT triangle, 17 closure phases were constructed from data taken on March 21, 2012 and reported in Ref. [14]. Where visibility amplitudes provide a measure of the "power" in an image at a given spatial scale, closure phases are particularly sensitive to asymmetry; e.g., a point-symmetric image has identically zero closure phases. These are consistent with a constant closure phase of 0 • , with typical uncertainties of 10 • .

V. RESULTS

Predicted images are shown in Fig. 1, corresponding to the synchrotron intensity at 230 GHz from DM. The upper panels correspond to a DM spike with γ sp = 7/3, whereas the maps in the lower panels are computed for the no-spike case, for which we consider a standard Navarro-Frenk-White (NFW) DM profile with power-law index γ = 1. We assume annihilation of 10 GeV DM particles into b b, 3 and a magnetic field of 10 G, for a static BH (left panels) and a maximally rotating BH (right panels).

The photon ring, i.e. the bright ring of radius ∼ 20 µas that surrounds the darker shadow of the BH, is clearly visible in the simulations for all DM models we consider, although in practice in the absence of a spike the signal is too weak to be FIG. 1. Simulated maps of the synchrotron intensity at 230 GHz from a spike of 10 GeV DM annihilating into b b, accounting for the strong gravitational lensing induced by the central BH, for a Schwarzschild BH (left panel) and a maximally rotating BH (right panel), in the presence (upper panels) and absence (lower panels) of a spike in the DM profile. Note that considering the wide range of intensities, we use different color scales, but with the same dynamic range spanning three orders of magnitude for comparison. The angular coordinates ξ and η correspond to the directions respectively perpendicular and parallel to the spin of the BH. For the spike cases, the slope of the DM spike is γ sp = 7/3, and the annihilation cross-sections correspond to the best fit to EHT observations (see text for details), namely 7.4 × 10 -31 cm 3 s -1 for the Schwarzschild case and 3.1 × 10 -31 cm 3 s -1 for the maximally rotating case. In the absence of a spike, the intensity is computed for the thermal s-wave cross-section of 3 × 10 -26 cm 3 s -1 . For all the simulated maps the magnetic field is 10 G.

detectable with the EHT, as discussed in the following. The presence of a photon ring introduces small-scale structure into the signal, readily observable with the EHT on long baselines. For a static BH the shadow is exactly circular. For all but the most rapidly rotating BHs it is also very nearly circular [31,32]. For a maximally rotating Kerr BH viewed from the equatorial plane the photon ring is flattened in the direction aligned with the BH spin.

At scales above 25 µas the DM spike-induced emission produces a diffuse synchrotron halo whose intensity falls with radius as a power law with index ≈ 3.5. This is generic, occuring independently of the BH spin and is present even when gravitational lensing is ignored. The extended nature of this component ensures that it is subdominant on Earth-sized baselines. In the absence of a spike, the profile is much flatter, and falls with radius as a power law with index ≈ 1. Fig. 1 also illustrates the fact that the intensity is significantly enhanced in the presence of a DM spike with respect to the no-spike case. To better stress this enhancement, we show the maps for the spike case computed for very small annihilation crosssections of a few 10 -31 cm 3 s -1 -corresponding to the best fits to the EHT data, as discussed below-, while in the absence of a spike we use the thermal s-wave cross-section of 3 × 10 -26 cm 3 s -1 . Shown in the blue solid line in Fig. 2 is the visibility amplitude at 230 GHz as a function of baseline length for the current EHT triangle, for the simulated DM-induced synchrotron signal, computed with a cross-section that gives the best fit to the EHT measurements from Refs. [13,14]. The left panel corresponds to the Schwarzschild case and the right panel to the maximally rotating case. The DM-induced visibility amplitudes shown in Fig. 2 correspond to the intensity maps shown in Fig. 1. Note that no additional astrophysical component has been included in these.

A standard NFW cusp actually results in visibility amplitudes that are about eight orders of magnitude lower than the EHT data. Note that these no-spike visibility amplitudes are not shown in Fig. 2 for clarity. Therefore, the EHT is only sensitive to spiky profiles, which makes it a dedicated probe of such sharply peaked DM distributions. 4As shown in Fig. 2, a spike of annihilating DM gives a good fit to the EHT measurements of the visibility amplitudes, with best fit cross-sections and reduced chi-squaredχ 2 red = χ 2 /d.o.f., with a number of degrees of freedom (d.o.f.) equal to 160 (data points) minus one (cross-section)-given in Tables I andII. While the fits appear by eye to be quite good, the reduced chi-squareds coupled with the large number of degrees of freedom result in a p-value < 0.002, implying that some structural component is missing in our model. This is not, in itself, surprising given the extraordinary simplicity of the DM spike model and our neglect of the contributions from the observed larger-scale radio emission associated with the jet.

The morphology of the predicted visibility amplitudes is only weakly sensitive to changes in the DM mass, the annihilation channel or the magnetic field, resulting primarily in different best-fit cross-sections. A small increase in the ratio of visibilities on long baselines to those on short baselines for larger magnetic fields arises from the higher intrinsic synchrotron peak coupled with the gravitational redshift, which increases the 230 GHz emission near the horizon.

The visibility amplitudes are more sensitive to the characteristics of the BH. Rotating BHs exhibit larger relative visibility amplitudes on long baselines than static BHs as a consequence of the bright emission from comparatively smaller radii. As a result, the quality of the fit is much better for a Schwarzschild BH.

The closure phases for the DM-spike-induced signal are shown in Fig. 3, for the Schwarzschild case (solid line) and the maximally rotating case (dashed line). The symmetry of the simulated signal is insensitive to the various parameters so we do not need to specify them here. As expected, the closure phase is identically zero for a Schwarzschild BH, while it is slightly larger than zero for the maximally rotating case. In both cases, closure phases for the DM-induced emission are consistent with the low closure phases observed.

Small closure phases are also typical of astrophysical models on the Hawaii-CARMA-SMT triangle [14]. This is because of the near degenerate nature of the projected baseline triangle due to the comparatively short CARMA-SMT baseline; closure phases on trivial triangles (in which one baseline has zero length) vanish identically. However, the inclusion of a number of additional sites in the near future will result in many additional, open triangles for which the closure phases are likely to differ substantially from zero [14]. These will be instrumental to discriminating both between astrophysical and DM-dominated models.

VI. CONCLUSION

In this paper, we have demonstrated the potential of the EHT for DM searches. We have presented the first model of the DM-induced synchrotron emission in the close vicinity of a SMBH accounting for strong gravitational lensing effects. Our conclusions are the following.

• The synchrotron emission from DM spikes should be readily visible in EHT images of M87 if present. This remains true even for very small values of the annihilation cross-section. The resulting emission follows the structure of the DM spike, resulting in a synchrotron halo extending from horizon scales to roughly 100 µas.

• Within the spike emission, the silhouette of the BH is clearly visible on small scales for all spins and models we considered. This imparts small-scale structure on the image on scales of 50 µas, and contributes substantially to the visibilities on long baselines.

• DM spike emission provides an adequate fit to the existing horizon-scale structural constraints from the EHT. This necessarily ignores astrophysical contributions associated with the jet launching region. Such an additional astrophysical component is strongly motivated by the existence of extended emission at wavelengths of 3 mm and longer. Even within the 1.3 mm data, the fit quality suggests that the simple spike structures we present here are incomplete.

• Nevertheless, the limits on M87's flux and small-scale structure place corresponding constraints on the putative DM annihilation cross-sections. For a 10 GeV DM candidate this cross-section must be less than a few 10 -31 cm 3 s -1 , close to the characteristic cross-sections for p-wave-suppressed annihilation. The introduction of additional astrophysical components would decrease this limit further.

• EHT observations in the near future will include a number of additional stations, enabling the reconstruction of M87's image with substantially higher fidelity. As a result, the limits on the existence of DM spikes and the properties of DM candidates will similarly improve. Thus, the EHT opens a new, powerful path to probing the structure and features of DM in the centers of galaxies. Here we describe how we normalize the profile corresponding to a DM spike growing adiabatically from an initial power-

v f ∂r - 1 3r 2 ∂ ∂r r 2 v p ∂ f ∂p + 1 p 2 ∂ ∂p p 2 ṗ f = Q, ( B1
)
can be rewritten as:

∂ f ∂r + ṗad + ṗsyn v ∂ f ∂p = - 1 vp 2 ∂ ∂p p 2 ṗsyn f + Q v , ( B2
)
where ṗad is the momentum gain rate due to adiabatic compression in the advection process, and v(r) the velocity field of the accretion flow. The associated characteristic curves are obtained by solving the following differential equation:

dp dr = ṗad + ṗsyn v . ( B3
)
Generalizing the method of Ref. [16] to an arbitrary powerlaw profile for the magnetic field, B(r) = B 0 (r/R S ) -α/2 , solv- where p inj ≡ p inj (R inj ; r, p) is the injection momentum of an electron injected at R inj (≥ r) and arriving at r with momentum p. Using Eq. (B4) and expressing p inj as a function of p, we obtain, for α = 0:

p inj (R inj ; r, p) = p   k 0 R -1 2 S c R 3 2 inj p r R inj -1 + R inj r 1 2   -1
.

(B7) Note that the denominator of p inj can vanish and become negative, leading to nonphysical values of the injection momentum. This is related to the efficiency of the accretion flow and characterizes the region of the injection parameters (R inj , p inj ) corresponding to a given arrival point (r, p). In practice, p inj remains positive for R inj < R 0 inj where

R 0 inj = r + c k 0 p r R S -1 2 .
(B8)

1 DoelemanFIG. 2 .

 12 FIG.2. Visibility amplitude at 230 GHz as a function of baseline length. The blue solid line represents the synchrotron emission from a spike of 10 GeV DM annihilating into b b, with B = 10 G (solid), B = 10 2 G (dashed), and B = 10 3 G (dot-dashed), for a Schwarzschild BH (left panel) and a maximally rotating BH (right panel). The annihilation cross-sections correspond to the best fit to the EHT data from Refs.[13,14], given in Tables I and II for the Schwarzschild and maximally rotating cases respectively.

ing

  Eq. (B3) with initial condition p(R inj ) = p inj leads to p(r; R inj , p inj ) = p inj m e c) 2 .(B5)We consider the case α = 0, corresponding to a homogeneous magnetic field.The solution of the propagation equation in the ultrarelativistic regime is then given by f (r, p)

  ). A.E.B. and M.K. receive financial support from the Perimeter Institute for Theoretical Physics and the Natural Sciences and Engineering Research Council of Canada through a Discovery Grant. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. We acknowledge support by the Canadian Institute for Advanced Research for the participation of A.E.B. and J.S. at the annual meetings of the CIFAR cosmology program where this project was initiated.

		30		carried out in the ILP LABEX (ANR-10-LABX-63) and sup-
		20	DM spike, Schwarzschild BH DM spike, maximally rotating BH Akiyama et al. 2015	ported by French state funds managed by the ANR, within the Investissements d'Avenir programme (ANR-11-IDEX-0004-02
	Closure phase (deg)	10 0 10
		20
		t UTC (hr) 5 6 7 8 9 10 11 12 13 14 30
	FIG. 3. Closure phase as a function of universal time, for a
	Schwarzschild BH (solid) and a maximally rotating BH (dashed), for
	the current VLBI triangle between Arizona, California and Hawaii.
	Project No. 267117 (DARK) hosted by UPMC, and has been

Appendix A: Normalization of the DM profile

Note that M87 cannot be seen by the SPT.

The mass of the central BH is M BH = (6.4 ± 0.5) × 10 9 M [4], the corresponding Schwarzschild radius is R S = 6 × 10 -4 pc and the distance of M87 is d M87 ≈ 16 Mpc[5]. Throughout we adopt the higher stellar dynamical mass for M87; using the (3.5 ± 0.8) × 10 9 M value found by gas dynamical studies[6] would reduce the angular scales by roughly a factor of two througout.

We focus on the standard b quark channel for simplicity. Injection of electrons through a different channel would primarily result in a rescaling of the intensity at the frequency of interest, slightly changing the best fit crosssections we derive.

Moreover, the intensity profile for the no-spike case is much flatter than for a spiky profile, leading to a larger ratio of flux on long to short baselines, with no significant contribution on long baselines.

ACKNOWLEDGMENTS

This work has been supported by UPMC, CNRS and STFC. This research has also received support at IAP from the ERC We then use this value as an effective upper bound for the integral of Eq. (B6).

TABLE II. Best fit annihilation cross-section and reduced chi-squared χ 2 red = χ 2 /d.o.f. for various DM masses and magnetic field strengths, for the maximally rotating case. Values of the reduced chi-squared are given for illustration.

m DM = 10 GeV m DM = 10 2 GeV m DM = 10 3 GeV B = 10 G σv bf = 3.1 × 10 -31 cm 3 s -1 , χ 2 red = 6.5 σv bf = 1.2 × 10 -29 cm 3 s -1 , χ 2 red = 6.0 σv bf = 5.2 × 10 -28 cm 3 s -1 , χ 2 red = 5.8 B = 10 2 G σv bf = 2.9 × 10 -31 cm 3 s -1 , χ 2 red = 11 σv bf = 1.3 × 10 -29 cm 3 s -1 , χ 2 red = 11 σv bf = 5.6 × 10 -28 cm 3 s -1 , χ 2 red = 11 B = 10 3 G σv bf = 1.3 × 10 -30 cm 3 s -1 , χ 2 red = 12 σv bf = 5.6 × 10 -29 cm 3 s -1 , χ 2 red = 12 σv bf = 2.5 × 10 -27 cm 3 s -1 , χ 2 red = 12 law profile ρ 0 (r/r 0 ) -γ [19] at the center of the M87 galaxy:

where the saturation density determined by DM annihilations reads

with m DM and σv respectively the mass and annihilation cross-section of the DM particle, and we take t BH = 10 8 yr for the age of the BH. The spike profile reads

where

γ sp = (9 -2γ)/(4γ). We use M BH = 6.4 × 10 9 M for the mass of the BH [4], the corresponding Schwarzschild radius is R S = 6 × 10 -4 pc, and we take α γ = 0.1. We fix r 0 = 20 kpc for the halo (as for the Milky Way), and we must then determine the normalization ρ 0 . We choose ρ 0 in such a way that the profile is compatible with both the total mass of the galaxy and the mass enclosed within the radius of influence of the BH, of order 10 5 R S . We thus follow the procedure described in Ref. [33]: the DM mass within the region that is relevant for the determination of the BH mass, typically within R 0 = 10 5 R S , must be smaller than the uncertainty on the BH mass ∆M BH . ρ 0 is thus obtained by solving the following equation:

with ∆M BH = 5 × 10 8 M . Considering the complex dependence of ρ on ρ 0 , we use the fact that the mass is dominated by the contribution from r R S , i.e., typically r > R min = O(100R S ). In this regime we have ρ ∼ ρ sp (r). can also factorize the dependence on ρ 0 in

3-γ , and we finally obtain

We take γ = 1, which corresponds to the NFW profile. The corresponding spike has a power-law index of γ sp = 7/3. Numerically, we get ρ 0 ≈ 2.5 GeV cm -3 for γ = 1. Finally, the total mass within 50 kpc is ∼ 4 × 10 12 M , compatible with the value derived from observations, 6 × 10 12 M [34].

In practice, the saturation radius r sat -for which ρ sat = ρ(r sat )-is smaller than the Schwarzschild radius of the BH for all values of the DM mass and annihilation cross-section of interest here, so that the DM profile reads more simply

Note that one usually assumes that the DM profile vanishes below 4R S (or 2R S from the full relativistic calculation for a static BH [35]) due to DM particles captured by the BH.

Here, for simplicity, to study the potential of the EHT for probing very steep power-law density profiles, we consider a DM spike that goes all the way down to the horizon of the BH, i.e. R S for a Schwarzschild BH and R S /2 for a maximally rotating BH. However, this simplification has a negligible impact on our results.

Appendix B: Solving the cosmic-ray equation in the presence of an advection flow towards the central BH

The propagation equation of electrons and positrons in the presence of advection and synchrotron losses,