R Kalvig

E Jedryka

P Aleshkevych

M Wojcik

W Bednarski

M Petit

L Michez

Ferromagnetic resonance in Mn 5 Ge 3 epitaxial films with weak stripe domain structure

Extensive X-band and Q-band FMR experiments have been performed in the Mn5Ge3 epitaxial films with thicknesses varying between 4.5 and 68 nm. FMR signals were recorded in the temperature range between 15 and 295 K, at different orientations of magnetic field with respect to the film plane. In addition to the acoustic FMR mode with well defined resonance field, originating from inside the magnetic domains, a broad absorption line has been observed at low fields and attributed to the unresolved spectrum of FMR modes having the origin in flux closure caps. The FMR results have been discussed in the context of the domain structure computed with the use of OOMMF micromagnetic calculations and giving good agreement with the experimental hysteresis curves. From the Q-band experiments, where the FMR signal is observed in magnetically saturated sample, the uniaxial anisotropy constant in films with different thicknesses has been determined as a function of temperature. This FMR study provides the evidence that the strong uniaxial anisotropy observed in epitaxial thin films of Mn5Ge3 leads to the formation of a stripe domain structure above 25 nm, in agreement with the published reports on magnetization studies in these films. It also eliminates a possible confusion that may arise from previously published FMR studies on films grown with the same method, which led their authors to conclude that the shape anisotropy can force the magnetization to the in-plane orientaion in this thickness range and even above it.

I. INTRODUCTION

Ferromagnetic compound Mn 5 Ge 3 is of interest for spintronic applications, due to considerable spin polarization (42%) and high Curie temperature (296 K) which can be further increased (up to 445 K) by the addition of carbon [1,2]. Importantly, this compound is compatible with the mainstream silicon technology and can be epitaxially grown on Ge, maintaining its magnetic properties and metallic conductivity [3,4]. Moreover, it displays strong magnetocrystalline anisotropy with hexagonal caxis as an easy direction, which makes it prospective for spintronic and magnetic recording applications [5].

It has been shown that this compound can be successfully grown by Solid Phase Epitaxy (SPE), consisting in a room temperature Mn deposition on the Ge(111) substrate followed by thermal annealing. Good quality epitaxial films can be obtained with (0001) hexagonal c-axis parallel to the Ge(111) direction [6]. A detailed analysis of magnetization reversal process in the SPE-prepared thin films of Mn 5 Ge 3 revealed that in the thickness range between 10 nm and 20 nm, magnetization direction at low temperatures undergoes reorientation from the inplane to out-of-plane direction (i.e. along the hexagonal c-axis) and at 25 nm a clear signature of a multidomain structure with perpendicular orientation was observed [5,7].

Meanwhile, another study in Mn 5 Ge 3 films prepared with the same technique [8], has reported-based on the analysis of Ferromagnetic Resonance (FMR) experiments-that in spite of the demonstrated uniaxial anisotropy perpendicular to the film plane, magnetization is oriented in sample plane, within the investi-gated range of temperatures (30-290 K), even in samples as thick as 30 nm. This evident discrepancy requires a closer examination, especially since the FMR experiments reported in [8] were performed at 9.08 GHz (Xband), and the corresponding FMR signals were observed at magnetic field values below sample saturation, therefore some kind of domain structure can be expected. To resolve these discrepancies and to get a deeper insight into the magnetic properties of epitaxial thin films prepared by SPE method, we have undertaken an extensive FMR study in the Q-band, where the FMR resonance takes place in a fully saturated sample, as well as in the X-band (in presence of a domain structure). Mn 5 Ge 3 films with thicknesses varying between 4.5 and 68 nm have been studied as a function of temperature and sample orientation. Our results show unambiguously that in the investigated film thickness range between 18 and 68 nm the uniaxial anisotropy prevails over the shape anisotropy leading to the stripe domain structure with closure caps containing a significant fraction of magnetic moments.

Analysis of our FMR data, supported by a simulation of the magnetic moments distribution using the OOMMF software package [9], made it possible to identify particular fragments of the domain structure giving rise to the resonance modes observed at X-band: on one hand an acoustic FMR mode with well-defined resonance field originating from magnetization inside the stripe domains, and on the other hand-a broad absorption spectrum observed at low fields, attributed to flux closure caps. The Q-band experiments, where the uniform precession mode is observed in magnetically saturated sample, made it possible to determine the uniaxial anisotropy constant as a function of temperature in films with different thicknesses.

II. EXPERIMENTAL DETAILS

Thin films of Mn 5 Ge 3 with thickness of 4.5 nm, 9 nm, 18 nm and 68 nm have been prepared on the Ge(111) substrates (500 µm thick) by solid phase epitaxy (SPE) in molecular beam epitaxy (MBE) chamber with a base pressure less than 10 -10 Torr. The SPE preparation consists of a room temperature deposition of Mn layer followed by thermal annealing at 720 K to activate interdiffusion between the deposited metal and the germanium buffer layer. The as prepared Mn 5 Ge 3 films have been capped with amorphous Ge to prevent oxidation. Further details on sample preparation are given in ref [6,10]. The quality of crystal growth has been monitored in situ using reflection high-energy electron diffraction (RHEED). The obtained samples have been characterized by XRD technique-the scan clearly indicates two narrow reflections (0002) and (0004) evidencing good quality Mn 5 Ge 3 phase. The XRD results confirm also that c-axis of Mn 5 Ge 3 is parallel to Ge(111) [7].

Magnetic characterization was performed using a Quantum Design superconducting quantum interference device magnetometer (SQUID). In order to understand better the magnetic properties and the evolution of a domain structure postulated from these measurements, micromagnetic simulations of the magnetic moment distribution within the sample have been performed with the use of the OOMMF software package [9] and compared with the experimental data.

Ferromagnetic resonance experiments at X-band (9.38 GHz) were performed by means of a commercial Bruker EMX electron spin resonance spectrometer with a rectangular microwave cavity working in the TE 102 mode and a standard phase sensitive detection. External magnetic field was modulated at 100 kHz. The temperature has been varied in the range 15-295 K with the use of an Oxford helium-flow cryostat. The samples were mounted on a quartz holder and the measurements were performed with magnetic field applied in the plane of a sample and perpendicular to it, whereas microwave field was applied parallel to the sample plane and perpendicular to the DC bias field. Due to a strong uniaxial anisotropy of the studied compounds, the resonance field in parallel orientation of the DC magnetic field was below the saturation field in the X-band, therefore additional measurements have been performed in the Q-band (34 GHz) with the use of BrukerElexSys E500 spectrometer equipped with flexline ER 50106QT/W resonator and ER 4118CF cryostat operating at 3.8-295 K. Temperature was controlled and stabilized with the Oxford Instruments temperature controller ITC 503S. FMR spectra were recorded as the first derivative of the microwave power (about 2 mW) absorption versus the external magnetic field in the range of 0-17 kOe. Magnetic field was modulated with the fre-quency of 100 kHz and modulation amplitude of the order of 10 Oe. The respective configuration of the pumping field versus DC bias magnetic field was the same as in case of the X-band experiment. Manual homemade goniometer was used in the studies of angular dependences of the resonance fields.

III. RESULTS AND DISCUSSION

A. Magnetization reversal and micromagnetic OOMMF simulations of the domain structure A thorough study of magnetization reversal at low temperatures in Mn 5 Ge 3 thin films used in this study has been reported in Ref. [5,7], as a function of film thickness. In case of magnetic field applied in the film plane, hysteresis loops recorded in very thin films have a square shape indicating that the magnetization follows the external field, and above 16 nm the loops become canted. As an example, in Fig. 1 we present the magnetic hysteresis loop measured at 15 K in 25 nm film. Magnetization reaches its saturation value of 1450 emu/cm 3 for the in-plane fields higher than 6.3 kOe. We also note a coercive field of 500 Oe and remanence magnetization corresponding to about 20% of saturation value, indicating the presence of an in-plane component of magnetization. In the same figure we plot the theoretical curve computed for a 30 nm thick film using the OOMMF package, as described below. Magnetization curves recorded from the same sample with magnetic field applied perpendicular to the film plane, are shown in Fig. 2. Magnetic saturation in the out-of-plane configuration is reached around 10.3 kOe.

With decreasing magnetic field the system reveals an inertia: magnetization maintains its saturation value, and then at around 8.1 kOe it suddenly drops down. This is a well-known effect due to the fact that the oppositely magnetized domain has to reach a certain minimum dimension in order to create a stable domain configuration. Therefore the magnetization curve in the "fielddown" direction does not follow the same values as in the "field-up" direction, creating a small loop just below the saturation field. This loop is better visible in the thicker samples, as shown in Fig. 2 inset, where we present the blow-up of the high field detail of the hysteresis loop recorded from the 68 nm thick film (note that the bubble nucleation region in case of 68 nm film is between 9.1 and 9.7 kOe). The curves corresponding to other films from the same series of samples, revealing the small loop feature, have been presented in the Fig. 3 of Ref. [5]. Similar small loops have also been reported in other films with stripe domain structure, e.g. ( 0001)hcp cobalt films [11] or thin films of Co 2 MnGa Heusler alloy [12]. The material parameters obtained from magnetization studies were used in Ref. [7] to determine the ratio between the uniaxial perpendicular magnetocrystalline anisotropy and the shape anisotropy of the thin film, described by the so called quality factor, defined as Q=K u /(2πM s 2 ), which was found in the studied materials to be around 0.6 at 15 K [7]. Thin films of the low Q-factor materials (Q<1), do not develop an open stripe domain structure, as is the case of materials with Q>1, but rather tend to form a flux closure caps at the film surface. Therefore in addition to the wide Bloch-type domain walls within the film volume, one can expect at the film surface a fraction of magnetic moments with the in-plane orientation, separated by the Neel-type domain walls. Interpretation of the FMR data in such a complex domain structure requires a detailed consideration of magnetic moment distribution. To this end, we performed micromagnetic calculations with the use of freely available software package OOMMF [9]. Simulated volume was a cuboid with dimensions 150 nm x 150 nm x 30 nm, whereas the discretization cell was cube of side 1.5 nm. Exchange stiffness parameter A was taken as A=1•10 -7 erg/cm from Ref. [5] and uniaxial anisotropy constant K u =5.66•10 6 erg/cm 3 was taken from the Q-band FMR experiment as described in the next section. Saturation magnetization was estimated as 1450 emu/cm 3 , based on the experimental hysteresis loops. Calculations were carried out by minimizing the total energy in each discretization cell, while the external magnetic field was applied with 200 Oe step, from 0 field to 10 kOe in sample plane (direction x on Fig. 3), and then in the opposite direction (-x) to -10 kOe and back to 10 kOe. Similar procedure was applied for magnetic field perpendicular to sample plane (along z-axis).

