Emilie Carré
email: emilie_carre@hotmail.com

Emmanuel Cantais

Olivier Darbin

Jean-Pierre Terrier

Michel Lonjon

Bruno Palmier

Jean-Jacques Risso
email: j.j.risso@imnssa.net

Technical Aspects of an Impact Acceleration Traumatic Brain Injury Rat Model with Potential Suitability for both Microdialysis and P ti O 2 Monitoring

Keywords: Traumatic Brain Injury, Microdialysis, P ti O 2, Weight-drop, Impact Acceleration, Diffuse Axonal Injury, Neuroprotection, Methodology

This report describes technical adaptations of a traumatic brain injury model -largely inspired by Marmarou -in order to monitor microdialysis data and P ti O 2 (brain tissue oxygen) before, during and after injury. We particularly focalize on our model requirements which allows us to re-create some drastic pathological characteristics experienced by severely head-injured patients: impact on a closed skull, no ventilation immediately after impact, presence of diffuse axonal injuries and secondary brain insults from systemic origin… We notably give priority to minimize anaesthesia duration in order to tend to banish any neuroprotection.

Our new model will henceforth allow a better understanding of neurochemical and biochemical alterations resulting from traumatic brain injury, using microdialysis and P ti O 2 techniques already monitored in our Intensive Care Unit. Studies on efficiency and therapeutic window of neuroprotective pharmacological molecules are now conceivable to ameliorate severe head-injury treatment.

Introduction

In order to explain the microdialysis and P ti O 2 (brain tissue oxygen) data monitored on severely head-injured patients from our Intensive Care Unit and to study the specific effects of anaesthesia and/or neuroprotection, a traumatic brain injury (TBI) animal model, similar to human head injury, is indispensable.

According to [START_REF] Bullock | Current status of neuroprotection trials for traumatic brain injury: Lessons from animal models and clinical studies[END_REF], "individual animal models rarely, if ever, model the entire spectrum of pathological characteristics observed in the patient population with severe head injuries". So the majority of drugs that have shown a neuroprotective effect on animal have usually few effects on severely head-injured patients. For this reason we have chosen to give priority to the concordance between our traumatized animals and head-injured patients characteristics.

The Marmarou impact acceleration model, commonly called "weight-drop", is a wellknown model of traumatic brain injury. This model, precisely described in two complementary articles [START_REF] Marmarou | A new model of diffuse brain injury in rats. Part I : Pathophysiology and biomechanics[END_REF][START_REF] Foda | A new model of diffuse brain injury in rats. Part II : Morphological characterization[END_REF], compiles some of the most important characteristics of human head-injury. Because of its closed skull impact, this model is more particularly in agreement with the cases of falls or road accidents, which are the most frequent situations in our Intensive Care Unit. Moreover this model was fully attested to produce diffuse axonal injury similar to that described in man -these diffuse axonal injuries are detected in more than 90% of fatal head injured patients [START_REF] Gentleman | Axonal injury : a universal consequence of fatal closed head injury[END_REF].

As far as we know, microdialysis and P ti O 2 monitoring have never been envisaged before and during impact in this type of model. So we report, in this article, some technical aspects of the TBI model adaptations to this intracerebral concomitant monitoring before, during and after impact. We particularly focalize on the requirements of our model -especially the limitation of neuroprotection -in order to recreate the most drastic conditions inherent in human head injury. 6

Materials and Methods

Ethics

The experimental protocol was approved by the local ethics committee and by the French Ministry of Defence.

Specificities of the impact acceleration TBI model

Preparation of intracerebral guides and P ti O 2 probe

Before any surgery, 2 intracerebral guides (CMA, Phymep, Paris), the first devoted to microdialysis and the second devoted to P ti O 2 , are bonded together with a specific angle, using Araldite glue (Bostik S.A., France). This angle should be determined so that sensitive areas of both microdialysis and P ti O 2 probes membranes will be very close but not directly in contact.

Because no introducer seems to be available for P ti O 2 , a modified dummy cannula is glued on the P ti O 2 probe so that the probe will perfectly fit to the guide.

Surgery

Two weeks before the experimentation, a male Sprague-Dawley rat weighting between 400 and 450g (OFA strain, Iffa Credo, France) is anaesthetised with sodium pentobarbital (60mg/kg intraperitoneally) before surgery.

The intracerebral microdialysis guide -bonded to the P ti O 2 guide -is stereotaxically implanted in the striatum, 4 mm higher than the precise location of the microdialysis probe membrane, according to the atlas of [START_REF] Paxinos | The rat brain in stereotaxic coordinates[END_REF]: coordinates relative to Lambda in mm A:+9.7; L:+3; H:+6.5. The two guides, and a support for fixing to the swivel, are fixed on the very anterior part of the skull with two screws and dental cement (Dentalon plus, Heraeus Kulzer, Germany). Immediately back, the 10mm-diameter impact site on the skull must absolutely stay clear of cement. The skin above this impact site must be closed by stitches in order to be preserved from air contact.

After surgery, the rat is allowed to recover in individual cage, under a 12-hours light/dark cycle (light on from 07h00 to 19h00), with free access to food and water for the two weeks before the experiment.

Microdialysis and P ti O 2 monitoring

Twelve hours prior to the experiment, the rat is placed in the experimental plexiglass cage in order to adapt to his new environment.

The day of the injury, the dummy cannula of the straight guide is replaced by the microdialysis probe (CMA12, 4mm-length, 0.5mm-diameter, Phymep, France). This probe is perfused by an artificial CSF (in mM ; NaCl:147 ; KCl:2.7 ; CaCl 2 :1.2 ; MgCl 2 :0.85) at the rate of 1µL/min. The modified P ti O 2 probe (LICOX CC1.R, 4mmsensitive area length, 0.5mm-diameter, Integra NeuroSciences, France) is inserted in the slope guide and connected to the LICOX CMP instrument (Integra NeuroSciences, France).

Because rat needs to be in freely-moving, a microdialysis liquid swivel has been especially adapted in our lab in order to allow concomitant passing of both electrical P ti O 2 signal and liquid microdialysis samples.

A first phase of parameters stabilisation (at least 2 hours) should be observed. Then basal levels of PtiO 2 and microdialysis parameters can be monitored.

Induction of Traumatic Brain Injury

Preliminary experiments lead us to induce Traumatic Brain Injury under transient 3%isoflurane anaesthesia (less than 10 min).

Stitches above site impact are removed and a 10mm-diameter 3mm-thick metallic disk -designed to protect against skull fracture -is placed directly in contact with the clear part of the skull. A cylindrical metallic 430g-weight is dropped from two meters through a metal tube onto the disk (Fig. 1). Weight rebound is prevented using an automated system especially designed in our lab.

A velocity sensor, designed in our lab as well, certifies reproducibility of the weight acceleration at the very moment of the impact.

After impact, the rat is allowed to return to his plexiglass cage, in order to recover from anaesthesia and to continue the recordings in conscious freely-moving conditions.

Evidences of traumatic brain injury

Behavioural test

Four days after injury, an adhesive-removal somatosensory test, modified from [START_REF] Schallert | CNS plasticity and assessment of forelimb sensorimotor outcome in unilateral rat models of stroke, cortical ablation, parkinsonism and spinal cord injury[END_REF], was used to reveal evidences of TBI-induced behavioural lesions.

A standardised adhesive stimulus was attached to one of the rat forelimb. Rats removed the stimulus using their teeth. The latency of stimulus contact and removal was recorded and compared between 6 injured and 6 sham-lesioned rats (non parametric C1 Fisher-Yates-Terry test).

Brain Fixation and Histopathological Preparation

At least 24h after impact, rats were deeply anaesthetised with an intraperitoneal injection of sodium pentobarbital. The chest was rapidly opened, a catheter was introduced into the ascending aorta, and the right atrium was incised. Firstly, 200 mL of heparin saline (1000U.I. of heparin in saline) were perfused through the catheter, at a rate of 25mL/min. Secondly, 400 to 500mL of fixative (4% formaldehyde, 3% acid acetic, in saline) were perfused at the same rate. The brain was carefully removed and stored in a fixative (4% formaldehyde in saline) for at least 24 hours. The brain was finally placed in 10%-formaldehyde for 12 hours before gross examination. Brain coronal and sagittal sections were embedded in paraffin. Sections 5 µm-thick were cut

Results and Discussion

Surgery

Microdialysis and P ti O 2 guides -included in cement -take up a large space on the rat head (Fig. 2). On account of the little place left, the impact site is quite closed to the cement. Despite we wanted to collect microdialysis data and P ti O 2 really nearby the impact site, such high closeness was firstly difficult to use: at the moment of the impact, the cement was sometimes pulled out from the skull because of the force of the impact.

In order to limit the incidence of such event, manual dexterity -when strongly interlocking skull and cement -must be carefully optimized, as one of the crucial point of our model. During the 2 weeks delay between implantation and experiment, the impact site behind the cement must be preserved of air contact in order to not weaken the bone structure before impact. So the first rats with an impact site skull exposed to air systematically experienced a severe skull fracture at the impact site. The use of stitches to close the skin above the impact area suppresses the incidence of skull fracture. Moreover, Piper reports important sliding frictions in the initial Marmarou impact acceleration model, probably because of metal/plexiglass contact. But our results show few differences between theoretic velocity and recorded velocity (respectively 6.26m/s vs. 6.06m/s, i.e. less than 3.3%), that is to say virtually no sliding frictions. A good suitability between the 2m-tube and the weight shapes may explain these limited frictions, but the main reason is that our tube is no longer made with plexiglass but with metal, and metal/metal contact is said to induce only few sliding frictions.

Reproducibility of the impact

In the initial Marmarou impact acceleration model, after the impact, the rat foam bed is commonly gently pushed away in a lateral direction to prevent a second impact. This is particularly hazardous because a good coordination between weight dropping and foam bed moving is indispensable, otherwise the rat could be one more time injured. This situation is absolutely inconceivable in our model: various microdialysis and P ti O 2 connections take up so many space on the rat head that a so fast moving is impossible without a risk of connections damage. Our automated system preventing weightrebound allows us to prevent from such problems because weight is caught when moving up in the metal tube after the impact. No weight-deceleration or problem of reproducibility were noticed using this automated system.

Evidences of traumatic brain injury.

Modifications from the initial Marmarou model did not seem to induce any alterations of the characteristics of the traumatic brain injury pattern.

For example, in our model, immediately after injury, all rats suffered from respiratory distress such as apneas. This is in accordance with [START_REF] Marmarou | A new model of diffuse brain injury in rats. Part I : Pathophysiology and biomechanics[END_REF], who

observed, on all rats surviving the 450g-2m impact, apneas immediately after injury and a reduction in respiratory rate of 20% for up to 30 min postinjury. After a period of 20 to 30 min, we no longer observed rat death. This is again in accordance with Marmarou who noted that, after 30 min postinjury, respiration in his animals gradually recovered and was not significantly different from his control rats by 2 hours postinjury.

In accordance with [START_REF] Foda | A new model of diffuse brain injury in rats. Part II : Morphological characterization[END_REF], our injured-rat recovery from anaesthesia was delayed for several minutes (around 5 min for injured rats vs. <1 min for sham-lesioned rats). [START_REF] Schallert | CNS plasticity and assessment of forelimb sensorimotor outcome in unilateral rat models of stroke, cortical ablation, parkinsonism and spinal cord injury[END_REF] described a bilateral somatosensory test that is useful in studies of loss and/or recovery of function following Central Nervous System injury. We simplified this test in an unilateral somatosensory test (removal of an adhesive label attached only to one rat forelimb) in order to detect behavioural TBI-induced alterations. Injured rats were actually behaviourally altered 4 days after injury: the latency of adhesive-removal was drastically increased for injured rats (1h 6min vs. 8min

for sham lesioned rats, alpha<0.01, non parametric C1 Fisher-Yates-Terry test). [START_REF] Foda | A new model of diffuse brain injury in rats. Part II : Morphological characterization[END_REF] observed in his model that massive diffuse axonal swelling and other forms of brain oedema reach a maximum after 24 hours. He also noted that diffuse axonal injuries, in the form of retraction balls, continue to be visible, although smaller in size and number, until the 10 th day after trauma. That is why we decided to realize our histopathological preparations at least 24 hours after the impact. We observed on our cerebral histological sections some of Foda's histological characteristics (1994): pink shrunken neurons, considered as a sign of neuronal death, and diffuse axonal injuries in the form of retraction balls (Fig. 3). The pink shrunken neurons were surrounded by perineuronal vacuolations. In association with these injured neurons, we noted pericapillary brain oedema, and several capillaries were congested with red blood cells. underline that recovery from anaesthesia must be as short as possible. We also noted for control rats that longer is the isoflurane anaesthesia period, longer is the recovery from anaesthesia ; this situation is aggravated for injured rats. This also implies that the anaesthesia duration must be minimized, in order to reduce the risk of mortality incidence.

