La page ci-présente existe en version notebook téléchargeable grâce au bouton Bouton (choisir le format .ipynb). On rappelle qu’l faut ensuite l’enregistrer dans un répertoire adéquat sur votre ordinateur (capa_num par exemple dans votre répertoire personnel) puis lancer Jupyter Notebook depuis Anaconda pour accéder au notebook, le modifier et exécutez les cellules de code adéquates.

6. Etude d’un mouvement à force centrale#

But : obtenir des trajectoires d’un point matériel soumis à un champ de force centrale conservatif

6.1. Position du problème#

On se propose d’étudier l’expérience de Rutherford où des particules alpha (charge \(+2e\)) d’énergie cinétique \(E_c\) bombardent des noyau d’or \(Z = 79\) supposés fixes dans un référentiel galiléen. On peut supposer que chaque particule ne rencontre qu’un seul noyau d’or qui est responsable de sa déviation.

Dans la modélisation habituelle, les particules sont initialement à l’infini et l’interaction entre le noyau et la particule est coulombienne :

Nous allons :

  • Etablir une condition sur la distance initiale pour pouvoir considérer que les particules sont initialement sur une trajectoire de diffusion.

  • Déterminer et tracer la trajectoire des particules alpha pour plusieurs paramètres d’impact \(b\) à énergie cinétique fixée.

  • Déterminer l’angle de diffusion en fonction du paramètre d’impact \(b\) et vérifier la cohérence avec la loi attendue.

  • Simuler un bombardement homogène et vérifier que la loi sur la section efficace attendue est correcte.

  • Etudier la dépendance de la précédente loi à l’énergie cinétique des particules alpha.

6.2. Modélisation du problème.#

On suppose la particules alpha initialement à une distance \(r_i\) suffisamment grande pour qu’elle puisse presque être considérée comme à l’infini (on ne peut pas la placer à l’infini pour une simulation numérique !). Le paramètre d’impact \(b\) représenté sur la figure est supposé connu. Elle possède initialement une vitesse \(v_0\) orientée comme sur le schéma.

Conditions initiales

6.2.1. Considération sur les unités#

Comme dans le cas de l’étude de la vibration moléculaire, on va se placer dans un système d’unités évitant des puissances de 10 trop grandes. On va donc travailler avec comme unité de référence :

  • la charge élémentaire devient \(e = 1 \rm{q_a}\)

  • la masse atomique \(m = 1.7 \times 10^{-27} \rm{kg}\) devient \(m = 1 \rm{m_a}\)

  • la distance \(r = 1.0 \times 10^{-14} \rm{m}\) (taille typique du noyau) devient \(r = 1 \rm{r_a}\)

  • l’énergie \(E = 1.6 \times 10^{-19} \rm{J}\) (1eV) devient \(E = 1 \rm{E_a}\)

Dans ces conditions, les données numériques utiles deviennent :

  • Masse des particules alpha : \(m = 4 \rm{m_a}\)

  • Energie cinétique initiale \(E_{c0} = 5.3 \times 10^6 \rm{E_a}\) (cette valeur sera susceptible de changer ensuite).

  • Constante \({e^2 \over 4 \pi \epsilon_0} = 1.4 \times 10^5 \rm{E_a.r_a} \)

  • un pas de temps de \(h = 1 t_a\) correspondra à \(1.0 \times 10^{-20} \rm{s}\)

6.2.2. Mise en équation#

Il s’agit d’un problème classique de force centrale coulombienne, la conservation du moment cinétique et le principe fondamental de la dynamique permettent de se ramener à un problème à trois inconnues :

\[\begin{split} {\rm{d} \over \rm{d}t} \begin{pmatrix} \theta \\ r \\ \dot r \end{pmatrix}(t) = \begin{pmatrix} {bv_0 \over r(t)^2} \\ \dot r(t) \\ {b^2 v_0^2 \over r(t)^3} + {2Ze^2 \over 4 \pi \epsilon_0 m r(t)^2} \end{pmatrix} \end{split}\]

6.2.3. Schéma d’Euler#

On utlisera directement la fonction native odeint de la bibliothèque scipy.integrate.

6.2.4. Conditions initiales#

