Martin Rosalie
email: martin.rosalie@uni.lu

Templates and subtemplates of Rössler attractors from a bifurcation diagram

Keywords: numbers: 05.45.-a Rössler system, bifurcation diagram, template, subtemplate

Templates and subtemplates of Rössler attractors from a bifurcation diagram

Introduction

Since 1976, the Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF] is well know for its simplicity (three differential equations with only one non linear term) and its dynamical richness producing chaos. Used as a basic system to demonstrate various properties of dynamical systems, this system is still a source inspiration for researchers. This system has been widely explored with several tools. The main goal of this paper is to extend the use of the topological characterization method to several chaotic attractors. We introduce a way to use the template as a global description that contains various attractors templates of the Rössler system.

In this paper, we will study this system in a parameter space to highlight common properties of neighbours attractors in this space. Castro et al. [START_REF] Castro | Characterization of the Rössler system in parameter space[END_REF] study of the parameter space of this system using Lyapunov exponents reflects its dynamics (stable, chaotic or trajectory diverging). The maps are built varying parameters a and c of the system. These maps display fractal structure and illustrate period doubling cascades. This principle is also employed by Barrio et al. [3] for the three parameters of the Rössler system. Their analysis of local and global bifurcations of limit cycles permits to a have a better understanding of the parameter space. Additionally, attractors with equilibrium points associated to their first return maps are also plotted to illustrate the various dynamics of this system, including coexistence of attractors. To enlarge this overview of recent work on the bifurcations and dynamics on the Rössler system, Genesio et al. [START_REF] Genesio | A global qualitative view of bifurcations and dynamics in the Rössler system[END_REF] use the first-order harmonic balance technique to study fold, flip and Neimark-Sacker bifurcations in the whole parameter space. Finally, the recent work of Sprott & Li [START_REF] Sprott | Asymmetric bistability in the Rössler system International Journal of Bifurcation and Chaos submitted[END_REF] introduces another way to reach coexisting attractors in addition to the cases identified by Barrio et al. [START_REF] Barrio | Qualitative analysis of the Rössler equations: Bifurcations of limit cycles and chaotic attractors[END_REF].

In this paper we study a bifurcation diagram of the Rössler system exhibiting various dynamics using topological properties of attractors. We use the topological characterization method based on the topological properties of the attractor's periodic orbits [START_REF] Gilmore | The Topology of Chaos[END_REF]. The purpose is not only to obtain templates of chaotic attractors but also to find common points or properties as it as already been shown for this system by Letellier et al. [START_REF] Letellier | Unstable periodic orbits and templates of the Rössler system : toward a systematic topological characterization[END_REF] for a "funnel attractor". In this particular case, a linking matrix describes the template depending on the number of branches. In this paper, we will explore a bifurcation diagram and show that only one template contains all the templates of attractors as subtemplates.

This paper is organized as follow. The first part introduces the Sprott & Li [START_REF] Sprott | Asymmetric bistability in the Rössler system International Journal of Bifurcation and Chaos submitted[END_REF] work with their bifurcation diagram. The second part details the topological characterization method; eight attractors are studied and their templates are obtained. Then we prove that the eight templates are subtemplates of a unique template. It describes the topological structure of all the attractors of the entire bifurcation diagram. Finally we provide a partition of the bifurcation diagram giving the symbolic dynamic associated with the unique template depending on the bifurcation parameter.

Bifurcation diagram

Barrio et al. [START_REF] Barrio | Qualitative analysis of the Rössler equations: Bifurcations of limit cycles and chaotic attractors[END_REF] highlight the fact that a Rössler system can have two coexisting attractors as solutions for a set of parameters. Sprott & Li [START_REF] Sprott | Asymmetric bistability in the Rössler system International Journal of Bifurcation and Chaos submitted[END_REF] parametrize the Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF] with the parameter α

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ẋ = -y -z ẏ = x + ay ż = b + z(x -c) with ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ a = 0.2 + 0.09α b = 0.2 -0.06α c = 5.7 -1.18α (1)
in order to explore bifurcations. When α = 1, two attractors coexist in the phase space. We reproduce their bifurcation diagram when α varies. The value of the fixed points of the system are

S ± = x ± = c ± √ c 2 -4ab 2
y ± = -c ∓ √ c 2 -4ab 2a z ± = c ± √ c 2 -4ab 2a .
(

) 2
The bifurcation diagram Fig. 1 is obtained using the following Poincaré section

P ≡ (y n , -z n ) ∈ R 2 -x n = -x - (3)
where x -is the x value of the fixed point S -(see [START_REF] Letellier | Toward a general procedure for extracting templates from chaotic attractors bounded by high genus torus[END_REF] for details on this Poincaré section). The uses of this Poincaré section explains why Fig. 1 is similar to FIG. 5 This diagram indicates parameter for the Rössler system where the solution is a limit cycle or chaotic. This diagram exhibits a doubling period cascade that is a classical route to chaos for -2 < α < 0.2. This is followed by a chaotic puff (α = 0.5) and by various regimes (banded chaos and almost fully developed chaos). We choose representative values of α where one or two attractors are solutions of the system

A α = -0.25 E 1 α = 1 B α = 0.5 E 2 α = 1 C α = 0.78 F α = 1.135 D α = 0.86 G α = 1.22 . (4)
We analyse these attractors using topological characterization method in order to obtain a generic description of the attractors while α is varied.

Topological characterization

The main purpose of the topological characterization method is to build a template using topological properties of periodic orbits. The template has been introduced by [START_REF] Birman | Knotted periodic orbits in dynamical systems -I[END_REF][START_REF] Franks | Entropy and knots[END_REF] further to the works of Poincaré [START_REF] Poincaré | Les méthodes nouvelles de la mécanique céleste[END_REF]. According to Ghrist et al. [START_REF] Ghrist | Knots and links in three-dimensional flows Books[END_REF], a template is a compact branched two-manifold with boundary and smooth expansive semiflow built locally from two type of charts: joining and splitting. Each charts carries a semiflow, endowing the template with an expanding semiflow, and the gluing maps between charts must respect the semiflow and act linearly on the edges. This topological characterization method is detailed by Gilmore & Lefranc [START_REF] Gilmore | Topological analysis of chaotic dynamical systems[END_REF][START_REF] Gilmore | The Topology of Chaos[END_REF]. Recently, we detail additional conventions to obtain templates that can be compared and sorted [START_REF] Letellier | Systematic template extraction from chaotic attractors: I-Genusone attractors with an inversion symmetry[END_REF][START_REF] Letellier | Systematic template extraction from chaotic attractors: II. Genus-one attractors with multiple unimodal folding mechanisms[END_REF]. We start with a brief description of the method. As the trajectories are chaotic, they are unpredictable in a long term behavior. But attractors have a time invariant global structure where its orbits compose its skeleton. The purpose of the method is to use the topological properties of these orbits to describe the structure of the attractor. We provide a sum up of this method with eight steps (including our conventions):

(i) Display the attractor with a clockwise flow;

(ii) Find the bounding torus;

(iii) Build a Poincaré section;

The first step permits to ensure that the study will be carried out to the respect of conventions: clockwise flow having a clockwise toroidal boundary and described by a clockwise template. This clockwise convention ensures us to describe template with a unique linking matrix, a keystone to work only with linking matrices [START_REF] Letellier | Systematic template extraction from chaotic attractors: I-Genusone attractors with an inversion symmetry[END_REF]. The toroidal boundary give a global structure that permit to classify attractors. For a given toroidal boundary, a typical Poincaré section is associated according to the Tsankov & Gilmore theory [START_REF] Tsankov | Strange attractors are classified by bounding tori[END_REF]. This Poincaré section contains one or more components to permit an effective discretization of trajectories and consequently an efficient partition of the attractor.

(iv) Compute the first return map and define a symbolic dynamic; (v) Extract and encode periodic orbits;

(vi) Compute numerically the linking numbers between couple of orbits;

The first return map details how two consecutive crossings of a trajectory through the Poincaré section are related and permits to associate a symbol to each point. It permits to define a partition of the attractor and a symbolic dynamic. Associated symbols depend on the parity of the slope (even for the increasing one, odd for the other). Up to this point, periodic orbits structuring the attractor are extracted and encoded using this symbolic dynamic. The linking number between a pair of orbits is a topological invariant indicating how orbits are wind one around another. In this paper, we use the orientation convention of Fig. 2. The linking numbers are equal to the half-sum of the oriented crossings in a regular plane projection. These crossings are obtained numerically by following the path of one orbit and counting the crosses with the other orbit. This algorithm is implemented by splitting orbits into continuous part having either positive or negative derivative of the x value. This permits to simplify the sign identification of the crossing knowing the orientation of orbits.

Convention Permutations Torsions

+1

-1 positive negative positive negative The final steps concern the template:

(vii) Propose a template;

(viii) Validate the template with the theoretical computation of linking numbers.

The template is clockwise oriented. The template of an attractor bounded by genus one torus is defined by a unique linking matrix. This matrix describes how branches are torn and permuted. We use the Melvin & Tufillaro [START_REF] Melvin | Templates and framed braids[END_REF] standard insertion convention: when the branches stretch and squeeze, the left to the right order of the branches corresponds to the bottom to top order. The diagonal elements of the matrix indicate the torsions of branches and off diagonal elements give the permutations between two branches. Finally, to validate a template, we use a procedure introduced by Le Sceller et al. [START_REF] Sceller | Algebraic evaluation of linking numbers of unstable periodic orbits in chaotic attractors[END_REF] that permits to compute linking numbers theoretically from a linking matrix. Linking numbers obtained theoretically with this method have to correspond with those obtained numerically at the step (vi) to validate the template. The challenge of this procedure resides in the step (vii) because it is non trivial to find a template whose theoretical linking numbers correspond to the numerically computed linking numbers.