A

The choice of shortest transient isoflurane anaesthesia was also dictated by the fact that we wanted to banish any neuroprotection in our model. Actually, isoflurane-related neuroprotection is still a matter for debate: although isoflurane anaesthesia was commonly said to protect CNS tissue against cerebral injuries such as trauma or ischemia, some authors believe that isoflurane anaesthesia only delays but does not prevent cerebral infarction in rats subjected to focal ischemia [START_REF] Warner | Isoflurane neuroprotection: a passing fantasy, again ?[END_REF][START_REF] Kawaguchi | Isoflurane delays but does not prevent cerebral infarction in rats subjected to focal ischemia[END_REF]. Whether isoflurane can induce neuroprotective effects or not, we must anaesthetize our animals according to the respect for Bioethics. For this reason, we have chosen, as a precaution, to limit isoflurane anaesthesia duration in order to minimize neuroprotective effects. Moreover we have observed on control rats that 3%isoflurane anaesthesia (same results for halothane anaesthesia) for more than 15 to 20 min may result in metabolic alterations up to the 30 hours following anaesthesia (unpublished data). Because our model was adapted in order to notably study metabolic and neurochemical alterations resulting from traumatic brain injury, such anaesthesiainduced alterations might be detrimental to our monitoring. All these observations led us to use isoflurane only transiently and less than 5 to 10 min (necessary time to induce anaesthesia and to induce traumatic brain injury) in order to avoid both rat death and metabolic alterations, as well as to limit neuroprotection. [START_REF] Marmarou | A new model of diffuse brain injury in rats. Part I : Pathophysiology and biomechanics[END_REF] , 1999 ;Brain Trauma Foundation, 1996), mechanical ventilation is a neuroprotective treatment notably against secondary cerebral hypoxia related to systemic hypoxemia. So, when ventilating, the rats might be both protected against secondary-induced hypoxia and against some form of diffuse axonal injuries. That is why we have chosen to banish ventilation in our model in order to stay in the most drastic conditions, in agreement with human head-injury.

Fifteen seconds after the impact, [START_REF] Marmarou | A new model of diffuse brain injury in rats. Part I : Pathophysiology and biomechanics[END_REF] noticed an increase of blood pressure, followed by a period of hypotension, which is a typical situation observed in head-injured patients (French Society for Anaesthesia and Intensive Therapy editors, 1999 ;Brain Trauma Foundation, 1996). Such periods of hypotension are commonly said to induce secondary brain insults from systemic origin -particularly cerebral ischemia -because of a cerebral blood flow decrease. [START_REF] Bullock | Current status of neuroprotection trials for traumatic brain injury: Lessons from animal models and clinical studies[END_REF] criticize most TBI animal models essentially because they do not include the secondary insults that are typically observed in severe human head injuries, such as hypoxia or hypotension. In order to study the true alterations resulting from traumatic brain injury, we needed such model with inherent secondary brain insults, which is closer to clinical head-injured patients situation immediately after injury.

Conclusion

Some of the most drastic conditions experienced by severely head-injured patients (i.e.

  Fig.1: Schematic representation of the impact acceleration device. metal 2m-tube 430g-weight foam bed cement and impact site velocity sensor system preventing weight-rebound (Figure freely adapted fromMarmarou, 1994)

Fig. 2 :

 2 Fig.2: Schematic representation of the intracerebral probes and of the impact site microdialysis guide and probe P ti O 2 guide and probe dental screw dental cement support for fixing to swivel skull 10mm-diameter metallic disk : impact site

  [START_REF] Piper | Monitoring weight drop velocity and foam stiffness as an aid to quality control of a rodent model of impact acceleration neurotrauma[END_REF] affirmed, about the initial Marmarou impact acceleration model, that the weight-drop velocity of a 450g-weight dropped from two meters through a plexiglass tube can vary by as much as 40%. The velocity control tests performed on repeated dropping in our model show only less than 1% variation of velocity (mean velocity 6.06 m/s, SEM 0.0079 m/s), which means a perfect reproducibility of our model.

Fig. 3 :

 3 Fig.3: Photomicrographs of brain sections of an injured rat, 72 hours after injury. A (HES X100): Retraction balls and many pink shrunken neurons can be observed. B (HES X400): Capillaries are congested with red blood cells. Pink shrunken neurons are surrounded by perineuronal vacuolation.

  impact on a closed skull, absence of ventilation immediately after impact, presence of diffuse axonal injuries and of cerebral secondary brain insults from systemic origin…) have been successfully reproduced in our model. Because this model allows us to monitor both microdialysis and P ti O 2 parameters -more and more commonly monitored in Intensive Care Unit -ongoing experiments concerning post-traumatic neurochemical and biochemical mechanisms are carried out in our Lab. Moreover studies focused on influence of neuroprotective pharmacological molecules are now conceivable to ameliorate severely head-injured patients treatment.

  observed that non-ventilated rats that did not survive experienced apnea lasting for up to 20 seconds immediately after impact, with a gradual slowing of respiration until death. He related his high mortality rate in spontaneously breathing animals to a transient central respiratory dysfunction that was readily reversible by respiratory support: 8.7% of mortality on intubated mechanically ventilated rats vs. 58.6% on non-intubated non-ventilated. Moreover he observed that the severity of posttraumatic clinical observations and pathological changes was similar in ventilated or non-ventilated animals, which mainly means that diffuse axonal injuries are beyond the use of mechanical ventilation. Nevertheless diffuse axonal injuries, typically TBI-induced, are sometimes said to be aggravated by ischemia and/or hypoxia

(cited by

Povlishock and Chritman, 1995)

. In accordance with the majority of guidelines on severely head-injured patients management (French Society for Anaesthesia and Intensive Therapy editors

Acknowledgements:

This work was supported by a DGA Grant (French Ministry of Defence): PEA n°000805.

Special thanks to Jean-Michel Salon, Pascal Schirck, Bruno Schmid and Boualem Zouani for mechanical conceptions and to Josy Magréau for histology.

with a rotary microtome, stained with haematoxylin-eosin-safran (HES), and examined under light microscopy.

Anaesthesia and neuroprotection

In order to select, in our model, the most appropriate type of anaesthesia during the impact, preliminary experiments were performed: 2 protocols of anaesthesia were compared. In the first group, 30 rats were submitted to the protocol of Traumatic Brain Injury (as described in the chapter "Material and Methods -Induction of Traumatic Brain Injury"). In the second group, 10 rats were submitted to the same protocol except for the 3%-isoflurane anaesthesia which was replaced by an intraperitoneal injection of sodium pentobarbital (60mg/kg) 30 minutes before the impact . Our results show an appreciable difference in the mortality of these injured rats according to the way they were anaesthetized during the impact: 100% of death under sodium pentobarbital anaesthesia (injection) vs. only 33% under transient isoflurane anaesthesia (inhalation). [START_REF] Marmarou | A new model of diffuse brain injury in rats. Part I : Pathophysiology and biomechanics[END_REF] did not observe such differences between his 2 types of anaesthesia: alpha-chloralose injection or isoflurane inhalation (respectively 44% of mortality vs. 58.6%, for a 450g/2m impact). In fact these surprising results should not be analysed according to the type of anaesthesia (injection or inhalation), but should be related to the anaesthesia duration (long-term/general or transient anaesthesia). Like this, two situations may be considered: firstly a high mortality under general long-term anaesthesia (100% in our model, and 44 to 58.6% for Marmarou), and a more moderate mortality under transient anaesthesia (33% in our model). This beneficial effect of transient anaesthesia (less than 10 min in our model) might be related to the quick recovery from transient anaesthesia because, as soon as they awake, they no longer suffer from apnea, maybe because they control again their breath. All transientlyanaesthetized rats that did not survive died before they awaked. Moreover, when rats were stimulated to awake, the corresponding mortality decreased. These results


Luis Basora
email: luis.basora@onera.fr

Jérôme Morio
email: jerome.morio@onera.fr

Corentin Mailhot
email: corentin.mailhot@aviation-civile.gouv.fr

A Trajectory Clustering Framework to Analyse Air Traffic Flows

This paper describes a framework to automatically identify air traffic flows from a set of trajectories by using a clustering algorithm. The framework offers two methods to cluster trajectories, each one using a different distance/similarity measure between trajectories. Results and performance characteristics of both methods are compared by applying them to real trajectories over a French Area Control Center. The framework can output statistics and figures for flow analysis and its use is facilitated by the relatively low number of parameters to be provided by the user. Its aim is to help support the SESAR vision of flow-centric operations by being integrated into Air Traffic Management tools, e.g. for airspace design/management or for analysis of traffic patterns in a free route environment.

I. INTRODUCTION

In Air Traffic Management (ATM), massive amounts of data are available to feed data-driven models for real-time decision support or to identify behaviours or patterns relevant to the operational performance of the system in post-event analysis. In particular, sources such as radars, Automatic Dependent Surveillance-Broadcast (ADS-B) or trajectory prediction models generate samples of data points representing trajectories that can then be clustered to identify traffic flows. Trajectory clustering algorithms could be useful in the context of the flow-centric operations vision of SESAR, as described in the European ATM Master Plan [START_REF] Undertaking | European ATM master plan[END_REF], by being integrated into tools supporting airspace design/management, complexity management, etc. This is particularly true in a free route environment, where the capacity to understand traffic flows is even more necessary as fixed routes will no longer structure the traffic.

Clustering is a widely used data analysis technique in the statistics/machine learning field. It is about grouping entities with similar characteristics together, so the notion of similarity/distance is essential to the problem. A number of clustering algorithms have been reported in the literature, such as kmeans [START_REF] Lloyd | Least squares quantization in pcm[END_REF], BIRCH [START_REF] Zhang | Birch: an efficient data clustering method for very large databases[END_REF], DBSCAN [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] and OPTICS [START_REF] Ankerst | Optics: ordering points to identify the clustering structure[END_REF], all of which are oriented towards the clustering of point data. Even though trajectories have a functional nature (curves), these algorithms can still be applied since trajectory data is available in the form of samples of data points.

For instance, in a recent paper [START_REF] Murc ¸a | Trajectory clustering and classification for characterization of air traffic flows[END_REF], DBSCAN shows promising results when applied to the characterisation of traffic flows based on recorded radar tracks in the terminal airspace of the New York Metro region. This algorithm has the capability to handle noise/outlier data and does not require the number of clusters to be provided as an input parameter. Previously, in [START_REF] Gariel | Trajectory clustering and an application to airspace monitoring[END_REF], another framework is described based on DBSCAN and k-means to analyse the patterns of traffic over the Northern California terminal area. In both of these studies, however, no distance/similarity measures between trajectories is provided, i.e. the clustering is based solely on the density of the individual trajectory points.

Alternatively, in [START_REF] Lee | Trajectory clustering: a partitionand-group framework[END_REF] and [START_REF] Nguyen | Identification of air traffic flow segments via incremental deterministic annealing clustering[END_REF], two different methods are presented to cluster trajectory segments rather than trajectory points. The former (TRACLUS) implements a variant of DB-SCAN to cluster the segments but has never been successfully applied to air traffic as far as we know. The latter is designed to identify air traffic flows in the National Airspace System (NAS), but it is limited to 2D and based on the development of specific algorithms for incremental clustering requiring a non-obvious parameter setting by the user.

Another interesting study that examines the problem of clustering air traffic trajectories is reported in [START_REF] Enriquez | Identifying temporally persistent flows in the terminal airspace via spectral clustering[END_REF] where a spectral clustering approach is used to consider the temporal characteristics of the air traffic flows in the US. However, the procedure to take into account the shape of the trajectories is not evident and the curvature of the clustering centroids might not always be very realistic from an operational point of view. Other approaches such as in [START_REF] El Mahrsi | Graph-based approaches to clustering network-constrained trajectory data[END_REF] relies on a graph structure like a road network, which is not adapted to ATM because of the direct routes often allocated by ATC for tactical reasons or when applied to a free route environment.

More recently, in order to overcome the limitations of some of these techniques, an approach based on entropy minimisation and Lie group modelling has been proposed in [START_REF] Puechmorel | Entropy minimizing curves with application to flight path design and clustering[END_REF]. Unfortunately, only partial results are available at the moment of writing this paper and the implementation of the algorithms is fairly complex.

In this paper, we propose a new framework for air traffic flow analysis based on an improved version of DBSCAN called Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN) [START_REF] Campello | Density-based clustering based on hierarchical density estimates[END_REF] which is able to manage clusters of different densities with a single input parameter. Two methods based on two different distance functions between trajectories, Euclidean Distance (ED) and Symmetrized Segment-Path Distance (SSPD) [START_REF] Guillouet | Apprentissage statistique: application au trafic routier à partir de données structurées et aux données massives[END_REF] can be selected by the user. These distances offer a trade-off between accuracy and runtime. In addition to the clustering capability, the framework provides several components to filter and interpolate trajectories, compute basic flow statistics and export clusters in Keyhole Markup Language (KML). The framework can be useful to analyse traffic patterns in a wide range of operational scenarios both in en-route and approach. In addition, its implementation and integration into ATM applications is facilitated by the public availability in several languages of some of the algorithms like HDBSCAN [START_REF] Mcinnes | hdbscan: Hierarchical density based clustering[END_REF]. This paper is organised as follows. Section II details the methodology and the different components of the framework. Section III presents the results and the comparison of the two methods by using DDR2 [16] trajectory data provided by Eurocontrol over the Reims Area Control Center (ACC). Finally, section IV summarizes and identifies some ideas for future work and potential applications.

