2 Méthodes et outils

Chapitre 2
Méthodes et outils

Au cours de ce chapitre, je présenterai les deux outils fondamentaux utilisés durant ma thèse. Un code numérique d’intégration des équations du mouvement, en application directe du formalisme hamiltonien présenté précédemment, et l’Analyse en Fréquence.

Dans une première partie, je discuterai d’abord des motivations qui m’ont conduit à me lancer dans l’écriture d’un code numérique. Ensuite, je donnerai une brève description des diﬀérentes philosophies et méthodes employées et expliquerai leurs domaines de validité. Je m’attarderai plus particulièrement sur l’avantage des nouvelles méthodes (schéma symplectique, algèbre de Lie) utilisées dans le monde des accélérateurs. Bien que ces méthodes soient aujourd’hui très répandues aux Etats-Unis, elles ont plus de mal à s’établir en Europe¹. Je présenterai alors une nouvelle classe d’intégrateurs symplectiques à pas tous positifs développée par Jacques Laskar (2000). Puis, je l’appliquerai aux Hamiltoniens décrivant chacun des principaux éléments magnétiques d’un accélérateur. Enﬁn, je donnerai une comparaison avec l’intégrateur de Ruth (1983), intégrateur symplectique le plus utilisé pour les machines à électrons.

Après avoir intégré les équations du mouvement, la dynamique des diﬀérentes machines à rayonnement synchrotron sera étudiée au moyen de l’Analyse en Fréquence. L’Analyse en Fréquence est l’outil central de mon travail de thèse. Je commencerai, dans la seconde partie de ce chapitre, par décrire la méthode sans noyer le lecteur dans des détails trop techniques. C’est un outil extrêmement puissant qui a été développé par Jacques Laskar initialement pour étudier la dynamique du Système Solaire. Avant de l’appliquer à la dynamique des accélérateurs (voir le chapitre 3), je présenterai quelques résultats pédagogiques sur le pendule rigide et l’application d’Hénon. Ces deux applications à un et deux degrés de liberté auront la particularité de permettre au lecteur de pouvoir lire plus simplement les cartes en fréquence des accélérateurs ; de plus, ils me permettront de présenter brièvement quelques résultats généraux sur les résonances et la stabilité des systèmes dynamiques.

2.1 Réalisation d’un intégrateur symplectique

2.1.1 Introduction

L’écriture d’un code d’intégration des équations du mouvement est arrivée tardivement dans mon travail de thèse. Dans le langage accélérateur, on parle souvent de code de tracking. De multiples raisons expliquent ce choix :

Pour commencer, de nombreux codes de tracking avaient déjà été écrits. Par exemple le programme MAD développé continuellement au CERN [50] est largement employé par la communauté des physiciens des accélérateurs. Cependant, généralement à chaque accélérateur correspond un ou plusieurs codes spéciﬁques. Ces codes sont modiﬁés, améliorés constamment et localement pour permettre de prendre en compte de nouveaux phénomènes physiques, de nouveaux éléments magnétiques. Très rapidement, le résultat peut être un code volumineux qui, souvent, est adapté à des tâches diverses et variées, allant de la conception d’un accélérateur à l’optimisation de sa dynamique et de ses performances. Le groupe théorie faisceau d’une machine peut être divisé schématiquement en deux sous-groupes : une personne spécialiste et développeur des codes utilisés par les membres du second groupe.

J’ai découvert le monde des accélérateurs au sein de l’équipe du Projet SOLEIL, il y a bientôt trois ans. Lorsque j’ai commencé à étudier la dynamique de SOLEIL, on m’a donné une description de la maille théorique de l’anneau de stockage et le manuel d’utilisation du logiciel BETA [93] pour faire le tracking. Le même processus s’est répété pour l’étude de l’Advanced Light Source avec le code DESPOT [42] (ici il n’y a pas de manuel d’utilisation).

Donc dans un premier temps, les codes d’intégrations se sont résumés à des boîtes noires « censées » donner les bons résultats. De plus, l’écriture d’un code de tracking est considérée comme une tache ardue par les utilisateurs. Mon directeur de thèse, astronome, prétendait le contraire après avoir travaillé avec les physiciens de l’ALS (1993). Initialement, je ne l’ai pas crû. Puis, en travaillant sur diﬀérentes sources de lumières (cf. chapitre 3), j’ai dû presque utiliser un code de tracking diﬀérent par machine. J’ai été confronté à des problèmes de modélisation, à des incompatibilités entre codes. Plus préoccupant, les codes peuvent donner des résultats complètement diﬀérents aux grandes amplitudes du mouvement et parfois même pour les paramètres linéaires. Le temps de calcul peut également être du simple au double suivant les méthodes utilisées.

Enﬁn, à ces motivations doivent encore s’en ajouter deux autres. Premièrement à ma connaissance, il n’existe pas dans la communauté française² de code de tracking optimisé et adapté à l’étude de la dynamique à long terme des accélérateurs d’électrons. Deuxièmement, l’équipe Astronomie et Systèmes Dynamiques est un lieu regroupant toutes les compétences nécessaires pour l’écriture d’un tel intégrateur. Ces dernières années, plusieurs intégrateurs y ont d’ailleurs été développés pour la Mécanique Céleste.

Le but principal de ce travail est donc d’obtenir un code de calcul eﬃcace permettant d’intégrer le mouvement d’une particule dans un élément magnétique puis dans l’accélérateur complet. L’intégrateur doit être à la fois rapide pour réaliser des études à long terme (plusieurs millions de tours), précis pour étudier ﬁnement la dynamique et valide à grande amplitude pour prendre correctement en compte l’eﬀet des éléments multipolaires. Dans la première version du code, la modélisation des défauts (déplacements, rotations des éléments magnétiques) ne nous intéresse pas ; cependant leur implémentation ne poserait aucune diﬃculté.

Je vais commencer par décrire les deux grandes approches selon lesquelles les codes numériques de tracking sont écrits. Mais tout d’abord, déﬁnissons les notions d’information locale et globale dans un accélérateur (cf. Forest et Hirata, 1992 et Forest, 1998).

Une quantité est appelée locale si elle dérive d’un élément magnétique indépendamment de sa position dans l’anneau et de l’anneau lui-même : intégration d’une particule à travers l’élément, le système de coordonnées de référence, la déviation de l’orbite nominale, le champ magnétique ou l’enveloppe du faisceau.

Une quantité est dite globale si elle ne peut être déduite qu’à partir de la connaissance de l’anneau complet : l’orbite fermée, les nombres d’ondes, les fonctions de Twiss ou l’émittance d’équilibre.

2.1.2 Ancienne approche

Dans l’ancienne approche, on déﬁnit un seul et unique système de référence (coordonnées curvilignes), le champ électromagnétique est alors développé dans ce système pour obtenir l’expression, soit d’un Hamiltonien global, soit d’équations newtoniennes globales. Les équations du mouvement sont exprimées de manière à permettre simultanément le tracking et le calcul des quantités globales caractérisant l’anneau. Bien que possible, cette approche devient très compliquée avec l’augmentation de la complexité et de la taille d’un anneau de stockage.

Les solutions ( $x$ ) du mouvement sont obtenues en eﬀectuant un développement de Taylor des équations du mouvement autour d’une orbite dite de référence — typiquement l’orbite fermée qui doit préalablement avoir été calculée — ; cependant son obtention est pénible à partir de l’ordre 3. L’expression générale entre le vecteur position initial ( $i x$ ) et ﬁnal ( $f x$ ) s’écrit :

∑ ∑ ∑ xfj = 𝒦j + ℛjkxik + 𝒯 jkmxikxim + 𝒰jkmnxikximxin + 𝒪(|x i|4) k km kmn

(2.1)

où $𝒦j$ représente l’orbite fermée, $ℛjk$ l’optique linéaire, $𝒯 jkm, 𝒰jkmn$ , les aberrations optiques d’ordre un et deux (voir Brown, 1982). L’optique linéaire est obtenue en tronquant la série de Taylor au premier ordre. Les éléments magnétiques sont modélisés par des matrices, puis selon les besoins, on inclut les termes d’ordres supérieurs. Généralement, le développement est arrêté au second ordre pour permettre d’inclure les éléments multipolaires qui sont presque toujours exprimés sous forme de lentilles minces, en particulier les inévitables hexapôles (voir par exemple les logiciels BETA [93], TRANSPORT [50]).

Cette approche est historiquement la première a avoir été développée (Courant et Snyder, 1958 et Bruck, 1966). Initialement, elle été vouée à la conception des accélérateurs, au tracking sur quelques tours, aﬁn de calculer les grandeurs linéaires et d’analyser la stabilité du mouvement au premier ordre. Une des hypothèses fondamentales utilisées est l’approximation des petits angles ou approximation dite des grandes machines, hypothèse qui rend les dipôles et quadripôles linéaires (et les sections droites : « éléments sans champ »).

Cette méthode a l’avantage de permettre d’obtenir simplement et rapidement les informations globales de l’anneau (e.g. fonctions de dispersion, orbite fermée, nombres d’ondes, chromaticités, fonctions de Twiss). Autre avantage, le temps de calcul de l’ouverture dynamique (paramètre essentiel pour l’optimisation) est suﬃsamment faible (quelques minutes) pour permettre de tester un grand nombre d’optiques de la machine³.

Je voudrais insister sur le fait qu’il n’existe pas de logiciel de calcul adapté à tous les besoins rencontrés dans le monde des accélérateurs : chaque code de calcul est écrit dans un cadre très précis (développement au voisinage de l’orbite fermée, approximation des grandes machines, conception de machines, optimisation d’optiques), cadre qui trop souvent est oublié par leurs utilisateurs. De plus, il faut ajouter que les performances exigées aujourd’hui pour les accélérateurs ne sont plus du tout les mêmes qu’il y a vingt ou trente ans. A l’époque, les machines étaient très linéaires, avaient de faibles acceptances en énergie et dynamique ; aujourd’hui, pour une machine comme SOLEIL, on exige des performances extrêmes — acceptance en énergie de six pour cent, grande ouverture dynamique, durée de vie de plus de dix heures, haute brillance, faible émittance, nombreux dispositifs d’insertion — et tout cela avec des composantes multipolaires de plus en plus fortes, rendant l’optimisation de plus en plus délicate du fait des nombreuses résonances induites.

Après cette parenthèse, nous convenons qu’il serait dangereux d’utiliser un code hors de son cadre. En particulier, si l’utilisation d’un développement de Taylor est adapté pour les LINAC (que la particule ne traverse qu’une seule fois), il serait illusoire de vouloir utiliser aujourd’hui cette approche pour une étude ﬁne de la dynamique et pour des études à long terme, et qui plus est, à grande amplitude (cf. infra Super-ACO, une machine à faible rayon de courbure). De par la non symplecticité de la troncature de la série de Taylor (en général), lorsque la trajectoire d’une particule est intégrée sur plusieurs milliers de tours pour des électrons ou plusieurs dizaines de millions de tours pour des protons, on observe des amortissements artiﬁciels conduisant à des orbites spirales et à la présence de résonances factices.

2.1.3 Approche moderne

La méthode dite moderne consiste à découpler la simulation de l’analyse, i.e. séparer la construction de l’application de premier retour et de l’intégrateur (problèmes locaux) de l’analyse des grandeurs caractéristiques déduites de l’application de premier retour (problèmes globaux) ; cette séparation doit avoir lieu autant du point de vue de la conception que du point de vue de la programmation.

2.1.3.1 Présentation

La construction d’un code de tracking est un problème entièrement local : chaque élément est déﬁni par un Hamiltonien local dont le système de référence est choisi à partir de considérations purement géométriques (formes et symétries du champ magnétique). Suivant le type d’élément magnétique telles ou telles considérations seront privilégiées, considérations qui peuvent être contradictoires d’un élément à l’autre. Il faut donc trouver une méthode pour relier les éléments les uns avec les autres : c’est l’approche LEGO si l’on reprend la terminologie de Forest et Hirata (1992). Comme nous l’avons déjà dit, le choix du système de coordonnées est entièrement libre ; cependant deux grandes familles de blocs sont à distinguer : le bloc courbe et le bloc rectiligne.

