Particule chargée dans un champ magnétique dipolaire

1. Introduction

Les particules chargées émises par le Soleil (protons et électrons) se propagent en suivant les lignes du champ magnétique solaire. Ce flux de particules constitue le vent solaire. Lorsqu'elles rencontrent le champ magnétique terrestre, certaines peuvent être piégées et rester confinées dans la magnétosphère, en particulier dans les ceintures de Van Allen.

La simulation suivante permet de retrouver les caractéristiques du mouvement des particules dans les ceintures de Van Allen. On suppose que le champ magnétique terrestre est assimilable à celui d'un dipôle et que les particules chargées sont non relativistes.

2. Équations du mouvement

On se place en coordonnées sphériques, l'axe polaire z étant l'axe du dipôle magnétique, de moment M. Le champ magnétique est :

$$\overrightarrow{B}=\frac{{\mu }_{0}M}{4\pi }\frac{(2cos\theta {\overrightarrow{e}}_r+sin\theta {\overrightarrow{e}}_{\theta })}{r^3}$$

Dans le cas non relativiste, l'équation du mouvement d'une particule de masse m et de charge q dans ce champ est :

$$m\frac{d\overrightarrow{v}}{dt}=q\overrightarrow{v}\wedge \overrightarrow{B}$$

On considère le cas d'une charge positive. Soit R une distance qui servira d'unité de longueur pour les calculs. On définit la pulsation cyclotronique :

$$\Omega =\frac{q{B}_0}{m}=\frac{q{\mu }_{0}M}{4\pi m{R}^3}$$

En prenant l'inverse de cette pulsation comme unité de temps, l'équation du mouvement devient :

$$\frac{d\overrightarrow{v}}{dt}=\overrightarrow{v}\wedge (2cos\theta {\overrightarrow{e}}_r+sin\theta {\overrightarrow{e}}_{\theta })$$

Les équations différentielles correspondantes en coordonnées sphériques sont :

$$\begin{array}{cc}\hfill & \frac{d^{2}r}{d{t}^2}-r(\frac{d\phi }{dt}{)}^2{sin}^2\theta -r(\frac{d\theta }{dt}{)}^2=-\frac{{sin}^2\theta }{r^2}\frac{d\phi }{dt}\hfill \\ \hfill & r\frac{d^2\theta }{d{t}^2}-r(\frac{d\phi }{dt}{)}^{2}sin\theta cos\theta +2\frac{dr}{dt}\frac{d\theta }{dt}=\frac{2}{r^2}cos\theta sin\theta \frac{d\phi }{dt}\hfill \\ \hfill & r\frac{d^2\phi }{d{t}^2}sin\theta +2r\frac{d\phi }{dt}\frac{d\theta }{dt}cos\theta +2\frac{dr}{dt}\frac{d\phi }{dt}sin\theta =\frac{1}{r^3}(\frac{dr}{dt}sin\theta -2r\frac{d\theta }{dt}cos\theta )\hfill \end{array}$$

3. Résolution numérique

On définit tout d'abord un système différentiel du premier ordre :

systeme={r'[t]==dr[t],
dr'[t]==-Sin[theta[t]]^2/r[t]^2*dphi[t]+r[t]*dphi[t]^2*Sin[theta[t]]^2+r[t]*dtheta[t]^2,
theta'[t]==dtheta[t],
dtheta'[t]==(2*Cos[theta[t]]*Sin[theta[t]]*dphi[t]/r[t]^2+r[t]*dphi[t]^2*Sin[theta[t]]
            *Cos[theta[t]]-2*dr[t]*dtheta[t])/r[t],
phi'[t]==dphi[t],
dphi'[t]==((dr[t]*Sin[theta[t]]-r[t]*dtheta[t]*2*Cos[theta[t]])/r[t]^3-2*r[t]*dphi[t]
            *dtheta[t]*Cos[theta[t]]-2*dr[t]*dphi[t]*Sin[theta[t]])/(r[t]*Sin[theta[t]])}
            

Comme premier exemple, considérons une particule située initialement dans le plan équatorial, avec une dérivée $\frac{d\phi }{dt}$ initialement non nulle :

initial={r[0]==1,dr[0]==0,theta[0]==Pi/2,dtheta[0]==0,phi[0]==0,dphi[0]==0.1}
tmax=200
solution=NDSolve[Join[systeme,initial],{r,dr,theta,dtheta,phi,dphi},{t,0,tmax},
                Method->{'ExplicitRungeKutta',DifferenceOrder->5},
                MaxSteps->Infinity,AccuracyGoal->6]
            

On trace r, theta et phi en fonction du temps :

Plot[r[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,2}},AxesLabel->{'t','r'}]