II. THE FRAMEWORK

The framework is designed to satisfy the following requirements:

1) It shall be possible to filter trajectories geographically, by airspace, by altitude and in time.

2) It shall be possible to cluster in 2 or 3 dimensions, optionally taking the heading into account as well.

3) The clustering algorithm shall allow for variable tradeoff between computing time and quality of the resulting clusters so as to be useful in scenarios with either large or small sets of trajectories. 4) The clustering algorithm shall work within an area where traffic density distribution is not uniform. 5) The clustering algorithm shall work with noise data and be able to identify outliers. 6) The clustering algorithm should require a reduced number of input parameters with clear ATM operational significance. 7) For each cluster it shall be possible to compute the centroid in addition to the set of associated trajectories. 8) For each cluster it shall be possible to compute statistics such as the flow rate, flight distribution per origin/destination, average distance and heading of the cluster trajectories. Figure 1 shows the framework architecture with the two clustering methods in the center and the pre-processing and post-processing steps. The user can choose the method and the associated parameters. The clustering method based on ED works faster than SSPD at the possible expense of a lower quality clustering result. The specific steps of the framework depend on the selected distance and are described in the following sub-sections.

The framework is implemented in Python 3 with the following Python libraries: scikit-learn 0.18.2, RDP 0.8, Hdbscan 0.8.8, Pyproj 1.9.5.1, Planar 0.4. The HDBSCAN library is highly optimised, but otherwise no parallelism, multi-threading or code optimisation has been used so far in areas such the SSPD implementation where considerable gains of performance can be expected.

A. Preprocessing

Trajectory datasets contain 3D position samples (latitude, longitude, altitude) for each trajectory and possibly other information like speed and heading. Latitude/longitude coordinates are projected in order to facilitate distance computation. Speed is not used by the clustering algorithms and heading is optional in the ED-based clustering method to ensure flows in opposite directions are separated into different clusters. If heading information h is required but not part of the dataset, it is calculated for each trajectory segment from geographic coordinates [(lat1, lon1), (lat2, lon2)] and ∆lon = lon2 -lon1 by: h = (arctan(a) + 360) % 360

(1) where a = cos(lat2) sin(∆lon) cos(lat1) sin(lat2) -sin(lat1) cos(lat2) cos(∆lon)

It is important to note that position samples can be irregularly distributed along the trajectory as shown in Figure 2. Many data points can be removed as they are redundant, whereas other positions need to be interpolated.

We use the Ramer-Douglas-Peucker (RDP) [START_REF] Ramer | An iterative procedure for the polygonal approximation of plane curves[END_REF] [18] algorithm to remove redundant trajectory information, but only in the case of SSPD where considerable runtime gains can be expected. Further details on the application of this algorithm are given in section II-B2.

Regarding the interpolation method, there are several methods available from simple linear interpolation to more sophisticated polynomial or spline interpolation. In the framework, the Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) [START_REF] Fritsch | Monotone piecewise cubic interpolation[END_REF] [20] method is used as it better takes into account the operational reality of air traffic trajectories. As shown in the example of the vertical profile of a flight in Figure 3, linear interpolation is too abrupt, whereas the order 3 spline ensures the continuity but with oscillations (overshoot). With PCHIP, however, we get a smoother curve much more representative of a true flight trajectory.

Finally, as for the filtering capability, in addition to a an altitude range or a time window, the framework allows the user to chose either a bounding box around the geographic area to be studied or a set of airspaces (e.g. ACC or sector). In the latter case, the convex hull is computed for the set of airspaces.

B. Clustering

Once the preprocessing phase is complete and trajectories have been properly prepared, we can apply the HDBSCAN algorithm to them to identify the traffic flows. Compared to DBSCAN, HDBSCAN presents several improvements. First, it only requires the user to define the minimum cluster size, i.e. the minimum number of points (trajectories in our case) to form a cluster. Secondly, HDBSCAN works better than DBSCAN for data with varying density, which is usually the case for air traffic trajectories. Thirdly, HDBSCAN prevents the "bridge effect" between two clusters because of a single or a few data points in the middle of the two clusters (potentially gluing them together) by considering these points as noise.

The choice of an appropriate distance/similarity measure between trajectories is as important as the choice of the clustering algorithm. There are quite a few functions to compute trajectory distances reported in the literature, each one offering a different trade-off between computation time, accuracy and sensitivity to noise data, from the simple ED to the more sophisticated and complex ones like Fréchet and Hausdorff. As our framework is to be used in a wide range of use case scenarios, it cannot be based on a single metric. Also, we need a similarity measure that accounts for both the physical distance and the global shape between trajectories, i.e. time or other trajectory characteristics will not be considered.

We decided to choose the ED function since it is computationally fast and so potentially useful for scenarios involving a large number of trajectories, such as the ones covering wide geographic areas or large time windows. For scenarios where accuracy is more important and execution time less of an issue, SSPD distance should be used instead.

HDBSCAN labels each trajectory with a cluster number (outliers are labelled with a -1) and measures the strength of the cluster membership for each trajectory in the cluster with a probability. In this framework, a cluster is defined by three attributes: a cluster number, the list of trajectories belonging to this cluster (with probability equal to 1), and a representative trajectory or centroid (this last attribute is computed after the execution of HDBSCAN).

In the following subsections, the specificities of the application of HDBSCAN with each distance are detailed.

1) ED-based Clustering: The first step in this method is to build a matrix where each row represents a whole trajectory built from the sequence of trajectory points. This requires all the trajectories to have exactly the same number of points, which is indeed the case as we interpolate the trajectories according to the number of points chosen by the user. This is done at least once in the pre-processing phase before filtering, and again if the trajectories are clipped as a result of a filter being applied.

For instance, in the case of a 3D clustering with heading, N trajectories and n points per trajectory, the matrix would be defined as:   

x 11 y 11 z 11 h 11 x 12 . . . x 1n y 1n z 1n h 1n x 1n . . . . . . . . .

x N 1 y N 1 z N 1 h N 1 x N 2 . . . x N n y N n z N n h N n x N n    (2)
where x ij y ij is the projected position, z ij the altitude (in meters) and h ij the heading computed as in [START_REF] Undertaking | European ATM master plan[END_REF]. In order for these elements to have the same weight, we standardise each value in the matrix by subtracting the mean and dividing by the variance.

Once these steps are concluded, we can apply the HDB-SCAN algorithm with the following two parameters: the minimum number of trajectories per cluster (user parameter) and the clustering matrix built as in [START_REF] Lloyd | Least squares quantization in pcm[END_REF]. ED is computed from the clustering matrix by the HDBSCAN itself as it is part of the standard distances already implemented in the algorithm.

The last step is to calculate the centroid of each cluster as the mean of the trajectories of the cluster.

2) SSPD-based Clustering: From a mathematical point of view, SSPD should be considered as a similarity measure rather than a distance or metric, because although it is symmetric, the triangle inequality property is not satisfied. Compared to the Fréchet distance, SSPD does not consider the speed of the trajectory since it is purely geometrical. On the other hand, if two trajectories are similar during the en-route phase because they share the same air routes and diverge only in the final phase of the flight at the terminal area, distance would be over-estimated by both Hausdorff and Fréchet. This is not the case with SSPD, because it better takes into account the global difference in the shape of both trajectories. Also, it presents a good trade-off between the simpler and faster ED and the more complex and time-consuming Hausdorff and Fréchet distances. However, with the SSPD distance it may be difficult to separate trajectories that are geographically close, similar in shape and length but having opposite directions. This may be less of an issue when the clustering is performed in 3D, as air routes for instance are designed to vertically separate flows in opposite directions.

After the pre-processing phase where trajectories are projected, interpolated and finally filtered for the area of interest, there is no need to interpolate the filtered trajectories as opposed to in the ED case. This is because SSPD does not require the trajectories to have the same number of points. However, the SSPD distance being much slower to compute than the ED, the RDP algorithm is applied first to remove redundant trajectory points.

The RDP algorithm can be particularly effective when applied to the en-route phase of the trajectory where aircraft follow long great circle segments linking the flight plan waypoints. Nevertheless, even during the cruise phase, aircraft do turn to follow jet routes and take advantage of favourable winds, avoid hotspot areas, etc. Therefore, we need to use this algorithm carefully enough so as not to exclude too many (if any) turning points. RDP accepts a threshold parameter allowing for the simplification, to a greater or lesser degree, of the trajectory, which can be specified by the user as a framework parameter. With a threshold of 200 meters we notice already an important reduction of redundant points for some trajectories like in Figure 4.

Once RDP has been applied to the trajectories and the distance matrix is computed, we can run HDBSCAN to perform the clustering. Afterwards, we finally obtain the cluster centroids by choosing the trajectory in the cluster that minimises the distance with the other trajectories in the cluster. Therefore, with this method, centroids are true trajectories and not virtual trajectories like with ED.

C. Post-processing

The objective of this phase is to generate from the clustering results all the necessary outputs to enable the user to analyse the structure of the traffic represented by the clusters. In particular, the framework can generate:

1) 3D and 2D plots of the centroids indicating visually the flow intensity and direction. 2) Histograms to show the distribution of trajectories per cluster. 3) A set of statistics per cluster: The objective is to identify the flows over Reims ACC and for that we need to specify the minimum number of trajectories (minimum cluster size) for a flow to be considered as such. This is the main clustering parameter to be defined by the user, and its value is highly application-dependent. Our purpose being purely the assessment of the algorithms, we set it up so that major flows can be identified.

The computer used in the evaluations is an Intel R Xeon Quad-Core 2.80GHz processor with 6GB of memory.

A. ED-based clustering results

This first experiment covers the complete volume of traffic in Reims ACC and is performed with the parameters specified in the first column of Table I.

We want to identify flows of at least 50 trajectories. Even though increasing the number of points after interpolation may improve the clustering accuracy, we have found 100 points to be a reasonable value. We also filter out any trajectory shorter than 20 NM, which reduces the number of trajectories to 9307.

PCHIP interpolation has the longest computation time, taking 99 seconds against 36 seconds for HDBSCAN. In total, outputs are generated in about 3 minutes, which is acceptable for a scenario with more than 9000 trajectories.

Figure 5 shows the 2D and 3D views of the 38 cluster centroids, where the thickness of the centroid lines is proportional to the intensity of the flow as indicated in the flight distribution per cluster figure. The algorithm has identified both en-route and terminal flows, the latter ones mainly around Paris Charles-de-Gaulle or Paris-Orly (Paris coordinates are 49 • , 2.3 • ) terminal area. We can use the KML output in Figure 6 to better visualise the centroids for the flows departing from and approaching Paris.