Le bloc rectiligne est un bloc avec deux faces parallèles séparées d’une distance $L$ , sur chacune desquelles est attaché un système local de coordonnées cartésiennes (Fig. 2.1-a). Les repères sont les mêmes pour les deux faces et sont reliés par une ligne droite. Ce type de bloc peut être utilisé pour la description des quadripôles, hexapôles, 2n-pôles, wiggler.
Le bloc courbe est un bloc pour lequel les axes verticaux $(Oy )$ des repères entrant et sortant sont parallèles alors que les axes horizontaux $(Ox )$ se coupent suivant un angle $θ$ (Fig. 2.1-b). Les origines des deux repères sont reliées par un arc de cercle de rayon $ρ$ . Un tel bloc a pour premier but de courber la trajectoire d’une particule d’un angle d’environ $θ$ degrés. Ce bloc sera plus particulièrement adapté à la description d’aimants de courbure (aimant dipolaire secteur ou droit).

pict

FIG. 2.1:

Schéma succinct d’un bloc rectiligne (a) et courbe (b) de longueur $L$ . Les axes de coordonnées sont portés par les faces d’entrée et de sortie des blocs — inspiré de Forest et Hirata (1992) —

2.1.3.2 Le ﬂot d’un système

Le but des blocs de base est de permettre de construire et calculer l’application qui transporte une particule de la face d’entrée (i) à la face de sortie (f) ; cette application, $ℳ$ , est généralement appelée application de transfert ou ﬂot de l’élément :

xf = ℳi →f xi

(2.2)

Remarquons qu’un élément magnétique pourra être composé de plusieurs blocs élémentaires. L’élément physique sera donc déﬁni par le couple bloc (rectiligne, courbe) et un modèle qui donne l’application de transfert à travers l’élément.

Donnons dès à présent, deux propriétés essentielles du ﬂot. Soit deux éléments notés (1) et (2), et $ℳ 0→1$ est le ﬂot associé à (1) et $ℳ 1→2$ celui associé à (2) alors le ﬂot résultant pour le système (1–2) est simplement l’application composée :

ℳ0 →2 = ℳ1 →2 ℳ0 →1

(2.3)

Pour un système autonome qui possède une périodicité (l’anneau), on est amené à calculer le ﬂot de l’anneau complet $𝒜$ sur $n$ tours de circonférence $2πR$ , alors :

n 𝒜2nπR →2(n−1)πR ...𝒜2 πR→4πR 𝒜0→2 πR = (𝒜0 →2πR)

(2.4)

Cette approche est bien locale, chaque élément est modélisé par un Hamiltonien exprimé dans les coordonnées naturelles (pour son intégration).

2.1.3.3 Intégrateur symplectique ab ovo

Nous allons introduire quelques notions d’algèbre de Lie et déﬁnir l’intégration symplectique (voir Lectures on nonlinear orbit dynamics de A. Dragt (1982) pour une introduction à l’algèbre de Lie appliquée au monde des accélérateurs).

Soit $ℋ$ le Hamiltonien d’un système autonome à $n$ degrés de liberté et $x = (pj,qj)j=1..n$ le vecteur de coordonnées et moments généralisés et $s$ , la coordonnée indépendante. Alors les équations de Hamilton sont :

dp ∂ ℋ dq ∂ℋ --j-= − ----, --j = ---- pour j = 1..n ds ∂qj ds ∂pj

(2.5)

Si on déﬁnit les crochets de Poisson de deux fonctions f et g des variables $x$ par :

d´ef∑ ∂f ∂g ∂f ∂g {f,g} = --------− -------- j ∂pj ∂qj ∂qj ∂pj

(2.6)

alors les équations 2.5 peuvent se réécrirent de manière plus concise :

dx --- = {ℋ, x} = Lℋx ds

(2.7)

avec $Lℋ$ l’opérateur diﬀérentiel déﬁni par $d´ef L ℋf = {ℋ, f }$ . On peut alors intégrer formellement les équations du mouvement pour des conditions initiales données $i x$ :

∑ sn xf = --Lnℋx i d=´efesLℋx i n≥0 n!

(2.8)

En pratique, la diﬃculté est l’évaluation du terme $esLℋxi$ de l’équation 2.8, car bien qu’étant une série convergente, la convergence est souvent lente et nécessite l’évaluation d’un grand nombre de termes dépendant du pas de temps $s$ et de la forme de $ℋ$ . Le principe de base utilisé pour intégrer l’équation 2.7 consiste à approximer l’opérateur $sLℋ e$ de manière symplectique.

Formellement, la solution exacte du système 2.8 est symplectique : le ﬂot hamiltonien, transportant le vecteur de l’espace des phases $x$ du temps $0$ au temps $s$ le long d’une trajectoire du système d’énergie $ℋ$ , conserve exactement la 2-forme symplectique⁴ :

f f i i dp ∧ dq = dp ∧ dq

(2.9)

Une méthode d’intégration vériﬁant cette propriété est appelée intégrateur symplectique. Un développement de Taylor à l’ordre $n$ ne vériﬁe généralement ni la condition de symplecticité 2.9 ni la conservation de l’énergie. Remarquons qu’idéalement, on souhaiterait également que l’énergie soit exactement conservée, cependant Ge et Marsden (1988) ont montré qu’il n’existe pas en général de tel schéma d’intégration pour un Hamiltonien non intégrable. Cependant, l’intégrateur symplectique conserve exactement le Hamiltonien $ℋ˜$ d’un système perturbé voisin de celui décrit par le Hamiltonien $ℋ$ initial :

˜ ∑ k ℋ = ℋ + s Hk k≥1

(2.10)

où $Hk$ est fonction des dérivés d’ordre $j ≤ k$ de $ℋ$ (Yoshida, 1990a, 1990b). Par déﬁnition, un intégrateur symplectique est d’ordre $n$ si :

˜ n ℋ = ℋ + 𝒪 (s )

(2.11)

2.1.3.4 Méthodes de construction

Brièvement, nous pouvons diviser les méthodes d’intégration en deux grands groupes (Yoshida, 1993) :

les méthodes implicites reposant sur
– un schéma de Runge-Kutta implicite (Sanz-Serna, 1988) ;
– l’utilisation de fonctions génératrices (Channel et Scovel, 1990).
les méthodes explicites
– Ruth, le premier, a introduit l’utilisation d’intégrateur symplectique (1983), a initialement utilisé des transformations canoniques successives pour réaliser des intégrateurs jusqu’à l’ordre 4 ;
– Neri (1988) a reformulé le problème en utilisant l’algèbre de Lie ;
– Yoshida (1990) a donné un schéma général pour construire un intégrateur d’ordre $2n+ 2$ à partir d’un intégrateur d’ordre $2n$ ;
– Abell (1995) a introduit l’utilisation d’applications de Crémona (applications polynomiales symplectiques).

Nous nous sommes intéressés plus particulièrement à un schéma explicite que nous détaillerons un peu plus loin. Tout d’abord, déﬁnissons un petit peu la philosophie de construction d’un intégrateur. Elle repose sur cinq points essentiels :

Le Hamiltonien est local ;
Il est strictement impossible de connaître le Hamiltonien réel de l’élément ou de l’accélérateur de par les erreurs de construction, les erreurs aléatoires (champs, déplacements, marées, vibrations) ;
Nous nous intéressons à la modélisation de la dynamique la plus proche de la dynamique réelle de l’accélérateur (nonlinéarités, stabilité) ;
Si la solution exacte d’un Hamiltonien est connue, alors elle est symplectique, sinon le sus-dit Hamiltonien est découpé en parties séparément intégrables ;
Si besoin, l’ordre de l’intégrateur peut être augmenté à volonté (aux dépens du temps de calcul).

A la suite de Forest (1998, passim), nous considérons que l’application générée par un intégrateur doit toujours être vue comme un nouveau modèle de l’aimant. C’est pourquoi les paramètres linéaires doivent — si besoin — être réajustés aﬁn d’obtenir avec l’intégrateur symplectique une machine dont les nombres d’ondes et paramètres linéaires correspondent à un intégrateur « matrice de transfert-kick ». Nous ne rentrerons pas dans les discussions des partisans de l’ancienne école, qui prétendent que leur méthode d’intégration est plus exacte que celle fondée sur un schéma symplectique. Dans la première, les éléments sont dits épais car représentés matriciellement (de longueur non négligeable), dans la dernière, les éléments sont parfois interprétés comme « coupés en rondelles » suivant un schéma « section droite-lentille mince–section droite » (cf. Peggs et Talman (1986), Schachinger et Talman (1987)).

Nous allons maintenant présenter les deux principales approches rencontrées dans la plupart des codes modernes dans le monde des accélérateurs.

Décomposition en produit d’applications symplectiques : La première approche consiste à écrire le ﬂot hamiltonien sous la forme d’un produit d’applications symplectiques (Dragt et Finn, 1976, Dragt et al., 1988, Dragt, 1996), i.e. sous la forme :

f i sL i sL sL sL sL i x = ℳx = e ℋx = e f1e f2e f3e f4 ...x

(2.12)

où les fonctions $fk$ sont des polynômes homogènes de degré $k$ déterminées à partir du Hamiltonien $ℋ$ du système considéré. Le facteur $sLf2 e$ représente la partie linéaire de la transformation usuellement écrite par une matrice de transfert $ℛ$ (cf. Eq. 2.1). L’application $esLf3$ décrit les termes nonlinéaires d’ordre deux (et parfois supérieurs) correspondant aux coeﬃcients $𝒯 jkm$ , $𝒰 jkmn$ , etc. du développement de Taylor 2.1 (Dragt et Forest, 1986 – Dragt et Abell, 1996 – Dragt, 1996). Le produit d’applications symplectiques 2.12 peut être tronqué sans aucune diﬃculté à n’importe quel ordre $n$ :

f i sLf1 sLf sLf i x = ℳx = e e 2 ...e n+1x

(2.13)

Par construction la troncature est symplectique. Formellement, cette procédure est appliquée pour chaque élément. Il est aussi possible de calculer l’application de transfert $𝒜$ pour l’anneau complet, on obtient alors l’application de premier retour qu’il suﬃt d’itérer pour obtenir un logiciel de tracking symplectique d’ordre $n$ .

Une diﬃculté pratique apparaît, car l’application de transfert $ℳ$ n’est pas intégrable pour $k$ supérieur à 2. Cette diﬃculté est levée par l’utilisation d’une fonction génératrice $F(qi,pf)$ :

i f i f i f F(q ,p ) = q ⋅ p + P (q ,p )

(2.14)

où $P$ est un polynôme à déterminer. On obtient alors la transformation canonique des variables :

