Charlène Esmieu

Maylis Orio

Stéphane Ménage

Stéphane Torelli

Influence of the copper coordination spheres on the N 2 Or activity of a mixed-valent copper complex containing a {Cu 2 S} core

come

Influence of the copper coordination spheres on the

N2Or activity of a mixed-valent copper complex containing a Cu2S core

INTRODUCTION

Nitrous oxide (N2O) is an aggressive pollutant due to its major greenhouse effect and global warming potential (GWP, 300 times that of CO2). 1 Besides natural emissions, anthropogenic N2O sources mainly rely on agricultural activities, nitric / adipic acids factories and transportations. Independently, Health issues have been connected with long-term N2O exposures. [START_REF] Sanders | Biologic Effects of Nitrous OxideA Mechanistic and Toxicologic Review[END_REF] Among the identified effects or consequences, one can cite the inactivation on methionine synthase [START_REF] Deacon | Inactivation of Methionine Synthase by Nitrous Oxide[END_REF] , as well as genetic [START_REF] Hoerauf | Genetic damage in operating room personnel exposed to isoflurane and nitrous oxide[END_REF] or cerebrocortical damages. [START_REF] Jevtovic-Todorovic | Ketamine potentiates cerebrocortical damage induced by the common anaesthetic agent nitrous oxide in adult rats[END_REF] Consequently, mitigation of the N2O levels using efficient remediation catalysts or elaboration of selective sensors constitute important scientific breakthroughs.

Despite its chemical inertness, N2O reduction is a thermodynamically favorable process (N2O + 2H + + 2e → N2 + H2O, E° = 1.35 V @ pH 7.0 vs NHE) that Nature masters with the so-called copper-containing enzyme nitrous oxide reductase (N2Or) during the denitrification pathway. [START_REF] Pauleta | Nitrous oxide reductase[END_REF] This reactivity is remarkable since the 59 kcal/mol (gas phase) energy barrier for the N-O cleavage are overcome under physiological conditions. N2Or is a homodimer containing two different copper motifs: a dinuclear CuA site acting as electron mediator and a tetranuclear Cu4S cluster for N2O reduction. To date, two different structures are known for the Cu4S active site, depending on purification methods. The first, CuZ * , is isolated under aerobic conditions while the second, CuZ, is obtained in reductive medium (Chart 1). [START_REF] Brown | Revisiting the catalytic CuZ cluster of nitrous oxide (N2O) reductase -Evidence of a bridging inorganic sulfur[END_REF] The main difference lies in the open CuI-CuIV edge with coordinated hydroxide/water or sulfide ions. Deep experimental investigations combined with computational methods provided a major knowledge to clarify the reaction mechanism. [START_REF] Gorelsky | Mechanism of N2O Reduction by the Tetranuclear CuZ Cluster of Nitrous Oxide Reductase[END_REF] Indeed, a better understanding of how N2Or activates and reduces N2O would contribute to the perspectives described above. The involvement of a Cuz° intermediate with a µ,1-3 coordination mode for N2O at the CuI-(µ,S)-CuIV edge, together with the participation of a hydrogen bond with a neighboring lysine residue to favor the N-O bond cleavage, has been proposed. [START_REF] Johnston | Spectroscopic Definition of the CuZ° Intermediate in Turnover of Nitrous Oxide Reductase and Molecular Insight into the Catalytic Mechanism[END_REF] In our group, we aim at developing bio-inspired Cu2S-containing complexes representative of the truncated active Cu4S native cluster to propose structure/activity correlations for N2O activation/reduction. In the literature, few examples of Cux-Sy (x = 2-4 and y = 1, 2) architectures have been reported for stoichiometric N2O reduction at low temperature (-80°C). [START_REF] Bar-Nahum | Reduction of Nitrous Oxide to Dinitrogen by a Mixed Valent Tricopper-Disulfido Cluster[END_REF] We already demonstrated that mixed-valent (MV) dicopper (I,II) cores embedded in a thiophenolate-containing ligand can reduce N2O under single turnover conditions depending on the metal ions environments at room temperature and under N2O atmospheric pressure. Indeed, while [1], 11c [3] 11a and [4] 11a are inactive, [2] 11b exhibits a unique N2Or activity (Chart 1). These results point out (i) the need of exchangeable position(s) for the reactivity and (ii) the importance of one of the exogenous coordinated ligand since one position has to be occupied by a water molecule (Cu-OH2 motif).

Here, we demonstrate that the coordination sphere of the Cu-X part adjacent to Cu-OH2 within a given dinuclear complex significantly affects the kinetic at a single turnover level. Indeed, the introduction of an acetonitrile ligand for

RESULTS AND DISCUSSION

The UV-vis/NIR spectrum recorded upon dissolution of crystalline [3] in acetone lead to an intense purple solution with features at 1205 nm (ε = 760 M -1 .cm -1 ), 775 nm (ε = 855 M -1 .cm - 1 ), 570 nm (ε = 450 M -1 .cm -1 ) and 430nm (ε = 1095 M -1 .cm -1 ) that are significantly different compared to the spectrum recorded in acetonitrile used for its preparation (Figure 1). In particular, the presence of transitions in the NIR region was already observed for fully delocalized [1] and[2] (1115 nm (ε = 1190 M -1 .cm -1 ) and 1255 nm (ε = 690 M -1 .cm -1 ), respectively, Figure S1). This result indicates clear structure changes upon dissolution in that case, with the formation of a new species denoted [5] that needs to be investigated. Indeed, titration with 1,4-diazabicyclo[2.2.2] (DABCO) up to one molar equivalent resulted in drastic changes on the absorption spectrum (Figure S2). As a control experiment, no changes occurred upon addition of DABCO to an acetonitrile solution of [4] (Figure S3) that corresponds to solvated- [3] in acetonitrile. 11a Consequently, a structure for [5] in which one copper would be coordinated by an acetonitrile exogenous ligand and the other one by a water molecule is a rational hypothesis (Cu-OH2 converted into Cu-OH upon deprotonation). This new structure was successfully DFT-optimized (Figure 2, left) and showed, as for [1], [2] and [3], the presence of a Cu-Cu bond (2.7 Å) that was further confirmed by Natural Bond Order analysis (NBO, Figure 2, middle). Each copper is pentacoordinated with two nitrogen atoms from the ligand (N1/N2 for Cu1 and N3/N4 for Cu2), the bridging thiophenolate (S1), one vicinal copper ion and either one oxygen (O1) or an extra nitrogen atom (N5) from a water or an acetonitrile molecule, respectively. For both, the geometry is best described as distorted tetragonal (τCu1 = 0.28 and τCu2 = 0.27). 12 Mülliken population analysis indicates closely distributed spin densities between Cu1 (0.24), Cu2 (0.26) and S1 (0.25) that account for 75% of the total. The remaining (25%) is spread over the pyridine rings (Figure S4 for the spin density plot). The localized Singly Occupied Molecular Orbital (SOMO, Figure 2, right) displays 56% Cu and 20% S character and shows a σ-antibonding interaction between the Cu 3dz 2 orbitals and the S 3px orbital. A non-negligible degree of covalence for the Cu-S bond is expected due to the contribution of the S atom in the SOMO. For comparison, the distributed spin densities calculated for [2] are 0.27 (Cu1), 0.27 (Cu2) and 0.23 (S), and the spreading between Cu's and S (77%). This indicates that chemical modification at one metal center does not affect the spin distribution.

The calculated electronic excitations at 1200 nm, 763 nm, 693 nm and 496 nm obtained by TDDFT for optimized [5] are now in good agreement with the experiment (multi-peak analysis, 1222 nm, 776 nm, 630 nm, 560 nm, Figure S5). The nature of the transitions is finally assigned to intervalence (IVCT, CuS→CuS, NIR) and ligand to metal (LMCT) charge transfer processes (Figure S6).

Examination of [5] by EPR spectroscopy (Figure 3 and Figure S7) evidenced an interesting temperature dependence behavior. Up to 120K, the spectra are indicative of a valence-localized S = ½ system (g1 = 2.201, g2 = 2.066, g3 = 2.036, giso = 2.101; A1Cu = 545 MHz, A2Cu = 64 MHz, A3Cu = 10 MHz and Aiso = 206 MHz) with a four lines pattern (ICu = 3/2, 2nI+1 = 4 with n = 1).

Going to 292K results in drastic changes with a more complicated multi-lines signature that resembles the one previously detected for fully delocalized [1] and [2]. [5] can thus be considered as a unique example of Cu2S-containing MV species belonging to the class II category according to the Robin and Day classification. [START_REF] Day | Mixed valence: origins and developments[END_REF] Independently, the spectrum recorded upon addition of one molar eq. of DABCO at 10 K also differs from that of [5], confirming the presence of a coordinating water ligand. The parameters (g1 = 2.224, g2 = 2.074, g3 = 2.042 and giso = 2.113; A1Cu = 582 MHz, A2Cu = 98 MHz, A3Cu = 10 MHz, Aiso = 230 MHz, Figure S8) are again in that case close to the analog generated with [2] under the same conditions (g1 = 2.180, g2 = 2.050, g3 = 2.012, giso = 2.08; A1Cu = 535 MHz, A2Cu = 150 MHz, A3Cu = 10 MHz, Aiso = 232 MHz, Figure S9). These parameters correlate well with the computed ones performed on DFT optimized structures considering a Cu II -OH structure for both deprotonated [2] and [5] (Figures S10 andS11). Given the localized character, no significant difference would have been expected since in both cases the Cu II -OH motifs have identical chemical environments provided by the ligand. with exclusive N-coordinating atoms and (or not) Cu-Cu bond (Figure S12). The electronic communication within the cores that relies on the Cu2S triangle motif is thus highly sensitive to the environment. The anodic system is the most affected, with higher redox potential in presence of the metal-metal bond that tends to stabilize the MV state over a larger potential range. While the two main absorption bands at 1205 nm and 775 nm vanished within 10 minutes, a new feature is detected at 415 nm. The presence of an isosbestic point at 390 nm indicates the consumption of the starting complex and its clean conversion into a new species without side reactions. This behavior was already evidences for [2], unique Cu2S MV reported at that time for single turnover N2O reduction at room temperature. An EPR and ESI-MS follow-up of the reaction gave additional information about redox character of the reaction (Figure S13). Indeed, the axial S= ½ signal characteristic of [5] gradually disappeared concomitantly with N2O bubbling and a +2 net charge fragment consistent with the integrity of the organic framework, two Cu(II) ions and the presence of an hydroxo ligand (Figure S14) are observed. These results are thus consistent with the formation of a dinuclear Cu(II) complex with an hydroxo bridge that confers the required antiferromagnetic coupling to silent the EPR signal. Additionally, examination of the headspace gas gave clear evidence about N2 formation and confirms N2O reduction (Figure S15). A reaction rate order for [5] of 0.77 ± 0.02 under pseudo-first order conditions was determined (Figure S16), suggesting (as for [2]) the involvement of a reaction intermediate that consumes [5] along the reaction (since a two-electron process is required for the reduction). Finally, the reaction rate was calculated to be 83 ± 3 min -1 whereas 8.0 ± 0.1 min -1 was obtained for [2] (Figure S17). This result is important form a structure/activity correlation point of view. Since the aquo species are the only active ones, the coordination sphere of the adjacent Cu center has a non-negligible effect on the reaction rate and is now a new parameter to be considered. In addition, for [2], a stabilizing H-bond between the water and triflate ligands could also explain such differences in reactivity. This result is of main interest for establishing a rationale for N2O activation with these bio-inspired systems. Indeed, not only the Cu-OH2 part is essential, but also the coordination sphere of the adjacent Cu center influences the reaction. For [5] the 3NCuSCuN2O1 motif reminds that of the dissymmetric CuICuIV edge at CuZ * with exclusively nitrogen atoms in addition to the oxygen form water. This chemical environment could be the common determinant for an efficient reactivity.

