Comprendre : Page pour l'impression

Dans cette section on répond à la question suivante: pour un système planétaire donné, quelles seront les observations ? Dans la section précédente on a introduit la fonction $f: \theta \in \mathbf{R}^n,t \in \mathbf{R}^m , \epsilon \in X^p \rightarrow y \in \mathbf{R}^m$ , où désigne un espace de variables aléatoires, qui donne des observations en fonction d'instants d'observation , d'erreurs alétoires $\epsilon$ et de paramètre du modèle $\theta$ incluant masse de la planète, de l'étoile, paramètres de l'orbites, etc. qui seront précisés. En l'occurrence, l'erreur est bien représentée par un "bruit additif", c'est à dire les observations sont de la forme $y=f(\theta,t) + \epsilon$ , où désigne une fonction modélisant la physique du phénomène observée, celle qui donnerait les observations "parfaites", sans aucune source d'erreur aléatoire. L'erreur $\epsilon$ modélise tous les phénomènes aléatoires intervenant dans les mesures. Selon la distinction adoptée ici, les paramètres du modèle peuvent inclure une source indésirable mais non aléatoire.

L'objet de ce chapitre est de construire cette fonction f_p avec une modélisation physique. La modélisation la plus précise, utilisant tout ce que l'on sait de la physique serait ici inutilement complexe, car compte tenu des erreurs de mesures la différence avec certains modèles plus simples serait si faible qu'il serait impossible de la détecter avec un niveau de confiance acceptable. Le niveau de précision du modèle présenté ici est couramment utilisé par les observateurs.

Remarque: En pratique, les algorithmes d'estimation de paramètres sont des algorithmes d'optimisation, qui nécessitent tous de donner un ou plusieurs points de départs pour la recherche. Lorsque des modèles plus complexes sont utilisés, les estimations précises peuvent être utilisées comme tels.

On présentera d'abord une modélisation physique de l'objet étudié, puis d'autres effets physiques à prendre en compte. Les bruits instrumentaux seront traités dans le chapitre suivant. Les effets physiques "indépendants de l'observateur" qui seront pris en compte sont:

Le mouvement de l'étoile dû aux planètes.
Le mouvement propre: mouvement de translation du barycentre du système {Etoile+Planètes}
Des effets stellaires variables (Ondes de pression, granularité, activité)

Ensuite, on modélise "comment" la lumière émise par l'étoile nous parvient, ce qui amènera à présenter:

La prise en compte du mouvement de la Terre (vitesse de celle-ci dans le référentiel barycentrique du système solaire pour les vitesses radiales, parallaxe et réduction des observations en astrométrie.
Les perturbations atmosphériques
Les variations de parallaxe et de mouvement propre.

Problème à deux corps

Problème à deux corps newtonien

Equation de Newton

Les planètes et l'étoile ont un certain volume. Cependant, on peut montrer que lorsque la distance entre deux corps en intéraction gravitationnelle augmente, leur comportement se rapproche de plus en plus de celui de deux points matériels. On néglige aussi les effets relativistes, de sorte que l'on est ramené à la modélisation de l'intéraction de deux points matériels, dont la résolution va suivre. Considérons deux points matériels (planète) et (pour "star" de sorte à éviter un conflit de notation avec l'anomalie excentrique ) de masses respectives et et l'origine d'un repère Galiléen. Le principe fondamental de la dynamique s'écrit: $m \frac{d^2\overrightarrow{OP}}{dt^2} = -GmM \frac{\overrightarrow{S P}}{S P^3}$ où est la constante universelle de gravitation.. En posant $\overrightarrow{r} = \overrightarrow{SP}$ et $\mu = G(M+m)$ on obtient:

$\frac{d^2 \overrightarrow{r}}{dt^2} = -G(m+M) \frac{ \overrightarrow{r}}{r^3} = -\mu \frac{\overrightarrow{r}}{r^3}$

On va montrer que la solution de cette équation décrit une conique plane, une ellipse dans le cas des planètes liées à une étoile. On suppose que le système est isolé, on peut donc choisir comme origine du repère le barycentre du système, ce qui permet d'obtenir facilement $\overrightarrow{OP}$ et $\overrightarrow{OS}$ par la relation $m\overrightarrow{OP} + M\overrightarrow{OS} = \overrightarrow{0}$ .

On suppose que le système est isolé, donc le mouvement du barycentre du système {Etoile+planètes} est rectiligne uniforme dans un référentiel galiléen (par exemple le référentiel barycentrique du système solaire).

Quantités conservées

L'équation de Newton est une équation différentielle de degré deux sur des vecteurs de $\mathbf{R}^3$ . Pour la résoudre il faut trouver six quantités conservées au cours du mouvement (ou intégrales premières du mouvement) indépendantes. En l'occurrence on peut facilement montrer que deux vecteurs sont conservés au cours du mouvement (ce qui fait deux fois trois composantes, on a bien six scalaires conservés). Notons $\overrightarrow{r}} = r \overrightarrow{u}$ , où $\overrightarrow{u}$ est un vecteur unitaire (et donc est la norme de $\overrightarrow{r}$ ) et $\dot{\overrightarrow{a}}$ la dérivée par rapport au temps d'un vecteur $\overrightarrow{a}$ . On a la conservation du moment cinétique par unité de masse $\overrightarrow{G} = \overrightarrow{r} \wedge \dot{\overrightarrow{r}}$ et du vecteur excentricité $\overrightarrow{P} = \frac{ \dot{\overrightarrow{r}} \wedge \overrightarrow{G} }{\mu} - \overrightarrow{u}$ au cours du mouvement c'est à dire leur dérivée temporelle est nulle

En effet, soit $\overrightarrow{r}(t)$ une solution de l'équation de Newton. Alors $\frac{d \overrightarrow{G}}{dt}(t) = \frac{d( \overrightarrow{r}(t) \wedge \dot{\overrightarrow{r}(t)})}{dt} = \dot{\overightarrow{r}}(t) \wedge \dot{\overightarrow{r}}(t) + \overightarrow{r}(t) \wedge \ddot{\overightarrow{r}(t)} =0$ car d'après l'équation de Newton, $\ddot{\overrightarrow{r}$ est colinéaire à $\overrightarrow{r}$ . D'autre part, $\frac{d\overrightarrow{P}}{dt}(t) = \frac{\ddot{\overrightarrow{r}}(t) \wedge \overrightarrow{G}(t) }{\mu} - \dot{\overrightarrow{u}}(t)$ , or

$\ddot{\overrightarrow{r}} \wedge \overrightarrow{G} = - \frac{\mu}{r^3} \overrightarrow{r} \wedge (\overrightarrow{r} \wedge \dot{\overrightarrow{r}}) = - \frac{\mu}{r^3} \overrightarrow{r} \wedge (\overrightarrow{r} \wedge r \dot{\overrightarrow{u}}) = - \frac{\mu}{r^3} r \overrightarrow{u} \wedge (r \overrightarrow{u} \wedge r \dot{\overrightarrow{u}}) = \mu \dot{\overrightarrow{u}}$

Donc $\frac{d\overrightarrow{P}}{dt} = \overrightarrow{0}$ .

Géométrie de l'orbite

Les quantités conservées définies page précédente permettent de donner une description géométrique de l'évolution de $\overrightarrow{r}$ .

On déduit de la conservation du moment cinétique que le mouvement est plan. En effet, $\overrightarrow{r}\cdot\overrightarrow{G} = \overrightarrow{r} \cdot (\overrightarrow{r} \wedge \dot{ \overrightarrow{r}) }$ et comme $\overrightarrow{r} \cdot (\overrightarrow{r} \wedge \dot{ \overrightarrow{r}) }$ est un produit mixte, il est invariant par permutation circulaire: $\overrightarrow{r} \cdot (\overrightarrow{r} \wedge \dot{ \overrightarrow{r}) } = \dot{\overrightarrow{r}} \cdot (\overrightarrow{r} \wedge \overrightarrow{r}) } = \overrightarrow{0}$ car le produit vectoriel de deux vecteurs colinéaires est nul. Le vecteur $\overrightarrow{r}$ est orthogonal à $\overrightarrow{G}$ à tout instant, autrement dit le mouvement est dans un plan orthogonal à $\overrightarrow{G}$ .

Notons l'angle entre $\overrightarrow{r}$ et $\overrightarrow{P}$ et posons $e = \|\overrightarrow{P} \|$ . Alors comme $\overrightarrow{u}$ est unitaire, $\overrightarrow{P} \cdot \overrightarrow{u} = e \cos v$ d'autre part en remplaçant $\overrightarrow{P}$ par sa définition, on a $\overrightarrow{P} \cdot \overrightarrow{u} = \frac{\overrightrarrow{r} \cdot ( \dot{\overrightarrow{r}} \wedge \overrightarrow{G} )}{r \mu} -1$ et comme $\overrightarrow{r} \cdot ( \dot{\overrightarrow{r}} \wedge \overrightarrow{G} )$ est un produit mixte, $\overrightarrow{r} \cdot ( \dot{\overrightarrow{r}} \wedge \overrightarrow{G} ) = \overrightarrow{G} \cdot ( \overrightarrow{r} \wedge \dot{\overrightarrow{r}} ) = \overrightarrow{G} \cdot \overrightarrow{G} = G^2$ . où est la norme de $\overrightarrow{G}$ . On obtient alors en fonction de :

$r(v) = \frac{p}{1+e \cos v}$

Où $p = \frac{G^2}{\mu}$ , qui est l'équation polaire d'une conique du plan. Cette équation donne une paramétrisation de la solution en fonction de , appelée anomalie vraie. Cependant, nous voulons exprimer la solution en fonction du temps, l'objet de la page suivante est d'exhiber une relation entre et le temps. On sait que lorsque e < 1 , il s'agit de l'équation d'une ellipse. On peut montrer géométriquement que p = a(1-e^2) où est le demi-grand axe de l'ellipse. Si e = 1 ou e>1 , la trajectoire est respectivement parabolique ou hyperbolique. Dans ces deux cas le mouvement n'est pas borné, il concernerait une planète en phase d'éjection, événement dont l'observation est très improbable et indiscernable d'une planète à très longue période ou du mouvement propre sur les données actuelles.

L'orbite est dans un plan perpendiculaire à $\overrightarrow{G}$ , et le vecteur $\overrightarrow{P}$ est parallèle à $\overrightarrow{FA}$ . Notons $\overrightarrow{I}= \frac{\overrightarrow{P}}{\|\overrightarrow{P} \|}$ et $\overrightarrow{K}= \frac{\overrightarrow{G}}{\|\overrightarrow{G} \|}$ . On introduit un vecteur $\overrightarrow{J}$ de sorte que $(\overrightarrow{I}, \overrightarrow{J}, \overrightarrow{K})$ forme une base orthonormale.

Géométrie du mouvement elliptique

Le mouvement est dans le plan $(\overrightarrow{P}, \overrightarrow{Y})$

Loi des aires

Nous avons montré que la solution au problème des deux corps est plane, et peut s'exprimer en fonction de l'anomalie vraie , angle entre le vecteur excentricité et le vecteur $\overrightarrow{r}$ .

$r(v) = \frac{a(1-e^2)}{1+e \cos v}$

Afin d'exprimer en fonction du temps on va introduire successivement deux variables, et appelées respectivement 'anomalie moyenne et 'anomalie excentrique.

Loi des Aires

Avec les notations précédentes, $\overrightarrow{u}(v) = \cos v \overrightarrow{I} + \sin v \overrightarrow{J}$ , donc $\dot{\overrightarrow{u}}(v) = -\sin v \overrightarrow{I} + \cos v \overrightarrow{J}$ en remplaçant dans, $\overrightarrow{G} = r \overrightarrow{u} \wedge (\dot{r}\overrightarrow{u}} + r \dot{\overrightarrow{u}})$ on obtient:

$r^2 \frac{dv}{dt}=G$

Cette équation s'appelle la Loi des aires et signifie que $\overrightarrow{r}$ balaie les aires à vitesse constant . En particulier l'aire totale de l'ellipse est parcourue en un certain temps fixe: le mouvement est périodique. L'aire de l'ellipse vaut:

$\pi a b = \pi a^2 \sqrt{1-e^2} = \int_{0}^{2\pi} \frac{1}{2} r^2 dv = \int_{0}^T \frac{G}{2}dt = \frac{T}{2} \sqrt{\mu a (1-e^2)}$

On retrouve bien la troisième loi de Kepler : $n^2 a^3 = \mu$ avec $n = \frac{2 \pi}{T}$ . On définit alors l'anomalie moyenne par M = n(t-t_0) où t_0 est le temps de passage au périastre ( sur la figure), qui est proportionnelle à l'aire FAC (voir figure).

Exercices

Conservation de l'énergie

Difficulté : ☆☆

Question 1)

Pour établir l'équation du mouvement, les conservations du moment cinétique et du vecteur eccentricité suffisent. On a une autre intégrale du mouvement: l'énergie. Considérons deux points matériels et de masses respectives et . On pose $\overrightarrow{r} = \overrightarrow{SP}$ , $r = \|\overrightarrow{r}\|$ et $\mu = \mathcal{G}(m+M)$ où $\mathcal{G}$ est la constante universelle de gravitation.

