Nesrine Ben Hassine
email: nesrine.ben-hassine@inria.fr

Dana Marinca
email: dana.marinca@uvsq.fr

Pascale Minet
email: pascale.minet@inria.fr

Dominique Barth
email: dominique.barth@uvsq.fr

Machine Learning and Popularity Prediction of a Video Content

Keywords: CDN, popularity prediction, YouTube contens, phase profile I

In Content Delivery Networks (CDNs), the knowledge of the popularity of video contents helps the manager to take efficient decisions about which video contents should be cached near the end users and also about the duplication degree of each video to satisfy the end user Quality of Experience. This paper focuses on predicting the popularity of video contents, in terms of number of solicitations. For that purpose, different software entities, called experts compute popularity values of video contents. Each expert uses its own prediction method. The expert prediction accuracy is evaluated by a loss function as the discrepancy between the prediction value and the real number of solicitations. The simulations based on real YouTube traces show that the accuracy of prediction is improved by splitting the video content profile in contiguous phases. Different prediction methods are compared and also different phase change-points detection methods are evaluated in order to identify the method (or method parameters) minimizing the cumulated discrepancy compared to real solicitations of video contents.

INTRODUCTION

A. Context

Due to the increasing use of the Internet services, there has been an important growth in network traffic. This results in decreases in service quality. To overcome the limitation of the Internet, Content Delivery Networks (CDNs) provide a method for a more efficient delivery of contents to end users. They distribute contents geographically; heavily requested contents are cached at various locations, close to users who request them the most. Thereby, it is crucial to determine whether a content is requested sufficiently often so that it is worth caching it. The user request is redirected to the nearest location storing the solicited content. Moreover, cache management is a challenging task. Given its finite storage capacity and the dynamicity of the demand, each stored content has to be continuously monitored so that it is pushed out from the cache once it is no longer required. These issues require an in-depth analysis of the dynamics of contents requests.

In the Global Internet phenomena report 2014 [START_REF]Global Internet phenomena report[END_REF], streaming traffic remains the major source of Internet traffic in the world. Understanding the dynamics of the number of requests for streaming traffic is consequently both critical and indispensable for a better resource management. Streaming contents are mainly formed of video contents and they are basically originated from Video on Demand (VoD) platforms and User Generated Content (UGC) systems such as YouTube.

To better manage the enormous number of video contents, Youtube maintains various categories of video contents. The most popular video categories are cartoon, music, sports, movies, news... where some contents may belong to several categories. Hence, this classification helps people find the contents interesting them more easily.

The behavior of video contents over time reveals distinct patterns of popularity evolution. A huge number of contents are viewed only a few times. Others are largely required. There are contents that are popular for a long period of time whereas others are much solicited for few days or even hours. Hence, it is possible to estimate the future popularity of content by examining the types of popularity growth behavior that contents display over time. Therefore, demands can have different behaviors depending on the content and the time interval.

As a matter of fact, popularity evolution pattern is a series of phases. Each phase may be very different from the others, whereas it retains an homogeneous behavior inside. Let us define the concept of phase: Definition 1: A phase can be defined as an interval of time during which a measured metric remains relatively stable. Thus, if a phase is properly identified, the prediction within a phase should be more accurate. That is what we want to prove in this paper.

The paper is organized as follows. In Section I-B, we describe different methods to detect phase changes. In Section II, we give the problem statement and define the experts used to predict popularity of video contents. Results obtained from real traces are presented and explained in Section III, where we use different methods to detect phase changes. Finally, we conclude in Section IV and list some further work.

B. Related Work

Let us consider a time series of data characterizing a phenomenon monitored. Usually, it is possible to identify subsets of contiguous data where a same model can be applied to represent the phenomenon considered. It is then of the utmost importance to identify the time points at which changes occur in the data collected.

A large number of change-points detection methods has been proposed to tackle this problem. Most of the studies consider a data set consisting of observations at discrete times 𝑋 1 . . 𝑋 𝑛 spaced at uniform time intervals, let us assume we can partition the data into an unknown number 𝐾 of partitions, 𝑆 1 . . 𝑆 𝑘 . One of the most difficult issues is estimating the number of change points. Bai and Perron [START_REF] Bai1 | Computation and analysis of multiple structural change models[END_REF] propose a method with an a priori fixed number of change-point s. Let 𝑚 be the maximum number of change-point s. For each 𝑘 ≤ 𝑚, a change-point location is the one ensuring the smallest within-segment square sums. To select the best number of change-points, they pick the partition that minimizes the Bayesian Information Criterion (BIC). The criterion of selection could be SIC (Schwarz Information Criterion) or AIC (Akaike Information Criterion)

In [START_REF] Taylor | Change-point analysis: A powerful tool for detecting changes[END_REF], Taylor developed a change point detection method to detect changes in time series (CUmulative SUM). This approach is based on the mean-shift model written as 𝑋 𝑖 = 𝑋 � + 𝜀 𝑖 . 𝑋 � is the simple average and 𝜀 𝑖 is the residual term associated with the 𝑖 𝑡ℎ observation. Residuals are assumed to be independent and identically distributed with mean zero. The cumulative sum is not the cumulative sum of the values but the sum of differences between the values and the average. The change-point is detected through searching for the maximum of absolute CUSUM. Once a change-point is detected, the time series is splitted into two segments and the procedure is repeated.

Structural Change Model (SCM) is another mean-shift model based method [START_REF] Taylor | Change-point analysis: A powerful tool for detecting changes[END_REF]. Initially, they split time series into two segments at any location 𝑘. Let the sum of squared residuals be defined as the mean square error estimator.

𝑀𝑆𝑀(𝑘) = ∑ (𝑋 𝑖 -𝑋 � 1 ) 2 + 𝑘 𝑖=1 ∑ (𝑋 𝑖 -𝑋 � 2 ) 2 𝑛 𝑖=𝑘+1
Where

𝑋 � 1 = ∑ 𝑋 𝑖 𝑘 𝑖=1 𝑘 , 𝑋 � 2 = ∑ 𝑋 𝑖 𝑛 𝑖=𝑘+1 𝑛-𝑘
. The change-point is found by minimizing the sum of squares residuals and located at 𝑘 � = 𝑎𝑎𝑎𝑚𝑖𝑛 1≤𝑘≤𝑁 𝑀𝑆𝑀(𝑘).

As CUSUM and SCM, Circular Binary Segmentation (CBS) [5][6] needs the whole data to perform its analysis. Each time a change-point is detected, data are splitted into two segments at the change-point location. The method stops when no change-point can be found in any segment. The location of change-point is the one maximizing the absolute value of the tstatic, where t-static is a statistic test given by

𝑇 𝑖,𝑗 = �𝑌 � 𝑖,𝑗 -𝑍 � 𝑖,𝑗 � * 𝑠 𝑖𝑗 2
Where 𝑌 � 𝑖,𝑗 is the mean over �𝑋 𝑖+1 . . 𝑋 𝑗 �, 𝑍 ̅ 𝑖,𝑗 is the mean over �𝑋 1 . . 𝑋 𝑖+1 , 𝑋 𝑗+1 . . 𝑋 𝑛 � and 𝑠 𝑖𝑗 2 is the corresponding mean squared error.

Instead of estimating that a change-point has occurred, the Bayesian method computes the probability of a hypothesis. It expresses uncertainty about the number and location of change-point s. It computes the probability of a change-point at each location. In the Bayesian analysis, the posterior probability measures how likely we have a change-point at t. It is the weight given to hypothesis ℎ where ℎ = {𝑡ℎ𝑒𝑎𝑒 𝑖𝑠 𝑎 𝑐ℎ𝑎𝑛𝑎𝑒𝑎𝑎𝑖𝑛𝑡 𝑎𝑡 𝑡}. This probability is computed using the Bayes' theorem given by

𝑃(ℎ|𝑀) = 𝑃(𝑀|ℎ) * 𝑃(ℎ) 𝑃(𝑀)
where E is an evidence. Once we obtain the posterior probability, it is used to compute the posterior prediction distribution. This latter is a weighted average of the predictions of each individual hypothesis.

All the methods mentioned above are off-line and could be applied only once all the data are collected. In this paper, we compare an off-line method based on the Bayesian approach with two on-line methods that we will present in Section II-D.

II. THEORETICAL FRAMEWORK

A. Problem Statement

The goal of this paper is to predict the number of solicitations of individual video contents. We use real data sets of number of solicitations of video contents extracted from YouTube platform.

We adopt a machine learning approach and propose a model using several software entities called experts. Each expert computes and predicts future solicitation number for individual video contents using its own computation method (for example Basic, DES) with its own parameter set (e.g. size of its observation window, smoothing parameter).

We use a set of experts 𝑀, each expert 𝑀 𝑖 ∈ 𝑀 computes at time 𝑀 the value 𝑎 𝑖,𝑡+1 , representing the predicted number of solicitation of a given video, using its own prediction strategy. The goal of each expert is to predict a solicitation number as close as possible of the real solicitation number 𝑦 𝑡+1 . We will describe two prediction strategies in Section II-B. The accuracy of experts prediction is evaluated by the cumulative loss, as defined in Section II-C.

Focusing on the evolution of the number of solicitations of a video content, called video profile, we observe the existence of phases. A phase is characterized by a similar evolution of the number of solicitations as shows the video content profile depicted in Figure 1. The video contents considered belong to the Music and Cartoon categories. Figures on the right represent the number of solicitations as a function of time. For each profile, we represent also the cumulated number of solicitations (Figure 1). Tics represent phase changes.

The goal of this study is to answer the following questions:

• Can we improve the prediction accuracy through the cooperation of several experts for a given video content instead of one expert prediction?

• Does the use of phases improve the accuracy of predictions?

• Has the phases detection method an impact on the quality of predictions?

This work was supported by ANR project NetLearn -ANR-13-INFR-004

B. Experts

We recall that each expert prediction concerns an individual video content. At a given time 𝑡, each expert 𝑀 𝑖 ∈ 𝑀 computes and predicts the value 𝑎 𝑖,𝑡+1 , which will be compared to the real value at the time 𝑡 + 1. We define hereafter two concepts: (a) the observation window 𝑎𝑜, a past time interval ending at time 𝑡, containing past consecutive values of the number of solicitations for the considered video and (b) the prediction window 𝑎𝑜, starting at time 𝑡, representing the future time interval at the end of which the accuracy of the prediction 𝑎 𝑖,𝑡+1 will be evaluated.

1) DES expert:

DES experts, Double Exponential Smoothing, [START_REF] Bianchides | Prediction, Learning, and Games[END_REF] have been introduced to cope with trends (e.g. increase) in the data to predict. A DES expert applies two exponential smoothing:

• first on the value 𝑦 𝑡 observed at time 𝑡 when it computes 𝑆 𝑖,𝑡 , . We have: 𝑆 𝑖,𝑡 ′ = 𝛼𝑦 𝑡 + (1 -𝛼)𝑆 𝑖,𝑡-1 ′ , where 𝛼 is the smoothing factor, with 0 < 𝛼 < 1.