CONCLUSION

Experiments

to prepare (i) new complexes combining various Cu-X spheres and the Cu-OH2 motif to confirm these result and (ii) derivatives lacking de Cu-Cu bond but still having the Cu-OH2 motif, which would to help in evaluating the influence of this structural feature on the reactivity and currently performed.

[ 5 ]

 5 (generated by dissolving [3] in acetone and fully characterized in this study) in place of a triflate ion in [2] results in a 10 times faster reduction process. A nitrogen-rich environment, as present in native N2Or, thus appears more prone to facilitate the reaction. Chart 1. Representation of the N2Or active sites and chemical structures of MV [1], [2], [3], [4] and [5]. The main differences within Cu2S cores are highlighted in red; OTf -stands for trifluoromethanesulfonate anion.

Figure 1 .

 1 Figure 1. UV-vis/NIR spectra of [3] recorded in acetone (black solid, corresponding to [5]) and acetonitrile (black dashed, corresponding to [4]) at 298 K.

Figure 2 .

 2 Figure 2. (left) Optimized structure of [5]; selected interatomic distances (Å) and angles(°): Cu1-Cu2, 2.703; Cu1-S1, 2.236; Cu1-N1, 2.083; Cu1-N2, 1.977; Cu1-N5, 2.029; Cu2-S1, 2.197; Cu2-N3, 2.059; Cu2-N4, 1.953; Cu2-O1, 2.158; Cu1-S-Cu2, 75.15 (see Table S1 in the Supporting Information for additional metric details); (middle) occupied natural orbital relevant to the Cu-Cu bond; (right) localized SOMO.

Figure 3 .

 3 Figure 3. X-band EPR spectra of [5] recorded at 10 K, 120 K and 292 K in acetone.

Figure 4 .

 4 Figure 4. Cyclic voltammogram of [5] in acetone with 0.1 M Bu4N(ClO4) as supporting electrolyte and glassy carbon as working electrode. The curve corresponds to the n+1 scan at 100 mV.s -1 .

Figure 5 .

 5 Figure 5. UV-vis monitoring of the conversion of [5] (solid black to red) upon N2O bubbling in acetone (0 to 10 min) and (insert) full UV-vis/NIR spectra of [5] (solid black) and the final species (solid red).

  In summary, a new dinuclear MV copper complex[5] containing an opened Cu2S core with acetonitrile and water molecules as exogenous ligands was generated upon dissolution of a bis(CH3CN) precursor[3] in acetone and fully characterized. The DFT-optimized structure is analogous to that of the other aqua parent [2], with the presence of a Cu-Cu bond. The electronic structure was probed by TDDFT methods and the NIR band connected to IVCT. [5] has an original (unique in this series and not observed for[2] with the (H2O)Cu-Cu(OTf) motif.) EPR temperature dependence behavior with a localized valence at least between 10 K and 120 K and a delocalized character at 292 K The electrochemical behavior of[5] reminds that of[2] and differs from the other members of the series where the water ligand and (or) the Cu-Cu bond (is) are missing ([3] and[4]). Even more interesting was the detection of a N2Or activity that led us to attest for the importance of a Cu-OH2 moiety when targeting such a reactivity, and corroborates the results obtained for[2]. Finally, the kinetic data indicate (i) a rate order under pseudo-first order conditions close to 1 in agreement with the participation of a reaction intermediate that consumes[5] to promote the two-electron reduction process and (ii) a reaction rate 10 times faster compared to[2].

ACKNOWLEDGMENTS

The authors gratefully acknowledge the Labex ARCANE and the CBH-EUR-GS (ANR-17-EURE-0003) fund for their partial financial support via the CNRS, CEA, and the University of Grenoble-Alpes. We are grateful to Colette Lebrun for her precious assistance in recording ESI-mass spectra and to Bernard Sartor for technical assistance.

ASSOCIATED CONTENT

Supporting Information

Full descriptions of materials and methods, DFT-optimized structure of [5], electronic structure of

Notes

The authors declare no competing financial interests


A Kemeny

F Colombet

T Denoual

How to Avoid Simulation Sickness in Virtual Environments during User Displacement

Driving simulation (DS) and Virtual Reality (VR) share the same technologies for visualization and 3D vision and may use the same technics for head movement tracking. They experience also similar difficulties when rendering the displacements of the observer in virtual environments, especially when these displacements are carried out using driver commands, including steering wheels, joysticks and nomad devices. High values for transport delay, the time lag between the action and the corresponding rendering cues and/or visual-vestibular conflict, due to the discrepancies perceived by the human visual and vestibular systems when driving or displacing using a control device, induces the socalled simulation sickness.

While the visual transport delay can be efficiently reduced using high frequency frame rate, the visual-vestibular conflict is inherent to VR, when not using motion platforms. In order to study the impact of displacements on simulation sickness, we have tested various driving scenarios in Renault's 5-sided ultra-high resolution CAVE. First results indicate that low speed displacements with longitudinal and lateral accelerations under a given perception thresholds are well accepted by a large number of users and relatively high values are only accepted by experienced users and induce VR induced symptoms and effects (VRISE) for novice users, with a worst case scenario corresponding to rotational displacements. These results will be used for optimization technics at Arts et Métiers ParisTech for motion sickness reduction in virtual environments for industrial, research, educational or gaming applications.

INTRODUCTION

Driving simulation and Virtual Reality (VR) are closely connected from the very beginning [START_REF] Kemeny | From Driving Simulation to Virtual Reality[END_REF]. Their history starts for both in the 1960-s and very quickly they are using both computer generated imagery (CGI), are interactive with sensory feedback and provide physical and mental immersion for an efficient use [START_REF] Shermann | Understanding Virtual Reality[END_REF].

Driving simulation (DS) and Virtual Reality (VR) technologies use the very same characteristics for visualization and may use similar technologies for head movement tracking and high end 3D vision. They also share similar difficulties in rendering movements of the observer in the virtual environments. The visual-vestibular conflict, due to the discrepancies perceived by the human visual and vestibular systems, induce the so-called simulation sickness in DS or cyber sickness in VR, when driving or displacing using a control device (ex. Joystick). Another cause for simulation sickness is the transport delay, the delay between the action and the corresponding rendering cues.

Another similarity between driving simulation and VR is the need for correct scale 1:1 perception. Correct perception of speed and acceleration in driving simulation is crucial for automotive experiments for example Advances Driver Aid Systems (ADAS) as vehicle behavior has to be simulated correctly and anywhere where the correct mental workload is an issue as real immersion and driver attention have a strong impact on it. Correct perception of distances and object size is crucial using HMDs or CAVEs, especially as their use is frequently involving digital mockup validation for design, architecture or interior and exterior lighting.

Today, the advents of high resolution 4K digital display technology allows near eye resolution stereoscopic 3D walls and integrate them in high performance CAVEs. High performance CAVEs now can be used for vehicle ergonomics, styling, interior lighting and perceived quality. The first CAVE in France, built in 2001 at Arts et Métiers ParisTech, is a 4 sided CAVE with a modifiable geometry with now traditional display technology. The latest one is Renault's 70M 3D pixel 5 sides CAVE with 4K x 4K walls and floor and with a cluster of 20 PCs. New application domains, such as automotive AR design, will bring combined features of VR and driving simulation technics, including CAVE like display system equipped driving simulators.

PERCEPTION OF MOVEMENT

Simulation sickness

If driving simulation and VR are using the same technics they share the same or similar difficulties in rendering movements of the observer in the virtual environments. The worst case is when the observer/driver is turning while driving or moving using a control device (ex. Joystick) in a CAVE or in a static driving simulator. This is due to the visuo-vestibular conflict, due to the discrepancies perceived by the human visual and vestibular systems induce the socalled simulation sickness [START_REF] Kennedy | Use of a motion history questionnaire to predict simulator sickness[END_REF] as the observer/driver may see rotations not perceived by the inner ears.

Another cause is the transport delay, the delay between the action and the corresponding rendering cues. This transport delay depends on the acquisition of observer/driver actions and movements (especially head movements), image frequency rate and motion rendering and may create strong perceptual discrepancies [START_REF] Kemeny | Recent developments in visuo-vestibular restitution of self-motion in driving simulation[END_REF]. When observing an image computed in function of head movements, small delays may induce nausea as the Vestibular Ocular Reflex (VOR), which is responsible of eye fixations, is very quick, about 20ms [START_REF] Berthoz | The Brain's Sense of Movement[END_REF]. When driving, these effects may be smaller than for passengers or passive observers as the efferent copies of information given by steering wheel, pedals or other control devices help the driver to anticipate the movement [START_REF] Dagdelen | Analysis of the visual compensation in the Renault driving Simulator[END_REF].

To avoid simulation sickness, transport delay should be limited to 50 ms in driving simulation [START_REF] Allen | Driving simulation, mechanization and application[END_REF], 20 ms using HMDs and CAVEs or head slaved or HMD based driving simulators because of the rapidity of the VOR reflex.

Other factors may play a role, as expectations based on previous experiences [START_REF] Reason | Motion Sickness[END_REF] or any other physiological parameters which influences postural stability. One way to deal with simulation sickness is to provide an acceptable face value of motion, for example by reducing sensory conflicts generating generic proprioceptive cues.

Perception of distance, speed and acceleration at scale 1:1

Correct perception of speed and acceleration in driving simulation is crucial for automotive experiments for Advances Driver Aid System (ADAS) experiments as vehicle behavior has to be simulated correctly and anywhere where the correct mental workload is an issue as real immersion and driver attention is depending on it [START_REF] Barthou | Driver trust and reliance on a navigation system: effect of graphical display[END_REF]. Speed perception is depending on the visual field of view, in part because of the role of peripheral field of view in motion perception by the human visual system and partly because of the corresponding generated immersion level and probably on other parameters, such as steering wheel force feedback and overall CGI quality [START_REF] Panerai | Speed and safety distance control in truck driving: comparison of simulation and real-world environment[END_REF]. The correct perception of accelerations requires the use of motion platform, which makes driving in a CAVE difficult, especially with rear projection floor, allowing only limited payload.