The number of outliers (32%) is high, but it depends on the size of the cluster (in our case 50) as shown in Figure 7. Also, it can be explained by the fact that the chosen day was the one with the largest volume of traffic for 2015 in Reims. This particularly high traffic density may have induced a higher complexity and required exceptional measures to resolve conflicts or hotspots.

B. SSPD-based clustering results

In this case, we use the SSPD-based clustering method to analyse the flows in the LFEEXR sector, with the parameters in the second column of Table I.

After filtering out the flights not entering this sector, a total number of 349 trajectories remain. The RDP algorithm, with a threshold of 200 meters, is next applied to remove redundant trajectory points. Then the SSPD distance is calculated, taking 44 seconds. Adding the 98 seconds of interpolation time for the initial set of 9442 trajectories, the total computation time is approx. 2 minutes and 30 seconds, which is relatively high compared with the 3 minutes for the over 9000 trajectories of the first experiment.

Five clusters are identified with 25% of outliers. We can see the five centroids and the number of flights per cluster in Seventh SESAR Innovation Days, 28 th -30 th November 2017 8. It has to be noted that with SSPD a cluster can have trajectories with opposite directions because of the alreadymentioned limitation of this distance and the fact that the clustering is performed in 2D. Therefore, the direction displayed for the cluster centroids may be the opposite direction of some of the trajectories in the cluster.

Longitude (°) Latitude (°) 2°W 2°W 1°W 0°1°E 2°E 3°E 4°E 5°E 6°E 7°E 8°E
In Figure 9, we superimpose the identified clusters centroids in red over the air routes in green to make sure that both are consistent in terms at least of geographical location. For instance, cluster number 2 is a northwards flow located in the east of Paris matching perfectly with air route UM733.

Longitude (°) Latitude (°) 0°1°E 2°E 3°E 4°E 5°E 6°E 7°E 47°N 48°N 49°N 50°N 51°N 52°N 52°N 0 1 2 3 4 5 centroids 0 1 2 3 4
Cluster number

C. ED versus SSPD comparison

The last experiment consists of comparing both clustering methods by applying them to the same use case. The goal is to identify the major flows in 2D (minimum of 50 flights per day) in Reims ACC between FL300 and FL500. Clustering is performed in 2D because it is more convenient for visual verification, as controllers are used to a 2D view of the traffic. II.

Parameters are shown in the third column of Table I and the results in Table

There is a significant difference in time performance between the two methods, with ED-based clustering being almost 300 times faster than SSPD-based clustering (about 5 minutes versus more than 24 hours). Additionally, ED-based clustering produces about 10% less outliers, which can be explained by the fact that SSPD better takes into account the differences in shape and physical distance between trajectories.

As for the quality of the clustering, ED performs worse in some cases like in cluster number 8 (see Figure 10). Even if the trajectories in this cluster are all southbound, it would be more appropriate to separate them into different clusters to better account for the diversity of headings. However, the ED distance is not sensitive enough to operate this separation and the only alternative would have been to decrease the minimum number of trajectories per cluster, which would result in an even a greater number of clusters. With SSPD, we obtain fewer and more homogeneous clusters, but the number of outliers (43 %) is considerable.

In order to further analyse the identified clusters, Table III gives the statistics for two of the flows identified by the SSPD method. Thus, 5% of the flights in cluster 7 follow the Paris Charles-de-Gaulle (LFPG) -Heathrow airport (EGLL) route, whereas 39% of the flights in cluster 8 follow the Nice (LFMN) -Paris Orly (LFPO) route. It may not at first be evident as to why cluster 8 overlaps cluster 7, rather than being merged into one single cluster. In fact both clusters contain a majority of northbound flights from south-east France and Italy. However, most of the flights in cluster 8 are for Paris (89% LFPO), whereas in cluster 7 the main destination is London-Gatwick (17% EGKK), therefore there are two distinct routes that are represented by these two clusters.

Figure 11 displays the centroids and their intensities (number of trajectories) to compare the results of both methods. It is not difficult to match the major flows, e.g. ED flows 3 (green), 10 (pale yellow), 1 (blue) correspond to SSPD flows 4 (salmon), 7 (orange), 5 (red). However, except for cases like the anomalous cluster ED 8 and the corresponding SSPD 14, SSPD clusters have a higher number of trajectories as this distance is unable to discriminate similar trajectories with opposite directions. For instance, we can observe that clusters 17 (red), 23 (brown), 25 (blue) identified by the ED method have been replaced by the single and bigger SSPD cluster 19 (dark orange).

D. Verification with planned trajectories

All of the previous experiments are based on the executed trajectories only (M3 in DDR2). However, the difference between the planned and executed traffic may be significant for days with an important volume of traffic as in our case, so it may be interesting to check how well the clustering With the executed traffic, we have already seen that some of the clusters correspond to air routes in a previous experiment (see Figure 9). In the case of the planned trajectories, which are based on ARN and have not been modified by ATC in the tactical phase, the matching should be even more evident if our clusters are consistent with the operational reality. We use the same parameters as those for ED (first column in Table I), i.e. with a minimum number of trajectories equal to 50. Even with the ED method, some of the obtained clusters are quite accurate e.g. clusters 4 and 20 in Figure 12.

These two clusters have 74 and 84 trajectories respectively. In each cluster, all trajectories are almost completely overlapping since no deviations from the flight plans were introduced by ATC. In Figure 13 it can be observed that sections of the published air routes UM133 and UM728 over Reims ACC (in green) match perfectly with the two clusters (in black). More generally, most of the clusters correspond to sections or links between sections of air routes. It is also interesting to note that even though the two routes are geographically close, the clustering is accurate enough to properly identify two separate flows.

IV. CONCLUSION

In this article, we have presented two methods based on the HDBSCAN clustering algorithm to identify air traffic flows from a set of trajectories. The choice of the method to use depends on the volume of trajectories to be processed and the desired accuracy. For a precise clustering, over a limited area, such a sector or small family of sectors, SSPD-based clustering is the best-adapted, except if we want to make sure that similar 2D flows in opposite directions are properly separated. Otherwise, ED-based clustering is probably the best option as SSPD may be prohibitively slow for some applications unless parallel computing for SSPD is implemented.

Our results have shown that we can match the obtained clusters to existing operational air routes. For the planned trajectories, clusters fit the route structure over Reims as expected, which reinforces our confidence in the framework results. On the other hand, outliers are in general quite high and further analysis is needed to characterise them. The best way to validate the framework would be to use it in a real application and have the operational experts (Flow Management Position or controllers) check that the identified flows make sense.

In particular, we are considering using the framework within the scope of the SESAR PJ08 project (Advanced Airspace Management) [START_REF][END_REF], where sector configurations could be optimised by minimizing flow cuts. In addition, we need to add further analytic capabilities in order to better understand the contribution of each flow, as well as that of the outliers, to the traffic complexity in a sector.

We would like to extend the framework by adding other clustering algorithms. A segment approach such as the one in TRACLUS seems particularly promising, where a trajectory could potentially be associated not only to a single cluster, but to a sequence of clusters, to better explain the different flight phases. Thus, in the en-route phase, the trajectories sharing the same air routes could have their en-route segments clustered, independently of whether or not they share the same departure/arrival procedures.

Figure 1 .

 1 Figure 1. Trajectory Clustering Framework process for both the ED and SSPD-based distance methods. The framework automatically identifies the flows from a set of trajectories by applying the HDBSCAN clustering algorithm and generates statistics for flow analysis.

Figure 2 .Figure 3 .

 23 Figure 2. Trajectory discretisation example

Figure 4 .

 4 Figure 4. Application of Ramer-Douglas-Peucker (RDP) algorithm example: points reduced from 100 to 15.

•

  average length/altitude/heading, • flow rate, • flight distribution per origin/destination pair, origin and destination. 4) A KML file to display in detail the clusters as well as other related ATM structures such as the airspaces and the air routes. III. RESULTS We assess the performance of both clustering methods by applying them to one day of traffic (26 June 2015) over the Reims ACC area. The dataset provided by Eurocontrol in DDR2 format contains 9442 trajectories. A bounding box with coordinates [(46.7 • , 1.328889 • ), (51.116667 • , 8.218611 • )] is defined around the Reims area to filter the trajectories geographically and different evaluations are performed on several altitude intervals. The coordinates are projected with Lambert-93 which is the official projection in Metropolitan France.

Figure 5 .

 5 Figure 5. ED-based clustering results

Figure 6 .

 6 Figure 6. ED-based clustering -Centroids for Paris flows

Figure

  Figure 8. It has to be noted that with SSPD a cluster can have trajectories with opposite directions because of the alreadymentioned limitation of this distance and the fact that the clustering is performed in 2D. Therefore, the direction displayed for the cluster centroids may be the opposite direction of some of the trajectories in the cluster.

Figure 8 .Figure 9 .

 89 Figure 8. SSPD-based clustering results for LFEEXR sector

Figure 10 .

 10 Figure 10. Cluster number 8 generated by the ED-based method

Figure 11 .Figure 12 .

 1112 Figure 11. ED versus SSPD clustering results

Figure 13 .

 13 Figure 13. Comparison between routes UM133 and UM728 (in green) and clusters 4 and 20 generated by the ED-based method (in black)

TABLE III

 III

	SSPD CLUSTERING STATISTICS EXAMPLE
	Cluster	7	8
	Avg. Heading ( • )	325	336
	Avg. Alt (ft)	36199	32753
	Avg. Length (NM)	255	75
	Flights/hour	14.7	2.5
	Main orig/dest	LFPG-EGLL	LFMN-LFPO
		5% (17)	39% (24)
	Main origin	LFPG	LFMN
		16% (56)	43% (26)
	Main destination	EGKK	LFPO
		17% (59)	89% (54)
	Main A/C type	A319	A320
		22% (76)	62% (38)

Seventh SESAR Innovation Days, 28 th -30 th November 2017

ACKNOWLEDGMENT

The authors would like to thank Sébastien Aubry, Judicaël Bedouet, Thomas Dubot, Antoine Joulia and Xavier Olive from ONERA and Pascal Lezaud from ENAC for their contributions to the present work.


Sandrine Grellier
email: sandrine.grellier@univ-orleans.fr

GENERIC COLOURFUL TORI AND INVERSE SPECTRAL TRANSFORM FOR HANKEL OPERATORS

Keywords: . 2010 Mathematics Subject Classification. 35B65, 37K15, 47B35 Cubic Szegő equation, action-angle variables, integrable system, Hankel operator, spectral analysis. P

This paper explores the regularity properties of an inverse spectral transform for Hilbert-Schmidt Hankel operators on the unit disc. This spectral transform plays the role of actionangles variables for an integrable infinite dimensional Hamiltonian system -the cubic Szegő equation. We investigate the regularity of functions on the tori supporting the dynamics of this system, in connection with some wave turbulence phenomenon, discovered in a previous work and due to relative small gaps between the actions. We revisit this phenomenon by proving that generic smooth functions and a G δ dense set of irregular functions do coexist on the same torus. On the other hand, we establish some uniform analytic regularity for tori corresponding to rapidly decreasing actions which satisfy some specific property ruling out the phenomenon of small gaps.

1. Introduction 1.1. The cubic Szegő equation. This paper explores the properties of some inverse spectral transformation related to an integrable infinite dimensional Hamiltonian system. Introduced in [START_REF] Gérard | The cubic Szegő equation[END_REF], the cubic Szegő equation reads [START_REF] Baxter | A norm inequality for a "finite-section" Wiener Hopf equation[END_REF] i∂ t u = Π(|u| 2 u) ,

where u = u(t, x) is a function defined for (t, x) ∈ R × T, T := R/2πZ, such that, for every t ∈ R, u(t, .) belongs to the Hardy space L 2 + (T) of L 2 functions v on T with only nonnegative Fourier modes, c n e inx .

∀n < 0 , v(n) = 0 .
It has been proved in [START_REF] Gérard | The cubic Szegő equation[END_REF] that ( 1) is globally well posed on Sobolev spaces H s + (T) := H s (T) ∩ L 2 + (T) for all s ≥ 1 2 , with conservation of the H 1 2 norm. Recall that, for elements in L 2 + (T), the H s Sobolev norm reads