• second on the value of 𝑆 𝑖,𝑡 ′ when it computes 𝑆 𝑖,𝑡 ′′ , using the same smoothing factor 𝛼. We have: 𝑆 𝑖,𝑡 ′′ = 𝛼𝑆 𝑖,𝑡 ′ + (1 -𝛼)𝑆 𝑖,𝑡-1

′′

, with 𝑆 𝑖,0 ′ = 𝑆 𝑖,0 ′′ = 𝑦 0 .

The value 𝑎 𝑖,𝑡+1 predicted by any Double Exponential Smoothing (DES) expert 𝑀 𝑖 at time 𝑡 + 1 is given by: 𝑎 𝑖,𝑡+1 = 𝐿 𝑖,𝑡 + 𝑇 𝑖,𝑡 , where 𝐿 𝑖,𝑡 = 2𝑆 𝑖,𝑡 ′ -𝑆 𝑖,𝑡 ′′ denotes the estimated level and 𝑇 𝑖,𝑡 = 𝛼 1-𝛼 (𝑆 𝑖,𝑡 ′ -𝑆 𝑖,𝑡 ′′ ) the estimated trend.

A DES expert has two parameters: 𝛼 the smoothing factor and ow the size of the observation window, which can be different for different experts.

2) Basic expert:

The basic expert is the expert that at time 𝑡 + 1 predicts 𝑎 𝑖,𝑡+1 = 𝑦 𝑡 + (𝑦 𝑡 -𝑦 𝑡-1 ) = 2𝑦 𝑡 -𝑦 𝑡-1 . In other words, it predicts the last known value plus its increment with the previous value. Such an expert is very simple to implement

C. Best expert

Informally, the best expert is the one providing the most accurate predictions. This accuracy can be evaluated by the loss incurred by this expert.

Definition 2: For any expert 𝑀 𝑖 , we define its instantaneous loss at time 𝑡 as 𝑙 𝑖,𝑡 = �𝑦 𝑡 -𝑎 𝑖,𝑡 � and its cumulated loss at time 𝑡 as 𝐿 𝑖,𝑡 = ∑ 𝑙 𝑖,𝑘 𝑡 𝑘=0 . Definition 3: For any set of experts 𝑀, we define the best expert as the expert ∈ 𝑀 that minimizes the cumulated loss. Notice that the expert that minimizes the cumulated loss also minimizes the average instantaneous loss. We now propose to extend the concept of best expert to the best expert on a given phase as follows:

Definition 4: For any given phase, for any set of experts 𝑀, we define the best expert for this phase as the expert ∈ 𝑀 that minimizes the cumulated loss on the phase.

D. Phase Detection

In this paper we distinguish several methods to detect phase changes.

• Detection by a sudden increase of the instantaneous loss of the Basic expert. This detection is simple to implement and can easily be done on-line.

• Detection by applying Bayesian inference method. A phase change is detected by applying Bayesian inference method, implemented in R tool. This detection requires the knowledge of the whole profile and consequently is done off-line only. We use the module implemented in the R tool.

• Detection of a contiguous period. Each phase has a constant duration of k time units, with k an integer greater than 1. This detection is very simple but may not correspond to a change in the profile of solicitations. Such detection is well adapted to profiles presenting periodic behaviors.

This paper studies the impact of phase detection method on the prediction accuracy.

E. The Prediction Method

We consider several video profiles belonging to different content categories. In an off-line analysis we run a lot of DES experts analyzing the considered profiles and predicting the number of solicitations. Among the identified experts we determine a subset of DES experts, having a cumulated loss less than 102% of the cumulated loss of the best DES expert on the entire profiles. These experts are characterized by particular values of the parameters 𝛼 and 𝑎𝑜. This simulation helps us to identify a short list of DES experts. Furthermore, we consider the profile segmentation according to the different methods specified in Section II-D. For each change-point detection method, the cumulated loss of each expert is evaluated at the end of each phase. The best expert on each phase is the expert having the smallest cumulated loss. At the end of the profile, the prediction cumulated loss represents the sum over all phases of the cumulated loss of the best expert on each phase.

III. SIMULATION RESULTS

A. Determination of the Best DES Experts

We considered 50 video content solicitation profiles collected from YouTube plate-forme, belonging to different video categories (society, music, films, children, etc.). To determine the best values of the smoothing factor 𝛼 and the observation window 𝑎𝑜 for the DES experts, we run many simulations on collected profiles. We succeeded to identify a list of six DES experts, the best ones, which are given in Table I. The parameter 𝑎𝑜𝑚𝑎𝑜 corresponds to the duration of the considered phase. The next simulations are based only on the selected subset of the best experts presented in the previous Table.

B. Phase Change Detection

We now focus on the detection of phase changes. We apply the three methods described in Section II-D. For the detection of a contiguous period, we use 7 days, 14 days and 21 days.

C. Prediction Evaluation 1) No phase detection:

In the first series of experiments, we consider video contents profiles without distinguishing any phase. We test the Basic expert described in Section II.B.2 and the best DES experts given in Table I. Results obtained on real traces show that for different profiles the associated best experts could be different. In Figure 2, the video content considered belongs to the category Music. For this content, the expert minimizing the cumulated loss is the DES expert with 𝛼 = 0.9999 whatever the size of the observation window. The Basic expert has close results. We now focus on a video content of category Cartoon. For this content, the Basic expert drastically outperforms the DES experts, as depicted in Figure 3. Notice that all the DES Brown experts have approximately the same cumulated loss. Thereby, the Basic expert has very good results in both cases.

2) Detection of Phase Changes:

In the second series of experiments, we want to answer the following question: can we improve the accuracy of prediction by detecting phase changes? Tested video profiles are segmented using the methods given in Section II-D. Each method has its detected points which differ from the output of other methods. Results are depicted in Figure 4 for a cartoon video content. We will compare the cumulated loss obtained by these methods with the cumulated loss obtained when ignoring the profiles phases, denoted as "No change-point detection". More precisely, two cumulated losses are compared over each profile:

(1) Using "No change-point detection" method, we compute the cumulated loss of the best expert at the end of the profile;

(2) Using each change-point detection method, we compute the sum over all profile phases of the cumulated loss of the best expert per phase. At the end of each phase, the cumulated loss of each expert is reinitialized and another best expert can replace the previous one. As depicted in Figures 5 and6, we observe that regardless of the method used, all change-point detection methods reduce the cumulated loss obtained without detecting phases. This can be explained by the expert changes at the end of each phase. These simulation results show that the change-point detection has a strong impact on the accuracy of the prediction. As shown in Figures 5 and6, Basic based change-point detection provides the largest cumulated loss compared to other change-point detection methods. Nevertheless, in both figures, the smallest cumulated loss is obtained when the change-point detection method uses a period of 7 consecutive days. Earlier studies [START_REF] Hongliang | Understanding user behavior in large-scale video-on-demand systems[END_REF] have shown that users' access follows a pattern over a week. This explains the results obtained when it is assumed that there exit a change every week. Based on this assumption, the cumulated loss could be further reduced if the end of the interval coincides with the weekend.

Thus, the best method to detect phase changes is the method considering fixed periods of 7 consecutive days. These results justify an approach based on detecting phase changes in video profiles evolution.

3) Behavior of experts near a phase change:

In the third series of experiments, we evaluate the behavior of experts near phase changes. For instance, in Figure 7, we compare the behavior of Basic and DES experts with and without change-point detection. We assume that we detect a change-point every period of 7 consecutive days. The prediction of experts with change-point detection deviates slightly from the real value. However, experts without changepoint detection make large errors at the phase change, especially the Basic expert. We also observe that within a phase all selected experts provide predictions close to the real value.

IV. CONCLUSION

Performances of Content Delivery Networks are usually measured by the hit ratio of user requests. The main goal of the cache management policy is to maximize the hit ratio over CDN. In this paper, we propose to predict the popularity of video contents, in terms of the number of solicitations. This popularity can be used to decide which contents to be cached. The best expert predicting the number of solicitations is identified as the expert minimizing the cumulated loss over the profile. This expert provides the prediction the closest to the real value. To improve the accuracy of predictions, we defined different phases over each profile. The best expert per phase outperforms the best expert on the whole video content profile. However a false detection of a phase change (false positive detection) may decrease the quality of the prediction.

  a) Music content: original profile. b) Music content: cumulative profile. c) Cartoon content: original profile. d) Cartoon content: cumulative profile.

Fig. 1 .

 1 Fig. 1. Phases in the number of solicitations of a video content.

Fig. 2 .

 2 Fig. 2. The total loss of selected experts over the whole music contents

Fig. 3 .

 3 Fig.3. The total loss of selected experts over the whole cartoon contents.

Fig. 4 .

 4 Fig. 4. Cartoon content: Detected change-points using different methods.

Fig. 5 .

 5 Fig. 5. Music content: The minimum cumulated loss with no change-point detection vs the cumulated loss with different change-point detection methods.

Fig. 6 .

 6 Fig. 6. Cartoon content: The minimum cumulated loss with no change-point detection vs the cumulated loss with different change-point detection methods.

Fig. 7 .

 7 Fig. 7. Cartoon content: Prediction based on change-point detection using fixed period of 7 days.

TABLE I .

 I BEST DES EXPERTS.

	𝛼	𝑎𝑜𝑚𝑎𝑜
	0.8	𝑎𝑜 = 7 ou 𝑎𝑜𝑚𝑎𝑜
	0.93	𝑎𝑜 = 7 ou 𝑎𝑜𝑚𝑎𝑜
	0.9999	𝑎𝑜 = 7 ou 𝑎𝑜𝑚𝑎𝑜


Fabio Martínez

Antoine Manzanera

Eduardo Romero
email: edromero@unal.edu.co

Automatic Analysis and Characterization of the hummingbird wings motion using dense optical flow features

Keywords: wing hummingbird flight, dense optical flow, hummingbird yaw turn analysis, Kinematic hummingbird patterns

come

Automatic analysis and characterization of the hummingbird wings motion using dense optical flow features

Introduction

Emulation of flying animal movements may facilitate the design of smarter micro air vehicles (MAVs), allowing them to perform more complex maneuvers and expend less energy [START_REF] De Margerie | Artificial evolution of the morphology and kinematics in a flapping-wing mini uav Bioinspir[END_REF]. Among these MAVs, hovering machines have succeeded in executing elaborated movements within very limited spaces [START_REF] De Margerie | Artificial evolution of the morphology and kinematics in a flapping-wing mini uav Bioinspir[END_REF]. Designs of such new hovering MAVs may be inspired in animals such as the Hummingbird [START_REF] Altshuler | The ecological and evolutionary interface of hummingbird flight physiology[END_REF] [START_REF] Keennon | Development of the nano hummingbird: a tailless flapping wing micro air vehicle 50th[END_REF], a bird species able to hover in mid-air and flight backwards [START_REF] Sapir | Backward flight in hummingbirds employs unique kinematic adjustments and entails low metabolic cost[END_REF]. Their maneuverability is based on their ability to easily turn in very restricted spaces, a skill developed to access flowers with sugar rich nectars textcolorMyColor and also to catch insects. A proper design of flapping-wing flight machines requires a precise comprehension of the many interactions between kinematic and morphological wing parameters. Temporal wing deformation patterns constitute the fundamental variables to estimate the nonlinear relationships between the different wing segments [START_REF] Shyy | Recent progress in flapping wing aerodynamics and aeroelasticity[END_REF]. Moreover, rotations and velocities are related not only with kinematic information but also with forces and torques produced by muscles and tendons [START_REF] Douglas | Neuromuscular control of wingbeat kinematics in Anna's hummingbirds (Calypte anna)[END_REF], [START_REF] De Margerie | Artificial evolution of the morphology and kinematics in a flapping-wing mini uav Bioinspir[END_REF]. Nevertheless, in spite of extensive studies, many of the morphological, local dynamic, kinematic and even aerodynamic details are poorly understood [START_REF] De Margerie | Artificial evolution of the morphology and kinematics in a flapping-wing mini uav Bioinspir[END_REF] because of the still limited tools used to capture flight patterns, like for instance flapping frequency, stroke amplitude, angular velocity and even spatio-temporal wing deformation.