Les conditions initiales ne sont pas directement \(\theta, r, \dot r\) mais :

\[\begin{split} \begin{pmatrix} b \\ r(t=0) = r_i \\ Ec(t=0) = E_{c0} \end{pmatrix} \end{split}\]

avec \(v_0\) suivant \(-\vec e_x\).

Il faudra donc déterminer \(\theta, r, \dot r\) à partir de ces contraintes avant de commencer l’intégration.

Pour qu’on puisse considérer que la particules est à l’infini et proche de l’axe Ox, il faut deux conditions sur \(r_i\):

  • \(r_i \gg r_{min}\) avec \(r_{min}\) le rayon limite à \(t=0\) pour lequel la particule alpha est placée dans un état de diffusion.

  • \(r_i \gg b\)

On choisira donc \(r_i = \max(1000 r_{min}, 1000 b)\).

6.2.5. Temps d’observation#

En supposant \(r_i\) suffisamment grand, le temps que mettra la particule pour atteindre à nouveau la distance \(r_i\) dans un état de diffusion est de l’ordre de \({2r_i \over v_0}\). On pourra se servir de cette estimer pour choisir un pas de temps donnant approximativement N points sur la trajectoire.

Pour se donner un peu de marge on choisira un un durée totale de simulation égale \({4r_i \over v_0}\) et un nombre de points de simulation \(N = 10000\).

6.3. Bibliothèques utiles#

import numpy as np
import matplotlib.pyplot as plt
import numpy.random as rd
from scipy.integrate import odeint

6.4. Travail sur les conditions initiales#

Exercice 1 :

  1. Trouver par le calcul une condition sur \(ri\) pour que la particule alpha, possédant une énergie cinétique \(E_{c0}\) soit dans un état de diffusion. On note \(r_{min}\) la valeur minimale trouvée. Pour pouvoir négliger l’énergie potentielle et donc l’interaction à l’instant initiale, on prend \(ri = 1000 r_{min}\) comme condition initiale.

  2. En vous aidant de la figure et du calcul précédent, écrire une fonction CI(b, Ec0) qui prend comme argument les conditions initiales b et Ec0 et qui renvoient les conditions initiales \(\theta(t=0), r(t=0), \dot r(t=0)\) sous forme d’un vecteur numpy. On considèrera que \(b\) est de signe quelconque et que \(v_0\) est positif pour un vecteur vitesse suivant \(-\vec e_x\).

6.5. Détermination de la trajectoire#

L’implémentation du schéma d’Euler sera ici un peut particulière : on ne va pas arrêter l’intégration à un temps fixé mais lorsque la particule aura atteint à nouveau la distance \(r_i\). On pourra alors considérer qu’elle est à nouveau “à l’infini”.

Vous ne connaîtrez donc pas à l’avance la taille des ensembles de points tk(temps), thetak, rk et rpointk.

Exercice 2 :

  1. Ecrire la fonction F(Y, t, C) introduite dans le schéma d’Euler explicite. Attention, Y est ici un vecteur à 3 coordonnées.

  2. Ecrire la fonction deviation(b, Ec0, N) qui doit, à partir du paramètre d’impact b, de l’énergie cinétique initiale Ec0 et du nombre de points de calculs N :

    • calculer les conditions initiales

    • créer un vecteur temps contenant les instants où on veut estimer les coordonnées polaires.

    • utiliser la fonction odeint pour réaliser l’intégration numérique

    • renvoyer le veceur temps et le tableau du résultat de l’intégration

  3. Ecrire une fonction pol_to_cart(r, theta) qui prend comme arguments deux vecteurs r et theta correspondant à des coordonnées polaires et qui renvoie les vecteurs de coordonnées x et y cartésiennes associées.

  4. Pour la valeur d’énergie cinétique donnée dans l’énoncé, déterminer puis tracer la trajectoire pour des paramètres d’impacts \(b\) allant de 0.01 (inclus) à 1000 (inclus, 16 valeurs, on utilisera logspace pour répartir les valeurs sur une échelle logarithmique) et commenter l’allure des trajectoires et l’angle de déviation observé. On choisira N = 10000 points. Est-ce qu’on attend ?