Montrer que $\frac{d^2\overrightarrow{r}}{dt^2} = -\mu \frac{\overrightarrow{r}}{r^3}$
Calculer l'énergie potentielle associée à la force de gravitation
On note $v = \left\| \frac{d \overrightarrow{r}}{dt} \right\|$ . Déduire du théorème de l'énergie mécanique que $\frac{1}{2} v^2 - \frac{\mu}{r}$ est une quantité constante au cours du mouvement
Quel est le point de l'orbite où la vitesse de l'étoile est la plus élevée ? la moins élevée ?

Equation de Kepler

L'équation de Kepler est une relation entre l'anomalie excentrique et l'anomalie moyenne, cette page présente un moyen de l'établir. L'aire $\mathcal{A}(FAC)$ est proportionnelle à l'anomalie moyenne .

$\mathcal{A} (AFC) = \frac{M}{2 \pi} \pi a^2 \sqrt{1-e^2} = \frac{1}{2} a^2 \sqrt{1-e^2} M$

L'ellipse de la trajectoire est obtenue par une affinité sur l'axe $\overrightarrow{Y}$ de rapport $\frac{b}{a} = \sqrt{1-e^2}$ . Donc

$\mathcal{A}(AFC') = \frac{\mathcal{A}(AFC)}{\sqrt{1-e^2}} = \frac{a^2 M}{2}$

Par ailleurs en notant l'angle $\widehat{AOC'}$

