Acessibilidade / Reportar erro

O Escurecimento de limbo via modelagem de trânsitos e inferência bayesiana

The limb darkening via transit modeling and bayesian inference

Resumos

Missões espaciais de busca por exoplanetas exigem um tratamento do escurecimento de limbo com maior acurácia e eficiência. Neste contexto, os coeficientes de escurecimento de limbo produzem efeitos sutis nas curvas de luz e por isso são frequentemente degenerados com outros parâmetros. Nesse trabalho, estudamos a estrela KIC 5357901 (Kepler-425) de tipo espectral K0, com temperatura efetiva de 5170K, a partir da análise da curva de luz obtida pelo satélite Kepler da NASA. A escolha dessa estrela foi baseada na multiplicidade de trânsitos. Produzimos uma modelagem do trânsito do exoplaneta Kepler-425b e aplicamos métodos de inferência Bayesiana para obter um intervalo de credibilidade baseado em um algoritmo de Monte Carlo via cadeia de Markov (MCMC). Obtivemos para os coeficientes de escurecimento de limbo da estrela Kepler-425 (KIC 5357901) os valores u1=0,559-0.070+0.071 e u2=0,104-0.117+0.120 pela lei quadrática. Esta lei quadrática fornece maior alcance e melhor ajuste de dados em comparação com as outras leis. Nossas determinações são consistentes com aquelas obtidas anteriormente com diferente parametrização. Nosso resultado com um número menor de iterações por cadeia apresenta acurácia similar, e da ordem de 10-2 para u1, além explorar de forma mais focada a região de interesse dos parâmetros.

Palavras chave:
Exoplanetas; Escurecimento de Limbo; Método do Trânsito Planetário; Fotometria


Space missions to search for exoplanets by photometry have been improved in terms of precision and require a more accurate and efficient treatment of limb darkening coefficients, and this seems to be a bottleneck in modeling the transits of discovered planets. Limb darkening coefficients often have subtle effects on light curves and are therefore often degenerate with other modeling parameters. In this work, we study the K0-type star KIC 5357901 (Kepler-425), with effective temperature of 5170K, from light curve obtained by NASA’s Kepler satellite. We drove photometric modeling of the Kepler-425b exoplanet and we applied a Bayesian inference method to obtain a credible interval based on a Monte Carlo Markov chain (MCMC) algorithm. We obtained the Kepler-425 limb darkening coefficients values of u1=0.559-0.070+0.071 and u2=0.104-0.117+0.120 according to the quadratic law. These determinations are consistent with values obtained previously with different parameterization. Our result based on a smaller number of iterations per chain presents similar accuracy, of about 10-2 for u1. In addition, we explored in a more focused way the region of the parameter of interest.

Keywords:
Exoplanets; Limb Darkening; Planetary Transit Method; Photometry


1. Introdução

A detecção em 1995 de 51 Pegasi b [1[1] M. Mayor e D. Queloz, Nature 378, 355 (1995).], o primeiro planeta fora do Sistema Solar (exoplaneta) orbitando uma estrela do tipo solar da sequência principal, impulsionou a busca por exoplanetas e este tópico se tornou uma importante realidade em muitos centros de pesquisas. Este evento culminou no Prêmio Nobel de Física em 2019 para os astrônomos Michael Mayor e Didier Queloz. No caso de 51 Pegasi b, o método de detecção utilizado foi o de velocidade radial [1[1] M. Mayor e D. Queloz, Nature 378, 355 (1995).], porém a evolução da instrumentação empregada para coleta de dados e, consequentemente, a chegada de medidas mais precisas e satélites especializados na busca por exoplanetas, fizeram do método do trânsito planetário (o qual será descrito nas próximas secções) o mais promissor em termos de busca intensiva por exoplanetas [2[2] D.G. Koch, W.J. Borucki, G. Basri, N.M. Batalha, T.M. Brown, D. Caldwell, J. Christensen-Dalsgaard, W.D. Cochran, E. DeVore, E.W. Dunham et al., The Astrophysical Journal 713, L79 (2010).]. Já do ponto de vista instrumental, os telescópios espaciais contribuíram significativamente para o avanço desse tipo de pesquisa graças ao grande número de sistemas detectados, e particularmente, os telescópios espaciais CoRoT, acrônimo de Convection, Rotation and Transit da European Space Agency (ESA) [3[3] M. Auvergne, P. Bodin, L. Boisnard, J.T. Buey, S. Chaintreuil, G. Epstein, M. Jouret, T. Lam-Trong, P. Levacher, A. Magnan et al., Astronomy and Astrophysics 506, 411 (2009).] e Kepler, [2[2] D.G. Koch, W.J. Borucki, G. Basri, N.M. Batalha, T.M. Brown, D. Caldwell, J. Christensen-Dalsgaard, W.D. Cochran, E. DeVore, E.W. Dunham et al., The Astrophysical Journal 713, L79 (2010).] da NASA (National Aeronautics and Space Administration ) escreveram uma importante página na história da busca por exoplanetas. Estas duas missões espaciais foram idealizadas inicialmente para detectar planetas do tipo Terra, ou seja, planetas rochosos e com raio e semi-eixo maior semelhantes ao sistema Sol – Terra. Uma vez detectados, estes exoplanetas poderiam ser avaliados em termos de condições de habitabilidade, que se dá quando um determinado sistema possui características compatíveis com a manutenção da vida como a conhecemos [4[4] R.K. Kopparapu, R. Ramirez, J.F. Kasting, V. Eymet, T.D. Robinson, S. Mahadevan, R.C. Terrien, S. Domagal-Goldman, V. Meadows e R. Deshpande, The Astrophysical Journal 765, 131 (2013).]. Contudo, análises dos planetas detectados mostram que a maior parte destes apresentam um raio orbital muito maior que o raio orbital da Terra e estão mais próximos da sua estrela hospedeira, quando comparado com o sistema Sol – Terra. Em parte, isso é explicado por um viés intrínseco do método observacional e que será abordado neste estudo.

A missão espacial Kepler da NASA foi lançada em 2009 e culminou na detecção de mais de 2.700 exoplanetas.1 1 https://exoplanetarchive.ipac.caltech.edu/docs/counts_detail.html Este satélite coletou dados em duas cadências diferentes: curta (1 ponto a cada 58 segundos) e longa (1 ponto a cada 29,42 minutos). Coletando dados fotométricos na faixa visível do espectro eletromagnético, o telescópio funcionou totalmente operante até meados de 2013 quando seu sistema giroscópios começou a apresentar defeitos, fazendo com que o instrumento tivesse o seu movimento comprometido [5[5] S.B. Howell, C. Sobeck, M. Haas, M. Still, T. Barclay, F. Mullally, J. Troeltzsch, S. Aigrain, S.T. Bryson, D. Caldwell et al., Publications of the Astronomical Society of the Pacific 126, 398 (2014).]. Inicialmente, o satélite Kepler observava uma região de interesse estreita e bem definida no céu, porém com a falha do giroscópio, a missão foi redirecionada para observar diferentes regiões do céu devido à pouca possibilidade de apontamento causada pelo giroscópio defeituoso. A missão K2, como foi batizada, apesar das suas limitações, gerou ainda muitos resultados interessantes, como por exemplo, a detecção de novos exoplanetas [5[5] S.B. Howell, C. Sobeck, M. Haas, M. Still, T. Barclay, F. Mullally, J. Troeltzsch, S. Aigrain, S.T. Bryson, D. Caldwell et al., Publications of the Astronomical Society of the Pacific 126, 398 (2014).] em diversas regiões da Galáxia. Em seguida, aconteceu o lançamento do satélite TESS [6[6] G.R. Ricker, J.N. Winn, R. Vanderspek, D.W. Latham, G.Á. Bakos, J.L. Bean, Z.K. Berta-Thompson, T.M. Brown, L. Buchhave, N.R. Butler et al., Journal of Astronomical Telescopes, Instruments, and Systems 1, 014003 (2015).] que atualmente se encontra em voo e fornece medidas fotométricas para milhares de estrelas brilhantes e provenientes de várias posições da nossa Galáxia, observando quase 85% do céu.