Correct perception of distances and object size is also crucial using HMDs or CAVEs (Figure 2), especially as their use is frequently involving digital mockup validation for design, architecture or interior and exterior lighting. Early VR validation studies showed already differences in distance and scale 1:1 perception between VR and real life conditions [START_REF] Loomis | Visual perception of egocentric distance in real and virtual environments[END_REF]. One of the most important visual depth cues is motion parallax and its role was shown real, virtual environments [START_REF] Panerai | Speed and safety distance control in truck driving: comparison of simulation and real-world environment[END_REF] [13] [START_REF] Kemeny | Evaluating perception in driving simulation experiments[END_REF]. Nevertheless, the precise role of binocular convergence [START_REF] Paillé | Stereoscopic Stimuli are not used in Absolute Distance Evaluation to Proximal Objects in Multi-Cue Virtual Environment[END_REF][16] and cognitive factors [START_REF] Glennerster | Humans ignore motion and stereo cues in favor of a fictional stable world[END_REF] are still to be further clarified both for scale 1:1 and distance perception.

User paradigm

Using Renault's high-performance CAVE, users require for a number of use-cases, especially for vehicle architecture and visual ergonomics studies, driving in a digital mock-up in a virtual environment, which is usually a duplicate of an existing real and well-known environment. In order to ensure an efficient use of the CAVE for these experiments and based in previous observations in the use of driving simulation experiments [START_REF] Kemeny | Evaluating perception in driving simulation experiments[END_REF], [START_REF] Kemeny | Evaluating perception in driving simulation experiments[END_REF], [START_REF] Kemeny | Evaluating perception in driving simulation experiments[END_REF], user experimentation initiated to identify the limits of possible vehicle movements, in order to keep the users in acceptable and comfortable driving conditions. A specific experimentation seemed useful in regard of the visual performances of Renault's high performance CAVE in terms of resolution, frequency, field-of-view, etc. (see section 3.1), which may have an impact on the visuo-vestibular effect. Furthermore, as in the majority of use-cases drivers are passive, simulation sickness effects may be reinforced [START_REF] Kemeny | Evaluating perception in driving simulation experiments[END_REF].

For the planned experiment, low and high level speed and acceleration thresholds were induced, thanks to the experimental driving scenarios settled using the SCANeR Studio © driving simulation software [START_REF] Kemeny | A Cooperative Driving Simulator[END_REF] to determine what amplitude of movements is acceptable when driving in a CAVE. Our hypothesis is that the higher the vehicle accelerations, the higher the visuo-vestibular conflict and the more the simulation sickness. In other terms, if the vehicle acceleration is under or close to the vestibular perception threshold, the visuo-vestibular conflict will be non-existent or weak enough to avoid Virtual Reality Induced Symptoms and Effects (VRISE).

Two types of vehicle movements were distinguished in the experiment: translations and rotations. Indeed, physical translational movements and rotational movements are perceived by the different organs in the vestibular system, the otoliths for linear and the semicircular canals for rotations [START_REF] Berthoz | The Brain's Sense of Movement[END_REF].

PROTOCOL

Experimental apparatus

Renault had recently put in use its high performance IRIS (Immersive Room and Interactive System) CAVE™, a 5 sided virtual reality room with a combined resolution of 70 M pixels, distributed over sixteen 4K projectors and two 2K projector as well as an additional 3D HD collaborative power wall. Images of the studied vehicle are displayed in real time thanks to a cluster of 20 HP Z800 computers with 24 Go RAM and 40 nVidia Quadro 6000 graphics boards. The experiment was conducted using Renault's SCANeR Studio © real time driving simulation software. The image generation frequency with SCANeR in IRIS was 60 Hz with a very detailed town environment. Head tracking was performed with an ART 8-camera infrared tracking system running at 60Hz. A driver station with controlled and measurable seat, steering wheel and driver body position was used for passive driving sessions.

General purpose

The aim of the experiment was to determine the impact of two factors on motion sickness: longitudinal accelerations and yaw rotations. Thus we built 2 scenarii to study these factors, in which the subjects are in passive driving conditions. Subjects were split into two groups, one group per scenario. And for each study, we compared two conditions: a lowlevel and a high-level stimulus. So every subject had to perform the two conditions of his or her case. In order to have no interference between the conditions, we asked subjects to perform the two conditions on two separate days. In order to avoid order effects, half of the subjects of each group performed the low-level stimulus condition first, and the other half performed the high-level stimulus condition first.

Longitudinal acceleration study

The first scenario takes place on a straight road (Figure 2). The longitudinal speed of the vehicle was controlled to generate a succession of longitudinal accelerations and decelerations. The complete motion of the controlled vehicle is obtained thanks to MADA, Renault's vehicle dynamics model.

Rotational (yaw) motion study

For the second scenario, a succession of straight and curved roads (curvature radius of 30 m, Figure 4) was used. For the two levels of stimulus, the same road was used but at different speeds. Thus the induced yaw speed and acceleration are different. Figure 5 presents these yaw speed and acceleration, given by our vehicle dynamics model (MADA). Table 2 presents the two levels of stimulus for this second scenario with 2°/s 2 and 13°/s 2 .

Table 2 -Levels of stimulus for the yaw rotation study

The vehicle yaw rotation speed can be simplified by the formula = , where V is the longitudinal speed and RC the curvature radius. As the curvature radius of the trajectory is the same in both conditions and only the longitudinal speed varies (12 km/h and 40 km/h), this is why we can observe that the rotational speeds are proportional between the 2 conditions (Figure 5, left).

Measures