$\mathcal{A}(AOC') = \frac{1}{2} a^2 E = \mathcal{A}(FOC') + \mathcal{A}(AFC')$

L'aire du triangle FOC' s'obtient facilement car $HC' = HC / \sqrt{1-e^2}$

$\mathcal{A}(FOC') = \frac{1}{2} a^2e \sin E$

On a finalement l'équation de Kepler

$E - e\sin E = M$

Cette équation est "transcendante", en conqéquence il n'existe pas d'expression analytique de en fonction de . Cependant, on peut développer en puissances de .

Anomalies

Représentation de l'anomalie vraie

, excentrique

et moyenne

Solutions en fonction du temps

Il reste à trouver une relation entre et , ce qui s'obtient aisément en exprimant la position (X,Y) de dans le plan $(\overrightarrow{I}, \overrightarrow{J})$ :

$\left\{ \begin{array}{l l l} X =& r \cos v =& a(\cos E - e) \\ Y =& r \sin v =& a \sqrt{1-e^2} \sin E\end{array}$

On en déduit les relations utiles: $\begin{array}{l l} \cos v = \frac{\cos E - e}{1 - e \cos E}, \quad & \sin v = \frac{\sqrt{1-e^2} \sin E}{1 - e \cos E} \end{array}$

A ce point, rappelons que l'objectif est d'exprimer la position et la vitesse de l'étoile comme des observables. En exprimant en fonction de puis en fonction de on a la position du corps à un instant quelconque. En ce qui concerne l'astrométrie, on peut s'arrêter aux équations ci-dessus et simplement faire un changement de référentiel, c'ést à dire exprimer la projection de $X(t)\overrightarrow{I} + Y(t)\overrightarrow{J}$ sur la sphère céleste et choisir un jeu de paramètres $\theta$ à ajuster. Pour les vitesses radiales nous avons besoin de $\frac{dX}{dt} = \frac{dX}{dE} \frac{dE}{dM} \frac{dM}{dt}$ . En dérivant l'équation de Kepler par rapport à on obtient:

$\frac{dE}{dM} = \frac{a}{r}$

D'où:

$\left\{ \begin{array}{l l l} \dot{X} =& -na\frac{\sin E}{1-e \cos E} &= -\frac{na}{\sqrt{1-e^2}} \sin v \\ \dot{Y} =& na\sqrt{1-e^2}\frac{\cos E}{1-e \cos E} &= \frac{na}{\sqrt{1-e^2}} (e+\cos v) \end{array}$

En pratique, il n'est pas nécessaire de calculer en fonction de . Les expressions ci-dessus suffisent. Pour mémoire, en remplaçant $\cos v$ dans la formule trigonométrique $\tan^2 \frac{v}{2} = \frac{1-\cos v}{1+\cos v}$ , on obtient:

$\tan \frac{v}{2} = \sqrt{\frac{1+e}{1-e}} \tan \frac{E}{2}$

Changement de référentiel

Les observations sont disponibles dans un référentiel $(\overrightarrow{i}, \overrightarrow{j}, \overrightarrow{k})$ où $\overrightarrow{k}$ est la direction d'observation et $(\overrightarrow{i}, \overrightarrow{j})$ sont choisis de sorte que le repère est orthonormé direct. Le repère $(\overrightarrow{I}, \overrightarrow{J}, \overrightarrow{K})$ est lui aussi orthormé direct. La matrice de passage du repère orbital au repère d'observation est donc une rotation, que l'on décompose en trois rotations dont les angles ont des noms usuels.

Une rotation d'angle $\omega$ d'axe $\overrightarrow{K}$ (parallèle au moment cinétique), cet angle est appelé argument du périastre. Le repère obtenu par cette rotation est noté $(\overrightarrow{I'}, \overrightarrow{J'}, \overrightarrow{K})$ .
Une rotation d'angle d'axe $\overrightarrow{I'}$ , c'est l' inclinaison
Une rotation d'angle $\Omega$ d'axe $\overrightarrow{k}$ , c'est l'ascension droite au noeud ascendant

$\left( \begin{array}{cc} x & \dot{x} \\ y & \dot{y} \\ z & \dot{z} \end{array} \right) = \mathcal{R}_3(\Omega) \mathcal{R}_1(i) \mathcal{R}_3(\omega) \left( \begin{array}{cc} X & \dot{X} \\ Y & \dot{Y} \\ Z & \dot{Z} \end{array} \right)$

Où:

$\mathcal{R}_1(\theta) = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos \theta & -\sin \theta \\ 0 & \sin \theta & \cos \theta \end{array} \right) \quad ; \quad \mathcal{R}_3(\theta) = \left( \begin{array}{ccc} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array} \right)$

Soit le barycentre du système {planète+étoile}, on note $\mathcal{R}_3(\Omega) \mathcal{R}_1(i) \mathcal{R}_3(\omega) = \left( \begin{array}{ccc} B & G & \star \\ A & F & \star \\ C & D & \star \end{array} \right)$ où A, B, F, G sont appelées les constantes de Thiele-Innes. On conserve la notation classique pour ces constantes, qui sautent quelques lettres de l'alphabet pour une raison inconnue des auteurs. La notation C,D n'est en revanche qu'une convention pour ce cours et ne se trouve pas spécialement dans la littérature. On ne donne pas de nom particulier aux éléments de la dernière colonne de la matrice car étant donné que Z=0 , ils n'apparaissent jamais dans les calculs

Paramètres d'un mouvement à deux corps newtonien

Le plan de référence est ici le plan d'observation.

Paramètres orbitaux classiques

On peut caractériser l'orbite par les éléments suivants:

: le demi-grand axe
: l'excentricité
: la période
$\omega$ : l'argument du périastre
: l'inclinaison
$\Omega$ : l'ascension droite au noeud ascendant

Ces éléments donnent la géométrie de l'orbite. Pour déterminer la position de la planète à un instant donné, il faut de plus connaître l'instant de son passage au périastre t_p . On peut alors calculer $M = 2\pi \frac{t-t_p}{P}$ connaisant on peut calculer l'anomalie excentrique par l'équation de Kepler, puis l'anomalie vraie . On en déduit la position sur l'ellipse par l'équation donnée page "Loi des aires", $r(v) = \frac{a(1-e^2)}{1+e \cos v}$ . Enfin, la position sur l'orbite est donnée par les rotations explicitées ci-dessus.

Conclusion pour le problème à deux corps

Dans un référentiel galiléen quelconque de centre , le mouvement de l'étoile en fonction du temps peut s'écrire:

$\begin{array}{ccc}\overrightarrow{OS}(t) &= \left( \begin{array}{c} x(t) \\ y(t) \\ z(t) \end{array}\right) &= \overrightarrow{OB}_0 + \alpha t \overrightarrow{u} + \overrightarrow{BS}(t) \\ \dot{\overrightarrow{OS}}(t) &= \left( \begin{array}{c} \dot{x}(t) \\ \dot{y}(t) \\ \dot{z} (t) \end{array}\right) &= \alpha \overrightarrow{u} + \dot{\overrightarrow{BS}}(t) \end{array}$

Où $\overrightarrow{OB}_0$ est la position du barycentre du système {Etoile, Planète} à t=0 , $\overrightarrow{u}$ est la direction du mouvement de et $\alpha$ son module. $\overrightarrow{BS}(t) = - \frac{m}{m+M} \overrightarrow{r}$ où $\overrightarrow{r}$ vérifiel'équation de Newton $\frac{d^2 \overrightarrow{r}}{dt^2} = -\mu \frac{\overrightarrow{r}}{r^3}$ , qui est une équation différentielle de degré deux sur l'espace. Lorsque $\mu$ est fixé et la position et la vitesse à t=0 , $(\overrightarrow{r}_0, \dot{\overrightarrow{r}}_0)$ sont connus, la position et la vitesse sont données par le flot: $\begin{array}{cccc } \Phi_{\mu}:& \mathbb{R} \times \mathbb{R}^3 \times \mathbb{R}^3 & \rightarrow & \mathbb{R}^3 \times \mathbb{R}^3 \\ & (t,\overrightarrow{r}_0, \dot{\overrightarrow{r}}_0)& \rightarrow &(\overrightarrow{r}(t), \dot{\overrightarrow{r}}(t))\end{array}$ . En d'autres termes, les sept paramètres $(\mu, X_0, Y_0, Z_0, \dot{X}_0, \dot{Y}_0,\dot{Z}_0)$ définissent une orbite de manière univoque. On peut faire un changement de variables pour décrire l'orbite par un autre jeu de sept paramètres, par exemple $a, e, P, \omega, i, \Omega$ et t_p , qui sont les paramètres classiques présentés dans la section précédente. La plupart des auteurs les utilisent pour ajuster le mouvement des planètes, mais ils ont l'inconvénient d'être très sensibles aux erreurs pour de faibles excentricités et inclinaisons. Pour palier à ce problème on définit:

$\varpi = \omega + \Omega$ : la longitude du périhélie
$k = e \cos \varpi$
$h = e\sin \varpi$
: le demi-grand axe
: l'inclinaison
$\Omega$ : l'ascension droite au noeud ascendant
$\lambda = M + \varpi$ : la longitude

Ces éléments sont des fonctions continues de l'inclinaison et de l'excentricité.

Exercices

Auteur: Nathan Hara

Bonnes et mauvaises configurations d'observation

Difficulté : ☆☆

Question 1)

On définit un repère (x,y,z) où la direction est selon la droite observateur - étoile. Avec les notations des pages précédentes, pour quelles valeurs de l'inclinaison les observations de vitesses radiales ne donnent aucune information sur le système ? lequelles rendent la configuration optimale ? Même question pour les angles et $\Omega$ en astrométrie.

Auteur: Nathan Hara

Qu'est-ce que les vitesses radiales permettent de mesurer ?

Difficulté : ☆☆☆

Question 1)

On définit un repère (x,y,z) où la direction est selon la droite observateur - étoile. On note $\left( \begin{array}{c} \dot{X} \\ \dot{Y} \\ 0 \end{array} \right)$ les composantes du mouvement dans le plan orbital. Calculer la composante du mouvement de l'étoile selon la direction d'observation. On mettra les observations sous la forme $K ( \cos(\omega+\nu) + e\cos(\omega))$ , où est une constante à déterminer. Quels paramètres de l'orbite peut-on observer par vitesse radiale ? Montrer en particulier qu'on ne peut pas déterminer la masse exacte de la planète.

Modèle de la trajectoire

Les expressions établies jusqu'ici permettent de paramétrer le mouvement d'une étoile autour de laquelle une planète orbite, en faisant l'hypothèse que les deux corps sont ponctuels, forment un système isolé. Nous allons utiliser ces expressions pour donner un modèle général, lorsque n_p planètes orbitent autour de l'étoile.

En notant $\overrightarrow{r}_0$ et respectivement la position et la masse de l'étoile et $\overrightarrow{r_k}, m_k$ , $k = \mathbf{1}..n_p$ les positions et masses des n_p planètes , les équations de la mécanique (classique) dans le référentiel barycentrique de ce système sont:

$\begin{array}{cccccc} \ddot{\overrightarrow{r}_0} & = &-G\sum\limits_{k=1}^{n_p} m_k \frac{\overrightarrow{r_0}-\overrightarrow{r_k}}{\|\overrightarrow{r_0}-\overrightarrow{r_k}\|^3} & &=& -GM\left(\sum\limits_{k=1}^{n_p} \frac{m_k}{M} \frac{\overrightarrow{r_0}-\overrightarrow{r_k}}{\|\overrightarrow{r_0}-\overrightarrow{r_k}\|^3 \right) \\ \vdots & & & & & \\ \ddot{\overrightarrow{r}_j} & = &-GM\frac{\overrightarrow{r_j}-\overrightarrow{r_0}}{\|\overrightarrow{r_j}-\overrightarrow{r_0}\|^3}& -G\sum\limits_{k=1, k\not=j}^{n_p} m_k \frac{\overrightarrow{r_j}-\overrightarrow{r_k}}{\|\overrightarrow{r_j}-\overrightarrow{r_k}\|^3}& = &-GM \left(\frac{\overrightarrow{r_j}-\overrightarrow{r_0}}{\|\overrightarrow{r_j}-\overrightarrow{r_0}\|^3} - \sum\limits_{k=1,k\not=j}^{n_p} \frac{m_k}{M} \frac{\overrightarrow{r_j}-\overrightarrow{r_k}}{\|\overrightarrow{r_j}-\overrightarrow{r_k}\|^3} \right) \\ \vdots & & & & & \\ \ddot{\overrightarrow{r}_{n_p}} & = &-GM\frac{\overrightarrow{r_{n_p}}-\overrightarrow{r_0}}{\|\overrightarrow{r_{n_p}}-\overrightarrow{r_0}\|^3}& -G\sum\limits_{k=1}^{n_p-1} m_k \frac{\overrightarrow{r_{n_p}}-\overrightarrow{r_k}}{\|\overrightarrow{r_{n_p}}-\overrightarrow{r_k}\|^3} &=& -GM \left(\frac{\overrightarrow{r_{n_p}}-\overrightarrow{r_0}}{\|\overrightarrow{r_{n_p}}-\overrightarrow{r_0}\|^3} - \sum\limits_{k=1}^{n_p-1} \frac{m_k}{M} \frac{\overrightarrow{r_{n_p}}-\overrightarrow{r_k}}{\|\overrightarrow{r_{n_p}}-\overrightarrow{r_k}\|^3} \right) \end{array}$

En négligeant tous les termes du type $\frac{m_k}{M} \frac{\overrightarrow{r_j} - \overrightarrow{r_k}}{\|\overrightarrow{r_j} - \overrightarrow{r_k} \|^3}$ , c'est à dire l'intéraction entre les planètes, on obtient:

$\begin{array}{cccc} \ddot{\overrightarrow{r}_0} & =& -GM\left(\sum\limits_{k=1}^{n_p} \frac{m_k}{M} \frac{\overrightarrow{r_0}-\overrightarrow{r_k}}{\|\overrightarrow{r_0}-\overrightarrow{r_k}\|^3 \right) = -\sum\limits_{k=1}^{n_p} \frac{m_k}{M} \ddot{\overrightarrow{r_k}} \\ \vdots & & & \\ \ddot{\overrightarrow{r}_j} & = &-GM\frac{\overrightarrow{r_j}-\overrightarrow{r_0}}{\|\overrightarrow{r_j}-\overrightarrow{r_0}\|^3}& \\ \vdots & & & \\ \ddot{\overrightarrow{r}_{n_p}} & = &-GM\frac{\overrightarrow{r_{n_p}}-\overrightarrow{r_0}}{\|\overrightarrow{r_{n_p}}-\overrightarrow{r_0}\|^3}& \end{array}$

En résolvant les n_p problèmes à deux corps associés à chacune des planètes, par $M \overrightarrow{r}_0 + \sum\limits_{k=1}^{n_p} m_k \overrightarrow{r_k} = \overrightarrow{0}$ (dont la première équation du système ci-dessus est la dérivée seconde) on obtient le mouvement de l'étoile. Le modèle de trajectoire complet dans un référentiel galiléen est:

$\begin{array}{ccccccc}\overrightarrow{OS}(t) &=& \left( \begin{array}{c} x(t) \\ y(t) \\ z(t) \end{array}\right) &=& \overrightarrow{OB}_0 + \alpha t \overrightarrow{u} + \overrightarrow{BS}(t) & =& \overrightarrow{OB}_0 + \alpha t \overrightarrow{u} - \frac{1}{M}\sum\limits_{k=1}^{n_p} m_k \overrightarrow{r_k} \\ \dot{\overrightarrow{OS}}(t) &= &\left( \begin{array}{c} \dot{x}(t) \\ \dot{y}(t) \\ \dot{z} (t) \end{array}\right) &=& \alpha \overrightarrow{u} + \dot{\overrightarrow{BS}}(t)&=& \alpha \overrightarrow{u} - \frac{1}{M}\sum\limits_{k=1}^{n_p} m_k \dot{\overrightarrow{r_k}} \end{array}$

Trajectoire observée depuis le barycentre du système solaire

Le modèle précédent donne le mouvement de l'étoile observée dans un référentiel galiléen, en particulier le référentiel barycentrique du système solaire $(O,\overrightarrow{I},\overrightarrow{J},\overrightarrow{K})$ , que l'on munit d'un repère de coordonnées sphériques $(O, \rho, \alpha, \delta)$ repérant l'étoile . Supposons que l'observateur est situé en . A un instant il mesure:

soit la vitesse radiale $\dot{\|\overrightarrow{OS}\|}(t) = \dot{\rho}(t)$
soit la position sur la sphère céleste $(\alpha(t), \delta(t))$ dans le cas de mesures astrométriques

On peut facilement montrer avec un développement limité qu'au premier ordre en $\frac{\| \overrightarrow{BS} \|}{\| \overrightarrow{OS}\|}$ , $\rho(t) = \rho_B(t)+ \overrightarrow{e_{\rho}} \cdot \overrightarrow{BS}(t)$ , $\alpha(t) \cos \delta(t) = \alpha_B(t)\cos \delta(t)+ \frac{1}{\rho_B} \overrightarrow{e_{\alpha}} \cdot \overrightarrow{BS}(t)$ , de même $\delta(t) = \delta_B(t)+ \frac{1}{\rho_B} \overrightarrow{e_{\delta}} \cdot \overrightarrow{BS}(t)$ (à faire en exercice). Supposons que n_p planètes indéxées par gravitent autour de l'étoile en et $(\overrightarrow{I_k}, \overrightarrow{J_k}, \overrightarrow{K_k})$ le repère orbital de la planète . On projette le mouvement $\left( \begin{array}{c} X_k(t) \\ Y_k(t) \\ 0\end{array} \right)$ dans le repère orthonormé direct $(\overrightarrow{e_{\rho}}, \overrightarrow{e_{\alpha}}, \overrightarrow{e_{\delta}})$ associé aux coordonnées sphériques $(O, \rho_B, \alpha_B, \delta_B)$ (voir figure). Avec les notations de la page Changement de référentiel:

$\Large \left( \begin{array}{c} \alpha(t) \cos \delta_B \\ \delta(t) \\ \rho(t) \end{array} \right) = \left( \begin{array}{c} \alpha_B(t)\cos \delta_B \\ \delta_B(t) \\ \rho_B(t) \end{array} \right) - \frac{1}{M} \sum\limits_{k=1}^{n_p} m_k \left( \begin{array}{ccc} \frac{A_k}{\rho} & \frac{F_k}{\rho} & \star \\ \frac{B_k}{\rho} & \frac{G_k}{\rho} & \star \\ C_k & D_k & \star \end{array} \right) \left( \begin{array}{c} X_k(t) \\ Y_k(t) \\ 0 \end{array} \right)$

Observation depuis le barycentre du système solaire

Position de l'étoile en coordonnées sphériques dans le référérentiel barycentrique du système solaire et dans le référentiel translaté au centre de masse de la Terre. Les observations sont disponibles dans le référentiel terrestre, représenté en rouge, qui dépend de l'instant de mesure

Trajectoire observée depuis la Terre

Le modèle précédent donne l'évolution de l'étoile dans le repère barycentrique du système solaire (RBSS). Or les observations sont disponibles depuis la Terre. Si est la position de l'observateur, l'étoile cible on a $\overrightarrow{TS} = \overrightarrow{TO} + \overrightarrow{OS}$

La mesure de vitesse radiale est celle de $\dot{\overrightarrow{TS}}$ , il faut donc soustraire des observations $\dot{\overrightarrow{TO}}$
En astrométrie on observe la projection du mouvement de l'étoile sur la sphère céleste à un certain instant. Non seulement il faut estimer $\overrightarrow{TO}$ comme dans le cas des vitesses radiales, mais aussi s'assurer que toutes les observations sont exprimées dans le même référentiel (par exemple le référenciel barycentrique translaté, sur la figure) alors que chaque observation est dans un référentiel dépendant de l'observation (en rouge sur la figure). L'expression des mesures prises à des instants différents selon une même échelle se pose aussi pour les vitesses radiales, mais passe par l'étalonnage de l'instrument de mesure.

La dérivée temporelle des vecteur est toujours définie par rapport à un référentiel. Ici la notation $\dot{a}$ désigne la dérivée temporelle par rapport au RBSS.

La détermination de la position du centre de masse de la Terre par rapport au barycentre du système solaire est un sujet à part entière. La trajectoire d'un corps céleste au cours du temps dans un référentiel donné est appelée une éphéméride. Les principaux laboratoires de calcul des éphémérides sont le JPL (NASA) et l'IMCCE (Observatoire de Paris). Les liens envoient sur les générateurs en lignes d'éphémérides respectifs des deux laboratoires. .

Vitesses radiales

La vitesse de l'observateur par rapport au barycentre du système solaire peut se décomposer en $\dot{\overrightarrow{OT}} = \dot{\overrightarrow{OC}} + \dot{\overrightarrow{CT}}$ . La partie $\overrightarrow{OC}$ est donné par les éphémérides, comme expliqué plus haut. La correction de la deuxième partie est essentielle: l'observateur parcourt deux fois le rayon terrestre en une nuit, ce qui donne une vitesse d'environ 300 m/s. Si M(t) désigne la matrice de changement de repère entre le référentiel lié à la Terre et le référentiel barycentrique du système solaire, on a $\dot{\overrightarrow{CT}} = \dot{M}(t) \overrightarrow{CT}_{R_t}$ où $\overrightarrow{CT}_{R_t}$ désigne la position de l'observateur dans le référentiel terrestre.

La détermination de est aussi un sujet à part entière, appelé "rotation de la Terre". Les paramètres de rotations officiels sont donnés par le Service de rotation de la Terre au SYRTE (Observatoire de Paris).

Parallaxe (astrométrie)

En ce qui concerne l'astrométrie, il faut exprimer la relation entre les positions mesurées et la position dans le RBSS, ce qui peut se décomposer en deux étapes: passer du référentiel terrestre au RBS translaté au centre de masse de la Terre (passer du référentiel rouge au référentiel noir à droite sur la figure), puis passer du référentiel translaté au RBS. Comme on le verra plus tard, il est aussi possible de se passer de cette étape en prenant un champ contenant des étoiles de référence.

Notons (x,y,z) la position de la Terre dans le RBS. En exprimant $\overrightarrow{OS}$ dans le RBSS de deux manières on obtient une relation entre les coordonnées $(\rho_B, \alpha_B, \delta_B)$ et $(\rho_T, \alpha_T, \delta_T)$ : $\begin{array}{lll} \rho_T \cos \delta_T \cos \alpha_T & = & \rho_T \cos \delta_S \cos \alpha_S -x\\ \rho_T \cos \delta_T \sin \alpha_T & = & \rho_T \cos \delta_S \sin \alpha_S - y \\ \rho_T \sin \delta_T & = & \rho_B \sin \delta_B -z\end{array}$

Lorsque est suffisamment loin, ces expressions différenciées au voisinage de $(\rho_B, \alpha_B, \delta_B)$ donnent au premier ordre en $\Delta \alpha = \alpha_T - \alpha_B$ , $\Delta \delta = \delta_T - \delta_B$ et $\Delta \rho = \rho_T - \rho_B$

$\begin{array}{lll} \Delta \alpha \cos \delta &=& \frac{x}{\rho_B} \sin \alpha - \frac{y}{\rho_B} \cos \alpha \\ \Delta \delta &=& \left( \frac{x}{\rho_B} \cos \alpha + \frac{y}{\rho_B} \sin \alpha \left) \sin \delta - \frac{z}{\rho_B} \cos \delta \\ \Delta \rho &= &-x \cos \alpha \cos \delta - y\sin \alpha \cos \delta + z \sin \delta \end{array}$

On appelle la quantité $\varphi = \frac{1}{\rho_B} \approx \frac{1}{\rho}$ la parallaxe de l'étoile (en général notée $\pi$ ou $\varpi$ , qui sont des symboles déjà utilisés dans le cours). Elle considérée comme constante au cours des observations et est ajustée aux observations en astrométrie.

Remarque: la procédure de changement de référentiel passe par des changements d'échelle de temps (UTC, UT1, TDB...) qui ne seront pas détaillés ici.

Changement de coordonnées

Accélération de perspective et changement de parallaxe (astrométrie)

Accélération de perspective

L'accélération de perspective est un effet purement géométrique, qui tient à la définition du mouvement propre d'une étoile. Le mouvement propre est en effet défini comme $\mu = \frac{V\sin \theta}{\rho}$ où est la vitesse du système observé, $\rho$ la distance entre l'observateur et le système, $\theta$ est l'angle entre la ligne de visée et la vitesse du système observé. Nous supposons que les mesures astrométriques sont sur un plan. Cependant, on mesure la projection du mouvement sur une sphère, donc $\theta$ n'est pas constant au cours du temps. En dérivant par rapport au temps, comme est constant:

$\frac{d \mu }{dt} = - \frac{V}{\rho^2} \sin \theta \frac{d\rho}{dt}+\frac{V}{\rho} \cos \theta \frac{d\theta}{dt}$

On voit sur le triangle OSP que $\frac{d\theta}{dt} = - \mu$ et $V \cos\theta =V_r=\frac{d\rho}{dt}$ . D'où:

$\frac{d\mu}{dt} = - \frac{\mu}{\rho} \frac{d\rho}{dt} - \frac{\mu}{\rho} V_r = -\frac{2\mu}{\rho} V_r$

En pratique on ajuste un terme quadratique $\eta$ , où le mouvement angulaire sur la sphère céleste est donné par $\mu (t) = \mu_0 t + \eta t^2$ .

Variation de parallaxe

La variation de parallaxe vaut: $\dot{\varphi}= -\frac{1}{\rho^2} \frac{d\rho}{dt} = -\frac{1}{\rho^2} V_r= - \varphi^2 V_r$ . Comme elle est du deuxième ordre en $\varphi$ , elle n'est en général pas prise en compte.

Changement de coordonnées

Définition de $\theta$

Etoile

L'ajustement des paramètres orbitaux permet de connaître $\mu = G(m+M)$ ou bien $\mu \sin i$ dans le cas des vitesses radiales. Il est impossible de distinguer les masses et séparément a priori. Cependant, on peut mesurer la masse de l'étoile de sorte à lever l'indetermination.

La modélisation des étoiles permet de distinguer les effets sur le spectre de la variation de leur flux lumineux de la présence de compagnons planétaires. Ces variations sont modélisées par des variables aléatoires suivant une certaine loi de probabilité, dont l'amplitude varie de plusieurs ordres de grandeur selon le type d'étoile. Plus l'étoile est active, plus l'amplitude du mouvement dû à la planète doit être grande pour distinguer la planète du bruit. En particulier certains types d'étoiles sont trop actives pour pouvoir détecter des potentielles super-Terres compagnon.

La description physique des étoiles est complexe car de nombreux phénomènes, tous interdépendants, ont lieu: convection, radiation, magnétisme... Les modèles utilisés en détection par vitesses radiales comprennent trois phénomènes

Les ondes de pression se propageant dans les couches extrérieures de l'étoile entrainent des compressions-dilatations de celles-ci, ce qui change la vitesse radiale moyenne observée par rapport à une photosphère statique.Cet effet affecte principalement les observations de vitesses radiales.
La granulation: phénomène convectif (c'est à dire lié au déplacement de matière dans l'étoile)
L'activité de l'étoile: la présence de régions sombres ou lumineuses à la période de rotation de l'étoile

L'impact de ces effets sur les mesures n'est pas simple à quantifier. Selon le type d'étoile et les instants d'observations, l'effet peut être très variable. Certains auteurs estiment les incertitudes par des simulations numérique, d'autres par des modèles théoriques. Nous donnons une description qualitative de ces phénomènes et une modélisation possible, mais il faut garder à l'esprit que c'est un sujet de recherche ouvert. Dans le cas des vitesses radiales, ces bruits sont classiquement modélisées par un processus stochastique d'une certaine densité spectrale de puissance, notion importante en statistique, qui est définie dans cette page.

Masse de l'étoile

La mesure de la masse de l'étoile peut se faire de deux manières

Par ajustement aux observations de modèles théoriques d'intérieurs stellaires.
Par mesure du mouvement dans un système binaire

Dans le premier cas, si la distance à l'étoile est connue (par exemple par mesure de parallaxe), on peut mesurer sa luminosité intrinsèque (si on ne connaît pas la distance on ne mesure évidemment que la luminosité apparente). Par son spectre, on peut mesurer sa température effective $T_{eff}$ . Des modèles d'intérieurs stellaires permettent ensuite d'évaluer la masse. Cette estimation peut être rafinée avec un modèle d'atmosphère stellaire. On peut alors avoir la gravité à la surface de l'étoile . Comme $L \approx T_{eff}^4 R^2$ et $g \approx \frac{M}{R^2}$ , on peut avoir une estimation de la masse.

Dans le cas des étoiles binaires (systèmes de deux étoiles), le spectre présente des raies des deux étoiles. Le mouvement de ces raies se fait à la même fréquence, mais dans des directions opposées (lorsqu'une étoile approche l'autre s'éloigne). La période de ces mouvements est liée à la masse du système par l'équation de Kepler. L'amplitude relative de ces mouvements permet de déterminer la masse des deux étoiles séparément. Comme dans le cas de la détermination des orbites des planètes, la masse n'est connue qu'à un faceteur $\sin i$ près. Pour lever cette indetermination, il faut déterminer l'inclinaison de l'orbite par rapport à l'observateur . Si on observe des eclipses (une binaire passe devant l'autre), $\approx 0$ . Si le système n'est pas dans cette configuration, on peut séparer angulairement les deux étoiles par des techniques d'interférométrie.

Processus stochastique et Densité spectrale de puissance

Processus stochastique

La notion de densité spectrale de puissance (DSP) n'est pas simple à définir, cependant très utilisée dans la littérature de traitement du signal. Nous donnons une définition mathématique pour qu'il n'y ait pas d'ambiguités mais compte tenu de la sophistication des notions introduites, le lecteur pourra se référer à la description qualitative suivante.

La densité spectrale de puissance est une propriété relative à plusieurs variables aléatoires. Les familles de variables aléatoires peuvent par exemple représenter des mesures sur lesquelles on a une incertitude. A chaque instant de mesure on associe une variable alétoire qui a une certaine densité de probabilité. En physique théorique ou en économie, on rencontre des processus stochastiques continus - typiquement le mouvement brownien, qui représente des mouvements d'atomes ou des fluctuations de prix. Formellement, un processus stochastique est une famille de variables aléatoires indexées par un ensemble totalement ordonné , toutes définies sur le même espace de probabilité ( $\forall t \in T, X_t \in (\Omega,\mathcal{F},\mathbb{P}) )$ . Dans ce cours on aura seulement besoin de $T=\mathbb{R}$ ou $\mathbb{N}$ . On note $\mathbb{E}$ l'espérance mathématique.

Dans le cas général, la densité de probabilité de la variable aléatoire $X_{t}$ (pour $t \in \mathbb{R}$ ) dépend des valeurs prises à d'autres "instants" par les autres variables aléatoires. En particulier on peut s'intéresser à une éventuelle probabilité de périodicité. Par exemple si on modélise un nombre de ventes de vêtement par jour, on verra des ventes plus importantes au moment des soldes (à peu près tous les six mois). La densité spectrale de puissance est un outil qui permet de visualiser ce genre de périodicité. Dans la section suivante, on voit que si on prend une famille de variables aléatoires certaines, c'est à dire que X_t vaut une certaine valeur réelle x(t) avec la probabilité 1, la DSP en une fréquence $\omega$ est égale à $|\hat{x}(\omega)|^2$ , où $\hat{x}$ est la transformée de Fourier de . Si maintenant X_t est une variable alétoire, la DSP sera la "transformée de Fourier typique" d'une réalisation de X_t .

Pour définir cette notion mathématiquement, on doit d'abord introduire les notions de convergences et intégrales en moyenne quadratique. Pour plus de précision le lecteur peut se référer au cours de Timo Koski à KTH.

Intégrale en moyenne quadratique (Mean Square Integral)

Rappelons d'abord que si Y_1, ..., Y_n sont des variables aléatoires et $\phi: \mathbb{R}^n\rightarrow \mathbb{R}^m$ une fonction mesurable alors $\phi(Y_1,...,Y_n)$ est une variable aléatoire. En particulier, si $\lambda$ est un scalaire, $\lambda Y_1$ et Y_1 + Y_2 +...+Y_n sont des variables aléatoires. Pour un rappel sur les variables aléatoires, voir le cours de Didier Pelat.

Soit $(\Omega,\mathcal{F},\mathbb{P})$ un espace de probabilités, on dit que la suite de variables aléatoires $(Y_n)_{n \in \mathbb{N}}$ telle que $E\{Y_n\} < + \infty$ , définies sur cet espace converge en moyenne quadratique si et seulement si:

$\lim\limits_{n \rightarrow \infty} \mathbb{E}\{|Y_n-Y |^2\}=0$

Soit $(X_t)_{\mathbf{t \in \mathbb{R}}}$ un processus stochastique continu ( $T=\mathbb{R}$ ) tel que chacune des variables aléatoires X_t a une espérance finie ( $E\{X_t\} < + \infty$ ). L'intégrale en moyenne quadratique du processus $(X_t)_{t \in T}$ sur l'intervalle [a,b] est définie comme la limite en moyenne quadratique (lorsqu'elle existe) de:

$\sum\limits_{k=1}^n X_{t_k} (t_k-t_{k-1}})$

Pour a=t_0 < t_1 <... <t_n=b et $\lim\limts_{n\rightarrow\infty}\max (t_k-t_{k-1}) = 0$ . On la note alors $\int_b^b X_t dt$ .

On définit alors la densité spectrale de puissance comme:

$P_X(\omega) = \lim\limits_{T\rightarrow \infty} \mathbb{E} \{|\hat{x}_T(\omega)|^2\}$

Où $\hat{x}_T(\omega) = \frac{1}{\sqrt{T}} \int_0^T X_t e^{-i \omega t} dt$

Cette définition un peu complexe peut être vue comme une généralisation de la transformée de Fourier à des processus stochastiques. En effet, lorsque le processus X_t est telle que X_t = x(t) avec une probabilité 1, l'intégrale en moyenne quadratique se comporte comme l'intérale de Riemann, alors $P(\omega)$ est le carré du module de la transformée de Fourier de la fonction réelle d'une variable réelle x(t) . Dans le cas où les X_t sont aléatoire, est le carré de la transformée de Fourier "en moyenne" des réalisations de X_t. Par exemple si X_t modélise une tension mesurée au cours du temps dans une expérience d'électronique réalisée un grand nombre de fois, donnant profils de tension u_k(t) à l'expérience (des réalisations du processus stochastique U_t ), la moyenne des carrés du module des transformées de Fourier des u_k(t) notée $\tilde{P}_n(\omega)$ sera approximativement égal à $P_U(\omega)$ . Si le nombre d'expérience tend vers l'infini $\tilde{P_n} \rightarrow P_U$ en norme 2.

Dans le cas d'un processus stationnaire discret ( $T= \mathbb{Z}$ ), on peut directement définir $\hat{x}_n(\omega)=\frac{1}{\sqrt{2n+1}} \sum\limits_{k=-n}^n X_ke^{-i\omega t_k}$ et $P_X(\omega) = \lim\limits_{n\rightarrow \infty} \mathbb{E} \{|\hat{x}_n(\omega)|^2\}$ .

Processus stationnaires au sens large

La densité spectrale de puissance a une définition plus simple lorsque le processus est stationnaire, c'est à dire lorsque le processus $(X_t)_{t\in T}$ vérifié:

Il existe $\mu \in \mathbb{R}$ tel que $\forall t \in T, \mathbb{E}\{X_t\} =\mu$
La fonction d'autocorrélation $\begin{array}{cccc}C: &\mathbb{R} \times \mathbb{R}: & \rightarrow & \mathbb{R} \\ &t,s & \rightarrow& \mathbb{E}\{X^{\star}_s X_t \} \end{array}$ ne dépend que de , i. e. $\forall t,s, t',s'$ tels que , on a . On note appelle alors habituellement la fonction d'autocorrélation , telle que

Les processus stationnaires modélisent des phénomènes qui ont une certaine invariance dans le temps, en particulier la covariance ne dépend pas de de manière absolue, mais de manière relative à un autre instant. Dans ce cas, la densité spectrale de puissance est égale au carré du module de la transformée de Fourier de la fonction . L'équivalence avec la définition de la densité spectrale de puissance donnée plus haut est établie par le théorème de Wiener-Khinchin.

Granulation

Le phénomène de granulation est lié à la convection du gaz dans l'étoile. La lumière rayonnée par le gaz chaud remontant à la surface va vers l'observateur, la longueur d'onde reçue est donc décalée vers le bleu. En rayonnant, le gaz se refroidit, puis repart vers le centre de l'étoile. Etant moins chaud, il émet moins de lumière, si bien que la lumière est globalement décalée vers le bleu. Ce phénomène est variable dans le temps, donc le décalage vers le bleu aussi. Cette variation peut apparaître dans le spectre et créer des fréquences parasites. La nature aléatoire du phénomène fait que même après ajustement, il reste un bruit résiduel. Pour une étoile de type solaire, il est de l'ordre de 0.5 - 1 m/s sur une observation.

La granulation est en général modélisée par un processus stochastique dits de bruits en créneaui (popcorn noise ou burst noise en anglais). Il s'agit de processus stochastiques pouvant prendre deux valeurs, par exemple -1 ou 1 avec une probabilité de changement suivant une loi de Poisson (loi exponentielle). Si à la valeur passe de 1 à -1, la densité de probabilité pour que la valeur passe à 1 à $t+ \Delta t$ est $e^{-\frac{t}{\tau}}$ où $\tau$ est un réel positif. Pour les vitesses radiales, la densité spectrale de puissance de ces bruits peut être modélisée par::

$P(\nu) = \frac{4 \sigma ^2 \tau}{1+(2\pi\nu\tau)^2}}$

Cette modélisation, due à Harvey (1985) a depuis été revue et d'autres densités spectrales de puissances ont été proposée à partir de simulations 3d de convection au sein d'une étoile. En pratique, le bruit dû à la granulation apparaitra comme un signal périodique de l'ordre de cinq minutes. Cependant, on observe aussi des phénomènes appelés meso-granulation et super-granulation sur des échelles de temps plus longues. La contribution totale de ces bruits est:

$P(\nu) = \frac{4 \sigma_g ^2 \tau_g}{1+(2\pi\nu\tau_g)^2}}+\frac{4 \sigma_{mg} ^2 \tau_{mg}}{1+(2\pi\nu\tau_{mg})^2}}+\frac{4 \sigma_{sg} ^2 \tau_{sg}}{1+(2\pi\nu\tau_{sg})^2}}$

Où les indices g, mg et sg se réfèrent respectivement à la granulation, la méso granulation et la super-granulation. En anticipant un peu sur le troisième chapitre, lorsque ces bruits sont pris en compte, les valeurs des $\sigma$ et $\tau$ sont ajustés sur le spectre de puissance du signal.

Activité magnétique

La formation d'arcs de champ magnétique à la surface de l'étoile inhibe le mouvement des particules, donc réduit la température et provoque donc des tâches sombres. Cet effet à des effets à court termes (à la fréquence de rotation de l'étoile $\approx 1$ mois), et à plus long terme à travers des cycles d'activité magnétique.

A court trerme, la tache introduit une dissymétrie entre la partie de l'étoile tournant vers l'observateur, et la partie s'en éloignant, ce qui engendre un décalage du spectre mesuré. D'autre part, la tache engendre un déplacement du photocentre de l'étoile périodique, pouvant être confondu avec la présence d'une planète. Pour éviter ces confusions, on estime la période de rotation de l'étoile par spectroscopie, et on ajuste des sinusoïdes à cette période et ses premières harmoniques.

L'activité magnétique peut se mesurer à travers divers indicateurs, dont on analyse les corrélations.

Toujours selon le modèle de Harvey (1985), la densité spectrale de puissance du bruit de vitesse radiale induit par une tache solaire est:

$P(\nu) = \frac{4 \sigma_{am} ^2 \tau_{am}}{1+(2\pi\nu\tau_{am})^2}}$

Effet de la rotation de l'étoile sur le spectre

La lumière provenant de la moitié de l'étoile ayant un mouvement vers l'observateur est décalé vers les hautes fréquences (vers le bleu). L'autre moitié est décalée vers les basses fréquences (vers le rouge). Dans l'hypothèse où l'étoile est sphérique, et a une luminosité identique partout sur sa surface, le décalage vers le rouge et celui vers le bleu ne fait qu'élargir les raies spectrales. Si une tache est présente, ici sur la partie bleue, la symétrie est brisée et le déficit de lumière entraine un décalage du spectre vers le rouge.

Oscillations radiales (vitesses radiales)

Des inhomogénéités de densité dues aux mouvements convectifs font que des ondes mécaniques se propagent au sein des étoiles. Certaines de ces ondes sont radiales, ce qui provoque un mouvement d'ensemble de la photosphère qui a une signature sur le décalage du spectre mesuré. L'étude de ces ondes est un domaine de la physique stellaire appelé "astérosismologie". Etant donné que la théorie est accessible au niveau licence, nous en donnons des principes généraux.

La théorie procède comme suit: on écrit localement 1) l'équation du mouvement linéarisée au premier ordre au voisinage d'un état d'équilibre, 2) la conservation de la masse ou équation de continuité, 3) l'équation de Poisson, liant le potentiel gravitationnel et la densité, 4) Le premier principe de la thermodynamique. On néglige l'effet du champ magnétique. En général, on fait l'hypothèse que le terme de transfert thermique dans le 1er principe est nul.