{ i d´ef ∂F- f -∂P- i f p = ∂qi = p + ∂qi(q ,p ) qf d=´ef ∂Ff = qi + -∂Pf(qi,pf ) ∂p ∂p

(2.15)

Le polynôme $P$ est déterminé de manière à ce que les applications 2.13 et 2.15 soient égales jusqu’à l’ordre $n$ : il s’exprime en fonction des générateurs $fk$ (pour plus de détail le lecteur pourra se rapporter par exemple au guide des méthodes physiques de logiciel MAD [58] ou au manuel de MARYLIE [29]).

On remarque que la première équation de l’expression 2.15 est implicite en $pf$ . Elle est résolue numériquement en utilisant une méthode de Newton. La deuxième équation peut ensuite être calculée.

Cette approche est suivie, par exemple, par les logiciels MARYLIE [29] et MAD [58]. Une approche voisine consiste à approximer $ℳ$ par des fonctions de Crémona (voir par exemple Dragt et Abell, 1996).

Intégration directe : La seconde approche, que nous avons retenue, repose sur l’intégration directe du Hamiltonien d’un élément de l’anneau. Ce choix se justiﬁe par notre volonté de réaliser un code d’intégration optimal en termes de vitesse et de précision. Cette méthode est simple à mettre en œuvre et découple deux problèmes distincts, à savoir la simulation et l’analyse globale de l’anneau. Si le Hamiltonien est intégrable et si la solution analytique se prête bien au calcul numérique, nous la conserverons. Dans le cas contraire, nous décomposons le Hamiltonien en deux parties séparément intégrables, $ℋ = A + εB$ où $ε$ peut être vu comme un paramètre de perturbation ; puis nous utilisons le théorème de Baker-Campbell-Hausdorﬀ (BCH) qui stipule que l’on peut écrire formellement :

esLAesLεB = esL˜ℋ

(2.16)

avec le Hamiltonien formel :

s s2 ℋ˜ = A + εB + -{A, εB } + ---({A {A,εB }} + {εB {εB, A }}) + ... 2 12

(2.17)

Ainsi l’intégrateur symplectique du premier ordre le plus simple est :

sL sL sL e ℋ = e Ae εB + 𝒪 (s ε)

(2.18)

Avec l’aide du théorème BCH, des intégrateurs ( $Sn$ ) d’ordre quelconque $n$ peuvent être construits :

∏n Sn (s) = ecisLAedisLεB i=1

(2.19)

où les coeﬃcients $(ci,di)i=1..n$ sont déterminées pour que le reste soit d’ordre $n$ . En particulier, les deux relations suivantes doivent être vériﬁées :

∑n ∑n ci = 1 et di = 1 i=1 i=1

(2.20)

Dans la suite, on se restreint à la classe des intégrateurs symétriques, i.e. $S−n1(s) = Sn (− s)$ . Par construction ces intégrateurs sont de reste paire $𝒪 (s2nε)$ . Le célèbre intégrateur d’ordre 2 appelé dans la littérature leapfrog integrator (« saute-mouton ») s’écrit (Ruth, 1983) :

S2 = ec1sLAed1sLεBec1sLA

(2.21)

avec $c=1- 1 2$ et $d = 1 1$ .

Une manière simple pour obtenir des intégrateurs d’ordre plus élevé repose sur le travail de Yoshida (1990). L’idée principale est de composer le précédent intégrateur d’ordre 2 pour obtenir un intégrateur symétrique d’ordre plus élevé. Donc, on pourra écrire un intégrateur d’ordre 4 comme :

S4(s) = S2(as)S2(bs)S2(as)

(2.22)

avec pour condition de faire un pas total de longueur $s$ et d’annuler le terme d’ordre 2, soit :

{ ( --1-- 2a + b = 1 { a = 2−213 2a3 + b3 = 0 ⇐⇒ ( -213-- b = − 2−213

(2.23)

De manière générale, Yoshida montre qu’un intégrateur d’ordre $2n + 2$ peut être construit à partir d’un intégrateur d’ordre $2n$ en utilisant le schéma symplectique :

S2n+2(s) = S2n(as)S2n(bs)S2n(as)

(2.24)

Immédiatement, on obtient : $a = ---1-√-- 2− 2n+12$ et $- b = − -2n+1√2√-- 2− 2n+12$ .

Les deux intégrateurs les plus fréquemment utilisés en physique des accélérateurs circulaires⁵ sont :

– l’intégrateur d’ordre 4 ( $S 4$ ) — initialement établi par Ruth (1983, puis reformulé par Forest et Ruth, 1990) — que l’on peut aussi écrire :

d1sLA c2sLεB d2sLA c3sLεB d2sLA c2sLεB d1sLA S4(s) = e e e e e e e

(2.25)

avec $c2= 11+α-$ , $c3 = (α − 1)c2$ , $d1 = c22$ , $d2 = αd1$ et $1 α = 1 − 23$ .

– l’intégrateur d’ordre 6 ( $S6$ ) (Yoshida, 1990) :

S6(s) = S2(ds)S2(cs)S2(bs)S2(as)S2(bs)S2(cs)S2 (ds )

(2.26)

avec trois jeux de solutions possibles (cf. Tab. 2.1) et avec toujours un grand pas négatif.


	Solution 1	Solution 2	Solution 3

a	-0.117767998417887E1	-0.213228522200144E1	+0.152886228424922E2
b	+0.235573213359357E0	+0.426068187079180E2	-0.214403531630539E1
c	+0.784513610477560E0	+0.143984816797678E1	+0.144778256239930E1

d	1-2(a+b+c)

TAB. 2.1:

Trois jeux de coeﬃcients pour construire un intégrateur symplectique d’ordre 6 en utilisant la méthode de Yoshida. Dans chacun des cas, il existe un grand pas négatif provenant de la relation 2.20 : $d + 2 (a + b + c) = 1$ .

2.1.4 Intégrateurs utilisés

La classe d’intégrateurs que nous avons utilisée est présentée dans l’article de Laskar et Robutel (2000) auquel le lecteur pourra se référer pour plus de détails.

Si l’on distingue $A$ et $B$ dans l’expression générale du Hamiltonien, les intégrateurs symétriques peuvent appartenir à l’une des deux classes $𝒮 𝒜 ℬ𝒜k$ et $𝒮ℬ 𝒜 ℬk$ déﬁnies par :

_2n	: e^c₁sL_Ae^d₁sL_εB…e^d_nsL_εBe^c_n+1sL_Ae^d_nsL_εB…e^d₁sL_εBe^c₁sL_A
_2n+1	: e^c₁sL_Ae^d₁sL_εB…e^c_n+1sL_Ae^d_n+1sL_εBe^c_n+1sL_A…e^d₁sL_εBe^c₁sL_A
_2n	: e^d₁sL_εBe^c₂sL_Ae^d₂sL_εB…e^d_nsL_εBe^c_n+1sL_Ae^d_nsL_εB…e^d₂sL_εBe^c₂sL_Ae^d₁sL_εB	(2.27)
_2n+1	: e^d₁sL_εBe^c₂sL_Ae^d₂sL_εB…e^c_n+1sL_Ae^d_n+1sL_εBe^c_n+1sL_A…e^d₂sL_εBe^c₂sL_Ae^d₁sL_εB

Par exemple l’intégrateur leapfrog (Eq. 2.21) appartient à la classe $𝒮 𝒜 ℬ𝒜1$ . C’est un intégrateur d’ordre 2 avec pour Hamiltonien formel $ℋ˜ = A + εB + 𝒪 (s2ε)$ . L’intégrateur d’ordre 4 de Forest et Ruth (Eq. 2.25) appartient à la classe $𝒮 𝒜 ℬ𝒜3$ avec $ℋ˜ = A + εB + 𝒪 (s4ε)$ . Nous pouvons remarquer la présence de deux pas négatifs pour cet intégrateur, ce qui rend la valeur absolue de chaque pas intermédiaire relativement grande pour un pas total de 1 : $d1 ≈ 0.6756$ , $d2≈−0.1756$ , $c2 ≈ 1.3512$ et $c3 ≈ − 1.7024$ . Il s’ensuit que pour de grands pas d’intégration, la méthode perd de son eﬃcacité (coût élevé, instabilités numériques). Suzuki (1991) a démontré qu’il est impossible de réaliser un intégrateur symplectique d’ordre $n ≥ 3$ avec des pas tous positifs. Cependant, le problème des pas négatifs peut être partiellement solutionné.

Jusqu’à présent, nous n’avons pas pris en compte l’existence du petit paramètre $ε$ . La méthode retenue consiste à déterminer les coeﬃcients $(cj,dj)$ des intégrateurs 2.27 pour avoir un reste d’ordre $𝒪 (snε + s2ε2)$ et non plus d’ordre $𝒪 (sn ε)$ .

Par exemple pour un intégrateur de classe $𝒮𝒜 ℬ 𝒜2$ (ordre quatre), on a :

SABA = ec1sLAed1sL εBec2sLAed1sLεB ec1sLA 2

(2.28)

l’unique solution pour les coeﬃcients est :

d1 = 1, c1 = 1(1 − c2) et c2 = √1-- 2 2 3

avec

( ) 2 2 1 c1 4 ℋ˜ = A + εB + s ε − ---+ -- {{A, B}, B} + 𝒪 (s ε) ◟-------24----4---◝◜------------------◞ 𝒪(s4ε+s2ε2)

(2.29)

De manière similaire pour un intégrateur de classe $𝒮ℬ 𝒜 ℬ2$ ,

d1sL εB c2sLA d2sLεB c2sLA d1sLεB SBAB2 = e e e e e

(2.30)

les coeﬃcients positifs solutions sont l’unique triplet :

1- 2- 1- d1 = 6 , d2 = 3 et c2 = 2

Ce sont ces deux intégrateurs qui ont été retenus pour l’écriture du code de tracking.

En fait, dans le cas particulier où $A$ est quadratique en les impulsions et $B$ ne dépend que des positions, il est possible d’améliorer encore la méthode en introduisant un correcteur $𝒞$ déﬁni par (voir Laskar et Robutel, 2000) :

𝒞 = e−s3ε2c2L{{A,B},B}

(2.31)

où le coeﬃcient du correcteur $c$ est déterminé pour annuler le terme d’ordre $2 2 𝒪 (s ε )$ (voir par exemple l’expression 2.29 pour l’intégrateur $SABA2$ ). On notera que le correcteur introduit un pas négatif, mais qu’il est d’autant plus petit que l’ordre de la méthode est élevé (Tab 2.2).


n	$SABAn$	$SBABn$

1	$1∕12$	$− 1∕24$
2	$√ -- (2 − 3)∕24$	$1∕72$

TAB. 2.2:

Coeﬃcient $c$ du correcteur pour les intégrateurs $SABAn$ et $SABAn$ (extrait de Laskar et Robutel, 2000)

Typiquement, le schéma symplectique avec correcteur s’écrit, par exemple pour un intégrateur $SABA2$ :

SABAC2 = e−s3ε2c2L{{A,B},B}SABA2e −s3ε2c2L{{A,B },B}

(2.32)

L’intégrateur avec correcteur est encore symétrique et son reste est d’ordre $n 4 2 𝒪(s ε + s ε )$ .

Les intégrateurs que j’utiliserai sont d’ordre deux et quatre. Typiquement, un intégrateur d’ordre 4 sera utilisé pour intégrer les éléments de type dipôles, quadripôles alors qu’un intégrateur d’ordre $2$ suﬃra pour l’intégration des multipôles individuels.

Nous montrerons qu’un intégrateur d’ordre deux avec correcteur est plus précis d’un ordre de grandeur que l’intégrateur de Ruth (cf. infra). Je vais maintenant présenter le Hamiltonien local pour chacun des principaux éléments magnétiques et les approximations réalisées.

2.1.5 Intégration des éléments parfaits

Nous avons vu page §, l’expression la plus générale du Hamiltonien pour un accélérateur (cf. Eq. 1.35) exprimée en fonction des trois paires de variables canoniques $(x,px)$ , $(y,py)$ et $(l,δ)$ . Cette expression dépend explicitement de la longitude $s$ prise comme variable indépendante.

Dans toute la suite, sauf mention explicite, j’ai choisi de travailler sur la forme quadratique en les impulsions en réalisant un développement limité au premier ordre de la racine carrée. Pour mémoire, je rappelle son expression :

p2+ p2 ℋ(x, y,l,px,py,δ) = (1 + hx )-x----y − hx(1 + δ) − eAˆs 2(1 + δ)

(2.33)

La première approximation réalisée a consisté à développer la racine carrée de l’expression 1.35, terme nonlinéaire qui caractérise le fait que les particules sont relativistes. Ce développement permet d’obtenir une expression polynomiale dans les variables accélérateurs. Cette approximation est à la base de nombreux code de tracking. Cependant, il est intéressant de noter que le premier terme négligé est le terme octupolaire,

(p2 + p2)2 Δℋ = --x----y-- 4(1 + δ)3

(2.34)

c’est un terme purement cinématique⁶ qui doit néanmoins être pris en compte comme correction quand l’émittance est grande et lorsque la fonction de Twiss $γ$ est importante (voir Papaphilippou et Abell, 2000).

Dans ce cas, il suﬃra soit d’introduire des termes d’ordre supérieur du développement limité de la racine carrée, soit de garder la racine carrée. Ce qui a l’inconvénient majeur de compliquer les calculs ; de plus, un intégrateur du type $𝒮 𝒜 ℬ𝒜n ou 𝒮 ℬ𝒜 ℬn$ avec correcteur ne pourra plus être implémenté, car le correcteur ${{A, B },B }$ n’est en général plus intégrable.

La deuxième approximation fondamentale concerne l’expression générale du champ magnétique (cf. Eq. 1.46 page §) : on suppose que le champ magnétique est constant à l’intérieur d’un élément et nul à l’extérieur, si bien que le Hamiltonien devient autonome pour un élément donné. En réalité, la transition est plus douce. En anglais, on parle d’approximation hard edge, ce qui signiﬁe que les éléments ont un proﬁl magnétique rectangulaire. Cette approximation ne sera plus suﬃsante pour les machines de faible rayon de courbure. Pour prendre en compte les champs de fuite, une solution simple consiste à compléter la description des éléments ; nous en discuterons plus en détail au moment d’aborder les éléments dipolaires.

Les équations du mouvement sont données par les équations de Hamilton :

( dx- ∂H- ( dpx- ∂H- |{ ds = ∂px |{ ds = − ∂x dyds- = ∂∂Hpy dpdys- = − ∂∂Hy- |( dl- ∂H- |( dδ- ∂H- ds = ∂δ ds = − ∂l

(2.35)

Etrangement, bien que certains éléments (dipôle, quadripôles) soient complètement intégrables à la suite des deux approximations précédentes, l’intégration est généralement réalisée en utilisant une solution approchée au moyen d’un intégrateur symplectique. Une des raisons est que cette démarche est nécessaire si l’on désire obtenir l’application de premier retour de l’anneau et utiliser des méthodes automatiques d’algèbre diﬀérentielle (e.g. le DA-Package appelé plus tard Truncated Power Series Algebra — TPSA — développé en FORTRAN77 par Berz en 1989, voir aussi son livre : Modern Map Methods in Particle Beam Physics, 1999 et plus récemment The Full Polymorphic Package qui est une extension écrite en FORTRAN90 par Forest [47]).

Au cours de ce travail, je n’étais pas intéressé — dans un premier temps — par le calcul de l’application de premier retour (analyse globale de l’anneau), c’est pourquoi je conserve la solution exacte si elle existe — sauf temps prohibitifs de calculs, instabilités numériques, rendant l’intégration symplectique plus eﬃcace — De plus, je reste convaincu qu’il est très souvent plus avantageux en termes de temps de calcul d’écrire un code de tracking construit sur le Hamiltonien local de chacun des éléments (au moins pour les machines à rayonnement synchrotron).

2.1.5.1 Section droite

Description et Hamiltonien : La section droite, drift en anglais, est l’élément le plus simple à modéliser dans un accélérateur, car c’est un élément sans champ magnétique ( $As = 0$ ). Elle est caractérisée par un seul paramètre : sa longueur notée $L$ . En coordonnées rectangulaires ( $h = 0$ ), son Hamiltonien se réduit à la forme (cf. Eq. 2.33) :

p2+ p2 ℋ (x,y,l,px,py,δ) = --x----y 2(1 + δ)

(2.36)

Les variables $(x,y,l)$ sont cycliques⁷, les équations du mouvement sont :

( dx px ( dp |{ ds = 1+-δ |{ dxs- = 0 dy = py- dpy- = 0 |( ddsl 1+ δp2x+p2y |( ddsδ ds = − 2(1+δ)2 ds = 0

(2.37)

Une seule hypothèse de calcul : La linéarité des équations est obtenue par le développement limité du terme cinématique (cf. supra) que j’appelle pour la suite approximation des petits angles.

Intégration des équations : Le système 2.37 est complètement intégrable de solution :

(| f i -pix- ( f i |{ x = x + 1+iδs |{ px = p x yf = yi + -py-s pf = pi ||( f i 1(p+ixδ)2+ (piy)2 |( yf y l = l − -2(1+-δ)2--s δ = δ

(2.38)

où les exposants $i$ et $f$ désignent les coordonnées canoniques respectivement à l’entrée et à la sortie de la section droite de longueur $L = s$ .

Prise en compte des termes négligés : Si la racine carrée est conservée, les équations du mouvement sont nonlinéaires. Il est néanmoins encore aisé d’intégrer les équations du mouvement en coordonnées rectangulaires ou curvilignes (cf. annexe A.1 et annexe A.3 expression A.36). Il est ainsi possible de vériﬁer la validité des approximations réalisées.

2.1.5.2 Aimant de courbure parfait

Un dipôle simple est caractérisé par un rayon de courbure constant $ρc$ et une longueur $L$ ; en moyenne, il courbe la trajectoire d’une particule d’un angle $θ = L∕ρc$ . Son Hamiltonien s’exprime naturellement en coordonnées curvilignes en utilisant les équations 2.33 et 1.51 :

p2x + p2y x2 ℋ (x,y,l,px,py,δ) = (1 + hx) --------− hx (1 + δ) + h (x +---) 2(1 + δ ) 2ρc

(2.39)

Par la suite, je ne distinguerai plus le rayon de courbure $ρ c$ de l’élément et celui introduit par le système de coordonnées curvilignes ( $−1 ρ = h$ ). Le terme quadratique en $x$ de l’expression 2.39 traduit une focalisation horizontale purement géométrique du dipôle.

Pour intégrer le Hamiltonien 2.39, je vais me placer dans un cadre un peu plus général. En eﬀet, dans certaines sources de lumière, comme à l’ALS, il existe des dipôles dits combinés, car ils comprennent en plus du terme dipolaire, une composante quadripolaire.

A) Dipôle combiné

Description et Hamiltonien : Le Hamiltonien d’un dipôle combiné est déduit des équations 2.33, 1.51 et 1.53 :

ℋ(x,y,l,p_x,p_y,δ)	= + hx − hx(1 + δ) + h(x + h) + (x² − y²)
	= + hx − hδx + h² + (x² − y²)	(2.40)

où le coeﬃcient $b2$ caractérise le gradient quadripolaire. Dans la plupart des codes de tracking (e.g. BETA [93], DESPOT [42], TRACY [11]), on fait de plus l’approximation dite des grandes machines en négligeant le terme encadré dans l’équation ci-dessus⁸, de type hexapolaire, qui contribue fortement à la chromaticité des petites machines (Dragt, 1982 et Forest, 1998). Si bien qu’au ﬁnal, le Hamiltonien du dipôle combiné est simplement :

p2+ p2 x2 b ℋ (x,y,l,px,py,δ) = --x----y − hδx + h2 ---+ -2(x2 − y2) 2(1 + δ) 2 2

(2.41)

On en déduit les équations du mouvement :

( dx -px- ( dpx 2 |{ ds = 1+δ |{ -ds = hδ − (h + b2)x dy = -py- dpy = b2y |( ddsl 1+δp2x+p2y- |( ddδs ds = − 2(1+δ)2 − hx ds = 0

(2.42)

Hypothèses de calcul : les approximations (a) des grandes machines, (b) des petits angles et (c) hard edge pour le champ magnétique ont été supposées pour exprimer le Hamiltonien du dipôle combiné sous la forme 2.41.

Intégration exacte des équations : Les équations du mouvement 2.42 sont complètement intégrables (cf. équations de deux oscillateurs harmoniques découplés). La solution exacte est donnée en annexe A.2 (p. §) en termes de fonctions trigonométriques circulaires et hyperboliques.

Intégration approchée des équations : Il est également possible d’utiliser un schéma symplectique pour intégrer le Hamiltonien 2.41 que l’on décompose en deux parties séparément intégrables⁹ :

p2+ p2 x2 b ℋ = A + B avec A = --x----y et B = − h δx + h2--+ -2(x2 − y2 ) 2(1 + δ) 2 2

(2.43)

Comme annoncé dans la section précédente, ce découpage du Hamiltonien va nous permettre d’utiliser un intégrateur de classe $𝒮 𝒜 ℬ𝒜n$ ou de classe $𝒮ℬ 𝒜 ℬn$ . Quelque soit le type d’intégrateur choisi, il suﬃt de savoir évaluer les deux opérateurs $esLA$ et $esLB$ . En fait le calcul est presque immédiat, dans chacun des cas, il y a au plus un seul crochet de Poisson à calculer. On obtient les applications respectivement pour $A$ et $B$ :

( i ( ||{ xf = xi + p1x+δs |{ pfx = pix sLA f i piy- f i e :|| y = y + 1+iδ2s i2 i |( py = py ( lf = li − (px2)(1++(δp)y2)s − hxis − h p1x+δ s22 δf = δ

(2.44)

( f i ( f i 2 i |{ x = x |{ px = p x − ((b2 + h )x − hδ)s esLB : yf = yi pfy = piy + b2yis |( f i |( f l = l δ = δ

(2.45)

Nous avons également vu qu’il est possible d’améliorer l’intégrateur en introduisant un correcteur (cf. Eq. 2.31). Il suﬃt de savoir calculer et évaluer le double crochet de Poisson $C={{A,B },B }$ :

1 ( ) {{A, B}, B} = ----- (α + kx )2 + b22y2 1 + δ

(2.46)

avec $α=− δh$ et $2 k = b2 + h$ . L’application déﬁnie pour le correcteur est alors :

{ f i { f i 2k(α+kxi) esL {{A,B },B} : x = x px = px − 1+ δ s yf = yi pf = pi − 2b22yis y y 1+δ

(2.47)

avec $s=− s3c2$ conformément à l’équation 2.31.

Correction due aux coins de l’aimant :

pict

FIG. 2.2:

Schéma d’un aimant secteur (a) et d’un aimant à faces parallèles (b). Les angles des coins d’entrée et de sortie sont ici égaux : $θ θe = θs = 2$ .

Pour être rigoureux, l’expression du Hamiltonien du dipôle combiné 2.41 décrit un aimant secteur ou aimant à faces tournées (cf. la déﬁnition du bloc courbe page § et le schéma 2.2-a). Souvent, les aimants d’un accélérateur sont à faces parallèles (e.g. à l’ALS ou Super-ACO), il est donc nécessaire de corriger les eﬀets de bord introduits par les coins des aimants (cf. Fig. 2.2-b). Au premier ordre, ces eﬀets se modélisent par un champ quadripolaire focalisant horizontalement et défocalisant verticalement¹⁰. L’application entre l’entrée et la sortie d’un aimant à faces parallèles est simplement — au premier ordre — la composition des applications d’un aimant secteur, $ℳsecteur$ , et des coins d’entrée, $𝒯 e$ , et sortie, $𝒯 s$ , $xf = 𝒯 sℳsecteur𝒯 exi$ (Forest et al., 1994), soit pour un angle d’entrée ou de sortie de l’aimant $θu$ :

{ f i i 𝒯 (θ ) : px = px + hx tan θu avec u = e,s u u pfy = piy − hyitan θu

(2.48)

Notons que cet eﬀet sera prépondérant pour les petites machines, car il est proportionnel à la courbure $h$ .

Champs de fuite :

pict

FIG. 2.3:

Schéma du proﬁl magnétique longitudinal d’un aimant de longueur $L$ . En approximation hard-edge, le champ magnétique est constant dans l’élément et nul à l’extérieur. En réalité, le champ magnétique décroît jusqu’à une valeur nulle de part et d’autre sur une longueur $λ$ : on parle de champ de fuite.

Pour les machines à faible rayon de courbure comme Super-ACO, un eﬀet supplémentaire est encore à ajouter dans le modèle : les champs de fuite. Ce phénomène est simplement lié au fait que le champ magnétique décroît sur les bords de l’aimant et n’est pas tout à fait nul juste à l’extérieur du dipôle (l’approximation hard-edge doit être complétée). Au premier ordre, les champs de fuite induisent une composante quadripolaire verticale proportionnelle à la courbure $h$ de l’anneau. Le déphasage vertical¹¹ ( $ψu$ ) est donné par la formule (Brown, 1982, pp. 116–117) :

2 1 +-sin-θu ψu = − 𝒦hg cosθu avec u = e,s

(2.49)

avec $g$ la distance entre les deux pôles de l’aimant, $h$ la courbure, $θu$ l’angle d’entrée (e) ou de sortie (s) et $𝒦$ , l’intégrale :

∫ +∞ By (s)[B0 − By(s)] 𝒦 = ----------2-------ds −∞ gB 0

(2.50)

où $B(s) y$ est l’amplitude du champs de fuite sur le plan moyen à la longitude $s$ mesurée depuis l’entrée de l’aimant et $B0$ est la valeur asymptotique de $By (s)$ dans l’aimant (cf. Fig.2.3).

L’application symplectique $ℱ u$ du premier ordre incluant les coins et le champ de fuite de l’aimant :

{ pf = pi+ hxitan θ ℱ u : xf xi i u avec u = e,s py = py − hy tan(θu − ψu )

(2.51)

B) Dipôle simple

Pour un dipôle simple, il suﬃt de poser $b2 = 0$ dans les formules établies pour le dipôle combiné. En annexe A.2, les équations du mouvement sont intégrées pour diﬀérentes approximations et méthodes d’intégration : (a) sans approximations des grandes machines et des petits angles en géométrie curviligne et rectangulaire et (b) en prenant en compte le terme des petits machines.

2.1.5.3 Quadripôle droit

Description et Hamiltonien : Un quadripôle droit est un élément magnétique dont la vocation est de focaliser la particule qui le traverse ; il a une longueur $L$ et un gradient magnétique $b2=K$ ; son Hamiltonien s’exprime naturellement en géométrie cartésienne à partir des expressions 2.33 et 1.53 :

p2+ p2 K ℋ (x,y,l,px,py,δ) = --x----y + --(x2 − y2) 2(1 + δ) 2

(2.52)

Les équations du mouvement du quadripôle sont ( $l$ est cyclique) :

( dx px- ( dpx- |{ ds = 1+ δ |{ ds = − Kx ddys = p1+yδ dpy- = Ky |( dl -p2x+p2y |( ddsδ ds = − 2(1+δ)2 ds = 0

(2.53)

Hypothèses de calcul : L’expression du Hamiltonien 2.52 est établie dans (a) l’approximation des petits angles et (b) l’approximation hard-edge.

Intégration « exacte » : Les équations du mouvement 2.53 s’intègrent directement pour une longueur $s$ et pour $K > 0$ , on obtient :

(f i --1--- i ( f i i |{x = cos(ωs)x + ω(1+ δ) sin(ωs )px |{ px = − ω(1 + δ)sin(ωs )x + cos(ωs)px yf = cosh(ωs)yi + --1---sinh (ωs)piy pfy = ω (1 + δ) sinh (ωs )yi + cosh (ωs)piy |(f i ω(1+δ) |( f l = l + Δl δ = δ

(2.54)

avec $∘ ---- -K- ω= 1+δ$ , $s = L$ et

Δl =	ω (xⁱ)² −ω (yⁱ)²
	− (p_xⁱ)² − (p_yⁱ)²	(2.55)
	xⁱp_xⁱ −yⁱp_yⁱ

Cette application de transfert décrit un quadripôle focalisant dans le plan horizontal et défocalisant dans le plan vertical. Dans le cas $K < 0$ , il suﬃt de poser $K = − K$ et d’inverser les fonctions circulaires et hyperboliques dans les équations 2.54. Le quadripôle est alors focalisant dans le plan vertical.

Intégrateur symplectique : Les solutions s’écrivent comme celles du dipôle combiné (cf. Eq. 2.44 à 2.47) en posant une courbure nulle, i.e. $h = 0$ .

2.1.5.4 Hexapôle parfait droit

Description et Hamiltonien : Les hexapôles sont inévitablement introduits pour corriger la chromaticité. Le Hamiltonien d’un hexapôle de force $b3 = S$ et de longueur $L$ s’exprime en utilisant les équations 2.33 et 1.53 en coordonnées rectangulaires par :

p2 + p2 S ℋ (x,y, l,px,py, δ) = -x----y-+ --- (x3 − 3xy2 ) ◟2(1 +◝◜-δ)◞ ◟3----◝◜ -----◞ A(px,py,δ) B(x,y)

(2.56)

Les équations du mouvement sont ( $l$ est cyclique) :

( ( |{ ddxs = 1p+xδ |{ dpx = − S (x2 − y2) dy -py ddpsy | ds = 1+δ2 2 | ds = 2Sxy ( ddls = − p2x(1++pδy)2 ( ddδs = 0

(2.57)

Hypothèses de calcul : L’expression du Hamiltonien 2.56 est établie dans (a) l’approximation des grandes machines, (b) l’approximation hard-edge. Le schéma d’intégration que nous allons présenter pour l’hexapôle peut simplement se généraliser pour modéliser un multipôle droit ou tourné quelconque. Par exemple, un $2n$ -pôles droit de force $bn$ est modélisable par un Hamiltonien de la forme :

p2 + p2 bn ℋ (x,y, l,px,py, δ) = -x----y-+ --Re ((x + jy)n) 2(1 + δ) n

(2.58)

Intégration approchée des équations : La partie $A (px,py,δ)$ contient le Hamiltonien d’une section droite dont la solution a déjà été calculée (cf. Eq. 2.38). Il ne reste plus qu’à intégrer la partie $B(x, y)$ contenant le potentiel hexapolaire, on obtient l’application :

( f i ( f i 2 2 |{ x = x |{ px = px − S(x − y ) s esLBxi : yf = yi pf = pi + 2Sxy s |( f i |( yf y l = l δ = δ

(2.59)

avec $s=L$ . Usuellement, l’hexapôle est modélisé par une lentille mince de longueur inﬁnitésimale, e.g. $−6 L=10$ m. Dans ce cas, il suﬃra de prendre un intégrateur symplectique d’ordre 2. Si l’hexapôle est modélisé par une lentille épaisse, un intégrateur d’ordre 4 est utilisé comme pour les éléments magnétiques précédents.

2.1.6 Etude comparative

2.1.6.1 Introduction

Un code d’intégration numérique a été écrit en langage FORTRAN90. Ce code prend en compte l’ensemble des modèles d’éléments magnétiques précédemment présentés. Ce programme est écrit de manière modulaire aﬁn de pouvoir traiter les petites et grandes machines ; il est possible d’introduire de nouveaux éléments magnétiques (e.g. les dispositifs d’insertion) ; le code est dédié et optimisé pour le tracking. Les intégrateurs avec ou sans correcteur de type $𝒮𝒜 ℬ 𝒜n$ et $𝒮ℬ 𝒜ℬn$ ont été programmés pour $n = 1, 2$ et $3$ . Des intégrateurs d’ordre plus élévé pourront sans diﬃculté être introduits par la suite. Pour les comparaisons internes, l’intégrateur d’ordre 4 de Forest et Ruth a également été programmé.

Diﬀérents types de comparaisons ont été eﬀectués, en particulier avec les codes DESPOT et MAD respectivement sur une maille de l’ALS (grande machine) et de Super-ACO (petite machine). Les calculs ont été réalisés en double précision sur une station DIGITAL PWS 433 AU (EV56 à 433 MHz).

2.1.6.2 Propriétés

Précision : Nous avons vériﬁé la précision de l’intégration pour chaque élément magnétique modélisé. Le Hamiltonien est une intégrale première du mouvement. Nous présentons quelques résultats pour le dipôle combiné dans l’approximation des grandes machines. Son Hamiltonien est rappelé :

p2x + p2y x2 b2 ℋ (x, y,l,px,py,δ) = --------− hδx + h2--- + --(x2 − y2 ) 2◟(1◝+◜-δ)◞ ◟--------2-◝◜ -2--------◞ A B

avec comme courbure $h = 0.2015 m −1$ , gradient quadripolaire $b2 = − 4 × 10−3m − 2$ et longueur $L=0.86m$ . Le dipôle est ainsi focalisant dans les deux plans.

Tout d’abord, il est nécessaire de remarquer que le Hamiltonien s’apparente à celui de deux oscillateurs. Les parties $A$ et $B$ sont comparables i.e. le petit paramètre $ε$ est voisin de l’unité. Donc pour cet élément, il sera inutile de vouloir utiliser un intégrateur d’ordre supérieur à 4 puisque la nouvelle classe d’intégrateurs n’a d’intérêt qui si $ε$ est un petit paramètre.

Les erreurs sur l’énergie du dipôle combiné sont calculées sur un temps d’intégration de 1 000 passages à travers l’aimant (Fig. 2.4). Pour l’intégration exacte (cf. annexe A.2.4, page §), l’erreur est aléatoire avec une dérive numérique linéaire du nombre d’itérations comme attendu (la précision machine est $2.26 × 10−16$ ).

Les trois schémas symplectiques choisis sont le schéma de Ruth et le schéma $SABA 2$ et $SBAB2$ avec correcteur. Le dipôle de longueur $L$ est intégré en trois étapes d’intégration ( $k=3$ ). La meilleure des méthodes est l’intégrateur $SABA2$ qui est plus d’un ordre de grandeur plus précis que le schéma de Ruth (Fig. 2.4).

pict

FIG. 2.4:

Comparaison à coût constant des intégrateurs $SABA2$ , $SBAB2$ avec correcteur et Ruth, pour un dipôle combiné. Erreur sur l’énergie $ΔE ∕E$ en fonction du nombre $N$ de passages dans l’aimant. Pour les schémas symplectiques l’erreur est bornée et périodique.

pict

FIG. 2.5:

Erreur relative sur l’énergie en échelle logarithmique en fonction de la taille du pas d’intégration $log10(s)$ . Les intégrateurs $SABA 2$ , $SBAB 2$ avec correcteur sont plus précis que l’intégrateur de Forest et Ruth respectivement d’un ordre et d’un demi ordre de grandeur.

Une étude systématique de la précision des méthodes a également été réalisée. Les intégrateurs $SABAn$ et $SBABn$ avec correcteur sont comparés à coût constant (i.e. à nombre égal d’évaluations) avec le schéma de Ruth. L’erreur relative des intégrateurs symplectiques est présentée pour $n = 2$ en fonction du nombre d’étapes d’intégration pour un aimant de longueur $L$ (Fig. 2.5). En échelle logarithmique, la pente des droites est 4, ce qui correspond à l’ordre des méthodes utilisées. La méthode $SABA 2$ avec correcteur est plus précise que le schéma de Ruth d’un ordre de grandeur.

Ce résultat va pouvoir être utilisé pour diminuer, à précision constante, le temps d’intégration des élément magnétique, soit de réduire le temps de tracking. Généralement, une à deux étapes d’intégration sont gagnées par comparaison avec l’intégrateur de Forest et Ruth.

Déphasage :


	k=4		k=6		k=8

Méthode	$N$	$νx$	$N$	$νx$	$N$	$νx$

« Exact »	-	0.185450	-	0.185450	-	0.185450
Ruth	5 441	0.185267	27 278	0.185414	87 163	0.185439
$SABA2$	330 503	0.185453	1 664 687	0.185451	5 256 904	0.185450

TAB. 2.3:

Comparaison du nombre $N$ nécessaire d’itérations pour obtenir une solution déphasée de $2π$ par rapport à l’intégration de référence d’une maille de l’ALS. Le déphasage et le nombre d’ondes sont exprimés pour diﬀérentes valeurs du nombre $k$ d’étapes d’intégration. Le schéma $SABA2$ est à chaque fois plus précis de plus de un à deux ordres de grandeur.

Un autre point remarquable de cette nouvelle classe d’intégrateurs est sa faible dérive en phase à nombre égal d’itérations. En comparant la solution numérique obtenue par le schéma symplectique par rapport à la solution exacte, on note que le déphasage introduit par la méthode $SABA2$ est très faible par rapport à l’intégrateur de Ruth et Forest (voir Fig. 2.6 pour le dipôle combiné et Fig. 2.7 pour une maille complète de l’ALS). Par exemple, le tableau 2.3 donne le nombre $N$ de tours de l’ALS pour obtenir un déphasage $Δ$ de $2π$ pour diﬀérentes valeurs du nombre d’étapes $k$ d’intégration des dipôles et des quadripôles de l’anneau¹². Une loi d’échelle entre le temps d’intégration $T = N × L$ et le nombre d’étapes $k$ d’intégration peut être établie. Si $L$ est la longueur de l’élément intégré $N$ fois, et si $k$ est le nombre d’étapes pour intégrer l’aimant. Alors le pas de l’intégrateur pour chaque étape est $L∕k$ et on établit :

( ) k2 4 N2 --- = --- k1 N1

où $N1$ , $N2$ correspondent au nombre de fois qu’il faut intégrer l’aimant de longueur $L$ , respectivement en $k1$ et $k2$ étapes, pour obtenir un déphasage $Δ = xsymp(T ) − xexact(T ) = 2π$ entre le schéma symplectique ( $xsymp$ ) et exact ( $xexact$ ). En eﬀet, en se souvenant que la méthode est d’ordre 4 et en ne gardant que les termes de plus bas degré, on peut écrire :

( L )4 Δ = |xsymp (T ) − xexact(T)| = T × -- = Cte k

Notons que le déphasage n’est pas une grandeur fondamentale pour déterminer si un intégrateur numérique est meilleur qu’un autre. Nous avons vu précédemment qu’un intégrateur symplectique peut être vu comme un modèle d’un élément magnétique. Ainsi pour retrouver les bons nombres d’ondes, il suﬃrait d’ajuster les conditions initiales.

pict

FIG. 2.6:

Déphasage $Δ$ des solutions du dipôle combiné introduit par les schémas $SABAC2$ et de Ruth par rapport à la solution « exacte ». Le déphasage est inférieur de deux ordres de grandeur pour la nouvelle classe d’intégrateur.

pict

FIG. 2.7:

Déphasage introduit par les schémas symplectiques pour l’intégration sur $N$ tours de l’ALS. Les quadripôles et dipôles sont intégrés en $k = 4$ étapes. Dès 5 441 tours, le schéma de Ruth introduit un déphasage complet (remarquer la diﬀérence d’échelle verticale pour les deux méthodes).

Jusqu’à présent, l’analyse de la précision des méthodes symplectiques n’a pas été réalisée. Cependant, une remarque concernant la taille des pas d’intégration peut être faite. Pour cela, décrivons succinctement le schéma d’intégration numérique d’un élément de longueur $L = 1$ avec l’intégrateur $SABA1$ dont l’expression est rappelée¹³ :

xf = ec1LA ed1LB ec1LAxi ◟-◝◜i-◞ ◟----◝◜-xc1-◞ xi d1

(2.60)

avec $2c1= 1$ et $d1 = 1$ . Dans ce cas l’intégration est réalisée en deux pas suivant $A$ et un pas suivant $B$ . Si l’intégration est réalisée en $k$ étapes, le schéma d’intégration 2.60 est composé $k$ fois avec un pas $L = L∕k$ .

pict

FIG. 2.8:

Schéma pour l’intégration des parties $A$ et $B$ d’un Hamiltonien avec (a) l’intégrateur de Forest et Ruth et (b) l’intégrateur $SABA1$ .

Ce processus est illustré pour les intégrateurs $SABA2$ et de Forest et Ruth (schéma 2.8 pour $k=1$ ). Pour l’intégrateur $SABA2$ les pas, tous positifs, sont petits, plus ﬁns sur les bords de l’élément et plus grands en son centre (cf. condition 2.20). L’intégrateur de Forest et Ruth est caractérisé par la présence de grandes valeurs de pas. Pour intégrer $A$ , on fait deux pas positifs $2×0.6756 = 1.3512$ et deux pas négatifs $2 × − 0.1724 = − 0.3512$ ; pour intégrer $B$ , deux pas positifs $2× 1.3512 = 2.7024$ et un grand pas négatif $− 1.7024$ .

La valeur des pas d’intégration a probablement une incidence sur la taille du reste des méthodes symplectiques et donc sur leur précision et phase. Des études plus complètes pourront être entreprises pour préciser ces phénomènes.

Conclusion : Dans les approximations réalisées, il est nécessaire de trouver un compromis entre la précision et le temps de calcul. L’utilisateur peut selon les besoins modiﬁer la précision de calcul en jouant sur le nombre d’étapes d’intégration d’un élément de longueur $L$ . La nouvelle classe d’intégrateurs utilisée permet de d’obtenir un intégrateur plus eﬃcace que le schéma classique de Forest et Ruth à précision de calcul égale.

2.2 Analyse en Fréquence

2.2.1 Introduction

Longtemps, le mouvement du Système Solaire a semblé être d’une régularité parfaite. Pour étudier sa dynamique globale, Jacques Laskar a développé l’Analyse en Fréquence qui lui a permis de mettre en évidence du chaos dans le Système Solaire (Laskar, 1988, 1990). Cette méthode numérique utilise une méthode de Fourier raﬃnée ; elle est adaptée pour l’analyse de la stabilité des orbites d’un système conservatif. Cette méthode a ensuite été appliquée à l’étude da la stabilité de l’obliquité des planètes puis de la Terre (Laskar et Robutel, 1993a et 1993b). Rapidement, elle a été utilisée pour caractériser la dynamique d’autres systèmes : application standard (Laskar, Froeschlé et Celleti, 1992 – Carletti et Laskar, 2000), dynamique galactique (Papaphilippou et Laskar, 1996 et 1998), dynamique des accélérateurs (Dumas et Laskar, 1993 – Laskar et Robin, 1996 – Papaphilippou, 1999).

L’Analyse en Fréquence permet de caractériser ﬁnement la dynamique d’un système à $n$ degrés de liberté, est particulièrement adaptée aux systèmes à 3 degrés de liberté. Pour des orbites régulières, la méthode donne une représentation analytique des solutions alors que pour des orbites non régulières, elle fournit des critères permettant de caractériser le degré de nonlinéarité du système.

2.2.2 Application fréquence

Pour l’écriture de cette partie, je me suis inspiré des articles de J. Laskar décrivant de manière exhaustive la méthode d’Analyse en Fréquence (1992, 1993, 1994 et 1999).

Pour décrire la méthode d’Analyse en Fréquence, nous allons d’abord nous placer dans le cas général d’un système Hamiltonien autonome à $n$ degrés de liberté écrit sous la forme :

ℋ = ℋ0 + εℋ1

(2.61)

où $ℋ0$ représente la partie intégrable du mouvement et $ℋ1$ un terme de perturbation — dans le cadre de la dynamique d’un accélérateur, nous avons vu au cours du chapitre 1 (page §) que $ℋ0$ décrit le mouvement bétatron et $ℋ1$ décrit par exemple les défauts multipolaires —

Supposons qu’en l’absence de perturbation ( $ε = 0$ ), le système puisse s’écrire en variables actions-angles $(I,θ)$ . Alors les équations du mouvement se réduisent à :

{ dθk ∂ℋ0- { dt = ∂Ik = νk(I) ⇒ θk(t) = νk(I)t + θ0k avec k = 1,...,n ddItk= 0 Ik(t) = Ik0

(2.62)

Les solutions 2.62 peuvent s’écrire sous forme complexe $jθ zk = I0ke k$ , soit :

jνkt zk(t) = z0ke , avec z0k = zk(0 ) et k = 1,...,n

(2.63)

Les orbites 2.63 sont conﬁnées sur des tores de dimensions $n$ , produits de cercles de rayons $I 0k$ parcourus aux fréquences $ν(I ) = (ν (I ),...,ν (I )) 0 1 0 n 0$ (cf. Fig. 2.9 pour $n = 2$ ).

pict

FIG. 2.9:

Espace des phases d’un système à deux degrés de liberté en coordonnées actions-angles $(I ,I ,θ ,θ ) 1 2 1 2$ . Les orbites ( $x i$ ) sont conﬁnées sur des tores de dimension 2.

S’il existe une bijection entre les actions et les fréquences, i.e. si la condition de non-dégénérescence suivante est vériﬁée :

( ) ( 2 ) det ∂ν-(I)- = det ∂-ℋ0-(I)- ⁄= 0 ∂I ∂I2

(2.64)

alors le mouvement peut être décrit de manière équivalente par les actions $Ik$ ou les fréquences $νk$ ; on déﬁnit alors l’application, $F$ , dite application fréquence :

F : (I1,...,In) ↦→ (ν1,...,νn )

(2.65)

Si de plus, on se place sur une surface d’énergie constante ( $ℋ0 (I ) = Σ$ ), alors seulement $(n−1)$ actions sont indépendantes et l’application fréquence 2.65 devient (si par exemple, $νn⁄=0$ ) :

F : (I ,...,I ) ↦→ (ν1,..., νn−1) 1 n− 1 νn νn

(2.66)

Le problème qui nous préoccupe est de savoir ce qu’il advient de ces orbites sous l’eﬀet d’une perturbation faible. La réponse est fournie par le théorème KAM (Kolmogorov-Arnold-Moser, voir par exemple les références incluses in Laskar, 1999) : sous des conditions très générales, pour une perturbation suﬃsamment faible, la plupart des tores du système non perturbé ( $ℋ0$ ) subsistent mais sont déformés. Ce sont les tores dont les vecteurs fréquences vériﬁent la condition dite diophantienne :

-Cε-- n | < m, ν > | ≥ |m |τ, τ > n et ∀ m ∈ ℤ

(2.67)

où $Cε$ est une constante dépendant de $ε$ et de $ℋ0$ : ils sont appelés tores KAM. Les tores de $ℋ0$ pour lesquels $< m, ν >= 0$ sont appelés tores résonants et sont détruits d’après le théorème KAM. Entre les tores KAM, les orbites sont en général chaotiques. Cependant, il est encore possible de construire une application fréquence (voir Laskar, 1999). Les solutions KAM peuvent alors s’exprimer sous la forme :

∑ zk(t) = zk0ejνkt + akmiejωmit, mi ∈ ℤn mi

(2.68)

Les coeﬃcients $ωmi$ sont combinaison linéaire de $n$ fréquences indépendantes appelées fréquences fondamentales du système, i.e. $ωmi = < mi, ν >= mi1ν1 + ...+ mn νin$ .

La méthode d’Analyse en Fréquence repose sur la construction de l’application fréquence numérique, $FT$ , en recherchant une décomposition quasi-périodique d’une trajectoire sur un temps d’intégration ﬁni $T$ et en utilisant une technique de Fourier raﬃnée (Laskar 1988, 1993 et 1999), i.e. sous la forme d’un nombre ﬁni $N$ de termes :

N jνkt ∑ k j<mi,ν>t n zk(t) = zk0e + amie , mi ∈ ℤ i=1

(2.69)

où les coeﬃcients $akmi$ sont ordonnés par amplitude décroissante¹⁴ et $ν$ est le vecteur de fréquences fondamentales.

A énergie ﬁxée, une trajectoire de l’espace des phases vit sur une surface d’énergie de dimension $2n − 1$ . Habituellement, on observe le mouvement dans une surface de Poincaré, i.e. que l’on enregistre une trajectoire discrète en ﬁxant un des angles $θ k$ , par exemple $θn=0$ modulo $2π$ (cf. Fig. 2.10). Les fréquences fondamentales déduites du signal quasi-périodique 2.69 sont exprimées par rapport à la section de Poincaré, donc par exemple dans notre cas $νn = 1$ . Le mouvement est donc restreint sur un espace de dimension $2n − 2$ .

pict

FIG. 2.10:

Section de Poincaré (a) : à chaque tour de l’accélérateur, la trajectoire discrète de la particule est enregistrée dans le plan $θ = 0 n$ . (b) Espace des phases $√--- √---- (xk=2Ik cosθk, pk = 2Ik sin θk)$ au temps (i).

Pour l’Analyse en Fréquence, une condition plus stricte est adoptée : on ﬁxe les $n$ angles $θk=θ0k$ . Dans ce cas, le mouvement est étudié dans un espace de dimension $n − 1$ . Il est caractérisé par la donnée de $n − 1$ actions.

En ne conservant que les fréquences fondamentales de la décomposition quasi-périodique (Eq. 2.69) obtenue après intégration sur un intervalle de temps $[τ,τ + T ]$ , on construit l’application fréquence numérique :

F_θ₀^T : ℝ × ℝⁿ⁻¹	−→ℝⁿ⁻¹
(τ,I₁,…,I_n−1)	(ν₁,…,ν_n−1)	(2.70)

2.2.3 Propriétés de l’application fréquence numérique

L’étude des propriétés de l’application $F T$ va nous donner des informations sur le système dynamique (résonances, nonlinéarités, chaos).

L’image $T F θ0(τ,I)$ est indépendante du choix dans angles initiaux $θ0$ (cf. sur un tore, les fréquences ne dépendent pas de la phase) et $F T$ est une fonction régulière pour une trajectoire KAM. Inversement, la non-régularité de $F T(τ,I)$ nous permettra de mettre en évidence les phénomènes résonants et les zones chaotiques.
La convergence de la méthode est rapide et suit une loi en $1T4$ en utilisant une fenêtre de Hanning alors qu’une transformée de Fourier rapide (FFT) converge en $1- T$ (Laskar, 1999). En eﬀet, l’algorithme de la FFT présuppose que le signal recherché est périodique, alors que l’Analyse en Fréquence cherche à reconstruire une plus grande classe de signaux, à savoir les signaux quasi-périodiques de la forme (2.69) ; ce point fondamental rend la méthode plus apte à l’étude d’un système dynamique. Cette propriété permet de la calculer sur un temps d’intégration $T$ court et de faire de la prédiction à long terme.
L’application $T F θ0$ est invariante par translation dans le temps $τ$ pour une trajectoire régulière (trajectoire KAM). A l’inverse, la variation avec le temps des fréquences, appelée diﬀusion, pourra être utilisée comme indice de stabilité et de chaoticité de la trajectoire.
La diﬀusion des orbites nous donnera de précieuses informations sur la stabilité globale de la dynamique du système étudié. Pour un système à deux degrés de liberté, l’espace des phases est de dimension 4, les surfaces d’énergie de dimension 3 et l’espace des fréquences de dimension 1 (cf. réduction par la conservation de l’énergie). Une orbite KAM est représentée par un point dans l’espace des fréquences et vit sur un tore de dimension 2. Une orbite non régulière n’est pas un point mais une courbe révélant la diﬀusion entre les orbites régulières (cf. Fig. 2.11-a). Il est clair que dans l’espace des fréquences, un point (dimension 0) sur une courbe de dimension 1 empêche tout mouvement régulier d’avoir de grandes excursions en amplitude. De manière équivalente, en se plaçant dans l’espace des phases, un tore de dimension 2 sépare une surface d’énergie (dimension 3) en une partie intérieure et une partie extérieure : c’est donc une barrière topologique à la diﬀusion.
Pour un système à trois degrés de liberté, la situation n’est plus du tout la même. Un point n’arrête pas la diﬀusion dans un plan (cf. Fig. 2.11-b). Vus dans l’espace des phases, les tores de dimension 3 ne sont plus des frontières topologiques à la diﬀusion pour une surface d’énergie de dimension 5. Cependant il peut être montré que la diﬀusion est quasi-nulle au voisinage des orbites régulières. Au voisinage d’une trajectoire KAM, la diﬀusion décroît exponentiellement si bien qu’un grand nombre de trajectoires régulières joue le rôle de barrière eﬀective à la diﬀusion (résultat de Morbidelli et Gorgilli, 1995).

pict

FIG. 2.11:

Diﬀusion pour un système à deux degrés de liberté (a), e.g. l’application d’Hénon (cf. infra) : la diﬀusion est conﬁnée par les trajectoires régulières qui sont des points dans l’espace des fréquences. Pour trois degrés de liberté (b), e.g. la dynamique transverse d’un accélérateur, la diﬀusion n’est plus conﬁnée par les trajectoires régulières ; la diﬀusion est cependant faible en leur voisinage (Laskar, 1994)

2.2.4 Applications préliminaires

Avant d’appliquer l’Analyse en Fréquence à l’étude de la dynamique des accélérateurs, deux applications simples vont être présentées succinctement. Le pendule mécanique qui modélise la dynamique au voisinage d’une résonance. Puis, l’application d’Hénon qui est une application nonlinéaire et me permettra d’introduire quelques notions supplémentaires sur l’étude de la dynamique des systèmes.

2.2.4.1 Le pendule

Le Hamiltonien d’un pendule rigide peut s’écrire sous la forme :

2 ℋ = p--− a cosq, a ∈ ℝ+ 2

(2.71)

où $q$ est l’angle de rotation du pendule et $p$ sa dérivée temporelle. Le pendule est un système dynamique à un degré de liberté complètement intégrable qui permet de modéliser la dynamique au voisinage d’une résonance (voir par exemple, Chirikov, 1979). Son portrait de phase $(p,q)$ est tracé pour diﬀérentes valeurs de l’énergie (Fig. 2.12-a). Deux régimes peuvent être mis en évidence :

le régime de libration ( $ℋ < a$ ) où le pendule eﬀectue des oscillations autour de l’origine de période :
$∫ arccosℋa- dq T(ℋ ) = ∘--------------- − arccosℋa 2(ℋ + acos q)$ (2.72)
le régime de circulation où l’énergie du pendule ( $ℋ > a$ ) est assez grande pour lui permettre d’eﬀectuer des mouvements de rotation de période :
$∫ 2π dq T (ℋ) = ∘--------------- 0 2 (ℋ + a cosq)$ (2.73)

pict pict pict

FIG. 2.12:

Espace des phases $(q,p)$ du pendule rigide (a) et courbes en fréquence $ν (q)$ . Au voisinage du point hyperbolique (b), la courbe en fréquence présente une singularité logarithmique. Au voisinage du point elliptique (c), la courbe en fréquence est nulle (plateau). Après échantillonnage en moment $p$ , les projections (d et e) des courbes en fréquences sont identiques pour le régime de circulation.

Le pendule admet deux points ﬁxes correspondant aux positions d’équilibre stable et instable : l’origine $(p,q) = (0,0)$ est un point ﬁxe elliptique et $(p,q) = (0,π )$ est un point ﬁxe hyperbolique.

La courbe en fréquence $ν (p) = 2π- π T$ est tracée pour $q = π$ , i.e. passant par le point hyperbolique et y admet une singularité logarithmique. La période de rotation (Eq. 2.73) est alors (voir par exemple Tabor, 1988 p. 11 sqq. et Laskar, 1993) :

2 ( π ∘ ---4a--) T(ℋ ) = ∘----------sn --, -------- , p2 + 4a 2 p2 + 4a

(2.74)

avec $sn(x,k)$ , la fonction sinus elliptique de Jacobi :

∫ x dφ sn(x,k) = ∘-------------- 0 1 − k2 sin2 φ

(2.75)

La courbe en fréquence $ν0(p)$ calculée pour $q = 0$ est identiquement nulle pour $2 p < 4a$ , car la fréquence de circulation est nulle. Par contre au-delà, elle vaut $∘ -------- ν0(p) = νπ( p2 − 4a)$ .

Ces deux courbes sont fondamentales, car elles décrivent le comportement de l’application fréquence au voisinage d’une résonance (Laskar, 1993). La courbe $ν (p) π$ décrit son comportement lorsqu’un point hyperbolique est traversé (Fig. 2.12-b), la courbe $ν0(p)$ lorsqu’une île de résonance est traversée (Fig. 2.12-c). Si l’on projette la courbe $ν (q )$ sur l’axe des fréquences, pour un échantillonnage uniforme en $p$ , l’allure de la courbe en fréquence ne dépend pas du choix de $q0$ dans le régime libration (cf. courbes 2.12-d et 2.12-e et propriété (1) p. §). Par contre l’échantillonnage est plus ﬁn au voisinage de île de libration qu’au voisinage du point hyperbolique, car la variation de la fréquence est plus rapide (cf. singularité logarithmique).

2.2.4.2 L’application d’Hénon

L’application d’Hénon (Hénon et Heiles, 1964 et Hénon, 1969) est la plus simple des applications non triviales symplectiques polynomiales à deux degrés de liberté (Bazzani et al., p. 77) :

(q ′) ( cosω sin ω ) ( q ) ′ = 2 p − sinω cosω p + εq

(2.76)

où $p$ et $q$ sont des variables canoniques et $ω$ est le nombre de rotation. Elle dérive du Hamiltonien dépendant explicitement du temps $s$ :

ω ε ℋ (p,q;s) = --(q2 + p2) +--q3δ2π(s) 2 3

(2.77)

avec $δ 2π$ la fonction de Dirac $2π$ -périodisée.

Malgré sa simplicité, l’application d’Hénon contient une grande partie de la dynamique nonlinéaire que nous rencontrerons dans les chapitres suivants. La partie quadratique du Hamiltonien 2.77 modélise le mouvement bétatron perturbé par un hexapôle situé en $s = 0$ (cf. Wiedemann, Tome II ou Lee, 1998). Je ne l’introduis ici que pédagogiquement pour illustrer les propriétés de l’application fréquence énoncées dans la section 2.2.3 et également pour discuter des résonances induites par une perturbation de type hexapolaire (pour plus de détails, voir le travail de Bazzani, Todesco, Turchetti et Servizi, 1994).

Dans ce cas, l’application fréquence est simplement déﬁnie par :

F^T : ℝ	→ ℝ
q	ν	(2.78)

Les portraits de phases de l’application d’Hénon sont tracés pour diﬀérentes valeurs de la fréquence linéaire $ω$ et pour un nombre d’itérations $T = 2000$ (Fig. 2.14) et mettent respectivement en évidence des résonances d’ordre 4, 5, 6 et 7 (cf. grands îlots de résonance). Pour chacun des cas, la courbe en fréquence $ν (q)$ pour $p = 0$ et $q > 0$ est adjointe. C’est une courbe monotone de l’amplitude $p$ excepté au voisinage des résonances. L’étude de la régularité de $ν(q)$ révèlent également des résonances d’ordre plus élevé.

Similairement au cas du pendule, lorsque $ν (q )$ traverse une île, on observe un plateau, et au voisinage d’un point hyperbolique une discontinuité de l’application fréquence. L’application d’Hénon est suﬃsamment simple pour pouvoir calculer analytiquement la variation de $ν$ avec $q$ . Pour cela, on écrit une forme normale de l’application d’Hénon au voisinage de l’origine qui est un point ﬁxe elliptique (voir Bazzani et al., 1994). Les largeurs de résonance peuvent aussi être calculées et dépendent à la fois de $q$ et de $ε$ . A l’origine, on retrouve la fréquence linéaire, i.e. $ν(0)= ω$ .

pict

FIG. 2.13:

Perturbation d’un système : destruction d’un tore rationnel de fréquence $ν = 1 4$ en quatre points ﬁxes elliptiques et quatre points ﬁxes hyperboliques. Les tores irrationels sont faiblement déformés. Au voisinage de chaque point elliptique, le raisonnement peut être réappliqué.

pict

FIG. 2.14:

Portraits de phase $(q,p)$ et courbes en fréquence associées $ν(q)$ pour $p = 0$ de l’application d’Hénon. La fréquence linéaire $ω$ a été choisie dans chacun des cas pour mettre en évidence des résonances d’ordre 4 (a, $ω = 0.251$ ), 5 (b, $ω = 0.205$ ), 6 (c, $ω = 0.180$ ) et 7 (d, $ω = 0.150$ ). Des résonances d’ordre plus élevé apparaissent et limitent l’ouverture dynamique.

L’apparition de résonances de tout ordre sous l’action d’une perturbation, me permet d’introduire un résultat général sur la dynamique d’un système Hamiltonien. En eﬀet, d’après le théorème KAM (cf. page §), nous avons vu que la majorité des trajectoires du système linéaire subsistent sous l’eﬀet d’une faible perturbation¹⁵. Par contre, tous les tores résonants sont détruits. La condition de résonance pour deux degrés de liberté est :

p < m, ν >= 0 ⇐ ⇒ qν − p = 0 ⇐ ⇒ ν = --, m = (p,q) ∈ ℤ × ℤ ⋆ q

(2.79)

Ces orbites sont $q$ -périodiques. Quel est le destin des tores résonants sous l’action d’une perturbation ? Le théorème des points ﬁxes de Poincaré-Birkhoﬀ (voir par exemple, Arnold et Avez, 1968) stipule que le cercle d’orbites périodiques de fréquence rationnelle $p ν = q$ pour le système non perturbé dégénère en nombre paire ( $2kq$ ) points ﬁxes : $kq$ points ﬁxes elliptiques alternés avec $kq$ point ﬁxes hyperboliques. Dans chacune des îles, ce schéma peut être appliqué de nouveau au voisinage des orbites elliptiques (cf. Fig. 2.13).

pict

FIG. 2.15:

Ouverture dynamique ( $qmax$ ) associée à l’application d’Hénon en fonction de la fréquence linéaire $ω$ normalisée par $2π$ . L’ouverture dynamique est limitée par la résonance entière pour $0 < ω < 0.013$ . Au-delà, elle est contrainte par des résonances d’ordre plus élevé. L’ouverture dynamique est nulle pour $ω = 1$ et $1∕3$ ; elle est inﬁnie en $ω = 0.5$ .

La notion de « dangerosité » d’une résonance pour la dynamique est un problème à part entière. Par exemple, les résonances $1$ et $1 ∕3$ conduisent toujours à des instabilités quelque soit l’amplitude d’oscillation initiale $q0$ : l’ouverture dynamique est nulle ; pour $ω=1∕2 2π$ , la dynamique est toujours stable : l’ouverture dynamique est inﬁnie. Pour toutes les autres résonances, il existe une amplitude maximale non nulle ( $q max$ ) au-deçà de laquelle le mouvement est stable. Au voisinage des points hyperboliques le mouvement est chaotique mais borné. Par contre aux grandes amplitudes ( $q > qmax$ ), les orbites résonantes sont suﬃsamment excitées pour détruire toutes les trajectoires. La dernière trajectoire fermée¹⁶ déﬁnit l’ouverture dynamique qui est fortement marquée par les résonances (cf. Fig. 2.14 et Fig. 2.15).

Ces résultats me permettent de rappeler, qu’un 2n-pôles peut engendrer des résonances de tout ordre. En eﬀet, il est parfois aﬃrmé, par exemple, que les hexapôles ne peuvent générer que des résonances d’ordre un et trois. Cette idée fausse provient du fait que si l’on applique une théorie de perturbation du premier ordre, seules des résonances d’ordre inférieur à $n$ apparaissent (voir par exemple Hagedorn, 1957, Schoch, 1958 ou Guignard 1978) ; de plus, les calculs ne sont en général pas faits aux ordres supérieurs (voir par exemple, Bazzani et al., 1994, Bengtsson, 1988).

2.2.5 Accélérateurs : carte en fréquence

2.2.5.1 Application fréquence

La première application de l’Analyse en Fréquence à l’étude de la dynamique des accélérateurs a été eﬀectuée sur des modèles de l’Advanced Light Source (Dumas et Laskar, 1993, Laskar et Robin, 1996).

Nous avons vu au cours du premier chapitre que la dynamique transverse d’un accélérateur est modélisée par un système à (2+1) degrés de liberté, où le troisième degré de liberté provient de la dépendance en longitude $s$ du Hamiltonien global.

L’Analyse en Fréquence se trouve être un outil particulièrement adapté à l’étude d’un tel système. En eﬀet, en adoptant les positions et moments canoniques usuelles, $(x,y,s,px,py,ps)$ , deux des conditions initiales peuvent être ﬁxées, par exemple les moments $p0x = p0y = 0$ , pour un jeu de conditions initiales $(x ,y ) 0 0$ , la trajectoire de phase est intégrée numériquement et l’on enregistre la trajectoire phase discrète $(k) (k) (k) (k) (x (t),y (t),px (t),py (t))$ sur une section de Poincaré pour une longitude $s$ donnée, typiquement $s = 0$ à chaque tour $k$ pendant un temps d’intégration $T$ (Fig. 2.16). L’approximation quasi-périodique du signal $zu(t) = u (t) + iv(t)$ est recherchée sous la forme (cf. page §) :

∑N zu(t) = aueiνut + amkei<mk,ν>t, k=1

(2.80)

avec $u=x (resp. y) et v = px$ (resp. $py$ ), $ν = (νx,νy,1)$ le vecteur de fréquences fondamentales qui correspond aux deux nombres d’ondes et à la fréquence de premier retour sur la surface de Poincaré ; $mk = (m1k, m2k, m3k)$ un multi-indice et $< mk, ν >= m1k νx + m2kνy + m3k$ .

pict

FIG. 2.16:

Section de Poincaré en $s = 0$ pour un accélérateur de particules et trajectoire de phase $(x,p ) x$ enregistrée aux temps discrets $i$ .

L’application fréquence est alors de dimension deux :

T 2 2 F : ℝ → ℝ (2.81) (x, y) ↦→ (νx, νy)

Donc, si l’on choisit les positions $(x0, y0)$ sur une grille de conditions initiales balayant par exemple l’ensemble de l’ouverture physique de l’accélérateur, on va pouvoir construire une carte en fréquence. Toute la puissance de l’Analyse en Fréquence pour un système à 3 degrés de liberté repose sur le fait qu’une carte en fréquence dévoile la dynamique de manière exhaustive. L’étude des cartes en fréquence que je présenterai dans le chapitre suivant révélera de manière naturelle les principales résonances contraignant la dynamique transverse du faisceau.

2.2.5.2 Principe de construction

Nous allons décrire la construction de la carte en fréquence de l’Advanced Light Source, carte qui a la particularité d’être très simple à lire.

Pour chacune des conditions initiales $(x0,y0)$ , la trajectoire de la particule est intégrée numériquement sur $2×$ 1 000 tours consécutifs. Ce nombre de tours est justiﬁé par le temps d’amortissement qui, pour l’ALS, est dix fois supérieur.

Si la particule a survécu (i.e. sa trajectoire n’a pas divergé), ses conditions initiales $(x0,y0)$ appartiennent à l’ouverture dynamique. Les 1 000 premiers tours de l’intégration sont utilisés pour calculer avec l’Analyse en Fréquence les nombres d’ondes $(ν(x1),ν(1y))$ . On a ainsi obtenu le premier point de l’ouverture dynamique de l’ALS et de sa carte en fréquence.

Les 1 000 tours suivants servent à calculer de nouveau les nombres d’ondes $(ν (2x),ν(y2))$ . Si l’orbite n’est pas régulière les deux jeux de fréquences vont diﬀérer, car l’orbite aura diﬀusé. On déﬁnit alors un indice de stabilité $D$ :

( ∘ (----------)----(----------)--) D = log ν(2)− ν(1) 2 + ν(2)− ν(1) 2 10 x x y y

(2.82)

La diﬀusion sera codée par des niveaux de gris ou des couleurs (Annexe B) à la fois sur la carte en fréquence et sur l’ouverture dynamique. Cet indice se révèle être un excellent critère de stabilité à long terme (voir le travail de Dumas et Laskar, 1993).

Ce processus de construction est réitéré pour chaque point de la grille de conditions initiales. Puisque la première contribution au glissement des nombres d’ondes est quadratique, $νx=ν0x+ kxxx2$ (cf. Eq. 1.76 page §), le choix le plus judicieux comme pas d’incrémentation des conditions initiales est un pas suivant une loi en racine carrée (i.e. à choisir l’action). Ce qui nous permet d’obtenir une détermination plus grossière de l’ouverture dynamique au voisinage du point de fonctionnement mais beaucoup plus ﬁne sur les bords de l’ouverture dynamique, i.e. là où précisément les résonances limitent la stabilité du faisceau.

La carte en fréquence et l’ouverture dynamique de l’ALS sont données par la ﬁgure 2.17. Le point de fonctionnement $(νx, νy ) = (14.25, 8.18 )$ est le coin supérieur droit de la carte. Le bord supérieur (respectivement inférieur) de la carte correspond à la variation des nombres d’ondes avec l’amplitude pour $νy(y)$ à $x = 0$ (respectivement $νx (x )$ à $y=0$ ).

Trois types de zones peuvent être distinguées sur la carte en fréquence (Fig. 2.17) :

les zones régulières: (zones A) où les points de l’espace des fréquence sont régulièrement espacés et caractérisés par une faible diﬀusion. Le mouvement est typiquement un mouvement bétatron régulier.
les résonances: (zones B) qui sont les droites de pente rationnelle, $pνx + qνy + r = 0$ . Si une résonance est traversée au voisinage d’une zone elliptique, on observe une droite dans l’espace des fréquences (ce phénomène est renforcé pour une résonance isolée). Si une résonance est traversée au voisinage d’une zone hyperbolique, il y a désertion de points autour de la résonance.
les zones chaotiques: (zones C) où toute structure est perdue avec une diﬀusion élevée et qui peut conduire à un comportement fortement nonlinéaire voire chaotique, par exemple aux grandes amplitudes et sur les bords de l’ouverture dynamique.

FIG. 2.17:

Carte en fréquence (a) et ouverture dynamique (b) calculées pour une maille parfaite de l’ALS pour la surface de section $s = 0$ ( $β = 11.3 m x$ et $β = 4.0 m y$ ). Trois types de zones sont observés : zones régulières (A), les résonances (B) et les régions chaotiques (C). Une résonance est révélée diﬀéremment suivant qu’elle est traversée au voisinage des régions hyperboliques ou elliptiques. La diﬀusion des orbites est codée suivant une échelle logarithmique en niveaux de gris. L’échantillonnage dans l’ouverture dynamique suit une loi en racine carrée de la distance à l’origine.

Nous obtenons une bijection entre l’espace des conﬁgurations et l’espace des fréquences. Les résonances peuvent être directement localisées sur l’ouverture dynamique. Dans cet espace (Fig. 2.17-b), les régions complètement instables correspondent à celles sans points ; les résonances ne sont pas des droites mais des « courbes » dont l’épaisseur est proportionnelle à la largeur de la résonance.

Pour une zone hyperbolique, la résonance est identiﬁée comme une courbe dont la diﬀusion est localement plus élevée.
Pour une région elliptique, les îles de la résonance correspondent à grande amplitude non pas à une courbe mais à un ensemble de conditions initiales (voir Fig. 2.17-b). Toutes les particules capturées dans l’île de résonance ont les mêmes fréquences. Localement, il n’y a plus de bijection entre l’espace des fréquences et l’espace des conﬁgurations. Ce phénomène est particulièrement marqué pour les résonances isolées [sur la carte en fréquence] dont l’amplitude est très grande.

[next] [prev] [prev-tail] [front] [up]

Chapitre 2Méthodes et outils

2.1 Réalisation d’un intégrateur symplectique

2.1.1 Introduction

2.1.2 Ancienne approche

2.1.3 Approche moderne

2.1.3.1 Présentation

2.1.3.2 Le ﬂot d’un système

2.1.3.3 Intégrateur symplectique ab ovo

2.1.3.4 Méthodes de construction

2.1.4 Intégrateurs utilisés

2.1.5 Intégration des éléments parfaits

2.1.5.1 Section droite

2.1.5.2 Aimant de courbure parfait

2.1.5.3 Quadripôle droit

2.1.5.4 Hexapôle parfait droit

2.1.6 Etude comparative

2.1.6.1 Introduction

2.1.6.2 Propriétés

2.2 Analyse en Fréquence

2.2.1 Introduction

2.2.2 Application fréquence

2.2.3 Propriétés de l’application fréquence numérique

2.2.4 Applications préliminaires

2.2.4.1 Le pendule

2.2.4.2 L’application d’Hénon

2.2.5 Accélérateurs : carte en fréquence

2.2.5.1 Application fréquence

2.2.5.2 Principe de construction

Chapitre 2
Méthodes et outils