Before performing the driving sessions, subjects had to read information notices and complete questionnaire about themselves and their use of virtual reality and driving simulation tools. To measure the motion sickness, the Simulator Sickness Questionnaire (SSQ [START_REF] Kennedy | Simulator Sickness Questionnaire: an enhanced method for quantifying simulator sickness[END_REF] was used. Subjects had to fill it before and after the driving to compare the results.

Postural stability was also measured, as some studies show possible relation with motion sickness [START_REF] Riccio | An ecological theory of motion sickness and postural instability[END_REF]. The Techno-Concept Stabilotest platform is a professional stabile-metric platform, used by podiatrists and allowing recording the center of gravity projection on the floor. The system can provide the area of the smallest ellipse surrounding the trajectory of the center of gravity. Two records were performed: one before and one after the driving session. The evolution of the ellipse area between the subjects was consequently measured for every subject. Each record lasted 50 seconds during which subjects were asked to stand as immobile as possible while keeping their arms along their body, keeping their eyes open and fixing a constant point in front of them.

RESULTS

Longitudinal acceleration study

Figure 6 presents the results of the longitudinal acceleration study. In both graphs, blue bars correspond to the low-level stimulus (0.44 m/s²) and red bars to the high-level stimulus (2.2 m/s²). The left graph presents the postural instability evolution (in percentage) between before and after the driving. The graph on the right presents the SSQ scores. We can notice that the majority of subjects weren't sick in any condition. On the other hand, subjects who experienced VRISE, the trend is opposite to our initial hypothesis: subjects have higher SSQ values with low speed values than with higher values Z=-1.84, p<0.06. A possible explanation is that VRISE is also induced by incoherencies induced by the transport delay, the lag between head motions and the corresponding computed images, which is more observable at low speed.

These results suggest that most subjects are accepting relatively high longitudinal accelerations, up to 0.2g, for short time periods, that is during a few seconds. On the contrary other incoherencies, due to high transport delays for example may have more severe effects on driver sickness or uneasiness, especially at low speeds.

There was no difference in postural instability between the two conditions. Furthermore, the measured postural instability values are not correlated to the SSQ scores. The ANOVA performed on the SSQ Scores showed a significant effect of rotation speed (F(2,14)=4.57, p<.05), the higher acceleration levels induced significantly higher SSQ values, showing a strong VRISE effect with the higher rotational accelerations.

Rotational (yaw) motion study

There was no difference in postural instability between the two conditions for rotational motions either. Furthermore, the measured postural instability values are not correlated to the SSQ scores. With subjects split into two groups, expert and novice subjects according to personal questionnaire answers, the ANOVA performed on the SSQ Scores showed a strong VRISE (F(2,7)=4.69, p=0.05) for novice subjects with the higher rotational motions, while experienced subjects have not experienced almost any simulation sickness (see Figure 8).

CONCLUSION

Results show that for rotational motion naïve subjects experienced significantly higher SSQ values when the acceleration values were higher (up to 13°/s 2 ) comparatively to low levels (2°/s 2 ) and with the lower level motions they have not experienced significant VRISE (nevertheless still two subjects experienced a slight uneasiness without a disturbing driving effect) and often with zero SSQ values. These naïve subjects seemed to experience an easy ride with the low, about 2°/s 2 values, though were uneasy with the higher, 13°/s 2 values.

Results also show that the group of experienced subjects had low SSQ values in the two different sessions, independently on the acceleration values and had always zero or very low SSQ values. Thus it would seem that with experienced subjects, such as gamers or VR specialists, these rotational VRISE are low but for naïve subjects low rotational motion levels are required to avoid significant VRISE.

In addition we can also observe that SSQ and postural stability values are not correlated in our experiments. Previous studies have already shown that the relationship between motion sickness and postural stability may vary according to motion characteristics [START_REF] Bos | Nuancing the relationship between motion sickness and postural stability[END_REF], especially in driving situations where even negative correlation may be observed [START_REF] Reed-Jones | The relationship between postural stability and virtual environment adaptation[END_REF], well in accordance with our findings.

Figure 1 -

 1 Figure 1 -Renault's High Performance CAVE with controlled driver station

Figure 3

 3 Figure 3 presents the obtained vehicle acceleration profile. This profile was repeated 8 times by subject for a total duration of 128 seconds.

Figure 2 -Figure 3 -

 23 Figure 2 -Screenshot of the straight road used for the first scenario (longitudinal acceleration study)

Figure 4 -Figure 5 -

 45 Figure 4 -Presentation of the curved road used for the second scenario (yaw motion study)

Figure 6 -

 6 Figure 6 -SSQ scores and postural instability evolution for the two longitudinal acceleration levels

Figure 7

 7 Figure 7 presents the results for the rotational (yaw) motion study. In both graphs, blue bars correspond to the low-level stimulus and red bars to the high-level stimulus. The left graph presents the postural instability evolution (in percentage) between before and after the driving. The graph on the right presents the SSQ scores.

Figure 7 -

 7 Figure 7 -Postural instability evolution and SSQ scores for the two rotational (yaw) motion levels

Figure 8 -

 8 Figure 8 -SSQ scores for the two rotational (yaw) motion levels -expert subjects on the left and novice subjects on the right

Table 1

 1

										15
										10
			_r
	-20				---Travel at 12 km /h			-15					---Travel at 12 km /h
					-Travel at 40 km /h								-Travel at 40 km /h
	-250	100	200	300	400	500	600	700	800	-20	100	200	300	400	500	600	700	800
				Travelled Distance [m]							Travelled Distance [m]

presents the two proposed levels of stimulus for this first scenario, generating 0.04g and 0.2g acceleration levels:

Table 1 -Levels of stimulus for the longitudinal acceleration study

 1


Louis Guillot
email: guillot@lma.cnrs-mrs.fr

Bruno Cochelin
email: bruno.cochelin@centrale-marseille.fr

Christophe Vergez
email: vergez@lma.cnrs-mrs.fr

A Taylor series based continuation method for solutions of dynamical systems

Keywords: Nonlinear dynamics, Numerical continuation, Quadratic recast, Asymptotic Numerical Method, Taylor series, Dynamical systems

This paper describes a generic Taylor series based continuation method, the so-called Asymptotic Numerical Method, to compute the bifurcation diagrams of nonlinear systems. The key point of this approach is the quadratic recast of the equations as it allows to treat in the same way a wide range of dynamical systems and their solutions. Implicit Differential-Algebraic Equations, forced or autonomous, possibly with time-delay or fractional order derivatives are handled in the same framework. The static, periodic and quasi-periodic solutions can be continued as well as transient solutions.

Introduction

Numerical continuation methods are commonly used in different scientific communities. In this paper, one Taylor series based continuation method, also called the Asymptotic Numerical Method (ANM), is applied to various solutions of different types of dynamical systems. When the dynamical system of interest is an Ordinary Differential Equations (ODE) system, there are a lot of methods available in the literature. The continuation of equilibrium is a common possibility of [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF][START_REF] Doedel | Numerical computation of periodic solution branches and oscillatory dynamics of the stirred tank reactor with a yields b yields c reactions[END_REF][START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF]. The detection and the continuation of bifurcation points has also been investigated in [START_REF] Dhooge | Matcont: a matlab package for numerical bifurcation analysis of odes[END_REF][START_REF] Doedel | Auto: A program for the automatic bifurcation analysis of autonomous systems[END_REF][START_REF] Khibnik | Continuation techniques and interactive software for bifurcation analysis of odes and iterated maps[END_REF]. These tools often allow the continuation of periodic orbits of ODE and of their bifurcations. More recently, some numerical tools were generalized to the continuation of equilibrium and periodic orbits of Delay Differential Equations [START_REF] Engelborghs | Numerical bifurcation analysis of delay differential equations using dde-biftool[END_REF] and of neutral type DDE [START_REF] Barton | Collocation schemes for periodic solutions of neutral delay differential equations[END_REF]. The methods developed for the representation of periodic solutions can be time-domain methods e.g piecewise polynomial orthogonal collocation [START_REF] Doedel | Computation of periodic solution bifurcations in odes using bordered systems[END_REF] or shooting [START_REF] Peeters | Nonlinear normal modes, part ii: Toward a practical computation using numerical continuation techniques[END_REF][START_REF] Seydel | From equilibrium to chaos: practical bifurcation and stability analysis[END_REF], mixed time and frequency methods e.g trigonometric collocation [START_REF] Jean | Periodic response investigation of large order non-linear rotordynamic systems using collocation[END_REF], alternating frequency-time methods [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF][START_REF] Farner | Contribution to harmonic balance calculations of self-sustained periodic oscillations with focus on single-reed instruments[END_REF], and frequency domain method with the full harmonic balance method [START_REF] Krylov | Introduction to non-linear mechanics[END_REF][START_REF] Nakhla | A piecewise harmonic balance technique for determination of periodic response of nonlinear systems[END_REF]. In the last few years, progress has been achieved on the continuation of quasi-periodic solutions of ODE, using different methods. In [START_REF] Schilder | Continuation of quasi-periodic invariant tori[END_REF] the invariant tori are continued when in [START_REF] Fontanela | Computation of quasi-periodic localised vibrations in nonlinear cyclic and symmetric structures using harmonic balance methods[END_REF][START_REF] Zhou | A variable-coefficient harmonic balance method for the prediction of quasi-periodic response in nonlinear systems[END_REF], a two dimensional alternating frequency-time scheme is developed. In [START_REF] Guillot | Continuation of quasi-periodic solutions with two-frequency harmonic balance method[END_REF], the system of equations is recast quadratically and a full quasi-periodic harmonic balance method is applied. Most of the methods presented above are designed for steady states regimes. However, they generally include the possibility of continuing solutions of boundary value problems and initial value problems if they are formulated as solutions of algebraic systems.

Commonly, the continuation tools use predictor-corrector methods for which the solutions on the branch are computed from (a) previous solution(s) and a predictor method, generally of first order. Then, a corrector is applied to the predicted solution to ensure that the corrected solution belongs to the solution branch. On the contrary, in this paper the continuation method is the Asymptotic Numerical Method (ANM) which is based on a high-order Taylor series expansion of the solution branch. It can be seen as a high-order predictor method which most of the time does not require a correction step. In a way, this is similar to homotopy techniques which are presented in [START_REF] Liao | Beyond perturbation: introduction to the homotopy analysis method[END_REF]. The representation as a Taylor series gives a continuous representation of the solution branch, as opposed to predictor-corrector methods. Three main approaches to compute the Taylor series coefficients of the solution branch can be considered. The first one is to do the calculation by hand for a specific set of equations [START_REF] Boumédiène | Nonlinear forced vibration of damped plates by an asymptotic numerical method[END_REF][START_REF] Claude | Eigensolutions to a vibroacoustic interior coupled problem with a perturbation method[END_REF]. The second one uses automatic differentiation tools as in [START_REF] Charpentier | Towards a full higher order ad-based continuation and bifurcation framework[END_REF][START_REF] Charpentier | The diamant approach for an efficient automatic differentiation of the asymptotic numerical method[END_REF]. The third one relies on the quadratic recast of the equations. This latest approach is the one developed by the authors for its genericity and its performance [START_REF] Guillot | A generic and efficient taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF]. It has first been described in [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF][START_REF] Cochelin | Asymptotic-numerical methods and pade approximants for non-linear elastic structures[END_REF][START_REF] Cochelin | The asymptotic-numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF] in the case of equilibrium branches of structures and then extended to periodic solutions of ODE using the Harmonic Balance Method [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. The generalization to the computation of quasi-periodic orbits has been done in [START_REF] Guillot | Continuation of quasi-periodic solutions with two-frequency harmonic balance method[END_REF]. The types of dynamical systems that can be treated in this framework is limited only by the need of a quadratic recast of the equations. This is far less restrictive than it could seem at first sight as will be shown throughout this paper. Periodic solutions of fractional order systems can be continued in this framework [START_REF] Vigué | Continuation of periodic solutions for systems with fractional derivatives[END_REF]. Some very recent work [START_REF] Guillot | A generic and efficient taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF] has shown that the periodic solutions of delay differential equations (possibly of neutral type) could also be investigated with this method.

This article generalizes the approach given in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. The first section briefly recalls the content of [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] and gives the generalization to a large class of systems and to a large class of their solutions. The remaining sections treat some examples of applications in detail. For clear brevity reasons, the examples treated do not cover all the combinations of solution type and system type. In the second section, the solutions sought for are all timeperiodic. The periodic solutions of five types of dynamical systems, namely Ordinary Differential Equations (ODE), Differential Algebraic Equations (DAE), Implicit Differential Equations (IDE), Delay Differential Equations (DDE) and fractional order equations are investigated. It shows the versatility of our method in terms of dynamical system type. In the third section, a system of ODE is fixed. Its equilibrium, its periodic solutions, its quasi-periodic solutions and some of its transient regimes are continued. It shows the versatility of our method in terms of solution type.

Elements of theory

The quadratic framework and the Asymptotic Numerical Method

This section consists on a brief recall of the quadratic framework that is used by the continuation method developed by the authors. This is a key idea that allows a low time computation, an automatic and exact Jacobian matrix of the nonlinear system considered and a generic implementation which minimizes problem dependent implementation.

The so called Asymptotic Numerical Method (ANM) first described in [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF][START_REF] Cochelin | Asymptotic-numerical methods and pade approximants for non-linear elastic structures[END_REF][START_REF] Cochelin | The asymptotic-numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF]] is based on the numerical continuation of algebraic systems of the form

R(u, λ) = 0 (1)
where u ∈ R n , λ ∈ R is the continuation parameter and R : R n × R → R n is an analytic function of its arguments. From now on, the vector of all the unknowns U = (u, λ) is introduced to simplify the notations. Let U 0 be a regular solution1 of the system (1). Let U 1 be a tangent vector at U 0 . The classical arc-length parameter [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF] a = U t 1 (U 0 -U) is introduced and the solution branch around U = U 0 is written as a power series2 with respect to a

U(a) = U 0 + aU 1 + a 2 U 2 + a 3 U 3 + • • • + a N U N .
(

) 2
This writing of the solution branch as a high-order power series is the base of the ANM. Several methods have been developed to compute the coefficients U p , p ≥ 1 of the series. The approach presented in this paper consists on a quadratic recast of the equations (1) with the use of auxiliary variables U a :

R f (U, U a ) = R f (U f ) = C + L(U f ) + Q(U f , U f ) = 0 (3)
where

U f ∈ R n f +1
, n f ≥ n and C, L and Q are respectively a constant, a linear and a quadratic operator with values in R n f . The subscript f stands for full and the subscript a stands for auxiliary. Once this recast is given, the auxiliary variables U a are expanded in Taylor series (2) in the same manner as for U. The expansions are then introduced in the quadratic equations (3) and the terms of same order in a are collected. It results in N linear systems of size n f to be solved (sharing the same matrix to be inverted). The systems can be solved with the help of an efficient condensation algorithm which finally result in the Taylor series (2) sought for. The interested reader is referred to [START_REF] Guillot | A generic and efficient taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF] for all the details on the practical computation of the series. It shall be emphasized here that once the quadratic recast is available, all the operations necessary to compute the Taylor series developments (2) can be automatized. Several examples of quadratic recasts are given throughout the paper. The next two subsections 1.1.1 and 1.1.2 give a hint of how generic this framework is through two examples. Firstly, it is shown how to recast polynomials and rational fractions quadratically. Secondly, a general method for the quadratic recast of transcendental functions is given. In the appendix C, the recurrence relations to compute the Taylor series development for these two examples are detailed.

Rational nonlinearity

The first step is to be able to recast rational fractions quadratically. To show the general idea of the method, the variable y = 1

x + x 3 is recast quadratically. For this purpose, two auxiliary variables u = 1