Perturbations atmosphériques

Vitesses radiales

L'atmosphère peut influencer la mesure de vitesses radiales de deux manières:

Par atténuation de l'intensité, ce qui augmente le bruit de photon
Par la modification du contenu spectral du signal

La présence de lignes spectrales d'émission ou d'absorption est difficile à corriger. C'est pourquoi on ne considère que des plages de fréquences où l'intensité des raies atmosphériques est inférieure à 1/10000 de l'intensité de la cible. De plus, la réponse de l'atmosphère dépend de la longueur d'onde. Le spectre obtenu est pondéré de sorte à corriger ces inhomogénéités.

Astrométrie

Avant de parvenir au télescope, la lumière issue de l'étoile traverse l'atmosphère. Les turbulences aux hautes altitudes créent des inhomogénéités de densité, donc d'indice optique, qui distordent le front d'onde de la lumière pénétrant dans l'atmosphère. La modélisation classique de ce phénomène passe par l'introduction de "Structure functions" qui modélisent la densité de probabilité des champs de vitesses, densité etc.

Afin de ne pas surcharger le cours, nous ne donnerons qu'un modèle très simple de ce phénomène. Les mouvements de turbulence forment des "cellules" de composition homogène et d'une taille typique d_0 . Vue du dessus, l'atmosphère se comporte approximativement comme un tableau de pixels de taille d_0 dont chaque cellule introduit une déviation angulaire du faisceau incident. Cette déviation est de l'ordre de $\approx \frac{\lambda}{d_0}$ . Dans le cas limite où une cellule a un indice 0 et les cellules avoisinantes ont un indice 1 (sont opaques), on retrouve l'ordre de grandeur de la diffraction. Les fronts d'ondes issus de ces différents pixels arrivent au télescope avec des angles différents, donc donnent chacun une image différente appelé "speckle". L'union de ces speckle forme une tache lumineuse élargie de taille environ égale à $1.22\frac{\lambda}{d_0}$ .