v 2 H s = ∞ n=0 (1 + n) 2s |v(n)| 2 .
Furthermore, it turns out that (1) enjoys an unexpected Lax pair structure, discovered in [START_REF] Gérard | The cubic Szegő equation[END_REF] and studied in [START_REF] Gérard | Invariant Tori for the cubic Szegő equation[END_REF], [START_REF] Gérard | An explicit formula for the cubic Szegő equation[END_REF], [START_REF] Gérard | The cubic Szegő equation and Hankel operators[END_REF]. More precisely, consider, for every u ∈ H

1 2
+ (T), the Hankel operator H u : L 2 + (T) → L 2 + (T) defined as H u (h) = Π(uh) . Notice that H u is an antilinear realization of the Hankel matrix Γ û, where, for every sequence α = (α n ) n≥0 of complex numbers, Γ α denotes the operator on 2 (Z + ) given by the infinite matrix (α n+p ) n,p≥0 . Indeed, if F denotes the Fourier transform v → v between L 2 + (T) and 2 (Z + ), it easy to check that

F H u F -1 = Γ û • C ,
where C denotes the complex conjugation. The Lax pair identity then reads as follows, see [START_REF] Gérard | The cubic Szegő equation[END_REF]. If s > 1 2 and u is a H s + solution of (1), then

dH u dt = [B u , H u ] ,
where B u is a linear antiselfadjoint operator depending on u. As a consequence, there exists a one parameter family U (t) of unitary operators on L 2 + (T) such that ∀t ∈ R , H u(t) = U (t)H u(0) U (t) * .

In particular, H 2 u(t) = U (t)H 2 u(0) U (t) * . Notice that H 2 u is a linear positive operator on L 2 + (T), and that

F H 2 u F -1 = Γ ûΓ * û ,
thus H 2 u is a trace class operator as soon as u ∈ H 1 2

+ (T), with

Tr(

H 2 u ) = ∞ n=0 (1 + n)|v(n)| 2 = v 2 H 1 2 .
Consequently, apart from 0, the spectrum of H 2 u is made of eigenvalues, which are conservation laws of [START_REF] Baxter | A norm inequality for a "finite-section" Wiener Hopf equation[END_REF]. In fact, a second Lax pair for (1) holds [START_REF] Gérard | Invariant Tori for the cubic Szegő equation[END_REF], which concerns the operator K u := S * H u = H u S = H S * u , where S denotes the shift operator on L 2 + (T), namely multiplication by e ix . Operator K u is also a Hankel operator, F K u F -1 = Γû • C , where Γα denotes the shifted Hankel matrix (α n+p+1 ) n,p≥0 . Again, it is possible to prove that ∀t ∈ R , K u(t) = V (t)K u(0) V (t) * , for some one parameter family V (t) of unitary operators on L 2 + (T), and consequently that the eigenvalues of K 2 u are conservation laws of [START_REF] Baxter | A norm inequality for a "finite-section" Wiener Hopf equation[END_REF]. Denote by (ρ 2 j ) j≥1 the positive eigenvalues of H 2 u and by (σ 2 k ) k≥1 the positive eigenvalues of K 2 u , so that the ρ j 's -resp. the σ k 's -are the singular values of Γ û -resp. of Γû . In view of the identity

K 2 u = H 2 u -(.|u) L 2
u and of the min-max theorem, there holds the following interlacement property,

ρ 1 ≥ σ 1 ≥ ρ 2 ≥ σ 2 ≥ . . . 1.2. The spectral transform. If u belongs to a dense G δ subset H 1 2 +,gen (T) of H 1 2 + (T), one can establish -see [4] -that ρ 1 > σ 1 > ρ 2 > σ 2 > . . . We set s 2j-1 = ρ j , s 2k = σ k , j, k ≥ 1 . The s r 's are called the singular values of the pair (H u , K u ). Of course ∞ r=1 s 2 r = Tr(H 2 u ) + Tr(K 2 u ) = ∞ n=0 (1 + 2n)|û(n)| 2 < ∞ .
Conversely, given a square summable -strictly -decreasing sequence (s r ) r≥1 of positive numbers, the set of u ∈ H

1 2
+ such that the s r 's are the singular values of the pair (H u , K u ), in the above sense, is an infinite

dimensional torus T ((s r ) r≥1 ) [4] of H 1 2
+ (T). This torus is parametrised by the following explicit representation [START_REF] Gérard | The cubic Szegő equation and Hankel operators[END_REF], where we classically identify functions of L 2 + (T) with holomorphic functions u = u(z) on the unit disc such that sup

r<1 2π 0 |u(re ix )| 2 dx < ∞ .
The current element of the infinite dimensional torus T ((s r ) r≥1 ) is then given by [START_REF] Böttcher | Toeplitz Matrices, Asymptotic Linear Algebra and Functional Analysis[END_REF]

u(z) = lim N →∞ C N (z) -1 (1 N ), 1 N , |z| < 1 , 1 N :=   1 . . . 1   ∈ C N , where (3) C N (z) := s 2j-1 e iψ 2j-1 -zs 2k e iψ 2k s 2 2j-1 -s 2 2k 1≤j,k≤N
. and (ψ r ) r≥1 ∈ T ∞ is an arbitrary sequence of angles. Furthermore, the evolution of the new variables (s r , ψ r ) r≥1 through the dynamics of ( 1) is given by

ds r dt = 0 , dψ r dt = s 2 r , r = 1, 2, . . .
A natural question is then the description of the regularity of u in these new variables. A first type of answer to this question is provided by results due to Peller and Semmes -see e.g. [START_REF] Peller | Hankel Operators and their applications Springer Monographs in Mathematics[END_REF]-, which characterise the Schatten classes

r≥1 s p r < ∞ , 0 < p < ∞ , in terms of the Besov spaces ∞ j=0 2 j T |∆ j u| p dx < ∞ ,
where (∆ j u) j≥0 denotes the dyadic blocks of u. In particular, if u is smooth, then (s r ) r≥1 satisfies

∞ r=1 s p r < ∞ , ∀p < ∞ .
However, the latter condition is far from being sufficient to control high regularity of u. In fact, Sobolev regularity H s for s > 1 2 cannot be easily described by the variables (s r , ψ r ) r≥1 , as shown by the following result.

Theorem 1 ([8]

). There exists a dense G δ subset of initial data in

C ∞ + (T) := s H s + (T)
such that the corresponding solutions of (1) satisfies, for every s > 1 2 ,

∀M ≥ 1 , lim sup t→∞ u(t) H s |t| M = +∞ lim inf t→∞ u(t) H s < ∞ .
In other words, in the (s r , ψ r ) r≥1 representation, the size of the high Sobolev norms may strongly depend on the angles (ψ r ) r≥1 . The goal of this paper is to investigate this phenomenon in more detail. 1.3. Overview of the results. Our first result claims that generic smooth functions u are located on a torus T ((s r ) r≥1 ) containing also very singular functions.

Theorem 2. There exists a dense G δ subset G of C ∞ + (T) such that every element u of G belongs to H 1 2
+,gen (T), and the infinite dimensional torus T ((s r ) r≥1 ) passing through u has a dense G δ subset -for the H 1 2 topology-which is disjoint of H s for every s > 1 2 . Theorem 2 states that, on the tori T ((s r ) r≥1 ) passing through generic smooth functions, the regularity changes dramatically from C ∞ to the outside of H s for every s > 1 2 . Of course, this result can be seen as a natural extension of Theorem 1 recalled above, of which we use the weaker form that tori T ((s r ) r≥1 ) passing through generic smooth functions are unbounded in H s for every s > 1 2 . However, in order to find singular functions on these tori , we combine it with a structure property of these tori, which we think has its own interest.

Lemma 1. Let s > 1 2 and let (s r ) r≥1 be a square summable decreasing sequence of positive numbers such that the numbers {s 2 r , r ≥ 1} are linearly independent on Q. Then we have the following alternative.

• Either T ((s r ) r≥1 ) is a bounded subset of H s , • Or T ((s r ) r≥1 ) \ H s is a dense G δ subset of T ((s r ) r≥1 ) for the H 1 2
topology. The point of Theorem 2 is that, even for fast decaying singular values (s r ), the regularity of u may be spoiled by the relative smallness of the gaps s r -s r+1 with respect to s r . In fact, if

u N (z) = C N (z) -1 (1 N ), 1 N , ψ r = 0, r = 1, 2, . . . ,
with the notation introduced above, then, using the positivity property of the Hankel matrices Γ ûN and Γû N equivalent to ψ r = 0 for all rsee [START_REF] Gérard | Inverse spectral problems for Hankel operators[END_REF], [START_REF] Gérard | An inverse problem for self-adjoint positive Hankel operators[END_REF]-, we prove in the appendix that (4)

u N C 1 (T) ≥ N j=1 s 2j-1 s 2j s 2j-1 -s 2j .
It is then easy to find fast decaying sequences (s r ) such that the above right hand side tends to infinity as N goes to infinity, which implies that (u N ) is unbounded in C 1 (T). However, at this stage we do not know how to conclude that u is not in C 1 (T).

Our two other results state some uniform analytic regularity for tori T ((s r ) r≥1 ) where the sequence (s r ) r≥1 satisfies some specific property ruling out the phenomenon of small gaps. Theorem 3. For every ρ > 0, there exists δ 0 > 0 such that, for any δ ∈ (0, δ 0 ), if ∀r ≥ 1 , s r+1 ≤ δs r , all functions u ∈ T ((s r ) r≥1 ) are holomorphic and uniformly bounded in the disc |z| < 1 + ρ. Consequently, for any initial datum corresponding to some of these functions, the solution of the cubic Szegő equation (1) is analytic in the disc of radius 1 + ρ for all time, and is uniformly bounded in this disc. In particular, the trajectory is bounded in C ∞ (T).

Theorem 3 applies in particular to geometric sequences s r = e -rh for h > 0 large enough. Our last result explores in more detail the case of geometric sequences s r = e -rh , where h > 0 is arbitrary. In this case, we still obtain some uniform analytic regularity, but with a constraint on the angles ψ r . Theorem 4. Let h > 0 and θ ∈ R. Assume (s r ) is given by s r = e -rh and (ψ r ) by ψ r = rθh. Then there exists ρ > 0 such that the corresponding elements of T ((s r ) r≥1 ) are holomorphic and uniformly bounded in the disc |z| < 1 + ρ.

We do not know whether or not geometric tori are embedded into the space of analytic functions on T. What we are able to prove is that, for transcendental γ , we have the following alternative,

• Either there exists ρ > 0 such that every element of T ((γ r ) r≥1 ) is holomorphic on the disc |z| < 1 + ρ, with a uniform bound, • Or the non-analytic elements of T ((γ r ) r≥1 ) form a dense G δ subset of T ((γ r ) r≥1 ) for the H 1 2 topology. This is a special case of an extension of Lemma 1 to analytic regularity (see Lemma 4).

1.4. Open problems. In view of the above theorems, the most natural open question is certainly to decide whether Theorem 3 can be generalised to any parameter δ < 1. In particular, as we questioned above, if 0 < γ < 1, is it true that the infinite dimensional torus T ((γ r ) r≥1 ) is included into the space of analytic functions on T ?

Another question connected to Theorem 2 relies on estimate (4). Assuming that

∞ j=1 s 2j-1 s 2j s 2j-1 -s 2j = ∞ ,
can one infer that the function u ∈ T ((s r ) r≥1 ) characterised by ψ r = 0, r = 1, 2, . . . , is not C 1 on T ? In view of Lemma 1, this would imply, if moreover the s 2 r are linearly independent on Q, that most of the points on this torus would be singular -say, not in H 2 . Then it would be interesting to draw the consequences of this property for long term behaviour of solutions of the cubic Szegő equation on this torus. 1.5. Organisation of the paper. The proof of Theorem 2 is provided in Section 2 after reducing to Lemma 1 and Theorem 1. The proof of Lemma 1 combines a Baire category argument and some elementary ergodic argument for the cubic Szegő flow. Section 3 is devoted to the proof of Theorem 3, which is based on brute force estimates on matrices C N (z). In Section 4, we prove Theorem 4 by a different approach relying on the theory of Toeplitz operators and a theorem by Baxter which reduces our analysis to proving that the restriction to T of a meromorphic function given by an explicit series, has no zero and has index 0, which can be realised using some elementary complex analysis and the Poisson summation formula. Finally, the estimate ( 4) is derived in Appendix A from an explicit calculation using Cauchy matrices, in the spirit of [START_REF] Gérard | The cubic Szegő equation and Hankel operators[END_REF] and of [START_REF] Gérard | Weighted model spaces and Schmidt subpaces of Hankel operators[END_REF].