Magnetic hysteresis loops computed in the in-plane and out-of-plane configuration reproduce correctly the features of experimental curves, as shown in Fig. 1 and Fig. 2, respectively. The magnetization curve calculated for field applied in sample plane (Fig. 1) is almost identical to the experimental data, with saturation field around 6.3 kOe and remanence magnetization of about 20% of saturation value. The curve computed for field oriented perpendicular to sample plane reveals a different value of domain nucleation field and much broader loops around saturation than the experimental data recorded for 25 nm film. This is probably due to the fact that in a real material there may be a certain distribution of demagnetizing fields due to the miniscule variations in surface roughness and to the influence of a capping layer. Nevertheless, the general character of the magnetization reversal is reproduced, proving that the calculated magnetic structure represents closely the magnetic moment distribution in a real material.

Fig. 3 and Fig. 4 visualizes the magnetic domain structure resulting from the micromagnetic calculations, as a function of magnetic field applied along the x-axis and z-axis, respectively (the coordinate system used is shown in panels A). The respective panel columns visualize the magnetic moment distribution in 30 nm Mn 5 Ge 3 film in the three projections: (1) Top view at sample surface, (2) x-y cross section at half-thickness and (3) y-z cross section. The rows A, B, C, D, E correspond to characteristic field values, as indicated on the hysteresis loops in Fig. 1 and Fig. 2. Point A in Fig. 1 corresponds to initial magnetic equilibrium state in this material before turning the magnetic field on. In panels A of Fig. 3 we note a clearly visible domain structure with two oppositely oriented domains, up and down along the c-axis separated by the Bloch-type walls. White color and blue arrows indicate magnetic moments pointing in the "+z" direction (domains "up"), whereas black color and red arrows describe domains oriented "down" in the "-z" direction. Green background and black arrows represent fractions of the sample volume where magnetic moments are oriented in sample plane. Flux closure caps on the film surface as well as Bloch wall centers in the sample volume are the main contributors to the in-plane magnetization (remanence around zero field). With increasing in-plane field value, i.e. going towards point B in Fig. 1, the magnetization inside domains tilts towards the field direction, forming a characteristic stripe domain pattern, until a full saturation in the film plane (point C) is reached. Upon subsequent decreasing the field value, the stripe domain pattern is recovered and it is preserved when the field is switched off (point D). Application of magnetic field in the opposite direction reverts the spins on the surface as shown in the panel E1 in Fig. 3, and magnetic moments inside the domains also begin to tilt towards the "-x" direction. The evolution of a domain structure in case of magnetic field applied in the direction perpendicular to the film is illustrated in Fig. 4, using the same convention of colors. Starting from the same domain pattern (panel A), which corresponds to the initial point A, we observe a progressive shrinking of the domains with magnetization oriented opposite to the applied field. At point B (5 kOe) only two cylindrical domains with magnetization in the "-z" direction remain, and at point C the sample is fully saturated in the "+z" direction. Point D is the characteristic point on the hysteresis curve where a sudden drop of magnetization from its saturation value takes place. At this field value we note the onset of a number of well-defined bubble domains magnetized in "-z" direction (panel D). Finally, at point E, in the absence of magnetic field, the equilibrium between "up" and "down" domains is recovered. In conclusion, the micromagnetic calculations performed in Mn 5 Ge 3 films with thickness of 30 nm reproduce correctly the experimental magnetization curves and provide the evidence for presence of a stripe domain structure typical for the low Q-factor materials, at least above 25 nm. Considering the systematic evolution of hysteresis curves as a function of film thickness reported in Ref. [5] one can extend this conclusion at least up to the thickest sample in this study (68 nm) where regular stripe domain structure is expected.

B. FMR studies

Fig. 5 presents the absorption derivative X-band FMR spectra recorded at 100 K from films with thicknesses between 4.5 and 68 nm, for the two orientations of DC magnetic field: in the film plane and perpendicular to it. All the films reveal two characteristic spectrum components: a narrow, well resolved resonance line and a very broad absorption signal at low field values. In case of 4.5 nm film this broad signal represents almost entire spectrum intensity, the narrow mode is merely marked.

A remarkable observation is that in case of 68 nm sample, the narrow FMR signal in the out-of-plane configuration could only be observed when the magnetic field was swept downwards after the sample was saturated in the out-of-plane direction. When starting from a demagnetized state and increasing the field all the way up to saturation, only the broad absorption at low fields was recorded, as shown in Fig. 6 where we present the X-band FMR spectra recorded in this sample at different temperatures, while sweeping the field upwards (red dashed line) and downwards (black line).

FIG. 5. Absorption derivative FMR spectra recorded from

Mn 5 Ge 3 films with thickness between 4.5 nm and 68 nm. The experiment was performed at X-band (9.38 GHz) at 100 K with DC field applied parallel (left) and perpendicular (right) to the film plane.

To understand better the origin of the observed FMR signals, we performed additional experiments at 34 GHz (Q-band). The absorption derivative FMR spectra recorded at different temperatures from the same 68 nm thick Mn 5 Ge 3 film at Q-band are presented in Fig. 7. Similar to the X-band experiments, the spectra recorded at Q-band consist of a narrow, well resolved resonance line and a very broad absorption signal at low field values. In case of Q-band experiment with the DC field applied perpendicular to the film plane, the resonance field of the narrow line in 68 nm sample is beyond the available field range at low temperatures-it could however be retrieved at 295 K, as shown in Fig. 7.

To identify the origin of the narrow resonance mode it is useful to consider the position of the resonance field as a function of film thickness. Fig. 8 summarizes the resonance field values of the well resolved FMR line (X-band and Q-band experiment) in case of the external magnetic field oriented in-plane and out-of-plane at 15 K. It also shows the respective saturation field values determined from the magnetization measurements. Vertical dashed lines separate three thickness regions where magnetization reversal studies indicate its equilibrium orientation in the film plane (below 10 nm), perpendicular to it (above 20 nm), and the reorientation region inbetween [7].

It is well established that the FMR spectra in thin magnetic films with stripe domains can have a complex structure, consisting of a number of resonance modes. In the simplest case of high-Q material, the resonances in the adjacent "up" and "down" domains, separated by the Bloch walls, are coupled via the dipolar demagnetization fields, splitting the FMR signal in two branches: in-phase mode ω + (denoted also as acoustic) and out-of-phase ω - mode (optical), the latter observed at higher field values [13,14]. Theoretical calculations of Vukadinovic et al. [15] have shown the dependence of the ω + and ω - frequencies in high-Q materials on magnetic parameters such as quality factor, demagnetization, and damping parameter, giving a satisfactory agreement with the experimental results on bismuth-substituted single crystal garnet films. In addition to these two basic FMR modes a third absorption mode, due to resonance of domain wall oscillations (DWR) has been reported in the FMR experiments in garnet thin films [16,17]. This mode is excited by the r.f. field component parallel to magnetization within the domain and can overlap with the FMR modes. An exhaustive review of FMR studies in high- Q films and hybridized FMR-DWR oscillations has been given by Vukadinovic [18].

In the low-Q magnetic films with perpendicular anisotropy the situation is more complex: below the critical thickness the material displays the in-plane orientation of magnetization, forming large in-plane domains separated by the Neel-type walls. Above a critical thickness these films reveal stripe domain pattern, which is referred to as "weak stripe domains", characterized by the presence of flux closure caps. The dynamic response to high frequency magnetic excitations of such low-Q magnetic thin films has been studied by FMR in hexagonal Co [19,20], FePd [18,21], Co 2 MnGa [12] and by zerofield microwave susceptibility in several magnetic alloys, e.g. CoNbZr [22], CoFeZr [23], CoFeB [24] CoZrTa [25], and permalloy [26].

These studies, involving micromagnetic calculations and comparison with the experimental results, have shown that a number of different FMR modes could be excited, depending on the respective orientations of the pumping field and the DC bias field enforcing stripe orientation [19][20][21]27]. In the particular configuration of microwave field in plane of the sample and perpendicular to the DC bias field (denoted as configuration "X") the selective excitation of two kinds of FMR modes in Co thin films has been demonstrated by Ebels [19,20]. One is the acoustic FMR mode originating from magnetization inside the domains and the second set consists of two resonance modes related with the flux closure caps (cap breathing mode) [19,20]. Vukadinovic presented the extended model, calculating the theoretical in-plane FMR spectra from the frequency and magnetic field dependence of magnetic susceptibility. These calculations, applied to FePd films with Q=0.4, confirmed that in the "X" configuration, the FMR spectra are dominated by the acoustic mode, having its continuation in the uniform mode above saturation [21].

The micromagnetic simulations of the domain structure presented in the previous section show that our samples (Q=0.58) display a similar weak stripe domain pattern as Co films (Q≈0.5) studied in [19,20] and FePd films (Q=0.4) studied in [18,21]. Our experimental configuration (microwave field in plane of the sample and always perpendicular to the DC bias field) corresponds to the configuration "X" in Ref. [19,21] Therefore, by analogy to the line assignment presented by Ebels [19], we attribute the broad absorption signal that we systematically observe at low fields in the X-band and Q-band experiments to the unresolved spectrum of excitation modes of magnetic moments located in different regions of flux closure caps. The volume involved in flux closure caps in our Mn 5 Ge 3 films is clearly larger than in case of hexagonal Co, resulting in a broad absorption, rather than clearly resolved two peaks reported in [19,20]. A similar broad absorption signal was observed at low fields in the FMR experiments on thin films of Heusler alloy Co 2 MnGa with stripe domain structure [12]. It is interesting to note that this signal has not been reported in previous FMR studies on Mn 5 Ge 3 thin films [8,28], even though they present FMR results in the unsaturated state.

The second component of our FMR spectra presented in Fig. 6 (X-band experiment in 68 nm Mn 5 Ge 3 film)i.e. a narrow, well-resolved resonance line-is interpreted as the acoustic mode of resonance excitations inside the domains, in agreement with the Ebels and Vukadinovic models for the "X" pumping configuration. This assignment is supported by the interesting feature presented in Fig. 6: the FMR signal is observed only when sweeping the DC field down after reaching an out-of-plane saturation, and the resonance field corresponding to this line (9.3 kOe) coincides with the domain nucleation region on hysteresis loop in this sample (inset in the Fig. 2). The domain structure corresponding to the point D (shown in the Fig. 4) consists of a number of well-defined cylindrical domains with magnetization opposite to the external field. This result shows clearly that the observation of this resonance mode is conditioned by the onset of a stable structure of oppositely oriented domains. The acoustic mode corresponds to the situation where the magnetization in neighboring, oppositely magnetized domains rotates in phase, coupled via the surface charges. The "optical" FMR mode is not observed in the present experiment, since it requires pumping along the DC field direction ("Y" configuration in Ebels notation).