Attractor A

In this section we will detail the previously described procedure step by step for attractor A.

(i) Display the attractor with a clockwise flow;

We propose to make a rotation of the attractor around the y-axis. Displaying the attractor in the phase space (-x, y), the flow evolves clockwise (Fig. 3a).

(ii) Find the bounding torus;

The attractor is bounded by a genus one torus: a surface with only one hole. Consequently, a Poincaré section with one-component is required.

(iii) Build a Poincaré section;

We build our Poincaré section using x -(2)

P ≡ (y n , -z n ) ∈ R 2 -x n = -x -. (5)
This Poincaré section is a half-plan transversal to the flow as illustrated in grey Fig. 3a. (iv) Compute the first return map and define a symbolic dynamic;

To compute the first return map, we first normalize the intersection value of the flow through the Poincaré section: ρ n . This value is oriented from the inside to the outside (Fig. 3a). Then the first return map is obtained by plotting ρ n+1 versus ρ n (Fig. 3b). This return map is the classical unimodal map made of an increasing branch followed by a decreasing one. This first return map indicates that the classical "horseshoe" mechanism generate this chaotic attractor. The symbolic dynamic is defined as follow: "0" for the increasing branch and "1" for the decreasing one.

(v) Extract and encode periodic orbits;

Using, the first return map, we can locate periodic orbits that the flow visits while it covers the attractor. For instance, there is only one period one orbit because the bisector crosses the map once. We extract five orbits with a period lower than six: 1, 10, 1011, 10110 and 10111 (Fig. 4).

- (vi) Compute numerically the linking numbers;

We compute linking numbers between each pair of periodic orbits. The linking number between a pair of orbits is obtained numerically (Tab. 1).

Table 1. Linking numbers between pairs of orbits extracted from the chaotic attractor A.

(1) (10) (1011) (10110)

(10) -1 (1011) -2 -3 (10110) -2 -4 -8 (10111) -2 -4 -8 -10
(vii) Propose a template;

Using these linking numbers and the first return map structure, we propose the template (Fig. 5). This template is made of trivial part with a chaotic mechanism. The latter is composed by a splitting chart that separates continuously the flow into two branches. The left one encoded "0" permutes negatively over the right one encoded "1"; the latter have a negative torsion. After the torsion and permutation, branches stretch and squeeze to a branch line using the standard insertion convention. This template is thus described by the linking matrix

T (A) = 0 -1 -1 -1 ⟧ . (6)
In literature, two conventions are used to describe templates. One needs a matrix to describe torsions and permutations with an additional array to order the branches for the squeezing mechanism. The other only requires one matrix with respect to the Melvin & Tufillaro [START_REF] Melvin | Templates and framed braids[END_REF] standard insertion convention. For more details about the links between these two conventions, see [START_REF] Gilmore | Algorithms for concatenating templates Chaos[END_REF]. Here, the double line in the right bracket refers to the standard insertion convention indicating that the squeezing mechanism is included in the linking matrix. This use permits to simplify computation between algebraic description of templates using only one matrix. Here, with respect to this convention, the branch "0" is on the right side of the branch "1" after the permutation and torsion, thus branch "0" is over the branch "1" when the branches squeeze. (viii) Validate the template computing theoretically the linking numbers.

The theoretical calculus using the linking matrix and the orbits permits to obtain the same table of linking numbers. This validates the template of A defined by T (A) (Fig. 5).

Attractors B to G

In this section, we will only give the key steps for others attractors: B, C, D, E 1 , E 2 , F and G. We start with the Fig. 6 displaying these attractors for parameters (4) and with a clockwise flow evolution. We can observe that these attractors are made of fully developed chaos (Fig. 6acg) or banded chaos (Fig. 6bdf) or coexisting attractors of banded chaos (Fig. 6e). All these attractors are bounded by a genus one torus. We use the Poincaré sections: [START_REF] Letellier | Unstable periodic orbits and templates of the Rössler system : toward a systematic topological characterization[END_REF] for A, C, G, (8) for D, E 1 , E 2 , F and (9) for B:

P ≡ (y n , -z n ) ∈ R 2 -x n = -x -, (7)
P ≡ (y n , -z n ) ∈ R 2 -x n = -x -, -ẋn < 0, y < -7 , (8)
P ≡ (y n , -z n ) ∈ R 2 -x n = -x -, -ẋn < 0, y < -9 . (9)
We compute the first return maps to these Poincaré sections using a normalized variable ρ n oriented from the inside to the outside of the boundary. The eight return maps (Fig. 7) are multimodal with differential points and the number of their branches are two, three or four. We choose a symbolic dynamic for each first return map with respect to the slope orientation of branches

A 0 1 E 1 0 1 B 0 1 E 2 1 2 C 0 1 2 F 0 1 2 3 D 0 1 2 G 0 1 2 . ( 10
)
We extract a set of orbits for each attractors and numerically compute linking numbers between pairs of orbits (Tab. 2). Thus, we propose templates. They are validated using Le Sceller et al. procedure [START_REF] Sceller | Algebraic evaluation of linking numbers of unstable periodic orbits in chaotic attractors[END_REF]. All the results are summed up in the Tab. 2.

Only attractors C and G have the same template even if they have not the same orbits of period lower than five. The dynamic is not fully developed on each branch: some orbits are missing in both attractors compare to the full set of orbits. On the other hand, the two coexisting attractors have the same linking numbers but the orbits are encoded in two distinct ways. Their linking matrix are related as their template too. In fact, the orientation of the chaotic mechanism is the opposite one. This case permits to underline how orientation conventions are necessary to distinguish these two attractors.

Table 2. Linking numbers between pairs of orbits extracted from the attractors of Fig. 6 and their associated linking matrix describing their template.

Linking numbers

Linking matrix

A (1) (10) (1011) (10110) (10) -1 (1011) -2 -3 (10110) -2 -4 -8 (10111) -2 -4 -8 - 10
T (A) = 0 -1 -1 -1 ⟧ B (1) (10)
T (C) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 -1 -1 -1 -1 -1 -1 -1 0 D (1) (221) (211) (2221) (2211) (221) -4 (211) -4 -12 (2221) -5 -15 -15 (2211) -5 -15 -15 -20 (22110) -6 -18 -18
-24 -24

T (D) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -4 -4 -4 -4 -3 -3 -4 -3 -2 E 1 (1) (10) (1011) (101111)
T (E 1 ) = -2 -3 -3 -3 ⟧ E 2 (1)
T (E 2 ) = -3 -3 -3 -2 ⟧ F (3)
T (F) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -4 -4 -4 -4 -4 -3 -3 -3 -4 -3 -2 -3 -4 -3 -3 -3 G (1) (21) (2111) (2021) (21) -1 (2111) -2 -3 (2021) -2 -3 -6 (2110) -2 -3 -6 -6 T (G) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 -1 -1 -1 -1 -1 -1 -1 0 4. Subtemplate

Algebraical relation between linking matrix

In this section, we will use some algebraic relations between linking matrices already defined in our previous papers [START_REF] Letellier | Systematic template extraction from chaotic attractors: I-Genusone attractors with an inversion symmetry[END_REF][START_REF] Letellier | Systematic template extraction from chaotic attractors: II. Genus-one attractors with multiple unimodal folding mechanisms[END_REF]. Here, we provide an overview of these relations.

In the following description, a strip also nominates a branch of a branched manifold and a writhe is the number of positive crossings minus the number of negative crossings of an oriented link diagram. First of all, a linker is a synthesis of the relative organization of n strips: torsions and permutations in a planar projection (Fig. 2). A mixer is a linker ended by a joining chart that stretch and squeeze strips to a branch line. In the previous section, templates are only composed by one mixer defined by a linking matrix. We also define the concatenation of a torsion with a mixer and the concatenation of two mixers using the operator "+" in the following equations (see [START_REF] Letellier | Systematic template extraction from chaotic attractors: I-Genusone attractors with an inversion symmetry[END_REF][START_REF] Letellier | Systematic template extraction from chaotic attractors: II. Genus-one attractors with multiple unimodal folding mechanisms[END_REF] for more details). The concatenation of a torsion with a torsion is a torsion; the concatenation of a mixer with a torsion is a mixer and the concatenation of a mixer with a mixer is a mixer. We remind that X designates the attractor, T X its template and T (X ) the linking matrix that defines its template. Using these algebraical relations between mixers and torsions, we can link the mixers of previously studied attractors:

T (C) = T (G) T (E 1 ) = T (E 2 ) p -3 -3 -3 -2 ⟧ = -2 -3 -3 -3 ⟧ p T (B) = [-5] + T (A) -6 -6 -6 -5 ⟧ = -5 -5 -5 -5 + 0 -1 -1 -1 ⟧ T (E 2 ) = [-2] + T (A) -2 -3 -3 -3 ⟧ = -2 -2 -2 -2 + 0 -1 -1 -1 ⟧ (11)
T (E 1 ) = T (E 2 ) p means that it is the transposed matrix: thus mixer of E 2 is symmetric to the mixer of E 1 in a sens that the branch order is reversed. Then we obtain T (B) = [-5] + T (A) because the concatenation of an odd torsion before a mixer reverse the branch order. Finally, T (E 2 ) = [-2] + T (A) because the concatenation of an even torsion before a mixer do not reverse the branch order of the mixer.

Subtemplates

A subtemplate is defined as follow by Ghrist et al. [START_REF] Ghrist | Knots and links in three-dimensional flows Books[END_REF]: a subtemplate S of a template T , written S ⊂ T , is a topological subset of T which, equipped with the restriction of a semiflow of T to S, satisfies the definition of a template. For the eight attractors previously studied (A, B, C, D, E 1 , E 2 , F and G) we will demonstrate that their templates are subtemplates of the template of C made of one mixer defined by

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 -1 -1 -1 -1 -1 -1 -1 0 . ( 12
)
Using the linking matrix defining T A and T C , we directly find that T (A) is a subset of T (C) with the two first lines and columns