The melting pot property

In this section, we prove Theorem 2. First we reduce the proof to Lemma 1 by the following classical argument.

Lemma 2. The set of u ∈ C ∞ + (T) ∩ H 1 2 +gen (T) such that the squares s r (u) 2 , r ≥ 1 of the singular values s r (u) are linearly independent on Q, is a dense G δ subset of C ∞ + (T). Proof. From the proof of [4], Lemma 7, we already know that C ∞ + (T) ∩ H 1 2 +gen (T) is a dense G δ subset of C ∞ + (T).
In fact, we can slightly modify the proof as follows. For every N , consider the open subset O N made of functions u ∈ C ∞ + (T) such that the first singular values of H u and K u satisfy

ρ 1 (u) > σ 1 (u) > ρ 2 (u) > σ 2 (u) > . . . ρ N (u) > σ N (u) ,
and such that any non trivial linear combination of

ρ 1 (u) 2 , σ 1 (u) 2 , ρ 2 (u) 2 , σ 2 (u) 2 , . . . ρ N (u) 2 , σ N (u) 2
with integer coefficients in [-N, N ], is not zero. Approximating elements of C ∞ + (T) by rational functions, and using the inverse spectral theorem of [START_REF] Gérard | Invariant Tori for the cubic Szegő equation[END_REF] for rational functions, we easily obtain that O N is dense. The conclusion follows from Baire's theorem.

Intersecting the dense G δ subset of C ∞ + (T) provided by this lemma with the one provided by Theorem 1 -or its weaker form, saying that the corresponding Szegő trajectories are unbounded in every H s , s > 1 2 , -we observe that Theorem 2 is a consequence of Lemma 1, which we restate for the convenience of the reader. Lemma 3. Let (s r ) r≥1 be a square summable decreasing sequence of positive numbers such that the numbers {s 2 r , r ≥ 1} are linearly independent on Q and let s > 1 2 . Then we have the following alternative.

• Either T ((s r ) r≥1 ) is a bounded subset of H s , • Or T ((s r ) r≥1 ) \ H s is a dense G δ subset of T ((s r ) r≥1 ) for the H 1 2 topology.
Proof. Recall [START_REF] Gérard | Invariant Tori for the cubic Szegő equation[END_REF][START_REF] Gérard | The cubic Szegő equation and Hankel operators[END_REF] that, for the H 1 2 topology, T ((s r ) r≥1 ) is homeomorphic to the infinite dimensional tori T ∞ , endowed with the product topology, through the parametrisation given by ( 2) and (3). In particular, it is a compact metrizable space. For every s > 1 2 , the function

v H s = ∞ n=0 (1 + n) 2s |v(n)| 2 1 2
is lower semi-continuous on T ((s r ) r≥1 ). For every positive integer , consider

F = {v ∈ T ((s r ) r≥1 ) : v H s ≤ } . Then F is a closed subset of T ((s r ) r≥1
), and the complement of the union of the F is precisely T ((s r ) r≥1 ) \ H s . Hence, by the Baire theorem, either this set is a dense G δ subset of T ((s r ) r≥1 ), or there exists ≥ 1 such that F has a nonempty interior. Assume that some F has a nonempty interior, and let us show that T ((s r ) r≥1 ) is a bounded subset of H s . Let (ψ 0 r ) r≥1 ∈ T ∞ such that the corresponding point v 0 in T ((s r ) r≥1 ) lies in the interior of F . In view of the product topology on T ∞ , there exists some integer N ≥ 1 and some ε > 0 such that all the elements of T ((s r ) r≥1 ) corresponding to

ψ r ∈]ψ 0 r -ε, ψ 0 r + ε[ , r = 1, . . . , N
, form an open set U contained into F . At this stage we appeal to the number theoretic assumption on the s 2 r , which we use classically under the form that the trajectory {(ψ 0 r + ts 2 r ) r=1,...,N , t ∈ R} is dense into the torus T N . Since, as recalled in the introduction, this trajectory is precisely the projection of the trajectory of the cubic Szegő flow Φ t on the first N components, we infer that every element of T ((s r ) r≥1 ) is contained into some open set Φ t (U ). Since the cubic Szegő equation is wellposed on H s [3], we infer that T ((s r ) r≥1 ) is covered by the union of the interiors of the F m for m ≥ 1. By compactness, it is covered by a finite union, which precisely means that T ((s r ) r≥1 ) is bounded in H s .

As stated in the introduction for geometric sequences, the following analogous result holds in the analytic setting. Lemma 4. Let (s r ) r≥1 be a square summable decreasing sequence of positive numbers such that the numbers {s 2 r , r ≥ 1} are linearly independent on Q. Then we have the following alternative.

• Either there exists ρ > 0 such that every element of T ((s r ) r≥1 ) is holomorphic on the disc |z| < 1 + ρ, with a uniform bound, • Or the non-analytic elements of T ((s r ) r≥1 ) form a dense G δ subset of T ((s r ) r≥1 ) for the H

1 2 topology.
Proof. The proof is an adaptation of the preceding one (Lemma 1) to the analytic setting. As, from [START_REF] Gérard | On the radius of analyticity of solutions to the cubic Szegő equation[END_REF], the cubic Szegő equation propagates analyticity, the result follows from the Baire theorem applied to the closed sets

F := {v ∈ T ((s r ) r≥1 ) : ∞ n=0 e n |v(n)| ≤ } for ≥ 1.

Example of bounded analytic tori

In this section, we prove Theorem 3.

Proof. Let u ∈ T ((s r ) r≥1 ). Recall that u = lim N →∞ u N where u N (z) := C N (z) -1 1 N |1 N , C N (z) := s 2j-1 e iψ 2j-1 -zs 2k e iψ 2k s 2 2j-1 -s 2 2k 1≤j,k≤N and
1 N =   1 . . . 1   ∈ C N .
Our assumption is

s r+1 = ε r s r , r ≥ 1 ,
where the sequence (ε r ) r≥1 satisfies 0 < ε r ≤ δ for some δ < 1.

Our aim is to prove that, for δ sufficiently small, the functions u N are holomorphic and uniformly bounded in some disc of radius 1+ρ, where ρ > 0, independently of N . Our strategy is to use that C N (0) is related to a Cauchy matrix, and hence, that an explicit formula for its inverse is known. We write

C N (z) = C N (0) -z ĊN = C N (0)(I -zC N (0) -1 ĊN )
where

ĊN := s 2k e iψ 2k s 2 2j-1 -s 2 2k 1≤j,k≤N
and we establish the following lemma.

Lemma 5. For any 0 < δ < 1, there exists some constant C δ > 0 such that, for any N ≥ 1,

j,k (C N (0) -1 ) j,k ≤ C δ s 1 . (5)
There exists a universal constant A > 0 such that, for δ ∈ (0, 1 2 ) and for any N ≥ 1, (6)

C N (0) -1 ĊN 1 → 1 ≤ A δ .
Let us assume this lemma proved. Take ρ > 0, and choose δ 0 such that Aδ 0 (1 + ρ) ≤ 1 2 . Hence, for any δ ∈ (0, δ 0 ), from estimate ( 6),

(I -zC N (0) -1 ĊN )
is invertible for any z with |z| < 1 + ρ and its inverse R N (z) is analytic and has uniformly bounded norm for any z with |z| < 1 + ρ. Indeed, for any N ≥ 1, any z with |z| < 1 + ρ, by the Neumann series identity,

(7) R N (z) 1 → 1 ≤ ∞ k=0 |z| k C N (0) -1 ĊN k 1 → 1 ≤ ∞ k=0 2 -k ≤ 2 . Writing C N (z) -1 = (I -z(C N (0) -1 ĊN )) -1 C N (0) -1 = R N (z)(C N (0)) -1
we get

u N (z) = R N (z)C N (0) -1 (1 N ), 1 N .
Using ( 5) and ( 7), we conclude that the series defining u N converges uniformly for |z| < 1 + ρ. Hence u N is analytic and uniformly bounded in the disc of radius 1 + ρ. We infer that u is as well analytic in the disc of radius 1 + ρ and bounded on this disc. This completes the proof of Theorem 3, modulo Lemma 5.

3.1. Proof of Lemma 5. Notice that

C N (0) = diag s 2j-1 e iψ 2j-1 T , T := 1 s 2 2j-1 -s 2 2k 1≤j,k≤N .
Since T is a Cauchy matrix, its inverse is explicitly known, so the inverse of C N (0) is given by

C N (0) -1 = (-1) j+k+N α (N ) j β (N ) k s 2 2j-1 -s 2 2k 1 s 2j-1 e iψ 2j-1 1≤k,j≤N where α (N ) j := (s 2 2j-1 -s 2 2 ) <j (s 2 2 -1 -s 2 2j-1 ) >j (s 2 2j-1 -s 2 2 -1 ) β (N ) k := (s 2 2 -1 -s 2 2k ) <k (s 2 2 -s 2 2k ) >k (s 2 2k -s 2 2 )
.

In particular,

|α (N ) j | = <j s 2 2 s 2 2 -1 <j 1 -2j-2 r=2 ε 2 r 1 -2j-2 r=2 -1 ε 2 r >j (1 -2 -1 r=2j-1 ε 2 r ) (1 -2 -2 r=2j-1 ε 2 r ) s 2 2j-1 (1 -ε 2 2j-1 ) ≤ <j ε 2 2 -1 s 2 2j-1 ∞ m=1 (1 -δ 4m ) .
Indeed, in the first line above, the factors in the second product are bounded by 1, while, in the third product, the -factor is bounded by 1 1-δ 4( -j) . Similarly, we have

|β (N ) k | = <k s 2 2 -1 s 2 2 >k 1 -2 -2 r=2k ε 2 r 1 -2 -1 r=2k ε 2 r <k (1 -2k-1 r=2 -1 ε 2 r ) (1 -2k-1 r=2 ε 2 r ) s 2 2k-1 (1 -ε 2 2k-1 ) ≤ <k 1 ε 2 2 -1 s 2 2k-1 ∞ m=1 (1 -δ 4m ) Setting B δ = 1 ∞ m=1 (1 -δ 4m ) 2 , we obtain C N (0) -1 kj ≤ B δ s 2j-1 s 2 2k-1 |s 2 2j-1 -s 2 2k | <j ε 2 2 -1 <k 1 ε 2 2 -1 ≤ B δ        1 1-δ 4(k-j)+2 s 2 2k-1 s 2j-1 j≤ <k 1 ε 2 2 -1 if j < k 1 1-δ 2 s 2j-1 if j = k 1 1-δ 4(j-k-1)+2 s 2j-1 s 2 2k-1 s 2 2k k≤ <j ε 2 2 -1 if j > k.
To summarize,

C N (0) -1 kj ≤ B δ 1 -δ 2 s 2j-1    δ 2(k-j) if j < k 1 if j = k, k + 1 δ 2(j-k-1) if j > k + 1 In particular, it gives k,j (C N (0) -1 ) jk ≤ 2B δ (1 -δ 2 ) 2 j≤N s 2j-1 ≤ 2B δ s 1 (1 -δ 2 ) 3 . (8)
This proves estimate [START_REF] Gérard | Inverse spectral problems for Hankel operators[END_REF]. For the second estimate, one has to consider

C N (0) -1 ĊN 1 → 1 ≤ sup k |(C N (0) -1 ĊN ) k |. Recall that ĊN := s 2 e iψ 2 s 2 2j-1 -s 2 2 1≤j, ≤N .
In particular,

( ĊN ) j ≤ 1 1 -δ 2 s 2 s 2 2j-1 if j ≤ 1 s 2 if j ≥ + 1 As (C N (0) -1 ĊN ) k = j (C N (0) -1
) kj ( ĊN ) j , we get from the preceding estimate (8

) on |(C N (0) -1 ) kj |, • if k > (C N (0) -1 ĊN ) k ≤ B δ (1 -δ 2 ) 2 × × j≤ s 2 s 2j-1 j≤r≤k-1 ε 2 2r + +1≤j≤k s 2j-1 s 2 j≤r≤k-1 ε 2 2r + j≥k+1 s 2j-1 s 2 k+1≤r≤j-1 ε 2 2r-1 ≤ B δ (1 -δ 2 ) 2 j≤ δ 2( -j)+1 δ 2(k-j) + +1≤j≤k δ 2(j--1)+1 δ 2(k-j) + j≥k+1 δ 2(j--1)+1 δ 2(j-k-1) ≤ δ B δ (1 -δ 2 ) 2 2δ 2(k-) s≥0 δ 4s + +1≤j≤k δ 2(j--1) δ 2(k-j) As k, k≥ +1 +1≤j≤k δ 2(j--1) δ 2(k-j) = j≥ +1 k≥j δ 2(j--1) δ 2(k-j) = 1 (1 -δ 2 ) 2 one gets k, k> (C N (0) -1 ĊN ) k ≤ δB δ (1 -δ 2 ) 4 1 + 3δ 2 1 + δ 2 . • If k < (C N (0) -1 ĊN ) k ≤ B δ (1 -δ 2 ) 2 × × j≤k s 2 s 2j-1 j≤r≤k-1 ε 2 2r + k+1≤j≤ s 2 s 2j-1 k+1≤r≤j-1 ε 2 2r-1 + j≥ +1 s 2j-1 s 2 k+1≤r≤j-1 ε 2 2r-1 ≤ B δ (1 -δ 2 ) 2 j≤k δ 2( -j)+1 δ 2(k-j) + k+1≤j≤ δ 2( -j)+1 δ 2(j-k-1) + j≥ +1 δ 2(j--1)+1 δ 2(j-k-1) ≤ δB δ (1 -δ 2 ) 2 2δ 2( -k) s≥0 δ 4s + k+1≤j≤ δ 2( -j) δ 2(j-k-1)
and, as before, k, k<

(C N (0) -1 ĊN ) k ≤ δ B δ (1 -δ 2 ) 4 1 + 3δ 2 1 + δ 2 . • For k = , (C N (0) -1 ĊN ) kk ≤ B δ (1 -δ 2 ) 2 j≤ s 2 s 2j-1 j≤r≤k-1 ε 2 2r + j≥k+1 s 2j-1 s 2 k+1≤r≤j-1 ε 2 2r-1 ≤ B δ (1 -δ 2 ) 2 j≤ δ 2( -j)+1 δ 2(k-j) + j≥k+1 δ 2(j--1)+1 δ 2(j-k-1) ≤ B δ (1 -δ 2 ) 2 j≤k δ 4(k-j)+1 + j≥k+1 δ 4(j-k-1)+1 = 2 δB δ (1 -δ 2 ) 2 (1 -δ 4 )
.

Eventually, if δ ≤ 1 2 , say, we obtain, with a universal constant A,

(C N (0) -1 ĊN ) 1 → 1 ≤ sup k (C N (0) -1 ĊN ) k ≤ A δ.
This completes the proof of Lemma 5.

The totally geometric spectral data

In this section, we consider the totally geometric case and prove Theorem 4. Namely, for some fixed h > 0 and θ ∈ R, we consider the symbol u with spectral data (s r , ψ r ) with s r = e -rh and ψ r = rθh. In particular, s r+1 = s r e -h so that, for h sufficiently large, it becomes a particular case of sub-geometric spectral data treated in Theorem 3. However, the result here does not require any smallness on e -h . Our strategy here is to use Toeplitz operators and a stability result from [START_REF] Baxter | A norm inequality for a "finite-section" Wiener Hopf equation[END_REF].

4.1. Background on Toeplitz operators. Let us first introduce some basic notation. For a continuous function Φ on T, we denote by T (Φ) the Toeplitz operator of symbol Φ defined on L 2 + (T) by

T (Φ)(f ) = Π(Φf )
or equivalently the operator defined on 2 (N) by

(T (Φ)((a k ))) j := ∞ k=0 Φ(j -k)a k , j ∈ N.
For any integer N , we denote by T N (Φ) the truncated operator defined by

T N (Φ) := Π N T (Φ)Π N . Here Π N : 2 (N) → 2 (N) is the orthogonal projector : (x 0 , x 1 , x 2 , . . . ) → (x 0 , x 1 , x 2 , . . . , x N -1 ).
The operator T N corresponds to the N × N truncated Toeplitz matrix

( Φ(j -k)) 0≤j,k≤N -1 .
Recall that a sequence of N ×N matrices (A N ) N ≥1 is said to be stable if there is an N 0 such that the matrices A N are invertible for all N ≥ N 0 and sup

N ≥N 0 A -1 N 2 → 2 < ∞.
Theorem 5 ([1], [START_REF] Böttcher | Toeplitz Matrices, Asymptotic Linear Algebra and Functional Analysis[END_REF]). The sequence (T N (Φ)) N ≥1 is stable if and only if T (Φ) is invertible.

Let us emphasize that the operators are considered as operators acting on 2 (N) or L 2 + (T) so that the stability is evaluated in the 2 (N) norm. The characterization of the invertibility of Toeplitz operators is well known. We recall it for the convenience of the reader. Theorem 6. Let Φ be a continuous function on the unit circle. If Φ has index 0 and does not vanish on the circle, then T Φ is invertible on L 2 + (T). Under these hypotheses, Φ = e ϕ = Φ + Φ -with Φ + = e Π(ϕ) and Φ -= e (I-Π)(ϕ)

and the inverse of T Φ is given by T

Φ -1 + T Φ -1 - .
As an immediate consequence, one gets the following characterization of the stability of truncated Toeplitz operators.

Corollary 1. Let Φ be a continuous function on the unit circle. The sequence of truncated Toeplitz operators (T N (Φ)) is stable if and only if Φ has no zero on the unit circle and has index 0.

We are going to use this argument to prove Theorem 4.

Totally geometric spectral data and Toeplitz operators.

We claim that in the case of totally geometric spectral data, the explicit formula giving u N involves the inverse of a truncated Toeplitz operator. From direct computation, one has

C N (z) = ω 2j-1 -zω 2k |ω| 4j-2 -|ω| 4k 1≤j,k≤N = 1 ω 2j-1 1 -zω 2(k-j)+1 1 -|ω| 4(k-j)+2 1≤j,k≤N .
where ω = e -h (1-iθ) . In that case, if T N (z) and T N,r (z) denote the matrices

T N (z) = 1 -zω 2(k-j)+1 1 -|ω| 4(k-j)+2 1≤j,k≤N and
T N,r (z) = r k-j 1 -zω 2(k-j)+1 1 -|ω| 4(k-j)+2 1≤j,k≤N
we get from our explicit formula, for any r > 0,

u N (z) =< T N (z) -1 (ω 2j-1 ), 1 >=< T N,r (z) -1 (r -j ω 2j-1 ) 1≤j≤N ), (r k ) 1≤k≤N > .
We consider for |ζ| = r, |ω| 2 < r < 1, z ∈ C, the symbol

