René Lozi

Clarisse Fiol

Global Orbit Patterns for Dynamical Systems On Finite Sets

Keywords: dynamical systems, chaotic analysis, combinatorial dynamics, global orbit pattern, locally rigid functions PACS: 05.10-a, 05.45-a, 02.10.0x √

In this paper, the study of the global orbit pattern (gop) formed by all the periodic orbits of discrete dynamical systems on a finite set X allows us to describe precisely the behaviour of such systems. We can predict by means of closed formulas, the number of gop of the set of all the function from X to itself. We also explore, using the brute force of computers, some subsets of locally rigid functions on X, for which interesting patterns of periodic orbits are found.

INTRODUCTION

In some engineering applications such as chaotic encryption, chaotic maps have to exhibit required statistical and spectral properties close to those of random signals. There is a growing industrial interest to consider and study thoroughly the property of such map [10,[START_REF]US Patent 6,999,445 -Multiple access communication system using chaotic signals and method for generating and extracting chaotic signal[END_REF]12].

Very often, dynamical systems in several dimensions are obtained coupling 1dimension ones and their properties are strongly linked [START_REF] Lozi | The importance of strange attractors for industrial mathematics[END_REF].

Quasi-periodic or chaotic motion is frequently present in complicated dynamical systems whereas simple dynamical systems often involve only periodic motion. The most famous theorem in this field of research is the Sharkovskii's theorem, which addresses the existence of periodic orbits of continuous maps of the real line into itself. This theorem was once proved toward the year 1962 and published only two years after [START_REF] Sharkovskii | Coexistence of cycles of continuous mapping of the line into itself[END_REF].

Mathematical results concerning periodic orbits are often obtained for functions on real intervals. However, most of the time, as the complex behaviour of chaotic dynamical systems is not explicitly tractable, mathematicians have recourse to computer simulations. The main question which arises then is: does these numerical computations are reliable ?

As an example we report the results of some computer experiments on the orbit structure of the discrete maps on a finite set which arise when the logistic map is iterated "naively" on the computer.

Due to the discrete nature of floating points used by computers, there is a huge gap between these results and the theoretical results obtained when this map is considered on a real interval. This gap can be narrowed in some sense (i.e. avoiding the collapse of periodic orbits) in higher dimensions when ultra weak coupling is used [START_REF] Lozi | Giga-periodic Orbits for Weakly Coupled tent and Logistic Discretized Maps, Modern Mathematical Models, Methods and Algorithms for Real World Systems[END_REF][START_REF] Lozi | New Enhanced Chaotic Number Generators[END_REF].

Nowadays the claim is to understand precisely which periodic orbit can be observed numerically in such systems. In a first attempt we study in this paper the orbits generated by the iterations of a one-dimensional system on a finite set X N with a cardinal N. The final goal of a good understanding of the actual behaviour of dynamical systems acting on floating numbers (i.e. the numbers used by computers) will be only reached after this first step will be achieved.

On finite set, only periodic orbits can exist. For a given function we can compute all the orbits, all together they form a global orbit pattern. We formalise such a gop as the ordered set of periods when the initial value thumbs the finite set in the increasing order. We are able to predict, using closed formulas, the number of gop for the set F N of all the functions on X . We also explore by computer experiments special subsets of F N , such as sets of locally "rigid" functions which presents interesting patterns of gop.

This article is organized as follows : in the section "Computational divergences" we display some examples of such computational divergences for the logistic map in various ways of discretization. In the section "Pattern defined by all the orbits of a dynamical system" we introduce a new mathematical tool: the global orbit pattern, in order to describe more precisely the behaviour of dynamical systems on finite sets. In the section "Cardinal of classes" we give some closed formulas related to the cardinal of classes of gop of F N . In the section "Functions with local properties" we study the case of sets of functions with a kind of local "rigidity" versus their gop, in order to show the usefulness of these new tools.

COMPUTATIONAL DIVERGENCES

Discretized logistic map

As an example of collapsing effects which happen when using computers in numerical experiments, we presents the results of a sampling study in double precision of a discretization of the logistic map f 4 : [0, 1] → [0, 1] (see Fig. 1)

f 4 (x) = 4x(1 -x)
(1) and its associated dynamical system

x n+1 = 4x n (1 -x n )
(2) which has excellent ergodic properties on the real interval. There exists an unstable fixed point 0. The set 5-

P(x) = 1 π x(1 -x) (3)
However, in numerical computations using ordinary (IEEE-754) double precision numbers -so that the working interval contains of the order of 10 16 representable points -out of 1,000 randomly chosen initial points (see Table 1),

• 596, i.e., the majority, converged to the fixed point corresponding to the unstable fixed point {0} in equation 2, • 404 converged to a cycle of period 15,784,521. Thus, in this case at least, the very long-term behaviour of numerical orbits is, for a substantial fraction of initial points, in flagrant disagreement with the true behaviour of typical orbits of the original smooth logistic map.

In others numerical experiments we have performed, the computer working with fixed finite precision is able to represent finitely many points in the interval in question. It is probably good, for purposes of orientation, to think of the case where the representable points are uniformly spaced in the interval. The true logistic map is then approximated by a discretized map, sending the finite set of representable points in the interval to itself.

Describing the discretized mapping exactly is usually complicated, but it is roughly the mapping obtained by applying the exact smooth mapping to each of the discrete representable points and "rounding" the result to the nearest representable point. In our experiments uniformly spaced points in the interval with several order of discretization (ranging from 9 to 2,001 points) are involved. In each experiment the questions addressed are:

• how many periodic cycles are there and what are their periods ?

• how large are their respective basins of attraction, i.e. , for each periodic cycle, how many initial points give orbits with eventually land on the cycle in question ? On an another hand, for relatively coarse discretizations the orbit structure is determined completely, i.e., all the periodic cycles and the exact sizes of their basins of attraction are found. Some representative results are given in Tables 2 to 4. In theses tables, N denotes the order of the discretization, i.e., the representable points are the numbers, j N , with 0 ≤ j < N.

The Table 2 shows coexisting periodic orbits for the discretization with regular meshes of N = 9, 10 and 11 points. There are exactly 3, 2 and 2 cycles. The Table 3 shows coexisting periodic orbits for the discretization with regular meshes of N = 99, 100 and 101 points. There are exactly 2, 4 and 3 cycles. The Table 4 shows coexisting periodic orbits for the discretization with regular meshes of N = 1, 999, N = 2, 000 and N = 2, 001 points.

It seems that the computation of numerical approximations of the periodic orbits leads to unpredictable results.

Statistical properties

Many others examples could be given, but those given may serve to illustrate the intriguing character of the results: the outcomes proves to be extremely sensitive to the details of the experiment, but the results all have a similar flavour : a relatively small number of cycles attract near all orbits, and the lengths of these significant cycles are much larger than one but much smaller than the number of representable points.

In [START_REF] Diamond | Statistical laws for computational collapse of discretized chaotic mappings[END_REF], P. Diamond and A. Pokrovskii, suggest that statistical properties of the phenomenon of computational collapse of discretized chaotic mapping can be modelled by random mappings with an absorbing centre. The model gives results which are very much in line with computational experiments and there appears to be a type of universality summarised by an Arcsine law. The effects are discussed with special reference to the family of mappings

x n+1 = 1 -|1 -2x n | ℓ 0 ≤ x ≤ 1 1 ≤ ℓ ≤ 2 (4)
Computer experiments show close agreement with prediction of the model.

However these results are of statistical nature, they do not give accurate information on the exact nature of the orbits (e.g. length of the shortest one, of the greater one, size of their basin of attraction ...). It is why we consider the problem of computational discrepancies in an original way in the next section.

PATTERN DEFINED BY ALL THE ORBITS OF A DYNAMICAL SYSTEM

In this section in order to describe more precisely which kind of behaviour occurs in discretized dynamical systems on finite sets we conceive a new mathematical tool: the global orbit pattern of a function that is the set of the periods of every different orbits of the dynamical system associated to the function when the initial points are took in increasing order.

General definitions

For every x 0 ∈ X , let {x i } be the sequence of the orbit of the dynamical system associated to the function f which maps X onto X defined by

x i+1 = f (x i ) for i ≥ 0. ( 5
)
For convenience ∀x 0 ∈ X we denote

f 0 (x 0 ) = x 0 (6) and ∀p ≥ 1, ∀x 0 ∈ X , f p (x 0 ) = f • f • . . . • f p times (x 0 ). (7)
Hence

x i = f i (x 0 ). (8)
The orbit of x 0 under f is the set of points O(x 0 , f ) = { f i (x 0 ), i ≥ 0} = {x i , i ≥ 0}. The starting point x 0 for the orbit is called the initial value of the orbit. A point x is a fixed point of the map f if f (x) = x. A point x is a periodic point with period p if f p (x) = x and f k (x) = x for all k such that 0 ≤ k < p, p is called the order of x. If x is periodic of order p, then the orbit of x under f is the finite set {x, f (x), f 2 (x), . . ., f p-1 (x)}. We will call this set the periodic orbit of order p or a p-cycle. A fixed point is then a 1-cycle. The point x is an eventually periodic point of f with order p if there exists K > 0 such that ∀k ≥ K f k+p (x) = f k (x). ∀x ∈ X , we denote ω(x, f ) the order of x under f or simply ω(x) when the map f involved is obvious.

A subset T of X is invariant under f if f -1 (T ) = T . That is equivalent to say that T is invariant under f if and only if f (T ) ⊂ T and f -1 (T ) ⊂ T .
Notation ♯X is the cardinal of the finite set X .

Map on finite set

Along this paper, N is a non-zero integer and ♯A stands for the cardinal of any finite set A. In this article, we consider X as an ordered finite set with N elements. We denote it X N , it is isomorphic to the interval [[0, N -1]] ⊂ N. Then ♯X N = N. Let f be a map from X N into X N . We denote by F N the set of the maps from X N into X N . Clearly, F N is a finite set and ♯F N = N N elements. For all x ∈ X N , O(x, f ) is necessarily a finite set with at most N elements. Indeed, let us consider the sequence {x, f (x), f 2 (x), . . ., f N-1 (x), f N (x)} of the first N + 1 iterated points. Thanks to the Dirichlet's box principle, two elements are equals because X N has exactly N different values. Thus, every initial value of X N leads ultimately to a repeating cycle. More precisely, if x is a fixed point O(x, f ) is the unique element x and if x is a periodic point with order p, O(x, f ) has exactly p elements. In this case, the orbit of

x under f is the set O(x, f ) = {x, f (x), f 2 (x), . . . , f p-1 (x)}. If x is an eventually periodic point with order p, there exists K > 0 such that ∀k ≥ K f k+p (x) = f k (x). In this case, the orbit of x under f is the set O(x, f ) = {x, f (x), f 2 (x), . . . , f K (x), f K+1 (x), . . ., f K+p-1 (x)}.

Equivalence classes

Components

Let f ∈ F N . We consider on X N the relation ∼ defined by : ∀x,

x ′ ∈ X N , x ∼ x ′ ⇔ ∃k ∈ N such that f k (x) ∈ O(x ′ , f ). The relation ∼ is an equivalence relation on X N .
S N / ∼ is the collection of the equivalence classes that we will call components of X N under f which constitute a partition of X N . The number of components are given in [START_REF] Kruskal | The expected number of components under a random mapping function[END_REF]. Asymptotic properties of the number of cycles and components are studied in [START_REF] Ljuben | Limit theorem concerning random mapping patterns[END_REF]. For each component, we take as class representative element the least element of the component. The components will be written

T N (x 0 , f ), . . ., T N (x p f ,N , f ) where x i is the least element of T N (x i , f ).
By analogy with real dynamical systems, we can define attractive and repulsive components in discretized dynamical systems as follows.

Definition 1 A component is repulsive when it is a cycle. Otherwise, the component is attractive.

Remark

In other words, a component is attractive when the component contains at least an eventually periodic element. The corresponding cycle is strictly contained in an attractive component.

Examples are given in Tables 5, 6 and7. For instance, in Table 6, the fonction f has {2, 7} as period-2 orbit and {1, 2, 7, 9} as component which is attractive because 1 and 2 are eventually periodic elements.

Order of elements

Here are some remarks on the order of elements of components.

Remark

The order of every element of a component is the length of its inner cycle.

Definition 2 For all x ∈ X N , there exists i ∈ [[0, p f ,N ]] such that x belongs to the component T N (x i , f ). Then ω(x, f ) is equal to the order ω(x i , f ). Remark For all i ∈ [[0, p f ,N ]], T N (x i , f ) is an invariant subset of X N under f .
In the example given in Table 5, the order of the element 0 is 2, the order of the element 1 is 2, the order of the element 4 is 3. The elements 1 and 3 have the same order.

Definition of global orbit pattern

For each f ∈ F N , we can determine the components of X N under f . For each component, we determine the order of any element. Thus, for each f ∈ F N , we have a set of orders that we will denote Ω( f , N). Be given f , there exist p f ,N components and p f ,N representative elements such that x 0 < x 1 < . . . < x p f ,N .

For each f ∈ F N , the sequence [ω(x 0 ), ω(x 1 ), . . . , ω(x p f ,N ); f ] F N with x 0 < x 1 < . . . < x p f ,N will design the global orbit pattern of f ∈ F N .

We will write gop

( f ) = [ω(x 0 ), ω(x 1 ), . . ., ω(x p f ,N ); f ] F N .
When the reference to f ∈ F N is obvious, we will write shortly