Récapitulation

Les effets expliqués précédemment sont pour la plupart périodiques, donc peuvent être confondus avec le signal d'une planète. Le tableau suivant récapitule ces effets, leurs amplitudes et périodes typiques pour des étoiles de type solaire.

Effets physiques
Effet	Description	Amplitude (vitesse radiale)	Amplitude (astrométrie)	Echelle de temps
Mouvement des planètes	Le mouvement des planètes autour de l'étoile engendre un mouvement périodique	1 cm/s - 200 m/s	Terre: 0.3, 0.03 $\mu$ as Jupiter: 500, 20 $\mu$ as à 1 et 10 parsec resp.	de un jour à plusieurs centaines d'années
Mouvement propre	Mouvement rectiligne uniforme	Dizaines de kilomètres par seconde	10 - 1000 mas pour les étoiles observables	Mouvement non périodique
Taches	La présence de taches sombres ou lumineuses dues à l'activité magnétique provoque une disymétris entre la luminosité de la partie bleue et la partie rouge de l'étoile, engendrant un décalage du spectre ou du photocentre.	1 m/s	Avec ces trois effets, 0.5 - 10 $\mu$ UA	Période de rotation de l'étoile (quelques jours)
Activité magnétique (long terme)	Le nombre de tâche à la surface de l'étoile peut varier de quelques unes à plusieurs centaines. Par exemple le soleil a une périodicité de 11 ans. L'effet précédent est modulé par ces variations à long terme.	10 m/s		1 - 10 ans
Granulation	Le mouvement convectif à l'intérieur de l'étoile provoque un mouvement du gaz dans la photosphère	0.5 - 1 m/s		Granulation: quelques minutes, mesogranulation:
Oscillations radiales (p-modes)	La propagations d'ondes acoustiques dans le manteau de l'étoile entraîne une oscillation de celui-ci	0.5 - 1 m/s	inconnu	5 - 10 minutes
Système multiple	La présence d'autres étoiles orbitant autour de l'étoile cible engendre un mouvement décrit par les mêmes équations que celles utilisées pour les planètes. Les autres étoiles étant beaucoup plus massives que des planètes, l'effet est plus important.	1 - 30 km/s	0.1 - 1 as/an	10 - 100 ans

effets observationnels
Effet	Description	Amplitude (vitesse radiale)	Amplitude (astrométrie)	Echelle de temps
Mouvement de l'observateur dans le RBSS	Comme l'observateur est en mouvement dans un référentiel galiléen,	30 km/s (après soustraction, on a une erreur $\approx$ 0.5 m/s à cause des incertitudes sur la rotation de la Terre et les éphémérides)	0.01 - 0.1 as/an	Une année
Perturbations atmosphériques	La présence de raies spectrales atmosphériques peut perturber les observations par vitesses radiale, et les turbulences dévient les fronts d'ondes, engendrant un déplacement apparent des sources.	0.5 m/s	1 mas	Quelques minutes

Instrumentation & Observations

Objectifs

Jusqu'à présent, on a modélisé la lumière qui parvient à un observateur théorique situé au voisinage de la Terre. Dans cette section, on présentera l'aspect concret de l'observation, en particulier les instruments. En astrométrie on mesure la position sur le detecteur CCD sur le plan focal du télescope, pour les vitesses radiales on utilise un spectrographe dont l'entrée est située au point focal du télescope . Il est essentiel de bien comprendre le fonctionnement des instruments pour des raisons évidentes: ils constituent notre seule source d'information et leur étude permet de mieux modéliser les mesures, donc d'exploiter au mieux les données.

Dans le cas de l'astrométrie comme celui des vitesses radiales, deux conditions sont nécessaires pour obtenir des mesures suffisamment significatives:

La quantité de lumière collectée doit être suffisament grande pour obtenir le signal sur bruit désiré
Les erreurs instrumentales doivent être suffisamment petites pour avoir un rapport signal sur bruit acceptable.

En astrométrie comme en détection par vitesses radiales on place un dispositif sur le plan focal d'un télescope, respectivement des capteurs CCD et l'entrée d'un spectromètre (ou spectrographe). C'est pourquoi les principes généraux des télescopes seront présentés. On évoquera le fonctionnement des spectrographes utilisés pour les détections par mesures de vitesses radiales.

Spectrographe ELODIE

Cet appareil a permis la détection de 51 Pegasi b en 1995 par Michel Mayor et Didier Queloz à l'Observatoire de Haute Provence

Crédit : Observatoire de Haute Provence

Télescope

Un télescope est un appareil permettant de recueillir un rayonnement eléctromagnétique. Pour observer un rayonneme nt dans le visible ou dans l'infrarouge proche, on utilise des télescopes à miroir parabolique. Pour éviter les ambiguités, on définira le plan focal comme le plan perpendiculaire à l'axe optique (ici l'axe de symétrie du télescope) passant par le foyer, et on fait l'hypothèse que les rayons reçus font un angle faible avec l'axe optique. Dans ces conditions, La relation donnant la distance au point focal de l'image d'un rayon arrivant avec un angle $\alpha$ sur le plan focal est en première approximation $FF_1 \approx f\tan \alpha \approx f \alpha$ , où est la distance focale. Pour les angles faibles, on peut travailler avec une lentille équivalente au télescope, de même diamètre et distance focale. On va introduire trois notions de bases sur les télescopes: le champ, la résolution angulaire et la vitesse d'acquisition.

Le champ est la portion du ciel observée par le détecteur du télescope. Comme le détecteur est au plan focal, il est égal à $\tan \frac{s}{f} \approx \frac{s}{f}$

La résolution angulaireest l'angle minimal entre deux sources permettant de les séparer par le système de détection. Cette définition est vague, et Supposons qu'une source ponctuelle émettant à une longueur d'onde $\lambda$ soit placée en un point du plan focal . A cause du phénomène de diffraction, la lumière ne sera pas émise selon une direction unique, mais son énergie sera répartie sur certains angles centrés sur $\alpha = \frac{FP}{f}$ . Comme la lumière suit le même trajet dans les deux sens, ce détecteur reçoit de la lumière provenant de ces angles.

Le rapport du diamètre et de la distance focale du télescope donne la "vitesse" de l'instrument. En effet, considérons une source circulaire de taille angulaire $\alpha$ , son image sur le plan focal est un cercle d'aire $4\pi(f\alpha)^2$ . L'intensité observée est proportionnel à l'aire du télescope, donc l'énergie par unité de temps reçue est proportionnelle à $\left(\frac{d}{f}\right)^2$ . On définit l'ouverture du télescope par $R = \frac{f}{d}$ Le rapport signal sur bruit des mesures de CCD est égal à $\frac{1}{\mu \Delta t}$ où $\mu$ est le nombre moyen d'électron par unité de temps et $\Delta t$ est le temps d'intégration (voir Bruit de photon). On peut calculer le temps d'intégration minimal pour obtenir un certain signal sur bruit compte tenu de l'ouverture du télescope, de la taille de la source et de son intensité. Plus l'aire et grande, plus la distance focale est petite, et plus l'instrument collecte rapidement le nombre de photons requis.

Cependant, lorsque l'ouverture augmente, la résolution angulaire diminue.

Miroir parabolique

Les rayons perpendiculaires à l'axe de symétrie d'une parabole parviennent au point focal

.Ceux qui arrivent avec un certain angle arrivent légèrement décalés.

CCD

Un détecteur CCD est composé de trois partie: des électrodes en silicium, un isolant en SiO_2 et une jonction de semiconducteurs NP. Les électrodes en polysilicium sont reliées périodoquement par un fil conducteur (toutes les trois électrodes sur la figure) de sorte à créer un profil de potentiel électrique alternant puits et régions plates. Le principe d'un capteur CCD est le suivant (voir figure)

Lorsqu'un photon percute un atome du slicium dopé N, un électron est émis par effet photo-électrique.
Une jonction de semi-conducteurs dopés crée un champ électrique (exactement comme dans une diode). Comme la force qui s'exerce sur l'électron est $\overrightarrow{F} = - e \overrightarrow{E}$ où $\overrightarrow{E}$ est le champ électrique, est la charge de l'électron, et que le champ est orienté vers le bas l'électron est piégé sous l'isolant.
L'électron est confiné dans un puit de potentiel électrique, selon l'endroit où il a été émis
Ces trois premières étapes ont lieu pendant tout le temps d'intégration. Ensuite, le nombre d'électron par puit de potentiel est compté, voir (b) sur la figure. Les puits de potentiels sont progressivement décalés et le nombre d'électrons s'y trouvant est compté par puit.

Le contrôle de l'erreur induite par la CCD est primordial et peut s'avérer très complexe. Nous ne rentrerons pas dans ces considérations.

Schéma de principe d'un capteur CCD

Bruit de photon

Le bruit de photon est une limitation fondammentale en astronomie car il est dû à la nature de l'émission de la lumière. En effet, les instants d'émission de deux photons sont indépendants les uns des autres.On peut montrer mathématiquement qu'une suite d'évènement indépendants, sans mémoire et stationnaire est nécessairement un flux de Poisson. Dans notre cas, si on note le temps d'attente entre deux émission de photon, la probabilité que soit supérieure à t+s sachant que T>s est $e^{-\mu t}$ où $\mu$ est une constante.

Cette constante a une interprétation physique: plus elle est grande, plus $e^{-\mu t}$ diminue rapidement à donné, c'est à dire plus la probabilité que soit long diminue. On peut montrer $\mu$ est le nombre d'évènement par unité de temps moyen, proportionnelle à l'intensité. Rappelons que le détecteur capte des photons pendant le temps d'intégration $\Delta t$ . Le nombre de photons reçus pendant ce temps est une variable aléatoire suivant la loi:

$\text{Pr}\{N=n\} = \frac{\mu \Delta t}{n!} e^{-\mu \Delta t}$

Qui est appelée loi de Poisson. La moyenne et la variance d'une telle loi sont égales et valent $\mu \Delta t$ . En conséquence, l'écart-type vaut $\sqrt{\mu \Delta t}$ donc le rapport signal sur bruit (Signal-to-Noise Ratio) est:

$SNR = \frac{\text{ecart-type}}{moyenne} = \frac{\sqrt{\mu \Delta t}}{\mu \Delta t} = \frac{1}{\sqrt{\mu \Delta t}}$

Exercices

Auteur: Nathan Hara

Magnitude et temps d'observation

Difficulté : ☆☆

Question 1)