Φ(z, ζ) := ∈Z 1 -zω 2 +1 1 -|ω| 4 +2 ζ .
The transpose of the matrix

r k-j 1 -zω 2(k-j)+1 1 -|ω| 4(k-j)+2 j,k≥1
corresponds to the matrix of the Toeplitz operator of symbol

Φ(z, r•) : ζ → Φ(z, rζ).
We are going to prove the following result.

Proposition 1. There exist |ω| 2 < r < 1 and ρ > 0 such that the function ζ → Φ(z, rζ) has no zero and has index 0 on the the unit circle, for every z such that |z| < 1 + ρ.

Assuming this result proved, we obtain by Corollary 1 that, uniformly in z, |z| < 1 + ρ, T N,r (z) -1 2 → 2 is bounded (or more precisely the norm of its transpose is bounded). As |ω| 2 < r < 1, we obtain that the sequence (u

N (z)) N with u N (z) =< T N,r (z) -1 (r -j ω 2j-1 ) 1≤j≤N ), (r k ) 1≤k≤N >
is uniformly bounded and converge to u(z) for any z, |z| < 1 + ρ. We conclude as in the previous section. This ends the proof of Theorem 4.

It remains to prove Proposition 1, which is the object of the next subsections. As a preliminary, observe that, for

|ω| 2 < |ζ| < 1, γ = |ω| 2 , Φ(z, ζ) = F γ (ζ) -zωF γ (ζω 2 ), where (9) F γ (ζ) = Φ(0, ζ) = j∈Z ζ j 1 -γ 2j+1 , γ = |ω| 2 .
We collect some basic properties of function F γ in the following lemma. Lemma 6. The function F γ has a meromorphic extension in C \ {0} given by

(10) F γ (ζ) = ∈Z γ 1 -ζγ 2 .

Its only poles in

C \ {0} are the γ 2 's, ∈ Z and F γ (γ 2 +1 ) = 0, ∈ Z. Furthermore (11) F γ 1 ζ = -ζF γ (ζ), F γ ζ γ 2 = γF γ (ζ).
Proof. Let us give another expression of F γ . By assumption,

γ < |ζ| < 1 hence, |ζ| > γ 2 and F γ (ζ) = ∞ j=0 ζ j 1 -γ 2j+1 + ∞ j=0 ζ -j-1 1 -γ -2j-1 = ∞ j=0 ζ j ∞ =0 γ (2j+1) - ∞ j=0 ζ -j-1 γ 2j+1 1 -γ 2j+1 = ∞ =0 γ ∞ j=0 (ζγ 2 ) j -γζ -1 ∞ =0 γ ∞ j=0 (ζ -1 γ 2 +2 ) j = ∞ =0 γ 1 -ζγ 2 - ∞ =0 γ +1 ζ -γ 2 +2 Hence, (12)
F γ (ζ) = ∈Z γ 1 -ζγ 2 .
The other properties are elementary consequences of this equality.

Remark 1. Set γ = e -πτ , τ > 0. From the second identity (11), we observe that the meromorphic function

G τ (w) = e 2iπw F γ (e 2iπw ) 2 satisfies ∀λ ∈ Z + iτ Z , , G τ (w + λ) = G τ (w) ,
which means that G τ is an elliptic function relative to the lattice Z + iτ Z. Since G τ has only double poles at the lattice points, with singularity

1 (ζ -1) 2 ∼ - 1 4π 2 w 2
at w = 0, and since it cancels at points i τ 2 + Z + iτ Z, we infer that

G τ (w) = - 1 4π 2 P τ (w) -P τ i τ 2 ,
where

P τ (w) = 1 w 2 + λ =0 1 (w -λ) 2 - 1 λ 2
denotes the Weierstrass P function relative to the lattice Z + iτ Z. See e.g. [START_REF] Saks | Analytic functions[END_REF] 4.3. Ruling out the zeroes on the unit circle. In this section, we prove the following lemma. Lemma 7. There exists ρ > 0 such that Φ(z, ζ) does not vanish in a neighbourhood of the circle |ζ| = 1 for any z such that |z| ≤ 1 + ρ.

Lemma 7 is a consequence of the following result. Lemma 8. For every γ ∈ (0, 1),

γ 1/2 max |ζ|=γ |F γ (ζ)| < min |ζ|=1 |F γ (ζ)| .
Proof. First of all we rewrite both sides of the above inequality. If

ζ = e iθ , F γ (ζ) = ∞ k=0 γ k 1 -γ 2k e iθ + ∞ =1 γ - 1 -γ -2 e iθ = ∞ k=0 γ k 1 -γ 2k e iθ + ∞ =1 γ γ 2 -e iθ = 1 1 -e iθ + ∞ =1 γ (1 + γ 2 )(1 -e -iθ ) 1 + γ 4 -2γ 2 cos θ = (1 -e -iθ ) 1 2(1 -cos θ) + ∞ =1 γ (1 + γ 2 ) 1 + γ 4 -2γ 2 cos θ , hence |F γ (ζ)| = 1 2| sin(θ/2)| + 2| sin(θ/2)| ∞ =1 γ (1 + γ 2 ) 1 + γ 4 -2γ 2 cos θ = | sin(θ/2)| ∈Z γ (1 + γ 2 ) 1 + γ 4 -2γ 2 cos θ .
Similarly, if ζ = γe iϕ , we have