Nesse contexto, a elaboração de missões espaciais cada vez de melhor desempenho e telescópios mais precisos com detectores mais sensíveis, fazem parte dos desafios constantes da área de busca por exoplanetas. Além desta linha, a física estelar e os modelos de evolução estelar se beneficiam muito dos produtos observacionais gerados nestes últimos 25 anos de busca por exoplanetas. Este cenário deve continuar ainda pelos próximos anos e, neste sentido, a preparação para missões como o Planetary Transits and Oscillations of stars (PLATO) [7[7] H. Rauer, C. Catala, C. Aerts, T. Appourchaux, W. Benz, A. Brandeker, J. Christensen-Dalsgaard, M. Deleuil, L. Gizon, M.J. Goupil, Experimental Astronomy 38, 249 (2014).] da ESA, e o James Webb Space Telescope (JWST) [8[8] J.P. Gardner, J.C. Mather, M. Clampin, R. Doyon, M.A. Greenhouse, H.B. Hammel, J.B. Hutchings, P. Jakobsen, S.J. Lilly e K.S. Long, Space Science Reviews 123, 485 (2006).], da NASA, deve impulsionar ainda mais a investigação sobre a física da formação dos exoplanetas e alargar a cobertura estatística de parâmetros importantes para compreender melhor os efeitos de segunda ordem presentes em um trânsito planetário. Para dar um exemplo, a precisão fotométrica que era de uma parte por milhão nas últimas missões espaciais [2[2] D.G. Koch, W.J. Borucki, G. Basri, N.M. Batalha, T.M. Brown, D. Caldwell, J. Christensen-Dalsgaard, W.D. Cochran, E. DeVore, E.W. Dunham et al., The Astrophysical Journal 713, L79 (2010).], deve aumentar nos próximos anos e isto terá impacto considerável nos modelos de trânsito planetário [9[9] S. Csizmadia, T. Pasternacki, C. Dreyer, J. Cabrera, A. Erikson e H. Rauer, Astronomy and Astrophysics 549, A9 (2013).]. Entre outros parâmetros, a profundidade do decaimento no brilho normalizado observado da curva de luz que é, primariamente, determinado pela razão entre o raio do planeta e o raio da estrela, deve passar por um crescimento na acurácia de sua medida. A busca por uma modelização cuidadosa do trânsito planetário em curvas de luz de alta precisão gera, além da fração de raios, a inclinação azimutal do sistema e estes dependem dos coeficientes de escurecimento de limbo [10[10] D. Kipping e G. Bakos, The Astrophysical Journal 730, 50 (2011).], que podem mascarar resultados e gerar efeitos que cobrem possíveis detecções de planetas pequenos como a Terra orbitando estrelas parecidas com o Sol [9[9] S. Csizmadia, T. Pasternacki, C. Dreyer, J. Cabrera, A. Erikson e H. Rauer, Astronomy and Astrophysics 549, A9 (2013).].

Nosso estudo foca exatamente na medida destes coeficientes de escurecimento de limbo e, também, em um melhor entendimento do perfil de brilho estelar. Apresentamos ao leitor uma descrição física dos elementos necessários para a modelização de um trânsito planetário com base nos dados do satélite Kepler. Uma vez que a abordagem dos estudos de escurecimento de limbo com base no método do trânsito é um desafio [11[11] H.M. Müller, K.F. Huber, S. Czesla, U. Wolter e J.H.M.M. Schmitt, Astronomy and Astrophysics 560, A112 (2013).], percebemos que devido aos seus efeitos na curva de luz, pode-se confundir com efeitos provenientes de outros parâmetros. Isso decorre do fato de que o escurecimento de limbo afeta o formato em U do trânsito, logo ele afeta tanto o ponto de mínimo da curva de luz (de onde se calcula a fração de raios) quanto os pontos de entrada e saída do planeta no trânsito a partir do qual se calcula a inclinação do sistema, por exemplo. Ressaltamos que o tratamento do escurecimento de limbo interfere na precisão da caracterização de exoplanetas [12[12] H.A. Knutson, D. Charbonneau, R.W. Noyes, T.M. Brown e R.L. Gilliland, The Astrophysical Journal 655, 564 (2007).] e na modelização das manchas estelares via curva de luz [13[13] A.A. Souza e A. Valio, Revista Brasileira de Ensino de Física 41, e20180323 (2019).]. As leis de modelização dos coeficientes de escurecimento de limbo, mais utilizadas na literatura e aplicadas são formulações lineares [14[14] E.A. Milne, Monthly Notices of the Royal Astronomical Society 81, 361 (1921).], quadráticas [15[15] A. Manduca, R.A. Bell e B. Gustafsson, Astronomy and Astrophysics 61, 809 (1977).], raiz quadrada [16[16] J. Diaz-Cordoves e A. Gimenez, Astronomy and Astrophysics 259, 227 (1992).] e logarítmicas [17[17] D.A. Klinglesmith e S. Sobieski, The Astronomical Journal 75, 175 (1970).]. Porém, em qualquer descrição adotada uma variação nos valores dos coeficientes de limbo pode causar imprecisões nas medidas inferidas e erro na fração de raios, semi-eixo maior e até na duração do trânsito [12[12] H.A. Knutson, D. Charbonneau, R.W. Noyes, T.M. Brown e R.L. Gilliland, The Astrophysical Journal 655, 564 (2007).], contaminando assim, outros parâmetros do sistema. Isto inviabiliza avanços reais que foram feitos na instrumentação e detecção ao longo dos últimos anos. Além disso, essas formulações possuem uma dependência direta com o ângulo entre o observador e a normal para a superfície da estrela.

É neste contexto que se insere este estudo, onde realizamos uma busca pelos coeficientes de escurecimento de limbo seguindo prescrições já estabelecidas [15[15] A. Manduca, R.A. Bell e B. Gustafsson, Astronomy and Astrophysics 61, 809 (1977).] e com novas implementações de buscas de soluções através do método Monte Carlo via cadeias de Markov (do inglês, MCMC Markov Chain Monte Carlo) modificada. A metodologia foi aplicada para a estrela de tipo espetral K0, KIC 5357901 (Kepler-425), previamente observada pelo satélite Kepler que, além de possuir vários trânsitos (ou seja, possui um período orbital curto) ainda possui trânsitos grandes (fração de raios grande). Neste sentido, além da modelização dos trânsitos para este objeto, contribuímos para a determinação do seu perfil de brilho. De forma mais ampla, estabelecemos uma metodologia que pode ser aplicada para o estudo de uma larga base de estrelas já observadas pelo satélite.

2. O Escurecimento de Limbo

O escurecimento de limbo é um fenômeno onde se tem a diminuição gradual do brilho do disco de uma estrela conforme se observa desde o centro da estrela até a sua borda (vide Figura1). Tal fenômeno ocorre devido à distribuição de temperatura não uniforme observada em uma estrela. Além disso, esse efeito depende da temperatura efetiva da estrela, ou seja, é função do tipo espectral. A física fundamental associada ao escurecimento de limbo advém da lei de Stefan-Boltzmann, que associa a potência total irradiada por unidade de área superficial de um corpo negro (que por sua vez é um objeto que absorve e emite toda a radiação eletromagnética nele incidida, ou seja, nenhuma luz o atravessa ou é refletida) e é diretamente proporcional à quarta potência da sua temperatura termodinâmica (ou seja, sua temperatura medida na escala absoluta), como representado na equação

(1) j * = σ T 4 ,

Figura 1:
Fotosfera solar no dia 30 de janeiro de 2021, com escurecimento do limbo perceptível. Dados: NASA/SDO.

onde σ é a constante de proporcionalidade de Stefan-Boltzmann cujo valor aproximado é de 5,610-8 W m 2K 4. Essa lei é mais bem explorada para o caso do Sol no trabalho de P. A. D. Almeida J. & Gregorio-Hetem [18[18] P.A.D. Almeida e J. Gregorio-Hetem, Revista Brasileira de Ensino de Física 44, e202100405 (2022).].

As estrelas são boas aproximações para corpos negros devido a serem corpos muito opacos, ou seja, o material estelar é um ótimo absorvedor de radiação.