x and v = x 2 are introduced. The second is already quadratic while the first one can be written implicitly ux = 1. The case x = 0 is not an issue here as this value is not a solution of the equation ux = 1. The final quadratic recast of the variable y is finally

   y -u -xv = 0 ux -1 = 0 v -x 2 = 0 (4)
This procedure can easily be generalized to any rational fraction. The appendix C details the computation of the coefficients of the series (2) for this example.

Transcendental nonlinearity

The next fundamental step is to be able to recast functions like exp, cos, sin, arctan, . . . quadratically. The key is to use a differentiated form of these variables. In the book [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF], these recasts are discussed in more details and a table of recasts of common functions is given in the article [START_REF] Guillot | A generic and efficient taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF]. Hereafter, the example of sine function's recast is given as an highlight of the general approach.

Let us define the variable y = sin(x). Its differentiated form is dy = cos(x) dx. An auxiliary variable z = cos(x), called companion variable, is introduced with its differentiated form dz = -sin(x) dx. The final recast is y -sin(x) = 0 and dyzdx = 0 z -cos(x) = 0 and dz + ydx = 0

(5)

The original (not differentiated form) is used to compute the zero-th order coefficient of the series (2) and the differentiated form is used to compute higher order coefficients. This works easily as the differentiated form is quadratic with respect to the variables (x, y, z, dx, dy, dz). This form is very generic and can be computed for all the elementary transcendental functions [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF][START_REF] Guillot | A generic and efficient taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF]. It gives the possibility to recast a very wide range of nonlinearities quadratically. The appendix C details the computation of the coefficients of the series (2) for this example.

The classical ODE framework for periodic solution by HBM and ANM

In [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], the case of periodic solutions of Ordinary Differential Equations (ODE) was treated. Let us consider the following ODE system :

ẏ = f (t, y, λ) (6)
with y an unknown function of the time and λ the continuation parameter. f is an analytic function of its arguments and is time-periodic. As in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], the unknown function y is expanded in truncated Fourier series :

y(t) = Y 0 + H k=1 Y c,k cos(kωt) + H k=1 Y s,k sin(kωt) (7)
The Harmonic Balance Method (HBM) consists in introducing this ansatz of y in the equations [START_REF] Charpentier | Towards a full higher order ad-based continuation and bifurcation framework[END_REF] and to balance the 2H + 1 harmonic terms corresponding to y's Fourier development. This yields an algebraic system of 2H + 1 vectorial equations for 2H + 1 vectorial unknowns Y 0 , Y c,k , Y s,k . The angular frequency ω is an additional unknown which is determined either by the angular frequency of the forcing term if the system is forced or by a phase condition [START_REF] Seydel | From equilibrium to chaos: practical bifurcation and stability analysis[END_REF] if the system is autonomous.

A key point is that the algebraic system obtained after balancing the harmonics is automatically quadratic if the original system (6) has been recast quadratically before applying the HBM. Therefore, in the same manner as in the previous section and as explained in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], auxiliary variables y a are introduced to obtain from (6) a system of the form

d( ẏf ) = c 0 (t) + λc 1 (t) + l 0 (y f ) + λl 1 (y f ) + q(y f , y f ) (8)
where y f = (y, y a ) is the full vector of unknowns and c 0 ,c 1 ,l 0 ,l 1 ,d,q are respectively two constant operators, three linear operators and one bilinear operator. The form of equation [START_REF] Claude | Eigensolutions to a vibroacoustic interior coupled problem with a perturbation method[END_REF] has been chosen to take into account that the continuation parameter λ is not intended to be developed in Fourier series. In the reference [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] it is shown how to recover the constant, linear and bilinear operators of the algebraic system obtained after the harmonic balance from the operators of the system [START_REF] Claude | Eigensolutions to a vibroacoustic interior coupled problem with a perturbation method[END_REF]. It essentially gives the following algebraic system :

ωD(Y f ) = C 0 + λC 1 + L 0 (Y f ) + λL 1 (Y f ) + Q(Y f , Y f ) (9)
where the capital letters stand for frequency domain operators and vector of Fourier coefficients. It means that, once the equations ( 6) has been recast in the form ( 8), its periodic solutions can be continued with the ANM as explained in section 1.1.

Generalization of the quadratic recast of ODE to a large class of dynamical systems and the continuation of their solution using ANM

Let us consider a differential system of the general form :

f (t, y, ẏ, ÿ, y(t -τ ), D α (y), λ) = 0 ( 10
)
where y is an unknown function of the time, D α is the Weyl fractional order α derivative [START_REF] Vigué | Continuation of periodic solutions for systems with fractional derivatives[END_REF] operator and f is an analytic function of its arguments which is time-periodic. The main objective of this paper is to show that a quadratic recast of this equation is possible and to show how this allows many types of solutions to be continued (equilibrium, periodic or quasi-periodic, transients) within a unique framework. The first step is to introduce two auxiliary vectors y τ a (t) = y(tτ ) and y α a = D α (y) in order to isolate the time delays and the fractional order derivatives. These two terms are rewritten quadratically following the frequency-domain implementation developed in [START_REF] Vigué | Continuation of periodic solutions for systems with fractional derivatives[END_REF] and [START_REF] Guillot | Continuation of periodic solutions of various types of Delay Differential Equations using Asymptotic Numerical Method and Harmonic Balance Method[END_REF] respectively. The system [START_REF] Cochelin | Asymptotic-numerical methods and pade approximants for non-linear elastic structures[END_REF] can now be written