0 -1 -1 -1 ⟧ ⊂ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 -1 -1 -1 -1 -1 -1 -1 0 . ( 13
) 0 1 2 (a) T C (b) T A ⊂ T C Figure 8. Template of A is a subtemplate of the template of C.
The strip organization of T A are the same of the two first strips of T C . This means that T A is a subtemplate of T C : T A ⊂ T C . This is illustrated on Fig. 8 where we only display the mixers and not the trivial part of the template that link the bottom to top on the left side to have a clockwise flow as shown Fig. 5; this is also the case for the remainder of this article. We will use graphical representation of the templates and subtemplates because it details the relation between template and subtemplate. This representation combined with algebraical relations between linking matrices proves that a template is a subtemplate of a template. Attractors B,E 1 and E 2 display banded chaos because they are composed by several strips, or bands, with writhes. We start with the template T B . We know that this template can be considered with five negative torsions before a "horseshoe" mechanism [START_REF] Poincaré | Les méthodes nouvelles de la mécanique céleste[END_REF]. Thus we have to find a subtemplate that goes through a "horseshoe" mechanism and have five negative torsions. Letellier et al. [START_REF] Letellier | Topological invariants in perioddoubling cascades[END_REF] underline the fact that if a writhe is observed in an attractor bounded by a genus one torus, this is equivalent to two torsions by isotopy; the sign of the torsions is the same of the writhe one (see FIG. 3 and FIG. 4 of [START_REF] Letellier | Topological invariants in perioddoubling cascades[END_REF] for additional details). As a consequence, a subtemplate with n portions induces n writhes that are equivalent to 2n torsions by isotopy; the sign of these 2n torsions is the sign of the writhe that is also the sign of the permutations of subtemplate portions. We propose to build such a subtemplate where Fig. 9a is T B as the subtemplate of T C .

Banded chaos

We propose to use algebraical relations between linking matrices to validate subtemplates. Fig. 9b is the template of B shaped as a subtemplate of C. When we establish the template of B, we chose to use the Poincaré section (9) that corresponds to the portion far from the inside of the attractor; this portion is labeled A on Fig. 9b. We propose to concatenate its components with respect to their relative order: from A to B, then B to C and C to A that are respectively: a mixer, a strip without torsion and a negative torsion. We finally consider the transformation by isotopy that does not have an impact on the orientation of the strips. Thus we decide to concatenate the 2n torsions after the concatenation of the components. As illustrated Fig. 9b, there is two negative permutations (B to C over A to B and C to A over A to B) that are equivalent to 2 × 2 = 4 negative torsions (dashed circles of Fig. 9b). Consequently, we concatenate all these mixers and torsions

-1 -1 -1 0 ⟧ + [-1] + [-2] + [-2] = -6 -6 -6 -5 ⟧ = T (B) (14)
and obtain the linking matrix defining the template of B. Consequently T B is a subtemplate of T C . We now consider the two coexisting attractors E 1 and E 2 . They have a similar structure and coexist in the phase space for distinct initial conditions. The mixer of T C with three branches has also a symmetric structure: the middle of the second strip of the mixer T (C) is a reflecting symmetry axis where the left side is symmetric of the right side. To build the T E1 and T E2 as subtemplate of T C , we take this symmetry into account.

We propose these subtemplates as illustrated Fig. 10. We compute the concatenation of torsions and mixer of these figures. For T E1 , we have two parts: one is a mixer and the other is a strip without torsion. These parts permute negatively once in a writhe; it algebraically corresponds to a concatenation of two negative torsions. Thus, the linking matrix of such a subtemplate (Fig. 10a) is Similarly, the linking matrix of T E2 as a subtemplate (Fig. 10b) is

0 -1 -1 -1 ⟧ + [-2] = -2 -3 -3 -3 ⟧ = T (E 1 ) . ( 15
) (a) T E1 ⊂ T C (b) T E2 ⊂ T C (c) T E1 and T E2 coexisting
-1 -1 -1 0 ⟧ + [-2] = -3 -3 -3 -2 ⟧ = T (E 2 ) . (16)
These algebraical relations between template and subtemplate linking matrices with Fig. 10 prove that T E1 ⊂ T C and T E2 ⊂ T C . Moreover, these two subtemplates are symmetric one to the other by reflection and coexist in the template of C (Fig. 10c).

Concatenation of mixers

We now consider T F , the template of F, made of four strips. In a previous paper [START_REF] Letellier | Systematic template extraction from chaotic attractors: II. Genus-one attractors with multiple unimodal folding mechanisms[END_REF], we demonstrate that the concatenation of two mixers is a mixer and its number of strips is the product of the number of strips of each mixer. Thus, the concatenation of two mixers made of two branches is a mixer with four branches. Our hypothesis is that the four strips of T F are the result of this process. To validate it, we draw its subtemplate by splitting the template of C into the two symmetric part containing a mixer; we obtain the Fig. 11a. We decompose this subtemplate (Fig. 11b) in parts: the two parts contain a mixer and these parts permute negatively once. Thus, we concatenate a mixer before a mixer and two negative torsions (cf. 4.2.1). The first concatenation gives a mixer defined by a linking matrix; the algebraical relation necessary to obtain this matrix are detailed in [START_REF] Letellier | Systematic template extraction from chaotic attractors: II. Genus-one attractors with multiple unimodal folding mechanisms[END_REF]. The linking matrix of T F as a subtemplate (Fig. 11a) is:

-1 -1 -1 0 ⟧ + 0 -1 -1 -1 ⟧ + [-2] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 0 -1 -1 0 0 + 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 + -1 -1 -1 -1 -1 0 0 -1 -1 0 0 -1 -1 -1 -1 -1 + [-2] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -2 -2 -2 -2 -2 -1 -1 -1 -2 -1 0 -1 -2 -1 -1 -1 + [-2] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -4 -4 -4 -4 -4 -3 -3 -3 -4 -3 -2 -3 -4 -3 -3 -3 = T (F ) (17)
This algebraical relation between template and subtemplate linking matrices associated to Fig. 11 prove that T F ⊂ T C . We now consider the attractor D, we have T D ⊂ T F directly from their mixers

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -4 -4 -4 -4 -3 -3 -4 -3 -2 ⊂ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -4 -4 -4 -4 -4 -3 -3 -3 -4 -3 -2 -3 -4 -3 -3 -3 ; (18)
this is illustrated on Fig. 12. We previously obtain that T F ⊂ T C , thus we prove that T D ⊂ T C . Consequently, we prove that the six templates of attractors A, B, D, E 1 , E 2 and F are subtemplate of the template of the attractor C; it is a template with six subtemplates. We remind that G and C have the same template (T G = T C ).

A template for the whole bifurcation diagram

We obtain a unique template containing the eight templates of attractors. We now consider the whole bifurcation diagram (Fig. 1) and not only specific attractors. In this section, we will show that the template of C contains all attractors templates for any parameter value take from this bifurcation diagram. We use the Poincaré section (3) and build return maps using y n for α ∈] -2; 1.8[ when an attractor is solution. We associate a symbolic dynamic with the three symbols "0", "1", "2" of T C . Note that Barrio et al. [START_REF] Barrio | Topological changes in periodicity Hubs of dissipativity systems[END_REF] also use this process to study return maps of a Rössler system from a Lyapunov diagram. The authors display return maps with superstability curve and coexisting stable points. Here we prefer to collect extrema points to make a partition of the bifurcation diagram. 3). This partition give a symbolic dynamic with three symbols "0", "1", "2" depending on α.

In the diagram Fig. 13, we indicate the values of y n splitting the return maps into two or three parts. We remind the reader that we orientate application from the inside to the outside of the attractor. Thus, the branches are labelled with symbol number increasing while y n decrease. Fig. 13 reveals that the separator values are linear to α. We note y 0 1 (α) the value of y n that split branches "0" and "1" and y 1 2 (α) the value of y n that split branches "1" and "2". A linear regression gives y 0 1 (α) = 1.43638α -6.76016 y 1 2 (α) = 2.18237α -11.1289 .

Up to this point, if there is an attractor solution, its orbits can be encoded with the symbols depending on the previous equations. This also requires the use of the Poincaré section (3). For a given range of a bifurcation parameters (α ∈] -2; 1.8[), the parameters of the Rössler system depends on α: a(α), b(α) and c(α). The fixed points depend on the parameters and the Poincaré section depend on the fixed points. Thus, we obtain a Poincaré section and its partition [START_REF] Gilmore | Algorithms for concatenating templates Chaos[END_REF] depending on α while the template is defined by the linking matrix

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 -1 -1 -1 -1 -1 -1 -1 0 . ( 20
)
The main result is that the topological characterization of chaotic attractors can be extended as a description of various attractors whose parameters come from one bifurcation diagram. In this bifurcation diagram, we show that there are regimes where the chaotic mechanisms are topologically equivalent (T C = T G ), symmetric (T E1 and T E2 ) and they are a subset of the same chaotic mechanism. The point is that our work can help to understand the complex structure of attractors considering them as subtemplates of their neighbors (in term of bifurcation parameter). This also enlarge the possibility to use the topological characterization to describe more than an attractor, but an entire bifurcation diagram.

Conclusion

In this paper we study eight attractors of the Rössler system. The parameters values of these attractors come from a bifurcation diagram that exhibits various dynamics such as coexisting attractors. For each attractor we apply the topological characterization method that give us a template of the attractor. These templates detail the chaotic mechanism and these are only made of stretching and folding mechanism followed by a squeezing mechanism with two, three and four strips.

