Table des matières Python

Particule chargée dans un champ magnétique axial

1. Introduction

On considère un champ magnétique à symétrie axiale. L'axe de symétrie est Oz et on utilise les coordonnées cylindriques (r,θ,z). Le potentiel vecteur est :

A=Aθ(r,z)eθ(1)

Le champ magnétique est :

Br=-Aθz(2)Bz=1r(rAθ)r(3)

On suppose que Aθ et ses dérivées sont connus en tout point.

L'équation du mouvement d'une particule de charge q et de masse m est :

dvdt=-qBmv(4)

Une intégrale première est l'énergie mécanique, qui se réduit à l'énergie cinétique :

E=12mv2(5)

Dans le cas particulier d'un champ uniforme (Bz uniforme et Br=0), la particule décrit un mouvement cyclotronique de pulsation

Ω=-qBzm(6)

La projection de la trajectoire dans le plan équatorial (perpendiculaire à l'axe Oz) est circulaire, décrite à la pulsation Ω. On s'intéresse aussi à la trajectoire dans le plan méridien tournant, qui est définie par la courbe paramétrée (z(t),r(t)).

Ce document établit les équations différentielles qui permettent d'obtenir ces deux trajectoires dans le cas général, en vue de leur intégration numérique. Dans un premier temps, on établit les équations générales à l'aide du formalisme hamiltonien, puis on considère l'approximation paraxiale.

2. Équations de Hamilton

Les équations différentielles en coordonnées cylindriques peuvent être obtenues directement à partir de l'équation de Newton (4). Nous allons toutefois utiliser le formalisme Hamiltonien, qui a l'avantage de faire apparaitre simplement la conséquence de la symétrie axiale, c'est-à-dire l'invariance par rotation autour de l'axe Oz. De plus, le formalisme hamiltonien est nécessaire pour poser un problème quantique.

On commence par définir le lagrangien (fonction des coordonnées et de leur dérivée) :

L=12mv2+qvA(7)

Pour montrer que ce lagrangien conduit bien aux équations de Newton (3), on se place en coordonnées cartésiennes :

L(x,y,z,x˙,y˙,z˙)=12m(x˙2+y˙2+z˙2)+q(x˙Ax(x,y,z)+y˙Ay(x,y,z)+z˙ Az(x,y,z))(8)

L'équation de Lagrange pour la première coordonnée est :

Lx-ddtLx˙=0(9)

En développant cette équation avec l'expression ci-dessus du lagrangien, on obtient bien la composante sur x de l'équation de Newton. Il en est de même pour les deux autres composantes.

Le moment conjugué de la variable x est par définition :

px=Lx˙=mx˙+qAx(10)

On a donc :

p=mv+qA(11)

L'hamiltonien est construit à partir du lagrangien de manière à éliminer les dérivées les dérivées des coordonnées au profit des moments conjugués. Il est défini par :

H=x˙px+y˙py+z˙pz-L(12)

On obtient ainsi :

H(x,y,z,px,py,pz)=12m[(px-qAx(x,y,z))2+(py-qAy(x,y,z))2+(pz-qAz(x,y,z))2](13)

Les équations de Hamilton (équivalentes aux équations de Lagrange ou de Newton), sont alors (pour la coordonnée x seulement) :

dxdt=Hpx(14)dpxdt=-Hx(15)

Nous obtenons deux équations du premier ordre au lieu d'une seule équation du second ordre. Considérons à présent le cas des coordonnées cylindriques pour un champ magnétique à symétrie axiale. Le lagrangien est :

L(r,θ,z,r˙,θ˙,z˙)=12m(r˙2+r2θ˙2+z˙2)+qrθ˙Aθ(r,z)(16)

Les moments conjugués sont :

pr=Lr˙=mr˙(17)pθ=Lθ˙=mr2θ˙+qrAθ(18)pz=Lz˙=mz˙(19)

L'hamiltonien est :

H(r,z,pr,pθ,pz)=12m(r˙2+r2θ˙2+z˙2)(20)=12m[(pθr-qAθ(r,z))2+pr2+pz2](21)

La symétrie axiale conduit à un hamiltonien indépendant de θ. On a donc :

dpθdt=-Hθ=0(22)

Le moment angulaire pθ est donc une constante du mouvement. Ce résultat peut être retrouvé avec l'équation de Newton en projection sur θ. La connaissance de r et z permet donc de calculer la vitesse orthoradiale :

rθ˙=1m(pθr-qAθ)(23)

On obtient finalement les équations de Hamilton pour les variables r et z :

drdt=prm(24)dzdt=pzm(25)dprdt=1m(pθr-qAθ)(pθr2+qAθr)(26)dpzdt=1m(pθr-qAθ)qAθz(27)

On obtient ainsi un système de 4 équations du premier ordre adapté à l'intégration numérique, qui permet d'obtenir la trajectoire dans le plan méridien tournant. Pour obtenir la rotation du plan méridien, on ajoute l'équation différentielle suivante :

dθdt=1m1r(pθr-qAθ)(28)

La singularité du terme Aθ/r en r=0 n'est qu'apparente. On peut en effet écrire la même équation sous la forme :

dθdt=pθmr2-qm(Bz-Aθr)(29)

Pour le calcul numérique, on pose m=1 et on remplace q par le rapport α=q/m.

Si la particule se trouve initialement dans une zone de champ nul, ou très faible, avec une vitesse initiale contenue dans le plan méridien, alors pθ=0, ce qui simplifie les équations. Si la particule se trouve initialement sur l'axe, on a aussi pθ=0.

Si la position initiale (r0,z0) est hors de l'axe, avec une vitesse initiale dans le plan méridien, on a :

pθ=qr0Aθ(r0,z0)(30)

Dans ce cas, pθ est non nul est les équations différentielles présentent une singularité lorsque r=0. Lorsque r tend vers 0, la vitesse orthoradiale tend vers l'infini. Or l'énergie cinétique étant une constante du mouvement, cette situation ne peut se produire. Cela signifie que si pθ est non nul, la particule ne coupe jamais l'axe, et la singularité n'est pas rencontrée.

Si le champ magnétique est connu sans le potentiel vecteur, on ne peut pas utiliser le moment angulaire. Dans ce cas, on doit résoudre les équations suivantes, qui se déduisent directement de l'équation de Newton :

d2rdt2=r(dθdt)2+qmrdθdtBz(31)d2θdt2=1r(-2drdtdθdt+qm(dzdtBr-drdtBz))(32)d2zdt2=-qmrdθdtBr(33)

En raison du terme 1/r dans la deuxième, ces équations ne sont pas utilisables pour une trajectoire partant de l'axe.

3. Approximation paraxiale

3.a. Équations

L'approximation paraxiale est utilisée en optique électronique pour traiter le mouvement d'une particule qui reste voisine de l'axe ([1]). Au voisinage de l'axe, un développement limité à l'ordre 2 du potentiel vecteur donne :

Aθ(r,z)=r2Bz(0,z)(34)

Le champ magnétique à l'ordre 1 est donc :

Br(r,z)=-r2Bz(0,z)z(35)Bz(r,z)=Bz(0,z)(36)

Le principal intérêt de cette approximation est de permettre les calculs avec la seule connaissance du champ sur l'axe.

Les équations différentielles deviennent :

d2rdt2=pθ2m2r3-q24m2rBz2(37)d2zdt2=qmpθdBzdz(38)dθdt=pθmr2-qm12Bz(39)

Dans la deuxième équation, le terme d'ordre 2 a été négligé.

Dans le cas où le moment angulaire est nul pθ=0, la vitesse axiale vz est constante. On obtient alors l'équation différentielle pour la trajectoire r(z) :

d2rdz2=-q24m2vzrBz2(40)

Cette équation est linéaire : si r(z) est une trajectoire, λr(z) est aussi une trajectoire.

Dans le cas où le moment angulaire est non nul, les équations deviennent non linéaires.

3.b. Exemple avec moment angulaire nul

Le champ d'une spire circulaire de rayon a sur son axe est :

Bz(0,z)=B0(1+z2a2)32(41)

On s'intéresse à une trajectoire partant d'un point I de l'axe en z=z0, ce qui permet d'avoir pθ=0. On résout numériquement le système suivant :

drdt=pr(42)dprdt=-α24rBz(0,z)2(43)dθdt=-α2Bz(0,z)(44)

Avec de plus :

z=z0+vzt(45)

On fixe une valeur initiale de la vitesse radiale pr0. Pour des valeurs données de α=q/m et de vz, soit r0(z) la trajectoire obtenue pour une valeur pr0 de la vitesse radiale initiale. En raison de la linéarité des équations, une autre valeur initiale qr1 de la vitesse radiale conduit à la trajectoire suivante :

r1(z)=pr1pr0r0(z)(46)

Si ces trajectoires coupent l'axe Oz, elles le coupent au même point. Cela démontre le stigmatisme pour les points de l'axe.

import math
import numpy
from matplotlib.pyplot import *
import scipy.integrate

alpha=1.0
a1=alpha*alpha/4
a2=alpha/2
vz=1.0
z0=-5.0
def equation(y,t):
    z = z0+vz*t
    B = math.pow(1+z*z,-1.5)
    return [y[1],-a1*y[0]*B*B,-a2*B]
tmax=20.0
te=0.01                         
t = np.arange(start=0,stop=tmax,step=te)
atol=1e-6
rtol=1e-6
y0=[0,0.01,0]
y = scipy.integrate.odeint(equation,y0,t,rtol=rtol,atol=atol)
yt = y.transpose()
r=yt[0]
theta=yt[2]
z=z0+vz*t
figure(figsize=(10,5))
plot(z,r)
xlabel('z')
ylabel('r')
grid()
            
figAfigA.pdf

Les particules issues de I avec la même vitesse axiale vz convergent vers un point conjugué J à environ z=14 rayons de spire. La spire constitue donc une lentille convergente. Le stigmatisme est obtenu dans le cadre de l'approximation paraxiale, à condition que les particules aient la même vitesse axiale. Voyons l'évolution de l'angle du plan méridien :

figure(figsize=(10,5))
plot(z,theta)
xlabel('z')
ylabel('theta')
grid()
            
figBfigB.pdf

Le plan méridien a tourné de 1 radians.

Voyons ce qui se passe pour une vitesse axiale plus faible :

vz=0.8
y0=[0,0.01,0]
y = scipy.integrate.odeint(equation,y0,t,rtol=rtol,atol=atol)
yt = y.transpose()
r=yt[0]
theta=yt[2]
z=z0+vz*t
figure(figsize=(10,5))
plot(z,r)
xlabel('z')
ylabel('r')
grid()
            
figCfigC.pdf

Le point conjugué J est beaucoup plus proche, à z=5.

Si la vitesse axiale est plus grande :

vz=2.0
y0=[0,0.01,0]
y = scipy.integrate.odeint(equation,y0,t,rtol=rtol,atol=atol)
yt = y.transpose()
r=yt[0]
theta=yt[2]
z=z0+vz*t
figure(figsize=(10,5))
plot(z,r)
xlabel('z')
ylabel('r')
grid()
            
figDfigD.pdf

Dans ce cas, la trajectoire s'éloigne toujours de l'axe. On a donc l'équivalent d'une lentille divergente, avec formation d'une image J virtuelle.

Si on réduit la vitesse axiale jusqu'à une valeur très faible :

vz=0.001
tmax=10000.0
te=0.01
t = np.arange(start=0,stop=tmax,step=te)
y0=[0,0.001,0]
y = scipy.integrate.odeint(equation,y0,t,rtol=rtol,atol=atol)
yt = y.transpose()
r=yt[0]
theta=yt[2]
z=z0+vz*t
figure(figsize=(10,5))
plot(z,r)
xlabel('z')
ylabel('r')
grid()
            
figEfigE.pdf
figure(figsize=(10,5))
plot(z,theta)
xlabel('z')
ylabel('theta')
grid()
            
figFfigF.pdf

On obtient dans ce cas une trajectoire enroulée autour des lignes de champ, avec un très grand nombre de rotations du plan méridien. Au voisinage du centre de la spire, la trajectoire peut être assimilée à un cercle qui dérive lentement dans la direction axiale.

3.c. Exemple avec moment angulaire non nul

On reprend le champ magnétique de la spire. La particule se trouve initialement hors de l'axe au point (r0,z0) avec une vitesse orthoradiale nulle. Le moment angulaire est :

pθ=qr0Aθ(r0,z0)(47)

Il est non nul; il faut donc utiliser les équations complètes, qui sont non linéaires. Il faut remarquer que la non linéarité est d'autant plus grande que r est petit. Il faut donc s'attendre à des solutions complètement différentes de celles obtenues pour un moment angulaire nul.

On définit le système différentiel, qui comporte les constantes α=q/m et qθ :

def equation(y,t):
    r=y[0]
    r2 = r*r
    z=y[1]
    u = 1+z*z
    B = math.pow(u,-1.5)
    dB = -3*z*math.pow(u,-2.5)
    if ptheta==0:
        dpr = -alpha/4*r*B
        dpz = 0.0
        dtheta = -alpha/2*B
    else:
        dpr = ptheta/(r2*r)-alpha/4*r*B
        dpz = alpha*ptheta*dB
        dtheta = ptheta/r2-alpha/2*B
    return [y[2],y[3],dpr,dpz,dtheta]
                

La fonction suivante calcule et trace une trajectoire dans le plan méridien tournant :

def trajectoire(z0,r0,pz0,pr0,tmax):
    global ptheta
    B0 = math.pow(1+z0*z0,-1.5)
    ptheta = r0*r0/2*B0
    te = 0.01
    t = np.arange(start=0,stop=tmax,step=te)
    y = scipy.integrate.odeint(equation,[r0,z0,pr0,pz0,0.0],t,rtol=1e-5,atol=1e-5)
    yt = y.transpose()
    r = yt[0]
    z = yt[1]
    plot(z,r)
                

La fonction suivante trace l'angle du plan méridien en fonction de z :

def angle(z0,r0,pz0,pr0,tmax):
    global ptheta
    B0 = math.pow(1+z0*z0,-1.5)
    ptheta = r0*r0/2*B0
    te = 0.01
    t = np.arange(start=0,stop=tmax,step=te)
    y = scipy.integrate.odeint(equation,[r0,z0,pr0,pz0,0.0],t,rtol=1e-5,atol=1e-5)
    yt = y.transpose()
    z = yt[1]
    theta = yt[4]
    plot(z,theta)
                

La fonction suivante trace la projection de la trajectoire dans un plan perpendiculaire à l'axe (plan équatorial) :

def projection(z0,r0,pz0,pr0,tmax):
    global ptheta
    B0 = math.pow(1+z0*z0,-1.5)
    ptheta = r0*r0/2*B0
    te = 0.01
    t = np.arange(start=0,stop=tmax,step=te)
    y = scipy.integrate.odeint(equation,[r0,z0,pr0,pz0,0.0],t,rtol=1e-5,atol=1e-5)
    yt = y.transpose()
    r = yt[0]
    theta = yt[4]
    plot(r*numpy.cos(theta),r*numpy.sin(theta))
                

Voici un exemple, avec une trajectoire issue d'un point sur l'axe et une autre d'un point très proche de l'axe, avec les mêmes vitesses axiale et radiale initiales :

alpha = 1.0
figure(figsize=(8,6))
z0 = -10.0
trajectoire(z0,0.1,1.0,0.01,50.0)
trajectoire(z0,0.0,1.0,0.01,50.0)
axis([z0,10,-0.5,0.5])
xlabel('z')
ylabel('r')
grid()
                
figGfigG.pdf

Au départ, les trajectoires sont pratiquement rectilignes et parallèles car le champ magnétique est très faible. Un fois le plan de la spire franchi, on observe une trajectoire complètement différente dans le cas r0=0.1. Nous avons effectivement montré plus haut que dans le cas d'un moment angulaire non nul, r ne peut s'annuler. Voyons l'angle du plan méridien :

figure(figsize=(8,6))
angle(z0,0.1,1.0,0.01,50.0)
angle(z0,0.0,1.0,0.01,50.0)
xlabel('z')
ylabel('theta')
axis([z0,10,-1.0,0.0])
grid()
                
figHfigH.pdf

Les rotations sont très voisines. Projection des trajectoires sur le plan équatorial :

figure(figsize=(6,12))
subplot(2,1,1)
projection(z0,0.0,1.0,0.01,50.0)
xlabel('x')
ylabel('y')
grid()
axis([-0.5,0.5,-0.5,0.5])
subplot(2,1,2)
projection(z0,0.1,1.0,0.01,50.0)
xlabel('x')
ylabel('y')
axis([-0.5,0.5,-0.5,0.5])
grid()
                
figIfigI.pdf

Dans le cas r0=0, la trajectoire coupe l'axe. Dans le cas r0=0.1, la particule s'approche très près de l'axe puis s'en éloigne. Le point de rebroussement provient du terme non linéaire de l'équation différentielle radiale, qui se manifeste à très courte distance de l'axe.

Pour voir s'il y stigmatisme pour un point hors de l'axe, on considère des trajectoires avec des vitesses radiales différentes :

alpha = 1.0
figure(figsize=(8,6))
trajectoire(z0,0.1,1.0,0.005,50.0)
trajectoire(z0,0.1,1.0,0.01,50.0)
trajectoire(z0,0.1,1.0,0.02,50.0)
trajectoire(z0,0.1,1.0,0.04,50.0)
axis([z0,10,-0.5,0.5])
xlabel('z')
ylabel('r')
grid()
                
figJfigJ.pdf

Pour un point hors de l'axe, le stigmatisme n'est plus parfait mais seulement approché. C'est une conséquence de la non linéarité des équations. À l'inverse, lorsque le point est sur l'axe, les équations sont linéaires et le stigmatisme est parfait (dans la mesure ou l'approximation paraxiale est vérifiée).

Références
[1]  O. Klemperer, M.E. Barnett,  Electron Optics,  (Cambridge University Press, 1971)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.