Au deuxième siècle avant J.-C., l'astronome Hipparque a classé cinquante étoile par ordre de brillance en 6 catégories, les plus brillantes occupant la première. L'échelle de magniture apparente moderne est sous la forme $m = -2.5 \log_{10} \frac{F}{F_0}$ où est le flux lumineux reçu et un flux de référence. Le choix de 2.5 et $F_0 = 3.52 \times 10^{-23}$ W $m^{-2}$ $Hz^{-1}$ fait que les étoiles de classe 1 d'Hipparque aient une magnitude comprise entre 0 et 1, les classes 2 on une magnitude entre 1 et 2 etc. De ce fait, il semble naturel que les objets les moins lumineux visibles à l'oeil nu soient environ de magnitude 6. On considère un télescope de diamètre , on suppose qu'on observe dans le visible $\Delta \lambda = [400,800]$ nm, et que le flux lumineux reçu ne dépend pas de la longueur d'onde sur cette plage.

Sachant qu'un photon a une énergie $h \frac{c}{\lambda}$ où $\lambda$ est la longueur d'onde et la vitesse de la lumière, montrer que le flux de photon reçus $\mu$ (photons par seconde) est lié au flux par $\mu = \pi \frac{D^2}{4} \frac{\Delta \lambda}{hc} F$
Donner l'expression du temps d'intégration pour obtenir un rapport signal sur bruit en fonction de la magnitude observée, de et $\Delta \lambda$
Le bruit de photon est il plus grand pour les étoiles peu lumineuses ou très lumineuses ?

Effet Doppler

On considère une source d'ondes et un observateur. Selon leur vitesse relative, la fréquence reçue par l'observateur varie: c'est ce qu'on appelle l'effet Doppler. Dans le cadre de la physique classique et d'une onde harmonique (purement sinusoidale), on peut le calculer simplement. En effet, notons la distance entre la source et l'observateur, la vitesse de l'onde et. Par définition d'une onde harmonique de pulsation $\omega$ et d'amplitude , l'amplitude mesurée à et vaut

$A\cos \left( \omega\left(t-\frac{x}{c}\right) + \phi\right)$

Où $\phi$ donne en particulier la phase en t=x=0 . Si la distance varie avec le temps selon x=x_0+vt , alors en et on mesure

$A\cos \left( \omega\left(1-\frac{v}{c}\right)t - \omega\frac{x_0}{c} + \phi\right) = A\cos \left( \omega't + \phi'\right)$

Donc pour l'observateur, tout se passe comme s'il recevait une onde de pulsation $\omega' = \omega \left(1-\frac{v}{c} \right)$ ou de longueur d'onde $\lambda' = \frac{1}{1-\frac{v}{c} } \lambda \approx \left(1+\frac{v}{c} \right)\lambda$ et de phase $\phi' =- \omega\frac{x_0}{c} + \phi$ .

Cependant, lorsque est la vitesse de la lumière, la relativité générale donne une description beaucoup plus précise du phénomène. Nous ne présenterons pas le calcul menant à l'expression de l'effet Doppler dans ce cadre mais en donnons l'expression:

$\lambda' = \lambda \frac{1 + \frac{v}{c} \cos \theta}{\sqrt{1-\frac{v^2}{c^2}}}$

Où $\theta$ est l'angle entre la direction de la vitesse relative et la ligne de visée observateur-source. Remarquons que lorsque $\theta$ tend vers $\pi/2$ et v/c tend vers 0, l'expression tend vers l'expression classique.

Spectrographe

L'observation des vitesses radiales nécessite de mesurer des longueurs d'onde très précisément, pour cela on utilise des spectrographe. Les équipes américaines et européennes utilisent des appareils différents, mais dans les deux cas ce sont des spectrographes d'échelle. Le principe d'un tel instrument est d'observer simultanément plusieurs ordres élevés de diffraction à l'aide de deux diffractions successives. La lumière est d'abord diffractée par un premier réseau. Un dispositif, appelé "echelle grating" est placé à un certain angle (blazing angle) du premier réseau de sorte à recevoir des ordres élevés de la première diffraction, qui sont diffractés à nouveau.

Ce dispositif permet "d'étaler" le spectre de sorte qu'une rangée de détecteurs CCD reçoit des longueurs d'ondes très proches, ce qui permet une haute résolution spectrale. En contrepartie, l'énergie est elle aussi répartie, ce qui augmente le temps d'intégration nécessaire pour recevoir suffisamment de lumière pour obtenir un certain rapport signal sur bruit.

D'une mesure à l'autre, à cause de variations internes à l'instrument (température, pression), une longueur d'onde donnée peut se décaler. Comme les mesures doivent pouvoir être comparées entre elles; ce problème doit être résolu efficacement: il faut étalonner l'instrument. Sur ce point, les instruments européens et américains diffèrent. Pour les premiers: ELODIE, CORALIE, HARPS, HARPS-N, l'étalonnage se fait en observant simultanément l'étoile cible est une source dont le spectre est connu. ELODIE observe le ciel, HARPS une lampe thorium-argon calibrée et HARPS-N utilise deux calibrations: un spectre de Fabry-Perot et un "laser frequency comb" (un laser dont le spectre est constitué de raies régulièrement espacées). Les instruments américains font passer la lumière par une cavité contenant de l'iode, dont la position des raies d'absorption est connue. Dans les deux cas on peut comparer Les raies du spectre de référence et celles de l'étoile observée. Si leur déplacement est corrélé (elles se décalent simultanément), il est dû à l'instrument.

La résolution spectrale de ces spectrographes, c'est à dire le rapport $R=\frac{\lambda}{\delta \lambda}$ d'une longeur d'onde $\lambda$ et de la sa variation détectable par le dispositi f $\delta \lambda$ est de l'ordre de 100000. Des simulations numériques (Hatzes & Cochran 1992) ont montré que l'écart type sur la mesure finale de vitesse radiale $\sigma_{rv}$ vérifie:

$\sigma_{rv} = k I^{-\frac{1}{2}} {\Lambda}^{-\frac{1}{2}} R^{-1}$

Où est l'intensité reçue, $\Lambda$ est la plage de fréquences considérées et est une constante de proportionnalité. Comme certaines longueurs d'ondes jugées contaminées peuvent être exclues de certaines mesures, cette valeur varie d'une mesure à l'autre.

Principe du spectrographe d'échelle

Le rayon incident est d'abord diffracté sur une grille standard (std. grating), puis à nouveau diffracté par le réseau d'échelle. Les trois spectres finalement obtenus sont reçus par des capteurs CCD.

Crédit : "Echelle Principle" by Boris Považay (Cardiff University) - Own work. Licensed under CC BY-SA 2.5 via Wikimedia Commons

Cross-correlation Function

Nous avons expliqué que la mesure des vitesses radiales se fait par mesure du déplacement du spectre. Cependant, nous n'avons pas précisé comment ce déplacement était mesuré. Il y a a priori beaucoup d'estimateurs classiques. La méthode majoritairement adoptée repose sur une "cross correlation function" (CCF) ou fonction de corrélation.

Les étoiles ont certains types de spectres que l'on sait reconnaître. Pour une étoile donnée, on définit un "masque" $M(\lambda)$ valant 1 pour $\lambda$ correspondant à une raie d'absorption de l'étoile et 0 ailleurs. Le spectre de l'étoile observée $I(\lambda)$ est multiplié par ce masque décalé d'une valeur $\delta \lambda$ et on mesure $CCF(\delta \lambda) = \int_{\lambda_0+\delta \lambda}^{\lambda_1+\delta \lambda} M(\lambda + \delta \lambda) I(\lambda) d \lambda$ . Où $\lambda_0$ et $\lambda_1$ désignent les bornes inférieures et supérieures du spectre observé. Ensuite, une fonction gaussienne $Ae^{\frac{(\lambda- \mu_{\lambda})^2}{\sigma^2}}$ est ajustée sur la CCF. La valeur $\mu_\lambda$ correspondant au minimum de la fonction ajustée est prise comme valeur moyenne du déplacement. L'analyse consiste ensuite à comparer les $\mu_\lambda$ issus d'observations différentes.

La CCF n'est elle même pas symétrique. Ses propriétés d'assymétrie sont analysées, car elles sont significatives d'effets physiques.

Réduction - Problème inverse

Objectif

La modélisation physique du phénomène nous a permis d'obtenir un modèle $f(\theta,T)$ où $\theta \in \mathbb{R}^n$ sont les paramètres du modèle et $T \in \mathbb{R}^m$ sont les instants d'observation. De nombreuses sources d'erreurs ont aussi été listées. On se pose maintenant la question suivante: comment estimer les paramètres du modèle compte tenu des observations ?

Après avoir listé diverses sources de signal et de bruits de l'émission de la lumière à la valeur donnée par le déteteur, nous allons maintenant donner une expression finale au modèle. Cette expression n'est pas universelle, et selon l'étoile observée, la précision recherchée, d'autres formulations peuvent être préférables.

Après avoir établi ce modèle, on donnera quelques principes d'analyse statistique et des moyens algorithmiques pour estimer les paramètres du modèle. On verra en particulier que l'analyse "dans le domaine fréquenciel" est particulièrement importante.

Modèle final

Le modèle comportera une partie déterministe et une partie aléatoire. Pour l'astrométrie, la position $x(t) = \alpha(t) \cos \delta(t)$ $y(t)) = \delta(t)$ sur la sphère céleste est

$\begin{array}{ccc} x(t) = x_0 + \mu_x t + \eta_x t^2 + \sum\limits_{k=1}^{n_p} B_j X_j(t) + G_j Y_j(t) + \varphi \Pi_x(t) + S_x(t) + \Delta x_{atm}+ \epsilon_{S_x} + \epsilon_{x,atm} + \epsilon_x\\ y(t) = y_0 + \mu_y t + \eta_y t^2 + \sum\limits_{k=1}^{n_p} A_j X_j(t) + F_j Y_j(t) + \varphi \Pi_y(t) + S_y(t) + \Delta y_{atm} + \epsilon_{S_y}+ \epsilon_{y,atm} + \epsilon_y \end{array}$