The second part of this paper is dedicated to the proof that the eights templates are subtemplates of a unique template: the template of C. The main result here is that a template is no longer a tool to describe one attractor but also a set of neighbours attractors (in the parameter space). Thus for a bifurcation diagram, we build a partition using symbols of T C . This partition over the whole diagram gives a global structure of attractors for a range of parameters.

This better understanding of the structure of bifurcation diagram can help researchers that want to explore the behaviour of their system, especially if it exhibits chaotic properties. For instance Matthey et al. [START_REF] Matthey | Experimental study of limit cycle and chaotic controllers for the locomotion of centipede robots, Intelligent Robots and Systems[END_REF] design a robot using coupled Rössler oscillators to simulate its locomotion. A similar theoretical analysis of their system can provide various set of parameters with specific chaotic properties that might induce new locomotion pattern. A partitioned bifurcation diagram details the various non equivalent dynamical behavior of the system to find them.

This work on templates of attractors from a unique bifurcation diagram is a first step that can lead to a description of manifolds using templates. It is also a new way to classify templates grouping them as subtemplates and not claiming that there exist six new attractors for the Rössler system. To conclude, this work permits to apply the topological characterization method to a set of attractors from a bifurcation diagram. The partition of a bifurcation diagram associated with a unique template is a new tool to describe the global dynamical properties of a system while a parameter is varied. In future works, we will apply this method on attractors bounded by higher genus torus to highlight how symmetry breaking or template number of branches are related in a bifurcation diagram.

  of[START_REF] Sprott | Asymmetric bistability in the Rössler system International Journal of Bifurcation and Chaos submitted[END_REF]: we use y n and they use M which is a local maximum of x. Consequently, in both case, values close to zero correspond to value close to the center of the attractor and oppositely, high absolute values correspond to the outside boundary of the attractor.

Figure 1 .

 1 Figure 1. Bifurcation diagram when α varies: α increasing (red) and α decreasing (black). This figure is similar to the last part of the FIG. 5 of [5]. Using α increasing or decreasing replaces different initial conditions used in the original paper [5]. A, B, C, D, E 1 , E 2 , F and G refers to attractors solutions with the parameters indicated in (4).

Figure 2 .

 2 Figure 2. Convention representation of oriented crossings. The permutation between two branches is positive if the crossing generated is equal to +1, otherwise it is a negative permutation. We use the same convention for torsions.

Figure 3 .

 3 Figure 3. (a) Chaotic attractor A solution to the equation (1). Parameter value α = -0.25 with the initial conditions x = -1.25, y = -0.72 and z = -0.1. (b) First return map to the Poincaré section (5) using ρn (the arrow indicates the orientation).

Figure 4 .

 4 Figure 4. Orbits of the chaotic attractor A.

0 1 Figure 5 .

 15 Figure 5. T A : template of attractor A.

G for α = 1 . 22 Figure 6 .

 1226 Figure 6. Eight chaotic attractors for different values of α from the bifurcation diagram (Fig. 1)

  A for α = -0.25 (b) B for α = 0.5 (c) C for α = 0.78 D for α = 0.86 (e 1 ) E 1 for α = 1 (e 2 ) E 2 for α = 1 F for α = 1.135 (g) G for α = 1.22

Figure 7 .

 7 Figure 7. First return maps of the eight attractors of the Fig. 6.

Figure 9 .

 9 Figure 9. (a) Template of attractor B is a subtemplate of the template of C. (b) Template of attractor B shaped as a subtemplate of C. The arrow over A refers to the Poincaré section (9) with the ρn orientation.

Figure 10 .

 10 Figure 10. Coexisting templates T E 1 and T E 2 are subtemplates of the template of C.

Figure 11 .

 11 Figure 11. (a) The template of F is a subtemplate of the template of C. (b) T F with two mixers.

Figure 12 .

 12 Figure 12. Template of D is a subtemplate of the template of F and consequently, it is also a subtemlate of the template of C.

Figure 13 .

 13 Figure 13. Partition of the bifurcation diagram when α varies build using first return maps on yn of the Poincaré section (3). This partition give a symbolic dynamic with three symbols "0", "1", "2" depending on α.


Marvin Faix
email: marvin.faix@inria.fr.

Jorge Lobo
email: jlobo@isr.uc.pt.

Raphael Laurent
email: raphael.laurent@probayes.com.

Dominique Vaufreydaz
email: dominique.vaufreydaz@inria.fr.

Emmanuel Mazer
email: emmanuel.mazer@inria.fr.

Stochastic Bayesian Computation for Autonomous Robot Sensorimotor System

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Bayesian approaches have been successfully applied to autonomous robots when they must operate in environments not specially designed for them, where they must handle uncertainty and decide with incomplete data [START_REF] Thrun | Probabilistic Robotics[END_REF][START_REF] Bessière | Probabilistic reasoning and decision making in sensory-motor systems[END_REF][START_REF] João | Probabilistic approaches to robotic perception[END_REF]. One of the drawbacks of such approaches can be the demands on the computational power, and insights on this can be found in biology. We want to go beyond the standard current paradigm of exact hardware running software, and explore approximate computing hardware that can support the required computations with enough accuracy.

This work was developed in the scope of the European BAMBI FET project (Bottom-up Approaches to Machines dedicated to Bayesian Inference, www.bambi-fet.eu). BAMBI's goal is to rely on Bayesian theory, and an understanding of computing by living beings, to design a non Von Neumann bio-inspired probabilistic machine. All living systems process information with uncertainty, and this occurs at all levels, from microscopic molecular scale up to complex human perception. An interesting approach to model these processes is Bayesian modelling. This kind of modelling is already used in behaviour prediction and decision making where Bayesian inference proves powerful to carry out such computation [START_REF] Bessière | Probabilistic reasoning and decision making in sensory-motor systems[END_REF]. The study of information transmission in unicellular organisms, namely the Chlamydomonas reinhardtii, allows us to conjecture that probabilities are coded with a kind of binary telegraphic signals. Using a bitstream representation to encode probabilities, we design a bioinspired probabilistic machine based on these principles. For an initial version we rely on clocked bit streams to build a machine that takes into account the uncertainty of its inputs and can performs the computations with a small low power circuit. One key benefit is that due to the stochastic representation, the circuit is tolerant to sporadic logic gate faults. This has been increasingly relevant, as the limits of Moore's law are ever imposing, with faults occurring when power is lowered or clock speeds increased in current evermore minute transistors.

In this paper we present our compilation tool chain for a first Bayesian machine based on coding probabilities with stochastic bitstreams. This allows the automatic generation of the circuit implementation. The machine computes an exact inference with approximate computation. A fast high level simulation platform was developed, and the full systems was tested on the robotic classic problem of obstacle avoidance.

In the next two sections we present inspiring works from the literature, and we clarify what is, for us, a probabilistic machine. In section IV we describe our compilation toolchain, that starts from the Bayesian model description, in a high level language programming, and provides the final circuit in VHDL (a hardware description language from which circuits can be synthesised). The following section presents the sensor fusion application as well as simulation and tests on the robot, followed by conclusions.

Related Work

Many applications in Robotics already use probabilistic inference and Bayesian theory [START_REF] Thrun | Probabilistic Robotics[END_REF][START_REF] Bessière | Probabilistic reasoning and decision making in sensory-motor systems[END_REF][START_REF] João | Probabilistic approaches to robotic perception[END_REF]. In this work we use a specific Programming language named ProBT [START_REF] Kamel Mekhnacha | Exact and approximate inference in ProBT[END_REF][START_REF] Bessière | Bayesian programming[END_REF] but many languages dedicated to probabilistic model exist Figaro [START_REF] Pfeffer | Practical probabilistic programming[END_REF], Blog [START_REF] Milch | BLOG: Probabilistic models with unknown objects[END_REF]. In the same way as BAMBI project interests, other teams try to perform probabilistic inference in hardware [START_REF] Vigoda | Analog logic: Continuous-Time analog circuits for statistical signal processing[END_REF][START_REF] Kumar | Natively probabilistic computation[END_REF][START_REF] Michael | Stochastic architectures for probabilistic computation[END_REF]].

Vigoda's architectures [START_REF] Vigoda | Analog logic: Continuous-Time analog circuits for statistical signal processing[END_REF] use message passing algorithm [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF] to propagate uncertainty and compute inferences. Close to the temporal coding of probabilities, Vigoda represents the probabilities with analog values. These architectures are used to perform exact inference. At the opposite end Mansinghka tries to perform Bayesian computation with approximate inference using a sampling method [START_REF] Kumar | Natively probabilistic computation[END_REF]. Jonas [START_REF] Michael | Stochastic architectures for probabilistic computation[END_REF] followed this idea and designed samplers based on Markov Chain Monte Carlo algorithms to represent probability distributions.

In our approach, we are focused on the design of a probabilistic machine which computes exact inference with approximate computation, it means less precision. This is due to the stochastic arithmetic used, based on stochastic time coding to represent probabilities. The previous work of Von Neumann [START_REF] Von | Probabilistic logics and the synthesis of reliable organisms from unreliable components[END_REF] in 1956 and Gaines [START_REF] Gaines | Stochastic computing systems[END_REF] showed the interest of a such coding. Indeed the simplicity of basic operators is the main strength of this approach. However the computation time and the low accuracy have not allowed this way to compute to cope with the emergence of faster and better processors as the number of transistors on a chip grew at an astounding rate. The new interest of the community in these works is due to the current constraints with higher numbers of transistors on a chip, where fault tolerance, low power solutions, and robustness to noise are now important factors (see for example [START_REF] Qian | An architecture for fault-tolerant computation with stochastic logic[END_REF][START_REF] Li | A stochastic reconfigurable architecture for fault-tolerant computation with sequential logic[END_REF]).