Dessa forma, é importante destacar que a fotosfera apresenta o fenômeno físico de extinção que espalha e absorve parte da luz que a atravessa, assim como ocorre na atmosfera da Terra (no caso terrestre, a extinção causa dificuldades para fazer observações astronômicas próximo a regiões metropolitanas, por exemplo). Nessa conjuntura, faz sentido introduzirmos o conceito de profundidade óptica, dado por τν. A profundidade óptica descreve a fração de radiação que é absorvida ou emitida (espalhada) em um caminho de a até b, e é dada pela integral abaixo

(2) τ ν = a b κ ν ρ d x .

Veja que de forma geral, κν representa o coeficiente de absorção (que por sua vez é função da frequênciaν) e ρ representa a densidade do meio. De maneira equivalente, a profundidade óptica pode ser interpretada como uma medida da transparência de um meio, para uma frequência específica de luz ν.

Do ponto de vista da teoria de transferência radiativa pode-se considerar a profundidade óptica como uma medida da profundidade geométrica de uma atmosfera estelar, como pode-se ver representado esquematicamente pela Figura2. No limite superior da camada, a profundidade óptica é zero (ou seja, toda energia é transmitida), enquanto que no centro da estrela, ela tende a infinito (toda energia é absorvida ou espalhada). Enquanto que na profundidade óptica τ=1 constrói-se uma superfície teórica onde apenas 1/e da energia é transmitida, sem ser absorvida ou espalhada. Nesse sentido, podemos ver que a superfície τ=1, representada na Figura2, na região próxima ao centro, está associada a uma camada com a temperatura mais alta em relação à mesma superfície na extremidade da estrela. Nesse sentido, como a temperatura mais alta está relacionada a uma potência irradiada mais alta (pela lei de Stefan–Boltzmann), dá-se o fenômeno do escurecimento do limbo.

Figura 2:
Representação física do escurecimento de limbo. A linha tracejada indica a profundidade óptica τ=1. As setas representam a direção da linha de visada do observador. Enquanto as cores, nessa figura, representam de forma figurativa as temperaturas em forma de gradiente da mais quente (azul) para a mais fria (vermelho). Fonte: Muller (2016) [19[19] H.M. Muller, Limb-darkening Measurements on Exoplanet Host Stars and the Sun. Dissertação de Mestrado, Universität Hamburg, Hamburg (2015).].

3. O Trânsito Planetário

A base do nosso estudo reside na modelagem do trânsito planetário e este consiste na passagem do planeta pela linha de visada que conecta o observador até a estrela [20[20] K. Mandel e E. Agol, The Astrophysical Journal 580, L171 (2002).] como mostrado na Figura3. As ocultações periódicas geram variações também periódicas do brilho estelar. A condição geométrica do alinhamento observador – planeta – estrela é suficiente para tornar esse evento estatisticamente raro. Os três corpos estarão perfeitamente alinhados quando o ângulo de inclinação i for igual a π2 (Figura3). O ângulo complementar à i, e descrito por θ0 representa o limite superior para a detecção do sistema.

Figura 3:
Esquema representando a geometria de um sistema típico estrela – planeta – observador. É apresentado também o ângulo da linha da visada assim como o eixo de rotação da estrela.

Podemos calcular a probabilidade 𝒫 desse alinhamento considerando o caso onde a distância do observador seria muito maior que o raio da estrela (dR*). Tal probabilidade pode ser obtida através do ângulo formado pelo sistema observador-planeta-estrela, mais especificamente, a função densidade de probabilidade correspondente é proporcional a cosθ, o que implica que a probabilidade é dada por:

(3) 𝒫 = 0 θ 0 cos θ d θ 0 π / 2 cos θ d θ = sin θ 0 = R * a .

Isso significa que a probabilidade de um sistema exoplanetário qualquer estar alinhado com um observador é proporcional ao raio estelar e inversamente proporcional ao semi-eixo maior. Substituindo valores, para o caso onde quiséssemos, por exemplo, detectar o planeta com o semi-eixo maior de Júpiter orbitando uma estrela como o Sol, teríamos:

(4) 𝒫 = R * a = 0 , 00465 U.A. 5 , 2 U.A. 0 , 089 % ,

onde U.A. é a unidade astronômica de comprimento que equivale a, aproximadamente, à distância Terra – Sol, ou seja, cerca de 150 milhões de quilômetros.

Ou seja, estatisticamente, seria necessário observar 1/(0,089%)1120 sistemas como esse descrito acima para encontrar um deles alinhado de forma que a linha de visada do observador favorecesse uma ocultação do fluxo estelar, e assim possibilitasse, talvez, uma detecção pelo método do trânsito planetário. Uma vez que este evento ocorresse, sua detecção através das variações periódicas no brilho da estrela, poderia ser representada por medidas fotométricas do fluxo Fλ em um gráfico do tipo fluxo luminoso em função do tempo, como representado na Figura4. Essa forma gráfica é uma série temporal, conhecida também como curva de luz, e é a partir dela que fazemos as principais análises físicas do problema da modelização dos trânsitos planetários.

Figura 4:
Exemplo de uma curva de luz de um trânsito planetário. Na parte interna uma ampliação de um dos trânsitos e a medida da queda do fluxo causado pela passagem do planeta.

Uma vez apresentada a geometria do trânsito, assim como sua probabilidade de ocorrência, descreveremos abaixo a modelagem propriamente dita e separamos esta descrição em dois casos, o primeiro e mais simples, que considera a estrela como uma fonte uniforme, e o segundo, que considera o escurecimento de limbo.

3.1. Caso mais simples de fonte uniforme

O modelo a ser descrito foi primeiramente apresentado em [20[20] K. Mandel e E. Agol, The Astrophysical Journal 580, L171 (2002).] e é amplamente utilizado nos mais diversos projetos de exoplanetologia. Esse modelo, na sua forma mais simplificada, apresenta cinco parâmetros principais:

t 0 Tempo de conjunção inferior a / r * Raio da órbita normalizado pelo raio estelar P orb Período orbital do planeta i Inclinação do sistema. r p / r * ou p Raio do planeta normalizado pelo raio estelar

Essa forma mais simplificada do modelo tem como primeira premissa que a órbita do planeta é circular e com isso evita-se dois parâmetros, a excentricidade (e) e o argumento de longitude de periastro (ω) relacionado à órbita elíptica. Além disso, outra premissa usada nesta formulação é que a estrela não possui escurecimento de limbo, com isso, é como se a estrela tivesse o mesmo brilho por toda superfície projetada. Ambas essas premissas são incomuns (ou seja, não são conservadoras) e longe da realidade astrofísica. A estimativa de órbita circular, na verdade, pode ser aplicada mais comumente quando se tem conhecimento dessa característica do sistema através de outro método como, por exemplo, medidas de velocidade radial. Dessa forma, o percentual de fluxo da estrela não ocultado pelo planeta é Fe(p,z)=1-λe(p,z), onde z é a distância projetada entre o centro da estrela e o centro do planeta normalizado pelo raio estelar e z=d/r*, como mostrado na Figura5 e λe(p,z) é dado pela Equação (5).