In case of the Q-band experiments the resonance fields are clearly higher than the sample saturation field in both configurations the DC magnetic field. The resonance takes place in magnetically saturated state, all spins are aligned in the DC field direction and this line can be unambiguously interpreted as the uniform precession mode. In order to derive the resonance equations, the commonly accepted free energy approach, developed originally by Smit and Beljers [13] and extended by Baselgia et al. [29], was used here. The anisotropic part of magnetic free energy density F includes the Zeeman energy, the demagnetizing energy and uniaxial K u anisotropy energy density: F=-M H+(2πM 2 -K u )cos 2 θ. The resonant condition at a fixed microwave frequency ω is given by

ω γ 2 = 1 M 2 ∂ 2 F ∂θ 2 1 sin 2 θ ∂ 2 F ∂φ 2 + cosθ sinθ ∂F ∂θ - 1 sinθ ∂ 2 F ∂θ∂φ + cosθ sin 2 θ ∂F ∂φ 2 (1)
where θ and φ are the polar and azimuthal angle of magnetization, respectively. The θ is counted from the film normal. The resonance field is obtained by evaluating Eq.( 1) at the equilibrium position of M ( ∂F/∂θ=0 and ∂F/∂φ=0) for a given orientation of H. When magnetic field is applied perpendicular or parallel to the film surface, the Eq.( 1) is reduced to a simple form known as classical Kittel's FMR equations [30]:

ω γ = H ⊥ -4πM S + 2K u M S (2)
in case of magnetic field applied perpendicular to sample plane and

ω γ 2 = H H + 4πM S - 2K u M S (3)
in case of magnetic field applied in sample plane.

The physical quantities H , H ⊥ are resonance field values in the respective configurations, M S -saturation magnetization, K u -uniaxial anisotropy constant. In case of 68 nm thick sample, the resonance field at 295 K was H = 11.1 kOe and H ⊥ = 14.4 kOe, whereas saturation magnetization M S = 250 emu/cm 3 . Using these values we can determine the anisotropy constant of Mn 5 Ge 3 at 295 K to be: K u =11•10 4 erg/cm 3 . This result is comparable with the value previously determined in bulk Mn 5 Ge 3 samples (30•10 4 erg/cm 3 ) [31]. The difference may be related to a From Ref. [31] at 77 K the fact that these values refer to the temperature bordering the Curie point.

In the 9 nm and 18 nm thick films the uniform precession signal could be observed in the Q-band at all studied temperatures, making it possible to calculate the K u values as a function of temperature for these films, as shown in Table 1. The values obtained for the 9 nm and 18 nm films are very close to each other at all temperatures and almost identical to the reported value for bulk at 77 K [31]. The scatter of K u values observed at room temperature ist most probably related to the vicinity of the Curie point, as already mentioned. When rotating the DC magnetic field between the outof-plane and in-plane direction of magnetic field one can trace the position of the acoustic mode, as shown in Fig. 9, where we plot the results of X-band experiments performed in 68 nm sample at 15 K and at 295 K. We note that at 295 K the position of resonance line is above saturation field in case of the in-plane orientation (H sat= 2 kOe) [32], and thus the observed resonance is the uniform precession mode. At the same time, for orientations close to the film normal the resonance field is still below saturation (H sat=7 kOe) [32] and the observed line is an acoustic mode. This is in line with the theoretical predictions that the acoustic mode transforms into uniform precession in the magnetically saturated state [20]. An attempt to fit the experimental data at 295 K using the Eq. ( 1), gave an excellent result for the following free parameters: K u =9.1•10 4 erg/cm 3 and g-factor g=1.99. A slight difference between this result and the K u value obtained from the Q-band experiment (11•10 4 erg/cm 3 ) is attributed to the fact that in the out-of-plane configuration the sample is not fully saturated. Nevertheless, closeness of these values and the characteristic shape of the resonance field versus angle dependence confirm the correctness of our interpretation. In case of the data obtained at 15 K such fit would not have a physical meaning, since the Eq. ( 1) applies to the saturated state, and not to the acoustic mode. The solid line for the data at 15 K shown in Fig. 9 is only a guide to the eye, underlining the characteristic shape of the resonance field versus angle dependence.

The resonance field in the in-plane configuration of the acoustic mode is around 1 kOe at 15 K. The domain structure at this field is expected to resemble structures visualized in panel D in Fig. 3, i.e. domain magnetization is tilted in the field direction, and close to the surface it is oriented in plane. With decreasing film thickness the proportion of in-plane oriented moments increases, as revealed by the huge intensity of the broad absorption line (see Fig. 5), but the acoustic mode is still traceable, even in films as thin as 4.5 nm.

IV. CONCLUSION

Extensive X-band and Q-band FMR experiments have been carried out in the Mn 5 Ge 3 epitaxial films with thicknesses varying between 4.5 and 68 nm. The FMR signals were recorded in the temperature range between 15 and 295 K, at different orientations of magnetic field with respect to the film plane. In parallel, we performed the micromagnetic calculations of the domain structure of a Mn 5 Ge 3 film 30 nm thick as a function of the magnetic field applied in the film plane and perpendicular to it. The theoretically calculated magnetization curves are fully consistent with the experimental ones, proving the correctness of the simulated domain structure generated within the same model. Analysis of the experimental FMR spectra in the context of the theoreticaly simulated domain structure, made it possible to assign the observed resonance modes to particular fragments of the domain structure: on one hand magnetization precession inside the stripe domain gives rise to an acoustic FMR mode with well-defined resonance field. On the other hand, a relatively large volume of magnetization inside closure caps, gives rise to the broad absorption spectrum of unresolved FMR lines, observed at low fields. Our results show unambiguously that in the film thickness range above 18 nm the uniaxial anisotropy leads to the stripe domain structure with closure caps containing a significant fraction of magnetic moments, supporting the results reported in Ref. [5]. This conclusion negates the interpretation of FMR data in Mn 5 Ge 3 films presented in Ref. [8], where the authors claim that magnetization is entirely oriented in the film plane, and interpret the acoustic resonance mode as the uniform FMR mode. Our results stress the importance of combining the analysis of the FMR modes observed in magnetically non saturated materials with micromagnetic calculations and simulations of a domain structure. The effective uniaxial anisotropy constant in 18 nm film was found to vary between 5.66•10 6 erg/cm 3 at 15 K to 4.73•10 6 erg/cm 3 at 100 K in perfect agreement the reported values for the bulk material. The room temperature values of K u reveal a significant scatter which is attributed to the vicinity of the Curie point.

FIG. 1 .

 1 FIG. 1. (Color online) Solid-dotted line: experimental hysteresis loop recorded with SQUID magnetometer at 15 K in Mn 5 Ge 3 epitaxial film (25 nm), with magnetic field swept in the film plane. Solid line: theoretical curve computed with the use of OOMMF software. Letters A-E indicate characteristic points giving rise to the domain structure presented in diagrams in Fig.3.

FIG. 2 .

 2 FIG. 2. (Color online) Solid-dotted line: experimental hysteresis loop recorded with SQUID magnetometer at 15 K in Mn 5 Ge 3 epitaxial film (25 nm), with magnetic field swept out-of-plane. Solid line: theoretical curve computed with the use of OOMMF software. Letters A-E indicate characteristic points giving rise to the domain structure presented in diagrams in Fig. 4. The inset shows the blow-up of the high-field section of the experimental curve for 68 nm thick sample, indicating the field region corresponding to bubble nucleation.

FIG. 3 .

 3 FIG. 3. (Color online) Diagrams showing magnetic moment distribution in different parts of 30 nm Mn 5 Ge 3 film, computed with the use of OOMMF software package. Magnetic field applied in the film plane, rows (A,B,...E) correspond to characteristic points of hysteresis loop in Fig. 1. Columns represent different projections: [1]: Top view at sample surface, [2]: x-y cross section at half-thickness and [3]: y-z cross section. For details see the text.

FIG. 4 .

 4 FIG. 4. (Color online) Diagrams showing magnetic moment distribution in different parts of 30 nm Mn 5 Ge 3 film, computed with the use of OOMMF software package. Magnetic field applied perpendicular to the film plane, rows (A,B,...E) correspond to characteristic points of hysteresis loop in Fig. 2. Columns represent different projections: [1]: Top view at sample surface, [2]: x-y cross section at half-thickness and [3]: y-z cross section. For details see the text.

FIG. 6 .

 6 FIG. 6. (Color online) Absorption derivative FMR spectra recorded from the 68 nm thick Mn 5 Ge 3 film at X-band (9.38 GHz) at temperatures between 15 K and 295 K. Left column: magnetic field applied in the film plane, right: magnetic field applied perpendicular to the film plane. Red dashed line (shifted along y axis to ensure visibility): sweeping upwards the magnetic field (note the absence of narrow FMR mode); black line: sweep-down. For details see the text.

FIG. 7 .

 7 FIG. 7.Absorption derivative FMR spectra recorded from the 68 nm thick Mn 5 Ge 3 film at Q-band (34 GHz) at temperatures between 15 K and 295 K. Left column: magnetic field applied in the film plane, right: magnetic field applied perpendicular to the film plane.

FIG. 8 .

 8 FIG. 8. (Color online) Magnetic saturation field values (in-plane and out-of-plane) and the respective resonance field values in X-band and Q-band FMR experiment at 15 K, as a function of film thickness. Vertical dashed lines indicated the magnetization reorientation region reported in Ref. [7].

FIG. 9 .

 9 FIG. 9. (Color online) X-band FMR resonance field recorded in Mn 5 Ge 3 film 68 nm thick at 15 K and 295 K as a function of angle between the external magnetic field and the film normal (0 deg.field along normal, 90 deg.-field in sample plane). Blue dashed line for acoustic mode at 15 K serves as a guide to the eye. Solid black line at 295 K represents the fit as described in the text.

TABLE I .

 I Uniaxial anisotropy constant (erg/cm 3 ) as a function of temperature determined from the FMR Q-band experiment in a series of Mn5Ge3 films .

	Thickness [nm]		Temperature [K]
		15	50	100	200	295
	9	5.41•10 6 5.06•10 6 4.69•10 6 3.22•10 6 2.65•10 4
	18	5.66•10 6 5.27•10 6 4.73•10 6 3.16•10 6 2.74•10 4
	68					11•10 4
	Bulk			4.2•10 6a		30•10


Kevin Schmidmayer
email: kevinsch@caltech.edu

Antoine Marty

Fabien Petitpas

Eric Daniel

Eric Daniel Ecogen

ECOGEN, an open-source tool dedicated to multiphase compressible multiphysics flows