Où (x_0,y_0) est la position initiale de l'étoile, $\mu_x, \mu_y$ sont les composantes du mouvement propre $\eta_x, \eta_y$ sont des termes d'accélération de perspective, B_j, G_j, A_j, F_j sont les constantes de Thiele-Innes de la planète , X_j et Y_j sont ses coordonnées sur son plan orbital, $\varphi$ est la parallaxe, $\Pi_x, \Pi_y$ sont les coefficients parametrant le mouvement de la Terre, Les $\epsilon$ sont les bruits résiduels modélisés par des bruits gaussiens.: $\epsilon_{S_x}, \epsilon_{S_y}$ sont les bruits stellaires, $\epsilon_{x,atm}, \epsilon_{y,atm}$ sont les bruits atmosphériques et $\epsilon_x, \epsilon_y$ représentent des bruits instrumentaux.

$V(t) = \dot{\overrightarrow{OT}} + \mu_z + K\left( \cos(v+\omega) + \cos \omega \right) + S_z(t) + z_{atm}(t)+ \epsilon_{S} + \epsilon_{mes}$

Où $\dot{\overrightarrow{OT}$ est la vitesse de l'observateur dans le référentiel barycentrique du système solaire, $\mu_z$ est la composante du mouvement propre dans la direction radiale, S_z(t) est le signal stellaire dû à la granulation, aux oscillations et à l'activité, $\epsilon_S$ est le bruit stellaire résiduel et $\epsilon_{mes}$ le bruit associé à la mesure.

Dans les deux cas, les techniques de réduction de données visent à trouver des paramètres qui sont "plausibles", en l'occurrence, qui reproduisent les observations.

Traitement statistique

Rappelons qu'une variable aléatoire peut être vue comme un programme informatique qui délivre des valeurs suivant une certaine distribution de probabilité lorsqu'on lui demande. Le problème que nous posons maintenant est équivalent au suivant. Supposons qu'un ordinateur ait en mémoire des paramètres $\theta \in \mathbb{R}^p$ (nombre de planètes, leurs caractéristiques orbitales, la période de rotation de l'étoile etc.). S'il n'y avait pas de bruit, une mesure à l'instant t_k reviendrait à demander à l'ordinateur d'évaluer une fonction $f(\theta,t_k)$ . Nous connaissons t_k et (c'est l'un des modèles de la page précédente), mais nous ne connaissons pas $\theta$ . Notre but est de le déterminer à partir des mesures $y_k = f (\theta,t_k)$ . Ce principe est similaire à la résolution d'une énigme: quelqu'un connaît une information $\theta$ et nous essayons de la deviner en posant une question. Dans le cas sans bruit, celui qui pose l'énigme ne nous induit pas en erreur, mais cela ne veut pas dire que la résolution est facile !

Exemple: on veut trouver les paramètres d'une fonction affine du temps f(a,b,t) = at+c (ici $\theta = \left( \begin{array}{c} a \\c \end{array} \right)$ ). On évalue la valeur de en t_1 , on obtient y_1 = at_1+c : on a deux inconnues pour une équation, on ne peut pas résoudre. Si on a y_2 = at_2 + c , avec $t_2 \neq t_1$ , alors on a $a = \frac{y_2-y_1}{t_2-t_1}$ et c = y_1 - at_1 = y_2 - at_2 .

Notre cas est plus compliqué. Les valeurs que nous obtenons sont $y_k = f(\theta, t_k) + b_k$ où b_k est la réalisation d'une variable aléatoire $\epsilon_k$ . A l'appel numéro du programme, l'ordinateur fait appel à un autre programme qui délivre une variable alétoire selon une certaine loi. En l'occurrence, nous supposons cette loi gaussienne. Si nous faisons mesures, tout se passe comme si un programme principal évaluait la fonction $f(\theta, t_k)$ et programmes secondaires $\epsilon_k$ , k=1..m retournent chacun une valeur b_k . Si nous pouvions remonter le temps et faire les mesures plusieurs fois aux mêmes instants, on aurait des des vecteurs de mesures $y_1 = \left( \begin{array}{c} y_1 \\y_2 \\ \\y_m \end{array} \right) = f(\theta,T) + b_1$ , puis $y_2 = f(\theta,T) + b_2$ etc. (notez qu'ici b_1 et b_2 sont des vecteurs.

Comme nous ne connaissons que la loi suivie par les variables $\epsilon_k$ , ils nous est impossible de connaître les paramètres avec certitude. On leur attache une "erreur", qui quantifie l'incertitude que l'on a sur eux. En reprenant l'exemple précédent on mesure y_1 = at_1 + c +b_1 et y_2 = at_2 + c +b_2 . Si on estime et avec les mêmes formules, on fera une erreur $\frac{\sqrt{\sigma_1^2 + \sigma_2^2}}{t_1-t_2}$ sur et une erreur $\sqrt{\left(\frac{t_1}{t_2-t_1} \sigma_2\left)^2 + \right((1-\frac{t_1}{t_2-t_1}) \sigma_1 \left)^2}$ sur (admis).

Bruits

Dans le modèle $f(\theta,T) + \epsilon$ , le symbole $\epsilon$ désigne un bruit gaussien. Comme ils apparaissent constamment en détection de planètes extrasolaires et ailleurs, nous allons en donner quelques propriétés.

A une expérience donnée, $\epsilon$ prendra une valeur imprévisible. La probabilité que la valeur de soit comprise entre b_1 et b_2 est $\int_{b_1}^{b_2} f(x) dx$ où f(x) est la densité de probabilité de $\epsilon$ . Dire que $\epsilon$ est un bruit gaussien veut dire que sa densité est de la forme $f(x) = \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{(x-\mu)^2}{\sigma^2}}$ où $\mu$ et $\sigma$ sont des réels, qui sont égaux respectivement à la moyenne et à l'écart-type de $\epsilon$ . On note souvent $\epsilon \sim g(\mu, \sigma^2)$ , qui signifie " $\epsilon$ suit une loi gaussienne de moyenne $\mu$ et de variance $\sigma^2$ . Dans la plupart des cas, le bruit est de moyenne nulle (c'est le cas ici).

Dans le modèle, des bruits d'origines différentes s'additionnent. Sachant que le résidu de l'activité stellaire que nous n'avons pas ajusté $\epsilon_S$ et le bruit de mesure $\epsilon_{mes}$ suiven une certaine loi, quelle loi suivra $\epsilon = \epsilon_{S} + \epsilon_{mes}$ ? Nous pouvons déjà dire que la moyenne de $\epsilon$ sera égale à la somme des moyennes de $\epsilon_S$ et $\epsilon_{mes}$ car l'espérance est un opérateur linéaire. Peut-on dire plus ? Si ces bruits dépendaient l'un de l'autre, la réponse pourrait être complexe. En l'occurrence, la physique de l'étoile cible et les erreurs instrumentales sont totalement indépendantes. On peut montrer que dans ces conditions, la variance de $\epsilon$ est égale à la somme des variances de $\epsilon_S$ et $\epsilon_{mes}$ . Nous pouvons même aller plus loin car la somme de deux variables gaussiennes indépendante est une variable gaussienne. En résumé, $\epsilon \sim g\left(\mu_{S} +\mu_{mes}, \sigma_{S}^2+ \sigma_{mes}^2\right)$ en l'occurrence $\mu_S$ et $\mu_{mes}$ sont nulles.

Lorsqu'on dispose de plusieurs mesures, à l'expérience numéro on a un certain bruit b_k réalisation d'une variable $\epsilon_k$ de densité f_k . La plupart du temps, on fait l'hypothèse que les brutis $\epsilon_k$ sont indépendants, c'est à dire que la probabilité d'obtenir le bruit b_k à l'expérience ne dépend pas des valeurs prises aux expériences précédentes et suivantes. Lorsque ce n'est pas le cas on parle de bruits corrélés. Pour les caractériser, on utilise souvent leur densité spectrale de puissance. Un certain profil de densité spectrale correspond à une "couleur" du bruit.

A retenir: la somme de variables gaussienne indépendantes $\epsilon = \sum\limits_{k=1}^n \epsilon_k$ où $\epsilon_k \sim g\left( \mu_k, \sigma_k^2)$ est une variable gaussienne suivant la loi $g \left( \sum\limits_{k=1}^n \mu_k , \frac{1}{n} \sum\limits_{k=1}^n \sigma_k^2 \right)$ .

Signification statistique

Il esiste plusieurs outils pour s'assurer qu'une détection possible n'est pas due au bruit. L'un des plus utilisé est le test de signification (significance en anglais), qui consiste à calculer la probabilité d'avoir le signal observé "au moins aussi grand" s'il n'y avait en réalité que du bruit. Par exemple, qupposons que l'on veuille mesurer une quantité qui est perturbée par un bruit gaussien additif de moyenne nulle et d'écart-type $\sigma=1$ , donnant une mesure y = a+b .

Si il n'y avait en réalité pas de signal ( a=0 ), les mesures seraient uniquement dues au bruit. On imagine deux cas de figures:

1er cas: on mesure . La probabilité d'avoir obtenue une valeur au moins aussi grande par hasard, c'est à dire la probabilité que $|y| \geq 0.1$ sachant que est $\int_{-\infty}^{-0.1} e^{-x^2}dx + \int_{0.1}^{\infty} e^{-x^2}dx = 2 \int_{0.1}^{\infty} e^{-x^2}dx \approx 0.818$
2ème cas: on mesure . La probabilité pour que alors que est $\int_{-\infty}^{-10} e^{-x^2}dx + \int_{10}^{\infty} e^{-x^2}dx = 2 \int_{10}^{\infty} e^{-x^2}dx \approx 2.06 \times 10^{-9}$

Ces valeurs ont une interprétation: si on réalisait exactement le même type de mesure alors qu'il n'y a pas de signal, on observerait $|y| \geq 0.1$ dans 81.8 % des cas et $|y| \geq 10$ dans $2.06 \times 10^{-7}$ % des cas.Dans le premier cas, la probabilité d'avoir un signal aussi grand que celui que l'on a mesuré est grande. On ne peut pas assurer qu'un signal a été détecté. Par contre, dans le deuxième cas on serait dans un des deux cas sur un milliard où le signal serait dû au bruit. On peut alors dire qu'on a détecté un signal à $10 \;\sigma$ , car la valeur de l'écart-type du bruit est , sa moyenne est , donc on a $\frac{y}{\sigma-0} = \frac{10}{1} = 10$ .

Détecter un signal "à 10 sigmas" est un luxe que l'on peut rarement se payer. Les détections sont annoncées plutôt pour des valeurs de sigmas.

Remarque importante: on calcule la probabilité d'avoir les observations sachant qu'il n'y a pas de signal et non la probabilité d'avoir un signal sachant les observations qui est une quantité qui a davantage de sens. Le calcul de cette dernière quantité se fait dans le cadre du calcul bayésien, outil très puissant qui ne sera pas développé dans ce cours.

Exercices

Auteur: Nathan Hara

Espérance et variance de la moyenne empirique

Difficulté : ☆☆☆

Question 1)

On rappelle que l'espérance et la variance d'une variable alétoire de densité de probabilité sont données par $\mathbb{E}(X) = \int_{-\infty}^{+\infty} xf(x)dx$ et $\text{Var}(X) = \mathbb{E}\left((X-\mathbb{E}(X))^2\right)= \int_{-\infty}^{+\infty} (x-\mathbb{E}(X))^2f(x)dx$

Soit une variable aléatoire et un réel. Montrer que $\mathbb{E}(aX) = a \mathbb{E}(X)$ et $\text{Var}(aX) = a^2 \text{Var}(X)$
Soient et deux variables aléatoires indépendantes. Montrer que $\mathbb{E}(X+Y) = \mathbb{E}(X) + \mathbb{E}(Y)$ .
On définit la covariance par $\text{Cov}(X,Y) = \mathbb{E}\{ (X-\mathbb{E}(X)) (Y-\mathbb{E}(Y))\}$ . Montrer que $\text{Var}(X+Y) = \text{Var}(X) + \text{Var}(Y) + 2 \text{Cov}(X,Y)$
Nous allons maintenant voir un cas où la précision d'un estimateur est facilement calculable. On considère mesures entachées de bruits gaussiens d'une quantité fixe . Plus précisément, la mesure numéro est modélisée par une variable aléatoire $X_j = K + \epsilon_j$ où $\epsilon_j$ est un bruit gaussien de moyenne nulle et de variance $\sigma^2$ . On suppose que les mesures sont indépendantes, ce qui implique en particulier que $\text{Cov}(\epsilon_j,\epsilon_i) = 0$ pour $i \neq j$ . Pour obtenir une estimation de , on fait la moyenne empirique des expériences, c'est à dire $M = \frac{X_1 + X_2 +... X_m}{m}$ . Montrer que $\mathbb{E}(M) = K$ et $\text{Var}(M) = \frac{\sigma^2}{m}$
La précision de l'estimateur augmente-t-elle avec le nombre de mesures ?

Loi du chi 2

Nous nous sommes toujours ramenés à des modèles du type: modèle déterministe + bruit gaussien. Afin de vérifier que les observations sont compatibles avec le modèle, on étudie les résidus, définis comme "les observation - le modèle ajusté".