F γ (ζ) = ∞ k=0 γ k 1 -γ 2k+1 e iϕ + ∞ =0 γ --1 1 -γ -2 -1 e iϕ = ∞ k=0 γ k 1 -γ 2k+1 e iϕ + ∞ =0 γ γ 2 +1 -e iϕ = ∞ =0 γ (1 + γ 2 +1 )(1 -e -iϕ ) 1 + γ 4 +2 -2γ 2 +1 cos ϕ , so that |F γ (ζ)| = 2| sin(ϕ/2)| ∞ =0 γ (1 + γ 2 +1 ) 1 + γ 4 +2 -2γ 2 +1 cos ϕ = | sin(ϕ/2)| ∈Z γ (1 + γ 2 +1 ) 1 + γ 4 +2 -2γ 2 +1 cos ϕ .
Consequently,

min |ζ|=1 |F γ (ζ)| -γ 1/2 max |ζ|=γ |F γ (ζ)| = min θ∈T | sin(θ/2)| ∈Z γ (1 + γ 2 ) 1 + γ 4 -2γ 2 cos θ -max ϕ∈T | sin(ϕ/2)| ∈Z γ +1/2 (1 + γ 2 +1 ) 1 + γ 4 +2 -2γ 2 +1 cos ϕ Set, for x ∈ R, θ ∈ T \ {0}, f γ,θ (x) = | sin(θ/2)| γ x (1 + γ 2x ) 1 + γ 4x -2γ 2x cos θ .
Then we are reduced to proving that inf

θ∈T\{0} k∈Z f γ,θ (k) -sup ϕ∈T\{0} k∈Z f γ,ϕ k + 1 2 > 0 .
Applying the Poisson summation formula, we have

k∈Z f γ,θ (k) = n∈Z fγ,θ (2πn) , k∈Z f γ,ϕ k + 1 2 = n∈Z (-1) n fγ,ϕ (2πn), where fγ,θ (ξ) = | sin(θ/2)| R γ x (1 + γ 2x ) 1 + γ 4x -2γ 2x cos θ e -ixξ dx = | sin(θ/2)| | log γ| ∞ 0 (1 + t 2 )t -iξ/ log γ 1 + t 4 -2t 2 cos θ dt = | sin(θ/2)| 2| log γ| ∞ 0 (1 + y)y -iξ/2 log γ-1/2 1 + y 2 -2y cos θ dy ,
where we have set t = γ x , y = t 2 . We calculate the above integral by introducing the holomorphic function Finally, for θ, ϕ ∈ (0, 2π),

g(z) = | sin(θ/2)| 2| log γ| (1 + z)z -iξ/2 log γ-1/2 1 + z 2 -2z
k∈Z f γ,θ (k) - k∈Z f γ,ϕ k + 1 2 = π 2| log γ|   n∈Z cosh (π-θ)πn log γ cosh π 2 n log γ - n∈Z (-1) n cosh (π-ϕ)πn log γ cosh π 2 n log γ   = π | log γ|   ∞ n=1 cosh (π-θ)πn log γ cosh π 2 n log γ + ∞ n=1 (-1) n+1 cosh (π-ϕ)πn log γ cosh π 2 n log γ   ≥ π | log γ| ∞ n=1 1 cosh π 2 n log γ
since the second series is an alternating series of the form ∞ n=1 (-1) n+1 a n , with a n decaying to 0 as n → ∞.

Therefore min

|ζ|=1 |F γ (ζ)| -γ 1/2 max |ζ|=γ |F γ (ζ)| ≥ π | log γ| ∞ n=1 1 cosh π 2 n log γ > 0 .
The proof is complete. Lemma 8 implies that Φ has no zeroes for |ζ| = 1 and |z| ≤ 1. By continuity, it has no zeroes in a neighbourhood of this set. Hence Lemma 7 is proved.

f (C R ) (0) := 1 2iπ C R f (ζ) f (ζ) dζ.
In this section, we prove the following lemma.

Lemma 9. For any r < 1 sufficiently close to 1, the function

ζ → F γ (rζ)
has index zero on the unit circle.

Notice that Φ(0, rζ) = F γ (rζ). As the index is valued in Z and the map z → Φ(z, ζ) is smooth, Lemma 9 implies that the index of ζ → Φ(z, rζ) is zero for any z with |z| ≤ 1 + ρ as long as r is sufficiently close to 1. Proof. Proof of Lemma 9. We could use Remark 1 in order to reduce to properties of Weierstrass P function. However, for the convenience of the reader, we prefer to give a self-contained proof. Let us assume that R is chosen so that R = γ 2 , ∈ Z and F γ = 0 on C R . We consider the index of F γ on C R around 0

I(R) := Ind Fγ (C R ) (0) := 1 2iπ C R F γ (ζ) F γ (ζ) dζ.
Statement of Lemma 9 is equivalent to

I(1 -) := lim R→1 -I(R) = 0.
By definition, I is valued in Z and is continuous on the intervals corresponding to the circles avoiding the zeroes and the poles of F γ . From properties [START_REF] Gérard | Weighted model spaces and Schmidt subpaces of Hankel operators[END_REF], one has ( 13)

I(R) + I 1 R = -1, I(Rγ 2 ) = I(R).
In particular ( 14)

I(R) + I γ 2 R = -1 and (15) I(1 + ) = I((γ 2 ) + )
where I(r ± ) = lim t→r ± I(t). We are going to compute I((γ 2 ) + ) by another way, using the zeroes and the poles of F γ . Let us first collect some basic relations.

Let n be the number of zeroes in the annulus

{z ∈ C, γ < |z| < 1}.
Since there is no poles inside this annulus, one has The lemma below explicitly computes S(u). Let us give a simple proof of this formula, inspired from calculations in [START_REF] Gérard | Weighted model spaces and Schmidt subpaces of Hankel operators[END_REF]. Denote by x k , k = 1, . . . , N , the components of C (1) -1 (1). We have x k ρ + σ k .

Then Q has degree N , P has degree at most N -1 and P (ρ j ) = Q(ρ j ) , j = 1, . . . N .

Since Q -P is a unitary polynomial of degree N which cancels at ρ j , j = 1, . . . , N , we have

Q(ρ) -P (ρ) = N j=1
(ρ -ρ j ) .

jk := σ k ρ j + σ k ρ j -σ k i =j ρ i + σ k ρ j -ρ i =k ρ j + σ σ k -σ .
Multiplying and dividing µ Notice that, for every j = 1, . . . , N , R(ρ j ) = (-1) N -1 =k (ρ j + σ ) .

Since R has degree N -1, we infer R(σ) = (-1) N -1 =k (σ + σ ) , so that

N j=1 µ (N ) jk = σ k (-1) N -1 =k (σ k + σ ) =k (σ k -σ ) N i=1 ρ i + σ k ρ i -σ k ,
which is the claimed formula. The positivity of each term is an easy consequence of the inequalities ρ 1 > σ 1 > ρ 2 > σ 2 > . . .

As a consequence of Lemma 10, we retain the following inequality, obtained after discarding most of the factors bigger than 1 in each of the products.

Corollary 3. u C 1 ≥ N k=1 σ k (ρ k + σ k ) ρ k -σ k .
Notice that this implies inequality [START_REF] Gérard | Invariant Tori for the cubic Szegő equation[END_REF]. Unfortunately, at this stage we do not have arguments allowing to extend this inequality to non rational functions u, which would imply that u / ∈ C 1 if

∞ k=1 ρ k σ k ρ k -σ k = ∞ .

  ) e -inx dx 2π , n ∈ Z denotes the Fourier coefficient of v ∈ L 2 (T), and Π denotes the orthogonal projector from L 2 (T) onto L 2 + (T), Π n∈Z c n e inx = ∞ n=0

  cos θ , on the domain C \ R + , where the argument of z belongs to (0, 2π). Integrating on the contour R ε and making R → ∞, ε → 0, we obtain, by the residue theorem, assuming θ ∈ (0, 2π) with no loss of generality, fγ,θ (ξ) 1 + e πξ/ log γ = 2iπ Res g(z), z = e iθ + Res g(z), z = e -iθ = iπ sin(θ/2) | log γ| 2 cos(θ/2) 2i sin θ e θξ/2 log γ + 2 cos(θ/2) 2i sin θ e (2π-θ)ξ/2 log γ = π 2| log γ| e θξ/2 log γ + e (2π-θ)ξ/2 log γ .

4. 4 .

 4 Studying the index. Let us first recall the definition of the index. For 0 < R < ∞, we denote by C R the circle {z ∈ C, |z| = R}.Let f be a holomorphic function near C R , with no zero on C R . The index on C R around 0 of f is given by Ind

Corollary 2 .

 2 For any r < 1 sufficiently close to 1, the functionζ → Φ(z, rζ)has index zero for any z. This corollary will complete the proof of Proposition 1.

  (16) n = I(1 -) -I(γ + )From Equation (14) with R = 1 -and R = γ + , I(1 -) + I((γ 2 ) + ) = -1 and I(γ + ) + I(γ -) = -1.Substracting these equalities givesI(1 -) -I(γ + ) = I(γ -) -I((γ 2 ) + ) hence (17) n = I(γ -) -I((γ 2 ) + ).Denote by m the number of zeroes on C γ . As γ is a zero of F γ , m ≥ 1, and(18) m = I(γ + ) -I(γ -)since there is no pole on C γ . Denote by N the number of zeroes on C 1 then (19)I(1 + ) -I(1 -) = N -1since 1 is the only pole on C 1 . Now, we compute I((γ 2 ) + ). From equality (17),I((γ 2 ) + ) = I(γ -) -n = I(γ + ) -m -n from equality (18) = I(1 -) -n -m -n from equality (16) = I(1 + ) -(N -1) -m -2n from equality (19)Recalling equality (15), we conclude that N + 2n + m = 1 so n = 0 and N + m = 1. Since m ≥ 1, this implies N = 0 and m = 1. From Equation (19) and the equality I(1 + )+I(1 -) = -1 (Equation (14) with R = 1 + ), one concludes I(1 + ) = -1 and I(1 -) = 0 as required.Appendix A. A formula for the C 1 norm Let u ∈ L 2 + (T) be a rational function corresponding to the finite list of singular values ρ 1 > σ 1 > • • • > ρ N > σ N and angles ψ r = 0 for r = 1, . . . , 2N . Then we checked in[START_REF] Gérard | Invariant Tori for the cubic Szegő equation[END_REF],[START_REF] Gérard | Inverse spectral problems for Hankel operators[END_REF] that this cancellation of the angles precisely corresponds to the positivity of the operators Γ û and Γû on 2 (Z + ). The representation formula (2), (3) then reduces tou(z) = C N (z) -1 (1 N ), 1 N , with (20) C N (z) := s 2j-1 -s 2k z s 2 2j-1 -s 22k 1≤j,k≤N , Furthermore, the positivity of the Hankel matrices Γ û and Γû implies the positivity of the Fourier coefficitients of u, since, denoting by (e n ) n≥0 the canonical basis of 2 (Z + ), Γ ûe n , e n = û(2n) , Γû e n , e n = û(2n + 1) . Therefore the C 1 norm of u on T is given by S(u) := ∞ n=1 nû(n) .

σ k ρ 2 j -σ 2 k

 22 σ k ρ j -σ k =k σ k + σ σ -σ k ,where every term in the above sum is positive.Proof. We haveS(u) = u (1) = Ċ C (1) -1 (1), t C (1) -1 (1) , with C (1) := 1 ρ j + σ k 1≤j,k≤N, Ċ := 1≤j,k≤N .Notice that C (1) is a Cauchy matrix, so that the expression of C (1) -1 (1) is explicit. We have(21) C (1) -1 (1) = N j=1 (ρ j + σ k ) =k (σ k -σ ) 1≤k≤N.

  σ k ) , P (ρ) := Q(ρ) N k=1

.

  k -ρ j ) = (-1) N -1 N j=1 (σ k + ρ j ) . Since x k = P (-σ k ) Q (-σ k ) ,this yields (21). Similarly, we have(22) t C (1) -1 (1) = N =1 (ρ j + σ ) i =j (ρ j -ρ i ) 1≤j≤NComing back to the proof of Lemma 10, we have, in view of (21) and (

  jk by i =j (ρ i -σ k ), we have, for every k, k R(σ k ) =k (σ k -σ ) N i=1 ρ i + σ k ρ i -σ k , with R(σ) = N j=1 i =j ρ i -σ ρ j -ρ i =k (ρ j + σ ) .