Many questions about the flight mechanisms still remain open, e.g., how the hummingbird accomplishes yaw turns during the hovering flight or how the bird is able to keep turning once he starts this maneuver. Altshuler et al. [START_REF] Altshuler | Middleton Wingbeat kinematics and motor control of yaw turns in Anna's hummingbirds (Calypte anna)[END_REF], have recently published a study about the wingbeat kinematics and the motor control of yaw turns. This investigation correlated muscle EMG signals with kinematic trajectories and found differences among the different stroke flight phases and even asymmetry between the wings. In many cases however, the only available information about these turns comes from other species like the gliding or flight patterns observed in insects [START_REF] Hedrick | Low speed maneuvering flight of the rose-breasted cockatoo (eolophus roseicapillus) i. kinematic and neuromuscular control of turning[END_REF], [START_REF] Hedrick | Low speed maneuvering flight of the rose-breasted cockatoo (eolophus roseicapillus). ii. inertial and aerodynamic reorientation[END_REF], [START_REF] Warrick | Asymmetrical force production in the maneuvering flight of pigeons[END_REF], [START_REF] Warrick | Bird maneuvering flight: blurred bodies, clear heads[END_REF], [START_REF] Fry | The aerodynamics of free-flight maneuvers in drosophila[END_REF], [START_REF] Dawson | Forewing asymmetries during auditory avoidance in flying locusts The[END_REF]. However these animals are actually very different since they are entailed with simple wings that operate over a larger range of speeds [START_REF] Warrick | The aerodynamics of hummingbird flight[END_REF] and are supported by an outer skeleton. In addition, while birds roll their body to orient the wing forces, insects can easily turn by asymmetric left-right wingbeats [START_REF] Altshuler | Middleton Wingbeat kinematics and motor control of yaw turns in Anna's hummingbirds (Calypte anna)[END_REF], [START_REF] Hedrick | X Wingbeat time and the scaling of passive rotational damping in flapping flight[END_REF], [START_REF] Ristroph | Automated hull reconstruction motion tracking (HRMT) applied to sideways maneuvers of free-flying insects[END_REF], [START_REF] Fry | Dickinson The aerodynamics of free-flight maneuvers in Drosophila[END_REF] . The hummingbird flight has been broadly studied using customized analysis systems that require a set of reference points, somehow highlighted and superimposed on the animal structure. These salient anatomical points are usually drawn by an expert [START_REF] Ortega-Jimenez | Source Aerial shaking performance of wet Annaś hummingbirds[END_REF], adding new variables to the analysis since such delineation is a frame-per-frame task which inevitably leads to high inter and intraobserver variabilities. Smoothing and filtering strategies are often used to remove errors introduced during the digitazing process [START_REF] Walker | Estimating velocities and accelerations of animal locomotion: a simulation experiment comparing numerical differentiation algorithms[END_REF]. However, these strategies are sensible to the parameters of starting and require a fully knowledge about the signal nature, for instance, a higher sample rate allows the trade off of the signal. Any modification of the signal can lead to an alteration of the observational data. Furthermore, these tools introduce artifacts that alter the natural hummingbird motion gestures since they need markers that must be placed on the wing bird, markers that in addition can be easily occluded during a wing cycle [START_REF] Walker | Taylor Deformable wing kinematics in freeflying hoverflies[END_REF].Currently, a pyramidal Lukas-Kanade optical flow strategy have been proposed to estimate the deformation of birlike airfoil showing coherent results by using artificial structures. This strategy is however limited to follow large displacements since the optical flow strategy used is based on a local regularization, thereby assuming that the object has small displacement among consecutive frames [START_REF] Bansmer | P Experimental and numerical fluid-structure analysis of rigid and flexible flapping airfoils[END_REF].. Classical approaches may miss many important details of the real dynamic patterns since they rely completely in both an expert user and an invasive tracking strategy. Recently, some approximations have automatically estimated the dynamics of flight trajectories, most of them using machine-learning methods that require a considerable number of videos for training and finding global flight patterns [START_REF] Straw | Multi-camera Realtime 3D Tracking of Multiple Flying Animals[END_REF]. Overall, these kinematic methods model the wing as a stiff airfoil [START_REF] Wu | Computational aerodynamics of low Reynolds number plunging, pitching and flexible wings for MAV applications Bioinspiration and Biomimetics[END_REF], but in reality the wing is a very flexible structure, able to asymmetrically produce non uniform movements [START_REF] Wu | Computational aerodynamics of low Reynolds number plunging, pitching and flexible wings for MAV applications Bioinspiration and Biomimetics[END_REF]. Some works have extracted patterns to characterize the hovering hummingbird dynamics [START_REF] Warrick | Aerodynamics of the hovering hummingbird[END_REF], supporting the importance of performing analyses of the whole structure rather than following some isolated reference points. A very useful tool has been the Particle image velocimetry (PIV) [START_REF] Bomphrey | Advances in Animal Flight Aerodynamics Through Flow Measurement[END_REF], [START_REF] Adrian | Particle-Imaging Techniques for Experimental Fluid Mechanics Annual review of fluid mechanics[END_REF], a strategy that permits to estimate aerodynamic variables by indirect measurement of the resultant flow produced in the air by the wing motion. In practice the set of particles, near the wing boundary, is supposed to closely follow the animal motion. The whole procedure consists then in illuminating such particles with a pulsed light laser and capturing this motion with a camera. A flow field is obtained by cross-correlating consecutive image patterns formed by particles highlighted by the laser [START_REF] Bomphrey | Advances in Animal Flight Aerodynamics Through Flow Measurement[END_REF], [START_REF] Adrian | Particle-Imaging Techniques for Experimental Fluid Mechanics Annual review of fluid mechanics[END_REF]. A main advantage of this method is its relative non-invasiveness, but a limitation is that the particle illumination may produce interference phenomena near the wing boundary. Besides, the whole dynamics can not be estimated with the PIV method since the captured flow is related to the external motion and hence the analysis of forces and torques is only indirect. In other words, this method may improve understanding of the hummingbird movement from an aerodynamic standpoint, but not from a kinematic point of view, whose importance is mostly associated to the local wing deformation and temporal rotation. Lately, several researchers have introduced mechanical and mathematical wing models, aiming to exploit some basic physical relationships to simulate different types of wing motions [START_REF] Bos | Influence of wing kinematics on aerodynamic performance in hovering insect flight[END_REF], [START_REF] Liang | Computational Analysis of Hovering Hummingbird Flight In 48th[END_REF], [START_REF] Wang | Unsteady forces and flows in low Reynolds number hovering flight: two-dimensional computations vs robotic wing experiments[END_REF], [START_REF] Altshuler | Aerodynamic forces of revolving hummingbird wings and wing models[END_REF]. However, these models are very restrictive, for instance wings are represented as rigid structures and therefore only very linear approximations are possible, whereby the obtained conclusions are limited and the predictive model capabilities are largely reduced [START_REF] Shyy | Computational aerodynamics of low Reynolds number plunging, pitching and flexible wings for MAV applications[END_REF].

In this work, a novel description of the wing dynamics is introduced, providing previously unavailable information about rotational and deformation wing patterns. This spatiotemporal description is based on a set of velocity vectors obtained from a dense optical flow, computed from sequences captured with a 2D camera which is located perpendicular to the main axis of the bird. Four principal temporal estimations are computed: (1) the global angular acceleration of the hummingbird, (2) the angular acceleration profiles computed over different wing segments, (3) the variance of the spatial orientation of the velocity of the optical flow field for each of the two wings and (4) a distribution map of the local orientation variation within the wing. These four complementary measures allow not only to follow up the wing but also to determine the instants of maximal and minimal deformations, as well as the regions with larger deformation. This paper is organized as follows: Section 2 introduces the proposed method, section 2.4 describes the dataset, section 3 demonstrates the effectiveness of the method and section 4 presents a discussion of the proposed work. The last section concludes and presents some possible future works.

Materials and Methods

Overall, a dense optical flow provides a motion description by estimating the instantaneous displacement of each pixel between consecutive frames. Such vector of velocity optical flow allows to compute four complementary motion descriptors at each time step, all of them further described hereafter.

The Dense Optical Flow Computation

First, a dense optical flow is estimated by matching all the points from frame to frame, using a nearest neighbor search in a projected feature space [START_REF] Manzanera | Local Jet Feature Space Framework for Image Processing and Representation[END_REF]. The multidimensional feature space is herein obtained using the local jets [START_REF] Florack | The Gaussian scale-space paradigm and the multiscale local jet[END_REF], an approximation to the image spatial structure provided by derivatives at different orders and scales:

f σ ij = f * ∂ i+j Gσ ∂x i 1 ∂x j 2
, where G σ is the 2d Gaussian function with standard deviation σ. Each pixel is associated then to a feature vector defined as the collection f σ ij ; i + j ≤ r, σ ∈ S , where r is the derivation order and S = {σ 1 , . . . σ q } the selected scales. The set of features, associated to a frame, is represented in a kd-tree structure which facilitates efficient distance search [START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF]. Finally V t (x), the apparent velocity vector for pixel x at time t is calculated by coupling the feature vector associated to x to its nearest neighbor in the feature space calculated at time t + 1. Thanks to the multiscale representation, which implies the matching of that extended spatial structure, this method provides a spatially consistent dense optical flow, without using explicit spatial regularization as required in other methods [START_REF] Aubert | Computing optical flow via variational techniques[END_REF], [START_REF] Berthold | Determining optical flow[END_REF].

Hummingbird motion patterns

The quantified optical flow is not altered by the hummingbird spinning (at ∼ 30 : Hz), basically because this movement is parallel to the camera axis and therefore the optical flow is unable to capture, i.e., only those particles moving in the plane perpendicular to the camera can be tracked. Quantifying kinematic flight patterns and their structural relations is crucial for optimal design and development of flapping-wing flight prototypes. These complex relationships may be primarily determined from the analysis of variables like the angular acceleration and temporal wing deformation, herein estimated from the dense optical flow. Such variables are temporal functions, some of them being global scalar functions of t and other local spatial maps indexed by t. Local dynamic patterns are determined by firstly segmenting the calculated apparent velocity flow, finding the two largest connected regions that correspond to the right and left wing and denoted V rt and V lt . The isolated remaining flow areas are assigned to the nearest region, according to the spatial Euclidean distance.