Indications utiles :

  • Pour rappel, la signature est logspace(start, stop, Np) pour créer Np point entre \(10^{start}\) et \(10^{stop}\) inclus.

6.6. Paramètre d’impact et diffusion#

On va étudier la relation entre le paramètre d’impact et l’angle de diffusion (\(D\) sur la figure ci-après).

Angle de diffusion

On veut tester la relation théorique :

\[ b = {\beta \over \tan({D \over 2})} \]

avec \(\beta = {1 \over 2E_{c0}}{2Ze^2 \over 4\pi \epsilon_0}\)

Exercice 3 :

  1. Ecrire une fonction D_th(b, Ec0) qui renvoie l’angle de diffusion pour un paramètre d’impact b et une énergie cinétique initiale Ec0 (la particule étant initialement à l’infini.

  2. Ecrire une fonction D_sim(Y) qui, partant du vecteur Y renvoyé par odeint, renvoie une estimation de l’angle de diffusion.

  3. Représenter l’angle de diffusion en fonction b suivant la relation théorique et l’intégration numérique puis représenter \(\tan(D/2)\) en fonction de \(1/b\). Utiliser ces graphiques pour commenter la fiabilité de l’intégration numérique.

6.7. Statistique des angles de diffusion#

En pratique, on n’envoie pas une particule alpha mais un jet uniformément répartie dans le plan \(y0z\). Pour un point de départ \((y_0, z_0)\), l’invariance par rotation autour de l’axe Ox permet de conserver l’étude précédente en prenant \(b = \sqrt{y_0^2 + z_0^2}\).

Les détecteur de l’expérience de Rutherford, répartie sur une sphère autour de l’échantillons n’ont pas accès à \(D\) et \(b\) mais au nombre de particules qui arrivent dans la surface \(d\Omega\) (appelée angle solide pour une sphère de rayon 1) située entre \(D\) et \(D + dD\) (\(dD\) est au final la taille d’un capteur).

Section efficace

Sur la figure, l’angle de diffusion est noté \(\theta\) et non \(D\).

On peut assez facilement calculer théoriquement ce nombre de particule \(dn\) avec le flux de particules émis en remarquant qu’elle sont forcément passées par la surface circulaire comprise entre les rayons \(b\) (angle de diffusion \(D\)) et \(b+db\) (angle de diffusion \(D +dD\)) :

\[ dn = {\left(\beta \over 2\right)}^2 {1 \over \sin^4 (D/2)} Jd\Omega \]

où

  • \(J\) est le nombre de particules émis par unité de surface.

  • \(d\Omega = 2\pi \sin D dD\)

On se propose de réaliser une simulation (de Monte-Carlo) pour vérifier l’expression théorique de la section efficace différentielle :

\[ {dn \over Jd\Omega} \]

Le principe est :

  1. On tire aléatoirement (tirage uniforme) deux coordonnées \(y_0\) et \(z_0\).

  2. On calcule \(b\) et s’il est inférieur à un \(b_{max}\) fixé (on prendra \(b_{max} = 100\)), on déterminer l’angle de diffusion \(D\) associé par intégration.

  3. On détermine alors la statistique des angles de diffusion et on calcule la section efficace différentielle en fonction de \(D\).

On réalisera \(N = 10000\) tirages et on prendra 100 bins pour la statistique des angles des fréquences statistiques. De plus, on calcule \(J\) par :

\[ J = N / (\pi bmax^2) \]

Exercice 4 :

  1. Ecrire une fonction tir qui prend comme argument Ec0, b_max, N et nbins et qui renvoie les fréquences des occurences (vecteur) des angles de diffusion en respectant les règles imposées par la simulation précédente et les valeurs d’angle de diffusion associées à chaque fréquence (vecteur).

  2. Obtenir à partir de ces deux vecteurs la section efficace différentielle en fonction \(D\) et la tracé graphiquement. Tracer son expression théorique et comparer.

Indications utiles :
On utilisera la fonction histogram de la bibliothèque numpy dont la signature est :

hist, bins = np.histogram(a, bins=Nd)

et qui renvoie pour un vecteur de valeurs a le fréquences hist associées aux valeurs bins. Le nobmre de bacs de valeurs et donné par bins.