gop( f ) = [ω(x 0 ), ω(x 1 ), . . ., ω(x p f ,N )] N or gop( f ) = [ω(x 0 ), ω(x 1 ), . . ., ω(x p f ,N )] .
For example, the same gop associated to the functions given in Tables 5, 6 and7

is [2, 2, 1, 3] 11 .
Another example is given in Table 8. In that example, we have ω(0

) = 2, ω(3) = 1, ω(4) = 4.
F N is called G (F N ).
For example, for N = 5, the set

G (F 5 ) is {[1]; [1, 1]; [1, 1, 1]; [1, 2]; [1, 1, 1, 1]; [1, 1, 2]; [1, 3]; [1, 1, 1, 1, 1]; [1, 1, 1, 2]; [1, 1, 2, 1]; ; [1, 2, 1, 1]; [1, 2, 2]; [1, 3, 1]; [1, 4]; ; [2, 1]; [2, 1, 1]; [2, 2]; [2, 1, 1, 1]; [2, 1, 2]; [2, 2, 1] ; [3]; [3, 1]; [3, 1, 1]; [3, 2]; [4]; [4, 1]; [5]}.

Class of gop

We give the following definitions :

Definition 4 Let be A = [ω 1 , . . . , ω p ] N a gop.
Then the class of A, written A, is the set of all the functions f ∈ F N such that the global orbit pattern associated to f is A.

Then the modulus of

A is |A| = p ∑ k=1 ω k . Remark [ω 1 , . . . , ω p ] N ≤ N. Notation [ω k ] N means [ω, . . ., ω k times ] N and [ω k , ν m ] N means [ω, . . . , ω k times , ν, . . . , ν m times ] N .

Threshold functions

Ordering the discrete maps

Theorem 1 The sets F N and [[1, N N ]] are isomorphic. Proof We define the function φ from F N to [[1, N N ]] by : for each f ∈ F N , φ ( f ) is the integer n such that n = N-1 ∑ k=0 f (k)N N-1-k + 1. Then φ is well defined because n ∈ [[1, N N ]].
Let n be a given integer between 1 and N N . We convert n -1 in base N : there exists a unique N-tuple (a n-1,0 ; a n-1,1 ; . . .

; a n-1,N-1 ) ∈ [[0, N -1]] N such that n -1 N = N-1 ∑ i=0 a n-1,N-1-i N N-i-1
. We can thus define the map f n with : ∀i ∈ X N ,

f n (i) = a n-1,N-i-1 .
Then φ is one to one.

Remark This implies F N is totally ordered.

Definition 6 Let f ∈ F N . Then n = N-1 ∑ k=0 f (k)N N-1-k + 1 (9)
is called the rank of f .

Threshold functions

Be given a global orbit pattern A, we are exploring the class A.

Theorem 2 For every A ∈ G (F N ), the class A has a unique function with minimal rank.

Definition 7 For every class A ∈ G (F N ), the function defined by the previous theorem will be called the threshold function for the class A and will be denoted by

Tr(A) or Tr(A).

To prove the theorem, we need the following definition :

Definition 8 Let f ∈ F N be a function. Let p a non zero integer smaller than N. Let be x 1 , . . . , x p p consecutive elements of X N . Then x 1 , . . . , x p is a canonical p-cycle in relation to f if ∀ j ∈ [[1, p -1]] , f (x j ) = x j+1 and f (x p ) = x 1 .
Proof Let [ω 1 , . . . , ω p ] be a global orbit pattern of G (F N ). We construct a specific function f belonging to the class [ω 1 , . . . , ω p ] and we prove that the function so obtained is the smallest with respect to the order on F N . With the first ω 1 elements of

[[0, N -1]], that is the set of integers [[0, ω 1 -1]], we construct the canonical ω 1 -cycle : if ω 1 = 1, we define f (0) = 0, else f (0) = 1, f (1) = 2, . . ., f (ω 1 -2) = ω 1 -1, f (ω 1 -1) = 0. Then ∀ j ∈ [[ω 1 -1, ω 1 + N -s -1]], we define f ( j) = 0. Then with the next ω 2 integers [[ω 1 + N -s, ω 1 + N -s + ω 2 -1]
] we construct the canonical ω 2 -cycle. We keep going on constructing for all j ∈ [ [3, p]] the canonical ω j -cycle.

In consequence, we have found a function f belonging to the class [ω 1 , . . ., ω p ]. Assume there exists a function g ∈ F N belonging to the class of f such that g < f . Let I = {i ∈ [[0, N -1]] such that f (i) = 0}. As g < f , there exists i 0 ∈ I such that g(i 0 ) < f (i 0 ). There exists also j 0 such that i 0 ∈ ω j 0 . If f (i 0 ) = i 0 , then ω j 0 = 1, g(i 0 ) < i 0 and then g(i 0 ) / ∈ ω j 0 . Then the global orbit pattern of g doesn't contain anymore 1 as cycle. The global orbit pattern of g is different from the global orbit pattern of f . If f (i 0 ) = i 0 + 1, then g(i 0 ) ≤ i 0 . Either g(i 0 ) = i 0 and then the global orbit pattern of g is changed, or g(i 0 ) < i 0 and we are in the same situation as previously. Thus, in any case, the smallest function belonging to the class [ω 1 , . . ., ω p ] is the one constructed in the first part of the proof.

The proof of the theorem gives an algorithm of construction of the threshold function associated to a given gop. The threshold function associated to the gop [2 2 , 1, 3] 11 is explained in Table 9. Its rank is n = 25, 938, 474, 637.

cal 2-cycle 0 → 1 1 → 0 2 → 3 → 4 → 5 → 6 → 7 → 8 → 9 → 10 → Construction of the last canonical 3-cycle 0 → 1 1 → 0 2 → 3 → 4 → 5 → 6 → 7 → 8 → 9 9 → 10 10 → 8 Construction of the canonical 1- cycle 0 → 1 1 → 0 2 → 3 → 4 → 5 → 6 → 7 → 7 8 → 9 9 → 10 10 → 8 Construction of the canonical 2- cycle 0 → 1 1 → 0 2 → 3 → 4 → 5 → 6 6 → 5 7 → 7 8 → 9 9 → 10 10 → 8
Filling the remaining images with 0

0 → 1 1 → 0 2 → 0 3 → 0 4 → 0 5 → 6 6 → 5 7 → 7 8 → 9 9 → 10 10 → 8
Theorem 3 There are exactly 2 N -1 different global orbit patterns in F N .

That is

♯G (F N ) = 2 N -1. (10)
For example, for N = 4, ♯G (F 4 ) = 2 4 -1 = 15.

Proof

Let p an integer between 1 and N. Consider the set L(p, N) of p-tuples (α 1 , . . ., α p ) ∈ (N * ) p such that α 1 + . . . + α p ≤ N. We write L(N) = {L(p, N), p = 1 . . .N}. L(N) and G (F N ) have the same elements. Then

♯G (F N ) = p=N ∑ p=1 ♯L(p, N) = p=N ∑ p=1 N p = 2 N -1.

Ordering the global orbit patterns

We define an order relation on G (F N ).

Proposition 1 Let A and B be two global orbit patterns of G (F N ). We define the relation ≺ on the set G (F N ) by

A ≺ B iff Tr(A) < Tr(B)
Then the set (G (F N ), ≺) is totally ordered.

Proof As the order ≺ refers to the natural order of N, the proof is obvious.

Let r ≥ 1, p ≥ 1 be two integers. Let [ω 1 , . . ., ω p ] and [ω ′ 1 , . . . , ω ′ r ] be two global orbit patterns of G (F N ). For example, if p < r, in order to compare them, we admit that we can fill [ω 1 , . . . , ω p ] with r -p zeros and write [ω 1 , . . . , ω p ] = [ω 1 , . . . , ω p , 0, . . ., 0].

Proposition 2 Let r ≥ 1, p ≥ 1 be two integers such that p ≤ r. Let A = [ω 1 , . . . , ω p ] and B = [ω ′ 1 , . . . , ω ′ r ] be two global orbit patterns. • If r = p = 1 and ω 1 < ω ′ 1 then A ≺ B. • If r ≥ 2 then * If ω 1 < ω ′ 1 then A ≺ B. * If ω 1 = ω ′ 1 then there exists K ∈ [[2; r]] such that ω K = ω ′ K and ∀i < K, ω i = ω ′ i . • If |A| < |B|, then A ≺ B. • If |A| = |B|, then if ω K < ω ′ K then A ≺ B.
For example, for N = 5, the global orbit patterns are in increasing order :

[1] ≺ [1 2 ] ≺ [1 3 ] ≺ [1, 2] ≺ [1 4 ] ≺ [1 2 , 2] ≺ [1, 2, 1] ≺ [1, 3] ≺ [1 5 ] ≺ [1 3 , 2] ≺ [1 2 , 2, 1] ≺ [1 2 , 3] ≺ [1, 2, 1 2 ] ≺ [1, 2 2 ] ≺ [1, 3, 1] ≺ [1, 4] ≺ [2] ≺ [2, 1] ≺ [2, 1 2 ] ≺ [2 2 ] ≺ [2, 1 3 ] ≺ [2, 1, 2] ≺ [2 2 , 1] ≺ [2, 3] ≺ [3] ≺ [3, 1] ≺ [3, 1 2 ] ≺ [3, 2] ≺ [4] ≺ [4, 1] ≺ [5].

Algorithm for ordering the global orbit patterns : a pseudo-decimal order

The Table 10 gives a method for ordering the gop : indeed, we consider each gop as if each one represents a decimal number : we begin to order them in considering the first order ω 1 . Considering two gops A = [ω 1 , ..., ω p ] and

A ′ = [ω ′ 1 , ..., ω ′ r ], if ω 1 < ω ′ 1 , then A ≺ A ′ . For example, [2, 1, 2] ≺ [4, 1]. If ω 1 = ω ′ 1 and |A| -ω 1 < |A ′ | -ω ′ 1 , then A ≺ A ′ .
For example to compare the gop [START_REF] Diamond | Statistical laws for computational collapse of discretized chaotic mappings[END_REF][START_REF] Knuth | The Art Of Computer Programming[END_REF] and the gop [1 4 ], we say that the first order ω 1 stands for the unit digit -which is ω 1 = 1 here, then the decimal digits are respectively 0.2 and 0.111. We calculate for each of them the modulus-ω 1 : we find