Angular Acceleration patterns

Rotational wing patterns are crucial to understand complex flight mechanisms, like: an optimal angle of attack, a large downstroke span ratio and a high wingbeat frequency that together maintain the achieved equilibrium and facilitate complex maneuvers [START_REF] Dickinson | Wing Rotation and the Aerodynamic Basis of Insect Flight Science[END_REF], [START_REF] Tobalske | Biomechanics of bird flight[END_REF]. In this work, the angular acceleration was estimated by associating the orientation of the computed apparent velocity optical flow and the wing, under the assumption that the rotation axis is perpendicular to the optical axis of the camera [START_REF] Verri | Torre Mathematical properties of the two-dimensional motion field: from singular points to motion parameters[END_REF], [START_REF] Higgins | The Interpretation of a Moving Retinal Image Proceedings of the Royal Society of London[END_REF]. Two configurations were herein proposed to analyze the hummingbird wing motion from an angular acceleration standpoint: a global angular acceleration of both wings and a temporal map of regional angular accelerations, constructed by estimating the angular acceleration profiles for a set of wing segments which correspond to a wing anatomical partition.

The Global Angular acceleration for each wing: For a non-zero flow vector V, let φ(V) be the apparent velocity flow orientation, the wing flows V r t and V l t were followed, corresponding to the right and left wings, which are bounded within a bounding-box at a particular time t and oriented to the main direction of each wing.Then, an averaged spatial orientation for each wing φ(V r,l ) t is calculated as the linear combination of every orientation, weighted by its norm. Under the assumption that there exists a dominant orientation of the velocity field, the rate of change of the angular velocity is then computed as Ω t = dφ(V r,l )t dt . This variable describes the turn angle during the flight, and highlights the relaxation and effort phases.

Angular acceleration spatial profiles: An angular acceleration profile was herein obtained

Hence, left and right temporal series are constructed by associating a wing region to the closer region in the precedent frame, in terms of the minimal distance d r,l t+1 between the centroids of consecutive wing regions (x, y) s r,l t , as follows:

d r,l t+1 = (x, y) s r,l t -(xy) s r,l t+1 2 2
. Then, a local angular estimation was herein obtained from the temporal acceleration profile of s t r,l by dividing the whole region into a set of N segments perpendicular to the main direction, whereby s t r,l ∈ s {r,l} 0 , s {r,l} 1 , . . . , s {r,l} i , . . . , s {r,l} N -1 t (Figure 1 (d)). Finally, a local angular acceleration for the i -th segment s i t was computed and its temporal map Ω(s i ) was obtained. Such representation is:

Ω(s i (x, t)) = dφ(V t (x))) dt ; x : x ∈ s i (1)

Wing deformation patterns

Although hummingbird kinematic patterns are generally similar to the insect flight, the morphological structure of the wings are very different given the mass and structural relationship, furthermore if considering that the wing muscle mass is largely determined by the need of producing efficient mechanical power. This is why the impact of the wing deformation during the flight must be quantified for understanding how the wing geometry supports forces and produces motion. For instance, it has been documented [START_REF] Warrick | Aerodynamics of the hovering hummingbird[END_REF] that the hummingbird upstroke and downstroke are not symmetrical. Such observation is strongly related to the local loads supported along the wing.

Wing deformation is herein estimated by computing the variance of the local velocity along the time. This estimation is carried out by assuming that every flow vector should be equally distributed along a rigid wing, while deformable structures produce a non homogeneous distribution, whereby the local motion variance is proportional to the physical deformation.

Local map of the wing deformation: The wing flow and its averaged orientation are used to compute a per-frame deformation map, estimated as a squared difference between each non-zero flow vector and the average orientation:

Q l,r t (x) = φ(V t (x)) -O l,r t 2 2 (2)
This analysis provides a stress map of the wing, highlighting the forces shearing the wing 1 .

Global wing deformation: For each wing, a global deformation W l,r t is estimated ‡ by averaging the local wing deformation map Q l,r t (x) . This measure allows to detect the largest deformation instants during a flight stroke.

Dataset description