This paper presents a new multiphase flow code, cast under an open-source GNU license. The main characteristics of the different flow models are given, then the numerical method used is briefly presented: it includes temporal flow solvers, meshing features (like AMR technics), results visualization. Two examples of flows solutions are presented: the interaction of a high-speed flow with a droplet and the second concerns the attenuation of a propagating shock wave.

INTRODUCTION

Many numerical flow models have been developed devoted to rather different flows that concern multiphase flows (1), ( 2), (3), (4), [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF], [START_REF] Schmidmayer | A model and numerical method for compressible flows with capillary effects[END_REF]. It appeared necessary to develop a stable numerical platform that includes all those models without losing efficiency in programming and make easier the diffusion in the scientific community. ECOGEN is the tool devoted to this task. It has been developed and distributed under an open-source (GNU GPLv3 license). This code is available online at the following address: https://github.com/Matshishkapeu/ECOGEN In this paper we present the basic of the multiphase compressible flow models include in ECOGEN. The main features of the codeconsidering the temporal solver, the new AMR technics, the domain decompositionare present. Then, the results concerning two high speed flow results are given: the first one concerns the atomization of a liquid droplet and the second is about the attenuation of a propagating shock wave. In this last case, we also present comparisons with experimental data obtained in our shock tube lab.

FLOW MODELS

Several compressible multiphase flow models are included in ECOGEN. There are two general common features for these models which are:

-they are based on a diffuse interface method -they are hyperbolic (strictly or at least nearly).

The first property allows to treat a large variety of multiphase flows including mixtures in equilibrium state as well as interface flows (meaning that each phase is a pure phase from either side of an interface). A control volume of fluid is made of several phases (components). There is no interface reconstruction neither interface tracking. The interface is seen as a phase quantity, generally the volume fraction.

The compressibility property implies the use of equations of states (EOS). Various EOS are included in ECOGEN: perfect gas, stiffened gas. Next developments concern JWL, Mie-Grüneisen EOS and extension to tabular equations of state. Flows governed by Euler Equations are also included in ECOGEN (see paragraph 5.3).

The multiphase flows are basically described by the Kapila model [START_REF]Two-phase modeling of DDT in granular materials: Reduced equations[END_REF]. In this model, each phase evolves with the same pressure and with the same velocity, that is a mechanical equilibrium. Without any transfer between the phases this model is expressed as:

{ 𝜕𝛼 𝑘 𝜕𝑡 + 𝑢 ⃗ . ∇ ⃗ ⃗ 𝛼 𝑘 = 𝐾∇ ⃗ ⃗ . 𝑢 ⃗ 𝜕𝛼 𝑘 𝜌 𝑘 𝜕𝑡 + ∇ ⃗ ⃗ . (𝛼 𝑘 𝜌 𝑘 𝑢 ⃗ ) = 0 𝜕𝜌𝑢 ⃗ ⃗ 𝜕𝑡 + ∇ ⃗ ⃗ . (𝜌𝑢 ⃗ 𝑢 ⃗ + 𝑝𝐼 ̿ ) = 0 ⃗ 𝜕𝜌𝐸 𝜕𝑡 + 𝛻 ⃗ . [(𝜌𝐸 + 𝑝)𝑢 ⃗ ] = 0 (1)
In this system, the quantities with k index are related to k-th pure phase, with no index theses terms represent mixture quantities. 𝛼 𝑘 is the volume fraction, other quantities are written with usual notations. Notice that E is the total energy 𝐸 = 𝜀 + ‖𝑢 ⃗ ⃗ ‖ 2 2 .

In the volume fraction evolution equation, K represents the difference in the compressibility behaviors of each phase: it is linked to the acoustical impedance of pure phase. This term can be analytically expressed, or it can be computed by the means of relaxations terms in other models (8) (9) [START_REF] Zein | Modeling phase transition for compressible twophase flows applied to metastable liquids[END_REF].

ECOGEN also contains a mechanical and thermal equilibrium model. This last model is well designed for mixture of pure fluids problems, meaning that the approach is closed to Direct Numerical Simulation. The system of partial differential equations looks like the multi-components Euler equations:

𝜕𝛼 𝑘 𝜌 𝑘 𝜕𝑡 + ∇ ⃗ ⃗ . (𝛼 𝑘 𝜌 𝑘 𝑢 ⃗ ) = ∑ ṁ𝑘 𝑗 𝑗 𝜕𝜌𝑢 ⃗ ⃗ 𝜕𝑡 + ∇ ⃗ ⃗ . (𝜌𝑢 ⃗ 𝑢 ⃗ + 𝑝𝐼 ̿ ) = 0 ⃗ 𝜕𝜌𝐸 𝜕𝑡 + 𝛻 ⃗ . [(𝜌𝐸 + 𝑝)𝑢 ⃗ ] = 0 (2)
Even if the temperatures are the same for all the phases, this model allows the treatment of mass transfer between phase k and j (ṁ𝑘 𝑗 ) because the Gibbs free energies depend on the phase considered and are not in equilibrium state.

This large number of model in ECOGEN allows to treat a set of compressible flows. Some exchanges between phases can also be considered:

Mass transfer: the models included in ECOGEN enable the treatment of cavitating flows [START_REF]Diffuse interface models for high speed cavitating underwater systems[END_REF]. Combustion of a pure phase as well as detonation waves can be modelled [START_REF]Modelling detonation waves in condensed energetic materials: Multiphase CJ conditions and multidimensional computations[END_REF]. Change phase problems such as vaporization or condensation can also be solved.

Capillary effects: a new model proposed in (6) enables a conservative treatment of such problem. An extension of Model ( 1) is necessary which implies the treatment of the different physics by a numerical splitting. An example is provided in next section, devoted to the droplet fragmentation in a high-speed flow.

Thermal transfers: in ECOGEN are included heat transfers between the compressible phases as well as between a phase and a wall.

Viscosity effects: those effects may appear when momentum exchange occurs between phases and also as a dissipative term for pure phase flows.

CODE FEATURES

In this section are presented the main features of ECOGEN. These features involve numerical scheme for solving the flow equations as well as computing considerations.

The scheme is based on Godunov's method in a finite volume framework. Fluxes are solved thanks to exact or approximate Riemann solvers (HLLC). A MUSCL type extension provides high order accuracy in space and time.

The code is developed in C++. The parallelization is based on a geometrical domain decomposition and used MPI functions library to ensure the communications between the various domains.

The meshes used are structured or unstructured, can be done in 1-2-3D. ECOGEN can read gmesh mesh format for unstructured geometries but can also generate Cartesian grids.

In order to improve the quality of the solution around discontinuities an AMR technics has been developed. A new algorithm ensures efficient numerical solution for unsteady flows and the final code is much easier to develop than with previous existing methods [START_REF] Schmidmayer | Simplified Adaptive Mesh Refinement algorithm based on dual cell-boundary trees for multiphase compressible flows[END_REF]. AMR is only available for structured meshes.

The data visualization can be done with VTK tool, output files are provided under XML standard format.

RESULTS

Introduction

We present here a sample of results obtained with ECOGEN. The first example is the simulation of the atomization of a droplet in a high-speed flow including capillary effects. The second example is about the attenuation of a propagating shock wave in ducts: numerical results are compared with experimental ones.

Atomization of a droplet in a high-speed flow

A mathematical model for multiphase flows with capillary effects that ensures conservation of mass, energy and momentum is presented in reference [START_REF] Schmidmayer | A model and numerical method for compressible flows with capillary effects[END_REF]. This model also verifies the second law of thermodynamics for capillary flows. This new model is a wide extension of Kapila model, in which the supplemental physics is numerically solved thanks to an operator splitting involving three different stages.

The first one solves the equations of the multiphase flow with a pressure disequilibrium. In the second stage, the effects of surface tension are solved and finally, the pressure equilibrium between phases is recovered thanks to a relaxation procedure.

The surface tension is modelled via the Brackbill volume force (12) that required a color function 𝑐(OM ⃗⃗⃗⃗⃗⃗ , t) and this force is expressed via the capillary tensor to ensure a conservative form of the equations. It is noticeable that in the new approach, the color function 𝑐(OM ⃗⃗⃗⃗⃗⃗ , t) is no longer a phase quantity such as the volume or the mass fraction and thus the model is shown hyperbolic. To compute the capillary tensor the normal vector to the surface is required: it is computed via the gradient of 𝑐(OM ⃗⃗⃗⃗⃗⃗ , t). The components of the gradient of 𝑐(OM ⃗⃗⃗⃗⃗⃗ , t) obey a conservative equation which is solved.

A 3D numerical solution is now presented. The test case consists in the interaction of a water droplet and a gaseous high-speed flow (that may be a shock wave).

The liquid obeys a stiffened gas equation of state:

𝑝 = (𝛾 𝑙𝑖𝑞 -1)𝜌𝜀 -𝛾 𝑙𝑖𝑞 𝑃 ∞,𝑙𝑖𝑞 (4)
with 𝛾 𝑙𝑖𝑞 = 4.4 and 𝑃 ∞,𝑙𝑖𝑞 = 6.10 8 𝑃𝑎.

The air obeys a perfect gas law: 𝑝 = (𝛾 𝑎𝑖𝑟 -1)𝜌𝜀 with 𝛾 𝑎𝑖𝑟 = 1.4.

The test case is depicted in Figure 1 in which are given some flows data. An AMR mesh containing initially 250x50x50 cells for a physical domain of (250mm*50mm*50mm) is performed. The AMR method developed in ECOGEN is a new one (11) (13). 4 levels of refinement are used that leads to an equivalent non-AMR mesh size of 2.56*10 9 cells.

In Figure 2 is shown a view of the mixture density as well as the mesh. One can easily observe the mesh refinement around the droplet and at its rear (around the vortices). This view is plotted at time t=1ms. The initial flow velocity is around 152m/s: this velocity jump is due to an incoming shock wave with a pressure jump equal to 1.83. The high values of the density numerical schlieren obviously correspond to the liquid and appear in red in this plot. One can observe that the droplet is flattened with filaments at its side. These filaments will be torn away and form small droplet during the atomization process. The method can separate filaments from the 'mother droplet' as well as to form very small droplets.

These smaller droplets could be counted in Figure 3, that shows the liquid volume fraction at the same time. In Figure 4, are compared the CPU execution time if one considers 3 maximal levels of refinement or 4 levels versus the physical time. There is an order of magnitude but with only 3 levels of refinement, the simulation does not provide enough accurate results for this 3D test case, accuracy is considered here as the numerical diffusion at the liquid-gas interface.

Attenuation of a propagating shock wave

ECOGEN was used in support of an experimental study about the propagation and the attenuation of a shock wave through a 'Y' bifurcation ( 14) with a trap located in one of the branch (see Fig.

In view of the complexity of the present flow in the vicinity of the trap, we can reasonably consider the numerical calculation is in good agreement with the experimental observations. In Fig. 7 comparisons of pressure signals obtained experimentally and numerically for an incident shock wave Mach number of 1.44 are shown for a sensor located on the end-wall of the trapped branch. Fig. 7 shows that no significant difference (less than 5% on the recorded mean value) exists between numerical and experimental pressure traces. The small difference between the two signals is certainly due to the slightly error between the experimental and numerical sensor location. Note that the experimental value of the end-wall reflected pressure is slightly below the numerical one and can be explained by the non-ideal experimental conditions as calibration of sensors, parasitic losses in the device, small leaks or unperfect initial experimental conditions. The purpose of the numerical support in this experimental study is to complete the experimental investigation of the flow field in the device by adding some main physical mechanisms that only a numerical analysis could afford. As the velocity flow field behind the incident shock wave or the complete mapping of the pressure all along the device and not only at the sensor location. This analysis includes a numerical shadowgraph of the flow field behind the shock which highly depends on the incident Mach number. See Figure 8 for example.

Finally, a numerical parametric study on the size and the location of the trap in order to maximize the effects on the attenuation of this one on the end-wall reflected pressure. The numerical contribution in this case of investigation saves us time and a better design of the experimental setup. The time of calculation is about three hours for an area of 0.036m² made up of about 750,000 mesh cells with a labtop HP Zbook G3, 3 processors.

Figure 1 :

 1 Figure 1: Description of the 3D test case.

Figure 2 :

 2 Figure 2: Mesh and mixture density numerical schlieren at time t=1ms.

Figure 3 :

 3 Figure 3: Liquid volume fraction at time t=1ms.

Figure 4 :

 4 Figure 4: Computational time function of physical time for the 3D atomization tests.

  5.) to see what could be the impact on the attenuation of the end-wall reflected pressure. The incident shock wave Mach number was experimentally deduced from pressure records (PCB 113B26) taken at two different locations along the shock tube wall. In the present case the incident Mach number is 1.44 and moves from the left to the right.

Figure 5 :

 5 Figure 5: Picture of the shock tube with the 'Y' bifurcation with a trap device.The physical model solved with ECOGEN is based on a Euler equations solver. Numerical results are now compared to the experimental results both qualitatively (wave patterns in Fig.6) and quantitatively (pressure signals in Fig.7). Fig.6presents a comparison of schlieren pictures (experimental on the left versus numerical on the right) taken at same time. The schlieren variable calculated in the present work is the of the density gradient computed at each cell and visualized using open-source Paraview software. Specifically, the numerical schlieren is calculated here as follows.

Figure. 6 :

 6 Figure. 6: Experimental (left) and numerical (right) schlieren pictures showing the expansion of a planar shock wave through trapped duct for a Mach number equal to 1.44.

Figure. 7 :

 7 Figure. 7: Comparison of the evolution of the end-wall reflected pressure signal in the Y-shape for an incident shock wave Mach number of 1.44. The experimental (black line) and its numerical equivalent (red line).

Figure 8 :

 8 Figure 8: Numerical superposition of the schlieren density (bottom) and the pressure field (top) for Mach number of 1.44.


Yirmeyahu J Kaminski
email: jykaminski@gmail.com

Jean Lévine
email: jean.levine@mines-paristech.fr

François Ollivier
email: francois.ollivier@lix.polytechnique.fr

Intrinsic and Apparent Singularities in Differentially Flat Systems and Application to Global Motion Planning

Keywords: differential flatness, jets of infinite order, Lie-Bäcklund isomorphism, atlas, local chart, apparent and intrinsic singularity, global motion planning

Dans cet article, nous étudions les singularités des systèmes différentiellement plat dans la perpective de fournir des méthodes de planifications globales ou semi-globales : les sorties plates peuvent ne pas être définies globalement, ce qui peut empêcher la planification en dehors de leur domaine de définition, c'est-à-dire pour les singularités. Ces singularités sont de deux types : apparentes ou intrinsèques. Une définition rigoureuse est introduite en termes d'atlas et de cartes locales dans le cadre de la géométrie des jets d'ordre infini et des isomorphismes de Lie-Bäcklund. Nous donnons un résultat d'inclusion permettant le calcul effectif d'une partie de ces singularités intrinsèques. Enfin, nous montrons que nos résultats s'appliquent à la planification de trajectoire globale pour le fameux exemple de la voiture non holonome.

Introduction

Differential flatness has become a central concept in non-linear control theory for the past two decades. See [START_REF] Fliess | Flatness and defect of nonlinear systems : introductory theory and examples[END_REF][START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF], the overviews [START_REF] Ph | Flat systems, Plenary Lectures and Minicourses[END_REF][START_REF] Sira-Ramirez | Differentially Flat Systems[END_REF] and [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF] for a thoroughgoing presentation.

Consider a non-linear system on a smooth n-dimensional manifold X given by ẋ = f (x, u) (

where x ∈ X is the n-dimensional state vector and u ∈ R m the input or control vector, with m ≤ n to avoid trivial situations. We consider infinitely prolonged coordinates of the form (x, u) (x, u, u, ü, . . .)

∈ X × R m ∞ X × R m × R m × • • •
where the latter cartesian product is made of a countably infinite number of copies of R m . Roughly speaking, system (1) is said to be (differentially) flat 1 at a point (x 0 , u 0 ) (x 0 , u 0 , u0 , . . .) ∈ X × R m ∞ , if there exists an m-dimensional vector y = (y 1 , . . . , y m ) satisfying the following statements :

y is a smooth function of x, u and time derivatives of u up a to a finite order β = (β 1 , . . . , β m ), i.e. y = Ψ(x, u, u, . . . , u (β) ), where u (β) stands for (u

(β 1 ) 1 , . . . , u (β m )
m ) and where u

(β i ) i
is the β i th order time derivative of u i , i = 1, . . . , m, in a neighborhood of the point (x 0 , u 0 ) ; y and its successive time derivatives ẏ, ÿ, . . . are locally differentially independent in this neighborhood ; x and u are smooth functions of y and its time derivatives up to a finite order α = (α 1 , . . . , α m ), i.e. (x, u) = Φ(y, ẏ, . . . , y (α) ) in a neighborhood of the point (y 0 , ẏ0 , . . .) (Ψ(x 0 , u 0 , u0 , . . . , u (β) 0 ), Ψ(x 0 , u 0 , u0 , . . . , u (β+1) 0

), . . .).

Then the vector y is called flat output.

Note that the latter informal definition becomes rigorous if we regard the above defined functions Φ and Ψ as smooth functions over infinite order jet spaces endowed with the product topology 2 [START_REF] Krasil'shchik | Geometry of Jet Spaces and Nonlinear Partial Differential Equations[END_REF][START_REF] Zharinov | Geometric Aspects of Partial Differential Equations[END_REF][START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF][START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF]. They are then called Lie-Bäcklund isomorphisms and are inverse one of each other (see [START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF][START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF]). However, these functions may be defined on suitable neighborhoods that need not cover the whole space. We thus may want to know where such isomorphisms do not exist at all, a set that may be roughly qualified of intrinsically singular, thus motivating the present work : if two points are separated by such an intrinsic singularity, it is intuitively impossible to join them by a smooth curve satisfying the system differential equations and, thus, to globally solve the motion planning problem 3 . More precisely, the notions of apparent and intrinsic singularities are introduced thanks to the construction of an atlas, that we call Lie-Bäcklund atlas, where local charts are made of the open sets where the Lie-Bäcklund isomorphisms, defining the flat outputs, are non degenerated, in the spirit of [START_REF] Chang | Construction of an atlas for global flatness-based parameterization and dynamic feedback linearization of quadcopter dynamics[END_REF][START_REF] Chang | Global Chartwise Feedback Linearization of the Quadcopter with a Thrust Positivity Preserving Dynamic Extension[END_REF] where a comparable idea was applied to a quadcopter model. Intrinsic singularities are then defined as points where flat outputs fail to exist, i.e. that are contained in no above defined chart at all. Other types of singularities are called apparent, as they can be ruled out by switching to another flat output well defined in an intersecting chart. Our intrinsic singularity notion may be seen as a generalization of the one introduced in [START_REF] Li | Flat outputs of two-inputs driftless control systems[END_REF] in the particular case of two-input driftless systems such as cars with trailers, and restricted to the so-called x-flat outputs.

Our main result, apart from the above Lie-Bäcklund atlas and singularities definition, then concerns the inclusion of a remarkable and effectively computable set in the set of intrinsic singularities. Note that, since finitely computable necessary and sufficient conditions of non existence of flat output are not available in general [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF], an easily computable complete characterization of the set of intrinsic singularities is not presently known and it may be useful to label all or part of the singularities as intrinsic thanks to their membership of another set.

To briefly describe this result, we start from the necessary and sufficient conditions for the existence of local flat outputs of meromorphic systems of [START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF] 4 . It consists in firstly transforming the system (1) in the locally equivalent implicit form :

F(x, ẋ) = 0 ( 2
)
where F is assumed meromorphic, and introducing the operator τ, the trivial Cartan field on the manifold of global coordinates (x, ẋ, ẍ, . . .), given by

τ = ∑ n i=1 ∑ j≥0 x (j+1) i ∂ ∂x (j) i
. Then, we compute the diagonal or Smith-Jacobson decomposition [START_REF] Cohn | Free Rings and Their Relations[END_REF][START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF] of the following polynomial matrix :

P(F) = ∂F ∂x + ∂F ∂ ẋ τ (3)
a matrix that describes the variational system associated to [START_REF] Antritter | An efficient algorithm for checking hyperregularity of matrices[END_REF], and that lies in the ring of matrices whose entries are polynomials in the operator τ with meromorphic coefficients. We prove that the set of intrinsic singularities contains the set where P(F) is not hyper-regular (see [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF]). As a corollary, we deduce that if an equilibrium point is not first order controllable, then it is an intrinsic singularity.

These results are applied to the global motion planning problem of the well-known non-holonomic car, which is only used here as a benchmark in order to show how the classical and simple flatness-based motion planning space associated to the system, once the set of intrinsic singularities has been removed, can be joined by a system's trajectory, and thus that this set is connected by arcs.

4. Other approaches to flatness characterization may be found in [START_REF] Aranda-Bricaire | A linear algebraic framework for dynamic feedback linearization[END_REF][START_REF] Chetverikov | New flatness conditions for control systems[END_REF][START_REF] Antritter | Flatness characterization : two approaches[END_REF] methodology can be extended in presence of singularities. It is also meant to help the reader verifying that the introduced concepts, in the relatively arduous context of Lie-Bäcklund isomorphisms, are nevertheless intuitive and well suited to this situation. Note that different approaches, also leading to global results, have already been extensively developed in the context of non holonomic systems, based on controllability, Lie brackets of vector fields and piecewise trajectory generation by sinusoids [START_REF] Murray | Nonholonomic motion planning : steering using sinusoids[END_REF][START_REF] Jean | The car with n trailers. characterisation of the singular configurations[END_REF][START_REF] Chitour | A global steering method for nonholonomic systems[END_REF][START_REF] Jean | Control of Nonholonomic Systems : from Sub-Riemannian Geometry to Motion Planning[END_REF], or using Brockett-Coron stabilization results [START_REF] Brockett | Asymptotic stability and feedback stabilization[END_REF][START_REF] Coron | Global asymptotic stabilization for controllable systems without drift[END_REF]. However, though some particular nonholonomic systems, as the car example, happen to be flat, our approach applies to the class of flat systems which is different, including e.g. pendulum systems, unmanned aerial vehicles and many others that do not belong to the nonholonomic class (see [START_REF] Ph | Flat systems, Plenary Lectures and Minicourses[END_REF][START_REF] Sira-Ramirez | Differentially Flat Systems[END_REF][START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF][START_REF] Chang | Construction of an atlas for global flatness-based parameterization and dynamic feedback linearization of quadcopter dynamics[END_REF][START_REF] Chang | Global Chartwise Feedback Linearization of the Quadcopter with a Thrust Positivity Preserving Dynamic Extension[END_REF]).

Remark that, in the car example, the obtained intrinsic singularities are the same as the ones revealed in [START_REF] Murray | Nonholonomic motion planning : steering using sinusoids[END_REF][START_REF] Jean | The car with n trailers. characterisation of the singular configurations[END_REF][START_REF] Chitour | A global steering method for nonholonomic systems[END_REF][START_REF] Jean | Control of Nonholonomic Systems : from Sub-Riemannian Geometry to Motion Planning[END_REF] where first order controllability fails to hold, or, according to [START_REF] Brockett | Asymptotic stability and feedback stabilization[END_REF][START_REF] Coron | Global asymptotic stabilization for controllable systems without drift[END_REF], where stabilisation by continuous state feedback is impossible. However, the degree of generality of this coincidence is not presently known.

The paper is organized as follows. In section 2, we introduce the basic language of Lie-Bäcklund atlas and charts. Then this leads to a computational approach for calculating intrinsic singularities. In particular, their links with the hyper-singularity of the polynomial matrix (3) of the variational system are established in Proposition 2 and Theorem 1, and then specialized in Corollary 1 to the case of equilibrium points.

In section 3, we apply our results to the non holonomic car. We build an explicit Lie-Bäcklund atlas for this model, compute the set of intrinsic singularities and apply the atlas construction to trajectory planning where the route contains several apparent singularities and starts and ends at intrinsically singular points. Finally, conclusions are drawn in section 4.

Lie-Bäcklund Atlas, Apparent and Intrinsic Singularities

Recall from the introduction that we consider the controlled dynamical system in explicit form [START_REF] Antritter | On the computation of π-flat outputs for linear time-varying differential-delay systems[END_REF], where x evolves in some n-dimensional manifold X. The control input u lies in R m . Then the system can be seen as the zero set of ẋf (x, u) in TX × R m , where TX is the tangent bundle of X. From now on, we assume that the Jacobian matrix ∂ f ∂u (x, u) has rank m for every (x, u). Converting system (1) into its implicit form consists in eliminating the input u or, more precisely, in computing its image by the projection π from TX × R m onto TX to get the implicit relation [START_REF] Antritter | An efficient algorithm for checking hyperregularity of matrices[END_REF], where we assume that F :

(x, ẋ) ∈ TX → R n-m is a meromorphic function, with m ≤ n.
Following [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF], we embed the state space associated to (2) into a diffiety (see [START_REF] Zharinov | Geometric Aspects of Partial Differential Equations[END_REF]), i.e. into the manifold X X × R n ∞ , where we have denoted by R n ∞ the product of a countably infinite number of copies of R n , with coordinates

x (x, ẋ, ẍ, . . . , x (k) , . . .), endowed with the trivial Cartan field :

τ X n ∑ i=1 ∑ j≥0 x (j+1) i ∂ ∂x (i) i .
Note that τ X is such that the elementary relations τ X x (k) = x (k+1) hold for all k ∈ N. The integral curves of both ( 1) and ( 2) thus belong to the zero set of {F,

τ k X F | k ∈ N} in X.
However, there might exist points x = (x, ẋ, ẍ, . . . , x (k) , . . .) ∈ X such that the fiber π -1 (x, ẋ) above x is empty, i.e. such that there does not exist a u ∈ R m such that ẋf (x, u) = 0. We indeed naturally exclude such points. It is easily proven that the integral curves of ( 1) and (2) coincide on the set X 0 given by

X 0 = {x ∈ X | τ k X F(x) = 0, ∀k ∈ N} \ {x ∈ X | π -1 (x, ẋ) = ∅}.
Therefore, the system trajectories are uniquely defined by the triple (X, τ X , F) that we call the system from now on (see [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF]). Without loss of generality, we may consider that this system is restricted to X 0 .

In order to get rid of any reference to an explicit system, such as the complementary of the empty fibers of the projection π, we more generally assume that X 0 is an open dense subset 5 of {x ∈ X | τ k X F(x) = 0, ∀k ∈ N}. Let us recall the definitions of Lie-Bäcklund equivalence and local flatness for implicit systems ( [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF]) :

Consider two systems (X, τ X , F) and (Y,

τ Y , G) where Y Y × R q ∞
, Y being a q-dimensional smooth manifold, where q ∈ N is arbitrary, with global coordinates y (y, ẏ, . . .) and trivial Cartan field

τ Y ∑ q i=1 ∑ j≥0 y (j+1) i ∂ ∂y (j) i .
As before, we denote by

Y 0 an open dense subset of {y ∈ Y | τ k Y G(y) = 0, ∀k ∈ N}.
Definition 1. We say that (X, τ X , F) and (Y, τ Y , G) are Lie-Bäcklund equivalent at a pair of points (x 0 , y 0 ) ∈ X 0 × Y 0 if, and only if, (i) there exist neighborhoods X 0 of x 0 in X 0 , and Y 0 of y 0 in Y 0 , and a one-to-one mapping Φ = (ϕ 0 , ϕ 1 , . . .), meromorphic from Y 0 to X 0 , satisfying Φ(y 0 ) = x 0 and such that the restrictions of the trivial Cartan fields

τ Y Y 0 and τ X X 0 are Φ-related, namely Φ * τ Y Y 0 = τ X X 0 ;
(ii) there exists a one-to-one mapping Ψ = (ψ 0 , ψ 1 , . . .), meromorphic from X 0 to Y 0 , such that Ψ(x 0 ) = y 0 and

Ψ * τ X X 0 = τ Y Y 0 .
The mappings Φ and Ψ are called mutually inverse Lie-Bäcklund isomorphisms at (x 0 , y 0 ).

5.

As a consequence of the implicit function theorem, the set of points where the fibers are empty is the complement of an open dense subset of the set {x ∈

X | τ k X F(x) = 0, ∀k ∈ N}.
The two systems (X, τ X , F) and (Y, τ Y , G) are called locally L-B equivalent if they are L-B equivalent at every pair (x, Ψ(x)) = (Φ(y), y) of an open dense subset Z of X 0 × Y 0 , with Φ and Ψ mutually inverse Lie-Bäcklund isomorphisms on Z.

Accordingly, Definition 2. The system (X, τ X , F) is said (differentially) flat at x 0 if, and only if, it is Lie-Bäcklund equivalent to the trivial system (R m ∞ , τ, 0) at (x 0 , y 0 ) where τ is the trivial Cartan field on R m ∞ with global coordinates6 y = (y, ẏ, . . .

), i.e. τ = ∑ m i=1 ∑ j≥0 y (j+1) i ∂ ∂y (j) i
, and where 0 indicates that there is no differential equation to satisfy. In this case, we say that y, or Ψ by extension, is a local flat output, welldefined and invertible from a neighborhood of x 0 to a neighborhood of y 0 . Finally, the system (X,

τ X , F) is said locally (differentially) flat if it is flat at every point of an open dense subset Z of X 0 × R m ∞ .

Lie-Bäcklund Atlas

From now on, we assume that system (1), or equivalently [START_REF] Antritter | An efficient algorithm for checking hyperregularity of matrices[END_REF] or, also equivalently, system (X, τ X , F) is locally flat.

We now introduce the notion of a Lie-Bäcklund atlas for flat systems. It consists of a collection of charts on X 0 , that we call Lie-Bäcklund charts and atlas, and that will allow us to define a structure of infinite dimensional manifold on a subset of X 0 , that can be X 0 itself is some cases.

Definition 3.

(i) A Lie-Bäcklund chart on X 0 is the data of a pair (U, ψ) where

U is an open set of X 0 and ψ : U → R m ∞ a local flat output, with local inverse ϕ : V → U with V open subset of ψ(U) ⊂ R m ∞ . (ii) Two charts (U 1 , ψ 1
) and (U 2 , ψ 2 ) are said to be compatible if, and only if, the mapping

ψ 1 • ϕ 2 : ψ 2 (ϕ 1 (V 1 ) ∩ ϕ 2 (V 2 )) ⊂ R m ∞ → ψ 1 (ϕ 1 (V 1 ) ∩ ϕ 2 (V 2 )) ⊂ R m ∞
is a local Lie-Bäcklund isomorphism (with the same trivial Cartan field τ associated to both the source and the target) with local inverse ψ 2 • ϕ 1 , as long as

ϕ 1 (V 1 ) ∩ ϕ 2 (V 2 ) = ∅. (iii) An atlas A is a collection of compatible charts.

For a given atlas

A = (U i , ψ i ) i∈I , let U A be the union U A i∈I U i .
Here our definition differs from the usual concept of atlas in finite dimensional differential geometry, since, on the one hand, diffeomorphisms are replaced by Lie-Bäcklund isomorphisms and, on the other hand, we do not require that U A = X 0 . The reason for this difference is precisely related to our objective, i.e. identifying the essential singularities of differentially flat systems. This will become clear in the sequel.

Apparent and Intrinsic Flatness Singularities

It is clear from what precedes that if we are given two Lie-Bäcklund atlases, their union is again a Lie-Bäcklund atlas. Therefore the union of all charts that form every atlas is well-defined as well as its complement, which we call the set of intrinsic flatness singularities, as stated in the next definition. Definition 4. We say that a point in X 0 is an intrinsic flatness singularity if it is excluded from all charts of every Lie-Bäcklund atlas. Every other singular point, namely every point x ∈ U i for some chart (U i , ψ i ) but for which there exists another chart (U j , ψ j ), j = i, such that x ∈ U j , is called apparent.

Clearly, this notion does not depend on the choice of atlas and charts. The concrete meaning of this notion is that at points that are intrinsic singularities there is no flat output, i.e. the system is not flat at these points.

On the other hand, points that are apparent singularities are singular for a given set of flat outputs, but well defined points for another set of flat outputs.

Note, moreover, that obtaining atlases may be very difficult in general situations and a computable criterion to directly detect intrinsic singularities should be of great help. A simple result in this direction is presented in the following section 2.3.

Intrinsic Flatness Singularities and Hyper-regularity

The purpose of this section is to give a tractable sufficient condition of intrinsic singularity and an algorithm to effectively compute the associated points.

With the notations defined at the beginning of section 2, we next consider the variational equation, in polynomial form, of system (2) :

P(F)dx = 0, P(F) = ∂F ∂x + ∂F ∂ ẋ τ X (4)
where the entries of the (nm) × n matrix P(F) are polynomials in τ X with meromorphic functions on X as coefficients.

Recall that a square n × n polynomial matrix is said to be unimodular if it is invertible and if its inverse is also a matrix whose entries are polynomials in τ X with meromorphic functions on X as coefficients. It is of importance to remark that, according to the fact that the coefficients are meromorphic functions, they are, in general, only locally defined. This local dependence will be omitted unless explicitly needed.

The (nm) × n polynomial matrix P(F) is said hyper-regular if, and only if, there exists an (nm) × (nm) unimodular polynomial matrix V and an n × n unimodular polynomial matrix U such that

VP(F)U = I n-m 0 (n-m)×m . (5)
In fact, it has been proven in [START_REF] Antritter | An efficient algorithm for checking hyperregularity of matrices[END_REF] (see also [START_REF] Antritter | On the computation of π-flat outputs for linear time-varying differential-delay systems[END_REF]Proposition 1]), that the latter definition may be simplified as follows :

Proposition 1. The polynomial matrix P(F) is hyper-regular if, and only if, there exists an n × n unimodular polynomial matrix U such that

P(F)U = I n-m 0 (n-m)×m . ( 6
)
Démonstration. P(F) is hyper-regular if, and only if, there are matrices S, of size (nm) × (nm) and T of size n × n such that SP(F)T = I n-m 0 (n-m)×m . Thus, using the identity

I n-m 0 (n-m)×m = S -1 I n-m 0 (n-m)×m S 0 (n-m)×m 0 m×(n-m) I m
we get

I n-m 0 (n-m)×m = S -1 (SP(F)T) S 0 0 I m = P(F) T S 0 0 I m P(F)U
which proves [START_REF] Chang | Construction of an atlas for global flatness-based parameterization and dynamic feedback linearization of quadcopter dynamics[END_REF]. The converse is trivial

We say that P(F) is hyper-singular at a given point if, and only if, it is not hyper-regular at this point, i.e. if this point does not belong to any neighborhood where P(F) is hyper-regular or, in other words, if at this point no unimodular matrix U satisfying (6) exists.

Let us denote by S F the subset of X 0 where P(F) is hyper-singular. The following proposition clarifies some previous results of [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF] in the context of flat systems at a point : Proposition 2. If system (2) is flat at the point x 0 ∈ X 0 , then there exists a neighborhood V of x 0 where P(F) is hyper-regular.

Démonstration. Assume that system (2) is flat at the point x 0 ∈ X 0 . Then, denoting as before y (y, ẏ, ÿ, . . .) and x (x, ẋ, ẍ, . . .), by definition, there exists a neighborhood V of x 0 and a flat output y = Ψ(x) (Ψ 0 (x),

Ψ 1 (x), Ψ 2 (x), . . .) ∈ Ψ(V) ⊂ R m
∞ for all x ∈ V and conversely, x = Φ(y) (Φ 0 (y), Φ 1 (y), Φ 2 (y), . . .) for all y ∈ Ψ(V) such that F(Φ 0 (y), Φ 1 (y)) = F(Φ 0 (y), τΦ 0 (y)) ≡ 0.

Taking differentials, we show that dy is a flat output of the variational system. Considering the Jacobian matrix dΦ 0 (y) (resp. dΨ 0 (x)) of the 0th order component Φ 0 (resp. Ψ 0 ) of Φ (resp. Ψ), we denote by P(Φ 0 ) (resp. P(Ψ 0 )) its polynomial matrix form with respect to τ (resp. w.r.t. τ X ) (see [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF]).

Since dy = dΨ(x)dx and dx = dΦ(y)dy, we get that dx = P(Φ 0 )dy ∈ T * V, dy = P(Ψ 0 )dx ∈ T * Ψ(V), P(F)P(Φ 0 ) ≡ 0 and P(Φ 0 ) left-invertible, since P(Ψ 0 )P(Φ 0 ) = I m .

We next consider the Smith-Jacobson decomposition, or diagonal decomposition [START_REF] Cohn | Free Rings and Their Relations[END_REF]Chap. 8], of P(F) : there exists an (nm) × (nm) unimodular matrix W, an n × n unimodular matrix U and an (nm) × (nm) diagonal matrix ∆ such that WP(F)U = ∆ 0 . Partitionning U into U 1 U 2 , we indeed get WP(F)U 1 = ∆ and WP(F)U 2 = 0, or P(F)U 2 = 0 since W is unimodular. Thus, by elementary matrix algebra, taking account of the independence of the columns of both U 2 and P(Φ 0 ), one can choose U such that U 2 = P(Φ 0 ).