|[1, 2]| -1 = 2 and |[1 4 ]| -1 = 3, thus [1, 2] ≺ [1 4 ]. Finally, if ω 1 = ω ′ 1 and |A| -ω 1 = |A ′ | -ω ′
1 , then also we use the order of the decimal part. For example, [1 5 ] ≺ [1, 2, 1 2 because 1.1111 < 1.211. Applying this process, we have the sequence of the ordered gop for N = 4 given in the previous paragraph.

Gop

Modulus Modulus-ω 1 Gop Modulus Modulus-ω 1

[1] 1 0 [2] 2 0 [1 2 ] 2 1 [2, 1] 3 1 [1 3 ] 3 2 [2, 1 2 ] 4 2 [1, 2] 3 2 [2 2 ] 4 2 [1 4 ] 4 3 [2, 1 3 ] 5 3 [1 2 , 2] 4 3 [2, 1, 2] 5 3 [1, 2, 1] 4 3 [2 2 , 1] 5 3 [1, 3] 4 3 [2, 3] 5 3 [1 5 ] 5 4 [1 3 , 2] 5 4 [3] 3 0 [1 2 , 2, 1] 5 4 [3, 1] 4 1 [1 2 , 3] 5 4 [3, 1 2 ] 5 2 [1, 2, 1 2 ] 5 4 [3, 2] 5 2 [1, 2 2 ] 5 4 [1, 3, 1] 5 4 [4] 4 0 [1, 4] 5 4 [4, 1] 5 1 [5] 5 1
For example, for N = 5, we construct one branch of a tree with ω 1 = 1 (see Fig. 3) : each vertex is an ordered orbit, the modulus of the gop is written on the last edge. However, the sequence of ordered gop differs from the natural downward lecture of the tree and has to be done following the algorithm.

1 1[1] 1 2[1, 1] 1 3[1, 1, 1] 1 4[1, 1, 1, 1] 1 5[1, 1, 1, 1, 1] 2 5[1, 1, 1, 2] 2 4[1, 1, 2] 1 5[1, 1, 2, 1] 3 5[1, 1, 3] 2 3[1, 2] 1 4[1, 2, 1] 2 5[1, 2, 1, 1]

CARDINAL OF CLASSES

In this section we emphasize some closed formulas giving the cardinal of classes of gop. Recalling first the already known formula for the class [1 k ] N for which we give a detailed proof, we consider the case were the class possesses exactly one k-cycle, the case with only two cycles belonging to the class and finally the main general formula of any cycles with any length. We give rigorous proof of all. The general formula is very interesting in the sense that even using computer network it is impossible to check every function of F N when N is larger than 100.

Discrete maps with 1-cycle only

The theorem 4 gives the number of discrete maps of F N which have only fixed points and no cycles of length greater than one. This formula is explicit in [START_REF] Knuth | The Art Of Computer Programming[END_REF] and [START_REF] Purdom | Cycle length in a random function[END_REF]. A complete proof is given here in detail.

Theorem 4 Let k be an integer between 1 and N. The number of functions whose

global orbit pattern is [1 k ] N (i.e. belonging to the class [1 k ] N ) is N -1 N -k N N-k . That is ♯[1 k ] N = N -1 N -k N N-k . ( 11
)

Proof

Let k be a non-zero integer. Let f be a function of F N . There are N k possibilities to choose k fixed points. There remain Nk points. Let p be an integer between 1 and Nk. We assume that p points are directly connected to the k fixed points. For each of them, there are k manners to choose one fixed point. There are k p ways to connect directly p points to k fixed points. There remains Nkp points that we must connect to the p points. There are ♯[1 p ] N-k functions. Finally, the number of functions with k fixed points is

N k N-k ∑ p=1 k p ♯[1 p ] N-k
. We now prove recursively on N

for every 0 ≤ k ≤ N that ♯[1 k ] N = N -1 N -k N N-k . We have ♯[1] 1 = 1. The formula is true. We suppose that ∀k ≤ N ♯[1 k ] N = N -1 N -k N N-k .
Let X be a set with N + 1 elements. We look for the functions of F N+1 which have k fixed points. Thanks to the previous reasoning, we have

♯[1 k ] N+1 = N + 1 k N+1-k ∑ p=1 k p ♯[1 p ] N+1-k . ♯[1 k ] N+1 = N + 1 k N+1-k ∑ p=1 k p ♯[1 p ] N-(k-1) .
We use the recursion assumption.

♯[1 k ] N+1 = N + 1 k N+1-k ∑ p=1 k p N -k p -1 (N -k + 1) N-k+1-p . ♯[1 k ] N+1 = k N + 1 k N-k ∑ p=0 N -k p k p (N -k + 1) N-k-p . ♯[1 k ] N+1 = k N + 1 k (N + 1) N-k . ♯[1 k ] N+1 = N k -1 (N + 1) N-k+1 . ♯[1 k ] N+1 = N N -k + 1 (N + 1) N-k+1 . q.e.d.

Discrete maps with k-cycle

We look now for the number of functions with exactly one k-cycle.

Theorem 5 Let k be an integer between 1 and N. The number of functions whose

global orbit pattern is [k] N is ♯[1 k ] N × (k -1)!. i.e. ♯[k] N = ♯[1 k ] N × (k -1)!. ( 12
)

Proof

There are N k ways of choosing k elements among N. Then, there are (k -1)! choices for the image of those k elements in order to constitute a k-cycle by f . We must now count the number of ways of connecting directly or not the remaining Nk elements to the k-cycle. We established already this number which is equal to

N-k ∑ p=1 k p ♯[1 p ] N-k . Finally, we have ♯[k] N = (k -1)! N k N-k ∑ p=1 k p ♯[1 p ] N-k . That is, ♯[k] N = ♯[1 k ] N × (k -1)!. q.e.d.

Discrete maps with only two cycles

We give the number of functions with only two cycles.

Theorem 6 Let N ≥ 2. Let p and q be two non-zero integers such that p + q ≤ N. Then,

♯[p, q] N = ♯[1 p+q ] N (p + q -1)! q = (N -1)! N N-(p+q) (N -(p + q))! q . ( 13
)
Proof We consider a function f which belongs to the class [1 p+q ] N . We search the number of functions constructed from f whose gop is [p, q] N . From the p fixed points of f , we construct a p-cycle. Thus, there are p+q-1 p-1 ways to choose p -1 integers among the p + q -1 fixed points. Counting the first given fixed point of f , we have p points which allow to construct (p -1)! functions with a p-cycle. Then there remain q points which give (q -1)! different functions with a q-cycle. Finally, the number of functions whose gop is [p, q] N is : p+q-1 p-1 (p -1)!(q -1)! that is the formula (p+q-1)! q .

Remark

We notice that for all k non-zero integer such that k

≤ N -1, ♯[k, 1] N = ♯[k + 1] N .

General case : discrete maps with cycles of any length

We introduce now the main theorem of the section which gives the number of gop of discrete maps thanks to a closed formula.

Given a global orbit pattern α, the next theorem gives a formula which gives the number of functions which belong to α. Theorem 7 Let p ≥ 2 be an integer. Let [ω 1 , . . . , ω p ] N be a gop of G (F N ). Then,

♯[ω 1 , . . . , ω p ] N = ♯[1 ω 1 +...+ω p ] N (ω 1 + . . . + ω p -1)! ω p × (ω p-1 + ω p ) × . . . × (ω 2 + . . . + ω p ) (14) ♯[ω 1 , . . ., ω p ] N = (N -1)! N N-(ω 1 +...+ω p ) (N -(ω 1 + . . . + ω p ))! p ∏ k=2 ( p ∑ j=k ω j ) (15)

Proof

We consider a function f which belongs to [1 ω 1 +...+ω p ] N . We search the number of functions constructed from f whose gop is [ω 1 , . . . , ω p ] N . From the ω 1 fixed points of f , we construct a ω 1 -cycle. Thus, there are ω 1 +...+ω p -1 ω 1 -1 ways to choose ω 1 -1 integers among the ω 1 + . . . + ω p -1 fixed points. Counting the first given fixed point of f , we have ω 1 points which allow to construct (ω 1 -1)! functions with a ω 1 -cycle. Then, the first fixed point of f which has not be chosen for the ω 1 -cycle, will belong to the ω 2 -cycle. Thus, there are ω 2 +...+ω p -1 ω 2 -1 ways to choose ω 2 -1 integers among the ω 2 + . . . + ω p -1 fixed points. So we have ω 2 points which allow to construct (ω 2 -1)! functions with a ω 2 -cycle. We keep going on that way until there remain ω p fixed points which allow to construct (ω p -1)! functions with a ω p -cycle. Finally, we have constructed :

ω 1 +...+ω p -1 ω 1 -1 (ω 1 -1)! ω 2 +...+ω p -1 ω 2 -1 (ω 2 -1)! × . . . × ω p-1 +ω p -1 ω p-1 -1
(ω p-1 -1)!(ω p -1)! functions. We simplify and obtain the formula.

Corollary 1 Let p be a non-zero integer. Let [ω 1 , . . ., ω p ] N be a gop of G (F N ). We suppose that there exists j such that ω j ≥ 2. Let h be an integer between 1 and ω j -1. Then

♯[ω 1 , . . . , ω j , . . ., ω p ] N = ♯[ω 1 , . . . , ω j -h, h, ω j+1 , . . . , ω p ] N × (h + ω j+1 + . . . + ω p ). ( 16
) Proof ♯[ω 1 , . . . , ω j -h, h, ω j+1 , . . . , ω p ] N × (h + ω j+1 + . . . + ω p ) = ♯[1 ω 1 +...+ω p ] N × (ω 1 +...+ω p -1)!(h+ω j+1 +...+ω p )
ω p (ω p-1 +ω p )...(ω j+1 +...+ω p )(h+ω j+1 +...+ω p )(ω j +ω j+1 +...+ω p )×...×(ω 2 +...+ω p ) .

We simplify and we exactly obtain

♯[ω 1 , . . . , ω j -h, h, ω j+1 , . . ., ω p ] N × (h + ω j+1 + . . . + ω p ) = ♯[ω 1 , . . . , ω j , . . ., ω p ] N .
Examples :

♯[2 2 , 1, 3] 11 = 11, 180, 400. ♯[5, 2, 10, 8, 15, 2, 3

FUNCTIONS WITH LOCAL PROPERTIES

Locally rigid functions

Obviously it is not possible to transpose to the functions on finite sets the notions of continuity and derivability which play a dramatic role in mathematical analysis since several centuries. In fact the class C 0 (I) of the continuous functions on the real interval I is a very small subset of the set I R of all the functions on I. Hence by analogy to this fact and trying to mimic some others properties of continuous functions, we introduce some subsets of particular functions of F N , which have local properties such as locally bounded range in a sense we precise further. Limiting the range of the function in a neighbourhood of any point of the interval induces a kind of "rigidity" of the function, hence we call these functions locally rigid functions. In these subsets, the gop are found to be fully efficient in order to describe very precisely the dynamics of the orbits. We first consider the very simple subset L R 1,N of functions for which the difference between f (p) and f (p + 1) is drastically bounded. In next subsection we consider more sophisticated subsets.

We consider the set :

L R 1,N ={ f ∈F N such that ∀p, 0 ≤ p ≤ N -2, | f (p) -f (p + 1)| ≤ 1}. Orbits of L R 1,N Theorem 8 If f ∈ L R 1,N
then f has only periodic orbits of order 1 or 2.

Proof

We suppose that f ∈ L R 1,N has a 3-cycle. We denote (a; f (a); f 2 (a)) taking a the smallest value of the 3-cycle. If a < f (a) < f 2 (a) then there exist two non-zero integers e and e ′ such that f (a) = a + e and f 2 (a) = f (a) + e ′ . Thus, We can prove in the same way that the function f can't have either 3-cycle or greater order cycle than 3.

f 2 (a) -e ′ ≤ f 3 (a) ≤ f 2 (a) + e ′ . That is f (a) ≤ a ≤ f (a) + 2e ′ .

Numerical results and conjectures

We have done numerical studies of the G (L R 1,N ) for N = 1 to 16, using the brute force of a desktop computer (i.e. checking every function belonging to these sets).

The Tables 11,12, 13, 14, 15 and 16 show the sequences for L R 1,1 to L R 1,16 .

In theses Tables we display in the first column all the gop of G (L R 1,N ) for every value of N. For a given N, there are two columns; the left one displays the cardinal of every existing class of gop (-stands for non existing gop). Instead the second shows more regularity, displaying on the row of the gop Then we are able to formulate some statements which have not yet been proved.

Statement 1 ♯[1 k ] L R 1,N = ♯[1 k+1 ] L R 1,N+1 for k ≤ N + 1 2 . ( 17
)
f ∈ L R 1,1 , f ∈ L R 1,2 , f ∈ L R 1,3 , f ∈ L R
] - + - + 1 + [2, 1, 2] - + - + - + [2 3 ] - + 1 1 4 4 [2 3 , 1] - + - + 1 1 Statement 2 ♯[2 k ] L R 1,N = ♯[2 k+1 ] L R 1,N+2 for k ≤ N 2 . ( 18
) Statement 3 ♯[2 k ] L R 1,N = ♯[2 k , 1] L R 1,N+1 for 2k ≤ N ≤ 3k -1. ( 19
)
+ - + 1 + [2, 1, 2 2 ] - + - + - + [2 2 , 1, 2] - + - + - + [2 4 ] 1 1 4 4 18 18 [2 4 , 1] - + 1 1 4 4 [1, 2 4 ] - + - + - + [2, 1, 2 3 ] - + - + - + [2 2 , 1, 2 2 ] - + - + - + [2 3 , 1, 2] - + - + - + [2 5 ] - + - + 1 1 Statement 4 ♯[2 k ] L R 1,N = k+1 ∑ i=1 ♯[2, 2, . . .. . . , 1 i th , . . . , 2 k+1 orders ] L R 1,N for 2k + 1 ≤ N = k+1 ∑ i=1 ♯[2 i-1 , 1, 2 k-i+1 ] L R 1,N for 2k + 1 ≤ N ( 20
)
] - + - + - + [2 6 ] - + 1 1 4 4 [2 6 , 1] - + - + 1 1 Statement ♯[1 N-k+1 ] L R 1,N =    1 if k = 1 2 if k = 2 4 27 (k + 1) × 3 k for 3 ≤ k ≤ N+1 2 ( 21
)
[2 7 ] 1 1 4 4 [2 7 , 1] - + 1 1 Remark We call u k = ♯[1 N-k+1 ] L R 1,N .
For k > 2, then u k is the sequence A120926 On-line Encyclopedia of integer Sequences : it is the number of sequences where 0 is isolated in ternary words of length N written with {0, 1, 2}.

These statements show that first the set L R 1,N is an interesting set to be considered for dynamical systems and secondly the gop are fruitful in this study. However the set

L R 2,N ={ f ∈F N such that ∀p, 0 ≤ p ≤ N -2, | f (p) -f (p + 1)| ≤ 2}
is too much large to give comparable results. Then we introduce more sophisticated sets we call sets with locally bounded range which more or less correspond to an analogue of the discrete convolution product of the local variation of f with a compact support function -→ α t .

Orbits and patterns of locally rigid function sets

Consider now the set :

L R -→ α t ,q,N ={ f ∈F N such that ∀p, 0 ≤ p ≤ N -r -1, r=t ∑ r=1 α r | f (p) -f (p + r)| ≤ q} { f ∈ F N such that ∀p,t ≤ p ≤ N -1, r=t ∑ r=1 α r | f (p) -f (p -r)| ≤ q} for the vector -→ α t = (α 1 , α 2 , . . ., α t ) ∈ N t , for q ∈ N.
TABLE 17. Numerical study of the set L R α t ,q,N for N = 10, t = 5, α 1 = 20, α 2 = 9, α 3 = 5, α 4 = 2 and α 5 = 1, for q = 20, .. The functions belonging to these sets show a kind of "rigidity": the less is q, the more "rigid" is the function, this "rigidity" being modulated by the vector -→ α t . Furthermore, the maximal length of a periodic orbit increases with q, and so the number of gop ♯G (L R -→ α t ,q,N ) and the maximal modulus of the gop.

Remark

Using this generalized notation, one has : L R 1,n = L R 1,1,n and L R 2,n = L R 1,2,n .