Four annas hummingbirds (Calypte Anna's) were recorded. Each hummingbird was placed inside a flight chamber that contained a wooden perch in the corner and a feeder, in which a 1 ml syringe was mounted at the end of a metallic arm (0.22m × 0.19m × 0.15m). The feeder arm was connected to a stepper motor (MDrive 23 Plus, Intelligent Motion Systems, Inc.) placed at the centre of the cage roof. This setup allowed the bird to feed from the spinning feeder, while maintaining its centre of gravity in the same spot, executing a pure yaw turn. A white cardboard square was placed above the hummingbird to offer a contrasting background to the hummingbird body. The hummingbirds were trained to follow the feeder at 30 r.p.m. A stationary hovering flight was recorded using a high speed camera (Fastec Imaging, Troubleshooter) with spatial resolution of (640 × 480) and a temporal resolution of 1000 f ps. In this study, a set of videos recorded the hummingbird flight from a ventral view perspective (see Figure 1 ). A total of 91 flight cycles were herein considered. For each wing bird, an average of X wingbeats were computed independently. A complete description of the data herein used were already reported by Altshuler et al. who also captured muscle activation patterns with an electromiography device [START_REF] Altshuler | Middleton Wingbeat kinematics and motor control of yaw turns in Anna's hummingbirds (Calypte anna)[END_REF]. The electrical muscle activity was not herein considered since it was out of the scope of this paper, devoted to introduce a fully automatic method for analysis of the kinetic patterns.

Additionally, a controlled experiment was carried out by estimating the angular velocity of a rigid plate perpendicularly coupled to a servo-motor and recorded with a high speed ‡ where l stands for the left and r for the right camera, as shown in Figure X. The rotating mechanism used an encoder to control the angular velocity with a temporal resolution of 4096 pulses per revolution. This mechanism was captured with a high-speed camera, 102 f ps, and a spatial resolution of 640 × 480.

Statistical Analyses:

Statistical Analyses: The downstroke and upstroke were separated and a transition phase was defined as the interval between 40 % and 60 % of the normalized stroke. The t-student approach was used for pair comparisons between different variables, utilizing a mean value for each bird and for each flight phase. A first comparison was carried out between the angular acceleration computed with the proposed approach and manual delinations of the wing. Likewise, the t-student was used for comparison between the right and left wings of the fourt birds.

Results

The set of motion measures, herein proposed to analyze the hummingbird yaw turn flight, provides a valuable spatio-temporal wing description from local and global perspectives, highlighting the maximal and minimal moments of the angular acceleration and the wing deformation. The introduced analysis framework, allows a direct wing kinematic analysis, promoted by a dense optical flow that facilitates a spatio-temporal characterization of a particular wing region. In other words, this method tracks the whole wing structure and permits continuous local or regional analyses of the wing, much more representative of the wing dynamics than the few classical points of interest marked by an expert. An accurate estimation of the velocity with an optical flow field would facilitate the extraction of any motion pattern. However, a proper understanding the dynamics along the flight requires a balance between global and local information, a task for which a conventional optical flow is insufficient because it performs pixel matching at a very local level, while the quality of the global information depends on the spatial regularizer. In contrast, the multidimensional dense flow herein implemented, is more appropriate since matching is performed at regional level and no explicit spatial regularization is needed. A typical multidimensional dense flow estimation of the hummingbird flight is shown in Figure 6(b), where the apparent velocity vectors are observed to be mainly within the wings. However these multi scale optical flow advantages, the matching in large scales ( computed with large σ) may present some spatial inconsistencies regarding the location of the resultant vectors, for instance the vectors surrounding the wings, yet these vectors do represent actual dynamic patterns. For the experiments presented hereafter, the dense optical flow was parameterized as follows: 7 scales with a σ 0 = 0.3 and σ n+1 = 2σ n , while up to the second order derivatives were used, resulting in a descriptor vector of dimension 42.

For the controlled experiment, the rotating plate mechanism was triggered for different angular accelerations, ranging from 0.1 to 1 rad f rame 2 , according to the encoder sensor. Once the mechanism started from a pre-defined value, the camera recorded a sequence and a dense optical flow was computed. Figure 3 shows the dense optical flow calculated from a typical sequence of the rigid plate. The dense optical flow without any spatial restriction tracked the plate, clustering the regions with motions. Under a RMSE metric, the proposed approach achieved an averaged accuracy of 87.2 w.r.t to the mechanism sensor. Figure 2c shows the temporal deformation map Q t (x) as a function of time and space (x, t) obtained during a stroke cycle. The color scale represents the magnitude of the local deformation along the wing. Such deformation maps, estimated from the velocity fields, are totally coherent with the classical bio-mechanical descriptions of the musculoskeletal hummingbird structure: during the first stroke phase (corresponding to the [0 -20%] of the temporal serie) [START_REF] Douglas | Neuromuscular control of wingbeat kinematics in Anna's hummingbirds (Calypte anna)[END_REF], [START_REF] Altshuler | Middleton Wingbeat kinematics and motor control of yaw turns in Anna's hummingbirds (Calypte anna)[END_REF]. This action produces an initial wing deformation, basically controlled by the proximal wing, as illustrated in the first panel of the deformation map(Figure 2

(c)). During the intermediate stroke phase, the wing deformation is propagated to the distal segments, devoid of any skeletal structure, and whose maximal deformation occurs at about 60% of the stroke, as observed in the mid panel of the Figure 2(c). The final cycle is completed by backwardly extending the wings, as illustrated in the last panel of figure 2(c).

The estimated global angular acceleration Ω t correlates with the wing deformation and was calculated as the derivative of the averaged orientation of the apparent velocity flow for each wing. Every flight stroke, computed from each video, was normalized to be between 0 % and 100 % and then manually selected and phase-aligned by an expert. The pattern of a cycle was estimated by determining the correlation in the experimental group: for the whole set of 91 recorded sequences and for each instant t, a statistical boxplot was computed, within which the red line represents the median of the angular acceleration, the blue lower and upper boxes stand for the spread of data (interquartile range), the maximum and minimum values are the vertical dotted lines and the red asterisks stand for some outliers of the remaining data, as illustrated in Figure 3(a) and 3(b) for the two wings respectively. Of course, the sources of error are variable and several measurements were performed to highlight the origin of such variation, namely how much is due to variations among flights, between right and left wings and how much is caused by the new optical flow. The latter question was addressed by comparing the angular acceleration § obtained either from the optical flow or from a set of manual delineations of the wing. A t-student test showed no significant differences (p < 0.05) when comparing the manual and automatic estimations for the four test subjects, 91 flights and the two wings. Regarding the other part of the variability, the temporal normalized stroke follows a very § the angular acceleration is computed by deriving twice the trajectoty of the Center of Mass during the stroke.

similar pattern for each of the recorded wings, with minimal variance in the passive phases and larger differences for the active parts of the cycle. On the one hand, the mid part of the temporal normalized stroke shows a larger variance, the transition between the upstroke and downstroke phases. The observed acceleration peak, between the 40 % and 60 % percent of the normalized stroke, represents two complementary actions: an initial relaxation, when the wings move forwardly and join, followed by the muscle contraction required to separate them. During the transition between these two actions, the mean and standard deviation of the angular acceleration were computed for the right, 0.058 ± 0.469 rad f rame 2 , and left, -0.79 ± 0.526 rad f rame 2 , wings. A conventional t-student test showed no significant differences for α = 0.05. However when this analysis was performed individually, that is to say left and right wings of the same bird, significant differences appear (t-student test, p < 0.05) as illustrated in table 1.

On the other hand, the upstroke and downstroke phases correspond to the [10 % -40 %] and [60%-90%] of the temporal normalized stroke, respectively. In such intervals, the flight is basically passive and smooth because motion is mainly driven by the inertia. For the upstroke phase, the mean and standard deviations for the right and left wings were -0.03±0.291 rad f rame 2 and 0.143 ± 0.247 rad f rame 2 respectively, while the mean and standard deviations were for the right, -0.025±0.249 rad f rame 2 , and left, -0.015±0.298 rad f rame 2 , wings during the the downstroke, differences that were found to be not significant (t-student test and p < 0.05).

Table 1 shows the mean and standard deviation of the estimated angular acceleration, for each bird, stroke phase and specific wing. As expected, during the upstroke and downstroke, the angular acceleration and its standard deviation are relatively small for the left and right wings of the set of experimental birds. In contrast, during the transition phase the averaged angular acceleration and the standard deviation increase. Interestingly, significant differences were found in the transition phase between the right and left wing for the first and fourth bird.

Yet manual annotations are inevitably influenced by the inter subject variability and therefore they are prone to errors, such estimations were herein used as a reference frame of the flight patterns. A comparison of the computed kinematic paths was carried out w.r.t estimated motion trajectories from the expert segmentations. For doing so, each wing was delineated for the complete stroke sequence of every available flight. Then, the global angular acceleration was computed as the second derivative of the Center-of-Mass trajectory of every wing. The resultant mean and standard deviation of the angular acceleration obtained from manual delineations is plotted in Figure 3 (green line), and a close match with the automatic estimations is observed in this graph, a statement supported by a correlation factor of 0.84 ± 0.12.

A local angular acceleration profile was computed by dividing the wing into segments, the angular acceleration being regionally estimated for each segment during the flight stroke cycle. Figure 3 shows a typical map acceleration profile, during a cycle and using eight wing segments. Each computed map reveals further details about how the angular acceleration are distributed along the wing, from the proximal to distal regions along the stroke flight. As 1. The mean and standard deviation of the angular acceleration are presented for the right and left wing of each of the four birds. Such analysis is independently performed for each of the different stroke phases. During the upstroke and downstroke phases, the angular acceleration showed not significant differences among the different bird strokes (t-student test and p < 0.05) since the the flight is basically passive. In contrast, during the transition phase the averaged angular acceleration and standard deviation increased since the motion is highly dependent of the bird. The significant differences found in this phase (t-student) test, p < 0.05 are marked with the asterisk. expected, in the first phase of the stroke cycle, the acceleration profile is very similar for the proximal and distal parts. At about a 30 % of the cycle, this pattern changes and the distal part is fully responsible for the observed motion changes. The proximal wing simply leverages the distal motion whose amplitude varies in several mini-cycles to a larger extension at the end of the cycle. From such acceleration map, it was also possible to observe the non linear distribution of the angular acceleration in both axes: wing position and temporal cycle.

Finally, a global wing deformation, during a stroke cycle, was estimated from the whole experimental group W l t (resp. W r t ), and shown in Figure 5. In general, every series of the wing deformation follows a very similar pattern, with a variance whose amplitude is smaller than a 10 % of the spanned wing amplitude. As for other deformation estimations, a statistic boxplot was computed at each time of the cycle, using every recorded video, and again the red line represents the sample median while the spread of data is represented by the blue lower and upper quartiles. The larger stroke deformations correlate with the maximum wing effort moments, when motion abruptly rises and there exists a change of direction. An independent analysis for the three flight stroke phases was also performed for the global deformation patterns. As expected, during the downstroke and upstroke phases, when there exists a relatively relaxed state and motion is almost constant, the act like rigid structures with small deformations, principally at the wing borders. The pattern is systematically reproduced by the different wings, with a very small standar deviation (around 1.96 rad f rame

2

). On the other hand, during the transition stroke phase, the deformation increases, as well as the variance. During this phase, the deformation is large because of the rapid wing adaptation to the new change of acceleration direction.

This dense optical flow software was developed in house by the authors and was implemented by authors in c++. This particular dense optical flow has been already tested in other motion applications [START_REF] Manzanera | Local Jet Feature Space Framework for Image Processing and Representation[END_REF], [START_REF] Martínez | Eduardo Romero Representing activities with layers of velocity statistics for multiple human action recognition in surveillance applications[END_REF]

Discussion

This work introduced a novel approach to perform a complementary kinematic analysis of the hummingbird wing motion, using a multidimensional dense optical flow and without any prior information or any invasive sensor that could alter the natural flight motion gesture. This markerless strategy opens up the possibility of capturing additional spatiotemporal relationships since this is a much more complete representation of the wing dynamic than the classical points of interest that are marked by an expert along a recorded video sequence. While optical flow approaches such as Lukas-Kanade are relatively robustness to noise, the resulting flow basically tracks certain salient features, in particular the edges [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF], [START_REF] Bouguet | Pyramidal implementation of the affine lucas kanade feature tracker description of the algorithm Intel Corporation[END_REF]. Alternatively, global approaches, like Horn and Schunk, are limited to smooth vector fields that capture slow changes but are also more sensible to noise because the flow representation is obtained at the pixel level. A hummingbird flight analysis from this type of dense flows may lose relevant flight features, namely, during the transition between the upstroke and downstroke phases. In contrast, the optical flow herein applied extracts relevant information rather than the apparent location of individual pixels, and therefore ends up by being a suitable and robust tool to determine the fundamental motion wing patterns. textcolorMyColorThe herein proposed approach is able to compute the angular velocity and to estimate the wing deformation during the hummingbird yaw turning using a simple dense velocity field computed from video sequences. Since this particular dense optical flow is basically a computation of the spatial changes at multiple scales by using the different derivatives, this approach results to be very robust to the noise. In addition, the temporal correlation is used to remove different types of illumination noise or artifacts produced by the capturing procedure. The main limitation, in contrast, comes from the number of scales used in terms of computational time. The whole picture is reconstructed by means of a set of kinematic measurements computed from the velocity optical flow field, which was in due turn used to spatially segment and cluster the different wing flow patterns. Angular acceleration and temporal wing deformation measures were herein estimated from local and global perspectives over a set of 91 videos. Quantitative metrics and results showing a low variance of the computed kinematic patterns among the evaluated sample as showed in the variance analysis evaluation.

Flight analyses by means of physical approximations are often based on the study of some temporal sparse biomarkers, and therefore only provide temporal information of very simplified structures, or become statistically insignificant. These approximations in consequence result insufficient or unsuitable towards understanding complex maneuvers, especially when such motions are directly related with temporal non-linear structural deformation. Some computational models partially overcome such simplicity, allow for simulation of several flight patterns by subdividing the wing into simpler parts which are steered by independently parameterized physical models [START_REF] Liang | Computational Analysis of Hovering Hummingbird Flight In 48th[END_REF]. However, such simulations are usually unrealistic and far from being interpretable, they would hardly discover complicated relationships among forces, interactions, neuromotor commands and joint movements that might produce complex motions such as the hang up flight. In contrast, the proposed approach is able to describe convolved kinematic patterns that are spatially spread out along the wing, using exclusively the video information, without prior hypotheses. In this sense, the presented strategy not only can find out reliable motion patterns, but it also may reveal local deformation patterns of the wing .

Rotational patterns are crucial in flight analysis since they allow to study the source of circulation air which is straightforwardly related to the upward and downward forces during the stroke cycle. Furthermore, the angular acceleration, the angle of attack and the span ratio are responsible for controlling the upstroke circulation, the wing orientation and the maintenance of the forward flight [START_REF] Dickinson | Wing Rotation and the Aerodynamic Basis of Insect Flight Science[END_REF], [START_REF] Tobalske | Biomechanics of bird flight[END_REF]. The Global angular acceleration herein estimated presents large gradients during the upstroke and downstroke transitions, a result of the upward and downward wing forces which are responsible for the circulation changes. From a local angular acceleration standpoint, the proposed approach also allows to estimate a regional angular acceleration profile along the wing, thereby discovering not only the phases where the sources of circulation change but also determining which particular wing regions are more involved in the generation of these forces. As illustrated in Figure ??, the main angular acceleration modes are initially concentrated in the proximal wing regions because the shoulder and pectoral muscles drive the forward wing motion, a pattern that is displaced to the distal wing regions when the passive accumulated energy rules the downstroke phase. Although the proposed method is able to capture local motion patterns from the dense optical flow, the spatial accuracy is restricted by the number of scales and the standard deviation used to convolve each frame. For instance, when the scales are obtained by convolving Gaussians with large sigmas, results would be more robust to noise from a motion but less accurate spatially. A proper balance between the noise control and the spatial accuracy is part of the method tuning as to obtain coherent results.

On the other hand, the stunning deformations of the wing are responsible for much of the salient abilities of hummingbirds such as their stability and relative conservation of the orientation, even if they are forced to change constantly their acceleration patterns [START_REF] Ennos | The kinematics and aerodynamics of the free flight of some Diptera[END_REF], [START_REF] Lehmann | Elastic deformation and energy loss of flapping fly wings[END_REF]. In addition, localized wing deformations may produce relative point circulation vortices that induce a non linear pattern along the wing structure and therefore completely alter the hummingbird flight smoothness [START_REF] Hedrick | Estimates of circulation and gait change based on a three-dimensional kinematic analysis of flight in cockatiels (Nymphicus hollandicus) and ringed turtledoves (Streptopelia risoria)[END_REF]. In the proposed strategy, deformation wing patterns were locally estimated along the wing as the relationship between each local orientation w.r.t an average orientation, resulting in a wing map of deformation for every stroke. Main wing deformations are commonly localized along the proximal and distal segments while the largest deformation is observed when the wing acceleration is maximum, i.e., during the downstroke and upstroke transitions (see Figure 5).

Recent methodological advances have improved our comprehension of the most complex hummingbird maneuvers and their intrinsic mechanisms, for instance, aerodynamic studies using PIV have shown an asymmetric distribution of the lift forces during the stroke, probably caused by unknown micro-mechanisms like the stall and leading-edge vortices [START_REF] Warrick | Aerodynamics of the hovering hummingbird[END_REF]. New analysis strategies are then required to approach such local and non linear phenomena whose role is yet an open problem. The proposed approach herein introduced is a global-to-local analysis tool that complements the existent methods, with the benefit of methodological simplicity. Likewise, the two variables herein estimated, angular acceleration and wing deformation, were locally and globally computed, opening thereby the possibility of conducting more subtle studies that discover the wing regions that contribute to every step of the hummingbird yaw turn.

The kinematic patterns herein reported were captured via a single camera from only a ventral perspective, however this analysis can be directly extended to analyze the complementary dynamic patterns in the other planes. Likewise, the proposed approach can easily adapted to a 3D reconstruction by using multiple calibrated cameras. The proposed approach may also able to locally compute force and torque patterns at the different wing segments by directly relating the angular acceleration with the segmented anatomical structured. Finally, this method can be used to automatically segment the hummingbird flight strokes in sequences of video, highly reducing inter and intra expert variability. In addition, the proposed approach can be potentially applied even when the animal is moving along the scene, since once the dense optical flow tracks the main regions in motion, this bounded region may be followed and characterized.

Conclusions

This work has presented a novel approach to quantify the yaw turn hummingbird motion, using a dense optical flow computed from video sequences. A main contribution is the markerless flight analysis framework that estimates two relevant kinematic hummingbird variables i.e., the angular acceleration and wing deformation. The proposed method opens up an actual possibility of understanding more complex yaw turn phenomena. Future works include the study and development of a robust method that temporally relates the kinematic dense flow description with force measurements generated by flapping, under actual anatomical and mass parameters.

Figure 1 .

 1 Figure 1. In (a) is shown the ventral configuration perspective used to record the hummingbird flight. In (b) is illustrated the computation of the N wing segments using as reference a rotated bounding box

  (a) hummingbird flight video sequence (b) Dense optical flow computation during the hummingbird flight. The optical flow field in each frame it is represented by a set of green arrows. (c) Deformation map obtained from the orientation variance w.r.t to the average orientation in each frame. In this color scale the larger deformation is represented by the red color.

Figure 2 .

 2 Figure 2. Hummingbird wing motion analysis from a dense optical flow computation.

Figure 3 .

 3 Figure 3. In (a) the artificial rotating plate mechanism is shown using the proposed approach to test. In (b) a typical dense optical flow computation from the artificial mechanism is illustrated, notice the green arrows represent the speed vectors.

  (a) Angular acceleration for the right wing (b) Angular acceleration for the left wing.

Figure 4 .

 4 Figure 4. The statistical box plot sequences represent the estimated angular acceleration for both wings using the proposed approach within which the red line represents the median of the angular acceleration, the blue lower and upper boxes stand for the spread of data (interquartile range) and the maximum and minimum values are the vertical dotted lines. The green box plot sequence corresponds to the manual estimations performed by an expert.

Figure 5 .

 5 Figure 5. Local angular acceleration profile by the subdivision of the wing eight segments.

  (a) Global wing deformation for the right wing (b) Global wing deformation for the left wing.

Figure 6 .

 6 Figure 6. The statistical box plot sequences represent the estimated gobal wing deformationn for both wings using the proposed approach

  ± 0.06 0.8 ± 0.36 -0.1 ± 0.03 0.1 ± 0.05 0.3 ± 0.46* -3.3 ± 0.51*

	upstroke	downstroke	transition
	Right	Left	Right	Left	Right	Left
	B.1 0.2 ± 0.75	0.2 ± 0.2	0.1 ± 0.11 -0.2 ± 0.08 -0.4 ± 0.36* -0.1 ± 0.46*
	B.2 0.3 ± 0.06 -0.3 ± 0.03 -0.1 ± 0.46 0.1 ± 0.68	0.3 ± 0.62	0.1 ± 0.76
	B.3 0.1 ± 0.28 -0.1 ± 0.38 0.1± 0.38 -0.1 ± 0.36 -0.1 ± 0.44	0.2 ± 0.37
	B.4 -0.7 Table

Acknowledgements

The data used in the present investigation was kindly provided by the department of Zoology of the University of British Columbia. We wish to thank Dr Douglas Altshuler, Elsa Quicazan, and Dimitri Skandalis for their discussions that highly enriched the present work.


R Azouit

A Sarlette

P Rouchon

Convergence and adiabatic elimination for a driven dissipative quantum harmonic oscillator *

We prove that a harmonic oscillator driven by Lindblad dynamics where the typical drive and loss channels are two-photon processes instead of singlephoton ones, converges to a protected subspace spanned by two coherent states of opposite amplitude. We then characterize the slow dynamics induced by a perturbative single-photon loss on this protected subspace, by performing adiabatic elimination in the Lindbladian dynamics.

Introduction

The harmonic oscillator is a standard quantum system. It features coherent states, equivalent of classical harmonic oscillator amplitudes, whose coherent superpositions also called "cat states" feature genuinely quantum properties with no classical equivalent. In a recent paper [START_REF] Mirrahimi | Dynamically protected cat-qubits: a new paradigm for universal quantum computation[END_REF], our collaborators take advantage of this fact to propose an implementation of logical quantum bits as "cat states", with potential to realize universal quantum computation using standard technological elements as quantum gates. Their key contribution, besides the insight that the cat states are inherently insensitive to part of the typical perturbations, is the design of an "engineered reservoir" that stabilizes a subspace of such states in open loop.

The contribution of the present paper is to precisely establish, from the Lindblad master equation [START_REF] Carmichael | An Open Systems Approach to Quantum Optics[END_REF], the stabilization properties of this scheme both,

• in ideal situations with = 0, by proving in theorem 1 global convergence of the infinite-dimensional nominal model towards the target "protected subspace",

• in presence of a small but dominant decoherence source with > 0, by establishing the approximate slow dynamics as a reduced Lindblad master equation (15).

The first goal builds on a typical Lyapunov-LaSalle strategy, with additional care needed due to the infinite dimension. For the second goal, we resort to a separation of the quantum dynamics into fast and slow components. We then apply an adiabatic elimination of the fast system to deduce a good approximation of the dynamics on a slow manifold, with the state remaining -close to the protected subspace. Studying such perturbations is standard for quantum Hamiltonian systems, where regular perturbation theory can be routinely applied [START_REF] Sakurai | Modern quantum mechanics[END_REF], but the Lindbladian case with singular perturbations has attracted much less attention. In [START_REF] Mirrahimi | Singular perturbations and Lindblad-Kossakowski differential equations[END_REF] and similarly [START_REF] Reiter | Effective operator formalism for open quantum systems[END_REF] singular perturbations up to second order are applied to a system with N ground states and eliminating relaxing excited states. In [START_REF] Atkins | Approximate master equations for atom optics[END_REF][START_REF] Warszawski | Adiabatic elimination in compound quantum systems with feedback[END_REF] specific atom optics dynamics and an ancilla-mediated feedback are investigated with the standard approach of [START_REF] Carmichael | An Open Systems Approach to Quantum Optics[END_REF]. In [START_REF] Kessler | Generalized Schrieffer-Wolff formalism for dissipative systems[END_REF] the so-called Schrieffer-Wolff formalism is generalized to Lindbladian dynamics; its basic form requires inversion of the nominal dynamics operator, which is not too practical and which we circumvent here for the derivation of the reduced slow master equation ( 15). The paper is organized as follows. Section II describes the mathematical model of the dynamics to be studied. Section III provides the global convergence proof for an idealized model, and Section IV analyzes precisely how this allows to counter the dominant external perturbation, which is single-photon loss. We pedagogically present the corresponding slow/fast perturbative argument (Section IV.B) as the translation to quantum notation of the standard dynamical systems approach, recalled in Section IV.A. Finally, simulations illustrate the validity of our analysis in Section V.

Driven dissipative pairwise photon process

The underlying space of the quantum harmonic oscillator is a Hilbert space H of infinite dimension spanned by the Fock states {|n } n∈N . The annihilation operator a is defined by a|n = √ n|n -1 for any n ≥ 1 and a|0 = 0. Its Hermitian conjugate a † , verifies a † |n = √ n + 1|n + 1 , for any n ≥ 0; . We denote N = a † a the photon number operator, satisfying N|n = n|n , for any n ≥ 0. For any α ∈ C, a coherent state |α ∈ H is characterized by a|α = α|α and

|α = e -|α| 2 2 ∞ n=0 α n √ n! |n .
The coherent states |α are viewed as the quantum model for a "classical state" of complex amplitude α. Indeed, whereas the classical undamped harmonic oscillator |αe iωt up to an irrelevant phase. As is customary, in the rest of the paper we describe the system in a frame rotating at the oscillator frequency ω, for which the oscillator has stationary solutions |α . We call Schrödinger cat state or simply cat state the coherent superposition of two coherent states with opposite amplitudes,

|c ± α = |α ± | -α γ ± (1)
where [START_REF] Tarasov | Quantum Mechanics of Non-Hamiltonian and Dissipative Systems[END_REF]). The set of density operators (quantum states) corresponds to elements of K 1 (H) that are non-negative and of trace 1. They are usually denoted by ρ. We denote also by K f (H), the subspace of K 1 (H) of operators whose range is included in a vector space spanned by a finite number of Fock states:

γ ± = 2(1 ± e -2|α|
K f (H) =          finite f n,n |n n | : f n,n ∈ C, f n,n = f * n ,n          .
We consider the quantum harmonic oscillator interacting with its environment as described in [START_REF] Mirrahimi | Dynamically protected cat-qubits: a new paradigm for universal quantum computation[END_REF]. An external coherent driving field of amplitude u (assumed here real and strictly positive without loss of generality) is applied such that the oscillator can only exchange photons in pairs. Furthermore, the quantum system is built such that similarly the main dissipative process is a pairwise photon loss with rate κ > 0. Nevertheless, due to physical constraints, even if a single photon loss process is much less frequent than the previous one, it can't be totally omitted. For this reason, this term appears with a small coefficient 0 < κ in the following Lindblad master equation which governs the system dynamics:

d dt ρ = u[(a † ) 2 -a 2 , ρ] + κL a 2 (ρ) + L a (ρ)
where [•, •] stands for the commutator and where, for any linear operator A on H, the super-operator L A is given by

L A (ρ) = AρA † -( A † Aρ + ρ A † A)/2 .
Noting α = 2u/κ and L = a 2 -α 2 one can reformulate the previous equation as :

d dt ρ = κL L (ρ) + L a (ρ) (2)
It is shown in [START_REF] Mirrahimi | Dynamically protected cat-qubits: a new paradigm for universal quantum computation[END_REF] that, for = 0, the two-dimensional Hilbert space

H α = span |α , | -α
is a decoherence-free space: for = 0, any density operator ρ with support included in H α is a steady state, i.e., L L ( ρ) = 0. We will not investigate here the well-posedness of (2) and the associated strongly continuous semigroup of linear contraction on K 1 (H). Such issues can be investigated via theorem 3.1 of [START_REF] Davies | Quantum dynamical semigroups and the neutron diffusion equation[END_REF] ensuring the existence of minimal solutions since for any α, κ, > 0, the operator -κL † La † a is the infinitesimal generator of a strongly continuous one parameter contraction semigroup on H. This means that, in the sequel, we will always assume that the Cauchy problem (2) with an initial condition ρ 0 ∈ K 1 (H), positive semidefinite and of trace one, admits a solution in K 1 (H) defined for any t > 0. Moreover, this minimal solution will remain positive semidefinite. We will investigate here the time asymptotic regime for = 0 and then for 0 < κ.

3 Convergence of (2) for = 0

Lemma 1. For any quantum state ρ in K f (H) we have Tr LL L (ρ)L † ≤ -2Tr LρL † . Proof. From Tr LL L (ρ)L † = Tr LρL † [L † , L] and [L † , L] = [(a † ) 2 , a 2 ] = -4N- 2I, we have Tr LL L (ρ)L † = -2Tr LρL † (2N + I) .
We conclude since N and LρL † are positive semidefinite.

Modulo arguments proper to the infinite-dimensional setting, this basically means that V(ρ) = Tr LρL † is an exponential Lyapunov function for the system (2) with = 0. Lemma 2. For any ν ≥ 1 there exists µ ≥ 0 such that, for any quantum state ρ in

K f (H) we have Tr L L (ρ)N ν ≤ -ν Tr ρN ν ν+1 ν + µ Proof. Denote by L * L the adjoint of L L : L * L ( A) = L † AL -(L † L A + AL † L)/2
for any Hermitian operator A on H. Computations relying on the identity a f (N) = f (N + I)a for any function f , yield

L * L ( f (N)) = -N(N -1)( f (N) -f (N -2I)) + α 2 2 a 2 ( f (N) -f (N -2I)) + α 2 2 ( f (N) -f (N -2I))(a † ) 2 . Since Tr L L (ρ) f (N) = Tr ρL * L ( f (N)) , we have Tr L L (ρ) f (N) = -Tr ρN(N -1)( f (N) -f (N -2I)) + α 2 2 Tr ρa 2 ( f (N) -f (N -2I)) + α 2 2 Tr ρ( f (N) -f (N -2I))(a † ) 2 . Take f (x) = x ν and define g(x) = f (x) -f (x -2) for x ≥ 2; g(x) = f (x) for 2 > x ≥ 0; g(x) = f (0) for x < 0. From a 2 g(N) = a g(N -I)a g(N)
and Cauchy Schwartz inequality

|Tr ρa 2 g(N) | = Tr √ ρa g(N -I) a g(N) √ ρ ≤ Tr ρg(N -2I)(N + I) Tr ρg(N)N . Since g(x -2)(x + 1) ≤ g(x)(x + 3) for x ≥ 0, we have Tr ρa 2 g(N) ≤ Tr ρg(N)(N + 3I) . Thus 1 2 Tr ρ(a 2 g(N) + g(N)(a † ) 2 ) ≤ Tr ρg(N)(N + 3I) and Tr L L (ρ)N ν ≤ Tr ρg(N) -N 2 + (α 2 + 1)N + 3α 2 . Since (x ν -(x -2) ν ) -x 2 + (α 2 + 1)
x + 3α 2 is equivalent to -2νx ν+1 for large x, there exists µ > 0 such that for all x ≥ 0, 1. Take ν ≥ 1. Then there exists γ > 0 such that, for any initial quantum state ρ(0) = ρ 0 satisfying Tr ρ 0 N ν < +∞, we have for all t > 0, Tr ρ(t)N ν ≤ max γ, Tr ρ 0 N ν .

((x + 2) ν -x ν ) -x 2 + (α 2 + 1)x + 3α 2 ≤ -νx ν+1 + µ.
2. Assume that ρ 0 ∈ K f (H). Then, there exists a quantum state ρ with support in H α such that, for any ν ≥ 0, lim

t →+∞ Tr N ν 2 (ρ(t) -ρ)N ν 2 = 0.
The limit ρ depends on ρ 0 . It can be derived from ρ 0 with the four Hermitian, bounded and independent operators that are in the kernel of the adjoint superoperator L * L and given in [START_REF] Mirrahimi | Dynamically protected cat-qubits: a new paradigm for universal quantum computation[END_REF].

Proof. The first statement is a direct consequence of Lemma 2 and of the fact that quantum states element of K f (H) are a dense subset of the quantum states σ of K Banach space. From the first statement, we know that ρ 0 being an element of

K 1 ν (H), ρ(t) remains always in K 1 ν (H). Since ν > ν, the injection of K 1 ν (H) into K 1 ν (H) is compact: {ρ(t) | t ≥ 0} is precompact in K 1 ν (H). Denote by ρ ∈ K 1
ν (H) an adherent point of ρ(t) for t tending towards infinity. Since ν ≥ 2, ρ and ρ belong to the domain of L L . Lemma 1 implies Tr L ρL † = 0, i.e., the support of ρ is contained in the kernel of L, which coincides with H α . Moreover, the semigroup associated to the Lindblad master equation is a contraction for the trace distance: for two trajectories ρ 1 (t) and ρ 2 (t), t → Tr |ρ 1 (t) -ρ 2 (t)| is a nonincreasing function. Thus t → Tr |ρ(t) -ρ| is non-increasing since ρ is a steady state. Consequently the adherent point ρ is unique: ρ(t) converges towards ρ in

K 1 ν (H).
4 Reduced slow dynamics of (2)

We have proved in the previous section that the system converges toward the decoherence free subspace H α when we neglect the photon loss channel ( = 0). When 0 < 1, the center manifold theorem allows us to separate the system into fast and slow dynamics. The fast dynamics makes the system globally converge to a subspace close to H α . The slow dynamics approximate the behavior of "protected states" |c + α , |c - α . The present section is aimed at characterizing these dynamics to the first order in .

The fast/slow dynamics reduction from nonlinear systems theory, also known as singular perturbation theory, is useful despite the linearity of Lindbladian dynamics, because the very high dimension and often high degeneracy of open quantum systems makes matrix diagonalization impractical to apply.

For the sake of clarity, we first particularize the fast/slow dynamics reduction theory to linear systems of finite dimension with the standard notations. We propose then a reduction procedure via an adapted duality viewpoint with the characterization based on ( 9) and (10) here below. We apply then this characterization to the quantum system (2) and show that it facilitates the computations up to first order with respect to e.g. [START_REF] Kessler | Generalized Schrieffer-Wolff formalism for dissipative systems[END_REF].

Reducing a linear system to its slow dynamics

We here review the theory of geometric singular perturbation, which was mainly developed by Fenichel in [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] and surveyed by Jones in [START_REF] Jones | Geometric singular perturbation theory[END_REF], in a linear context. Consider the nominal linear system d dt x = Ax with x = (x 1 , x 2 ) ∈ R m × R n and converging globally to the m-dimensional subspace S = {x ∈ R m+n : x 2 ∈ R n = 0}. To this nominal dynamics we add an arbitrary perturbation matrix B of order 1. In matrix notation, the dynamics can be written in block form which yields:

d dt x 1 = A 1 x 2 + (B 1 x 2 + B 0 x 1 ) (3)
d dt x 2 = A 2 x 2 + (B 2 x 2 + B 3 x 1 ) .
The assumption that S is globally exponentially stable for = 0 corresponds to A 2 having all eigenvalues with strictly negative real parts, thus it is invertible. Hence we can define the regular change of variables

x1 = x 1 -A 1 A -1 2 x 2 x2 = x 2 (4)
which yields dynamics in Tikhonov normal form

d dt x1 = (B 0 -A 1 A -1 2 B 3 ) x1 (5)
+(B 0 + B 1 -A 1 A -1 2 (B 2 + B 3 )) x2 = f ( x1 , x2 ) . d dt x2 = A 2 x2 + (B 3 x1 + (B 2 + B 3 A 1 A -1 2 ) x2 ) (6) = g( x1 , x2 , ) .
The Tikhonov conditions for reducing the system by singular perturbations is that the first (slow) subsystem has eigenvalues going down as , while the second (fast) subsystem has eigenvalues bounded away from zero for = 0. These conditions are satisfied above. The Tikhonov theorem then allows the following reduction.

Proposition 1. The trajectories of the full system (5),(6) (with initial conditions satisfying g = 0) remain -close over at least a time of order 1/ , to the trajectories of the system restricted to the "slow submanifold" g( x1 , x2 , 0) = 0 and where dynamics are given by replacing x2 in f ( x1 , x2 ) by the solution of g( x1 , x2 , 0) = 0.

In our linear case, the slow manifold comes down to x2 = 0 and the slow dynamics trivially reduce to the first term in [START_REF] Jones | Geometric singular perturbation theory[END_REF]. Transforming back to the original coordinates the slow manifold corresponds just to x 2 = 0 and the dynamics are

d dt x 1 = (B 0 -A 1 A -1 2 B 3 )x 1 . (7)
The second term reflects the influence of the fast x 2 dynamics on the slow variable x 1 : by blindly setting x 2 = 0 in the original system (3) and neglecting its second line, we would miss this term and get an incorrect approximation.

Computing the corrective term A 1 A -1 2 B 3 by explicit inversion of A 2 can be a tedious task when the fast subsystem has a large dimension (in our quantum case, x 2 would rigorously be of infinite dimension). However if first integrals of the system with = 0 are known, these can facilitate the computations via the following dual viewpoint.

Consider a linear functional p T = (p T 1 , p T 2 ) ∈ R * m+n which is conserved by ẋ = Ax, i.e. satisfying ṗ = A T p = 0 or equivalently p T 1 A 1 x 2 + p T 2 A 2 x 2 = 0 for all x 2 ∈ R n . (The notation • T denotes matrix transpose.) Again using invertibility of A 2 , we see that p T actually satisfies

p T 1 A 1 A -1 2 y + p T 2 y = 0 ∀y ∈ R n . ( 8
)
Knowing m linearly independent functionals (p T (k)) k=1,...,m satisfying ( 8) is sufficient to fully characterize the corrective term in ( 7): we have

d dt x 1 = (B 0 + Q)x 1 (9)
with Q defined by the set of linear equations:

p T 1 (k)Q = p T 2 (k)B 3 , k = 1, 2, ..., m . (10)
4.2 Quantum system (2) with 0 < κ

We know apply the same procedure to our quantum system. The nominal ẋ = Ax corresponds to ρ = L L (ρ). We have shown in Section 3 that:

• S corresponds to a four-dimensional real subspace of Hermitian operators spanned by

|c + α c + α |, |c - α c - α |, |c + α c - α | + |c - α c + α |, and i(|c + α c - α | -|c - α c + α |). These correspond to the m = 4 coordinates of x 1 .
• The subspace S is globally asymptotically stable under ρ = L L (ρ). This corresponds to the invertibility condition on A 2 independently of .

We therefore introduce the projector

P c = |c + α c + α | + |c - α c - α | (11)
such that ρ s = P c ρP c [START_REF] Tarasov | Quantum Mechanics of Non-Hamiltonian and Dissipative Systems[END_REF] corresponds to the "slow" x 1 space of the previous section. The projection onto "x 2 space" is given by ρ f = ρ -P c ρP c .

Our goal is to compute the evolution of ρ s , which is the equivalent of [START_REF] Bris | Low-rank numerical approximations for highdimensional lindblad equations[END_REF]. For this we take advantage of the dual formulation ( 9), [START_REF] Reiter | Effective operator formalism for open quantum systems[END_REF]. The following procedure can in principle be applied to any perturbative dynamics, we here focus on L a as a physically relevant case.

The perturbative dynamics on the slow manifold features a first component, corresponding to B 0 , obtained simply by projection onto the slow manifold. The identities a|c ± = α γ ∓ γ ± |c ∓ quickly yield its explicit expression:

P c L a (ρ s )P c = α 2 L X (ρ s ) ( 13
)
where

X = γ + γ - |c + c -| + γ - γ + |c -c + | . ( 14
)
To compute the corrective term by duality, using the equivalent of ( 9), [START_REF] Reiter | Effective operator formalism for open quantum systems[END_REF], we need to identify m = 4 conserved functionals of the system, which are the fixed points of the dual nominal dynamics d dt ξ = L * L (ξ). Fortunately, those invariants are known for the particular operator L:

• One easily checks that ξ a = I the identity operator is in the kernel of any L * L .

• The parity operator ξ b = (-1) a † a is in the kernel of L * L because photons are exchanged by pairs.

• The appendix of [START_REF] Mirrahimi | Dynamically protected cat-qubits: a new paradigm for universal quantum computation[END_REF] gives two more operators ξ c and ξ d in terms of Bessel functions; one checks that they are linearly independent for finite α.

We will also use the following key property of the conserved quantities, which is specific to the structure of quantum Lindblad dynamics.

Lemma 3. Any Hermitian operator ξ in ker(L * L ) commutes with P c the orthogonal projector onto H α .

Proof. From 2L † ξL = L † Lξ + ξL † L and LP c = 0 = P c L † we have L † Lξ P c = 0 = P c ξL † L. Thus the Hermitian operator A = P c ξ + ξ P c satisfies L AL † = (L † L A + AL † L)/2, i.e. belongs to the kernel of L L . The support of A is thus included in H α and thus [P c , A] = 0. This implies that P c ξ = P c ξ P c = ξ P c .

Considering projections on respective subspaces, the equivalent of the B 3 term of the linear perturbation is given by L a (ρ s ) -P c L a (ρ s )P c . The equivalent of equation [START_REF] Reiter | Effective operator formalism for open quantum systems[END_REF] characterizes the Hermitian operator Q with support on H α as follows:

Tr P c ξ ν P c Q = Tr ξ ν -P c ξ ν P c L a (ρ s ) -P c L a (ρ s )P c = Tr ξ ν L a (ρ s ) -P c L a (ρ s )P c = 1 2 Tr ξ ν (P c -I)a † aρ s + ρ s a † a(P c -I) = 1 2 Tr ξ ν a † aρ s (P c -I) + (P c -I)ρ s a † a = 0
for any ν = a, b, c, d. From the first to the second line, P c ξ ν P c is readily dropped since P c is a projector. For the next one we have to write out L a and see that P c aρ s a † P c = aρ s a † for the particular perturbation operator a. The last line follows by using Lemma 3 and P c ρ s = ρ s . As a conclusion, we get that the corrective term is Q = 0 for our particular case. We can summarize these computations as follows.

The trajectories of system (2) with initial conditions having supports in H α remain -close over at least a time of order 1/ , to the trajectories of the "slow variable" ρ s which is a linear combination of |c

+ α c + α |, |c - α c - α |, |c + α c - α | + |c - α c + α |, and i(|c + α c - α | -|c - α c + α |).
The slow state ρ s follows the Lindblad dynamics:

d dt ρ s = α 2 L X (ρ s ) ( 15
)
where

X = γ + γ - |c + c -| + γ - γ + |c -c + | . (16)
In applications considering |c + α and |c - α as canonical states |0 , |1 of a logical qubit [START_REF] Mirrahimi | Dynamically protected cat-qubits: a new paradigm for universal quantum computation[END_REF], the operator X corresponds to a bit-flip in the limit γ + γ -→ 1 of large coherent amplitude α and to a decoherence to the vacuum |0 in the limit of Fock states |c + α = |n = 0 , |c - α = |n = 1 when α = 0. For all other cases, the qubit dynamics (15),( 16) corresponds on the canonical Bloch sphere to:

d dt x = -α 2 (γ 2 + -γ 2 -) 2 2γ 2 + γ 2 - x d dt y = -α 2 (γ 2 + + γ 2 -) 2 2γ 2 + γ 2 - y d dt z = -α 2 γ 4 + + γ 4 - γ 2 + γ 2 - z - γ 4 + -γ 4 - γ 4 + +γ 4 - .
This converges to x = 0 (slowly for γ + γ - 1), y = 0, and z =

γ 4 + -γ 4 - γ 4 + +γ 4 - (which is 0 for γ + γ - 1).

Numerical simulations

To illustrate the interest of this model reduction based on adiabatic elimination of the rapidly converging variables, we compare via numerical simulations trajectories of the complete system and of the reduced one. To do so, we use a numerical scheme which preserves the positiveness for the Lindblad equation and similar to one used in [START_REF] Bris | Low-rank numerical approximations for highdimensional lindblad equations[END_REF]. We choose the following values of the parameters : a time-step of 10 -3 , the decoherence strength κ = 1, u = 1/2 and = 0.01. With α = 1, the population of photons for n > n max = 40 in the coherent state |α is almost zero since less that 1 n max ! . Consequently, we truncate the infinite-dimensional Hilbert space to a numerical system space spanned by {|1 , |2 , . . . |40 } in the Fock basis. We denote by ρ the density matrix of the resulting system. It is worth stressing that the complete system is represented by an n max × n max matrix while the reduced system is represented by a 2 × 2 matrix (on the basis |c + α , |c - α ). Thus the computation is much faster on the second one.

We first take as initial condition the vacuum state, ρ 0 = |0 0|. The state of the reduced system, ρ s is then initialized at |c + α c + α | because both ρ 0 and |c + α c + α | are +1 eigenstates of the parity operator ξ b = (-1) a † a , which is a conserved quantity (see Section 4.2 and [8]). To compare the trajectories of (2) initialized at |0 0| and of (15) initialized at |c + α c + α |, we show in figure 1 the expectation values Tr ρσ z and Tr ρ s σ z of the operator

σ z = |c + α c + α | -|c - α c - α |, commonly denoted σ z .
After a transitional regime of typical duration 1/κ, one can see a strong similarity between Tr ρσ z and Tr ρ s σ z up to a constant offset. The value of this offset is of order . Furthermore, we plot the fidelity F(ρ, ρ s ) = tr √ ρ s ρ √ ρ s between ρ s and ρ. For better readability, figure 2 shows the logarithm of 1 minus the fidelity, i.e. of its deviation from the ideal value 1. This deviation quickly converges to an order 10 -4 , corresponding to 2 as expected. It then further decreases, incidentally, as both systems converge towards the unique equilibrium of the slow dynamics.

To emphasize the influence of γ + and γ -in (16), we add a simulation with the same parameters but with the following and same initial condition for the complete and reduced system:

ρ0 = 1 2 |c + α + |c - α c + α | + c - α |
Figure 3 shows that the expectation value of σ x = |c + α c - α | + |c - α c + α | slowly decreases over time, as expected from bit-flip dynamics. The slope of this decrease is approximated to ∼ 4% accuracy by the reduced dynamics. Moreover, figure 4 establishes that σ z does not remain zero. This is due to the fact that, with γ + > γ -, equation (15) "promotes" the population of |c + α c + α | over the population of |c - α c - α |, unlike a pure bit-flip.

The simulations thus confirm the validity of our approximation of the complete model by the reduced one.

CONCLUSIONS

We have rigorously proved convergence of a harmonic oscillator Lindblad dynamics with two-photon exchanges, to a protected subspace. We have also established the approximate slow dynamics on this protected subspace when a typical perturbation is added, and illustrated its validity in simulations. The methods used for this particular example are applicable to general Lindbladian dynamics.

The reduction by singular perturbations and adiabatic elimination is of course applicable in general to evaluate the remaining slow dynamics in quantum systems with (engineered) protected subspaces. Extension to k-photon processes a k with k > 2 can be addressed in the same way. The fact that the slow variable still follows a Lindbladian master equation may not be surprising but remains to be proved in the general case. The fact that the dynamics reduces to the orthogonal projection of the Lindbladian onto the protected subspace (i.e. B 0 without any correction due to B 3 , in the terms of Section 4.1) for the case examined here would be in agreement with the physicists' "quantum Zeno" viewpoint. However, under which formulation this viewpoint should be applied in the general case also remains to be rigorously characterized.

d

  dt x = -ωp, d dt p = ωx with x, p ∈ R has solutions x(t) = αe iωt , the quantum harmonic oscillator d dt |ψ = -iω(N + I/2)|ψ features solutions |ψ(t) = |α(t) =

≥

  Finally we get Tr L L (ρ)N ν ≤ -νTr ρN ν+1 + µ Since x ν ν+1 is concave, Tr ρN ν+1 ν ν+1 Tr ρN ν . Theorem 1. Consider a trajectory [0, +∞[ t → ρ(t) ∈ K 1 (H) of the master equation (2) with = 0, κ > 0 and α > 0. The following statements hold true:

  1 (H) with Tr σN ν finite. Tr ρ(t)N ν remains bounded since d dt Tr ρN ν = κTr L L (ρ)N ν ≤ -vκTr ρN ν ν+1 ν + κµ. Thus for Tr ρN ν ≥ λ = µ ν ν ν+1 , d dt Tr ρN ν ≤ 0. The second statement exploits the first one. We can assume ν ≥ 2. Take ν > ν. Denote by K 1 ν (H) the supspace of trace-class operators σ such that Tr N ν 2 σN ν 2 is finite. The space K 1 ν (H) with the norm σ ν = Tr |σ| + Tr N

  σz on the complete system σz on the reduced system

Figure 1 :

 1 Figure 1: Comparison of Tr σ z ρ and Tr σ z ρ s for ρ solution of the complete system (2) with = 1 100 , α = 1 and vacuum initial condition (truncation up to 40 photons), and for ρ s solution of the reduced system (15) with initial condition |c + α c + α |.

Figure 2 :

 2 Figure 2: log 10 (1 -F) where F is the fidelity between ρ s and ρ for the simulations of figure 1.

  σx on the complete system σx on the reduced system

Figure 3 :

 3 Figure 3: Comparison of Tr σ x ρ and Tr σ x ρ s for ρ solution of the complete system (2) with = 1 100 , α = 1 (truncation up to 40 photons), and for ρ s solution of the reduced system (15), with the same initial condition ρ(0) = ρ s (0) = 1 2 |c + α + |c - α

  σz on the complete system σz on the reduced system

Figure 4 :

 4 Figure 4: Comparison of Tr σ z ρ and Tr σ z ρ s where ρ and ρ s correspond to simulations of figure 3.

  2 ) is a normalization factor. If |α| 1 then we have γ + ≈ γ -≈ √ 2. Throughout this paper, we denote by K 1 (H) the set of trace-class operators on H, i.e. compact Hermitian operators on H whose eigenvalues (σ k ) k∈N satisfy k≥ 0 |σ k | < +∞. This K 1 (H) equipped with the trace-norm Tr | • | is a Banach space (see e.g.

ACKNOWLEDGMENT

The authors thank Mazyar Mirrahimi for many useful discussions.

* This work was partially supported by the Projet Blanc ANR-2011-BS01-017-01 EMAQS.