Ce projet a pour but de mettre en application tout votre nouveau savoir théorique pour l'appliquer à un cas pratique. Nous allons prendre des données réelles de transit, chercher si une planète s'y cache, et le cas échéant, déduire les paramètres de la planète associée, tels que son rayon et sa période. Vous comparerez alors vos résultats à ceux qui ont vraiment étaient publiés par l'équipe de recherche qui a traité ces données il y a quelques années.
Nous allons analyser des données provenant de SuperWASP qui est composé de 8 télescopes de 11cm chacun, ce qui permet de couvrir presque 500 degrés carrés sur le ciel (i.e. la surface d'environ 2500 pleines lunes à chaque prise). Il y a en fait un SuperWASP au Nord (La Palma) et un au Sud (Afrique du Sud) pour pouvoir couvrir tout le ciel (plus de détails ici ).
Voici les données historiques qui ont mené à la première détection de la planète autour de l'étoile WASP-80 en 2013. Dans ce fichier, vous trouverez les mesures photométriques de l'étoile WASP-80 (située dans la constellation de l'Aigle) à différentes dates. Les 3 colonnes indiquent respectivement la date du point photométrique (en BJD pour Barycentric Julian Date, qui donne le jour julien correspondant à l'observation), la magnitude de l'étoile mesurée à cette date, et l'erreur associée sur la magnitude.
Nous allons utiliser python pour ce mini-projet. Vous pouvez aussi utiliser votre langage de prédilection si Python ressemble à des hiéroglyphes pour vous.
Commencez par charger le fichier wasp80data.txt qui contient les observations de WASP-80 et faites un graphe représentant la magnitude en fonction du temps, avec les barres d'erreurs associées à chaque point. En python, cela donne :
import numpy as np
import matplotlib.pyplot as plt
bjd, mag, err = np.loadtxt("wasp80data.txt", delimiter=' ', skiprows=1, unpack=True)
plt.errorbar(bjd,mag,err)
On pourrait travailler avec les magnitudes, mais il est plus commun de travailler directement en flux. En effet, un transit produit une baisse du flux total qui provient de l'étoile alors que ça provoque une augmentation de la magnitude et donc on ne verrait pas une baisse mais une croissance de la magnitude lors d'un transit en magnitude. Vous pourrez essayer plus loin de suivre la même procédure en travaillant sur les magnitudes (pour constater que l'on arrive aux mêmes résultats) mais pour le moment on travaillera en flux.
Reproduire le même graphique que précédemment en unité de flux
Que voyons-nous ? Est-ce ce à quoi l'on s'attend ? Voyons-nous un transit ? Pas vraiment ! Ça ne ressemble pas exactement aux courbes théoriques que l'on a vues dans le cours avec une belle chute de flux quand la planète passe devant son étoile. Pour cela il va falloir travailler un peu les données.
Un œil averti se rendra compte que ces données on été pré-traitées. On voit que la magnitude oscille autour de zéro, la composante de magnitude due à l'étoile a donc été retirée. On voit aussi que les données semblent suivre une ligne horizontale, sans qu'il y ait de variations linéaires ou polynomiales. En réalité, quand on obtient les données brutes, c'est beaucoup moins propre que cela, et il faut retirer les effets qui ne sont pas dûs au transit mais plutôt à l'observation en tant que telle. Il peut y avoir des déviations du signal dues à l'instrumentation (par exemple à cause de variations thermiques), ou même des variations de l’atmosphère au cours de l'observation. Il y a différentes méthodes pour réajuster les données afin de corriger ces variations qui ne sont pas dues au transit. De manière générale et simplifiée, cela revient à trouver une fonction polynomiale qui suit au mieux les variations à long terme et de retirer cette fonction du signal. On appelle cela le detrending. Les données que l'on traite ici ont déjà subi ce detrending et peuvent être maintenant exploitées pour chercher un transit potentiel et trouver les paramètres de la planète le cas échéant.
Que faut-il faire pour essayer de voir un transit dans ces données ?
Ici, on a la magnitude (ou le flux) du système en fonction du temps. Pour pouvoir voir un transit, et c'est là toute la puissance de cette méthode de détection, il va falloir sommer les différents transits potentiels qu'il y a dans ces données. Pour cela, il faut trouver la bonne période sur laquelle sommer les données, et peut-être qu'alors avec une somme de plusieurs transits, on verra effectivement une courbe de lumière typique avec une baisse de flux quand la planète passe devant son étoile.
On peut essayer par exemple de sommer les signaux sur une période de 2 jours. En python cela donne :
from PyAstronomy.pyasl import foldAt
best_period=2#jours
phases = foldAt(bjd, best_period)
sortIndi = np.argsort(phases)
phases = phases[sortIndi]*best_period
mag = mag[sortIndi]
plt.plot(phases, mag)
On ne voit toujours rien. Idéalement, il faudrait pouvoir faire cela pour un très grand nombre de périodes et visualiser le résultat jusqu'à l'obtention potentielle d'une baisse de flux à un endroit de la courbe (i.e. voir le transit). C'est la méthode que l'on va utiliser, mais bien sûr, il est impossible de visualiser des millions de courbes à l'oeil et il faut trouver un bon indicateur que l'ordinateur puisse tester lui-même pour nous dire quand il y a effectivement une baisse de flux visible et donc un transit.
La méthode la plus couramment utilisée dans le domaine de la recherche est la méthode BLS (pour Box Least Square, développée par Kovacs et al. 2002). L'idée est de dire qu'un transit peut être en première approximation modélisé par une fonction porte (fonction rectangulaire). Un transit est en effet bien plus proche d'une fonction porte que d'un sinus ou cosinus car les variations lors de l'entrée ou de la sortie du transit qui peuvent être non rectangulaires ne représentent qu'une petite partie du signal et le reste du transit est plutôt plat. La décomposition en série de Fourier ne fournirait pas une fréquence donnée dominante mais plutôt un ensemble de fréquences alors que la fonction porte ne requiert qu'un terme pour bien modéliser le signal, d'où l'avantage.
On crée alors un modèle simple de fonction porte avec trois paramètres : sa profondeur, sa longueur (durée du transit) et sa phase (la date où le transit a lieu). On fait ensuite tourner le modèle sur un très grand nombre de périodes. Pour chaque période testée, l'algorithme BLS essaye de trouver les meilleurs paramètres (profondeur, longueur, phase) du modèle pour expliquer les observations. L'algorithme compare le meilleur modèle aux données et estime alors la vraisemblance (en log) du modèle en question (par exemple par une méthode des moindres carrés) et crée un périodogramme qui donne la vraisemblance du modèle en fonction de la période. On peut alors voir si certaines périodes paraissent plus vraisemblables que d'autres. Les pics dans le périodogramme indiquent la présence de planètes en transit ou du bruit non pris en compte dans le pré-traitement (detrending). Il faut faire attention aux faux-positifs avec cette méthode et les repérer est un champ actif de recherche dans lequel nous ne rentrerons pas en détails.
Trouvez la période du transit et sa profondeur dans les données de WASP-80.
Vérifiez ces valeurs en sommant les données initiales sur la période trouvée et en affichant le flux en fonction de la phase (vous pouvez utiliser le code de la question 1 ci-dessus pour sommer les différents transits sur une période donnée). Vous produirez aussi une deuxième courbe où vous moyennerez le signal dans des bins temporels d'environ 15 minutes pour avoir une courbe de transit plus lisse.
En faisant ce travail , on trouve le diagramme de la figure ci-jointe. :
Maintenant, on reconnaît une courbe de transit typique. On peut même zoomer sur la phase de transit pour obtenir la figure ci-jointe.
On voit alors que la durée du transit totale est inférieure à 0.1 jours, i.e. <2.4h. Pour obtenir la courbe rouge rebinnée en python, on peut écrire un code ad hoc ou procéder ainsi :
from scipy.stats import binned_statistic
bin_means = binned_statistic(phases, mag, bins=300)
bin_means2 = binned_statistic(phases, phases, bins=300)
plt.plot(bin_means2[0], bin_means[0],color='red',linewidth=5)
Déduisez maintenant le rayon de la planète et son demi-grand axe.
Il faudra d'abord trouver le rayon de l'étoile WASP-80 afin de pouvoir déterminer le rayon de la planète.
Le calcul du rayon de l'étoile que l'on vient de mener n'est qu'un calcul approximatif etTriaud et al. (2013) sont capable d'estimer ce rayon plus précisément à partir de leur jeu de données, ils trouvent , que l'on utilisera par la suite pour comparer à leurs résultats.
Maintenant on peut déduire le rayon de la planète. En supposant que l'orbite est circulaire, on peut aussi déduire facilement son rayon orbitale à partir de la 3ème loi de Kepler
On peut maintenant comparer aux valeurs trouvées dans le papier originel de Triaud et al. (2013) (tableau ci-joint).
On voit que nos valeurs pour , la période et le demi-grand axe sont très proches des valeurs originelles.
On peut aussi comparer les valeurs trouvées avec celles de l'encyclopédie exoplanet.eu.
Comme vu dans le cours, on peut faire une première vérification pour tester pour un faux positif, en comparant la durée du transit théorique avec sa durée réelle que l'on a obtenue à la question 3 de la page précédente.
Pour aller plus loin vous pouvez aussi essayer d'obtenir plus d'informations à partir des données, par exemple en utilisant, au lieu de la fonction porte, un vrai modèle de transit à la Mandel & Agol (2002) pour pouvoir déduire le paramètre d'impact et une possible excentricité en prenant en compte les effets d'assombrissement centre au bord (il y a une fonction python qui implémente ce modèle ).
Ce dossier(une fois décompressé, le fichier readme.txt décrit les différents fichiers) donne un ensemble de données qui permettent de faire d'autres analyses de cette planète:
Le dossier donne aussi des données Spitzer (infrarouge) du même système qui ont été utilisées dans Triaud et al. (2015). Dans ces données, on peut voir l'éclipse secondaire (en plus du transit) car le contraste dans l'infrarouge est plus favorable. Attention, il va falloir faire du detrending pour bien corriger les données brutes. Vous pouvez alors déduire la température de brillance de la planète. On s’attend à ce que la température d’équilibre et de brillance soient du même ordre, mais pas exactement les mêmes. On s'attend à ce que le côté jour de la planète soit légèrement plus chaud et le côté nuit légèrement en dessous de la température d’équilibre.
Pour avoir un jeu de données complet, on fournit aussi les observations en vélocimétrie radiale avec Coralie et HARPS pour ce même système afin de pouvoir aussi déterminer la masse de la planète et ainsi connaître sa densité (et donc pouvoir inférer sa composition).
Vous pouvez faire une courbe de la vitesse radiale en fonction de la phase (comme pour la courbe de transit) et voir ce qu'il se passe, en particulier pendant le transit. On peut, en effet, voir l'effet de Rossiter-McLaughlin dans les données HARPS, ce qui permet de contraindre le sens de rotation de l'orbite de la planète et l’orientation de l’orbite planétaire, à savoir l’angle entre le plan orbital et l’axe de rotation de l’étoile.
Au boulot !
pages_ind-transits/wasp89-donnees.html
La conversion d'un Flux F à une magnitude m est telle que
pages_ind-transits/wasp80-recherche.html
Il faut utiliser la méthode BLS décrite précédemment. Soit vous la codez vous-même, soit vous utilisez une fonction pré-existente qui intègre déjà cet algorithme BLS. Il y a de telles fonctions en python ou fortran par exemple.
Pour gagner du temps nous utiliserons la fonction python bls.eebls qui permet de faire tourner l'algorithme BLS en une ligne de code (il faut suivre l'installation du module ici ). On peut aussi essayer de réinventer la roue et de le développer soi-même pour comprendre tous les tenants et aboutissants de cette fonction.
Faisons appel à la fonction BLS python pour trouver notre transit. On appellera la fonction bls.eebls où l'on doit fournir les arguments qui suivent en entrée:
bls.eebls(bjd, flux, u, v, nf, fmin, df, nb, qmi, qma)
bjd est la date. Flux est le flux ou la magnitude. Puis les autres arguments sont fixés comme ceci :
On appelle ensuite la fonction avec u et v qui sont vides et servent à des calculs intermédiaires et on stocke le résultat dans la variable results :
results = bls.eebls(bjd, mag, u, v, nf, fmin, df, nb, qmi, qma)
On stocke les résultats dans différentes variables décrites ci-dessous
power, best_period, best_power, depth, q, in1, in2 = results
power est le spectre de puissance de dimension nf aux fréquences f = fmin + arange(nf) * df,
best_period est la meilleure période de transit trouvée par l'algorithme (dans la même unité que le temps mis en entrée),
best_power est la puissance (du spectre de puissance) pour cette meilleure période trouvée,
depth est la profondeur du transit pour la meilleure période,
q est la durée du transit par rapport à la période,
in1 est l'index du bin où le transit commence,
in2 est l'index du bin où le transit se termine.
Puis on crée le périodogramme résultant comme ceci : plt.plot(np.arange(fmin, fmin+nf*df, df), power)
Cela retourne le périodogramme de la figure ci-jointe.
On constate donc qu'il y a une fréquence orbitale privilégiée à environ 0.33 orbite par jour, i.e. 1/0.326=3.067 jours par orbite. L'algorithme retourne une période (best_period) de 3.0674846 jours et une profondeur de transit (depth) de 0.0313459 pour cette meilleure période.
pages_ind-transits/wasp80-planete.html
Nous allons estimer le rayon de l'étoile à partir de sa température et de sa distance. La température est estimée à partir de ses magnitudes dans la couleur bleue (B) et visible (V) grâce au graphe ci-joint.
Les données de l'étoile, sa distance et ses magnitudes B et V sont données dans le tableau venant de l'article de Triaud et al. ( ci-joint). Ces données sont aussi accessibles dans la base de données Simbad
Notez que Simbad donne la parallaxe de l'étoile en milliarcseconde d'arc. La distance exprimée en parsec est l'inverse de la parallaxe exprimée en seconde d'arc.
On utilise le fait que la luminosité d'une étoile est , que
, où
est la magnitude absolue de l'étoile
(i.e. la magnitude qu'aurait cette étoile si elle était à 10pc) et enfin que la magnitude absolue
et la magnitude visuelle,
sont liées à la distance
par
La magnitude absolue du Soleil est de 4.83 et sa température est de 5800K.
Dans ces calculs, on utilise différentes magnitudes dont les définitions sont données ici. La magnitude visuelle est notée ou
.
La différence de magnitude en bande B et en bande V égale à B-V=0.93. On trouve une température d'environ 4000K. On a affaire à une étoile de type K, plus petite que notre soleil (de type G2).
La distance est .
et donc
où l'index
est pour le Soleil et
pour l'étoile WASP-80. D'après la définition de la magnitude, on a aussi que
, où
est la magnitude absolue de l'étoile. Pour transformer la magnitude mesurée
(qui dépend de la distance de l'étoile) en magnitude absolue, on utilise
, soit
et donc on obtient
, soit
.
On utilise Tpla= 3.0674846 jours et la masse de l'étoile est dans Triaud et al. (2013), soit 0.57Ms.
Comme on connait la profondeur du transit de 0.0313459 et en utilisant la relation du cours sur taille de la planète, la profondeur = , on déduit que
.
Le demi-grand axe .
On aurait aussi pu deriver la masse en utilisant les relations entre masse et rayon ou entre masse et luminosité.
Comme on a vu dans le cours, la durée du transit est égale à : en supposant un paramètre d'impact nul. On obtient donc une durée de 0.075 jours ou 1.8 heures, ce qui est proche de ce que l'on voit sur le transit zoomé de la page précédente.