As an example, we explore numerically the case : N = 10, t = 5, α 1 = 20, α 2 = 9, α 3 = 5, α 4 = 2 and α 5 = 1, for q = 20, . . ., 142. The results are displayed in Table 17. In this Table "modulus" means the maximal modulus of the gop belonging to this set for the corresponding value of q in the row, "gop number" stands for ♯G (L R -→ α t ,q,N ) and "functions number" for ♯L R -→ α t ,q,N . One can point out that for the particular function -→ α t of the example; it is possible to find 10 intervals I 1 , I 2 , . . . , I 10 ⊂ N such that if q ∈ I r then there is no periodic orbit whose period is strictly greater than r, (e.g., I 6 = [[74, 87]]). Furthermore it is possible to split these intervals into subintervals I r,s in which ♯G L R -→ α t ,q,N is constant when q thumbs I r,s . This is not the case for ♯L R -→ α t ,q,N .

CONCLUSION

A discrete dynamical system associated to a function on finite ordered set X can only exhibit periodic orbits. However the number of the periods and the length of each are not easily predictable. We formalise such a gop as the ordered set of periods when the initial value thumbs X in the increasing order. We can predict by means of closed formulas, the number of gop of the set of all the function from X to itself. We also explore, using the brute force of computers, some subsets of locally rigid functions on X , for which interesting patterns of periodic orbits are found. Further study is needed to understand the behaviour of dynamical systems associated to functions belonging to these sets.

FIGURE 1 .

 1 FIGURE 1. Graph of the map f (x) = 4x(1x) on [0, 1]

FIGURE 2 .

 2 FIGURE 2. Invariant measure of the logistic map

Definition 5

 5 For example, for N = 11, the class of the gop [2, 2, 1, 3] 11 contains the following few of many functions defined in Tables 5, 6 and 7. The periodic orbit which are encountered have the same length nevertheless there are different. Let be A = [ω 1 , . . ., ω p ] N a gop.

FIGURE 3 .

 3 FIGURE 3. Branch of the tree for the construction of the gop with ω 1 = 1 on G (F 5 )

  And finally we have the relation a + e ≤ a which is impossible. If a < f 2 (a) < f (a) then there exist two non-zero integers e and e ′ such that f 2 (a) = a + e and f (a) = f 2 (a) + e ′ . Thus, f (a)e ≤ f 3 (a) ≤ f (a) + e. That is f (a)e ≤ a ≤ f (a) + e. But f (a)e = a + e ′ . And finally we have the relation a + e ′ ≤ a which is impossible.

  [2 k ] the sum of the cardinals of all the classes of the gop of the form [2, 2, . . .. . . , 1 i th , . . ., 2 k+1 orders ] which exist.

TABLE 1 .

 1 Coexisting

		periodic orbits found using 1,000 random
	initial points for double precision numbers
	Period	Orbit	Relative Basin size
	1 {0} (unstable fixed point) 15,784,521 Scattered over the interval	596 over 1,000 404 over 1,000

TABLE 2 .

 2 Coexisting periodic orbits for the discretization with regular meshes of N = 9, 10 and 11 points

	N Period	Orbit	Basin size
	9 9 9 10 10 11 11	1 1 1 1 2 1 4 {3, 8, 6, 9} {0} {6} {3, 7} {0} {3, 8} {0}	3 over 9 2 over 9 4 over 9 2 over 10 8 over 10 3 over 11 8 over 11

TABLE 3 .

 3 Coexisting periodic orbits for the discretization with regular meshes of N = 99, 100 and 101 points

	N Period	Orbit	Basin size
	99 99 100 100 100 100 101 101 101	1 10 {3, 11, 39, 93, 18, 58, 94, 15, 50, 97} {0} 1 {0} 1 {74} 6 {11, 39, 94, 18, 58, 96} 7 {7, 26, 76, 70, 82, 56, 97} 1 {0} 1 {75} 1 {16, 61, 95}	3 over 99 96 over 99 2 over 100 2 over 100 72 over 100 24 over 100 3 over 101 2 over 101 96 over 101

TABLE 4 .

 4 Coexisting periodic orbits for the discretization with regular meshes of N = 1, 999;2, 000 and 2, 001 points

	N Period	Orbit	Basin size
	1,999 1,999 1,999 2,000 2,000 2,000 2,000 2,000 2,001 2,001 2,001 2,001 2,001 2,001	1 4 8 {3; 11; 43; 168; 615; 1, 702; 1, 008; 1, 997} 1,006 over 1,999 3 over 1,999 {0} 990 over 1,999 {554; 1, 601; 1, 272; 1, 848} 1 2 over 2,000 {0} 1 14 over 2,000 {1, 499} 2 138 over 2,000 {691; 1, 808} 3 6 over 2,000 {276; 1, 221; 1, 900} 8 {3; 11; 43; 168; 615; 1, 703; 1, 008; 1, 998} 1,840 over 2,000 1 5 over 2,001 {0} 1 34 over 2,001 {1, 500} 2 92 over 2,001 {691; 1, 809} 608 over 2,001 8 {3; 11; 43; 168; 615; 1, 703; 1, 011; 1, 999} 18 263 over 2,001 {35; 137; 510; 1, 519; 1, 461; 1, 574; . . .} 25 {27; 106; 401; 1, 282; 1, 840; 588; . . .} 1,262 over 2,001

TABLE 5 .

 5 Orbits and components of a function belonging to F 11 with gop [2, 2, 1, 3] 11 .

	Function	orbit/component/nature
	→ 6 → 3 → 2 → 5 → 8 → 10 → 9 → 4 → 7 10 → 5 → 6	period-2 orbit : {6, 9} period-2 orbit : {5, 10} {1, 3, 5, 10} attractive attractive {0, 6, 9} fixed point : {2} repulsive {2} period-3 orbit : {4, 8, 7} {4, 8, 7} repulsive

TABLE 6 .

 6 Orbits

	Function	orbit/component/nature
	→ 4 → 2 → 7 → 3 → 8 → 10 → 5 10 → 6 → 1 → 4 → 2	period-2 orbit : {4, 8} period-2 orbit : {2, 7} fixed point : {3}	attractive {1, 2, 7, 9} attractive {0, 4, 8} {3} repulsive

and components of a function belonging to F 11 with gop [2, 2, 1, 3] 11 . period-3 orbit : {5, 10, 6} {5, 10, 6} repulsive

TABLE 7 .

 7 Orbits and components of a function belonging to F 11 with gop [2, 2, 1, 3] 11 .

	Function	orbit/component/nature
	0 → 9 1 → 6 2 → 4 3 → 7 4 → 10 5 → 3 6 → 1 7 → 5 8 → 2 10 → 10 9 → 0	period-2 orbit : {0, 9} period-2 orbit : {1, 6} fixed point : {10} period-3 orbit : {3, 7, 5}	repulsive repulsive {2, 4, 8, 10} attractive {0, 9} {1, 6} {3, 7, 5} repulsive

TABLE 8 .

 8 Orbits and components of a function belonging to F 8 with gop [2, 1, 4] 8 . The set of all the global orbit patterns of

	Function	orbit/component/nature
	0 → 1 1 → 0 2 → 0 3 → 3 4 → 5 5 → 6 7 → 4 6 → 7	period-2 orbit : {0, 1} fixed point : {3} period-4 orbit : {4, 5, 6, 7} {4, 5, 6, 7} repulsive attractive {0, 1, 2} repulsive {3}
	Definition 3

TABLE 9 .

 9 Algorithm for the threshold function construction for the gop [2 2 , 1, 3] 11 .

	First step	Second step	Third step	Fourth step	Fifth step
	Construction of
	the first canoni-

TABLE 10 .

 10 Ordered gop for N = 5 with modulus and modulus-ω 1

  ] 50 = 29, 775, 702, 147, 667, 389, 218, 762, 343, 520, 975, 006, 348, 329, 578, 044, 480, 000, 000, 000, 000, 000. ♯[5, 2, 10, 8, 15, 2, 3] 50 ∼ = 2.98 × 10 63 among the 8.88 × 10 84 functions of F 50 .

TABLE 11 .

 11 Numbering the locally rigid functions for

