Jean-Michel André
email: jean-michel.andre1@upmc.fr

Karine Le Guen

Philippe Jonnard

X-ray Quasi-Lamellar Etched Multilayers: Analysis by Coupled-Mode Theory

Keywords: X-ray mirrors, Multilayers, Gratings, Reflectivity, Coupled-mode theory

An etched multilayer (EM) is a two-dimensional structure obtained by etching a periodic multilayer following the profile of a grating. Generally a laminar grating shape, i.e. with a rectangular shape, is desired. However, we observe that different EM give different diffraction efficiency curves. We implement the coupled-mode theory to calculate the diffraction pattern of the EM as a function of the shape of the etching profile. Thus, we try to correlate the variability of the experimental reflectivity curves to various shapes of the profiles, rectangular and trapezoidal, and the corresponding simulated reflectivity curves.

Introduction

Multilayer gratings or etched multilayers (EMs), consisting in a periodic stack of bilayers with nanometric thickness laterally etched according to a periodic pattern, belong to this class of twodimensional photonic crystals that has opened up new ways for x-ray optics. Thus lamellar (or laminar) multilayer gratings can be used as efficient polychromators [1] and are also implemented as improved monochromators with respect to mere multilayer Bragg mirrors for the soft-x-ray domain [2]- [START_REF] Van Der Meer | Proc. SPIE[END_REF] . Recently, we have demonstrated that EMs can be used as dispersive element of curved-crystal soft x-ray spectrometer. In this case an improvement of the spectral resolution by a factor 4 was obtained in the C Kα emission range [START_REF] Jonnard | [END_REF] .

Calculations of the performances (efficiencies, band widths, …) of these devices have been initially done in the framework of the "rigorous" modal theories which give rise to timeconsuming computations and tricky numerical problems [10]- [12] . Among the latter are the numerical instabilities [13] [14] and the mathematical difficulties in solving the transcendental equations to find the modal wavenumbers [15] [16] . The level of difficulty is increasing as soon as the profile of the pattern is not lamellar, that is when the profile of the multilayer bar varies with the depth.

To overcome these difficulties, it has been recently proposed to appeal to the well-known coupled-waves approach; this latter has been successfully performed for the lamellar EMs [17]- [20] . In this paper, we show that it is possible to implement the coupled-mode theory (CMT) to calculate the diffraction pattern of a quasi-lamellar grating with different profiles of the multilayer bar and we adopt this approach to study the influence of this profile. Quasi-lamellar means that the shape slightly differs from the rectangular profile: trapeze (bottom wider than top) and inverse trapeze shapes are considered. A series of EMs has been prepared from Mo/B 4 C mirrors. Their diffraction pattern has been recorded around 182 eV (corresponding to the B K emission band) and compared to that originating from CMT calculations. The CMT approach has also been used in Ref. [18] to consider deviation of the multilayer period as a function of the depth and in Ref. [20] to consider passivation and contamination layers, non-rectangular profiles and sidewall scalloping of the multilayer bar. Here, we show that the comparison between experimental and simulated reflectivity curves can help to understand how the profile shape influences the optical performances of etched multilayers.

Description of the multilayer grating and geometry

We consider the diffraction by an EM: the scheme of the grating and the geometry are given in Figure 1. The EM consists in a multilayer stack built up with N bilayers of thickness d; each bilayer is made up with a heavy material 1 of thickness gd and polarizability ! ! , and a light material 2 of thickness (1-g)d and polarizability ! ! . We consider abrupt interfaces so that the polarizability ! ! of the bilayer is given by:

! ! = ! ! + ! ! (! ! -! ! ) (1)
where

! ! = g + ! ! !!! !"# !! ! ! ! ! (2)
with

! ! = !!!"# ! ! ! ! ! ! ! ! ! ! (3)
stands for the Heaviside unit step function. The stack is etched to obtain a grating of period D with multilayer bar of width ! ! ! at the depth z. We will consider the case of grating with symmetrical shape with respect to the zaxis, the profile of which is given by:

! !, ! = ! ! + 2 !"#(! ! ! ! ) ! ! cos ! ! ! ! ! !!! (4)
Then the function !(!) allows us to deal with a grating whose bar profile is arbitrary. In the usual case of a laminar grating, one has ! ! = Γ ! . The two-dimensionnal polarizability of the grating under consideration is given by:

! !, ! = ! ! ! !, ! (5)

Equations of coupled-mode theory

We recall here the main equations developed on the basis of the CMT to solve our problem.

More details are given in Ref. [19] . The electric field ! !, ! in the structure obeys the following propagation equation for the s-polarization case:

D GD d gd z x 5 ! ! ! !, !, ! + ! ! 1 -! !, ! ! !, !, ! = 0 (6)
solution of which can be expanded as:

! !, !, ! = ℰ ! ! exp ! ! !! ! !! !! !! exp ! ! ! ! (7)
with

! !" = ! !"#$ !"#$ + ! ! ! ! (8)
and

! ! = ! !"#$ !"#$ ( 9
)
where k is the wavenumber in vacuum and θ and ϕ the grazing and azimuthal angles of the incident radiation, respectively. Following the principle of the CMT, the n th order of the field ℰ ! ! is searched as a superposition of two counter-propagating waves whose respective amplitudes F and B vary with z, that is:

ℰ ! ! = ! ! ! exp ! ! !! ! + ! ! ! exp -! ! !! ! (10)
where