The TrueNorth project [START_REF] Steven K Esser | Cognitive computing systems: Algorithms and applications for networks of neurosynaptic cores[END_REF] from IBM is also interested in designing a non Von Neumann bio-inspired architecture. It proposes a neuromorphic system for neuronal networks and uses fixed-point arithmetic units to compute the output of a neuron knowing its inputs. In this point we differ with this approach because we think that the use of probabilities is a natural way to handle the uncertainty which shapes the world [START_REF] Edwin | Probability theory: the logic of science[END_REF].

3 What is a Bayesian Machine ?

In this part we specify what is a probabilistic machine dedicated to Bayesian inference in our context. We are focused on soft evidences (as opposed to hard evidences) which are probability distributions over an evidence. A Bayesian machine is a machine which resolves an inference problem. This machine takes soft evidences as inputs and outputs a probability distribution over the searched variables.

Soft Evidences

In the context of Bayesian model hard evidence could be defined as a deterministic observation whereas soft evidence represents uncertainty over known variables. In this case a soft evidence is represented by a probability distribution over the corresponding evidence variable. In our robotic sensor fusion application, soft evidences describe the precision of the sensor reading.

Structure of A Bayesian Machine

A Bayesian machine takes soft evidences as inputs. The outputs are probability distributions over searched variables. With this representation we ensure homogeneity between inputs and outputs; both are probability distributions. The third set of signals are model parameters of the joint distribution which are constant probability values. Each probability value is coded as a bitstream. Because bitstreams are drawn from a Bernoulli sequence, we use biased PRNG (Pseudo Random Number Generator) to generate the constant parameters of the model, a fast biased TRNG (True Random Number Generator) [START_REF] Cherkaoui | A very high speed true random number generator with entropy assessment[END_REF] could be used to have the randomness in the circuit.

Inside the probabilistic machine there are two different computation blocks which correspond to the sums and products to be computed for the exact inference. Figure 1 shows the structure of a probabilistic machine for a very simple example. We infer over variable M knowing two soft evidences P (D 1 ) and P (D 2 ), P (M ) is the prior of M , P (D 1 ) and P (D 2 ) are the constant parameters of the model.

Biased TRNGs

The Compilation ToolChain

In this part we explicit the compilation tool chain which starts from a Bayesian model described in a Bayesian programming language, ProBT [START_REF] Bessière | Bayesian programming[END_REF], and automatically designs the probabilistic machine which implements the inference computation over the model described. The electronic design is a VHDL file which can be used both with our simulation platform or to implement the circuit on a FPGA (Field Programmable Gate Array, a reconfigurable logic device).

Bayesian Program

A Bayesian machine can be defined using properties, notation and formalism of Bayesian programming [START_REF] Bessière | Bayesian programming[END_REF]. To design our machine we need to define the Bayesian model, the soft evidence inputs, the constant parameters of the model and the inference to compute. A Bayesian program is a program which describes the probabilistic link between variables. A general model is a joint probability distribution over a set of discrete and finite variables: P (S ∧ K ∧ F ) where S, K and F are respectively Searched, Known and Free variables of our model. Each variable could be a conjunction of variables, for example K = K 1 ∧ ... ∧ K i . The inputs of the machine are soft evidences, so we define the soft evidences over the Known variables K i , P (K i ) which is a probability distribution over K i . To do an inference, as for instance computing P (S| P (K 1 ) ∧ ... P (K i )), is equal to compute the following expression: P (S| P (K 1 ) ∧ ...

∧ P (K i )) = 1 Z K1 P (K 1 ) . . . Ki P (K i ) F P (S ∧ K ∧ F ),
where Z is the normalization constant:

Z = S   K 1 P (K1) . . . K i P (Ki) F P (S ∧ K ∧ F )   . (1)
ProBT generates a computation tree which is simplified thanks to the Successive Reduction Algorithm (SRA) [START_REF] Kamel Mekhnacha | Exact and approximate inference in ProBT[END_REF], itself similar to the sum-product algorithm [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF]. The simplification is possible thanks to the dependencies description between variables in the model defined in ProBT. The simplified computation tree is given as input to the compiler which will design automatically the machine.

Compilation

The compiler starts with the Bayesian model description written in ProBT and it generates the corresponding VHDL file which describes the dedicated probabilistic computer. In exact inference the expression to compute is composed of many sums and products. This is the fundamental point of using the temporal coding of probabilities. Sum and product operators are easy to implement and need few components using this kind of coding. In this case the terms stochastic sum and stochastic product are used. Every probability in our model is coded with a bitstream. Assume that two bitstreams BS 1 and BS 2 , encoding the probability values P 1 and P 2 , are decorrelated, the result of these two bitstreams in inputs of an AND gate is a bitstream BS p which encodes the product P 1 × P 2 . In the same way, just using an OR gate extended with a memory we can sum probabilities. These operation are non-normalised. Normalisation can be implemented thanks to a JK flip flop. A schema of these components is given in figure 2. More information about these components and their utilisation is available in [START_REF] Faix | Cognitive computation: A bayesian machine case study[END_REF].

Figure 2: The three components for product, sum and division in hardware using bitstream representation

First Bayesian Machines

We have previously validated our approach on two examples, see [START_REF] Faix | Cognitive computation: A bayesian machine case study[END_REF]. The first one is a very simple example to explicit the corresponding VHDL file from the computational tree given by ProBT. Let a joint probability distribution P (M ∧ D 1 ∧ D 2 ) over three variables, the searched variable of the Model M and two known Data D 1 and D 2 . The inference to compute the expression P (M | P (D 1 ) P (D 2 )):

P (M ) = 1 Z P (M )( D 1 P (D1)P (D1|M ))( D 2 P (D2)P (D2|M )), (2)
where P (D 1 ) and P (D 2 ) are the soft evidences. By nature of Bayesian inference, the circuit is duplicated several times. Indeed only the inputs and constants model signals change. Figure 3 shows a subblock of RTL (Register Transfer Level) description generated synthesising the VHDL code to compute a subpart of the expression 2 which is D1 P (D 1 )P (D 1 |m). The second example is a Bayesian filter as an application in telecommunication. The goal was to synchronise two Linear Feedback Shift Registers (LFSRs). The first LFSR emits bits corresponding to its state and the second one receives these bits through a noisy transmission chain. A mathematical study shows that inference of the emitter state is done by computing the following expression:

P (S T = s| PT (O T )) = 1 Z P (a(s))P Out(s) (O T ) , ( 3
)
where s is the state, PT (O T ) is the soft evidence over the received bit, a(s) is the ancestor state of s (see [START_REF] Faix | Cognitive computation: A bayesian machine case study[END_REF] for more detail). Figure 4 shows the corresponding architecture for the LFSR synchronisation. We have shown in [START_REF] Faix | Cognitive computation: A bayesian machine case study[END_REF] that the LFSR synchronisation circuit works with our architecture with transmission error less than 30%, table 1. The mathematical study of the LFSR synchronisation application drastically reduced the number of components needed to compute the inference. In the following section we present out target robotic application that provides a bigger challenge concerning the size of the Bayesian machine.

An Autonomous Sensorimotor System

We chose a classic robotic application to show our solution on a problem with a bigger number of variables and components. We next present the implemented autonomous sensorimotor system, going over the the objectives of the sensor fusion, defining the sensor model, the actual implementation of such an architecture and the results obtained.

Objectives

The goal is to allow an autonomous robot to adjust its trajectory to avoid obstacles given a distance estimate provided by two types of sensors. The infrared sensors and ultrasonic sensors are fused to increase precision over the distance estimation. The robot is shown in figure 5. The result of this inference is the probability distribution over the rotation velocity enabling the robot to avoid the obstacle while moving forward.

Sensor Model

Each sensor has is own model depending on the distance. We define three levels of distance: close, medium and far. The robot turns by increasing the velocity of the wheel corresponding to the opposite side of the direction wanted. We define five rotation velocities: speed on the left, half speed on the left, null, half speed on the right, speed on the right. We set an uniform distribution over the prior distance variables. Both infrared and ultrasonic sensors have 80% confidence on the reading knowing the distance. We know the rotation velocity knowing the distance model too. A piece of ProBT code which describes this system is given in appendix.

Fusion

We define our joint distribution model:

P (D 0 ∧ D 1 ∧ D 2 ∧ IR 0 ∧ IR 1 ∧ IR 2 ∧ U S 0 ∧ U S 1 U S 2 ∧ V ROT ) = P (V ROT |D 0 D 1 D 2 ) × P (D 0 ) × P (D 1 ) × P (D 2 ) P (IR 0 |D 0 ) × P (IR 1 |D 1 ) × P (IR 2 |D 2 ) P (U S 0 |D 0 ) × P (U S 1 |D 1 ) × P (U S 2 |D 2 ) (4)
This is the conjunction of variables defined in the sensor model. It means conjunction between the three priors variables, the three conditionals variables over IR (infrared) variable, the three conditionals variables over US (ultrasonic) variable and the conditional variable over rotation velocity knowing distances. Inference over this model corresponds to asking P (V ROT | PIR0 PIR1 PIR2 PUS0 PUS1 PUS2 ) and computing the following expression 5:

P (V ROT | PIR0 PIR1 PIR2 PUS0 PUS1 PUS2 ) = D2 P (D 2 ) × ( U S2 P (U S 2 )P (U S 2 |D 2 )) × ( IR2 P (IR 2 )P (IR 2 |D 2 ))× D1 P (D 1 ) × ( U S1 P (U S 1 )P (U S 1 |D 1 )) × ( IR1 P (IR 1 )P (IR 1 |D 1 ))× D0 P (D 0 ) × ( U S0 P (U S 0 )P (U S 0 |D 0 )) × ( IR0 P (IR 0 )P (IR 0 |D 0 )× P (V ROT |D 0 D 1 D 2 ) (5)

Circuit Implementation for Sensorimotor Fusion

The generated circuit involves about 2500 components and several thousands of signals. The circuit is duplicated five times, indeed the searched variable P (V ROT ) has five possible values P (V ROT 0 )...P (V ROT 4 ). So the circuit is fully parallel in terms of the cardinality of the searched variable P (V ROT i ). The important thing to observe is the number of stages in the circuit. The time dilution of probabilities, due to the many required products, decreases drastically the probabilistic value. So the bitstream representing the probability value is filled with a lot of 0. The bitstream size has to be very long to evaluate this value because of our temporal coding choice.

To use the robot platform, a simulation c++ code is generated from the VHDL code. It exactly corresponds to the VHDL code in terms of component and signals.

Simulation and Results

To simulate our circuit we used the c++ code which takes as inputs the sensors reading, computes the probability distributions and send the order for the rotation velocity. The given linear velocity is 0.2ms -1 because of the time computation of our simulator. Indeed for a bitstream size of 5000, the simulator needs 86ms to send an order. This computation time is linear with the bitstream length used to the computation in the simulator. The specification for this robot was 100ms between two orders with this linear velocity. So if we want to increase the linear velocity we need to make simulations faster or decrease the bitstream size. These results obviously show that a bigger bitstream size increases the precision over the rotation velocity probability distribution. The result over the variable V ROT0 is volatile because of the low probability value for this variable. A final validation with the robot avoiding various obstacles was done. In an indoor office scenario, the robot successfully avoided obstacles using the Bayesian machine implemented on the simulated circuit. A video is available at http://www.bambi-fet.eu/a-robot-controller-powered-with-the-bayesian-machine-1/, and image 5 shows some of the frames where you can see the robot and the obstacles put in its path, in this case the corridor walls.

Bit stream size

Tool V ROT0 V ROT1 V ROT2 V ROT3 V ROT4

Conclusions

We presented our work on Bayesian machine using stochastic arithmetic to compute exact Bayesian inference. These non Von Neumann architectures face the classic electronics constraints: fault tolerant circuit design, low power, energy efficient and robust circuit, especially because the probability coding intrinsically deals with these constraints. The temporal coding is one of the main strength of our approach, it relies with living beings model computation studied in BAMBI project. For us, a bio-inspired machine is a machine which takes probability distribution in inputs, computes with an arithmetic dedicated to the probabilistic model, transmits information with stochastic binary signals representing temporal coding of probabilities and outputs probability distribution.

Now we are able to transform any Bayesian program into a Bayesian machine, which computes a specific inference, described in VHDL language. This description can both be simulated with our specific high level simulation platform or synthesized and emulated on a FPGA. Our method has been tested in [START_REF] Faix | Cognitive computation: A bayesian machine case study[END_REF] and the automatic compilation toolchain presented was tested on a classic robotic application which is the sensor fusion enabling the robot to know the distance of an obstacle and avoid it. It is possible to reduce the circuit size by the number of states of the searched variable P (V ROT ) designing multiplexer in inputs of components which need directly a soft input or parameter constant value of the system. In this case multiplexer selects the right soft inputs and parameter constant values of the model for the corresponding P (V ROT i ) computed.

The main limitations of this work are the classic limits in exact Bayesian computation, as the number and dimension of variables increases, problems become intractable. Problems intractable with classic computations on standard processors are still intractable with our method. We are now working on an approximate inference compiler to scale some intractable Bayesian problems in exact inference.

P 1 Z 1 Z 1 ZFigure 1 :

 1111 Figure 1: The probabilistic machine: its inputs are the soft evidences (arrows incoming from the left) and its output is the probability distribution over the variable of interest, in this case M . The probability values of the internal parameters of the terms of the joint distribution are set by a set of TRNGs (True Random Number Generators)

Figure 3 :

 3 Figure 3: RTL (Register Transfer Level) description for the stochastic circuit computing D1 P (D 1 )P (D 1 |m)

Figure 4 :

 4 Figure 4: The architecture of a circuit implementation of a Bayesian filter using stochastic arithmetic for LFSR synchronisation.

Figure 5 :

 5 Figure 5: Robot performing obstacle avoidance using the Bayesian machine.

Table 1 :

 1 Results LFSR synchronisation

	Transmission Error	0%	10%	20%	30%
	Correctly
	Recognized	198	188	184	177
	States
	Iterations
	Before	2	57	59	64
	Synchronization

Table 2 :

 2 Comparison between theoretical result and simulation for different bitstream sizes

	without Bit stream	ProBT	0.005	0.012	0.068	0.293	0.622
	5000	simulator error	0 n.a	0 n.a	0.052 24%	0.21 28%	0.74 19%
	10000	simulator error	0 n.a	0 n.a	0.040 41%	0.320 9%	0.640 2.9%
	100000	simulator 0.004 error 25%	0.007 4%	0.064 6%	0.331 13%	0.594 4.5%
	1000000	simulator 0.003 error 41%	0.013 8%	0.068 0%	0.289 1.5%	0.628 1%

Acknowledgement

This work was made possible thanks to two grants from CNRS: Nano-Bayes and Defi-Bayes, and thanks to the EU collaborative FET Project BAMBI FP7-ICT-2013-C, project number 618024.

project BAMBI -Bottom-up Approaches to Machines dedicated to Bayesian Inference -FP7-ICT-2013-C, project number 618024 (www.bambifet.eu).


Stéphane Théveny

Frédéric André

Yves Elskens

On frequency and time domain models of traveling wave tubes

Keywords: traveling wave tube, wave-particle interaction, time-domain simulation, Gel'fand βtransform, envelope modulation PACS numbers: 84.40.Fe (microwave tubes) 52.35.Fp (Plasma: electrostatic waves and oscillations) 52.40.Mj (particle beam interaction in plasmas) 52.20.Dq (particle orbits)

We discuss the envelope modulation assumption of frequency-domain models of traveling wave tubes (TWTs) and test its consistency with the Maxwell equations. We compare the predictions of usual frequency-domain models with those of a new time domain model of the TWT.

I. INTRODUCTION

Two approaches are usually applied in describing the dynamics in a traveling wave tube (TWT) : (i) the time domain approach makes no assumption on the TWT working frequency, nor even on the shape of the signal wave, while (ii) the frequency domain approach implies that all variables of interest depend on time like e iωt . Socalled multi-frequency models can accommodate several frequencies, provided they are multiples of a fundamental frequency. This hampers the prediction of non-harmonic instabilities, as may result from nonlinearities or defects.

Current simulations for TWT design rely mainly on frequency domain models because of their fast numerical execution. However, the current need for higher power and gain brings these models to the limits of their reliability. Therefore, we develop a family of time domain models [START_REF] André | Hamiltonian description of self-consistent wave-particle dynamics in a periodic structure[END_REF][START_REF] André | Hamiltonian description of the electron and its radiated field in a periodic structure[END_REF][START_REF] André | First comparison of new TWT discrete model with existing models[END_REF][START_REF] Bernardi | Utilisation et amélioration du modèle discret d'excitation d'un guide d'onde périodique pour la simulation pratique du tube à onde progressive en domaine temporel[END_REF] inspired by [START_REF] Kuznetsov | On one form of excitation equations of a periodic waveguide[END_REF][START_REF] Ryskin | Non-stationary non-linear modeling of an electron beam interaction with a coupled cavity structure, i. Theory[END_REF], in order to provide more compact, accurate and complete descriptions of such regimes. This paper compares this time-domain approach and frequency models, as well as the foundations of the latter.

In section II, we revisit a central assumption of industrial frequency models [START_REF] Antonsen | christine : a multifrequency parametric simulation code for traveling wave tube amplifiers[END_REF][START_REF] Li | The beam-wave interaction in a Ka-band relativistic coupled-cavity TWT[END_REF][START_REF] Li | Theory and design of microwave-tube simulator suite[END_REF][START_REF] Waller | Modélisation numérique de l'interaction et diagnostic expérimental du faisceau d'électrons dans un tube à ondes progressives spatial[END_REF] and discuss its consistency with Maxwell equations when applied to threedimensional geometry. Then, we rederive the basic equation for the modulation amplitude in frequency models. In section III, we compare frequency-domain models to assess the importance of the selection of modes which must be made in order to run the model. In section IV, we recall the principles of our discrete model, which we compare to frequency-domain models in section V to assess their mutual consistency.

II. CONSISTENCY OF FREQUENCY MODELS

Frequency domain models used in the industry share a common representation of the electromagnetic field interacting with the electron beam. The "hot" field at frequency ω is simply the sum of "beamless" modes modulated by an envelope factor function of the axial position z only. One [START_REF] Waller | Modélisation numérique de l'interaction et diagnostic expérimental du faisceau d'électrons dans un tube à ondes progressives spatial[END_REF] or several [START_REF] Antonsen | christine : a multifrequency parametric simulation code for traveling wave tube amplifiers[END_REF][START_REF] Li | Theory and design of microwave-tube simulator suite[END_REF] modes can be used in this expansion. The term "mode" refers to the different propagating modes existing in the delay line at the considered frequency ω, each mode m having a real propagating constant β m as shown in Fig. 1. Consequently, we write in harmonic form

E c (r) = m C m (z)E m (r), (1)
H c (r) = m C m (z)H m (r). (2)
The envelope factors are the C m 's, and z is the longitudinal coordinate along the tube axis, with unit vector e z . In this paper, we call such models Cold Wave Amplification Models (CoWAMs). The real electric field is

1 2
E c (r)e iωt when only one frequency is considered. The case of several frequencies, in particular harmonics, can be treated by summing the expansions at each considered frequency. Hereafter, only one frequency is considered.

By definition, the cold fields E m and H m solve the homogeneous Maxwell equations in harmonic form, viz. the Helmholtz equation,

∇ ∧ E m = -iωµ 0 H m , (3)
∇ ∧ H m = iω 0 E m , (4)
with appropriate boundary conditions. In presence of a beam, the physical fields E(r) and H(r) must satisfy the full Maxwell equations, in particular Maxwell-Faraday for ∇ ∧ E. The electric field then reads where the space-charge field E sc is a gradient. Then the curl of the full field E(r)e iωt is only the curl of the circuit field (1),

E(r) = E sc (r) + E c (r), (5)
∇ ∧ E(r) = ∇ ∧ m C m (z)E m (r) = m dC m dz e z ∧ E m (r) + C m (z)∇ ∧ E m (r) = m dC m dz e z ∧ E m (r) -iC m (z) ω m µ 0 H m (r) , (6)
where the second equality follows from a vector analysis identity and the third equality results from (3). Finally, recalling (2) yields

∇ ∧ E c (r, t) + µ 0 ∂ t H c (r, t) = m dC m dz e z ∧ E m (r) e iωt .
(7) The Maxwell-Faraday equation states that the left-hand side must vanish. Therefore, the envelope form (1)-( 2) is consistent with the Maxwell equations only if the right hand side vanishes. If the electric field of each mode is purely longitudinal, this condition is met automatically.

However, generically the electric field is not purely longitudinal through the whole cross section of the structure. Otherwise the Poynting vector would be purely transverse and the considered mode would carry no power. This precludes the possibility to use a single mode in the expansion [START_REF] André | Hamiltonian description of self-consistent wave-particle dynamics in a periodic structure[END_REF]. Indeed, the only possibility to satisfy the Maxwell-Faraday equation with a single cold mode is that ∂ z C 1 = 0, i.e. that no amplification occurs. In practice, C 1 is slowly varying, so that ∂ z C 1 remains low enough for a good approximation.

The expansion on a single mode violates the Maxwell equations, so the question is now to determine whether they can be satisfied by an expansion on several modes.

To this end, consider the Maxwell-Ampère equation

∇ ∧ H = J(r, t) + 0 ∂ t E (8)
and recall that H = H c . Again, substituting (1)-( 2)-( 5) and using the same vector identity yields

m dC m dz e z ∧ H m e iωt = J + 0 ∂ t E sc . (9)
Denote by Jω (r) the Fourier component of J + 0 ∂ t E sc . Multiplying ( 9) with E * n and using a vector identity obtains

m dC m dz (H m ∧ E * n ) • e z = Jω • E * n . (10)
Similarly, dot-multiplying [START_REF] Doveil | Experimental observation of the nonlinear synchronization due to a single wave[END_REF] with H * n and using a vector identity yields

0 = m dC m dz (E m (r) ∧ H * n (r)) • e z ( 11
)
at every point r. Subtracting ( 11) from [START_REF] Li | Theory and design of microwave-tube simulator suite[END_REF] and integrating over a transverse, planar section S yields

m dC m dz S (H m ∧ E * n -E m ∧ H * n ) • e z d 2 r = S Jω • E * n d 2 r. (12)
For a periodic structure (as considered in mvtrad [START_REF] Waller | Modélisation numérique de l'interaction et diagnostic expérimental du faisceau d'électrons dans un tube à ondes progressives spatial[END_REF], christine [START_REF] Antonsen | christine : a multifrequency parametric simulation code for traveling wave tube amplifiers[END_REF] and bwis [START_REF] Li | Theory and design of microwave-tube simulator suite[END_REF]), the integral in the left hand side vanishes for m = n, thanks to the mode orthogonality relation

[6] 1 4 S (E m ∧ H * n + E * n ∧ H m ) • e z d 2 r = δ mn σ m P m ( 13
)
where P m is the absolute value of the (longitudinal) electromagnetic power flow in mode m (incorporating the normalization of the eigenfields E m and H m ) and σ m = ±1 according to the direction of the power flow. Equation ( 12) so reduces to

dC m dz = - σ m 4P m S Jω • E * m d 2 r. (14)
These are necessary, but not sufficient, conditions to satisfy the Maxwell equations. Equations ( 14) completely determine the amplitude of the wave inside the TWT, along with the dynamics of the electrons. At this point, we can make an important remark on which modes will be amplified or not. At the start of the amplification process, all modes have small amplitude and the TWT operates in the linear regime. We know that in this case the beam is carrying space charge waves [START_REF] Louisell | Coupled mode and parametric electronics[END_REF]. These waves have their own wave number which is given by the beam velocity v 0 in a first order approximation, β 0 = ω/v 0 . Therefore, only modes with β m close to β 0 will grow significantly in equation ( 14), because otherwise the right hand side is the integral of an oscillating function. This will be illustrated in the next section with a numerical example. Now, we compare the relative importance of these waves in term of power by calculating the total power flow from the Poynting vector. The same orthogonality theorem gives (with c.c. denoting complex conjugate)

P = 1 2 S E c ∧ H * c d 2 r + c.c. = 1 2 m |C m | 2 P m , (15)
i.e., the total power results from the power carried by each individual mode. Introducing the mode amplitude a m = C m √ P m [START_REF] Louisell | Coupled mode and parametric electronics[END_REF], the power of each mode is simply |a m | 2 and equation (15) becomes

da m dz = - σ m 2 √ 2 β m Z cm S Jω • E * m d 2 r. ( 16
)
where Z cm is the coupling impedance of mode m at the working frequency and E m is the electric field of that mode divided by its amplitude. With this final form of the interaction equation, we can see that modes with large coupling impedance will rapidly dominate over modes of low coupling impedance as we move toward the TWT output.

To conclude this part, we are left with a contradiction. We have just seen that only the modes synchronized with the beam can grow, as confirmed by dedicated experiments [START_REF] Doveil | Experimental observation of the nonlinear synchronization due to a single wave[END_REF]. In the usual situation of a practical TWT, only one mode satisfies this synchronization condition. Therefore, this mode alone rapidly dominates in the amplification process. But we have also seen that the Maxwell equations cannot be satisfied with a single mode. It follows that the envelope model (CoWAM) can only approximate the physics. How accurate this approximation is remains difficult to assess. One motivation to develop the discrete model, beside its time domain capabilities, is to lift these theoretical approximations. Based on our experience however, we expect that the approximation is minor, at least inside common ranges of parameters used for practical devices. This fact will also be illustrated hereafter with a numerical example.

III. RELEVANCE OF BACKWARD WAVE IN COWAM

In this numerical example, we assess the relative importance of the forward and backward modes for a standard helix tube. Such a tube comprises attenuating sections, where the fields are partly absorbed ; these attenuating sections are modeled in CoWAMs on adding a loss term in [START_REF] Waller | Modélisation numérique de l'interaction et diagnostic expérimental du faisceau d'électrons dans un tube à ondes progressives spatial[END_REF] to give where the spatial loss rate α is a function of position accounting for propagation losses and localized attenuator. The purpose of the present calculation is to compare the relative power of the different modes radiated by a given modulated electron beam. Consequently, the current density J is an input generated from a preliminary run of mvtrad. The same modulated current density is used to calculate the radiation on the forward mode alone or on the forward and backward modes together. The differences in the resulting electric field acting on the electron trajectories is not taken into account, i.e. these trajectories are frozen. We first compare the power radiated in the model taking a single mode into account (m = 1) versus taking also the backward wave into account (m = ±1). The full line on Fig. 2 results from integrating (17) given a modulation current computed with only one mode and computing the power gain as

dC 1 dz = - σ 1 4P 1 S E * 1 (z) • J(z) d 2 r -α (z)C 1 (z) (17)
P dB = 10 log 10 |C 1 | 2 C 2 0 ( 18
)
where C 0 = C 1 (0) is the envelope amplitude at the TWT inlet.

The broken line is obtained similarly for m = -1, using also mvtrad formally. Indeed, though mvtrad does not use the backward mode to compute the wave amplification, (17) can be integrated with the actual current density and the eigenfield E -1 , to check whether the beam modulations might be resonant with the backward mode. Fig. 2 shows that the backward mode is poorly coupled for this case, where the phase shift per cell βd is moderate : near the outlet, the backward mode reaches 5 dB while the direct mode culminates at 45 dB, and through the whole tube the backward mode is always at least 20 dB weaker than the direct mode. Quite obviously, the forward wave is significantly damped by the two attenuating sections, inserted to avoid multiple reflexions in each section. The attenuation of the backward wave is less visible because its power is not building up, the synchronization condition being not fulfilled. It results that the power radiated by the beam at any location is not negligible compared to the incoming power carried by the wave, including inside the attenuating section.

The role of backward waves is more relevant when the TWT is operated near the π mode because they synchronize with the beam at this operating frequency. This is confirmed by our simulations. For this simulation, we need a much higher frequency (out of the actual operation band), at which the coupling impedance is much smaller so that amplification occurs more slowly. Near the π mode, the backward wave couples efficiently with the slow space charge wave of the beam, generating the backward wave instability, where there is a risk of developing an oscillation. To stress the effect, our simulation is run without attenuation, and shows that the backward wave generated near the outlet of the TWT has an intensity exceeding significantly the input of the direct wave (here by about 8 dBm at z = 0). So, near the π mode indeed, the CoWAM must incorporate the backward wave, and a model like bwis will definitely be more accurate than mvtrad.

However, a CoWAM is unable to predict an instability if it is not foreseen by a genuine identification of the modes to be included in the run. Time-domain models will not suffer from this lack of prescience.

IV. DISCRETE MODEL

In contrast with the previous models, our model dimo is formulated in the time domain [START_REF] André | First comparison of new TWT discrete model with existing models[END_REF] and does not make use of the same expansion on cold fields that has been critically reviewed in section II. It implements the analysis presented in Ref. [START_REF] André | Hamiltonian description of self-consistent wave-particle dynamics in a periodic structure[END_REF] and can describe any delay lines : folded waveguides, helices, coupled cavities. . .

We outline here the derivation of the discrete model which has been presented in detail in the past.

To take advantage of the periodicity of the slow wave structure, one first determines the solutions E sβ (r), H sβ (r) to the Helmholtz equation in a periodic cell (denoted V 0 , typically 0 ≤ z ≤ d, with a given profile (x, y) ∈ S(z), e.g. x 2 + y 2 ≤ R 2 for a cylinder), with the Floquet boundary condition E β (r + de z ) = e -iβd E β (r) for the solenoidal (or circuit) fields E c , H c . For each propagation constant 0 ≤ β < 2π/d, these solutions form a basis with eigenfrequencies Ω sβ . Subscript s labels different eigenmodes meeting the same Floquet condition, i.e. the different bands (four of them are displayed in Fig. 1).

Physical functions G(r, t) are expressed in terms of these eigenmodes by first applying the Gel'fand βtransform [START_REF] Kuznetsov | On one form of excitation equations of a periodic waveguide[END_REF]

G(r + nde z , t) = d 2π 2π/d 0 G β (r, t) e -inβd dβ, ( 19
)
G β (r, t) = +∞ n=-∞ G(r + nde z , t) e inβd , ( 20
)
and expanding E β and H β on the Floquet divergence-free eigenbasis of the slow wave structure,

E β (r, t) = s V sβ (t)E sβ (r) -∇φ β , (21)
H β (r, t) = i s I sβ (t)H sβ (r), ( 22
)
where the first equation includes the irrotational spacecharge field contribution, described with the potential φ, and the imaginary unit is introduced so that final physical variables of the model are purely real numbers. This expansion is the essential difference between the discrete and the envelope models. The E sβ are the (vector) eigenfunctions of the Helmholtz equation, therefore they clearly constitute a basis on which any divergencefree vector field can be expanded at any time. The discrete model makes no approximation in the function choice. On the contrary, the expansion (1)-( 2) is performed on the propagating modes at a given frequency which apparently do not necessarily constitute a basis : at a frequency below the cut-off of a waveguide, for example, evanescent modes do exist although there are no propagating modes on which they could be expanded. Even at a frequency where propagation is possible, the modes with higher cut-off frequency are still possible in the form of evanescent modes but cannot be described by the modes propagating at this frequency.

The time-dependent coefficients V sβ , I sβ should generally not coincide (nor be merely proportional) as they will obey their own coupled evolution equations, whereas the envelopes C m in (1)-( 2) had to coincide to describe the same modulation of cold fields in CoWAMs.

The Maxwell equations for the field propagation translate into evolution equations for the Floquet coefficients V sβ , I sβ :

Vsβ + Ω sβ I sβ = - 1 N sβ V0 J β (r, t) • E * sβ (r) d 3 r,( 23
)
İsβ -Ω sβ V sβ = 0, (24)
where N sβ is the electromagnetic energy of mode (s, β) in a unit cell of the slow wave structure, and J is the beam electric current density. Imposing a constant ratio V sβ /I sβ in the discrete model (by analogy with the C m 's) would be an extra condition, making the set of equations for the coefficients overdetermined. Keeping linearly independent coefficients V sβ , I sβ enables the discrete model to satisfy Maxwell equations, in contrast with CoWAMs.

On introducing the field (closely related to the vector potential, see [START_REF] André | Hamiltonian description of self-consistent wave-particle dynamics in a periodic structure[END_REF])

F s,n (r) := 1 2π 2π 0 E sβ (r) N sβ e -inβd d(βd), (25)
these equations read in real space

Vsn + m Ω s,m I s,n-m = - V Z J(r, t) • F s,-n (r) d 3 r, (26) İsn -
m Ω s,m V s,n-m = 0, ( 27
)
where V Z is the full extent of the slow wave structure.

In real space and time variables, (21)-( 22) give the fields in terms of these coefficients as

E(r, t) = s,n V sn (t) E s,-n (r) -∇φ(r, t), (28)
H(r, t) = i s,n I sn (t) H s,-n (r), (29)
where φ(r, t) is the beam space charge potential.

The detailed modeling of the slow wave structures lies (i) in the frequency matrix Ω s , which is typically a band matrix Ω s,n,m = Ω s,n-m , symmetric (Ω s,n,m = Ω s,m,n by reciprocity condition), with a rather short range (Ω s,n-m = 0 if |n -m| > p with, say, p = 1 for coupled cavities [START_REF] Ryskin | Non-stationary non-linear modeling of an electron beam interaction with a coupled cavity structure, i. Theory[END_REF] and p ∼ 5 for a helix), and (ii) in the explicit functions E s,n (r), F s,n (r), H s,n (r), whose construction involves the coupling impedances Z sβ over the relevant bandwidth. As a first approximation, one mode s suffices to capture the physics of the TWT. In contrast with the local picture of envelope modulation (1)-( 2), the fields E s,n and F s,n should not be viewed as local to a cell n but rather may have quite long a range in terms of n to express how the cell couples with the beam.

The beam coupled with the wave is described by macro-electrons, with charge to mass ratio -η = -|e|/m e and mass m, position r k (t) and velocity v k = ṙk , injected at cathode potential -V K at a constant rate to match the physical current I beam at the electron gun. Their equation of motion then reads

d dt [(1-|v k | 2 /c 2 ) -1/2 v k ] = -η[E(r k , t)+µ 0 v k ∧H(r k , t)].
(30) The space-charge potential solves the Poisson equation with particles as sources [START_REF] Rowe | Nonlinear electron-wave interaction phenomena[END_REF] and boundary conditions fixed by the slow wave structure.

V. COMPARISON OF FREQUENCY AND TIME DOMAIN MODELS

We now compare the cold wave amplification model in frequency domain with the discrete model in time domain. Given the experimental characteristics (dispersion diagram, coupling impedance) of the tube, we interpolated the eigenmode electric field on the axis E β (z) and energy density N β , and constructed the matrix Ω n-m and interaction field F n (z) defined by the Gel'fand transform.

Once again, the current modulation J is obtained from a preliminary mvtrad run and is injected in equation (26). This procedure permits to compare the discrete model and CoWAM with exactly the same current modulation. Indeed, these two models differ only in how electromagnetic waves are radiated from the beam, not in the dynamics of electrons given by (30).

For mvtrad and the bwis-like model respectively, we solved ( 14) with m = 1 and m = ±1, as in section III. All three models are run lossless (viz. we insert no attenuation).

Fig. 4 displays the power in the first harmonic along the tube, normalized to the input power. All three models agree qualitatively, and rather well quantitatively, for the growth of the main harmonic. Because of their richer harmonic structure, bwis and dimo show more oscillations near the entrance to the tube, but these oscillations become negligible once the field becomes so intense that the first harmonic dominates. The end of the tube was not displayed on this figure, to zoom on the region where the contrast between models is larger. The shoulder near z = 80 mm is a saturation effect before the amplification resumes, and it is identically predicted by the three models, confirming their agreement.

Though these simulations were run without attenuation, dimo can accommodate these by an additional loss term in (26), which becomes Vsn = -V Z J(r, t)•F s,-n (r) d 3 rm Ω s,m I s,n-m -α sn V sn , (31) with a localized positive time-decay rate α sn = α (z)v g for z in cell n, with v g the group velocity at the tube operating frequency. Runs with attenuation confirm the agreement between all three models.

VI. CONCLUSION

We have first shown that the spatial envelope modulation assumed in CoWAMs, which made frequency models efficient for longitudinal, one-dimensional simulations, is inconsistent with the Maxwell equations in the general case. This encourages to search for alternative modelings, and the proposed discrete model is free from this failing. Moreover, since time-domain models are not restricted to a prescribed family of frequencies, they are good candidates for investigating nonlinear regimes and the appearance of unplanned resonances, such as a backward wave and drive induced oscillations.

Our second observations compare two one-dimensional CoWAMs and a simple time-domain model, dimo. On the one hand, we see that the relevant modes for a CoWAM can be reasonably predicted on the basis of their possible resonance with the beam (the backward wave is found negligible in fig. 2). On the other hand, we show that the time-domain simulation reproduces well the well-tested amplification regime where CoWAMs are reliable.

These results show the prospects opened by timedomain direct simulation using a compact discrete model with β-transformed basis fields, instead of a full electromagnetic model. An improved version of dimo is currently under development, taking advantage of the explicit hamiltonian nature of particle-wave dynamics.

FIG. 1 :

 1 FIG. 1: Typical dispersion diagram of a helix supported by dielectric rods inside a metallic cylinder. The helix can propagate six distinct modes at the chosen frequency.

FIG. 2 :

 2 FIG. 2: Power along the axis for the first harmonic. Full line : contribution from coefficient C+1. Broken line : contribution from coefficient C-1. Parameters : F = ω/(2π) = 11.5 GHz, tube length 150 mm, attenuating sections centered at 50 mm and 100 mm.

FIG. 3 :

 3 FIG. 3: Same as Fig. 2, near the π mode for β = π/d and no attenuation.

FIG. 4 :

 4 FIG. 4: First harmonic of the electric field along the axis. Thin solid red line : dimo. Blue dots : forward model mvtrad. Thick solid green line : forward and backward (bwis-like) model. Parameters : F = 11.5 GHz, tube length 150 mm, no attenuation. Ordinate scale shows the power rather than field amplitude.

Acknowledgment

The authors are pleased to thank P. Bernardi, F. Doveil and D. Minenna for fruitful discussions. S. Théveny was supported by a CIFRE doctoral grant.