TABLE 12 .

 12 Numbering the locally rigid functions forf ∈ L R 1,5 , f ∈ L R 1,6 , f ∈ L R 1,7 .

			1,4 .
	g.o.p.	N=1 N=1 N=2 N=2 N=3 N=3 N=4 N=4
	Total number		1		4		17	68
	[1]	1	+	2	+	7	+	+
	[1 2 ] [1 3 ] [1 4 ] [2]	----	+ + + +	1 --1	+ + + 1	4 1 -4	+ + + 4	+ + + 18
	[2, 1]	-	+	-	+	1	1	4
	[1, 2]	-	+	-	+	-	+	+
	[2 2 ]	-	+	-	+	-	+	1
	g.o.p.		N=5 N=5 N=6 N=6	N=7 N=7
	Total number		259		950		387
	[1]		95	+ 340	+ 1,193	+
	[1 2 ] [1 3 ] [1 4 ] [1 5 ] [1 6 ] [1 7 ] [2]		50 16 4 1 --70	+ 174 + 58 + 16 + 4 + 1 + -70 264 264 + + + + + +	600 204 60 16 4 1 952 952 + + + + + +
	[2, 1]		12	18	45	70	166 264
	[1, 2]		6	+	25	+	98	+
	[2 2 ] [2 2 , 1] [1, 2 2		4 1	4 1	18 4	18 4	70 17	70 18

TABLE 13 .

 13 Numbering the locally rigid functions forf ∈ L R 1,8 , f ∈ L R 1,9 , f ∈ L R 1,10 .

	g.o.p.	N=8	N=8	N=9	N=9	N=10	N=10
	Total number		11,814		40,503		13,6946
	[1]	4,116	+ 14,001	+ 47,064	+
	[1 2 ] [1 3 ] [1 4 ] [1 5 ] [1 6 ] [1 7 ] [1 8 ] [1 9 ] [1 10 ] [2]	2,038 700 214 60 16 4 1 --3,356	+ + + + + + + + + 3,356 11,580 11,580 39,364 6,852 + 22,806 2,366 + 7,896 742 + 2,520 216 + 754 60 + 216 16 + 60 4 + 16 1 + 4 -+ 1	+ + + + + + + + + 39,364
	[2, 1]	590	952	2,062	3,356	7,072	11,580
	[1, 2]	362	+	1,294	+	4,508	+
	[2 2 ] [2 2 , 1] [1, 2 2 ] [2, 1, 2]	264 62 6 2	264 70 + +	952 222 28 14	952 264 + +	3,356 770 113 69	3,356 952 + +
	[2 3 ] [2 3 , 1] [1, 2 3 ]	18 4 -	18 4	70 18	70 18	264 69	264 70

TABLE 14 .

 14 Numbering the locally rigid functions forf ∈ L R 1,11 , f ∈ L R 1,12 , f ∈ L R 1,13 .

	g.o.p.	N=11	N=11	N=12	N=12	N=13	N=13
	Total number		457,795		1,515,926		4,979,777
	[1]	156,629	+ 516,844	+ 1,693,073	+
	[1 2 ] [1 3 ] [1 4 ] [1 5 ] [1 6 ] [1 7 ] [1 8 ] [1 9 ] [1 10 ] [1 11 ] [1 12 ] [1 13 ] [2]	75,292 26,098 8,434 2,756 756 216 60 16 4 1 --132,104 132,104 438,846 + 246,762 + 85,556 + 27,904 + 8,658 + 2,590 + 756 + 216 + 60 + 16 + 4 + 1 + -	+ + + + + + + + + + + + 438,846 1,445,258 1,445,258 803,706 + 278,580 + 91,488 + 28,738 + 8,730 + 2,592 + 756 + 216 + 60 + 16 + 4 + 1 +
	[2, 1]	23,941	39,364	80,108	132,104	265,548	438,846
	[1, 2]	15,423	+	51,996	+	173,298	+
	[2 2 ] [2 2 , 1] [1, 2 2 ] [2, 1, 2]	11,580 2,634 429 293	11,580 3,356 + +	39,364 8,883 1,555 1,142	39,364 11,580 + +	132,104 29,659 5,478 4,227	132,104 39,364 + +
	[2 3 ] [2 3 , 1] [1, 2 3 ] [2, 1, ]	952 255 7 2	952 264 + +	3,356 899 35 16	3,356 952 + +	11,580 3,098 152 86	11,580 3,356 + +
	[2 2 , 1, 2] [2 4 ] [2 4 , 1] [1, 2 4 ] [2, 1, ]	-70 18 --	+ 70 18 + +	2 264 70 --	+ 264 70 + +	20 952 263 1 -	+ 952 264 + +
	[2 2 , 1, 2 ] [2 3 , 1, 2] [2 5 ] [2 5 , 1] [1, 2 5	--4 1	+ + 4 1	--18 4	+ + 18 4	--70 18	+ + 70 18

TABLE 15 .

 15 Numbering the locally rigid functions for f ∈ L R 1,14 , f ∈ L R 1,15 .

	g.o.p.	N=14	N=14	N=15	N=15
	Total number		16,246,924		52,694,573
	[1]	5,511,218	+ 17,841,247	+
	[1 2 ] [1 3 ] [1 4 ] [1 5 ] [1 6 ] [1 7 ] [1 8 ] [1 9 ] [1 10 ] [1 11 ] [1 12 ] [1 13 ] [1 14 ] [1 15 ] [2]	2,603,258 901,802 297,728 94,440 29,050 8,746 2,592 756 216 60 16 4 1 -4,725,220	+ + + + + + + + + + + + + + 4,725,220 15,352,392 15,352,392 8,391,360 + 2,904,592 + 962,888 + 307,848 + 95,676 + 29,140 + 8,748 + 2,592 + 756 + 216 + 60 + 16 + 4 + 1 +
	[2, 1]	873,149	1,445,258	2,851,350	+
	[1, 2]	572,109	+	1,873,870	+
	[2 2 ] [2 2 , 1] [1, 2 2 ] [2, 1, 2]	438,846 98,135 18,873 15,096	438,846 132,104 + +	1,445,258 322,310 63,967 52,569	1,445,258 438,846 + +
	[2 3 ] [2 3 , 1] [1, 2 3 ] [2, 1, ]	39,364 10,460 605 389	39,364 11,580 + +	132,104 34,845 2,282 1,596	132,104 39,364 + +
	[2 2 , 1, 2] [2 4 ] [2 4 , 1] [1, 2 4 ] [2, 1, ]	126 3,356 942 8 2	+ 3,356 952 + +	641 11,580 3,292 44 18	+ 11,580 3,356 + +
	[2 2 , 1, 2 ] [2 3 , 1, 2] [2 5 ] [2 5 , 1] [1, 2 5 ] [2 6 ] [2 6 , 1]	--264 70 -18 4	+ + 264 70 + 18 4	2 -952 264 -70 18	+ + 952 264 + 70 18

TABLE 16 .

 16 Numbering the locally rigid functions for f ∈ L R 1,16 .

	g.o.p.	N=16	N=16
	Total number		170,028,792
	[1]	57,477,542	+
	[1 2 ] [1 3 ] [1 4 ] [1 5 ] [1 6 ] [1 7 ] [1 8 ] [1 9 ] [1 10 ] [1 11 ] [1 12 ] [1 13 ] [1 14 ] [1 15 ] [1 16 ] [2]	26,932,398 9,314,088 3,097,650 996,764 312,456 96,096 29,158 8,748 2,592 756 216 60 16 4 1 49,610,818	+ + + + + + + + + + + + + + + 49,610,818
	[2, 1]	9,255,822	15,352,392
	[1, 2]	6,096,570	+
	[2 2 ] [2 2 , 1] [1, 2 2 ] [2, 1, 2]	4,725,220 1,051,686 213,975 179,597	4,725,220 1,445,258 + +
	[2 3 ] [2 3 , 1] [1, 2 3 ] [2, 1, ]	438,846 114,798 8,284 6,146	438,846 132,104 + +
	[2 2 , 1, 2] [2 4 ] [2 4 , 1] [1, 2 4 ] [2, 1, ]	2,876 39,364 11,246 204 106	+ 39,364 11,580 + +
	[2 2 , 1, 2 ] [2 3 , 1, 2] [2 5 ] [2 5 , 1] [1, 2 5 ] [2 6 ] [2 6 , 1] [2 7 ] [2 7 , 1] [2 8 ]	22 2 3,356 951 1 264 70 18 4 1	+ + 3,356 952 + 264 70 18 4 1


Cristel Chandre
email: cristel.chandre@cnrs.fr

Bradley A Shadwick
email: shadwick@unl.edu

Four-field Hamiltonian fluid closures of the one-dimensional Vlasov-Poisson equation

published or not. The documents may come

Four-field Hamiltonian fluid closures of the one-dimensional Vlasov-Poisson equation

I. INTRODUCTION

In order to simulate the dynamics of a plasma, there is a variety of models which are used according to the type of question and the level of detail in the description of the plasma. Most of these models can be categorized as kinetic or fluid, whether the dynamical field variables are functions of the phase-space coordinates (x, v) of the particles or just configuration space coordinates x. Compared to kinetic models, fluid models have the significant advantage to be defined in a dimensionally reduced space, which makes them particularly desirable from a computational viewpoint. The central question is how to define these fluid models from a parent kinetic model. There are plethora of methods to do this, some better suited than others depending on the specific problem at hand. For instance, some reductions rely on an assumption on the shape of the distribution function, [START_REF] Gosse | Two moment systems for computing multiphase semiclassical limits of the Schrödinger equation[END_REF][START_REF] Jin | Multi-phase computations of the semiclassical limit of the Schrödinger equation and related problems: Whitham vs Wigner[END_REF][START_REF] Shadwick | Hamiltonian description of lowtemperature relativistic plasmas[END_REF][START_REF] Fox | Higher-order quadrature-based moment methods for kinetic equations[END_REF][START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF] or introduce suitably designed dissipative terms. [START_REF] Hammett | Fluid moment models for Landau damping with application to the ion-temperature-gradient instability[END_REF][START_REF] Passot | A fluid description for Landau damping of dispersive MHD waves[END_REF][START_REF] Sarazin | Entropy production and collisionless fluid closure[END_REF] Here we follow a different route by requiring that the reduced fluid model preserves an important dynamical property of the parent model, namely its Hamiltonian structure. [START_REF] Perin | Higher-order Hamiltonian fluid reduction of Vlasov equation[END_REF][START_REF] Perin | Hamiltonian closures for fluid models with four moments by dimensional analysis[END_REF][START_REF] Perin | Hamiltonian fluid closures of the Vlasov-Ampère equations: From water-bags to N moment models[END_REF] Rather than being an additional constraint on the reduction, we will see that this requirement provides a way to perform the reduction, and precisely define the relevant closures leading to the definition of relevant Hamiltonian fluid model(s). In order to illustrate this point we consider the one-dimensional Vlasov-Poisson equation. This equation describes the evolution of the distribution function f (x, v, t) of charged particles (of charge e and mass m) in an electric field E(x, t):

∂f ∂t = -v ∂f ∂x - e E m ∂f ∂v ,
where E is the fluctuating part of the electric field E whose dynamics is given by ∂E ∂t = -4π j, and j = e vf dv is the current density. We assume periodic boundary conditions in x with period 2L x , so that the fluctuating part is defined as

E = E - 1 2L x Lx -Lx E dx.
We consider a fluid description obtained by using the first four fluid moments of the distribution function, more precisely, the density ρ(x, t), the fluid velocity u(x, t), the pressure P (x, t) and the heat flux q(x, t) defined by ρ = f dv, u = ρ -1 vf dv,

P = (v -u) 2 f dv, q = 1 2 (v -u) 3 f dv.
From the Vlasov-Poisson equation, we obtain the equations of motion for these moments and for the electric field:

∂ t ρ = -∂ x (ρu), (1a)
∂ t u = -u∂ x u - 1 ρ ∂ x P + e E m ( 1b
)
∂ t P = -u∂ x P -3P ∂ x u -2∂ x q, (1c)
∂ t q = -u∂ x q -4q∂ x u + 3P 2ρ ∂ x P - 1 2 ∂ x R, (1d)
∂ t E = -4πe ρu, (1e)
where R = (v -u) [START_REF] Fox | Higher-order quadrature-based moment methods for kinetic equations[END_REF] f dv, which is related to the kurtosis (in velocity) of the distribution function f . Here and in what follows, ∂ x and ∂ v denote the partial derivatives of a function of x and v with respect to x and v, respectively. In order to close the set of equations of motion, we need an equation of state of the form

R = R(ρ, u, P, q).
An example of closure is obtained by assuming a Gaussian distribution for f

(see Ref. 3), f (x, v, t) = ρ √ 2πσ 2 e -(v-u) 2 /(2σ 2 ) ,
which leads to R = 3P 2 /ρ, independent of u and q. One of the main problems of the Gaussian closure is that the resulting model breaks the original Hamiltonian structure of the parent model, the Vlasov-Poisson equation. [START_REF] De Guillebon | Hamiltonian structure of reduced fluid models for plasmas obtained from a kinetic description[END_REF] As a consequence, this closure introduces unphysical dissipation.

Based on the preservation of the Hamiltonian structure, another closure based on dimensional analysis was proposed in Ref. 10, namely

R = P 2 ρ + 4q 2 P . (2)
We notice that this closure depends explicitly on the asymmetries of the distribution function, measured by q, and is still independent of the fluid velocity u. However this closure has a fundamental drawback which is that homogeneous equilibria are all unstable. In order to see this, we linearize the equations of motion around one of such equilibria with q 0 = 0, u 0 = 0 and E 0 = 0, i.e., ρ = ρ 0 + δρ, u = δu, P = P 0 + δP , q = δq and E = δE. The linearized equations of motion for δX = (δρ, δu, δP, δq, δE) in Fourier space, i.e., for

δX = ∞ k=-∞ δX k e ikx , reduce to δX k = A δX k , (3)
where

A =           0 -ikρ 0 0 0 0 0 0 -ikρ -1 0 0 e/m 0 -3ikP 0 0 -2ik 0 ikρ -2 0 P 2 0 /2 0 ikρ -1 0 P 0 /2 0 0 0 -4πeρ 0 0 0 0           .
The matrix A does not have purely imaginary eigenvalues for

k 2 < k 2 c = 4πe 2 ρ 2 0 mP 0 ,
from which we conclude that all equilibria with q 0 = 0 are unstable.

Here we are looking for a closure which combines two important properties of the Vlasov-Poisson equation, namely, the stability of symmetric homogeneous equilibria, and its Hamiltonian structure.

We do not assume any particular form for the distribution function. Instead we solve the Jacobi identity in order to determine all possible R(ρ, u, P, q) for which this identity is satisfied. As a result, we unveil a one-parameter family of Hamiltonian fluid closures. We show that for these closures, the associated Poisson bracket has two Casimir invariants of the entropy type, i.e., two observables C of the form C = dx ρ Γ(ρ, P, q). These Casimir invariants provide normal variables in which the closure in parametric form is found to be polynomial. We then examine numerically some properties of the resulting Hamiltonian model in two cases: plasma oscillations and the two-stream instability.

II. DERIVATION OF THE FOUR-FIELD HAMILTONIAN CLOSURE

The one-dimensional Vlasov-Poisson equation has a Hamiltonian structure [START_REF] Morrison | Poisson brackets for fluids and plasmas[END_REF] (see also

Refs. 14 and 15 for a review), i.e., the equations of motion can be recast using a Hamiltonian and a Poisson bracket:

∂ t f = {f, H}, (4a)
∂ t E = {E, H}, (4b)
where

H[f, E] = dxdvf mv 2 2 + dx E 2 8π .
The Poisson bracket between two scalar functionals of f and E is given by

{F, G} = 1 m f ∂ x δF δf ∂ v δG δf -∂ v δF δf ∂ x δG δf -4πe δF δE ∂ v δG δf -∂ v δF δf δG δE dxdv, (5)
where δF δf and δF δE denote the functional derivatives of F with respect to f and E respectively. In particular, this bracket satisfies the Jacobi identity, i.e., {F, {G, H}} + {H, {F, G}} + {G, {H, F }} = 0, for all observables F , G and H. For simplicity of the notations and without loss of generality, we assume that m = 1.

Remark: Gauss's law is derived from a Casimir invariant of the bracket (5):

C[f, E] = ∂ x E -4πe dvf.
Here we consider a neutral plasma, i.e., such that the value of this Casimir invariant is -4πe which expresses the presence of a neutralizing background.

Regardless of the truncation, Eq. ( 1) can be recast in the following form (see Ref. 10 for more details):

∂ t X = {X, H},
where X = (ρ, u, S 2 = P/ρ 3 , S 3 = 2q/ρ

∂ x S i + α ij 1 ρ 2 δF δS i δG δS j + ∂ x 1 ρ δF δS i β ij 1 ρ δG δS j .( 6
)
The 2 × 2 matrices α = ∂ x γ and β are given by

γ =   2S 3 2S 4 -3S 2 2 3S 4 -6S 2 2 3S 5 -12S 2 S 3   , and
β =   4S 3 5S 4 -9S 2 2 5S 4 -9S 2 2 6S 5 -24S 2 S 3   ,
where S 4 = R/ρ 5 and S 5 is an arbitrary function of ρ, u, S 2 and S 3 . As a consequence, since the bracket is antisymmetric, the models are all conserving energy regardless of the closure S 4 = S 4 (ρ, u, S 2 , S 3 ) and S 5 = S 5 (ρ, u, S 2 , S 3 ). We notice that β = γ + γ T . This allows us to rewrite the Poisson bracket in a more antisymmetric way

{F, G} = dx δG δu ∂ x δF δρ - δF δu ∂ x δG δρ -4πe δG δu δF δE - δF δu δG δE - 1 ρ ∂ x S i δG δu δF δS i - δF δu δG δS i + 1 ρ δG δS i γ ij ∂ x 1 ρ δF δS j - 1 ρ δF δS i γ ij ∂ x 1 ρ δG δS j .
The Jacobi identity for the above bracket leads to the following constraints on the matrix γ:

(γ kn + γ nk ) ∂γ ij ∂S n = (γ jn + γ nj ) ∂γ ik ∂S n , ( 7a
)
∂γ in ∂S m ∂γ jk ∂S n = ∂γ jn ∂S m ∂γ ik ∂S n , (7b)
for all i, j, k, m (and repeated summation over n).

A. Explicit expression for the Hamiltonian closure

In Ref. 10, it was shown that in order for the bracket (6) to be Hamiltonian, the closures S 4 and S 5 needs to be of the form S 4 = S 4 (S 2 , S 3 ) and S 5 = S 5 (S 2 , S 3 ), i.e., they do not depend on ρ and u. The conditions (7) boil down to three constraints

6S 5 = 12S 2 S 3 + 4S 3 ∂S 4 ∂S 2 + (5S 4 -9S 2 2 ) ∂S 4 ∂S 3 , ∂S 5 ∂S 2 = 4S 3 + ∂S 4 ∂S 3 ∂S 4 ∂S 2 -3S 2 , ∂S 5 ∂S 3 = ∂S 4 ∂S 2 + ∂S 4 ∂S 3 2 .
Equivalently, a necessary and sufficient condition is that the closure function S 4 satisfies the following two coupled nonlinear partial differential equations

4S 3 ∂ 2 S 4 ∂S 2 2 -(9S 2 2 -5S 4 ) ∂ 2 S 4 ∂S 2 ∂S 3 - ∂S 4 ∂S 2 ∂S 4 ∂S 3 -12S 3 = 0, (8a)
4S 3 ∂ 2 S 4 ∂S 2 ∂S 3 -(9S 2 2 -5S 4 ) ∂ 2 S 4 ∂S 2 3 - ∂S 4 ∂S 3 2 -2 ∂S 4 ∂S 2 + 12S 2 = 0. (8b)
From these equations, we readily check that the Gaussian closure S 4 = 3S 2 2 is not a solution of these equations, which means that the Gaussian closure is not Hamiltonian. In addition, we check that the solution given by Eq. ( 2), corresponding to the dimensional analysis of Ref. 10, i.e., S 4 = S 2 2 + S 2 3 /S 2 , is the simplest solution. However, this is not an adequate solution since all homogeneous equilibria are always found to be unstable, as pointed out above. To solve Eqs. ( 8), we start by looking for solutions close to symmetric distributions, i.e.,

S 4 (S 2 , S 3 ) = f 0 (S 2 ) + S 2 3 f 1 (S 2 ) + O(S 4
3 ).

We insert this expansion in Eqs. ( 8) and consider their leading behavior near S 3 = 0. This lead to a set of two coupled ordinary differential equations

2f 0 -(9S 2 2 -5f 0 )f 1 -f 1 f 0 -6 = 0, -(9S 2 2 -5f 0 )f 1 -f 0 + 6S 2 = 0.
By combining these two equations, we obtain one single ordinary differential equation

f 0 (9S 2 2 -5f 0 ) + 2f 2 0 -18S 2 f 0 + 20f 0 = 0.
Near S 2 = 0, we look for solutions of the type

f 0 (S 2 ) = kS α 2 .
A possible solution is obviously the one obtained using the dimensional analysis 10 , i.e.,

f 0 (S 2 ) = S 2 2 .
In addition there is a less trivial family of solutions for α = 5/3. More generally, we look at solutions which can be expanded in Puiseux series

f 0 (S 2 ) = ∞ n=5 a n S n/3 2 .
We show that the only possible solutions are f 0 (S 2 ) = S 2 2 and

f 0 (S 2 ) = kS 5/3 2 ,
for any value of k. For practical purposes, we define κ = 5k/9. We notice that contrary to the solution provided by dimensional analysis, the second solution comes as a family parameterized by κ. The interesting feature is that this family extends to a Hamiltonian closure for arbitrary large values of S 3 . Indeed we are looking for a solution which can be expanded as

S 4 (S 2 , S 3 ) = ∞ n=0 f n (S 2 )S 2n 3 . (9)
Inserting this ansatz in Eq. ( 8) leads to a recurrence relation for the coefficients f n (S 2 ):

f 0 (S 2 ) = 9κ 5 S 5/3 2 , (10a)
f 1 (S 2 ) = κ -2S 1/3 2 3S 2 (κ -S 1/3 2 ) , (10b)
f n+1 = - 1 9(n + 1)(2n + 1)S 5/3 2 (κ -S 1/3 2 ) (4n -1)f n + n m=1 m(12m -7 -2n)f m f n+1-m , (10c)
and an addition constraint where f n has to satisfy

S 5/3 2 (κ -S 1/3 2 )(n + 1)f n+1 -(n + 1) κ 3 S 2/3 2 f n+1 + 2 9 f n + 1 9 n m=1 (6m -n -1)f m f n+1-m = 0, (11)
for all n ≥ 1. The first few terms are given by

f 2 (S 2 ) = 1 3 4 S 11/3 2 (κ -S 1/3 2 ) , f 3 (S 2 ) = 2 5κ -3S 1/3 2 3 8 S 19/3 2 (κ -S 1/3 2 ) 3 , f 4 (S 2 ) = 48κ 2 -61κS 1/3 2 + 18S 2/3 2 3 11 S 9 2 (κ -S 1/3 2 ) 5
.

The expression of other terms of the series expansion of S 4 can be obtained using a MAT-LAB [START_REF]MATLAB, Version 9.12.0[END_REF] code available at Ref. 17. We are not able to prove directly that for all n, the f n s obtained by the recursion relation (10) satisfy Eq. (11). However, we have checked that for n below 25, these conditions are satisfied using symbolic computations available from the MATLAB [START_REF]MATLAB, Version 9.12.0[END_REF] code. Beyond this value of 25, the symbolic computations are too complex to allow simplifications in a reasonable amount of time. By truncating the series (9), i.e., by considering

S 4 (S 2 , S 3 ) = nmax n=0 f n (S 2 )S 2n 3 ,
we have found that the Jacobi identity is satisfied up to orders S 2nmax 3 for the values of n max we have tested. This led us to conjecture that the limit n max → ∞ corresponds to a Hamiltonian closure. We notice that the closure is singular at

S (c) 2 = κ 3 , so this explicit closure S 4 = S 4 (S 2 , S 3 ) is valid only in the range S 2 ∈ [0, S (c) 2 [
. Remark: Scaling. We notice that the functions f n satisfy

f n (λ 2 S 2 ; λ 2/3 κ) = λ 4-6n f n (S 2 ; κ),
for all n ≥ 0. Therefore, we have a scaling relationship for S 4 :

S 4 (λ 2 S 2 , λ 3 S 3 ; λ 2/3 κ) = λ 4 S 4 (S 2 , S 3 ; κ).
A contour plot of S 4 in the plane (S 2 , S 3 ) is represented in Fig. 1 In particular one interesting feature is that the first order of the closure does not depend on ρ, i.e., R(ρ, P, q = 0) = 9κ 5 P 5/3 .

Remark: Relation between the kurtosis and the skewness.

B. Casimir invariants

A very interesting property of the noncanonical Poisson bracket (6) is that it possesses a number of Casimir invariants, i.e., observables C such that {C, F } = 0 for any other observable F . First we are looking for Casimir invariants of the entropy type, i.e., C = dxρΓ(S 2 , S 3 ).

The function Γ satisfies the following conditions:

β ij ∂ 2 Γ ∂S i ∂S n + ∂γ ij ∂S n ∂Γ ∂S i = 0, (12)
for all j, k and n in (2, 3) (and where we assumed implicit summation over the repeated index i). We assume that we have K solutions, denoted Γ k for k = 2, . . . , K. Using the property β = γ + γ T , we prove that the above-conditions are equivalent to

∂ ∂S n ∂Γ k ∂S i β ij ∂Γ l ∂S j = 0, (13)
for all n, k and l.

Using series expansions, we found two solutions to Eq. ( 13):

C 2 = ∞ n=0 dxρg n (S 2 )S 2n 3 , C 3 = ∞ n=0 dxρh n (S 2 )S 2n+1 3 ,
where the first elements in the series are:

g 0 (S 2 ) = S 1/3 2 , g 1 (S 2 ) = - 1
3 3 S 7/3 2 (κ -S 1/3 2 ) , g 2 (S 2 ) = -4κ + 3S 1/3 2 3 6 S 5 2 (κ -S 1/3 2 ) 3 , h 0 (S 2 ) = 1 3S 2 (κ -S 1/3 2 ) , h 1 (S 2 ) = 2 3 4 S 11/3 2 (κ -S 1/3 2 ) 2 , h 2 (S 2 ) = 2(5κ -4S 1/3 2 ) 3 7 S 19/3 2 (κ -S 1/3 2 ) 4 .
The functions h n and g n for n ≥ 1 are determined from the recurrence relations:

g n+1 = - 1 9(n + 1)(2n + 1)S 5/3 2 (κ -S 1/3 2 ) (4n + 1)g n + n m=0 m(6n + 4m + 1)f n+1-m g m , h n+1 = - 1 9(n + 1)(2n + 3)S 5/3 2 (κ -S 1/ 3
2 )

(4n + 3)h n + n m=0 (2m + 3n + 3)(2m + 1)f n+1-m h m ,
which are both obtained from Eq. ( 12) with j = 3 and n = 3.

These Casimir invariants allow us to define particularly relevant variables, referred to as normal variables, in which the Hamiltonian system is greatly simplified. We perform a local change of variables: (S 2 , S 3 ) → (Γ 2 , Γ 3 ), where

Γ 2 = ∞ n=0 g n (S 2 )S 2n 3 , Γ 3 = ∞ n=0 h n (S 2 )S 2n+1 3 .
The bracket (6) becomes

{F, G} = dx δG δu ∂ x δF δρ - δF δu ∂ x δG δρ -4πe δG δu δF δE - δF δu δG δE - ∂ x Γ i ρ δG δu δF δΓ i - δF δu δG δΓ i + 1 ρ ∂ x 1 ρ δF δΓ i βij δG δΓ j , ( 14
)
where β is a symmetric matrix whose elements are

βkl = ∂Γ k ∂S i β ij ∂Γ l ∂S j ,
with an implicit summation over repeated indices. From Eq. ( 13), we deduce that the matrix β is constant. As a consequence, the bracket ( 14) always satisfies the Jacobi identity.

Therefore the existence of two Casimir invariants of the entropy type for the bracket ( 6) is sufficient to ensure that it is a Poisson bracket. Note that we use the terminology Casimir invariant also for a bracket which is a priori not of the Poisson type. Using the expressions for S 3 = 0, the matrix β takes the very simple form

β =   0 1 1 0   .
In addition, the existence of two Casimir invariants of the entropy type ensures a third Casimir invariant:

C 1 = dx u - ρ 2 Γ k ( β-1 ) kl Γ l ,
which is equal to

C 1 = dx (u -ρΓ 2 Γ 3 ) .

Its expansion is given by

C 1 = dx u -ρ ∞ n=0 k n (S 2 )S 2n+1 3 ,
where

k n = n m=0 g n-m h m ,
for n ≥ 0, and the first elements of the series are given by

k 0 (S 2 ) = 1 S 2/3 2 (κ -S 1/3 2 ) , k 1 (S 2 ) = 1 3 3 S 10/3 2 (κ -S 1/3 2 ) 2 , k 2 (S 2 ) = 4κ -3S 1/3 2 3 6 S 6 2 (κ -S 1/3 2 ) 4
.

We notice that three Casimir invariants similar to C 1 , C 2 and C 3 (but of course, different) have been found for the Hamiltonian closure obtained using the dimensional analysis (see

Ref. 10).

The advantage of working in the variables Γ i instead of the variables S i is that the closure functions S 4 and S 5 are no longer present in the Poisson bracket. They are now in the Hamiltonian through the change of variables (S 2 , S 3 ) → (Γ 2 , Γ 3 ). If we truncate the closure functions S 4 and S 5 -a natural step since these functions are given as series in S 3 -the system remains Hamiltonian in the variables Γ i whereas if these truncations are performed in the bracket in the variables S i , the system would likely loose the Hamiltonian property.

C. Parametric expression for the Hamiltonian closure

There is another significant advantage to working with normal variables Γ i : What is not fully satisfactory with the variables S i is that the closure is given as a relatively complex expansion, and consequently we were not been able to check the Jacobi identity at all orders in the expansion. The origin of this complication is due to the search for an explicit closure function S 4 (S 2 , S 3 ), not to the search of a Hamiltonian closure per se. Here instead we are looking at a parametric expression of the closure, and we consider the normal variables as parameters of the closure. More precisely, we consider an arbitrary change of coordinates from some variables Γ i to variables S i :

S 2 = S 2 (Γ 2 , Γ 3 ), S 3 = S 3 (Γ 2 , Γ 3 ),
and the closure functions are given by

S 4 = S 4 (Γ 2 , Γ 3 ), S 5 = S 5 (Γ 2 , Γ 3 ).
We start with the bracket ( 14) which is a Poisson bracket since the matrix β is constant. The question of finding Hamiltonian closures is reformulated as follows: What are the functions S i for which the bracket ( 14) expressed in the variables S i is the original bracket ( 6)? The answer is given by two sets of equations

∂S i ∂Γ k β kl ∂S i ∂Γ l = β ij , (15)
∂ 2 S i ∂Γ n ∂Γ k β kl ∂S i ∂Γ l = ∂γ ij ∂Γ n , (16)
for all i, j and n. The first set of equations ( 15) defines parametrically the functions S 3 , S 4

and S 5 :

S 3 = 1 2 ∂S 2 ∂Γ 2 ∂S 2 ∂Γ 3 , S 4 = 9 5 S 2 2 + 1 5 ∂S 2 ∂Γ 2 ∂S 3 ∂Γ 3 + 1 5
∂S 2 ∂Γ 3 ∂S 3 ∂Γ 2 , S 5 = 4S 2 S 3 + 1 3 ∂S 3 ∂Γ 2 ∂S 3 ∂Γ 3 .
Once the function S 2 is specified, all of the other functions S i are uniquely determined by the above equations. By inverting the equations Γ i = Γ i (S 2 , S 3 ) or by solving one of the constraints ( 16), we obtain the following expression for S 2 (Γ 2 , Γ 3 ):

S 2 (Γ 2 , Γ 3 ) = Γ 3 2 + Γ 2 (κ -Γ 2 )Γ 2 3 . ( 17
)
Inserting this expression in the parametric equations for S 3 , S 4 and S 5 leads to the following expressions:

S 3 (Γ 2 , Γ 3 ) = Γ 2 Γ 3 (κ -Γ 2 ) 3Γ 2 2 + (κ -2Γ 2 )Γ 2 3 , (18a)
S 4 (Γ 2 , Γ 3 ) = 9κ 5 Γ 5 2 + 6Γ 3 2 (κ -Γ 2 ) 2 Γ 2 3 + Γ 2 (κ -Γ 2 ) κ 2 -3Γ 2 (κ -Γ 2 ) Γ 4 3 , (18b)
S 5 (Γ 2 , Γ 3 ) = 9κΓ 5 2 (κ -Γ 2 )Γ 3 + 10Γ 3 2 (κ -Γ 2 ) 3 Γ 3 3 + Γ 2 (κ -Γ 2 )(κ -2Γ 2 )(κ 2 -2κΓ 2 + 2Γ 2 2 )Γ 5 3 . ( 18c
)
We notice that the closure is no given as an infinite series. In particular, the functions S n for n = 2, 3, 4, 5 are polynomials in the two variables Γ 2 and Γ 3 , and the degree in Γ 3 is n and the degree in Γ 2 is n+1. Using Mathematica 23 , we have checked that the constraints ( 16) are all satisfied. The code is available at Ref. 17. The series expansion of the explicit closure S 4 = S 4 (S 2 , S 3 ) given in Eqs. ( 9) is obtained by inverting Eqs. ( 17) and (18a), and inserting them in Eq. (18b). As a consequence, this proves the Jacobi identity for the explicit closure

S 4 = S 4 (S 2 , S 3 ).
For S 2 to be positive, a necessary and sufficient condition is that κ > Γ 2 > 0 or if Γ 2 > κ,

Γ 2 3 < Γ 2 2 /(Γ 2 -κ).
This means that S 2 can take arbitrarily large values, provided that S 3 is not too large. We notice that the point (S 2 = κ, S 3 = 0) in Fig. 1 is obtained for Γ 2 = κ regardless of the value of Γ 3 .

In Fig. 2, we have represented the closure function S 4 given parametrically by Eqs. (18) for a selected range of parameters (Γ 2 , Γ 3 ). The surface gets more complicated, with more branches, as the range of (Γ 2 , Γ 3 ) is extended (see the Mathematica code available at Ref. 17). We notice that there is a central brighter patch where there is a single value of S 4 for a given (S 2 , S 3 ). It corresponds to the explicit closure S 4 = S 4 (S 2 , S 3 ) as depicted in Fig. 1. The Mathematica code is available at Ref. 17

D. Equations of motion

The Poisson bracket (14) becomes

{F, G} = dx δG δu ∂ x δF δρ - δF δu ∂ x δG δρ -4πe δG δu δF δE - δF δu δG δE - ∂ x Γ i ρ δG δu δF δΓ i - δF δu δG δΓ i + 1 ρ δG δΓ 2 ∂ x 1 ρ δF δΓ 3 - 1 ρ δF δΓ 2 ∂ x 1 ρ δG δΓ 3 ,
and the Hamiltonian is

H[ρ, u, Γ 2 , Γ 3 , E] = 1 2 dx ρu 2 + ρ 3 S 2 (Γ 2 , Γ 3 ) + E 2 4π ,
where S 2 is given by Eq. ( 17). The equations of motion are given by Ḟ = {F, H}:

∂ t ρ = -∂ x (ρu), (19a)
∂ t u = -u∂ x u - 1 ρ ∂ x ρ 3 S 2 + e E (19b)
∂ t Γ 2 = -u∂ x Γ 2 - 1 2ρ ∂ x ρ 2 ∂S 2 ∂Γ 3 , (19c)
∂ t Γ 3 = -u∂ x Γ 3 - 1 2ρ ∂ x ρ 2 ∂S 2 ∂Γ 2 , (19d)
∂ t E = -4πe ρu. (19e)
Remark 1: In the case of an external time-dependent electric field E 0 (x, t), the closure is identical. First we need to autonomize the bracket. For the Vlasov-Poisson equation, the variables are the fields f (x, v, t) and E 1 (x, t), together with t and K (K being the canonically conjugate variable to time t), such that the total electric field is

E = E 0 + E 1 .
The Hamiltonian is

H[f, E 1 , t, K] = dxdvf v 2 2 + dx E 2 1 + 2E 1 E 0 8π + K,
and the Poisson bracket

{F, G} = f ∂ x δF δf ∂ v δG δf -∂ v δF δf ∂ x δG δf -4πe δF δE 1 ∂ v δG δf -∂ v δF δf δG δE 1 dxdv + F t G K -F K G t .
For the reduced fluid equations, the Hamiltonian becomes

H[ρ, u, Γ 2 , Γ 3 , E 1 , t, K] = 1 2 dx ρu 2 + ρ 3 S 2 + E 2 1 + 2E 1 E 0 4π + K,
and the Poisson bracket

{F, G} = dx δG δu ∂ x δF δρ - δF δu ∂ x δG δρ -4πe δG δu δF δE 1 - δF δu δG δE 1 - ∂ x Γ i ρ δG δu δF δΓ i - δF δu δG δΓ i + 1 ρ δG δΓ 2 ∂ x 1 ρ δF δΓ 3 - 1 ρ δF δΓ 2 ∂ x 1 ρ δG δΓ 3 + F t G K -F K G t .
The equations of motion consists in changing E by E 0 + E 1 in the Vlasov equation and in the momentum equation, and replacing E by E 1 in the Ampère equation.

Remark 2: By rescaling the parameters Γ 2 and Γ 3 , and by rescaling the density ρ in the following way

Γ 2 = κ Γ (r) 2 , Γ 3 = √ κ Γ (r) 3 , ρ = κ -3/2 ρ (r) ,
the equations of motion (19a)-(19d) are not longer explicitly depending on κ. The parameter κ appears only in Ampère's equation or equivalently in Gauss' law. This means that the parameter of the closure κ can be viewed as the coupling parameter between the fluid part and the electrostatic part. The parameter κ can also be removed completely from the equations of motion by rescaling the charge and the electric field as

e = κ 3/4 e (r) , E = κ -3/4 E (r) .
As a consequence, the one-parameter family of Hamiltonian closures can be seen as a unique

Hamiltonian model, and the parameter κ is now in the initial condition.

E. Stability of the symmetric and homogeneous equilibria

We have found a one-parameter family of closures which fulfill the first requirement, namely, the resulting models are Hamiltonian. The second requirement is the stability of the equilibria q 0 = 0. The linearized equations of motion reduce to Eq. ( 3) with

A =           0 -ikρ 0 0 0 0 0 0 -ikρ -1 0 0 e/m 0 -3ikP 0 0 -2ik 0 0 0 -3ik(κP 2/3 0 -ρ -1 0 P 0 )/2 0 0 0 -4πeρ 0 0 0 0          
.

From the dispersion relation, we define

ω 2 0 = ω 2 p + 3κP 2/3 0 k 2 ,
where ω p = 4πe 2 ρ 0 /m is the plasma frequency. The eigenvalues of A are all purely imaginary if

ω 2 0 > ω 2 BG ,
where ω BG (k) is the Bohm-Gross dispersion relation given by

ω 2 BG = ω 2 p + 3 P 0 ρ 0 k 2 .
The non-zero eigenvalues of A are

iω = ± i √ 2 ω 2 0 ± ω 4 0 -4ω 2 p (ω 2 0 -ω 2 BG ) 1/2 .
Therefore the homogeneous equilibria are stable for ω 2 0 > ω 2 BG , which is equivalent to requiring that S 2 < S (c) 2 or Γ 2 < κ. In terms of the parameters of the equilibrium, this means that the pressure P 0 is such that P 1/3 0 /ρ 0 < κ. A crucial factor is that the closure R(ρ, P, q = 0) does not depend on ρ, and in this case, the necessary and sufficient condition for stability is

ω 2 p + k 2 ∂R ∂P > ω 2 BG .
We recall that R(ρ, P, 0) = ρ 5 S 4 P ρ 3 , 0 . The fractional exponent 5/3 in the closure comes from the requirement that R does not depend on ρ, ensuring the stability of the equilibria. More general cases for stability would be that at q = 0 ∂R ∂P > 3P ρ , ∂R ∂ρ ≤ 0, for all ρ > 0 and P > 0. However, these conditions do not ensure that the resulting model is Hamiltonian. As expected, the requirement that the model is Hamiltonian is more stringent than requiring that homogeneous equilibria are stable.

III. NUMERICAL APPLICATIONS

The objective of this section is not to offer a detailed comparison between the numerical implementation of the Hamiltonian fluid model and the one of the parent kinetic model.

The objective is more modest since we limit ourselves to a couple of illustrations of the Hamiltonian fluid model, demonstrating the feasibility and practicality of the fluid model, which could trigger further questions of a more practical nature than the ones we consider in what follows. We consider two applications, one where the fluid model leads to stable plasma oscillations and the other one where it is unstable. In all the simulations, we consider

a domain x ∈ [-L x , L x ] and v ∈ [-L v , L v ] with L v = 10.

A. Plasma oscillations

We consider the following initial distribution function, built from a skew-normal distribution,

f (x, v, 0) = 1 √ 2π (1 -A cos k x) 1 + erf αv √ 2 e -v 2 /2 ,
with A = 10 -4 , k = λ D /2 and α = 0.1 (where λ D is the Debye length). Here the velocities are in units of the thermal velocity v th = √ k B T . Given that the equilibrium has some initial fluid velocity, the Bohm-Gross dispersion relation is becomes

ω BG = ±ω p 1 ± u k ω p + 3 2 ρ 2 S 2 k 2 ω 2 p ± 2ρ 3 S 3 k 3 ω 3 p + O k 4 ω 4 p .
This is the same dispersion relation given by the fluid and the kinetic models. For the skew-normal equilibrium,

ρ = 1, u = α 2 π(1 + α 2 ) , S 2 = 1 - 2α 2 π(1 + α 2 ) , S 3 = α 3 (1 + α 2 ) 3/2 2 π 4 π -1 .
We consider the fluid model with κ = 1. Given the initial values of S 2 and S 3 , we compute the initial values for Γ 2 and Γ 3 . We represent the values of E(x, t) in Fig. 3 obtained with the fluid and the kinetic model. We notice some qualitative similarities between the kinetic and the fluid model, such as plasma oscillations. However, as expected, the fluid model does not capture the damping of the field (clearly visible for L x /λ D = 2π), which is a purely kinetic effect. For larger values of L x , i.e., L x /λ D = 3π the damping is reduced as expected, and the agreement between the kinetic and the fluid simulations is improved.

B. Two-stream instability

Next, we consider the two-stream instability with the initial distribution For this distribution, v th = √ 3 v 0 . To simplify comparison with the existing literature, we take λ = λ D / √ 3 and v 0 to be our length and velocity scales, respectively. We set A = 10 kinetic model, i.e., a growth rate of 0.25924553 ω p (which has been corrected for the effects of the spatial grid). We notice that both models display some similar features, such as the oscillations at the beginning. Also, the slope of the higher-order modes corresponds rather well, despite the fact that these modes are higher in amplitude for the fluid model.

f (x, v, 0) = (1 -A cos kx) v 2 e -v 2 /2v 2 0 √ 2π v 3 0 . ( a
The main discrepancy between both models occur when the amplitude of the field saturates, which is when the kinetic effects are predominant, and these cannot be described by the fluid model. In addition, all wavenumbers are unstable in the Hamiltonian fluid model while only the fundamental mode is unstable in the kinetic model (the higher harmonics are driven by the fundamental mode through nonlinear couplings). For both models, the initial electric field has the same initial amplitude. Nonetheless, the amplitude of the fundamental mode is slightly larger in the fluid model compared to the kinetic model (cf. the blue curves on the left panel of Fig. 4). This is due to differences in how the initial condition projects onto the system modes in the two models. In both cases, a linear analysis produces mode amplitudes that are in excellent agreement with the numerical results.

CONCLUSIONS

We

.FIG. 1 .

 1 FIG. 1. Contourplot of S 4 given by Eq. (9) as a function of S 2 and S 3 for κ = 1. The summation is computed up to S 30 3 . The vertical red line corresponds to the location of the singularity at S 2 = S (c) 2 . The MATLAB 16 code to compute symbolically the terms of the closure and obtain numerically this figure is available at Ref. 17.

  A scaling of kurtosis (related to S 4 ) with squared skewness (related to S 3 ) for plasma density fluctuations and sea-surface temperature fluctuations was found in Refs. 18-22. Using the Hamiltonian closure, this relation is found as the first two terms of the closure, i.e., S 4 = b + aS 2 3 + O(S 4 3 ) where a and b are functions of ρ and P .

FIG. 2 .

 2 FIG. 2. Parametric representation of S 4 given by Eqs. (18) as a function of S 2 and S 3 for κ = 1.

FIG. 3 .

 3 FIG. 3. Contour plot of E(x, t): Panel (a): Hamiltonian fluid model with κ = 1 and L x /λ D = 2π. Panel (b): one-dimensional Vlasov-Poisson equation with L x /λ D = 2π. Panel (c): Hamiltonian fluid model with κ = 1 and L x /λ D = 3π. Panel (d): one-dimensional Vlasov-Poisson equation with L x /λ D = 3π.

- 6 and k λ = 1 / 2 .FIG. 4 .

 6124 FIG. 4. The left panel shows magnitude of the Fourier modes E k (t) of the electric field E(x, t) as functions of time for k λ = 1/2 (blue curves), k λ = 1 (orange curves), k λ = 3/2 (green curves) and k λ = 2 (read curves). The continuous curves are for the kinetic model, and the dashed curves are for the Hamiltonian fluid model with κ = 1.30834. The right panel shows the amplitudes for k λ = 1/2 for Hamiltonian fluid model (dashed violet), kinetic (blue) and unstable mode from linear kinetic theory (red).

  have exhibited a one-parameter family of Hamiltonian fluid models with the first four fluid moments -fluid density, fluid velocity, pressure and heat flux -as a result of the reduction of the one-dimensional Vlasov-Poisson equation. The closure involves an equation for the kurtosis in velocity of the distribution function. In the course of the reduction to a Hamiltonian fluid model, we have identified some normal variables in which the closure expressed parametrically is found to be polynomial in the normal variables. Each reduced Hamiltonian fluid model possesses three Casimir invariants, two of the entropy type and one generalized velocity. We have shown that some of these models ensures the stability of symmetric homogeneous equilibria, depending on the parameter of the closure and the initial conditions.

ACKNOWLEDGMENTS CC acknowledges useful discussions with J. Féjoz and BAS acknowledges useful discussions with Frank M. Lee. This material is based upon work supported by the National Science Foundation under Grant No. DMS-1440140 while the authors were in residence at the Mathematical Sciences Research Institute in Berkeley, California, during the Fall 2018 semester. This work has been carried out within the framework of the French Federation for Magnetic Fusion Studies (FR-FCM). BAS was supported in part by the National Science Foundation under Contract No. PHY-1535678.

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.


Svenja Hermann

Pauline Butaud

Gaël Chevallier

A new Concept of Vibration Damper Based on Friction Induced Dissipation Between Magnetic Materials

Keywords: Vibration Damping, Magnetoactive Elastomers, Assembled Structure, Magnetic Contact Pressure, Friction

This paper presents a new concept for passive vibration mitigation in sandwich structures. The test structure in the present work is composed of permanently magnetized MagnetoActive Elastomers (MAEs) and ferromagnetic sheets representing the core and the skins of the sandwich respectively. The magnetic forces which develop between the core and the skins in the contact interfaces are sufficient to maintain the structure assembled without the need of glue. The sandwich structure is subjected to harmonic vibrations around a single resonance mode with different excitation magnitudes. A nonlinear evolution of the frequency response function's (FRF) amplitude are stated. This non-linearity induces on the one hand a decrease of the resonance frequency. On the other hand, the damping evolves in a non-monotonic way according to the excitation amplitude. This is very characteristic of frictional damping.

INTRODUCTION

The capacity of assembled structures to mitigate vibrations has been the subject of intensive research over the last years [START_REF] Brake | The mechanics of jointed structures: recent research and open challenges for developing predictive models for structural dynamics[END_REF]. Nowadays, the existence of reliable characterization methods allows the optimization of existing solutions as well as the testing of new concepts. In this work we propose a new type of assembled structure in which magnetoactive elastomers are used to mitigate vibrations efficiently by friction damping. Friction damping is a passive mechanism for vibration mitigation. It occurs in contact interfaces of vibrating assembled structures. When contact surfaces in the assembly are in relative motion, the friction partially opposes the movement which leads to the attenuation effect. Depending on the ratio of tangential stress and normal stress in a contact interface, the contact state can range from total stick to partial sliding to mainly sliding [START_REF] Scheel | Challenging an experimental nonlinear modal analysis method with a new strongly friction-damped structure[END_REF]. The tightening load establishes the normal stress which affects friction damping as well as the stiffness of the structure [START_REF] Teloli | Good practices for designing and experimental testing of dynamically excited jointed structures: The Orion beam[END_REF]. Vibrations generate tangential stresses and the contact state therefore depends on the vibration magnitude [START_REF] Bouaziz | Vibration reduction of an assembly by control of the tightening load[END_REF][START_REF] Vakilinejad | A comparison of robustness and performance of linear and nonlinear Lanchester dampers[END_REF]. This leads to a variation of the resonance frequency and a nonlinear evolution of the system's response when the contact state changes from total stick to mainly sliding.

An emerging type of mechanical joint in assembled structures is based on magnetically induced forces, where electromagnets [START_REF] Dai | Structural passive control on electromagnetic friction energy dissipation device[END_REF] or permanent magnets [START_REF] Choi | Smart damper using the combination of magnetic friction and precompressed rubber springs[END_REF][START_REF] Amjadian | Analytical modeling of a simple passive electromagnetic eddy current friction damper[END_REF] have been used to generate the normal tightening load. Permanent magnets are particularly interesting for this purpose as the tightening load is maintained without an external energy supply. However, the magnets add a considerable mass to the system and existing solutions are rather bulky. In this work, we present a concept for vibration mitigation that relies on a magnetic tightening load. The latter is generated by MagnetoActive Elastomers (MAEs) in a sandwich structure. Existing studies on vibration damping with MAE sandwich structures deal with the adaptive control of the material properties of non-permanently magnetized MAEs [START_REF] Hu | Experimental investigation of the vibration characteristics of a magnetorheological elastomer sandwich beam under non-homogeneous small magnetic fields[END_REF][START_REF] Szmidt | Adaptive Damping of a Double-Beam Structure Based on Magnetorheological Elastomer[END_REF][START_REF] Dyniewicz | Semi-active control of a sandwich beam partially filled with magnetorheological elastomer[END_REF]. In these studies, all contact interfaces are assumed to be full sticking. In contrast to the prior studies, we aim to explore the possibility of passive vibration attenuation by friction in sandwich structures with permanently magnetized MAE.

BACKGROUND

To demonstrate the concept, a sandwich structure is assembled from two steel skins, five MAE patches and an aluminum spacer (figure 1a). Our MAE is a permanently magnetized composite material with a residual magnetic flux density of 0.3 T [START_REF] Hermann | Magnetic and dynamic mechanical properties of a highly coercive MRE based on NdFeB particles and a stiff matrix[END_REF]. The sandwich is tested under clamped/free boundary conditions. The sandwich is designed so that the clamp does not influence the vibratory responses of the structure. The outer dimensions of the sandwich measure 22 cm x 3 cm x 0.26 cm.

In the sandwich, the ferromagnetic skins are attracted by the magnetized MAE. The magneto-mechanical interactions maintain the skins clamped to the core without additional means of fixation. The normal stress in the contact interfaces (figure 1b) is hence established by the magneto-mechanical interactions. When the structure is subjected to transverse vibrations, it is deformed and shear stresses appear in the contact interface. While the magnetic normal stresses remain constant, the tangential shear stresses vary with the excitation magnitude and hence determine the contact state.

In the experiment, the sandwich is clamped directly on an electromagnetic shaker (figure 1b). A reference signal is measured by an accelerometer placed on the base. The velocity of the sandwich beam is measured with a contactless laser vibrometer. A harmonic excitation is performed with different acceleration magnitudes around the resonance frequency to study the nonlinear response of the structure. For the harmonic excitation, stepped-sine signals with amplitudes between 1 m/s² and 90 m/s² are used in tests with increasing frequency ("f up") and decreasing frequency ("f down"). Figure 1c shows the experimentally obtained FRF amplitudes of the third resonance mode for the different tests.

ANALYSIS

The structure shows a nonlinear behavior with respect to the excitation magnitude: the resonance frequency decreases with increasing amplitude and the shape of the amplitude peaks is modified. For the lowest excitation magnitudes, the FRF is symmetric around the resonance frequency. When the excitation magnitude increases, the maximum value of the FRF decreases and the peak is distorted to the left. The decrease of the peak amplitude of the FRF shows that the damping increases with the acceleration magnitude from 1 m/s² up to 50 m/s². For a further increase of the excitation magnitude, the peak of the transfer function straightens, its shape becomes more symmetric. The inclination of the resonance peak and the nonlinear evolution of the peak amplitude are characteristic for systems which include friction damping. For small excitation magnitudes, the tangential stresses are small. The contact surfaces mainly stick together and the system behaves like a single beam structure. With increasing excitation magnitude up to 50 m/s², more parts of the contact interface start to slide as the shear stresses increase. The stiffness of the structure is therefore decreased and friction damping occurs. Above the threshold of 50 m/s², the shear stresses become much larger than the normal stresses. The contact surfaces mainly slide above each other which reduces the effect of friction damping. Since the nonlinear damping effect becomes smaller, the resonance peak straightens again.

A finite element model of the experiment has been implemented to obtain information about the local contact state. The results of the numerical simulations confirm the constant, relatively distributed contact pressure. Moreover, in studies with increasing excitation magnitudes, the local contact state evolves from mainly sticking to mainly sliding which confirms the experimental observations.

CONCLUSION

The present work outlines the concept of magnetically induced friction damping in sandwich structures which contain permanently magnetized magnetoactive elastomers. The sandwich structure studied in this work illustrates how to employ magnetoactive elastomers for friction damping. Compared to rigid magnets, the MAEs are more lightweight and can adapt to the deformation of the assembly. The main drawback of the concept is the added mass. A possible subject for future studies is the influence of the sandwich design on the FRF. Instead of using a full core, isolated patches could be placed at different locations in the sandwich for example. Harmonic excitations with different controlled amplitudes reveal the characteristic, nonlinear behavior of the sandwich structure and proof the existence of friction damping. The advantage of a shape adaptability is that is maintains the contact interface which favors friction damping. The local contact state has been analyzed by the help of a numerical model of the experiment. The simulation results indicate that the contact state evolves with the vibration amplitude which confirms the experimental results. The characteristic evolution of the FRF amplitudes clearly shows the presence of nonlinear damping. It is assumed presently that is results from friction. However, additional passive mechanisms of energy dissipation like viscoelastic damping and eddy current damping are probably activated in the sandwich structure. The contribution of these different passive damping mechanisms will be quantified in future studies.

Figure 1 :

 1 Figure 1: Components of the MAE sandwich structure (a), schematic side view of the experimental setup (b) and FRF amplitudes for the third bending mode of the structure with different excitation amplitudes in frequency-up ("f up") and frequency-down ("f down") tests (c).

ACKNOWLEDGEMENTS

This work has been funded by the French company Moving Magnet Technologies SA in collaboration with EIPHI Graduate School, ANR-17-EURE-0002.


R Boudot

M Abdel Hafiz

M Petersen

E Rubiola

A CPT-based Cs cell self-sustained microwave oscillator

Keywords: microwave oscillator, phase noise, AM noise, vapor cell

We report the implementation of a CPT-based Cs cell microwave feedback oscillator. Here, a high-Q atomic resonance initiated through CPT in a buffer-gas filled Cs vapor cell is detected with a fast photodiode, amplified and sustained in a loop oscillator configuration, comparable to those of optoelectronic oscillators (OEOs). The Cs cell oscillator currently demonstrates at 9.192 GHz a phase noise of -24, -82 and -112 dBrad 2 /Hz, for offset frequencies of 1, 100 and 10 kHz, respectively, in good agreement with the Leeson model. The AM noise was also measured, showing a comparable floor and 1/f dependence for offset frequencies lower than 200 Hz

I. INTRODUCTION

In a CPT-based vapor cell clock, the clock resonance signal is usually detected with a slow photodiode at the cell output as an increase of the light transmitted through the atomic medium. Alternatively, other studies in the early 2000s have demonstrated the possibility to directly detect opticallyinduced microwaves generated at the cell output with a fast photodiode and subsequently construct a self-oscillating atomic frequency reference [1][2][3][4]. Here, we aim to stimulate a renewed interest to the unusual latter concept with the demonstration of an all-optical microwave feedback oscillator using a Cs vapor cell resonator.

II. METHODS/RESULTS

Compared to previous studies [1][2][3][4], our experimental setup includes the use of a fibered Mach-Zehnder electro-optic modulator driven at 4.596 GHz that forces the request of a frequency divider in the loop. In addition, we report a simple and easy-to-understand phase noise characterization of the CPT-based feedback oscillator. The phase noise of the Cs cell atomic oscillator, in agreement with the Leeson model and shown in Fig. 1, is -24 dBrad 2 /Hz at 1 Hz, with a flicker frequency noise compatible with a stability of 8×10 -12 at 1 s, and currently degraded by the flicker phase noise of the EOM stage converted into frequency fluctuations of the oscillator loop. The noise floor of the oscillator is -112 dBrad /Hz and is limited by the low microwave power at the photodiode output.

We have also measured the AM noise of the oscillator, shown on Fig. 1, obtaining a spectrum well explained by the formalism developed in [START_REF] Rubiola | A generalization of the Leeson effect[END_REF]. These results were submitted for publication [START_REF] Boudot | [END_REF]. Our analysis suggests that the reduction of the residual noise of the sustaining loop might pave the way to stability levels at 1 s in the 4×10 -13 to 5×10 -13 range, competitive with those of state-of-the-art vapor cell clocks [5][6]. In that sense, the exploration and test of noise reduction techniques to reduce phase fluctuations of the sustaining loop is under progress. Latest results will be presented at the conference

FIG. 1 .

 1 FIG. 1. PM AND AM NOISE, MEASURED AT 9.192 GHZ, OF THE CS CELL FEEDBACK OSCILLATOR