WASP-80 : Les transits

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.

AideSolution

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 :

• u=np.empty(bjd.shape)
• v=np.empty(bjd.shape)
• nf=1000#nombre de bins de fréquences à tester
• fmin=1./(100.)#fréquence minimum à tester (100 jours)
• df=0.0005#taille des bins de fréquences
• nb=1000#nombre de bins à utiliser pour la courbe de lumière sommée
• qmi=0.001#durée de transit minimum à tester
• qma=0.2#durée de transit maximum à tester

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.

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)