La loi du $\chi^2$ est un outil commode pour étudier le comportement de plusieurs variables gaussiennes. Considérons d'abord une famille de variables aléatoires gaussiennes indépendantes $(\epsilon_k)_{k=1..m}$ , de moyenne nulle et de variance unité. On forme la quantité $\chi^2_m= \epsilon_1^2 + \epsilon_2^2 + ... + \epsilon_m^2$ . Comme les $\epsilon_k$ sont des variables aléatoires, les $\epsilon_k^2$ le sont aussi. La somme de variables alétoires étant toujours une variable alétoire, $\chi^2_m$ suit une certaine loi de probabilité. Dans l'analogie avec un programme informatique, la variable alétoire $\chi^2_m$ se comporte comme un programme qui appelle programmes générant une variables gaussiennes, puis additionne leurs carrés. Elle est appelée loi du $\chi^2$ à m degrés de liberté. On peut montrer qu'en moyenne une variable gaussienne au carré $\epsilon^2$ a une moyenne de 1. En conséquence, $\chi^2_m$ vaudra typiquement .

Pourquoi cette loi serait utile pour notre cas ? Si le modèle est bien ajusté, les résidus doivent se comporter comme un bruit gaussien. En supposant que les erreurs sont toutes indépendantes, de moyenne nulle et de variance unité, les résidus $r_k = y-f(\theta,t_k})$ en sont une réalisation. Donc $R^2=\sum\limits_{k=1}^m r_k^2$ est une réalisation d'une loi du $\chi^2$ à degrés de liberté. Si R^2 est de l'ordre de , le modèle est cohérent. Sinon, le modèle ou les paramètres ajustés sont à revoir.

En pratique, les erreurs $\epsilon__k$ ne sont évidemment pas de variance unité et parfois pas indépendantes. Par contre on peut à bon droit supposer qu'elles sont de moyenne nulle. Pour se ramener au cas précédent, on calcule non pas une réalisation de $\epsilon_1^2+...+\epsilon_m^2$ mais de $\chi^2_m'=\epsilon^T \Sigma^{-1} \epsilon$ où $\epsilon = \left( \begin{array}{c} \epsilon_1 \\ \epsilon_2 \\ \\ \epsilon_m \end{array} \right)$ , $\epsilon^T$ est sa transposée et $\Sigma$ est la matrice des variances-covariances de $\epsilon$ . Dans le cas où les $\epsilon_k$ sont indépendantes, la matrice des variances-covariances est diagonale, son -ème terme diagonal étant $\frac{1}{\sigma_k^2}$ , soit l'inverse de la variance de $\epsilon_k$ .

La vraisemblance

Etant donné des paramètres $\theta$ , le modèle global $f(\theta,T)+\epsilon$ est une variable aléatoire: il est somme d'une variable aléatoire valant $f(\theta,T)$ avec une probabilité 1 et d'un vecteur de variables aléatoires gaussienne $\epsilon$ . A ce titre, il a une certaine densité de probabilité que l'on note $L(y|\theta)$ . Le symbole | se lisant "sachant". La lettre L vient de Likelihood, qui veut dire vraisemblance en anglais. Il s'agit dans l'idée de la probabilité d'obtenir $y = f(\theta,T)+\epsilon$ pour une valeur de $\theta$ donnée..

La fonction $\theta \rightarrow L(y|\theta)$ est souvent appelée "fonction de vraisemblance". La valeur de $\theta$ maximisant $L(y|\theta)$ est appelé l'estimateur du maximum de vraisemblance. Il a de bonnes propritétés statistiques. En effet, on peut montrer que c'est un estimateur:

Non biaisé: lorsque le nombre de mesure tend vers l'infini, l'estimateur du maximum de vraisemblance tend vers la bonne valeur $\theta$
De variance minimale: parmi les estimateurs non biaisés, il a la variance la plus faible que l'on puisse espérer.

Dans notre cas, si les $\epsilon_k$ sont des variables indépendantes, leur densité de probabilité jointe est égale au produit de leurs densité de probabilité. $L(y|\theta) = \prod\limits_{k=1}^m g_k(y-f(\theta,t_k))$ où g_k est la densité de probabilité de la variable $\epsilon_k$ . De plus, si ces varibles sont gaussiennes et indépendantes, on a: $L(y|\theta) = \prod\limits_{k=1}^m \frac{1}{\sqrt{2\pi}\sigma_k} e^{-\frac{(y-f(\theta,t_k))^2}{\sigma_k^2}} = \frac{1}{(\sqrt{2\pi})^m \prod\limits_{k=1}^m \sigma_k} e^{-\sum\limits_{k=1}^m \frac{(y-f(\theta,t_k))^2}{\sigma_k^2}}}$

Méthode des moindres carrés

Pour des bruits gaussiens indépendants, maximiser la vraisemblance, revient à minimiser $\theta \rightarrow F(y,\theta)=-ln(L(y|\theta)) =\sum\limits_{k=1}^m \frac{(y-f(\theta,t_k))^2}{\sigma_k^2}}$ puisque $x \rightarrow e^{-x}$ est une fonction décroissante. La méthode consistant à minimiser $F(y,\theta)$ s'appelle la méthode des moindres carrés. C'est la méthode d'estimation de loin la plus utilisée dans tous les domaines. Elle est parfois utilisée quand les bruits ne sont pas gaussiens, mais il faut garder à l'esprit qu'elle n'a alors plus de propriétés statistiques sympathiques (sauf quand le modèle est linéaire en $\theta)$ .

Dans notre cas, les paramètres $\theta$ sont les éléments des orbites, les paramètres du bruit stellaire, du mouvement propre, etc. La fonction $F(y,\theta)$ a donc de nombreux paramètres, et trouver son minimum global est une tâche ardue qui fait l'objet d'une littérature très vaste.

Lorsque le modèle est linéaire en $\theta$ i. e. $f(T,\theta) = A \theta$ où est une certaine matrice dépendant des instants d'observation $T = (t_k)_{k=1..m}$ , l'ajustement est beaucoup plus simple car il a une solution explicite (voir mini-projet). On essaye de se ramener autant que possible à des ajustements linéaires. La plupart du temps, on estime les paramètre les uns après les autres, puis un ajustement global est réalisé. Une démarche classique consiste à:

Ajuster le mouvement propre et la parallaxe, puis soustraire ces ajustements du signal initial
Chercher les périodes éventuelles dans le signal à l'aide du périogogramme (voir page suivante)
Si on trouve une période, on ajuste une orbite, que l'on soustrait du signal
On recommence les étapes 2 et 3 jusqu'à ce que le périodogramme ne présente plus de pic significatif

Sur la figure, on représente les étapes d'un ajustement d'un signal astrométrique simulé (2 planètes, 45 observations). De gauche à droite et de haut en bas:

Signal brut en rouge En vert signal après correction des bruits instrumentaux et atmosphériques . Le mouvement propre (tendance linéaire) et la parallaxe (oscillation) dominent très nettement le signal.
En vert: le signal planétaire réel, sans aucun bruit $(s^{\star})$ . En rouge: signal estimé en soustrayant le mouvement propre, l'accélération de perspective et la parallaxe au signal .
En rouge: signal . En vert : orbite d'une planète ajustée sur .
En rouge : signal dont on a soustrait . En vert :orbite d'une planète ajustée sur .
En rouge: signal . En vert: , ajustement de deux planètes en partant des estimations et .
En vert, signal $(s^{\star})$ . On ajuste tous les paramètres simultanément sur en partant des estimations précédentes. Le signal en rouge représente le signal planétaire de cet ajustement.

Exemple d'ajustement en astrométrie (données simulées)

Six étapes d'ajustement aux données. Pour représenter la chronologie des mesures, on trace des segments reliant deux points de mesure consécutifs.

Périodogramme

Evaluer les résidus sur une grille d'un modèle à paramètres, où chacun d'eux peut prendre valeurs recquiert n^p évaluations, ce qui devient rapidement ingérable numériquement. Les planètes ont un mouvement périodique, donc il est raisonnable de checher des signaux périodiques dans le signal en ne faisant varier que la période du signal recherché. Pour des signaux échantillonnés à intervalles réguliers, on utilise la transformée de Fourier. Le périodogramme est un moyen de checher des signaux périodiques dans des données échantillonnées irrégulièrement. On les notera . Le périodogramme de Lomb-Scargle d'un signal $(y_k)_{k=1..m}$ échantillonné aux instants $(t_j)_{j=1..m}$ est défini comme suit pour une fréquence quelconque $\omega$ :

$P_y(\omega) = \frac{1}{2} \left( \frac{ \left( \sum\limits_{j=1}^m y_j \cos \omega(t_j-\tau) \right)^2}{ \sum\limits_{j=1}^m \cos^2 \omega(t_j-\tau)}} + \frac{ \left( \sum\limits_{j=1}^m y_j \sin \omega(t_j-\tau) \right)^2}{ \sum\limits_{j=1}^m \sin^2 \omega(t_j-\tau)}} \right)$ où $\tau$ vérifie:

$\tan 2 \omega \tau = \frac{\sum\limits_{j=1}^m \sin 2\omega t_j}{\sum\limits_{j=1}^m \cos 2\omega t_j}$

Cette expression est équivalente à $P_y(\omega) = A^2 +B^2$ où et sont les paramètres minimisant $\sum\limits_{j=1}^m (y_j-A\cos(\omega t_j) - B\sin(\omega t_j))^2$ . Le modèle $A,B, t \rightarrow A \cos \omega t + B \sin \omega t$ est linéaire en A, B , on a donc une solution explicite à la minimisation.

Le périodogramme a une propriété très intéressante: si le signal d'entrée est un bruit gaussien $\epsilon$ de variance unité, p_0 une valeur réelle fixée et $\omega$ une fréquence quelconque, la probabilité que $P_{\epsilon}(\omega)$ dépasse p_0 est $Pr\{P_\epsilon(\omega)>p_0\} = e^{-p_0}$ . En d'autres termes, la probabilité qu'une valeur du périodogramme à $\omega$ fixée soit "au moins aussi grand que p_0 " par hasard décroît exponentiellement. Supposons que l'on ait un signal y(t) = x(t) + b(t) où b(t) est un bruit gaussien de variance unité et nous trouvons un pic de taille p_0 , on calcule la probabilité de trouver un pic au moins aussi grand si le signal n'est composé que de bruit: $e^{-p_0}$ . Si cette valeur est petite, on pourra confirmer la détection d'un signal avec une erreur de fausse alarme de $e^{-p_0}$ . Ce procédé n'est autre qu'un test de signification statistique.

La figure montre un exemple de périodogramme. Il s'agit d'un périodogramme d'une des coordonnées d'un signal astrométrique simulé dont on a soustrait le mouvement propre et la parallaxe. En bleu, on représente un périodogramme idéal, sans bruit, avec 10000 observations. Le périodogramme représenté en rouge est lui calculé pour 45 observations. Le pic le plus haut correspond bien à une fréquence réelle. Par contre, le deuxième pic le plus important (à 0.37 rad/s) ne correspond pas à une sinusoïde. C'est ce qu'on appelle un alias de la fréquence principale.

En pratique, la variance du bruit n'est pas unitaire et dépend de l'instant de mesure. On peut corriger ce problème en minimisant un critère pondéré $\sum\limits_{j=1}^m \frac{(y_j-A\cos(\omega t_j) - B\sin(\omega t_j))^2}{\sigma_j^2}$ où $\sigma_j$ est la variance du bruit à la mesure .

Exemple de périodogramme

En rouge: périodogramme d'un système simulé à deux planètes après soustraction du mouvement propre et de la parallaxe (signal très peu bruité), pour 45 mesures. En bleu: périodogramme du même système sans bruit et avec 10000 observations.

Comprendre

Modélisation - Problème direct

Présentation du modèle

Problème à deux corps

Problème à deux corps newtonien

Equation de Newton

Quantités conservées

Géométrie de l'orbite

Loi des aires

Loi des Aires

Exercices

Conservation de l'énergie

Equation de Kepler

Equation de Kepler

Solutions en fonction du temps

Changement de référentiel

Paramètres orbitaux classiques

Conclusion pour le problème à deux corps

Exercices

Bonnes et mauvaises configurations d'observation

Qu'est-ce que les vitesses radiales permettent de mesurer ?

Modèle de la trajectoire

Trajectoire observée depuis le barycentre du système solaire

Trajectoire observée depuis la Terre

Vitesses radiales

Parallaxe (astrométrie)

Accélération de perspective et changement de parallaxe (astrométrie)

Accélération de perspective

Variation de parallaxe

Etoile

Etoile

Masse de l'étoile

Processus stochastique et Densité spectrale de puissance

Processus stochastique

Intégrale en moyenne quadratique (Mean Square Integral)

Processus stationnaires au sens large

Granulation

Activité magnétique

Oscillations radiales (vitesses radiales)

Perturbations atmosphériques

Vitesses radiales

Astrométrie

Récapitulation

Instrumentation & Observations

Objectifs

Télescope

CCD

Bruit de photon

Exercices

Magnitude et temps d'observation

Effet Doppler

Spectrographe

Cross-correlation Function

Réduction - Problème inverse

Objectif

Modèle final

Traitement statistique

Bruits

Signification statistique

Exercices

Espérance et variance de la moyenne empirique

Loi du chi 2

La vraisemblance

Méthode des moindres carrés

Périodogramme

Réponses aux exercices

Exercice 'Bonnes et mauvaises configurations d'observation'

Exercice 'Qu'est-ce que les vitesses radiales permettent de mesurer ?'

Exercice 'Magnitude et temps d'observation'