plotA.pdf

Plot[theta[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{'t','theta'}]

plotB.pdf

Plot[phi[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,2*Pi}},AxesLabel->{'t','phi'}]

plotC.pdf

On trace aussi la vitesse en fonction du temps pour vérifier qu'elle est constante :

Plot[Sqrt[dr[t]^2+(r[t]*dtheta[t])^2+(r[t]*Sin[theta[t]]*dphi[t])^2]/.solution,{t,0,tmax},
                PlotRange->{{0,tmax},{0,0.2}},AxesLabel->{'t','V'}]

plotD.pdf

Trajectoire dans l'espace :

ParametricPlot3D[{r[t]*Sin[theta[t]]*Cos[phi[t]],r[t]*Sin[theta[t]]*Sin[phi[t]],
                    r[t]*Cos[theta[t]]}/.solution,{t,0,tmax},PlotRange->{{-2,2},{-2,2},{-2,2}},
                    AxesLabel->{'x','y','z'},PlotPoints->1000]

plotE.pdf

On constate que la particule décrit un mouvement de rotation cyclotronique superposé à une dérive lente en longitude.

Pour obtenir un rayon cyclotronique 5 fois plus petit, on divise la vitesse initiale par 5. Ajoutons aussi une vitesse initiale perpendiculaire au plan équatorial :

initial={r[0]==1,dr[0]==0,theta[0]==Pi/2,dtheta[0]==0.02,phi[0]==0,dphi[0]==0.02}
tmax=1000
solution=NDSolve[Join[systeme,initial],{r,dr,theta,dtheta,phi,dphi},{t,0,tmax},
                Method->{'ExplicitRungeKutta',DifferenceOrder->5},
                MaxSteps->Infinity,AccuracyGoal->6]
            

Plot[r[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,2}},AxesLabel->{'t','r'}]

plotA1.pdf

Plot[theta[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{'t','theta'}]

plotB1.pdf

Plot[phi[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{'t','phi'}]

plotC1.pdf

Plot[Sqrt[dr[t]^2+(r[t]*dtheta[t])^2+(r[t]*Sin[theta[t]]*dphi[t])^2]/.solution,{t,0,tmax},
                    PlotRange->{{0,tmax},{0,0.1}},AxesLabel->{'t','V'}]

plotD1.pdf

ParametricPlot3D[{r[t]*Sin[theta[t]]*Cos[phi[t]],r[t]*Sin[theta[t]]*Sin[phi[t]],
                    r[t]*Cos[theta[t]]}/.solution,{t,0,tmax},PlotRange->{{-2,2},{-2,2},{-2,2}},
                    AxesLabel->{'x','y','z'},PlotPoints->1000]

plotE1.pdf

En superposition au mouvement cyclotronique (visible sur l'angle phi), la particule se dirige tout d'abord vers le pôle sud puis rebrousse chemin à phi=2, traverse à nouveau l'équateur et se dirige vers le pôle nord où elle rebrousse à phi=1,1. Le mouvement est donc constitué d'aller et retours entre les deux hémisphères, avec aussi une dérive lente en longitude (dérive de phi).

Augmentons la vitesse dphi[0] :

initial={r[0]==1,dr[0]==0,theta[0]==Pi/2,dtheta[0]==0.05,phi[0]==0,dphi[0]==0.02}
tmax=1000
solution=NDSolve[Join[systeme,initial],{r,dr,theta,dtheta,phi,dphi},{t,0,tmax},
                Method->{'ExplicitRungeKutta',DifferenceOrder->5},
                MaxSteps->Infinity,AccuracyGoal->6]
            

Plot[r[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,2}},AxesLabel->{'t','r'}]

plotA2.pdf

Plot[theta[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{'t','theta'}]

plotB2.pdf

Plot[phi[t]/.solution,{t,0,tmax},PlotRange->{{0,tmax},{0,Pi}},AxesLabel->{'t','phi'}]

plotC2.pdf

Plot[Sqrt[dr[t]^2+(r[t]*dtheta[t])^2+(r[t]*Sin[theta[t]]*dphi[t])^2]/.solution,{t,0,tmax},
                    PlotRange->{{0,tmax},{0,0.1}},AxesLabel->{'t','V'}]

plotD2.pdf

ParametricPlot3D[{r[t]*Sin[theta[t]]*Cos[phi[t]],r[t]*Sin[theta[t]]*Sin[phi[t]],
                        r[t]*Cos[theta[t]]}/.solution,{t,0,tmax},PlotRange->{{-2,2},{-2,2},{-2,2}},
                        AxesLabel->{'x','y','z'},PlotPoints->1000]

plotE2.pdf

Les points de retour sont à des latitudes plus élevées.