Following [START_REF] Fliess | A remark on Willems' trajectory characterization of linear controllability[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF] (see also [START_REF] Antritter | On the computation of π-flat outputs for linear time-varying differential-delay systems[END_REF] in a more general context), we introduce the free differential module K[dy] finitely generated by dy 1 , . . . , dy m over the ring K of meromorphic functions from X 0 to R and the differential quotient module H K[dx]/K[P(F)dx] where K[P(F)dx] is the differential module generated by the rows of P(F)dx. Taking an arbitrary non zero element z = (z 1 , . . . , z m ) in K[dy], and its image ξ = P(Φ 0 )z, we immediately get P(F)ξ = P(F)P(Φ 0 )z = 0 which proves that ξ is equivalent to zero in H. Since

U = U 1 P(Φ 0 ) is unimodular, it admits an inverse V = V 1 V 2
and thus

U 1 V 1 + P(Φ 0 )V 2 = I n .
Multiplying on the left by WP(F) and on the right by ξ, and using the relation

P(F)P(Φ 0 ) = 0, we get 0 = WP(F)ξ = WP(F)U 1 V 1 ξ + WP(F)P(Φ 0 )V 2 ξ = WP(F)U 1 V 1 ξ. Consequently, recalling that WP(F)U 1 = ∆, we have that ζ V 1 ξ = V 1 P(Φ 0 )z satisfies 0 = WP(F)U 1 ζ = ∆ζ.
Consequently, if the entries of the diagonal matrix ∆ contain at least one polynomial of degree larger than 0 with respect to τ, say δ i for some i = 1, . . . , nm, then ). Therefore, the entries of the matrix ∆ must belong to K, which implies that there exists a submatrix U ′ 1 such that U ′ U ′ 1 P(Φ 0 ) is unimodular and satisfies WP(F)U ′ = I n-m 0 , and thus, according to [START_REF] Antritter | An efficient algorithm for checking hyperregularity of matrices[END_REF] or [1, Proposition 1], that P(F) must be hyper-regular in the considered neighborhood.