f (t, y, ẏ, ÿ, y τ a , y α a , λ) = 0 y τ a (t) -y(t -τ ) = 0 y α a -D α (y) = 0 (11)
The two last equations are treated separately using appropriate methods (see [START_REF] Vigué | Continuation of periodic solutions for systems with fractional derivatives[END_REF] and [START_REF] Guillot | Continuation of periodic solutions of various types of Delay Differential Equations using Asymptotic Numerical Method and Harmonic Balance Method[END_REF]). The first equations are a set of Implicit Differential Algebraic Equations (IDAE). They are treated in the same manner as the ODE system of the previous section, following the idea of [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. Hence, the second step is to recast the system in a quadratic format with the help of auxiliary variables y a :

c 0 (t) + λc 1 (t) + λ 2 c 2 (t) + l 0 (y f ) + λl 1 (y f ) + q(y f , y f ) + d( ẏf ) + λd 1 ( ẏf ) + dd( ÿf ) = 0 ( 12
)
The form of these equations generalizes the form [START_REF] Claude | Eigensolutions to a vibroacoustic interior coupled problem with a perturbation method[END_REF]. c 0 ,c 1 ,c 2 ,l 0 ,l 1 ,d,d 1 ,dd,q are respectively three constant operators, five linear operators and one bilinear operator. Again, this form as been chosen to take into account that λ is constant in time.

Once this recast is available, the time-continuous system obtained ( 12) can be discretized automatically with different methods. The next two subsections 1.3.1 and 1.3.2 show how to obtain a quadratic algebraic system of the form (3) that can be continued with the ANM. Firstly, time domain discretization methods are discussed and secondly frequency domain methods.

Time-domain discretization methods

Without any requirement on the type of solution sought for, an integration scheme can be applied to the quadratic system [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. Given a set of consistent initial conditions3 y f (0), ẏf (0), the solution can be approximated on the grid {0, h, 2h, . . . , M h} (h > 0 is the integration step and M a given integer) for instance by a finite difference scheme with centered derivatives. For 0 ≤ k ≤ M -1, the equations are then

0 = c 0 (kh) + λc 1 (kh) + λ 2 c 2 (kh) + l 0 (y f (kh)) + λl 1 (y f (kh)) + q(y f (kh), y f (kh))+ d( y f ((k+1)h)-y f ((k-1)h 2h ) + λd 1 ( y f ((k+1)h)-y f ((k-1)h 2h ) + dd( y f ((k+1)h)-2y f (kh)+y f ((k-1)h h 2 ) (13)
Like the system (12), this system is quadratic in the unknowns y f (mh), 0 ≤ m ≤ M and can be continued with the Asymptotic Numerical Method. In this case, the transcendental functions can be treated exactly in the same manner as briefly explained in section 3.4.

It is clear here that many other schemes could be used. The last example of this paper 3.4 describes an implementation of a Newmark scheme. In appendix B, time-integrations of a pendulum are performed with a class of Euler schemes that cover the explicit Euler, the implicit midpoint and the implicit Euler schemes. A collocation method could also be used to solve the system of equations, as done in [START_REF] Karkar | A comparative study of the harmonic balance method and the orthogonal collocation method on stiff nonlinear systems[END_REF] for the continuation of periodic solutions.

Frequency-domain discretization methods

Following the idea of [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] recalled in section 1.2, the method used to compute periodic solutions is here a purely frequency-domain method. The solutions of the quadratic system [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] are expanded in Fourier series [START_REF] Charpentier | The diamant approach for an efficient automatic differentiation of the asymptotic numerical method[END_REF] and after balancing the harmonics, in the same manner as for the equations (9), the system obtained is :

0 = C 0 + λC 1 + λ 2 C 2 + L 0 (Y f ) + λL 1 (Y f ) + Q(Y f , Y f )+ ωD(Y f ) + λωD 1 (Y f ) + ω 2 DD(Y f ) (14)
The two last terms of this system are cubic in the unknowns (Y f , ω, λ) of the system. Hence, two additional auxiliary variables are introduced Λ = λω and Ω = ω 2 to obtain the final quadratic recast

0 = C 0 + λC 1 + λ 2 C 2 + L 0 (Y f ) + λL 1 (Y f ) + Q(Y f , Y f )+ ωD(Y f ) + ΛD 1 (Y f ) + ΩDD(Y f ) 0 = Ω -ω 2 0 = Λ -λω (15)
As opposed to the time-domain discretization method, the transcendental nonlinearities cannot be treated directly with the formalism recalled in section 1.1.2 that is detailed in the book [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF] and the article [START_REF] Guillot | A generic and efficient taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF]. The reference [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF] states clearly that it is sufficient to treat the case of an auxiliary variable v(t) = g(u(t)) where g is a transcendental function, which is defined with a differential system of the form :

v(0) = g(u(0)) v(t) = w(t) u(t) (16)
where the variables u,v and w are periodic functions and ġ(u(t)) = w(t). The final quadratic recast with U, V, W being the vectors of Fourier coefficients of u, v, w respectively is

v(0) = g(u(0)) and dv(0) = w(0)du(0) dL(V) = dQ(W, U) (17)
where dL and dQ are suitable linear and bilinear operators respectively. The details on this algebra are given in appendix A.

The generalization of this approach to quasi-periodic solutions is explained in [START_REF] Guillot | Continuation of quasi-periodic solutions with two-frequency harmonic balance method[END_REF]. The solutions are expressed as double truncated Fourier series of the form

y(t) = Y 0 + H1 k1=1 H2 k2=1 (Y c,k1,k2 cos(k 1 ω 1 t + k 2 ω 2 t) + Y s,k1,k2 sin(k 1 ω 1 t + k 2 ω 2 t)) . (18)
As there are two independent angular frequencies, an additional phase condition must be added to solve the algebraic system obtained after balancing the harmonics. For brevity the algebra is not detailed in this paper. The interested reader is referred to the article [START_REF] Guillot | Continuation of quasi-periodic solutions with two-frequency harmonic balance method[END_REF] or to the PhD thesis [START_REF] Vigué | Solutions périodiques et quasi-périodiques de systèmes dynamiques d'ordre entier ou fractionnaire : applications à la corde frottée[END_REF] for further details.

Periodic solutions of ODE, DAE, IDE, DDE and fractional order systems

In this section, periodic solutions of ODE, DAE, IDE, DDE and fractional order systems are continued.

Conservative and non-conservative systems are both treated. The functions appearing below are supposed analytic on their domain of definition. As the solutions investigated in this section are all periodic, all the functions of the time appearing are supposed periodic. The example of the pendulum is treated in three different formulations : solving the polar coordinates equations of motion, solving the Cartesian coordinates equations of motion and solving the equation of conservation of the energy. This illustrates the ability to cope with ODE, DAE and IDE in the same framework and to check that the results are consistent. A Van der Pol oscillator is then treated with different feedbacks, to show the application of the methods to fractional order systems and to delay differential equations.

Ordinary Differential Equations : Usual polar coordinates pendulum

The ODE are a very widespread type of equations that can be given in the following form:

ẏ = f (y, t), where y and f (y, t) ∈ R n (19)
The equations of the motions of the pendulum in polar coordinates are given by:

θ = φ φ = -sin(θ) -λφ (20)
where λφ is an unfolding term [START_REF] Munoz-Almaraz | Continuation of periodic orbits in conservative and hamiltonian systems[END_REF] added to compute the family of periodic solutions easily. Once the system solved, λ is found to be zero since no periodic solution exists otherwise. These equations are recast quadratically, in the form of equations ( 12) and ( 16) : These equations are treated automatically with the method proposed in section 1.3.2. The results are shown and discussed in section 2.4.

       θ = φ φ = -

Differential Algebraic Equations : Cartesian coordinates pendulum

The DAE are a type of equations where there are both differential equations and algebraic equations between the unknowns. It can be given in the following form:

ẏ1 = f (y 1 , y 2 , t), where y 1 , f (y 1 , y 2 , t) ∈ R n diff and y 2 , g(y 1 , y 2 , t) ∈ R n alg 0 = g(y 1 , y 2 , t) (22)
The equations of the motions of the pendulum in Cartesian coordinates are given by:

   ẍ + λ ẋ + T x = 0 ÿ + λ ẏ + T y -1 = 0 x 2 + y 2 -1 = 0 ( 23
)
where T is the tension of the string, x, y are respectively the horizontal and the vertical position of the pendulum. Similarly to the previous case, λ ẋ and λ ẏ are unfolding terms [START_REF] Munoz-Almaraz | Continuation of periodic orbits in conservative and hamiltonian systems[END_REF] added to compute the family of periodic solutions easily. Once the system solved, λ is found to be zero since no periodic solution exists otherwise. These equations are already quadratic according to equations [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] and are then treated automatically with the method proposed in section 1.3.2. The results are shown and discussed in section 2.4.

Implicit Differential Equations : Conservation of the Energy

The IDE are a type of equations where the derivatives cannot be expressed explicitly (as for the ODE ( 19)).

It can be written

f ( ẏ, y, t) = 0, where y, f ( ẏ, y, t) ∈ R n (24)
The potential energy of the pendulum is defined as -1 -cos(θ) so that its maximum value for periodic orbits is 0. The equation of conservation of the total energy of the pendulum is given by: θ2

2 -1 -cos(θ) = λ ( 25
)
where λ is the energy of the pendulum. Its quadratic recast is given by:

         φ 2 2 -1 -c = λ φ = θ s(0) = sin(θ(0)) and ṡ = c θ c(0) = cos(θ(0)) and ċ = -s θ ( 26
)
These equations are compatible with equations ( 12) and ( 16) and are then treated automatically with the method proposed in section 1.3.2. The results are shown and discussed in section 2.4.

Bifurcation diagram of the pendulum

The branch of periodic solutions of the pendulum is represented on the bifurcation diagrams in figure 1. The relative error ǫ between the theoretical value T th and the computed value T of the period of the oscillations is given by ǫ = T th -T T th

. The theoretical value is given by T th = 4K(sin( max(θ) 2

)

2 ) where K is the complete elliptic integral of the first kind (see [START_REF] Lévy | Précis élémentaire de la théorie des fonctions elliptiques avec tables numériques et applications[END_REF] or [START_REF] Marquis De Laplace | [END_REF] for what could be the first derivations of this formula). The relative error gets bigger while the energy gets closer to its maximum value 0. It is probably because the computation would need more than just a double precision setting and the solution more than H = 500 harmonics to be well represented in this area. The periodic solution represented in figure 2 shows a solution with an energy equal to -10 -14 , very close to its maximum value 0. We recall here that the reference of the potential energy of the system has been taken so that the maximal value of the total energy of a periodic solution is 0. Fig. 2 A solution of the pendulum with the energy equals to -10 -14 (that is very close to its maximum value 0) and H = 500 harmonics. It is almost a homoclinic orbit in the sense that it remains a very long time very close to the unstable equilibrium θ = π (or -π).

The three bifurcation diagrams for the three different formulations are superimposed on the top figure 1. It is not possible to distinguish the three bifurcation diagrams. On the bottom figure 1 the relative errors between the analytic reference for the period and the one computed by our method are shown. This error remains very small (under one percent) but shows that the three formulations do not give exactly the same results. This is probably due on the one hand to the numerical accuracy and on the other hand to the system solved that are not exactly the same.

The time computation for ten continuation steps is around 45 seconds for each system. The bifurcation diagrams obtained with ODE and DAE formulations require 50 continuation steps and the one with IDE formulation 75 continuation steps. The IDE formulation [START_REF] Jean | Periodic response investigation of large order non-linear rotordynamic systems using collocation[END_REF] has directly the energy as its continuation parameter. The DAE formulation [START_REF] Guillot | Continuation of periodic solutions of various types of Delay Differential Equations using Asymptotic Numerical Method and Harmonic Balance Method[END_REF] gives the tension of the string without post-processing as it is an unknown of the equations. The ODE formulation [START_REF] Fritzsche | From holomorphic functions to complex manifolds[END_REF] can be used to compute the stability of the solutions obtained. Most of the time, this is the only way to have access to the stability informations of a system as a lot more results are known for ODE than for other types of systems. It seems that the DAE formulation is more robust when the solution gets closer to the homoclinic orbit, that is the limit where the period of the oscillations tends to infinity (and the energy, as defined in [START_REF] Jean | Periodic response investigation of large order non-linear rotordynamic systems using collocation[END_REF], to zero). This can be seen through the relative error on the period of the oscillations on the top of figure 1 that is everywhere smaller than for the other systems solved. The reason behind this is not clear yet for the authors.

Delay Differential Equations and fractional order systems

Two other types of systems have been studied using ANM and HBM. The delay differential systems [START_REF] Guillot | Continuation of periodic solutions of various types of Delay Differential Equations using Asymptotic Numerical Method and Harmonic Balance Method[END_REF] and the systems with fractional order derivatives [START_REF] Vigué | Continuation of periodic solutions for systems with fractional derivatives[END_REF]. To summarize these type of systems, a Van der Pol oscillator with feedback is explored

ẍ -ε(1 -x 2 ) ẋ + x = -kεF (x) (27)
where F (x) = x τ or D α (x) with D α the operator of Weyl fractional order derivation of order α and x τ (t) = x(t-τ ). To recast these equations quadratically, three auxiliary variables y = ẋ, r = 1-x 2 and F nl = ry-kF (x) are introduced. This gives :

ẏ + x + εF nl = 0 y -ẋ = 0 r -1 + x 2 = 0 F nl -kF (x) + ry = 0 (28)
Since equations ( 28) is quadratic according to the framework defined by equations [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], the results are obtained automatically according to the process detailed in section 1.3.2. The figure 3 shows two bifurcation diagrams computed with the system [START_REF] Khibnik | Continuation techniques and interactive software for bifurcation analysis of odes and iterated maps[END_REF], where the amplitude of the solution x of ( 27) is represented with respect to the continuation parameter which is τ or α depending on whether F (x) = x τ or F (x) = D α (x) respectively. The two bifurcation diagrams are completely different. There are two Hopf bifurcations 4 for the fractional order system and it seems that there are two vertical asymptotes at α = -1 and α = 3. It is remarkable that the method still gives the results for α negative. The value α = -1 allows to solve integrodifferential equations with our method. The delay differential system shows only one Hopf bifurcation in the area represented in figure 3. Even when τ is negative, the method gives some results which would not be possible to obtain with standard solvers as this represents a feedback from the future : the system is not causal anymore. Amplitude of x F(x) = x(t-). Continuation with respect to . F(x) = D (x). Continuation with respect to . Equilibrium Fig. 3 The solution branches of ( 27) with ε = 2 and k = 1.5 and for F (x) = x(tτ ) in blue and F (x) = D α (x) in red. The continuation parameter is respectively τ and α. The blue point and the red point with abscissa 0 coincide since both systems are equal :

x(t -0) = x(t) = D 0 (x)(t).
3 Equilibrium, periodic, quasi-periodic and transients regimes of an ODE

This section shows the different types of solutions that can be obtained with our method. For the sake of clarity, an academic example consisting in a system of coupled Van der Pol oscillators is studied in detail. This example has been used to illustrate the quasi-periodic Harmonic Balance Method (QPHBM) in [START_REF] Guillot | Continuation of quasi-periodic solutions with two-frequency harmonic balance method[END_REF]. The equations are

ẋ1 + a 1 ẋ1 + Ω 2 1 x 1 = (b 1 λ -c 1 λ(x 1 + x 2 ) -d 1 λ(x 1 + x 2 ) 2 )( ẋ1 + ẋ2 ) (29) ẋ2 + a 2 ẋ2 + Ω 2 2 x 2 = (b 2 λ -c 2 λ(x 1 + x 2 ) -d 2 λ(x 1 + x 2 ) 2 )( ẋ1 + ẋ2 ) ( 30
)
The value of the dimensionless parameters are

Ω 1 = 1, a 1 = 0.01, b 1 = 0.5, c 1 = d 1 = 2 ; Ω 2 = 2.5, a 2 = 0.025, b 2 = 1, c 2 = d 2 = 4.
The first order quadratic recast is the one derived in [START_REF] Guillot | Continuation of quasi-periodic solutions with two-frequency harmonic balance method[END_REF] :

ẋ1 = y 1 ẏ1 = -a 1 y 1 -Ω 2 1 x 1 + b 1 λ(y 1 + y 2 ) -c 1 λw -d 1 λv ẋ2 = y 2 ẏ2 = -a 2 y 2 -Ω 2 2 x 2 + b 2 λ(y 1 + y 2 ) -c 2 λw -d 2 λv 0 = r -(x 1 + x 2 ) 2 0 = w -(x 1 + x 2 )(y 1 + y 2 ) 0 = v -(x 1 + x 2 ) × w (31)
Now these quadratic equations (in the format [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] previously defined) will be used to continue equilibrium, periodic solutions, quasi-periodic solutions and transient regimes of the system. For equilibrium or periodic solutions, the stable solutions are represented in solid lines and the unstable solutions in dotted lines. Fig. 4 A bifurcation diagram of the equilibrium of system [START_REF] Marquis De Laplace | [END_REF]. The Hopf bifurcations are represented by red dots.

Equilibrium

In the quadratic recast [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF], the derivatives are replaced by 0 to obtain the system to solve to find the equilibrium. This system is obtained (easily) automatically from [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF]. It is clear here that the equilibrium is just

x 1 = y 1 = x 2 = y 2 = r = w = v = 0.
The stability of this equilibrium requires the computation of the eigenvalues of the Jacobian matrix of the first order system. This computation is automatized with the help of the condensation operation, detailed in [START_REF] Guillot | A generic and efficient taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF]), and gives two Hopf bifurcations for λ = 0.02 and λ ≃ 0.025 as shown in figure 4.

Periodic solutions

The quadratic recast [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF] is used for periodic solutions as explained in section 1.3.2. From the Hopf bifurcations computed in the previous section, two periodic orbits arise that can be continued. The starting point is computed automatically using the eigenvectors and the eigenvalues of the Jacobian matrix at the bifurcation point [START_REF] Beyn | Numerical continuation, and computation of normal forms[END_REF]. The stability of the periodic orbit is given by the Floquet exponents of the systems, computed with Hill's method as explained in [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF] and [START_REF] Bentvelsen | Modal and stability analysis of structures in periodic elastic states: application to the ziegler column[END_REF]. With this system two Neimark-Sacker bifurcations are detected, one on each branch, for λ ≃ 0.03 on the branch arising from the second Hopf bifurcation and for λ ≃ 0.35 on the branch arising from the first Hopf bifurcation. The bifurcation diagram obtained is represented on figure 5.

Quasi-periodic solutions

The quadratic recast ( 31) is also used here for quasi-periodic solutions. They are represented by double Fourier series [START_REF] Farner | Contribution to harmonic balance calculations of self-sustained periodic oscillations with focus on single-reed instruments[END_REF], using the Quasi-Periodic Harmonic Balance Method (QPHBM) to write the equations [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF] in the frequency domain. The QPHBM has been explained in [START_REF] Guillot | Continuation of quasi-periodic solutions with two-frequency harmonic balance method[END_REF] using this same example. In this same reference, the reader will find the explanation on how to obtain an algebraic system automatically to continue the quasi-periodic solutions. The bifurcation diagram obtained is represented on figure 6.

Only one branch of quasi-periodic solutions that connects the two bifurcations (see figure 6) can be computed from the Neimark-Sacker bifurcations. The starting point is computed automatically from the eigenvectors of Hill's matrix and the Floquet exponents at the bifurcation point (not detailed here). The question of the stability of a quasi-periodic orbit is a tough problem that is not investigated here.

Transient regimes

In this section, the system of equations ( 31) is discretized on time samples : 0, h, 2h, . . . , T -h, T where T = 400 and the integration step is h = T N -1 , N ∈ N. With this definition the number of time samples is N which is set to 2000 in the applications. The algebraic system is obtained by writing a Newmark scheme [START_REF] Newmark | A method of computation for structural dynamics[END_REF] in the way explained in section 1.3.1, that is, for all n ≥ 1 :

1 (nh) = x 1 ((n -1)h) + hy 1 ((n -1)h) + h 2 (( 1 2 -β)ẍ 1 ((n -1)h) + β ẍ1 (nh) y 1 (nh) = y 1 ((n -1)h) + h(1 -γ)ẍ 1 ((n -1)h) + hγ ẍ1 (nh) x 2 (nh) = x 2 ((n -1)h) + hy 2 ((n -1)h) + h 2 (( 1 2 -β)ẍ 2 ((n -1)h) + β ẍ2 (nh) y 2 (nh) = y 2 ((n -1)h) + h(1 -γ)ẍ 2 ((n -1)h) + hγ ẍ2 (nh) r(nh) = (x 1 (nh) + x 2 (nh)) 2 w(nh) = (x 1 (nh) + x 2 (nh))(y 1 (nh) + y 2 (nh)) v(nh) = (x 1 (nh) + x 2 (nh)) × w(nh) (32)
β and γ are the parameters of the scheme and will be set to 1 4 and 1 2 respectively for the computations. The second order derivatives are given by :

ẍ1 (nh) = -w 2 1 x 1 (nh) -a 1 y 1 (nh) + (b 1 λ)(y 1 (nh) + y 2 (nh)) -(c 1 λ)w(nh) -(d 1 λ)v(nh) ẍ2 (nh) = -w 2 2 x 2 (nh) -a 2 y 2 (nh) + (b 2 λ)(y 1 (nh) + y 2 (nh)) -(c 2 λ)w(nh) -(d 2 λ)v(nh) (33)
Initial conditions are added to the system [START_REF] Leimkuhler | Simulating hamiltonian dynamics[END_REF]. They are set to x 1 (0) = x 2 (0) = y 2 (0) = 0 and y 1 (0) = 0.01. It should be noticed that the system (32) can be automatically obtained thanks to the quadratic recast [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF]. The transient solutions are continued with respect to λ. Instead of a bifurcation diagram, the figure 7 represents all the transient solutions along the solution branch in a waterfall plot. For small values of λ 0.05, the oscillations are either evanescent λ < 0.02 or very slowly growing to a limit cycle that would need a longer time of integration to be reached. Between λ ≃ 0.05 and λ ≃ 0.2144, the solutions tend to the first limit cycle, which has a frequency similar to the resonance frequency of the first oscillator, Ω 1 . For λ 0.2144, the solutions tend to the second limit cycle, which has a frequency similar to resonance frequency of the second oscillator, Ω 2 . When λ ≃ 0.2144, the transients become very long and look very much like the quasi-periodic regimes found in the previous section. However, as this is only a timeintegration algorithm, one should not forget that these results depend on the time of integration and on the integration step. When the number of time samples is increased to N = 4000, 6000, 8000, this quasi-periodic transition occurs for λ ≃ 0.1780, 0.1720, 0.1700 respectively.

Phase space representation

The different types of solutions continued in the previous sections are represented altogether in 2-dimensional and 3-dimensional projections of the 4-dimensional phase space (x 1 , y 1 , x 2 , y 2 ) in figure 8. Two transient trajectories are represented together with the two stable limit cycles and the quasi-periodic orbit. The two transients were computed with slightly different initial conditions : for both, x 2 (0) = y 1 (0) = y 2 (0) = 0. The transient that converges towards the limit cycle at Ω 1 (in blue, in the center of the 3-D projection) has been computed with x 1 (0) = 10 -3 and the other one, in red, that converges towards the limit cycle at Ω 2 , with x 1 (0) = 10 -2 . As the 3-D projection is in the space (x 2 , y 2 , y 1 ), the two transients start from the same point, but this is of course not the case in the full 4-D phase space.

All the solutions represented in figure 8 have been computed using a unique code developed by the authors. The value of the continuation parameter is λ = 0.2144. Fig. 8 These figures show three different projections of the phase space for λ = 0.2144. The quasi-periodic orbit is in pink while the two limit cycles associated with the pulsation Ω 1 and Ω 2 are in blue and in red respectively. The transients in dark blue and dark red converge respectively to the blue and red limite cycles.

Conclusion and perspectives

This article shows how the quadratic recast of any smooth nonlinear system can be applied to continue its equilibrium, periodic or quasi-periodic solutions as well as its transients with the Asymptotic Numerical Method (ANM). The approach presented has been shown to be applicable to ordinary differential equations, delay differential equations, implicit differential algebraic equations and fractional order equations. The solutions of these systems are continued using the same implementation of the ANM. Therefore, further developments of the ANM benefit directly to all the continuation problems handled, requiring very few modifications. The academic examples proposed have illustrated the genericity of the framework described. An object-oriented implementation of the method has been developed5 under the name Manlab. The code Manlab tries to unify the treatment of the large variety of solutions and systems presented in this paper. From a user point of view as well as from a developer point of view, there is a minimum of work to do to switch from a type of system to another and from a type of solution to another. The crucial point is to provide a quadratic recast of the dynamical systems, most of the other operations have been automatized or could be automatized.

The perspectives of this work are many. The global continuation of an integration scheme offers the possibility to study its properties from a numerical point of view as can be seen in appendix B. The attractors of a system can be sought for by continuation with respect to the initial conditions. Chaotic areas could be investigated using the sensibility of the continuation to the initial conditions. The high robustness of the con-tinuation with the ANM is clearly an advantage to treat in a unified framework a wide range of problems coming from either purely algebraic equations or time domain and frequency domain discretizations of continuous equations. The comparison of several standard methods of representation of the solutions of a system is possible, allowing global comparison between direct time integration, harmonic balance, orthogonal polynomial collocation, trigonometric collocation, alternating frequency time scheme, etc... in the same framework.

A Algebra of the HBM applied on transcendental functions

The case of an auxiliary variable v defined as v(t) = g(u(t)) is treated in detail. Let us suppose that the auxiliary variable w(t) = ġ(u(t)) has been rewritten quadratically 6 , the variable v can be defined as the unique solution of the differential system

v(0) = g(u(0)) v(t) = w(t) u(t) (34)
which is the same as [START_REF] Doedel | Numerical computation of periodic solution branches and oscillatory dynamics of the stirred tank reactor with a yields b yields c reactions[END_REF]. The auxiliary variables u 0 = u(0), v 0 = v(0) and w 0 = w(0) are introduced to recast the first equation quadratically with the system :

v 0 = g(u 0 ) dv 0 = w 0 du 0 . ( 35
)
This system is treated with the formalism explained in [START_REF] Guillot | A generic and efficient taylor series based continuation method using a quadratic recast of smooth nonlinear systems[END_REF] briefly recalled in section 1.1.2. For the second equation, let us define the complex truncated Fourier series of u,v and w by :

u(t) = H h=-H U h exp(ihωt) , v(t) = H h=-H V h exp(ihωt) and w(t) = H h=-H W h exp(ihωt) (36)
The equation v(t) = w(t) u(t) becomes :

H h=-H ihωV h exp(ihωt) =   H h=-H W h exp(ihωt)     H h=-H ihωU h exp(ihωt)   (37)
which can be written

H h=-H ihωV h exp(ihωt) = 2H h=-2H     k 1 +k 2 =h -H≤k 1 ,k 2 ≤H W k 1 ik 2 ωU k 2     exp(ihωt). (38)
The right-hand-side series is truncated at order H and the harmonics are balanced :

∀h ∈ -H, H , ihωV h = k 1 +k 2 =h -H≤k 1 ,k 2 ≤H W k 1 ik 2 ωU k 2 . (39)
ω is the angular frequency of the periodic solution and thus remains positive. It can be simplified to obtain

∀h ∈ -H, H , hV h = k 1 +k 2 =h -H≤k 1 ,k 2 ≤H W k 1 k 2 U k 2 . ( 40
)
This expression is quadratic with respect to the Fourier coefficients. Hence the transcendental variable v can be written with a quadratic formalism in the frequency domain in an automated way. The operators V → dL(V) and (U, W) → dQ(W, U) are respectively the left hand side and the right hand side of the equations [START_REF] Seydel | From equilibrium to chaos: practical bifurcation and stability analysis[END_REF]. U, V and W are the vectors of Fourier coefficients of u, v and w respectively. It is easy to see that dL and dQ are respectively linear and bilinear in their argument(s).

B Time-integration of the pendulum with a class of Euler schemes

The ODE equations of the pendulum in polar coordinates (20)

θ = φ φ = -sin(θ) (41)
Fig. 9 The top figure shows a waterfall of the trajectories of ( 42). The bottom figure shows the evolution of the discrete energy (in color) with respect to time (y-axis) for all the values of λ (x-axis). From λ = 0 to λ = 1, the scheme evolves from an explicit Euler to an implicit Euler. The value λ = 0.5 is plotted in black to enhance the trajectory obtained with this scheme, which is known to be symplectic [START_REF] Leimkuhler | Simulating hamiltonian dynamics[END_REF]. It is clear in the bottom figure that the scheme either creates or dissipates energy when λ is respectively smaller or greater than 0.5.

are solved using a generalized Euler scheme on the discretized grid {0, h, 2h, . . . , (N -1)h}. The standard notation θn = θ(nh) and φn = φ(nh) is used. The initial conditions are θ 0 = π 2 and φ 0 = 0. The discretized system is then

θ n+1 = θn + h((1 -λ)φn + λφ n+1 ) φ n+1 = φn -h sin((1 -λ)θn + λθ n+1 ) (42)
where the parameter λ is the continuation parameter and interpolates an explicit Euler scheme (for λ = 0) and an implicit Euler scheme (for λ = 1). The centered Euler scheme (λ = 1 2 ) is known to be a simple symplectic scheme [START_REF] Leimkuhler | Simulating hamiltonian dynamics[END_REF]. It is then well suited to solve these conservative equations. The auxiliary variables ϑn = (1λ)θn + λθ n+1 and Fn = sin(ϑn) are added to the system [START_REF] Vigué | Continuation of periodic solutions for systems with fractional derivatives[END_REF] with the companion variables Gn = cos(ϑn). The final quadratic recast with the differentiated forms of the auxiliary variables is

∀n ∈ 0, N -1 ,          θ n+1 = θn + h((1 -λ)φn + λφ n+1 ) φ n+1 = φn -hFn ϑn = (1 -λ)θn + λθ n+1
Fn = sin(ϑn) and dFn = Gndϑn Gn = cos(ϑn) and dGn = -Fndϑn

The top figure 9 shows the trajectories obtained for N = 500 and h = 2π 100 , that is over around 5 periods of the analytic solution. The bottom figure 9 shows the discrete energy Hn = φ 2 n 2 -1 -cos(θn) in color for 0 ≤ λ ≤ 1. It is clear that energy is created when λ < 0.5 while energy is dissipated when λ > 0.5. The trajectory obtained for λ = 0.5 is represented in black showing the stability of the energy over time, as expected for this symplectic scheme.

C Computation of Taylor series coefficients

Some examples (taken from the section 1.1) of practical computation of the Taylor series coefficients are shown. However, this procedure is completely automatized from the quadratic formulation of the system of equations. It is detailed here for pedagogical purposes. The following computations illustrates how easy and straightforward it is to derive Taylor series coefficients when a quadratic recast is given.

C.1 Rational case

The system of equations is :

   y -u -xv = 0 ux -1 = 0 v -x 2 = 0 (44)
Assume that (x 0 , y 0 , u 0 , v 0 ) is a solution of the equations. The n-th Taylor series coefficient of the unknown x is written xn, such that : Developing the products and equating the terms of same order in a finally gives, for n ≥ 0:

   yn -un -n i=0 x i v n-i = 0 n i=0 u i x n-i -1δ n=0 = 0 vn -n i=0 x i x n-i = 0 (47)
where δ n=0 is the Kronecker symbol, its value is 1 if and only if n = 0. This system for n = 0 gives

   y 0 -u 0 -x 0 v 0 = 0 u 0 x 0 -1 = 0 v 0 -x 2 0 = 0 (48)
which is verified from the assumption that (x 0 , y 0 , u 0 , v 0 ) is a solution of the equations (44).

For n ≥ 1, the terms of higher order are put apart to make the recurrence relations appear :

   yn -un -x 0 vn -xnv 0 = n-1 i=1 x i v n-i u 0 xn + unx 0 = -n-1 i=1 u i x n-i vn -2x 0 xn = n-1 i=1 x i x n-i (49)
which is equivalent to the linear system

J 0 ×     xn yn un vn     =   n-1 i=1 x i v n-i -n-1 i=1 u i x n-i n-1 i=1 x i x n-i   ( 50
)
where J 0 =   -v 0 1 -1 x 0 u 0 0 x 0 0 -2x 0 0 0 1   is the Jacobian matrix of the system at the solution point (x 0 , y 0 , u 0 , v 0 ). From the knowledge of the initial solution point (x 0 , y 0 , u 0 , v 0 ) it is then possible to deduce the whole Taylor series of the variables x, y, u and v.

C.2 Transcendental case

The system of equations and its differential form are given by : y -sin(x) = 0 and dyzdx = 0 z -cos(x) = 0 and dz + ydx = 0 (51)

Assume that (x 0 , y 0 , z 0 ) is a solution of the equations. Thus y 0 = sin(x 0 ) and z 0 = cos(x 0 ). Contrarily to the previous section, the differential form of the equations is used to compute the higher order coefficients of the Taylor series by writing Developing the products and equating the terms of same order in a finally gives, for n ≥ 0:

(n + 1)y n+1 -n i=0 z i (n + 1i)x n+1-i = 0 (n + 1)z n+1 + n i=0 y i (n + 1i)x n+1-i = 0 (55)

The terms of higher order are put apart and everything is divided by (n + 1) to make the recurrence relations appear :

y n+1 -z 0 x n+1 = n i=1 z i n+1-i n+1 x n+1-i z n+1 + y 0 x n+1 = -n i=1 y i n+1-i n+1 x n+1-i (56)
which is equivalent to the linear systems

J 0 ×   x n+1 y n+1 z n+1   = n i=1 z i (n + 1 -i)x n+1-i -n i=1 y i (n + 1 -i)x n+1-i (57)
where J 0 = -z 0 1 0 y 0 0 1 is the Jacobian matrix of the system at the solution point (x 0 , y 0 , z 0 ). From the knowledge of the initial solution point (x 0 , y 0 , z 0 ) it is then possible to deduce the whole Taylor series of the variables x, y and z.

  sλφ s(0) = sin(θ(0)) and ṡ = c θ c(0) = cos(θ(0)) and ċ = -s θ (21)

Fig. 1 A

 1 Fig.1A bifurcation diagram of the pendulum with H = 500. The top figure represents the amplitude of the oscillations with respect to the energy. The bottom figure shows the relative error between the theoretical value (see text) and the computed value of the period of the oscillations in logarithmic scales. In blue the ODE system (21) is solved, in dashed red the DAE system[START_REF] Guillot | Continuation of periodic solutions of various types of Delay Differential Equations using Asymptotic Numerical Method and Harmonic Balance Method[END_REF] and in dotted yellow the IDE system[START_REF] Karkar | A comparative study of the harmonic balance method and the orthogonal collocation method on stiff nonlinear systems[END_REF].

Fig. 5 AFig. 6 A

 56 Fig.5A bifurcation diagram of the periodic solutions of the system[START_REF] Marquis De Laplace | [END_REF]. The x-axis is in logscale to help to distinguish the two branches arising from H1 and H2, which are the first and the second Hopf bifurcations shown on figure4. NS1 and NS2 are the two Neimark-Sacker bifurcations detected.

  x

Fig. 7

 7 Fig.7Two waterfalls of the transients of[START_REF] Marquis De Laplace | [END_REF]. The left figure represents the motion of the first oscillator (of resonance frequency Ω 1 ) and the right figure the motion of the second oscillator (of resonance frequency Ω 2 ). For small values of λ, the transients vanish towards the stable Equilibrium. For larger values of λ, the transients first converge towards a periodic orbit at pulsation around Ω 1 , then there is a narrow transition zone where the transients become very long and finally, for large values of λ the transients converge towards a periodic orbit at a pulsation around Ω 2 .

  by their Taylor series in the equations gives :

  ∂a of the Taylor series of the unknown x with respect to a is : by the corresponding series in the equations (52) gives :   +∞ n=0 (n + 1)y n+1 a n -+∞ n=0 zna n +∞ n=0 (n + 1)x n+1 a n = 0 +∞ n=0 (n + 1)z n+1 a n + +∞ n=0 yna n +∞ n=0 (n + 1)x n+1 a n = 0 (54)

Regular means that the Jacobian matrix J at U 0 is of full rank.

This is possible thanks to the analytic version of the implicit function theorem[START_REF] Fritzsche | From holomorphic functions to complex manifolds[END_REF].

In the case of a system with delay or fractional order derivatives, it may be required to give the initial conditions over an interval of the form [-T, 0].

Following[START_REF] Vigué | Continuation of periodic solutions for systems with fractional derivatives[END_REF], we use here the name Hopf bifurcation. To our knowledge, this type of bifurcation has not been characterized theoretically but is analogous to a Hopf bifurcation point from a numerical point of view.

Manlab is available online on the dedicated website https://manlab.lma.cnrs-mrs.fr/ .

If this is not the case, an auxiliary variable z(t) = ẇ(t) is added and the procedure explained here is applied recursively.

Acknowledgements This work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the Investissements d'Avenir French Government program managed by the French National Research Agency (ANR). Conflict of Interest: The authors declare that they have no conflict of interest.