(5) λ e ( p , z ) = { 0 , , d > r p + r * . 1 π [ p 2 k 0 + k 1 + k 2 ] , | r * - r p | < d r p + r * . p 2 , , d r * - r p . 1 , , r p d + r * .

Figura 5:
Esquema de modelização do trânsito. Adaptado de [20[20] K. Mandel e E. Agol, The Astrophysical Journal 580, L171 (2002).].

onde k0=cos-1[(p2+z2-1)/2pz], k1=cos-1[(1-p2+z2)/2z] e k2=4z2-(1+z2-p2)24.

Percebe-se que a Equação (5) está dividida em alguns casos. Quando d>rp+r* o planeta está totalmente fora do trânsito. Em |r*-rp|<drp+r* o planeta está entrando ou saindo do trânsito, ou seja, uma parte do planeta oculta o brilho da estrela enquanto outra parte não. Já em dr*-rp o planeta está totalmente na frente da estrela. Toda a sua área oculta o brilho estelar, portanto, é nesse trecho que o fluxo ocultado é máximo. Por fim, em rpd+r* teríamos um caso em que o raio do planeta seria maior que o raio estelar e, assim, bloquearia todo o fluxo emitido pela estrela. Com isso, a função Fe(p,z) será a curva de luz. Para essas premissas, podemos identificar alguns parâmetros físicos diretamente da curva de luz. Por exemplo, a fração de raios (p ou rp/r*) é dada pela profundidade máxima de variação no brilho ΔF como indicado na Equação (6).

(6) Δ F F total - F transito F total = ( r p r * ) 2 .

No trabalho de Saeger & Mallén-Ornelas [21[21] S. Seager e G. Mallén-Ornelas, Scientific Frontiers in Research on Extrasolar Planets 294, 419 (2003).], foi mostrado que para esse caso mais simplificado com 5 parâmetros existe uma solução única para a curva de luz, dado que tenhamos a relação massa-raio da estrela no limite que M*Mp, juntamente com as outras premissas já apresentadas no início desta seção.

Tomando o limite que o raio da órbita é muito maior que o raio estelar (ar*), é possível chegar nas relações analíticas para os valores de duração total do trânsito (tT), e duração do trânsito na situação em que o planeta já está totalmente na frente da estrela (tF). Isto levando em conta o referencial do ponto de vista do observador, ou seja, em nossa parametrização d<(r*-rp). Representado na Equação (7).

(7) ( t F t T ) 2 = ( 1 - r p r * ) 2 - ( a r * cos i ) 2 ( 1 + r p r * ) 2 - ( a r * cos i ) 2 .

Por conveniência, definimos uma nova variável chamada de parâmetro de impacto (b) que tem um sentido geométrico por ser o semi-eixo maior normalizado pelo raio estelar projetado pela inclinação, como representado na Figura3. Assim, ele é definido pela Equação (8).

(8) b = a r * cos i .

Com o acréscimo do parâmetro b podemos reescrever a Equação (7) da forma.

(9) ( t F t T ) 2 = ( 1 - r p r * ) 2 - b 2 ( 1 + r p r * ) 2 - b 2 .

É importante frisar que as equações apresentadas nessa seção consideram órbitas circulares, ou seja, com excentricidade muito próxima de zero ou zero. Mesmo que seja do conhecimento da física que as órbitas devem seguir secções cônicas, para esse problema, podemos descartar as órbitas não circulares devido ao curto período orbital do planeta. Nesta configuração, é razoável que a dissipação de maré amorteça qualquer nível inicial de excentricidade orbital para um patamar irrisório.

4. Modelagem dos Trânsitos com Escurecimento de Limbo

Tecnicamente, no nível da operacionalização computacional deste trabalho, usamos a implementação das equações descritas acima através da biblioteca batman : BAsic Transit Model cAlculatioN in Python [22[22] L. Kreidberg, Publications of the Astronomical Society of the Pacific 127, 1161 (2015).] em Python, que oferece uma interface bastante amigável para simular curvas de luz, além de estar implementada em uma linguagem de mais baixo nível que o Python, portanto, trazendo uma vantagem no custo computacional da simulação. Para uma curva de luz com 100 pontos de trânsito, a biblioteca batman pode simular um milhão de curvas de luz em 30 segundos para o modelo de escurecimento de limbo quadrático com um único processador Intel Core i5 de 1,7 GHz [22[22] L. Kreidberg, Publications of the Astronomical Society of the Pacific 127, 1161 (2015).].

Para modelar o efeito do escurecimento é definida uma nova parametrização μ como sendo o cosseno de θ, onde θ é o ângulo entre o observador e o vetor normal a superfície projetada da estrela, como representado na Figura5. Na busca pela melhor modelização, foram propostas várias definições, com formulações lineares, exponenciais e até mesmo uma formulação não linear, entretanto, a lei que é mais utilizada e aceita na literatura para trabalhos com o método do trânsito, é a lei quadrática (vide os trabalhos de Espinoza & Jordán [23[23] N. Espinoza e A. Jordán, Monthly Notices of the Royal Astronomical Society 457, 3573 (2016).] e Kipping [24[24] D.M. Kipping, Monthly Notices of the Royal Astronomical Society 435, 2152 (2013).]), definida por Manduca et al. [15[15] A. Manduca, R.A. Bell e B. Gustafsson, Astronomy and Astrophysics 61, 809 (1977).], de acordo com a Equação (10):

(10) I ( μ ) I ( 1 ) = 1 - u 1 ( 1 - μ ) - u 2 ( 1 - μ ) 2 ,

onde I(1) é brilho no centro da estrela, μ=cosθ, sendo θ o ângulo entre o observador e a normal para a superfície da estrela e u1 e u2 são os coeficientes de escurecimento de limbo da estrela que entrarão como variáveis livres em nosso modelo.

Veja que esses coeficientes não podem ter quaisquer valores. Existem premissas físicas que devem ser respeitadas, dentre elas, a intensidade de brilho da estrela deve sempre diminuir (nunca aumentar) com o raio, ou seja, a derivada de I(θ) em relação a θ deve ser negativa e, em nossa parametrização, a derivada da Equação (10) em relação a μ deve ser positiva. Dessa forma, podemos chegar em um primeiro vínculo na Equação (11).

(11) d I ( μ ) d μ > 0 , u 1 + 2 u 2 ( 1 - μ ) > 0 .

Além disso, a estrela deve sempre ter brilho positivo, portanto, a Equação (10) deve ser sempre positiva. Com isso, podemos obter mais uma restrição com a Equação (12).

(12) I ( μ ) > 0 , u 1 ( 1 - μ ) + u 2 ( 1 - μ ) 2 < 1 .

A partir das Equações (11) e (12), ao analisar os pontos de contorno do problema, ou seja, em μ0 (θ=900) e em μ1 (θ=00), chegamos nos limites para u1 e u2, como representado na Equação (13).

(13) u 1 + u 2 < 1 , u 1 > 0 , u 1 + 2 u 2 > 0 .

Diante disso, com as três equações apresentadas em (13), obtemos um triângulo de interesse que contém os valores físicos para os coeficientes de escurecimento de limbo, como representado na Figura6. Entretanto, fazer uma amostragem desse triângulo pode não ser trivial. Para isso, vamos usar uma reparametrização de u1 e u2, com base na elegante ideia proveniente da computação gráfica apresentada por Kipping [24[24] D.M. Kipping, Monthly Notices of the Royal Astronomical Society 435, 2152 (2013).], que é uma forma inteligente de simplificar as restrições impostas pelas equações (13) sem acrescentar nenhuma variável livre às variáveis paramétricas q1 e q2 pela Equação (14).

(14) q 1 = ( u 1 + u 2 ) 2 , q 2 = u 1 2 ( u 1 + u 2 ) .

Figura 6:
Representação das equações (13), a área demarcada em cinza representa o espaço fisicamente válido entre u1 e u2.

Com isso, ao propor uma amostragem no intervalo entre zero e um, fazemos com que q1 e q2 tenham acesso a todo o espaço físico de possibilidade dos coeficientes de escurecimento de limbo. Com base no descrito, temos que a transformada inversa será determinada pela Equação (15).

(15) u 1 = 2 q 1 q 2 , u 2 = q 1 ( 1 - 2 q 2 ) .

5. Dados Fotométricos da Missão Kepler

A missão Kepler foi desenvolvida com foco nos exoplanetas e para responder à pergunta: “Qual é a frequência de ocorrência de outras Terras em nossa Galáxia?”. Em particular, a missão sondou a frequência estatística dos planetas do tamanho da Terra em zonas habitáveis de estrelas do tipo solar. O satélite foi lançado no dia 7 de março de 2009, em uma órbita heliocêntrica. Para nossa análise, usamos os dados do Kepler obtidos no Mikulski Archive for Space Telescopes (MAST).2 2 https://archive.stsci.edu/index.html Usamos ainda os dados de curta cadência (a cada 58 segundos) em detrimento da alta cadência (a cada 29,42 minutos) para obter uma descrição mais contínua das curvas de luz, além da fotometria de abertura corrigida PDCSAP (Pre-search Data Conditioning Simple Aperture Photometry) que tem correções de determinadas tendências instrumentais de longo prazo. Para uma análise ótima do escurecimento de limbo, faz-se necessário a escolha de um objeto com alta razão sinal-ruído, além de ter um período orbital curto para podermos ter o maior número de trânsitos possíveis. Nesse contexto, selecionamos para análise a curva de luz da estrela KIC 5357901. Esta é uma estrela de tipo espectral K0, com temperatura efetiva de cerca de 5170 K e uma luminosidade de 0,86L[25[25] G. Hébrard, A. Santerne, G. Montagnier, G. Bruno, M. Deleuil, M. Havel, J.M. Almenara, C. Damiani, S.C.C. Barros, A.S. Bonomo et al., Astronomy and Astrophysics 572, A93 (2014).] e hospeda o planeta Kepler-425b com período orbital de cerca de 3,8 dias.

Ao analisar os dados provenientes do satélite Kepler, vemos que as curvas de luz apresentam uma modulação que deve ser tratada antes de aplicar o modelo do trânsito planetário. Antes de normalizar, separamos a curva de luz em duas partes, uma parte que continha os trânsitos do exoplaneta e as partes sem trânsito (contínuo). É importante perceber que em um cenário teórico o contínuo deveria ser constante, mas pequenas modulações são observadas decorrentes de efeitos instrumentais, resíduos do processo de redução e também de atividade estelar, manchas solares. Fazemos uma normalização do contínuo usando mínimos quadrados com um polinômio de Chebyshev de primeira ordem e grau 90. O ajuste em questão está representado na Figura7. Além disso, para remover pontos fora da média, outliers, usamos um outro filtro que exclui todos os pontos a mais de 4σ de distância em relação à média do contínuo.

Figura 7:
Preto: Pseudocontínuo. Vermelho: Trânsitos. Azul: ajuste do contínuo com polinômio de Chebyshev.

Aqui vale ressaltar que essas medidas são fotométricas e que a normalização é feita em torno de um pseudocontínuo. Dessa forma, o fato do fluxo normalizado oscilar em torno de um é decorrente do método de normalização da própria missão Kepler e da junção dos quarters da curva de luz. No geral, ter o valor do fluxo normalizado em torno da unidade tem pouca influência na curva de luz uma vez que estamos interessados somente nos deeps (que são as quedas no brilho da estrela que sinalizam a passagem do planeta) da mesma.

6. Modelagem da Curva de Luz e Determinação dos Parâmetros de Trânsito

Para determinar os coeficientes de escurecimento de limbo é necessário primeiro um conhecimento preciso dos parâmetros de primeira ordem no modelo do trânsito, principalmente do período orbital (Porb ) e tempo do centro de trânsito (t0). Claramente, os valores imprecisos de Porb e t0 produzem impossibilidades na determinação correta dos parâmetros restantes na descrição do trânsito, pois tanto variações no Porb quanto no t0 deslocam horizontalmente (na coordenada temporal) o modelo da curva de luz do trânsito. Portanto, é essencial determinar Porb e t0 com a maior precisão possível. Como primeira etapa, usamos o modelo de trânsito para ajustar esses dois parâmetros. Neste passo optamos pelo uso do algoritmo downhill simplex, mantendo todos os outros parâmetros fixos. O algoritmo downhill simplex (ou algoritmo Nelder–Mead) é um método numérico de busca direta, e que é projetado para obter o mínimo ou máximo de uma função em um espaço multidimensional. É um método comumente aplicado a problemas de otimização não lineares, especialmente para os quais as derivadas podem não ser conhecidas e, portanto, o método de Newton-Raphson não pode ser utilizado. Depois de determinar Porb e t0, ajustamos os parâmetros restantes: a razão dos raios estrela-planeta, o ângulo de inclinação, o semi-eixo maior e os coeficientes de escurecimento de limbo quadráticos u1 e u2, mantendo dessa vez os valores de Porb e t0 já minimizados fixos, novamente utilizando o algoritmo downhill simplex. Essa minimização inicial é importante, pois usaremos esses resultados de ajuste como valores iniciais para nossa amostragem MCMC, como descrita a seguir na Seção7 7. Abordagem por Simulações Monte Carlo via Cadeias de Markov Métodos como o MCMC são cada vez mais frequentes em ciências guiadas por dados. Objetivamente, esses métodos são usados para obter amostragens da distribuição a posteriori da quantidade de interesse e são especialmente convenientes quando o problema em questão envolve muitos parâmetros livres ou quando os modelos são caros computacionalmente e as observações têm uma baixa razão sinal-ruído. Elementos muito frequentes em muitos ramos da física e em análises astrofísicas, como no caso desse trabalho. Fundamentalmente, uma cadeia de Markov é uma sucessão de estados, de forma que cada estado só é dependente do estado exatamente anterior a ele, já Monte Carlo implica que essa passagem de um estado para o outro vai ser determinada de forma não determinística. Ao tratar de distribuições a priori e a posteriori, estamos usando aqui o paradigma Bayesiano para estatística. Na estrutura Bayesiana, como todas as inferências são baseadas nas posterioris, existe uma base matemática bem estabelecida para quantificar as incertezas nos parâmetros do modelo, ao contrário dos métodos “frequencistas”, que normalmente resultam em estimativas pontuais dos parâmetros, já que o conceito de probabilidade é diferente de acordo com as interpretações “frequencistas” e Bayesianas. O principal aspecto da estatística na abordagem Bayesiana é o teorema de Bayes que pode ser derivado partindo de uma probabilidade conjunta p⁢(x,θ), reescrevendo a probabilidade de x pela integral da probabilidade marginal em θ, ou seja p⁢(x)=∫p⁢(x,θ)⁢dθ. Dessa forma, juntamente com a definição de probabilidade conjunta, (p⁢(x,θ)=p⁢(x)⁢p⁢(θ|x)), o teorema de Bayes, pode ser definido pela Equação (16). Sendo x os dados e θ o conjunto de parâmetros do modelo. (16) p ⁢ ( θ | x ) = p ⁢ ( x , θ ) p ⁢ ( x ) = p ⁢ ( θ ) ⁢ p ⁢ ( x | θ ) ∫ p ⁢ ( x , θ ) ⁢ d θ . Veja que o resultado da integral no denominador da Equação (16) é um escalar que não depende de θ. Assim, funciona como uma constante de normalização para a Equação, portanto, o teorema de Bayes assume outra forma mais usual numericamente pela Equação (17) (17) p ⁢ ( θ | x ) ∝ p ⁢ ( θ ) ⁢ p ⁢ ( x | θ ) , onde p⁢(θ|x) é a distribuição a posteriori a qual buscaremos amostrar, p⁢(x|θ) é a likelihood, ou verossimilhança, e p⁢(θ) é a distribuição a priori, que quantifica algum conhecimento que temos a priori sobre o sistema como, principalmente em nosso caso, os limites físicos para todos os parâmetros e a amostragem inteligente do escurecimento de limbo usando a parametrização descrita na Seção4. Essa amostragem inteligente consiste numa restrição na busca de parâmetros nos limites físicos do escurecimento de limbo, com a parametrização usada no trabalho de Kipping [24]. Dessa forma, não há gasto de tempo computacional em regiões não-físicas. O termo “inteligente” advém do fato de que os coeficientes de escurecimento de limbo não se limitam apenas à região física do problema. 7.1. Verossimilhança e distribuições a priori No contexto bayesiano, a verossimilhança quantifica a qualidade do ajuste de um modelo estatístico em relação a uma amostra de dados na qual queremos determinados valores de parâmetros desconhecidos. Se partimos de uma premissa que nossos dados apresentam barras de erros corretas, gaussianas e independentes, podemos expressar a verossimilhança como uma distribuição normal e, por estabilidade numérica, usaremos o logaritmo da verossimilhança conforme a Equação (18). (18) ln ⁡ p ⁢ ( x | θ ) = - 1 2 ⁢ ∑ j = 1 N [ ( x j - θ j ) 2 σ j 2 + ln ⁡ ( 2 ⁢ π ⁢ σ j 2 ) ] , onde N é o número de pontos nos dados e σj representa o erro do ponto j. Nesse trabalho, construímos um algoritmo para cobrir uma variação de 5 parâmetros, rp/r⋆, a, i, q1 e q2, enquanto mantivemos fixados Porb e t0. Manter Porb e t0 como parâmetros fixos é razoável pois, como são parâmetros que interferem muito no ajuste da curva de luz, temos estimativas desses parâmetros com pequenas barras de erro, enquanto que o Porb ainda pode ser obtido por outras técnicas além destas abordadas neste trabalho como, por exemplo, períodogramas de Lomb-Scargle [26]. Já para as distribuições a priori, pelo princípio da indiferença, considerando que outras informações não são suficientemente relevantes para tomarem forma de distribuições analíticas específicas, vamos usar prioris não informativas definindo uma distribuição uniforme normalizada por todo o espaço físico de cada parâmetro, facilitando assim, a exploração das cadeias pelo espaço físico rejeitando apenas os valores não físicos dos parâmetros. .

7. Abordagem por Simulações Monte Carlo via Cadeias de Markov

Métodos como o MCMC são cada vez mais frequentes em ciências guiadas por dados. Objetivamente, esses métodos são usados para obter amostragens da distribuição a posteriori da quantidade de interesse e são especialmente convenientes quando o problema em questão envolve muitos parâmetros livres ou quando os modelos são caros computacionalmente e as observações têm uma baixa razão sinal-ruído. Elementos muito frequentes em muitos ramos da física e em análises astrofísicas, como no caso desse trabalho.

Fundamentalmente, uma cadeia de Markov é uma sucessão de estados, de forma que cada estado só é dependente do estado exatamente anterior a ele, já Monte Carlo implica que essa passagem de um estado para o outro vai ser determinada de forma não determinística. Ao tratar de distribuições a priori e a posteriori, estamos usando aqui o paradigma Bayesiano para estatística. Na estrutura Bayesiana, como todas as inferências são baseadas nas posterioris, existe uma base matemática bem estabelecida para quantificar as incertezas nos parâmetros do modelo, ao contrário dos métodos “frequencistas”, que normalmente resultam em estimativas pontuais dos parâmetros, já que o conceito de probabilidade é diferente de acordo com as interpretações “frequencistas” e Bayesianas. O principal aspecto da estatística na abordagem Bayesiana é o teorema de Bayes que pode ser derivado partindo de uma probabilidade conjunta p(x,θ), reescrevendo a probabilidade de x pela integral da probabilidade marginal em θ, ou seja p(x)=p(x,θ)dθ. Dessa forma, juntamente com a definição de probabilidade conjunta, (p(x,θ)=p(x)p(θ|x)), o teorema de Bayes, pode ser definido pela Equação (16). Sendo x os dados e θ o conjunto de parâmetros do modelo.

(16) p ( θ | x ) = p ( x , θ ) p ( x ) = p ( θ ) p ( x | θ ) p ( x , θ ) d θ .

Veja que o resultado da integral no denominador da Equação (16) é um escalar que não depende de θ. Assim, funciona como uma constante de normalização para a Equação, portanto, o teorema de Bayes assume outra forma mais usual numericamente pela Equação (17)

(17) p ( θ | x ) p ( θ ) p ( x | θ ) ,

onde p(θ|x) é a distribuição a posteriori a qual buscaremos amostrar, p(x|θ) é a likelihood, ou verossimilhança, e p(θ) é a distribuição a priori, que quantifica algum conhecimento que temos a priori sobre o sistema como, principalmente em nosso caso, os limites físicos para todos os parâmetros e a amostragem inteligente do escurecimento de limbo usando a parametrização descrita na Seção4 4. Modelagem dos Trânsitos com Escurecimento de Limbo Tecnicamente, no nível da operacionalização computacional deste trabalho, usamos a implementação das equações descritas acima através da biblioteca batman : BAsic Transit Model cAlculatioN in Python [22] em Python, que oferece uma interface bastante amigável para simular curvas de luz, além de estar implementada em uma linguagem de mais baixo nível que o Python, portanto, trazendo uma vantagem no custo computacional da simulação. Para uma curva de luz com 100 pontos de trânsito, a biblioteca batman pode simular um milhão de curvas de luz em 30 segundos para o modelo de escurecimento de limbo quadrático com um único processador Intel Core i5 de 1,7 GHz [22]. Para modelar o efeito do escurecimento é definida uma nova parametrização μ como sendo o cosseno de θ, onde θ é o ângulo entre o observador e o vetor normal a superfície projetada da estrela, como representado na Figura5. Na busca pela melhor modelização, foram propostas várias definições, com formulações lineares, exponenciais e até mesmo uma formulação não linear, entretanto, a lei que é mais utilizada e aceita na literatura para trabalhos com o método do trânsito, é a lei quadrática (vide os trabalhos de Espinoza & Jordán [23] e Kipping [24]), definida por Manduca et al. [15], de acordo com a Equação (10): (10) I ⁢ ( μ ) I ⁢ ( 1 ) = 1 - u 1 ⁢ ( 1 - μ ) - u 2 ⁢ ( 1 - μ ) 2 , onde I⁢(1) é brilho no centro da estrela, μ=cos⁡θ, sendo θ o ângulo entre o observador e a normal para a superfície da estrela e u1 e u2 são os coeficientes de escurecimento de limbo da estrela que entrarão como variáveis livres em nosso modelo. Veja que esses coeficientes não podem ter quaisquer valores. Existem premissas físicas que devem ser respeitadas, dentre elas, a intensidade de brilho da estrela deve sempre diminuir (nunca aumentar) com o raio, ou seja, a derivada de I⁢(θ) em relação a θ deve ser negativa e, em nossa parametrização, a derivada da Equação (10) em relação a μ deve ser positiva. Dessa forma, podemos chegar em um primeiro vínculo na Equação (11). (11) d ⁢ I ⁢ ( μ ) d ⁢ μ > 0 , u 1 + 2 ⁢ u 2 ⁢ ( 1 - μ ) > 0 . Além disso, a estrela deve sempre ter brilho positivo, portanto, a Equação (10) deve ser sempre positiva. Com isso, podemos obter mais uma restrição com a Equação (12). (12) I ⁢ ( μ ) > 0 , u 1 ⁢ ( 1 - μ ) + u 2 ⁢ ( 1 - μ ) 2 < 1 . A partir das Equações (11) e (12), ao analisar os pontos de contorno do problema, ou seja, em μ→0 (θ=900) e em μ→1 (θ=00), chegamos nos limites para u1 e u2, como representado na Equação (13). (13) u 1 + u 2 < 1 , u 1 > 0 , u 1 + 2 ⁢ u 2 > 0 . Diante disso, com as três equações apresentadas em (13), obtemos um triângulo de interesse que contém os valores físicos para os coeficientes de escurecimento de limbo, como representado na Figura6. Entretanto, fazer uma amostragem desse triângulo pode não ser trivial. Para isso, vamos usar uma reparametrização de u1 e u2, com base na elegante ideia proveniente da computação gráfica apresentada por Kipping [24], que é uma forma inteligente de simplificar as restrições impostas pelas equações (13) sem acrescentar nenhuma variável livre às variáveis paramétricas q1 e q2 pela Equação (14). (14) q 1 = ( u 1 + u 2 ) 2 , q 2 = u 1 2 ⁢ ( u 1 + u 2 ) . Figura 6: Representação das equações (13), a área demarcada em cinza representa o espaço fisicamente válido entre u1 e u2. Com isso, ao propor uma amostragem no intervalo entre zero e um, fazemos com que q1 e q2 tenham acesso a todo o espaço físico de possibilidade dos coeficientes de escurecimento de limbo. Com base no descrito, temos que a transformada inversa será determinada pela Equação (15). (15) u 1 = 2 ⁢ q 1 ⁢ q 2 , u 2 = q 1 ⁢ ( 1 - 2 ⁢ q 2 ) . .

Essa amostragem inteligente consiste numa restrição na busca de parâmetros nos limites físicos do escurecimento de limbo, com a parametrização usada no trabalho de Kipping [24[24] D.M. Kipping, Monthly Notices of the Royal Astronomical Society 435, 2152 (2013).]. Dessa forma, não há gasto de tempo computacional em regiões não-físicas. O termo “inteligente” advém do fato de que os coeficientes de escurecimento de limbo não se limitam apenas à região física do problema.

7.1. Verossimilhança e distribuições a priori

No contexto bayesiano, a verossimilhança quantifica a qualidade do ajuste de um modelo estatístico em relação a uma amostra de dados na qual queremos determinados valores de parâmetros desconhecidos. Se partimos de uma premissa que nossos dados apresentam barras de erros corretas, gaussianas e independentes, podemos expressar a verossimilhança como uma distribuição normal e, por estabilidade numérica, usaremos o logaritmo da verossimilhança conforme a Equação (18).

(18) ln p ( x | θ ) = - 1 2 j = 1 N [ ( x j - θ j ) 2 σ j 2 + ln ( 2 π σ j 2 ) ] ,

onde N é o número de pontos nos dados e σj representa o erro do ponto j. Nesse trabalho, construímos um algoritmo para cobrir uma variação de 5 parâmetros, rp/r, a, i, q1 e q2, enquanto mantivemos fixados Porb e t0. Manter Porb e t0 como parâmetros fixos é razoável pois, como são parâmetros que interferem muito no ajuste da curva de luz, temos estimativas desses parâmetros com pequenas barras de erro, enquanto que o Porb ainda pode ser obtido por outras técnicas além destas abordadas neste trabalho como, por exemplo, períodogramas de Lomb-Scargle [26[26] N.R. Lomb, Astrophysics and Space Science 39, 447 (1976).]. Já para as distribuições a priori, pelo princípio da indiferença, considerando que outras informações não são suficientemente relevantes para tomarem forma de distribuições analíticas específicas, vamos usar prioris não informativas definindo uma distribuição uniforme normalizada por todo o espaço físico de cada parâmetro, facilitando assim, a exploração das cadeias pelo espaço físico rejeitando apenas os valores não físicos dos parâmetros.

8. Resultados

Determinamos o perfil de brilho do objeto KIC 5357901 com auxílio da técnica de MCMC com 200 walkers independentes e com 5000 iterações para amostrar distribuições à posteriori dos parâmetros do sistema. Utilizamos um descarte de 40% para as primeiras iterações burn-in.3 3 Também conhecido como warm-up e é referente ao período que as cadeias ainda não convergiram para a região de interesse e estão espalhadas em uma área maior do espaço de parâmetros. As cadeias foram iniciadas em ‘bolas’ gaussianas ao redor dos valores minimizados pelo algoritmo de “Simplex downhill” descrito na Seção6 6. Modelagem da Curva de Luz e Determinação dos Parâmetros de Trânsito Para determinar os coeficientes de escurecimento de limbo é necessário primeiro um conhecimento preciso dos parâmetros de primeira ordem no modelo do trânsito, principalmente do período orbital (Porb ) e tempo do centro de trânsito (t0). Claramente, os valores imprecisos de Porb e t0 produzem impossibilidades na determinação correta dos parâmetros restantes na descrição do trânsito, pois tanto variações no Porb quanto no t0 deslocam horizontalmente (na coordenada temporal) o modelo da curva de luz do trânsito. Portanto, é essencial determinar Porb e t0 com a maior precisão possível. Como primeira etapa, usamos o modelo de trânsito para ajustar esses dois parâmetros. Neste passo optamos pelo uso do algoritmo downhill simplex, mantendo todos os outros parâmetros fixos. O algoritmo downhill simplex (ou algoritmo Nelder–Mead) é um método numérico de busca direta, e que é projetado para obter o mínimo ou máximo de uma função em um espaço multidimensional. É um método comumente aplicado a problemas de otimização não lineares, especialmente para os quais as derivadas podem não ser conhecidas e, portanto, o método de Newton-Raphson não pode ser utilizado. Depois de determinar Porb e t0, ajustamos os parâmetros restantes: a razão dos raios estrela-planeta, o ângulo de inclinação, o semi-eixo maior e os coeficientes de escurecimento de limbo quadráticos u1 e u2, mantendo dessa vez os valores de Porb e t0 já minimizados fixos, novamente utilizando o algoritmo downhill simplex. Essa minimização inicial é importante, pois usaremos esses resultados de ajuste como valores iniciais para nossa amostragem MCMC, como descrita a seguir na Seção7. . Na obtenção dos resultados, o diagnóstico de convergência foi dado a partir da análise básica de três pontos. Primeiro, o gráfico de traçados (trace plot), que mostra como as cadeias exploraram o espaço de parâmetros. Segundo, analisamos a fração de aceitação média e por último avaliamos o tempo de auto-correlação médio. Este tempo de auto-correlação médio é, de forma sucinta, o “tempo” (em medidas de iterações) que a cadeia leva pra convergir. Já a fração de aceitação média é a taxa que as propostas geradas pelo Stretch Move estão sendo aceitas. Para o nosso problema, o resultado do tempo de auto-correlação médio das cadeias foi de 83.569 iterações e a fração de aceitação média foi de cerca de 51%. Sendo assim, comparando nosso resultado com valores descritos na literatura, percebe-se que estamos dentro do aceitável, onde o tempo de auto-correlação médio é várias vezes menor que o número de iterações totais para se obter uma boa amostragem da distribuição a posteriori. Já para a fração de aceitação média, aceita-se algo entre 25% e 75% para a convergência.

Um subproduto do MCMC é uma estimativa do ajuste da curva de luz usando as medianas das cadeias. Este ajuste, é representado na Figura8. O modelo aqui apresentado pode ser aplicado na íntegra para dados públicos de outros alvos observados pelo satélite Kepler. O resultado principal do presente trabalho é o processo de implementação que está representado na Figura9. Nesta figura podemos analisar as covariâncias entre os parâmetros e a amostragem de distribuição posteriori. A partir desta distribuição, no contexto Bayesiano, podemos fazer estimativas de erro dado por quantis de 1σ. Esses resultados fornecem as seguintes estimativas dos parâmetros de trânsito para Kepler-425b e intervalo 1σ: a/r*=11.61-0.0749+0.0796; rp/r*=0.1144-0.0007154+0.0006754; i=87.03-0.06839+0.07456; u1=0.5595-0.07034+0.07103; u2=0.1047-0.1177+0.1207.

Figura 8:
Painel superior: Preto: Curva de luz de KIC 5357901 dobrado na fase (diagrama da magnitude pela fase). Vermelho: modelo mais provável resultante do MCMC. Painel inferior: resíduo entre o modelo e os dados.
Figura 9:
Gráfico corner para KIC 5357901 representando as distribuições a posteriori conjuntas entre diferentes parâmetros, assim com suas distribuições a posteriori marginais. Os valores de q1 e q2 são associados com os coeficientes de escurecimento de limbo u1 e u2. As linhas tracejadas representam os quantis 0.16, 0.5 e 0.84, respectivamente. Gráfico feito utilizando a biblioteca corner.py [27[27] D. Foreman-Mackey, The Journal of Open Source Software 1, 24 (2016).].

Agradecimentos

G.Z.D. agradece ao Programa de Educação Tutorial (PET) pelo apoio financeiro e estrutural em todas as fases desta pesquisa. Agradece também ao Grupo de Estrutura e Evolução Estelar-DF/UFRN4 4 http://ge3.dfte.ufrn.br/ por todo o apoio acadêmico em todas as fases deste trabalho, em especial para Leandro de Almeida (Ted) e Eduardo Nunes. J. D. N. Jr é financiado parcialmente pelo CNPq com a bolsa PQ1. Tratamento, redução e distribuição dos dados do satélite Kepler é gentilmente realizado pela NASA, através do Mikulski Archive for Space Telescopes (MAST).

Apêndice

A. Implementação do MCMC

Existem diversas variações do MCMC que podem ser aplicadas dependendo das condições de contorno de cada contexto de dados e modelos. Um dos mais utilizados na literatura é o algorítmo de Metropolis-Hastings [28[28] W.K. Hastings, Biometrika 57, 97 (1970).], [29[29] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller e E. Teller, Journal of Chemical Physics 21, 1087 (1953).] baseado em um passeio aleatório, que consiste em gerar uma proposta de parâmetros baseada em uma distribuição normal centrada nos parâmetros atuais. Considerando distribuições a priori uniforme, e caso essa proposta tenha maior verossimilhança que a verosimilhança dos parâmetros atuais (ou seja, Lj/Lj-1>1) então ela será aceita, caso contrário, a proposta ainda pode ser aceita desde que Lj/Lj-1>q onde q é um número amostrado de uma distribuição uniforme entre 0 e 1, e Lj representa a verossimilhança na iteração j. Dessa forma o algorítmo consegue evitar ficar preso em mínimos locais e ter maior liberdade para explorar o espaço de parâmetros. Entretanto, nesse trabalho, vamos usar uma variação proposta por Goodman & Weare (2010) [30[30] J. Goodman e J. Weare, Communications in Applied Mathematics and Computational Science 5, 65 (2010).] conhecida como Stretch Move ou ainda Stretch Move paralelo, que mantém a maneira que o algorítmo de Metropolis-Hastings é proposto, porém evita mínimos locais enquanto altera a forma que é feita a proposta de parâmetros. Essa variação depende de implementarmos cadeias de Markov paralelas simultâneas, denominadas de walkers5 5 No contexto deste artigo, “walkers” são as caminhadas na busca de soluções explorando o espaço de parâmetros. Iremos manter aqui a terminologia inglesa. . Os walkers evoluem de forma que a direção da distribuição proposta para um determinado walker depende da posição de outros walkers do conjunto, como explícito no algorítimo (1). A magnitude do passo proposto depende de uma distribuição g(z), como representado pela Equação (A1), onde a é uma constante que pode ser ajustada para cada aplicação específica, mas o valor padrão usado na maioria dos casos [30[30] J. Goodman e J. Weare, Communications in Applied Mathematics and Computational Science 5, 65 (2010).] é a=2, e é esse o valor que usamos nesse trabalho. A Figura 10 traz uma representação visual de como esse algorítimo funciona.

(A1) g ( z ) { 1 z , se z [ 1 a , a ] 0 , caso contrário.

Para K walkers no conjundo S={Xk}, a proposta para um “walker” k se baseia em uma “walker” do conjunto complementar Sk={Xj,jk}.

Algoritmo 1:
Stretch Move
Figure 10
À esquerda, representação visual da direção proposta pelo “Stretch Move”. À direita, magnitudes máximas (Laranja) e mínimas (Amarelo) da proposta com a=2.

Nesse trabalho, especificamente com relação à implementação, contamos com a biblioteca EMCEE [31[31] D. Foreman-Mackey, D.W. Hogg, D. Lang e J. Goodman, Publications of the Astronomical Society of the Pacific 125, 306 (2013).] para auxílio no processo de amostragem por “Stretch Move”. Essa biblioteca oferece uma implementação completa em Python, tanto do amostrador quanto do MCMC e é amplamente usada nessa classe de problemas descritos na literatura.

Referências

  • [1]
    M. Mayor e D. Queloz, Nature 378, 355 (1995).
  • [2]
    D.G. Koch, W.J. Borucki, G. Basri, N.M. Batalha, T.M. Brown, D. Caldwell, J. Christensen-Dalsgaard, W.D. Cochran, E. DeVore, E.W. Dunham et al., The Astrophysical Journal 713, L79 (2010).
  • [3]
    M. Auvergne, P. Bodin, L. Boisnard, J.T. Buey, S. Chaintreuil, G. Epstein, M. Jouret, T. Lam-Trong, P. Levacher, A. Magnan et al., Astronomy and Astrophysics 506, 411 (2009).
  • [4]
    R.K. Kopparapu, R. Ramirez, J.F. Kasting, V. Eymet, T.D. Robinson, S. Mahadevan, R.C. Terrien, S. Domagal-Goldman, V. Meadows e R. Deshpande, The Astrophysical Journal 765, 131 (2013).
  • [5]
    S.B. Howell, C. Sobeck, M. Haas, M. Still, T. Barclay, F. Mullally, J. Troeltzsch, S. Aigrain, S.T. Bryson, D. Caldwell et al., Publications of the Astronomical Society of the Pacific 126, 398 (2014).
  • [6]
    G.R. Ricker, J.N. Winn, R. Vanderspek, D.W. Latham, G.Á. Bakos, J.L. Bean, Z.K. Berta-Thompson, T.M. Brown, L. Buchhave, N.R. Butler et al., Journal of Astronomical Telescopes, Instruments, and Systems 1, 014003 (2015).
  • [7]
    H. Rauer, C. Catala, C. Aerts, T. Appourchaux, W. Benz, A. Brandeker, J. Christensen-Dalsgaard, M. Deleuil, L. Gizon, M.J. Goupil, Experimental Astronomy 38, 249 (2014).
  • [8]
    J.P. Gardner, J.C. Mather, M. Clampin, R. Doyon, M.A. Greenhouse, H.B. Hammel, J.B. Hutchings, P. Jakobsen, S.J. Lilly e K.S. Long, Space Science Reviews 123, 485 (2006).
  • [9]
    S. Csizmadia, T. Pasternacki, C. Dreyer, J. Cabrera, A. Erikson e H. Rauer, Astronomy and Astrophysics 549, A9 (2013).
  • [10]
    D. Kipping e G. Bakos, The Astrophysical Journal 730, 50 (2011).
  • [11]
    H.M. Müller, K.F. Huber, S. Czesla, U. Wolter e J.H.M.M. Schmitt, Astronomy and Astrophysics 560, A112 (2013).
  • [12]
    H.A. Knutson, D. Charbonneau, R.W. Noyes, T.M. Brown e R.L. Gilliland, The Astrophysical Journal 655, 564 (2007).
  • [13]
    A.A. Souza e A. Valio, Revista Brasileira de Ensino de Física 41, e20180323 (2019).
  • [14]
    E.A. Milne, Monthly Notices of the Royal Astronomical Society 81, 361 (1921).
  • [15]
    A. Manduca, R.A. Bell e B. Gustafsson, Astronomy and Astrophysics 61, 809 (1977).
  • [16]
    J. Diaz-Cordoves e A. Gimenez, Astronomy and Astrophysics 259, 227 (1992).
  • [17]
    D.A. Klinglesmith e S. Sobieski, The Astronomical Journal 75, 175 (1970).
  • [18]
    P.A.D. Almeida e J. Gregorio-Hetem, Revista Brasileira de Ensino de Física 44, e202100405 (2022).
  • [19]
    H.M. Muller, Limb-darkening Measurements on Exoplanet Host Stars and the Sun Dissertação de Mestrado, Universität Hamburg, Hamburg (2015).
  • [20]
    K. Mandel e E. Agol, The Astrophysical Journal 580, L171 (2002).
  • [21]
    S. Seager e G. Mallén-Ornelas, Scientific Frontiers in Research on Extrasolar Planets 294, 419 (2003).
  • [22]
    L. Kreidberg, Publications of the Astronomical Society of the Pacific 127, 1161 (2015).
  • [23]
    N. Espinoza e A. Jordán, Monthly Notices of the Royal Astronomical Society 457, 3573 (2016).
  • [24]
    D.M. Kipping, Monthly Notices of the Royal Astronomical Society 435, 2152 (2013).
  • [25]
    G. Hébrard, A. Santerne, G. Montagnier, G. Bruno, M. Deleuil, M. Havel, J.M. Almenara, C. Damiani, S.C.C. Barros, A.S. Bonomo et al., Astronomy and Astrophysics 572, A93 (2014).
  • [26]
    N.R. Lomb, Astrophysics and Space Science 39, 447 (1976).
  • [27]
    D. Foreman-Mackey, The Journal of Open Source Software 1, 24 (2016).
  • [28]
    W.K. Hastings, Biometrika 57, 97 (1970).
  • [29]
    N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller e E. Teller, Journal of Chemical Physics 21, 1087 (1953).
  • [30]
    J. Goodman e J. Weare, Communications in Applied Mathematics and Computational Science 5, 65 (2010).
  • [31]
    D. Foreman-Mackey, D.W. Hogg, D. Lang e J. Goodman, Publications of the Astronomical Society of the Pacific 125, 306 (2013).

Datas de Publicação

  • Publicação nesta coleção
    16 Set 2022
  • Data do Fascículo
    2022

Histórico

  • Recebido
    16 Jun 2022
  • Revisado
    03 Ago 2022
  • Aceito
    04 Ago 2022
Sociedade Brasileira de Física Caixa Postal 66328, 05389-970 São Paulo SP - Brazil - São Paulo - SP - Brazil
E-mail: marcio@sbfisica.org.br