δ i ζ i = 0,

Remark 1. The above proof may be summarized by the following diagram of exact sequences :

0 -→ R m ∞ Φ -→ ←- Ψ X 0 F -→ 0 d ↓ ↓ d 0 -→ TR m ∞ dΦ -→ ←- dΨ TX 0 P(F) -→ 0 Since TR m ∞ , is isomorphic to the free differential module K[dy],
then TX 0 , that may also be seen as a differential module, is necessarily free. In other words, the kernel of P(F) must be equal to the image of TR m ∞ by the one-to-one linear map dΦ, thus sending a basis of TR m ∞ (flat outputs) to a basis of TX 0 .

Remark 2. Due to the Smith-Jacobson decomposition, the hyper-regularity property gives a practical row-reduction algorithm to compute S F (see [START_REF] Antritter | An efficient algorithm for checking hyperregularity of matrices[END_REF] and the car example in section 3.3). The hyper-singular set is then deduced by complementarity.

According to Proposition 2, it is clear that on S F , the system cannot be flat. We thus have the following straightforward result : Theorem 1. The set S F is contained in the set of flatness intrinsic singularities of the system.

In fact (see [START_REF] Fliess | A remark on Willems' trajectory characterization of linear controllability[END_REF][START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF]), S F corresponds to the points where the system is no more F-controllable, i.e. controllable in the sense of free modules, and therefore non flat (see [START_REF] Charlet | Sufficient conditions for dynamic state feedback linearization[END_REF][START_REF] Fliess | Flatness and defect of nonlinear systems : introductory theory and examples[END_REF][START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF][START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF]). As a consequence of this theorem, the points where the matrix P(F) is hyper-singular are automatically intrinisic singularities of the system.

Note that, at equilibrium points, F-controllability boils down to first order controllability, i.e. controllability of the tangent linear system. Corollary 1. The set made of equilibrium points that are not first order controllable is contained in the set of flatness intrinsic singularities of the system.

Applications : Route Planning For the Non Holonomic Car

In this section, we show on a specific example how the above carried out theoretical analysis applies.

Car Model

The car (kinematic) model is made of the following set of explicit differential equations (see e.g. [START_REF] Murray | Nonholonomic motion planning : steering using sinusoids[END_REF]) :

   ẋ = u cos θ ẏ = u sin θ θ = u l tan ϕ (7)
Details about the notations are given in the caption of figure 1. In explicit form, the system evolves in the manifold 1 where the variables are (x, y, θ, u, ϕ). For the sake of clarity, we note X 11 = R 2 × S 1 for the space of state variables (x, y, θ) and X 12 = R × S 1 for the space of control variables (u, ϕ). The tangent bundle of X 11 is denoted by TX 11 . This system can thus be seen as the zero set in TX 11 × X 12 of the following function :

X 1 = R 2 × S 1 × R × S
F(x, y, θ, ẋ, ẏ, θ, u, ϕ) =   ẋ -u cos θ ẏ -u sin θ θ -u l tan ϕ  
As in section 2 and again following [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF][START_REF] Lévine | On necessary and sufficient conditions for differential flatness[END_REF], we consider the local implicit representation of the system, obtained by projecting F on TX 11 by the canonical projection π : TX 11 × X 12 → TX 11 , which amounts to eliminating the controls. In this context, the dynamics [START_REF] Chang | Global Chartwise Feedback Linearization of the Quadcopter with a Thrust Positivity Preserving Dynamic Extension[END_REF] are locally equivalent to the zero set of the following function :

F(x, y, θ, ẋ, ẏ, θ) = ẋ sin θ -ẏ cos θ = 0. ( 8
)
Figure 1 -Car Model : the state vector is made of the coordinates (x, y) of the rear axle's center and of the angle θ between the car's axis and the x-axis. The controls are the speed u and the angle ϕ between the wheels' axis and the car's axis. The length l is the distance between the two axles.

We then embed the state space associated to (8) into the diffiety X = R 2 × S 1 × R 3 ∞ , endowed with the trivial Cartan field :

τ X = 3 ∑ i=1 ∑ j≥0 x (j+1) i ∂ ∂x (i) i
, where

we have set x 1 = x, x 2 = y and x 3 = θ.

The system trajectories now live in X 0 , the subset of {x ∈ X | τ k X F = 0, ∀k ∈ N}, where we have excluded the set Z {(x, y, θ, ẋ, ẏ, θ) ∈ TX 11 } | ẋ = ẏ = 0, θ = 0} of points of TX 11 where the fibers associated to π are empty, i.e. the points of TX 11 such that there does not exist u and ϕ such that F(x, ẋ) = 0 (see section 2). Thus