! !! = ! ! -(! !" ! + ! ! ! ) (11)
For the sake of uniqueness, the additional condition is required:

! ! ! exp ! ! !! ! + ! ! ! exp -! ! !! ! = 0 ( 12
)
where the dot indicates the derivative with respect to z. Substituting Eqs. ( 10) and ( 12) into Eqs. ( 6) and ( 7), it follows that the amplitudes satisfy the following system of differential equations:

!! ! z !" = - ! ! ! 2 ! !! ! ! ! !!! ! ! ! z exp +! ! !! -! !! ! !! !!!! + ! ! z exp -! ! !! + ! !! ! (13)
In this study, we assume that there is a single plane wave incident on the grating, which gives the following boundary conditions:

! ! 0 = ! !,! and ! ! (!) = 0 (14)
Let us note that this theory can be applied to the p-polarization case provided that one replaces in Eq. ( 6) the electric field by the modified magnetic field ! !, ! = ! !, ! / 1 -! !, ! and the polarizability ! !, ! by the following modified susceptibility [17] :

! !, ! = ! !, ! - ! ! ! ! ∇ ! ! !,! !!! !,! + ! ! ! ! ∇! !,! ! !!! !,! ! (15)

Comparison with experiment

The studied EM is designed as a narrow bandwidth monochromator for the spectral range around The reflectivity of the multilayer (prior etching) has been measured at the 0.154 nm wavelength. The analysis of the reflectivity curve shows that the multilayer period is 6.11 nm, the Mo thickness being 2.28 nm and the B 4 C one 3.83 nm. The interfacial roughness is limited to 0.4 nm rms. Some fits have been done with the period thickness varying along the stack. In these cases, the maximum possible drift is 0.0002 nm per period, that is to say a difference of only 0.04 nm between the top and the bottom of the multilayer.

The diffraction efficiency, or reflectance, of the EM at the zero th diffraction order of the grating and 1 st order of the multilayer has been measured at 182 eV versus the glancing angle [START_REF] Jonnard | [END_REF] at the BEAR beamline of the ELETTRA synchrotron facility. The reflectivity curves of six different samples, a priori prepared in the same conditions but not in the same batch, are presented in Figure 2. They are shifted vertically for sake of clarity. A common general shape is observed: a main peak close to an angle of 35° corresponding to the 1 st Bragg peak of the multilayer, surrounded by secondary peaks due to the coupling with the grating diffraction orders. However, all the curves are different: the maximum of reflectance varies between 0.075 and 0.150; the relative intensities of the secondary peaks change as well as their positions relative to the Bragg peak. For example, the intensity of the first feature toward the low angles varies from 0.034 to 0.054 and its position can move by about 0.1°, between 34.7 and 34.8°, leading this structure to be observed as a shoulder when too close from the Bragg peak (zero th diffraction order of the grating). It has been checked that the shape of the curves is the same at different location on the sample surface, that is to say the reproducibility problem is not due to a surface inhomogeneity introduced by the process of the multilayer. To gain in understanding this variability of the reflectivity curves, some samples have been cleaved in order to observe their cross section by scanning electron microscopy (SEM), see Figure 3. Let us note that the samples devoted to the SEM observation are not the one on which reflectivity measurements have been carried out.

First, it is observed that the grating period of 1.0 µm is respected. However, depending on the sample, different shapes of the multilayer bars are observed. The grating is not perfectly laminar: multilayer bars can present the shape of a trapeze or some roundness. Thus the G ratio value, even if close to 0.25, can vary from 0.15 to 0.3 from the top to the bottom of the multilayer. Two origins of this problem, coming from the etching process itself have been identified: a bad contact between the mask and the substrate during the photolithography step;

the difficulty to stop the RIE process just at the very border between the stack and the silicon substrate. Indeed, there is a thin Cr adhesion layer between the multilayer and the substrate, which is not sensitive to the etching process. Then, if the etching process is continuing after the layer is reached, ions will attack the foot of the multilayer bar and thus will lead to a nonrectangular shape. We believe that this variability of profile shapes can explain some part of the shape variations of the reflectivity curves. Thus, the developed CMT analysis is applied to study to what extent the profile of the multilayer bar influences the shape of the diffraction pattern of an EM. Five modes, i.e. 11 terms (-5, …, 0, …, +5) were taken into account in the calculations.

Figure 4 shows the EM efficiency in the range of the B K emission as a function of the glancing angle calculated for different profile shapes and for two mean G ratios, 0.18 and 0.25. In the first case the maximum of the simulated curves are normalized to the experimental reflectivity of one of the analysed samples. As explained above, the shape of the profile is taken into account through the depth variation of the G ratio. The mathematical functions describing the considered shapes are mentioned below in the case of G equal to 0.18:

-rectangular (lamellar grating): ! ! = 0.18;

-trapezoidal with the summit narrower than the basis: ! ! = 0.18 + 10 !! z (z in nm);

-trapezoidal with the summit broader than the basis: ! ! = 0.18 -10 !! z (z in nm).

These different shapes of the multilayer bars lead to small but sensitive reflectivity curve shapes.

The maximum of reflectivity varies by 0.01 and the intensities of the secondary structures by 0.02 (about 30% relative variation) in the case of the secondary structure toward the high angles of the main peak. The comparison of the experimental and calculated reflectivity curves shows that the reflectance is overestimated in the calculations by roughly a factor of 2.5, see Fig. 4(b) without normalization. A possible explanation can be the fact that the roughness at the multilayer interfaces as well as on the side of the multilayer bar is not taken into account. Now comparing the shape of the curves and the positions and relative intensities of the main and secondary maxima, it can be seen that none of the considered shapes is able to fully reproduce the experiment. For example the feature experimentally observed at 35.0° is only a faint shoulder in the simulations. The trapeze shape with G equal to 0.18 is the case where the position of the different structures is in better agreement with the experiment. In this case the intensity ranges from 0.011 to 0.018 in the simulations and is 0.025 experimentally. We have also made some calculations with an asymmetrical shape of the profile of the multilayer bar, that is to say a mix of rectangular shape for its ascending part (the part of the profile toward the lower x) and trapezoidal shape for its descending part (the part of the profile toward the higher x). In this case, the relative intensity of the shoulder at 35.0° increases and the relative intensities of the secondary peaks changes more prominently with respect to the symmetrical shapes. However, this does not seem compatible with the rather symmetrical shapes observed on the SEM images, see Fig. 3.

Conclusion

We have implemented the coupled-mode theory to calculate the diffraction pattern of an etched multilayer as a function of the shape of the etching profile. The simulations are not able to fully reproduce the features of the experimental reflectivity curve, however the positions and relative intensities of the secondary structures due to the coupling with the grating diffraction orders are sensitive to the shape of the multilayer bar and also to the ratio of the multilayer bar width to the grating period. To improve the comparison between experiments and simulations it does not seem that considering a variation of the multilayer period along the surface could be considered because we have verify the reproducibility of the reflectivity measurements at different positions of the surface. A drift of periodicity along the stack could also be envisaged but is very limited in our case, 0.04 nm at maximum over 200 periods, which could lead to some broadening and some reflectivity loss [18] , and could explain a part of the observed discrepancy.

Experimentally, we observed in the boron K range that the reflectivity curves from different Mo/B 4 C samples prepared a priori in the same conditions exhibit different shapes.

Fig. 1 .

 1 Fig.1. Scheme of an EM with the multilayer bars having a rectangular profile (laminar grating). d is the period of the multilayer stack while D is the period of the grating. gd is the thickness of the heavy material. G(z)D is the width of the multilayer bar at the depth z, constant in this particular case. The crossing of the x-and z-axes marks the origin (x=0; z=0).

  180 eV corresponding to the B K emission band. It consists in 200 Mo/B 4 C bilayers deposited onto an ultra-smooth silicon substrate, then etched following a two main steps process: first deep UV photolithography; second fluorine reactive ion etching (RIE). The aimed grating period D is 1 µm and G ratio is 0.25. The multilayer period d and g the ratio of the Mo thickness to this period are equal to 6.2 nm and 0.35 respectively. Thus the total thickness of the multilayer bar is 1240 nm. The size of the samples is 10 mm x 10 mm.

Fig. 2 .

 2 Fig. 2. Diffraction efficiency (reflectance) of six different etched multilayers. The curves are shifted vertically by 0.02 for sake of clarity.

Fig. 3 .

 3 Fig. 3. Scanning electron microscopy images of cross sections of two different etched multilayers.

Fig. 4 .

 4 Fig. 4. CMT calculated reflectivity curves in the B K emission range for different profiles of the multilayer bars with the G ratio of the multilayer bar width to the grating period equal to 0.18 (a) and 0.25 (b). The calculations with the 0.18 G ratio have been normalized to the experimental reflectivity value.

Acknowledgments

Incoatec is thanked for providing us with the multilayers. J.-R. Coudevylle and N. Isac from IEF, Orsay, France, are thanked for their help during the etching process. F. Delmotte, from Institut d'Optique, Palaiseau, France, is thanked for performing the reflectivity experiment at 0.154 nm and its fit.

This comes from small inhomogeneity from one sample to another during the preparation that leads to different shapes of the multilayer bars as shown by the SEM images.


Jaroslav Ne

Patrice Ossona

D E Mendez

A NOTE ON CIRCULAR CHROMATIC NUMBER OF GRAPHS WITH LARGE GIRTH AND SIMILAR PROBLEMS

Keywords:

In this short note, we extend the result of Galluccio, Goddyn, and Hell, which states that graphs of large girth excluding a minor are nearly bipartite. We also prove a similar result for the oriented chromatic number, from which follows in particular that graphs of large girth excluding a minor have oriented chromatic number at most 5, and for the pth chromatic number χp, from which follows in particular that graphs G of large girth excluding a minor have χp(G) ≤ p + 2.

Introduction

The circular chromatic number χ c (G) of a graph G, which is a refinement of its chromatic number, has received much attention recently (see [START_REF] Zhu | Recent developments in circular colouring of graphs[END_REF] for a survey on recent developments).

Recall that the circular chromatic number χ c (G) of a graph G is the infimum of rational numbers n k such that there is a mapping from the vertex set of G to Z n with the property that for adjacent vertices are mapped to elements at distance ≥ k.

In general, graphs of large girth can have arbitrary given circular chromatic number:

Theorem 1 (Nešetřil and Zhu [START_REF] Nešetřil | Construction of graphs with prescribed circular colorings[END_REF]). For any rational r ≥ 2 and any positive integers t, l, there is a graph G of girth at least l such that G has exactly t rcolorings, up to equivalence. In particular, for any r ≥ 2 and for any integer l, there is a graph G of girth at least l which is uniquely r-colorable and hence it has χ c (G) = r.

However, if restricted to special classes of graphs, large girth graphs may be forced to have small circular chromatic number. For instance: Theorem 2 (Galluccio, Goddyn, and Hell [START_REF] Galluccio | High-girth graphs avoiding a minor are nearly bipartite[END_REF]). For any integer n ≥ 4, for any ǫ > 0, there is an integer g such that every K n -minor free graph G with girth at least g has χ c (G) ≤ 2 + ǫ.

The aim of this paper is to extend Theorem 2 to other classes, all with a bounded expansion. Classes with bounded expansion have been introduced in [START_REF] Nešetřil | Linear time low tree-width partitions and algorithmic consequences[END_REF][START_REF]Grad and classes with bounded expansion I. decompositions[END_REF], and are based on the requirement for the graph invariant ∇ r (G) to be bounded in the class for each r.

Denote by V (G) and E(G) the vertex set and the edge set of G. Also denote by

|G| = |V (G)| (resp. G = |E(G)|) the order of G (resp. size).
Let G, H be graphs with V (H) = {v 1 , . . . , v h } and let r be an integer. A graph H is a shallow minor of a graph G at depth r, if there exists disjoint subsets A 1 , . . . , A h of V (G) such that

• the subgraph of G induced by A i is connected and as radius at most r, • if v i is adjacent to v j in H, then some vertex in A i is adjacent in G to some vertex in A j . We denote [START_REF]Grad and classes with bounded expansion I. decompositions[END_REF][START_REF]Sparsity (graphs, structures, and algorithms[END_REF] by G ▽ r the class of the (simple) graphs which are shallow minors of G at depth r, and we denote by ∇ r (G) the maximum density of a graph in G ▽ r, that is:

∇ r (G) = max H∈G ▽ r H |H| The expansion of a class C is the function Exp C : N → R ∪ {∞} defined by Exp C (r) = sup G∈G ∇ r (G). A class C has bounded expansion if Exp C (r) < ∞ for each value of r.
For instance, the class D of all graphs with maximum degree 3 has Exp D (r) = (3/2)2 r , while a class C has uniformly bounded Exp C (that is: Exp C (r) ≤ C for some constant C, independently of r) if and only if there is a graph F that is a minor of no graph in C. The expansion function Exp C is indeed non-decreasing and no other general constraint exists on the grow rate of expansion functions. In particular, for every non-decreasing function f : N → N with f (0) = 2 there exists a class C with Exp C (r) = f (r) (see [START_REF]Sparsity (graphs, structures, and algorithms[END_REF], Exercice 5.1).

However, when studying properties of a class, it is sometimes the case that the frontier between classes that verify the property and those that do not can be (at least approximately) expressed by means of a threshold expansion function (see for instance [START_REF] Dvořák | Small graph classes and bounded expansion[END_REF], where it is proved that every class C with expansion Exp C (r) ≤ c r 1/3-ǫ is small, but that some non-small class C exists, for which Exp C (r) ≤ 6•3 √ r log(r+e) ). In such a setting, Theorem 2 can be restated as the fact that if Exp C (r) ≤ C (for some constant C) then for any ǫ > 0, there is an integer g such that every G ∈ C with girth at least g has χ c (G) ≤ 2 + ǫ. However, such a statement cannot be extended to all classes with bounded expansion, and even to all classes C with exponential expansion: there exists 3-regular graphs with arbitrary large girth and circular chromatic number at least 7/3 [START_REF] Hatami | Random cubic graphs are not homomorphic to the cycle of size 7[END_REF], thus the class D of all graphs with maximum degree 3 (which is such that Exp D (r) = (3/2)2 r ) does not have the property that large girth implies a circular chromatic number arbitrarily close to 2.

We show in this paper (Corollary 1) that exponential expansion is indeed a threshold, in the sense that the condition log Exp C (r)/r = o(1) (as r → ∞) is sufficient to ensure the conclusion of Theorem 2. Hence a natural problem arises, to make this threshold more precise:

Problem 1. What is the maximum real integer c such that for every class C with lim r→∞ log Exp C (r) r ≤ c
and every positive real ǫ > 0 there is an integer g with the property that every graph G ∈ C with girth at least g has χ c (G) ≤ 2 + ǫ?

It follows from [START_REF] Hatami | Random cubic graphs are not homomorphic to the cycle of size 7[END_REF] and Corollary 1 that 0 ≤ c < log 2, and we conjecture c = 0. Other application of our approach are Theorem 5, which concerns the oriented chromatic number and which is the subject of Section 3, and Theorem 7,.which concerns the pth chromatic number χ p and which is the subject of Section 5.

Path-degeneracy

One of the main tools used to address homomorphism properties of graphs (or directed graphs) with high-girth is the notion of path-degeneracy (see for instance [START_REF] Galluccio | High-girth graphs avoiding a minor are nearly bipartite[END_REF][START_REF] Nešetřil | Colorings and girth of oriented planar graphs[END_REF].

Definition 1. A graph G is p-path degenerate if there is a sequence G = G 0 , G 1 , . . . , G t of subgraphs of G such that G t is a forest, and each G i (i > 0) is obtained from G i-1
by deleting the internal vertices of a path of length at least p, all of degree two.

We shall relate path-degeneracy of large girth graphs in a class C to expansion properties of the class C. First we recall a result that will needed for the proof.

Theorem 3 (Diestel and Rempel [START_REF] Diestel | Dense minors in graphs of large girth[END_REF]). Let d be an integer and let G be a graph. If

girth(G) > 6d + 3 and δ(G) ≥ 3 then ∇ 2d (G) ≥ 2 d
We now explicit the connection between path-degeneracy and expansion properties.

(G ′ ) ≥ 2 d . As (C ▽ (p -1)/2) ▽ 2d ⊆ C ▽ (2dp + (p -1)/2) it holds ∇ ⌈lp/3⌉ (G) ≥ ∇ 2dp+(p-1)/2 (G) ≥ ∇ 2d (G ′ ) ≥ 2 d . Hence it holds log ∇ r (G) r ≥ 1 p • d log 2 2d + 1 > 1 5p
, what contradicts our assumption.

Circular Chromatic Number of Graphs with Large Girth

We shall make use of the following property.

Lemma 2 (Bondy and Hell [START_REF] Bondy | A note on the star chromatic number[END_REF]). For every graph G and integer k, the following are equivalent:

• χ c (G) ≤ 2 + 1 k • G → C 2k+1
Existence of a homomorphism to an odd cycle is also linked to path-degenacy.

Lemma 3 (Galluccio, Goddyn, and Hell [START_REF] Galluccio | High-girth graphs avoiding a minor are nearly bipartite[END_REF]). A p-path degenerate graph G admits a homomorphism to any odd circuit of length at most p + 1.

In other words, a p-path degenerate graph G has

χ c (G) ≤ 2 + 1 ⌊p/2⌋
.

We deduce the following extension of Theorem 2. Then for every positive real ǫ > 0 there is an integer g such that every graph G ∈ C with girth at least g has χ c (G) ≤ 2 + ǫ.

Oriented Chromatic Number of Graphs with Large Girth

Recall that the oriented chromatic number χ o (G) of a (simple) graph G is the minimum k such that every orientation of G admits a homomorphism into some simple digraph with k vertices.

It was proved in [START_REF] Nešetřil | Colorings and girth of oriented planar graphs[END_REF] that there are planar graphs with arbitrarily large girth having oriented chromatic number 5, and that every planar graph with girth at least 16 has oriented chromatic number at most 5. The bound on the girth was reduced to 13 [START_REF] Borodin | Homomorphisms from sparse graphs with large girth[END_REF] (where the result generalized to graphs embeddable on the torus, or the Klein bottle) and then to 12 [START_REF] Borodin | Oriented 5-coloring of sparse plane graphs[END_REF]. In [START_REF] Borodin | Homomorphisms from sparse graphs with large girth[END_REF], it is proved that the considered highgirth graphs not only have oriented chromatic number 5, but that every orientation of these have a homomorphism to the same 5-vertex regular tournament C 2 5 . This tournament is a particular case of the (circular) digraphs C d n , which are defined as the digraphs with vertex set Z n and arcs (i, j) when j -i ∈ 1, . . . , d (mod n).

The following lemma relates p-path degeneracy and existence of homomorphism from a directed graph to C d n . Lemma 4 (Nešetřil, Raspaud, and Sopena [START_REF] Nešetřil | Colorings and girth of oriented planar graphs[END_REF]). Then end-vertices of every oriented path P of length at least (n -1)/(d -1) can be mapped by a homomorphism P → C d n to any pair of (not necessarily distinct) vertices of C d n . We can now relate expansion properties of a class to the oriented chromatic number of large girth graphs in the class. .

Then there is an integer g such that for every orientation G of a graph G in C with girth at least g there exists a homomorphism G → C d n .

Proof. Under the assumptions of the theorem there exists, according to Lemma 1 an integer g such that every graph G ∈ C with girth at least g is (n -1)/(d -1)path degenerate. For every orientation G of G, the existence of a homomorphism f : G → C d n the easily follows by an iterative use of Lemma 4. From this theorem immediately follows: Theorem 6. Let C be a class of graphs such that lim sup

r→∞ log Exp C (r) r ≤ 1 20
Then there is an integer g such that every G ∈ C with girth at least g has oriented chromatic number at most 5 (and actually are homomorphic to the tournament C 2 5 ). Proof. This follows from Theorem 5, by considering d = 2 and n = 5.

In particular, we have Corollary 2. Let C be a proper minor closed class of graphs. Then there is an integer g such that every G ∈ C with girth at least g has χ o (G) ≤ 5.

Centered Coloring of Graphs with Large Girth

Recall that the tree-depth td(G) of a graph G is the minimum height of a rooted forest F such that G is the subgraph of the closure of F . For positive integer p, the pth chromatic number χ p (G) is the least number of colors in a vertex coloring of G such that every i ≤ p colors induce a subgraph with tree-depth at most i. A p-centered coloring of a graph G is a vertex coloring such that, for any (induced) connected subgraph H , either some color c(H) appears exactly once in H, or H gets at least p colors. For detailed properties of tree-depth and χ p we refer the reader to [START_REF]Sparsity (graphs, structures, and algorithms[END_REF]. The χ p invariants are related to p-centered coloring as follows.

Lemma 5 (Nešetřil and Ossona de Mendez [START_REF]Tree depth, subgraph coloring and homomorphism bounds[END_REF]). For every graph G and every integer p, χ p (G) is the minimum number of colors in a (p + 1)-centered coloring of G. Lemma 6. Let G be a graph and let p be a positive integer. If G is 2p-path degenerate, then G has a p-centered coloring with at most p + 1 colors.

Proof. According to the definition of p-degeneracy, there is a sequence

G = G 0 ⊃ G 1 ⊃ • • • ⊃ G t of subgraphs of G such that G t is a forest, and each G i (i > 0) is obtained from G i-1 by
deleting the internal vertices of a path P i of length at least p, all of degree two. For each 0 < i ≤ t select an internal vertex v i of P i at distance at least p (in P i ) from both endvertices of P i . and orient all the edges of P i toward v i . Let C = {v i : 0 < i ≤ t}. Then G -C is a rooted forest, with edges oriented from the roots. Color the vertices G -C with p colors inductively as follows: each root is colored 0, and each other vertex u is colored as (c + 1) mod p, where c is the color assigned to the father of u. Note that this coloring is obviously a p-centered coloring of G -C. Now color the vertices in C with color p. It is easily checked that every path linking two vertices in C already got at least p -1 colors on the internal vertices by the above coloring of G -C. It follows that the every (p -1)-colored connected component of G is either included in G-C (hence has a uniquely colored vertex), or contains exactly one vertex in C, that is exactly one vertex colored p. Thus the defined coloring is a p-centered coloring of G.

The following theorem now easily follows. .

Then there is an integer g such that every graph G in C with girth at least g has χ p (G) ≤ p + 2.

Proof. According to Lemma 1 there is an integer g such that every graph G ∈ C with girth at least g is (2p + 2)-path degenerate thus, according to Lemma 6, such that χ p (G) ≤ p + 2.

In particular, we have Corollary 3. Let C be a proper minor closed class of graphs. Then for every integer p there is an integer g such that every G ∈ C with girth at least g has χ p (G) ≤ p+2.

As χ p (P n ) = p + 1 for sufficiently large n, it follows that the upper bound cannot be lowered under p + 1. The bounds χ 1 (G) ≤ 3 and χ 2 (G) ≤ 4 are optimal, as witnessed by a long odd cycle and the series-parallel graphs depicted below, respectively.

We leave as a problem the question whether the bound p + 2 is tight for every integer p.

Concluding Remarks

Remark that the property of a class that high girth graphs in the class are nearly bipartite or that they have oriented chromatic number at most 5 is not "topological" in the sense that satisfaction (or non-satisfaction) of this property is not preserved when the graphs in the class are subdivided: Consider an arbitrary class of graphs C. Construct the class C ′ as follows: for each graph G ∈ C (with |G| vertices) put in C ′ the |G|-subdivision of G. Then the class C ′ has the property that high girth graphs in the class are nearly bipartite and have oriented chromatic number at most 5, although this does not have to hold for C. In this sense, considering classes defined by a forbidden minor does not seem to be here optimal.

It is possible, from a class C, to construct the sequence (which we call resolution of C) of the classes of shallow minors of graphs in C at increasing depth

C ⊆ {G ▽ 0 : G ∈ C} ⊆ • • • ⊆ {G ▽ r : G ∈ C} ⊆ . . .
and we note that the results we obtained in this paper are expressed by means of the growth rate of the average degrees of the graphs in these classes, as r → ∞.

A related approach consists in considering the growth rate of the clique number of the graphs in these classes, as r → ∞. For a graph G and an integer r define ω r (G) = max{ω(H) : H ∈ G ▽ r}.

Then we have the following connection between the grow rate of ω r and the existence of sublinear vertex separator in a class of graphs: Then the following holds As ω r (G) ≤ 2∇ r (G) + 1, the classes C with lim r→∞ (log Exp C (r))/r = 0, which we considered in this note, satisfy the conditions of Theorem 8. However, the condition on ω r is weaker, as witnessed by classes of d-dimensional meshes with bounded aspect ratio (see [START_REF] Miller | Separators in two and three dimensions[END_REF] for a definition), that are such that sup G∈C ω r (G) grows polynomially with r (see [START_REF] Plotkin | Shallow excluded minors and improved graph decomposition[END_REF]).

We believe that the study of the growth rate of graph parameters on the resolution of a class can be useful to determine threshold where certain class properties stop being true.

Theorem 4 .≤ ǫ 10 .

 410 Let C be a class of graphs and let ǫ > 0. If it holds lim inf r→∞ log Exp C (r) r Then there is an integer g such that every graph G ∈ C with girth at least g has χ c (G) ≤ 2 + ǫ. Proof. Let p = ⌊2/ǫ⌋. Let g = [3r/p]. According to Lemma 1, every graph G ∈ C with girth at least g is p-degenerate. Thus, according to Lemma 3, every graph G ∈ C with girth at least g is such that χ c (G) ≤ 2 + ǫ. Corollary 1. Let C be a class of graphs such that lim r→∞ log Exp C (r) r = 0.

Theorem 5 .

 5 Let C be a class of graphs, and let n, d ∈ N. Assume it holds log Exp C (r) r ≤ d -1 5(n -1)

Theorem 7 .

 7 Let C be a class of graphs, and let p ∈ N. Assume it holds log Exp C (r)

Theorem 8 (

 8 [START_REF] Nešetřil | Linear time low tree-width partitions and algorithmic consequences[END_REF][START_REF]Sparsity (graphs, structures, and algorithms[END_REF]). Let C be a class of graphs such thatlim r→∞ log sup G∈C ω r (G) r = 0.Then the graphs of order n in C have vertex separators of size s(n) = o(n) which may be computed in time O(ns(n)) = o(n 2 ).Note that this threshold is sharp, in the sense that for every c > 0 there is a class C with limr→∞ log sup G∈C ω r (G) r ≤ cbut no sublinear vertex separators (consider the class C of k-subdivisions of cubic graphs).

  Lemma 1. Let C be a class of graphs, let p, g ∈ N, and let r = ⌈pg/3⌉. If it holds

	log Exp C (r) r	≤	1 5p	.

Then every graph G in C with girth at least g is p-path degenerate. Proof. Our proof follows similar lines as the proof ([6], Lemma 2.5). Assume for contradiction that C include graphs with girth at least g that are not p-path degenerate, and let G be a minimal such graph. The graph G is 2-connected, since a graph is p-path degenerate if all its blocks are. The graph G is not a circuit since, as l ≥ p + 1, the graph G would be p-path degenerate. Hence G is neither an edge or a circuit. Therefore, there exists a unique graph G ′ with minimum degree at least 3, which is homeomorphic to G. The graph G can be obtained from G ′ by replacing each edge e ∈ E(G ′ ) with a path P (e) of length at least one and at most (p -1). Thus G ′ ∈ C ▽ (p-1)/2. Let d = (g -3)/6. According to Theorem 3, if G ′ has girth greater than l then ∇ 2d

Supported by grant ERCCZ LL-1201 and CE-ITI, by the European Associated Laboratory "Structures in Combinatorics" (LEA STRUCO) P202/12/6061, and partially supported by ANR project Stint under reference ANR-13-BS02-0007.


Soizic Terrien
email: terrien@lma.cnrs-mrs.fr

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

Patricio De

La Cuadra
email: pcuadra@uc.cl

Benoît Fabre
email: benoit.fabre@upmc.fr

Is the jet-drive flute model able to produce modulated sounds like Flautas de Chinos ?

Flautas de chinos -prehispanic chilean flutes played during ritual celebrations in central Chile -are known to produce very particular beating sounds, the so-called sonido rajado. Some previous works have focused on the spectral analysis of these sounds, and on the input impedance of the complex resonator. However, the beating sounds origin remains to be investigated.

Throughout this paper, a comparison is provided between the characteristics of both the sound produced by flautas de chinos and a synthesis sound obtained through timedomain simulation of the jet-drive model for flute-like instruments. Jet-drive model appears to be able to produce quasiperiodic sounds similar to sonido rajado. Finally, the analysis of the system dynamics through numerical continuation methods allows to explore the production mechanism of these quasiperiodic regimes.

INTRODUCTION

The Flautas de Chinos are prehispanic chilean flutes still performed nowadays in ritual celebrations. The specificity of these intruments lies both in the production of very particular sounds (the so-called sonido rajado, literally "torn sounds") and in the complex form of their resonator (highlighted in the schematic representation of the instrument provided in figure 1), formed by two pipes in series with different diameters. Early works by Wright and Campbell [START_REF] Wright | Analysis of the sound of chilean pifilca flutes[END_REF] have focused on the characteristics of the produced sounds (in terms of fundamental frequency, waveform and spectrum). Blanc et al. [START_REF] Blanc | Acoustics of the "flautas de chinos[END_REF] then attempted to describe and model the complex resonator. However, the production mechanism of sonido rajado remains unexplained. Moreover, to the authors knowledge, the behaviour of such instruments has never been investigated through the study of a physical model of sound production. We do not develop in this paper an accurate, specific physical model of flautas de chinos, but we propose to evaluate the ability of the classical jet-drive model for flute-like instruments to produce regimes similar to sonido rajado. Moreover, such a study allows us to investigate the generation mechanism of those particular regimes.

The second section focuses on the description of the sound produced by flautas de chinos. In section 3, we recall the state-of-the-art model for flute-like instruments, whose capacity to produce modulated sounds is evaluated in section 4 using time-domain simulation. Finally, we examine in section 5 the origin of such regimes.

ANALYSIS OF SONIDO RAJADO

The "sonido rajado" produced by Flautas de Chinos seems to be the combination of a very loud sound, a rich harmonic structure, and an outstanding beating component [START_REF] Wright | Analysis of the sound of chilean pifilca flutes[END_REF][START_REF] Blanc | Acoustics of the "flautas de chinos[END_REF]. These beats can evoke the particular features of some quasiperiodic regimes, also produced for example in recorders or transverse flutes [START_REF] Coltman | Jet offset, harmonic content, and warble in the flute[END_REF][START_REF] Fletcher | Sound production by organ flue pipes[END_REF].

The sound produced by a musician playing a particular flauta de chinos named puntera -whose dimensions are specified in figure 1-has been recorded through a microphone located outside the instrument, and is represented in figure 2. One can note that this radiated sound shows strong envelope modulations, whose modulation frequency can be roughly approximate between 13 and 15 Hz.

To understand the production mechanism of these particular sounds, the knowledge of the regime nature is an essential information. However, the study of spectrum and waveform may be insufficient to distinguish, for example, a quasiperiodic regime from a chaotic regime. The computation of the Poincaré section (see for example [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF][START_REF] Gibiat | Phase space representations of acoustical musical signals[END_REF]) is a commonly used method that allows to reveal the nature of the regime. However, in the present case, both the small length of the recorded signal and its nonstationarity prevent the use of this tool.

Return times can be computed to estimate periodicity at different time scales in embedded time series. The socalled return times are the times for which the time series returns to its location in the pseudo phase space [START_REF] Zou | Exploring reccurences in quasiperiodic dynamical systems[END_REF]. When this is done by chosing any point of the time series as a reference point, an histogram can be computed, which represents the occurrence of the different return times. Strong peaks in the histogram indicate periodicity. According to Slater's theorem [START_REF] Zou | Exploring reccurences in quasiperiodic dynamical systems[END_REF], for a noise-free signal, a quasiperiodic regime is characterized by the presence of three peaks in the histogram, corresponding to three "return times". Two of them are respectively related to the "base frequency" and to the "modulation frequency", whereas the third is the sum of the two others. In experimental time series, noise is inevitable. Its influence has been studied on return times by Zou [START_REF] Zou | Exploring reccurences in quasiperiodic dynamical systems[END_REF] who concludes that other peaks are likely to occur, but the sum relationship still holds.

The histogram of the recorded sound, displayed in figure 3, shows two principal peaks, respectively for the return times t 1 = 0.0054 s and t 2 = 0.0767 s. It highlights the presence, in the recorded sound, of two incommensurate frequencies: the "base frequency" f 1 = 185 Hz and the "modulation frequency" f 2 = 13.04 Hz. All other return times of the histogram are either multiples of t 1 or combinations of the form t 2 + p • t 1 , with p an integer. The histogram structure thus seems to satisfy Slater's theorem. Indeed, in the case of a chaotic regime, other peaks should appear in the histogram, related to return times which are neither multiples nor linear combinations of t 1 and t 2 .

Therefore, the computation of return times allows us to conclude about the quasiperiodic nature of the sonido rajado produced by flautas de chinos.

MODEL OF FLUTE-LIKE INSTRUMENT

General mechanism of sound production

Although instruments such as transverse flutes, flue organ pipes or the prehispanic chilean flutes studied here may appear different from each other, they involve the same general mechanism of sound production. When the musician blows into the instrument, a naturally unstable jet is created at a channel exit. This channel is formed by the player's lips for transverse flutes and flautas de chinos, whereas it is a part of the instrument for recorders and flue organ pipes (see figure 1). This jet oscillates around an edge called "labium". The resulting interaction produces acoustic energy in the resonator, constituted by the air column contained in the instrument. In turn, the acoustic field thus created in the instrument perturbs the jet at the channel exit. This perturbation being convected and amplified from the channel exit toward the labium, it sustains the oscillation of the jet around the labium, and thus allows the emergence of self-sustained oscillations. This functioning is detailed in [START_REF] Chaigne | Acoustique des instruments de musique (Acoustics of musical instruments)[END_REF], whose main elements are recalled in sections 3.2 and 3.3.

Exciter

Hydrodynamics of the unstable jet

Once a permanent operating regime of the instrument is reached, the perturbation of the jet (the so-called receptivity) is provided by the acoustic field, at the channel exit. According to the empirical model proposed by de la Cuadra [START_REF] De La Cuadra | Visualization and analysis of jet oscillation under transverse acoustic perturbation[END_REF], the receptivity is represented as a transverse displacement η(0, t) of the jet at the channel exit:

η(0, t) = h U j v ac (t), (1)
where h is the channel height (i.e. the opening height of the musician's lips), v ac is the acoustic velocity in the pipe, and U j is the jet central velocity (directly related to the fact that the musician blows hard or not).

As the jet is naturally unstable, this perturbation is amplified while convecting along the jet, from the channel exit to the labium. Based on the work of Rayleigh concerning instability of an infinite jet [START_REF] Rayleigh | The theory of sound second edition[END_REF], this phenomenon can be described in a simplified way as an exponential amplification of the perturbation with respect to the convection distance x:

η(x) = η(0)e αix , (2)
Based on experimental works by de la Cuadra [START_REF] De La Cuadra | The sound of oscillating air jets: Physics, modeling and simulation in flute-like instruments[END_REF], an approximation of the amplification parameter α i through the empirical expression:

α i ≈ 0.4/h, (3)
seems reasonable regarding the range of the Strouhal number Str = f h Uj considered throughout the paper. Finally, we can express the jet transversal deflection at the labium:

η(W, t) = η(0, t -τ )e αiW = h U j v ac (t -τ )e αiW , ( 4
)
where W is the distance between the channel exit and the labium (corresponding to the length of the jet, highlighted in figure 1) and τ is the convection delay of the perturbation along the jet. Calling c p the convection velocity of this perturbation, τ is given by: τ = W cp . Both theoretical and experimental works [START_REF] De La Cuadra | Visualization and analysis of jet oscillation under transverse acoustic perturbation[END_REF][START_REF] Rayleigh | The theory of sound second edition[END_REF][START_REF] Nolle | Sinuous instability of a planar jet: propagation parameters and acoustic excitation[END_REF] have shown that in the general case 0.3U j ≤ c p ≤ 0.5U j , which leads to:

τ ≈ W 0.4U j (5)

Nonlinear jet-labium interaction

The interaction between the perturbed jet described above and the labium (see figure 1) makes the jet oscillate from one side to another of the labium. This oscillation causes a flow injection alternatively inside and outside the instrument. Based on the jet-drive model proposed by Coltman and followed by Verge [START_REF] Coltman | Jet drive mechanisms in edge tones and organ pipes[END_REF][START_REF] Verge | Sound production in recorder-like instruments. ii. a simulation model[END_REF], these two (localised) flow injection in phase opposition are modelised as a dipolar pressure source:

∆p = - ρδ d W H • dQ in dt , (6)
In this equation, ρ is the air density, H the width of the jet and we note, as proposed by Verge [START_REF] Verge | Jet oscillations and jet drive in recorder-like instruments[END_REF], δ d the effective distance between the two flow injection points. Based on both theoretical works and empirical estimation of parameters, Verge [START_REF] Verge | Jet oscillations and jet drive in recorder-like instruments[END_REF] proposed the approximation δ d ≈ 4 π (2hW ), which is adopted here. Q in represents the flow injection in the pipe [START_REF] Chaigne | Acoustique des instruments de musique (Acoustics of musical instruments)[END_REF], given by:

Q in = H y0-η(t) -∞ U (y)dy, (7)
In this expression, y 0 is the offset between the labium and the jet centerline, and U (y) the velocity profile of the jet, representing by a Bickley profile:

U (y) = U j sech 2 ( y b ), ( 8
)
where b is the half width of the jet, related (under certain assumptions detailed for example in [START_REF] Ségoufin | Experimental study of the influence of the mouth geometry on sound production in a recorderlike instrument: windway length and chamfers[END_REF]) to the height of the channel exit through b = 2h/5. By injecting equations ( 7) and ( 8) in equation ( 6), we finally obtain the expression of the pressure source that excites the resonator:

∆p src (t) = ρδ d bU j W • d dt tanh η(t) -y 0 b . (9)

Nonlinear losses

Between the channel exit and the labium of flue instruments, the presence of an important transversal flow induced by the acoustic field in the pipe can cause, for high acoustic velocities, vortex shedding at the labium [START_REF] Fabre | Vortex shedding in steady oscillation of a flue organ pipe[END_REF]. This phenomenon is modeled as nonlinear losses, and so represented by an additional nonlinear term ∆p los in equation ( 9):

∆p(t) =∆p src (t) + ∆p los (t) = ρδ d bU j W d dt tanh η(t) -y 0 b - ρ 2 v ac (t) α vc 2 sgn(v ac (t)) (10)
where α vc ≈ 0.6 is a vena contracta factor.

Resonator

The acoustical response of the air column contained in the pipe, excited by the pressure source described above, is represented through the input admittance Y in = V ac /∆P of this resonator. V ac and ∆P are respectively the frequencydomain expressions of the acoustic velocity at the pipe inlet and the pressure source. The use of a modal decomposition of Y in is interesting as it allows an independant control of the different resonance mode characteristics. In the frequency domain, Y in is thus represented as a sum of resonance modes:

Y in = m a m jω ω 2 m -ω 2 + jω ωm Qm , (11)
where a m , ω m and Q m are respectively the modal amplitude, the resonance pulsation and the quality factor of the m th resonance mode.

As the particularity of the flautas de chinos lies in the shape of their resonator, we use a modal decomposition of the measured admittance of the "puntera" flute whose sound is studied in section 2. The measure was realised thanks to an impedance sensor for wind instruments [START_REF] Macaluso | Trumpet with nearperfect harmonicity: Design and acoustic results[END_REF]. The corresponding modal parameters (displayed in table 1) of the first five resonance modes are estimated by fitting the admittance calculated using equation ( 11) to the measured one. One can note that the agreement between the original admittance and the fitted one, which are both displayed in figure 4, can be improved by increasing the number of resonance modes taken into account in the fit process. However, it considerabily increases the computation cost for the foregoing different resolution methods.

Numerical resolution of the model

The complexity of the model, due to its neutral nature and to the presence of nonlinear terms, imposes the use of numerical methods for solving the model. The neutral nature is related to the presence of a delayed derivative term, as discussed in [START_REF] Terrien | Numerical resolution of a physical model of flute-like instruments: comparison between different approaches[END_REF]. In this study, we compare the results provided by two complementary methods: a classical timedomain simulation scheme, and a calculation method of mode number periodic solutions coupled with a numerical continuation algorithm. Both are presented below and used in sections 4 and 5.

a m ω m Q m 1 11

Time-domain simulations

Time-domain simulations are performed in Matlab/Simulink, with a classical fourth-order Runge-Kutta method. Such an approach is interesting as it allows to access many kind of stable stationary regimes (for example static, periodic, quasiperiodic or chaotic) and to transients between these regimes. However, results are very sensitive to initial conditions and to numerical parameters such as the sampling frequency. Moreover, unstable solutions cannot be computed: a slight disturbance (noise or even round-off errors, for example) is enough to make the system move away from these solutions. Therefore, it can be difficult to access, with such a method, to a global knowledge of the different oscillation regimes, essential to understand the system dynamics.

Orthogonal collocation and numerical continuation

Orthogonal collocation method allows to compute periodic solutions of dynamical systems. Its principle is based on the discretization of a single (unknown) period of the solution x on n intervals. On each interval, the solution is approximated by a polynomial of degree d. Projecting the model equations on this set of representation points leads to an algebraic system, whose unknowns are the value of x at each discrete point, and the value of the period T .

This method, which allows to determine a solution x 0 for a set of model parameters λ 0 , is coupled to a numerical continuation algorithm. Starting from x 0 , it computes, using a Newton-Raphson prediction-correction method, the neighboring solution corresponding to a set of parameter λ 0 + ∆λ [START_REF] Krauskopf | Numerical continuation methods for dynamical systems[END_REF]. Proceeding by successive iterations, one finally access to a complete branch of periodic solutions.

The computation of the Floquet multipliers allows subsequently to determine the stability properties of each point of the branch [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF]. Indeed, according to Floquet theory, a periodic solution (i.e. a point of the branch) is stable until all its Floquet multipliers lies in the unit circle. When a periodic solution loses its stability, the resulting regime observed after this bifurcation point depends on which way the Floquet multipliers leave the unit circle at the bifurcation point. For example, a quasiperiodic regime may appear if two complex conjugate Floquet multipliers leave the unit circle [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF].

Giving access to the bifurcation diagram, which represents in an ideal case all the periodic and static solutions of the studied system, this method provides a more global knowledge of the system dynamics (coexistence of multiple solutions, unstable solutions ...), and thus permits an easier interpretation of phenomenon such as hysteresis or emergence of non periodic regimes (see for example [START_REF] Terrien | Numerical resolution of a physical model of flute-like instruments: comparison between different approaches[END_REF]).

The neutral nature of the model prevents the use of classical numerical continuation software (such as, for example, AUTO [START_REF] Doedel | AUTO: A program for automatic bifurcation analysis of autonomous systems[END_REF] or Manlab [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]). We use here DDE-Biftool [START_REF] Engelborghs | DDE Biftool: a Matlab package for bifurcation analysis of delay differential equations[END_REF], a software specifically developped for numerical continuation in delay differential equations, and its extension for neutral systems [START_REF] Barton | Collocation schemes for periodic solutions of neutral delay differential equations[END_REF].

TIME-DOMAIN SIMULATION: CHARACTERISTICS OF THE SYNTHESIS SOUND

Although we study here a very generic model of flute-like instruments, time-domain simulations lead to non periodic regimes which can perceptually recall the sonido rajado produced by flautas de chinos. Because of the few studies about parameter estimation in flautas de chinos, some of the model parameters are chosen arbitrarily close to those currently uses in the case of a model of an alto recorder [START_REF] Terrien | Numerical resolution of a physical model of flute-like instruments: comparison between different approaches[END_REF], except for the excitation window area W H and the labium offset y 0 which are here slightly higher than in the case of an alto recorder. As these parameters are fixed by the player in the case of flautas de chinos, the chosen values (W = 1cm, H = 6mm and y 0 = 0.2 mm) seem reasonable. Parameter values used throughout the paper are provided in table 2.

A first simulation is achieved using a constant jet velocity U j = 39 m/s, which would correspond to a case where the musician blows in a steady way. This value, which would be equivalent, using the stationary Bernoulli equation, to a mouth pressure of about 915 Pa, is of the same order of magnitude as the values measured in the player's mouth for the puntera flute (between 886 Pa and 1546 Pa).

The computed acoustic velocity at the pipe inlet v ac , represented in figure 5 as a function of the time, shows, as the real sound (see figure 2), envelope modulations. A Fourier analysis of the signal envelop provides the modulation frequency f mod = 8.6 Hz. One can note that this value is of the same order of magnitude as the modulation frequency of the real sound, f 2 = 13.04 Hz (see section 2), and coherent with the observations of Wright and Campbell [START_REF] Wright | Analysis of the sound of chilean pifilca flutes[END_REF]. In the case of time-domain simulation, this frequency strongly depends on the the jet velocity U j , and a slight shift of this parameter is enough to achieve the experimental value of 13.04 Hz. However, the fact that the modulation frequency is related to the jet velocity is a significant difference compared to the real instruments behaviour. Indeed, as observed by Wright and Campbell [START_REF] Wright | Analysis of the sound of chilean pifilca flutes[END_REF], the modulation frequency shows only little variations with the pressure in the musician's mouth. Although the physical model is not specific to the instrument studied here, and that some parameters can only be roughly estimated, the implementation of a measured admittance allows us to compare qualitatively, in figure 6, the spectrum of the synthesis sound with that of the real sound. Such a comparison shows that both the fundamental and modulation frequencies (particularly visible around 370 Hz) of the synthetized sound are of the same order of magnitude as those of the real sound. One can note that the amplitudes can not be compared: indeed, we consider the spectrum of the sound radiated by the real instrument whereas we consider, for the model, the internal acoustic field (which would correspond to a measure under the labium of the instrument). Moreover, to make easier the qualitative comparison of the frequencies, the two spectra have been normalized with respect to their own maximum.

As for the real sound, one can wonder about the quasiperiodic or chaotic nature of the synthesis sound. The study of the Poincaré section allows us to distinguish these different kinds of regime from each other. As shown on figure 7, the three-dimensional Poincaré section of the synthesis signal reveals a densely packed closed curve, which is characteristics of a two-period quasiperiodic solution (see for example [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF]), and thus allows to conclude about the quasiperiodic nature of the simulated signal. As a first conclusion, one can note that this study highlights the ability of the state-of-the-art physical model of flute-like instruments to produce regimes where some of their features such as quasiperiodic nature, fundamental and modulation frequencies, are similar to those of the sound produced by a flauta de chinos.

QUASIPERIODIC REGIME: GENERATION MECHANISM

As it provides a more global knowledge of the system dynamics, the study of the bifurcation diagram of the model studied, which ideally represents all its stable and unstable periodic solutions as a function of one of its parameters, has recently shown its interest in the understanding of the functioning of musical instruments [START_REF] Karkar | Oscillation threshold of a clarinet model: a numerical continuation approach[END_REF].

The bifurcation diagram of the model displayed in figure 8 represents the oscillation frequency of the different periodic solutions as a function of the jet velocity U j . Solutions are computed for a range of the jet velocity U j for which quasiperiodic regimes may occur in time-domain simulations (see section 4). It shows the existence, in this range of U j , of two periodic solutions branches: the first corresponds to the first register of the instrument (that is to say, to oscillations emerging from an instability at a frequency close to the first resonance frequency of the res-onator), and the second is related to the second register. Solid and dotted lines respectively represent stable and unstable parts of periodic solution branches. In the previous section, we highlighted that time-domain simulation of the model reveals the existence of a quasiperiodic regime for U j = 39m/s. For this value, the bifurcation diagram predicts that the first register is unstable, whereas the second register is stable.

In order to understand the origin of this quasiperiodic solution, we computed a second time-domain simulation, during which the jet velocity U j follows a decreasing ramp, from U j = 39m/s, to U j = 32m/s. The result, displayed in figure 9, shows a gentle transition from the quasiperiodic regime studied in the previous section, to a periodic regime. A Fourier analysis of this periodic regime provides an oscillation frequency of 188.4 Hz, which is very close to the first resonance frequency ω1 2π = 184.4 Hz, and thus allows to conclude that this periodic solution corresponds to the first register. Thereby, the confrontation of figures 8 and 9 suggests that the quasiperiodic solution is the result of the loss of stability of the first register (predicted by the bifurcation diagram at U j = 35.6m/s -see figure 8) through a direct Neimark-Sacker bifurcation [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF], .

If the numerical tools used here do not allow computa-tion and continuation of quasiperiodic solution branches, a further analysis of the stability properties of the periodic solutions permits to analyse the birth of quasiperiodic solutions (see section 3.4). The Floquet multipliers associated to the first register, around the bifurcation point (U j = 35.6 m/s) are shown in figure 10 for both the point of the branch corresponding to U j = 35.3 m/s (just before the loss of stability) and the point corresponding to U j = 35.6 m/s (just after the loss of stability). One can clearly observe that the loss of stability is the result of the crossing of the unit circle by two complex conjugate Floquet multipliers. This analysis not only confirms the generation of a quasiperiodic regime, but it also provides the value of the modulation frequency at the quasiperiodic regime threshold (i.e at U j = 35.6 m/s), through the computation of the Floquet exponents γ f , which are related to the Floquet multipliers µ f through [START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF]:

γ f = 1 T ln(µ f ) ( 12
)
where T is the period of the periodic solution at the bifurcation point. According to Floquet theory, the imaginary part of the Floquet exponent corresponds to the modulation pulsation ω mod . Here, such an analysis of the Floquet multipliers predicts a modulation frequency: f mod = ω mod /2π = 26.4 Hz. A further study of the synthesis signal envelop represented in figure 9 shows a threshold of the quasiperiodic regime at about U j = 34.2 m/s (see a zoom of figure 9 in figure 11). As in section 4, a Fourier analysis of this envelop signal provides the modulation frequency related to the quasiperiodic regime: f modsimu = 26.3 Hz. Although the quasiperiodic regime threshold seems to be underestimated in the time-domain simulation (which may be a consequence of both the sampling and the dynamics of the control parameter U j [START_REF] Bergeot | Prediction of the dynamic oscillation threshold in a clarinet model with a linearly increasing blowing pressure[END_REF]), the modulation frequency shows good agreement with those predicted by the Floquet analysis.

CONCLUSION

This study, rather than proposing a specific physical model of flautas de chinos, focuses on the ability of the jet-drive As a first conclusion, a study of the sound produced by a flauta de chinos has highlighted the quasiperiodic nature of sonido rajado.

Secondly, time domain simulations show that not only the jet-drive model can produce, as transverse flutes (for example), periodic solutions associated to the first or second register of the instrument, but also quasiperiodic regimes, which present similarities with the sonido rajado of prehispanic chilean flutes.

Moreover, an analysis of the system dynamics through its bifurcation diagram highlights the production mechanism of those quasiperiodic regimes. Comparison between the sound produced by flautas de chinos and the model dynamics thus suggests that the desired behaviour of such instruments would correspond to a quasiperiodic regime emerging from a direct Neimark-Sacker bifurcation of the periodic solution branch related to the first register of the instrument.

However, a quantitative comparison between experimental and synthesis sounds is impossible here as we study on the one hand the sound radiated by the real instrument whereas we consider, in the model, the internal acoustic field. But more importantly, both the fact that the precise values of different parameters of the model are unknown and the fact that different phenomena are not taken into account (for example, the jet is certainly turbulent [START_REF] Blanc | Acoustics of the "flautas de chinos[END_REF]) limit such a quantitative comparison. These elements probably explain why significant differences remain between the real sound and the synthesis one, such as, for quasiperiodic regimes, the variation of the modulation frequency with the jet velocity.

  Copyright: c 2013 Soizic Terrien et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 Unported License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Figure 1 .

 1 Figure 1. Schematic representation of both the resonator of the studied flauta de chinos and the excitation mechanism in flute-like instruments.

Figure 2 .

 2 Figure 2. Radiated sound produced by a flauta de chinos "puntera".

Figure 3 .

 3 Figure 3. Histogram representing the return times of the sound produced by the "puntera" flauta de chinos.

Figure 4 .

 4 Figure 4. Measured and fitted input admittance of a puntera prehispanic flute. Only the five first resonance modes taken into account in the model are represented.

Figure 5 .

 5 Figure 5. Acoustic velocity computed through a timedomain simulation of the model for a constant jet velocity U j = 39 m/s, represented as a function of the time. The sampling frequency used for computation is 22x44.1 kHz.

Figure 6 .Figure 7 .

 67 Figure 6. Normalized spectrums of both the sound produced by a "puntera" flauta de chinos (represented in figure 2) and the synthesis sound represented in figure 5.

Figure 8 .

 8 Figure 8. Bifurcation diagram of the studied system, with two periodic solution branches corresponding to the first and the second registers. Abscissa: jet velocity U j (m/s). Ordinate: oscillation frequency (Hz). The horizontal dotted lines correspond to the two first resonance frequencies of the resonator.

Figure 9 .

 9 Figure 9. Acoustic velocity computed through a timedomain simulation, for a linear decreasing ramp of the jet velocity U j . The sampling frequency used for simulation is 22x44.1 kHz.

Figure 10 .

 10 Figure 10. Floquet multipliers represented in the complex plane, for two points of the first register periodic solution branch. The first point is located just before the loss of stability of the first register (U j = 35.3 m/s) and the second one just after (U j = 35.6 m/s).

Figure 11 .

 11 Figure 11. Zoom of figure 9: acoustic velocity computed through a time-domain simulation, for a linear decreasing ramp of the jet velocity U j . The sampling frequency used for simulation is 22x44.1 kHz.

Table 1 .

 1 Fitted modal coefficients corresponding to the measured admittance of the puntera flute.

		.39 1156.7	26
	2	7.05 2342.8 34.4
	3	9.55 4796.4 50.7
	4	8.12 5943.4 52.9
	5	12.93 8418.9	58

Table 2 .

 2 Parameter values used for numerical resolutions of the model.

	Parameter Numerical value (S.I.)
	α i		0.4/h = 400
	δ d	4 π	(2hW ) = 0.0057
	b		2h/5 = 0.0004
	y 0		0.0002
	h		0.001
	W		0.01
	ρ		1.2
	c p		0.4U j
	α vc		0.6