X 0 {x ∈ X | τ k X F = 0, ∀k ∈ N} \ Z.

Lie-Bäcklund Atlas for the Car Model

We now define an atlas on X 0 by simply enumerating the charts, as in [START_REF] Chang | Construction of an atlas for global flatness-based parameterization and dynamic feedback linearization of quadcopter dynamics[END_REF][START_REF] Chang | Global Chartwise Feedback Linearization of the Quadcopter with a Thrust Positivity Preserving Dynamic Extension[END_REF] in the context of quadcopters. Each chart is defined on an open set associated to a local Lie-Bäcklund isomorphism ψ i from X 0 to R 2 ∞ with local inverse denoted by φ i : R 2 ∞ → X 0 . For simplicity's sake, we only define φ i by its three first components. The other ones are deduced by differentiation, i.e. by applying τ X to them an arbitrary number of times. A similar abuse of notation has been used for the definition of ψ i . A point in X 0 is denoted by x.

Over U 1

{ ẋ = 0}, we take y 1 = (x, y) = ψ 1 (x) and the inverse Lie-Bäcklund transform is given by

φ 1 =   x y tan -1 ( ẏ ẋ )   2.
Over U 2 { ẏ = 0}, we take y 2 = (x, y) = ψ 2 (x) and the inverse Lie-Bäcklund transform is given by

φ 2 =   x y cotan -1 ( ẋ ẏ )   3. Over U 3 { θ = 0}, we take y 3 = (θ, x sin θ -y cos θ) = ψ 3 (x).
Here for the sake of simplicity, we shall denote (z 1 , z 2 ) the components of y 3 . In that case the inverse Lie-Bäcklund transform is given by

φ 3 =    ż2 ż1 cos z 1 + z 2 sin z 1 ż2 ż1 sin z 1 -z 2 cos z 1 z 1    4.
Finally note that the above charts do not contain the set V = X 0 \ 3 i=1 U i = { ẋ = ẏ = θ = 0}, which corresponds to the set of equilibrium points of the system. Note that, by the definition of X 0 , ẋ = ẏ = 0 implies θ = 0. Therefore, V = X 0 \

3 i=1 U i = { ẋ = ẏ = 0}
One can check that for all i, j, Im(φ i ) ⊂ X 0 and that the ψ j • φ i 's satisfy the compatibility definition of section 2.1 on R 2 ∞ . Therefore we have indeed defined an atlas of 3 i=1 U i = X 0 \ { ẋ = ẏ = 0}. Among other things, this allows us to conclude that the car dynamics is globally controllable provided one avoids the singular set V, as illustrated in section 3. Note that at this level, we are not able to conclude that the set { ẋ = ẏ = 0} is an intrinsic flatness singularity since, according to definition 4 above, we still have to prove that no other atlas can contain this set, hence the importance of the next section based on the results of section 2.3.

Flat Outputs and Intrinsic Flatness Singularities of the Car Example

One first considers the differential of the implicit equation : Note that this direct computation, from the variational system, of the intrinsic singularity confirms that the atlas construction of section 2 was complete in the sense that adding more charts would not reduce the set of intrinsic singularities.

Remark 3. Let us stress that the intrinsic singularity obtained in section 3.2 and the planned trajectory of the next section 3 do not depend on the choice of atlas and charts. Another choice, using e.g. the formulas given in [START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF]Section 6.2.4] would be equally possible, leading to a similar construction. Remark 4. In this example, we could prove that S F is in fact equal to the set of intrinsic singularities of the system. Indeed, it would be most interesting to have an idea of the generality of this situation. However, examples where S F does not coincide with the set of flatness intrinsic singularities of the system are not presently known by the authors.

Route Planning

Next, we show how the previously built atlas can be used to control the car over a route along which there are several apparent and intrinsic singularities, as the one depicted in figure 2. This route has been defined in several steps. First, the way points A, C and following, up to K, were chosen in the (x, y)-plane to start from the equilibrium point A (intrinsic singularity) along the y-axis, which is an apparent singularity for y 1 (see section 3.2). The car accelerates up to B and then travels at constant speed up to C where it starts making a right turn up to D. The route between C and D has been designed by a univariate spline fitting in order to join the previous vertical line to the horizontal segment DE, an apparent singularity for y 2 . The next segment FG, after the arc EF, again designed by spline fitting, corresponds to a constant heading angle θ, an apparent singularity for y 3 . Finally, on the arc H J, the car speed remains constant and then The flat outputs x and y parametrized by time thus an intrinsic singularity.

The whole route has been parametrized, in a first step, by its arc length variable on the interval [0, L], with unit speed, in order to allow the design of an arbitrary speed profile over time (see Fig. 4).

The trajectory design is done according to the flatness-based method described in [START_REF] Rouchon | Flatness and motion planning : the car with n-trailers[END_REF][START_REF] Rouchon | Flatness, motion planning and trailer systems[END_REF] on each route section. The flat output used is y 2 on AC, y 1 on CE, indifferently y 1 or y 2 on EG, and y 1 on GK since the component y attains its minimum on this arc, thus with ẏ = 0.

The obtained speed profile of the car is shown in figure 3.

For the computation of ϕ, we exclude the end points where the speed vanishes and thus where ϕ is only asymptotically defined. See figure 5. Those points, which are indeed intrinsic singularities, can be approached as close as we want but exactly stopping on them with a prescribed orientation and bounded controls is impossible.

Concluding Remarks

In this paper, the concepts of intrinsic and apparent flatness singularities have been defined. These notions are of paramount importance for global trajectory planning, namely planning through apparent singularities, avoiding intrinsic singularities, with the possibility of approaching them as close as possible.

We have also shown that intrinsic singularities include a remarkable set, namely the points where the matrix P(F) of the variational system, which plays a major role in the process of flat output computation, is hyper-singular.

This analysis is illustrated by the global motion planning of a non holonomic car. In this context, we have exhibited an atlas of flat outputs and a complex trajectory safely passing through all possible charts of this atlas.

Note that this approach may be applied in the same way to other flat systems which do not belong to the class of nonholonomic systems. Moreover, it might be possible to extend it to the computation of the largest reachable set of a system.

  dF = d ẋ sin θ + ẋ cos θdθd ẏ cos θ + ẏ sin θdθ = ( ẋ cos θ + ẏ sin θ)dθ + sin θd ẋcos θd ẏ

  the plane (x,y) parametrized by arc length

Figure 2 -

 2 Figure 2 -Planned car route, parametrized by arc length.

Figure 3 -

 3 Figure3-The speed corresponding to the route depicted in figure2

Figure 4 -

 4 Figure 4 -The flat outputs parametrized first by arc length and then by time corresponding to the route depicted in figure 2

Figure 5 -

 5 Figure5-The angles θ and ϕ parametrized by time corresponding to the route depicted in figure2. For the computation of ϕ, the car length has been chosen equal to l = 2m.

  and since ζ i ∈ K[dy], we have proven that the non zero component ζ i is a torsion element of K[dy], thus leading to a contradiction with the fact that K[dy] is free (see e.g. [19, Theorem 7.3, Chap. III] or [11, Corollary 2.2, Chap. 8, Sec. 8.2]

The number of components of y must be equal to m (see[START_REF] Fliess | A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems[END_REF][START_REF] Lévine | Analysis and Control of Nonlinear Systems : A Flatness-based Approach[END_REF]).

Note that, if z is an arbitrary variable of the system, we have d ż = d(τ X z) = τ X dz, i.e. the exterior derivative d commutes with the Cartan field τ X , and the matrix P(F) reads as follows :

for all dx, dy, dθ that are differentials of the variables x, y, θ satisfying system [START_REF] Charlet | Sufficient conditions for dynamic state feedback linearization[END_REF]. Now in the context of the car system given by ( 8), we are ready to prove the following : Proposition 3. The intrinsic singular set of system [START_REF] Charlet | Sufficient conditions for dynamic state feedback linearization[END_REF], given by { ẋ = ẏ = 0}, is equal to S F . Démonstration. We compute the set where P(F) is not hyper-regular. Let us define A = ẋ cos θ + ẏ sin θ.

Up to a column permutation, P(F) reads [A, (sin θ)τ X , -(cos θ)τ X ]. Then the first column of U, say u 1 is u 1 = [1/A, 0, 0] t (the superscript t denotes the transposition operator). The second one u 2 is given by [P 0 , P 1 , P 2 ] t where P 0 , P 1 , P 2 are polynomials of τ X with deg (P 0 ) = 1 + max i=1,2 deg (P i ), such that AP 0 + (sin θ)τ X P 1 -(cos θ)τ X P 2 = 0, or P 0 = -1 A ((sin θ)τ X P 1 -(cos θ)τ X P 2 ). The third column u 3 is obtained in the same way :

such that the matrix

is unimodular. Therefore every decomposition exhibits at least one singularity defined by the vanishing of A. Moreover, it is readily seen that the following 0 degree choice

is singular if, and only if, A = 0. We thus conclude that P(F) is hyper-regular if and only if A = 0.

Finally, the equation A = ẋ cos θ + ẏ sin θ = 0, combined with F = ẋ sin θẏ cos θ = 0 leads to ẋ = ẏ = 0. We therefore have shown that S F = { ẋ = ẏ = 0}, in other words that the only obstruction to the hyper-regularity of P(F) is a flat output singularity, hence intrinsic according to Theorem 1.