EDP Sciences
Free Access
Issue
A&A
Volume 595, November 2016
Article Number A89
Number of page(s) 14
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201527699
Published online 03 November 2016

© ESO 2016

1. Introduction

Stellar activity is one of the main sources of uncertainty for planet detection and characterization. It causes the emergence of both cooler (therefore darker) and hotter (brighter) than average regions on the stellar photosphere. Such dark and bright spots, or activity features, cross the visible stellar disk as the star rotates, and therefore modulate the amount of flux emitted by the star. These spots are distributed in groups, and vary in size, temperature, and position on the stellar disk along an activity cycle.

In transit photometry, activity features can induce systematic errors in the determination of the planetary parameters. Czesla et al. (2009) showed that transit normalization is affected by non-occulted spots on the stellar disk, which leads to an overestimation of the transit depth. They and Silva-Valio et al. (2010) discussed how spots occulted during a transit act in the opposite way, producing an underestimation of the same parameter. Léger et al. (2009) showed that stellar activity can lead to an underestimation of the stellar density. Csizmadia et al. (2013) studied the effect of starspots on the estimate of limb-darkening coefficients. Alonso et al. (2009), Barros et al. (2013), and Oshagh et al. (2013b) showed that stellar activity can induce apparent transit timing variations (TTV), and can introduce errors in the determination of the transit duration as well.

Several approaches have been tried to disentangle the signal produced by a planet from the one which comes from activity features. The main attempts focus on the data reduction; that is, data are either corrected for the identified activity signatures, or the activity features are masked.

To correct for the errors in the derived planet-to-star radius ratio, Czesla et al. (2009) proposed adopting a different transit normalization technique to the standard one. The standard normalization consists in dividing each transit profile by a low-order polynomial fitted to the flux adjacent to both sides of the transit. With the normalization of Czesla et al. (2009), the out-of-transit flux modulations are taken into account. Moreover, they assumed that dark spots are dominant over bright ones, and proposed using the lower envelope of the deepest transits in order to recover a transit profile closer to the true one. With this approach on CoRoT-2 b, they found a ≃ 3% larger planet-to-star radius ratio than the one reported in the discovery paper (Alonso et al. 2008) where a standard approach was used.

Various methods have been developed in order to model stellar activity features. Some are based on analytic models (e.g. Budding 1977; Dorren 1987; Kipping 2012; Montalto et al. 2014, and references therein), while others make use of numerical techniques (e.g. Oshagh et al. 2013a; Dumusque et al. 2014, and references therein). The fitting techniques also differ, going from the division of the stellar surface in segments on which a χ2 minimization is performed (e.g. Huber et al. 2010), to maximum entropy regularization (Rodonò et al. 1995; Collier Cameron 1997; Lanza et al. 1998), and to modified Markov chain Monte Carlo (MCMC) algorithms to sample the spot parameter space (Tregloan-Reed et al. 2015).

In this work, we present a method to take into account the imprint of activity features on the transit parameters. We update an existing code for activity features modelling, and implement it into our Bayesian modelling package PASTIS. Then, we use it to model the light curve of CoRoT-2. The clearly visible activity pattern of the star of this system, both outside and inside the transits, has made it a widely studied case in the literature. Our goal is to take advantage of the modelling in order to more consistently determine the transit parameters.

The paper is structured as follows: in Sect. 2 we present our method; in Sect. 3 we describe the test case and the application of our method; in Sects. 4 and 5 the results are presented and discussed; and we conclude in Sect. 6.

2. Method

The modelling of activity features is known to be an ill-posed inversion problem because the involved parameters are highly degenerate. To reach convergence in the fit, some of the parameters are kept fixed, or only a part of the data is fitted (e.g. Lanza et al. 2009; Huber et al. 2010, and references therein). The computation time required by the modelling is another serious inconvenience. Numerical methods define a high-resolution two-dimensional grid of the stellar surface or of its projection onto a plane, and numerically integrate over two coordinates to obtain the light curve. They can deal with a large set of activity features configurations: in principle, they allow for the modelling of any shape for the features, flux profile, and limb-darkening law. The integration, however, is very expensive in terms of computation and time.

Analytic models, on the other hand, use analytic formulae to derive a synthetic light curve. They can be order-of-magnitudes faster, but usually require restricting constraints on the parameters to simplify the equations. The main restrictions are that they require simple circular shapes for the activity features, they assume a small size for the spots relative to the stellar surface, and they do not consider the overlapping of the features or of the features with a planet (Kipping 2012).

Semi-analytic models were employed in some cases in order to overcome these limitations (Béky et al. 2014). Alternatively, analytic models which allow for the activity features to overlap, as well as transit modelling, have recently been presented and implemented in freely available codes (e.g. KSint, Montalto et al. 2014). Their rapidity of execution allows them to be implemented in MCMC algorithms.

MCMC methods have already been proven effective in order to find best-fit values, uncertainties, correlations, and degeneracies for the photometric spot modelling problem (Croll 2006). Hence, we decided to exploit the additional transit modelling capabilities of an analytic starspot modelling code and the efficiency of MCMC simulations. We implemented the code KSint into the MCMC algorithm used by PASTIS (Díaz et al. 2014; Santerne et al. 2015). Hereafter, the combination of the codes will be referred to as KSint + PASTIS.

KSint models a light curve containing both planetary transits and activity features. The transits, modelled with the formalism of Pál (2012), are characterized by the planet-to-star radius ratio kr, orbital period Porb, orbital inclination ip1, eccentricity e, planet argument of pericentre ω, and mean anomaly M. The star is assigned an inclination angle i, a rotation period P, a density ρ, and quadratic limb-darkening coefficients ua and ub. The activity features are characterized by the same limb-darkening law as the star.

The features are assumed to be spherical caps. Each of them is described by four parameters: longitude λ, latitude φ, angular size α, and contrast f. The code models dark and bright spots. It should be noted that bright spots are not faculae, as their contrast does not change from the centre to the limb of the star, as happens for faculae. Bright spots, instead, can be modelled with the same limb-darkening law used for dark spots.

In the original version of the code, the time evolution of the features, which has been observed for many stars, is not included. The addition of parameters for feature evolution worsens the problems of correlation, degeneracy, and non-uniqueness of the solution. On the other hand, it allows longer parts of the light curve to be fitted.

thumbnail Fig. 1

Model light curve from KSint plotted over the data. The eight segments of the fit are divided by colour. The residuals are shown in the lower panels; the error bars are not shown for clarity. The larger amplitude of the residuals in correspondence with the transits is due to the full resolution kept for the transits. The out-of-transit binning is of 2016 s, and inside the transits 32 s.

Open with DEXTER

We therefore introduced a simple law for activity features evolution in KSint. Following the prescription of Kipping (2012), we used a linear variation of the angular size. The size parameter α was translated into the maximum size reached by a feature during its evolution, αmax. Then, four parameters were added to the description of every feature: 1) the time at which the maximum size is reached, tmax; 2) the time during which the feature keeps its maximum size, tlife; 3) the time of growth from zero to maximum size, ; and 4) the time of decay from maximum to zero size, .

The features of our model should be considered as representative of groups of features, rather than features taken individually. This allows large sizes and lifetimes to be used without losing physical meaning. PASTIS uses three other parameters to model a light curve: a normalized flux offset, a contamination term, and an instrumental jitter.

3. Application to CoRoT-2

3.1. Presentation of the system: previous studies

CoRoT-2 A is a G7V-type star observed during the LRc01 run of the CoRoT space telescope. It hosts the hot Jupiter CoRoT-2 b (Alonso et al. 2008), which has mass 3.31 ± 0.16MJ and radius 1.465 ± 0.029RJ. The orbit of the planet has a period of 1.74 days, and is almost aligned with the stellar equator (Bouchy et al. 2008). Its radius is about 0.3 RJ larger than expected for an irradiated hydrogen-helium planet of this mass. Models strive to explain a longer contraction time during the evolution of the planet, and allow for scenarios with a 30 to 40 Myr old, pre-main sequence host star (Guillot & Havel 2011). A detailed study of several age indicators favours a main sequence star with an age between 100 and 300 Myr, while the inflated radius of the planet can be explained by a transient tidal circularization and a corresponding tidal heating in the interior of the planet (Gillon et al. 2010). However, the system might actually be much older, and the radius anomaly of the planet can be explained by a stellar companion gravitationally bound to the system (Schröter et al. 2011). If the star were effectively older than observations suggest, its high level of activity may have been enhanced by its tidal interaction with the planetary companion (Poppenhaeger & Wolk 2014). The main characteristics of the system are listed in Table A.1.

The light curve, shown in Fig. 1, indicates that a varying fraction of the stellar surface – up to a few tens of percent – is covered by activity features (Lanza et al. 2009; Huber et al. 2010; Silva-Valio et al. 2010). Fits with a few and with several features have been performed on CoRoT-2, independently from the study of the planet.

Wolter et al. (2009) worked on a single transit to model an occulted spot. They constrained the size of the spot (between 4.5 and 10.5°) and its longitude with a precision of about . Lanza et al. (2009) fitted segments of light curve no longer than 3.2 days because the evolution of activity features was not allowed for in their models. They removed the transits and used both a three-spot model and a maximum-entropy regularization method, finding a stellar rotation period of ~4.5 days. They recovered two active longitudes on different hemispheres, and estimated the relative differential rotation of the star to be lower than ~1%. They also measured cyclic oscillations of the total area covered by active regions with a period of ≃ 29 days, and found that the contribution of faculae to the optical flux variations is significantly lower than in the present Sun.

Fröhlich et al. (2009) presented a Bayesian analytic spot model of the out-of-transit part of the light curve. They modelled three long-lived spots and recovered the same active longitudes as those found by Lanza et al. (2009). However, their model determines a differential rotation rate which exceeds by one order of magnitude the value found by these authors. Fröhlich et al. state that their assumption on spots longevity, not adopted by Lanza et al., is necessary to recover such a high value of differential rotation.

Huber et al. (2010) modelled the light curve, including transits, over two stellar rotations. They fitted the activity features on the whole light curve but not the transit parameters, and modelled the transits with a combination of the values published by Alonso et al. (2008) and Czesla et al. (2009). They managed to reproduce the photometric signal with a low-resolution surface model of 36 longitude strips, and found that the belt occulted by the planet (close to the stellar equator) is ~6% darker than the average remaining surface. This study was extended to the whole light curve by Huber et al. (2010), whose results are in agreement with the previous study. Significant indications of stellar differential rotation were found.

Silva-Valio et al. (2010) analysed the occulted features inside all transits, but excluded the out-of-transit part of the light curve from the fit. Up to nine spots per transit were needed for the fit. These authors found the size of spots to be between 0.2 and 0.7 planetary radii, and a spot coverage of 1020% in the belt transited by the planet. For the spots, they found contrasts between 0.3 and 0.8. Finally, they found a planet-to-star radius ratio of the deepest transit (assumed to be less affected by spots) of 0.172, i.e. 3% larger than the value found by Alonso et al. (2008).

Table 1

Data sets prepared for the modellings.

The occulted features were modelled in a spot map for every transit by Silva-Valio & Lanza (2011). For this map, 392 spots were used. The evolution timescale of the coverage of active regions in the transit chord was found to be between 9 and 53 days.

All these approaches are complementary and rely on specific assumptions to assess a different aspect of the light curve. However, none takes full advantage of the information encoded in the out-of-transit part of the light curve in order to properly correct the imprint of stellar activity on the transit profile. This was the motivation for our study: with our model, we aim at a more thorough exploration of the impact of the active regions on the transit parameters.

3.2. Data reduction

We used the light curve processed by the CoRoT N2 pipeline2. Only the part sampled every 32 s was used (~145 over 152 days of observation) to take advantage of the full resolution inside transits. The points classified as poor in quality, including those related to the South Atlantic Anomaly, were discarded. The data were first 3σ-clipped and the transits were identified with the ephemerides given by Alonso et al. (2008).

Five data sets were prepared, as summarized in Table 1.

Data set .

This data set was prepared as a reference to be compared to a standard transit analysis. The flux was divided by its median value. A second-order polynomial was fitted to the flux adjacent to the transits. For this, a window as long as a transit duration was considered on the sides of each transit. The flux was then divided by the resulting polynomial.

Data set .

This data set was prepared for comparison with a standard transit fit, where non-occulted features are taken into account in the data reduction phase, as done by Czesla et al. (2009). The normalization was performed following the prescription of these authors. For every bin i, we calculated the normalized flux zi as (1)where fi is the observed flux in bin i, ni is the value of the second-order polynomial fitted to the out-of-transit local continuum in the same bin, and p is the largest flux value in the light curve. This last parameter is assumed to be the least affected by dark spots.

Data set .

This data set was used for the fit of activity features (the label stands for “spots”) while fixing the transit parameters. To reduce the computation time, we followed Huber et al. (2010) and sampled the out-of-transit flux every 2016 s (33.6 min). This corresponds to about one-third of the orbital period of the CoRoT satellite (Auvergne et al. 2009). Given the stellar rotation period, this corresponds to a resolution of ~1.9°/bin, which is sufficient not to significantly affect the out-of-transit flux measurements. The transits were sampled every 160 s, i.e. with a resolution of ~0.1°/bin. The transits are ~0.09 days long, which means about 49 points were left for each transit. This sampling keeps enough information to force some features to be modelled inside the transits. The uncertainty of the binned data was calculated as the standard deviation of the points in each bin, divided by the square root of the number of points in the binning window. The average relative uncertainty on this data set is 3.9 × 10-4. After resampling, the flux was divided by its maximum value and no transit normalization was performed.

Data set .

This data set was prepared for active regions and transit fitting. It was prepared as data set , but full resolution was kept inside the transits. The average relative standard deviation on this data set is of 1.12 × 10-3.

Data set .

The out-of-transit model fitting data set was used to normalize the transits. The preparation of this data set is described in more detail in Sect. 3.3.

3.3. Modelling approaches

We performed two reference fits with standard approaches. Then, we performed three fits including activity modelling.

3.3.1. Reference fits

Fit : activity neglected.

Fit was performed on data set as a standard transit fit. We used a modified version of the EBOP code3 (Nelson & Davis 1972; Etzel 1981; Popper & Etzel 1981) included in the MCMC functionality of PASTIS. The priors for the free parameters are listed in Table 2. Because of the biases in our current knowledge of the radius distribution of planets, the incompleteness of the transit surveys, and other systematic effects, we followed Díaz et al. (2014) and used a Jeffreys prior for the planet-to-star radius ratio kr. We used a Jeffreys prior for the orbital separation-to-stellar radius ratio a/R, which is used by EBOP in place of the stellar density. A normal prior was also used for the orbital period Porb. To impose an isotropic distribution on orbit orientations, we set a sine prior for the orbit inclination ip. Uniform priors were used for the linear and the quadratic limb-darkening coefficients ua and ub. Following Alonso et al. (2008), we fixed e = 0, and consequently also ω = 0. The flux offset and the jitter were fitted on the whole light curve.

Table 2

Prior distributions used in the combined fit with KSint + PASTIS.

Alonso et al. (2008) used a 5.6 ± 0.3% contamination rate for CoRoT-2, as calculated before the CoRoT launch using a set of generic PSFs (contamination level 0, Llebaria & Guterman 2006). Gardes et al. (2012) updated this value to 8.81 ± 0.89% using a more realistic estimate of the CoRoT PSF (contamination level 1). We adopted this last value, and fitted the contamination rate using this prior. A larger contamination rate increases the variability amplitude of the star and therefore of the measured transit depth. We expect this latter to increase on average by (1−0.056)/(1−0.081) = 2.7%.

The same MCMC simulation was also performed with KSint instead of EBOP in order to verify its correct behaviour for a standard transit-only fit.

Ten chains were run for every MCMC set. The chains were then thinned according to their correlation length and merged into a single chain, giving the credible intervals for the parameters of that segment. To be considered robust, a merged chain was required to consist of at least a thousand uncorrelated points.

Fit : activity taken into account during normalization.

For this fit we used the same method as for fit , but on data set .

3.3.2. Fits that include activity

Because of the high level of degeneracy and correlation of the activity features parameters in the starspots modelling problem, different Markov chains do not converge to the same values for the parameters. We therefore divided the fit into two steps. First, only the parameters of the active regions were fitted with a MCMC run. This resulted in a local maximum-likelihood solution. Starting from this solution, we then used another MCMC run to sample the posterior distribution of a subset of key activity features parameters together with the transit parameters. Eventually, we performed a transit fit by taking advantage of the maximum-likelihood solution for the activity features.

Fit  of  active  regions  (Fit  ).

This fit takes into account the activity features only, and prepares the fit of the transits. The light curve modelling was carried out on the data set . The transit parameters (Porb,ρ,ip,kr,e,ω,ua,ub) were fixed to the values of Alonso et al. (2008). A uniform prior was imposed for the mean anomaly M. The priors indicated in Table 2 were used. The quasi-perpendicularity between the sky-projected stellar spin axis and the planetary orbit found by Bouchy et al. (2008) allowed us to model the system as if these two were perpendicular. Alonso et al. (2008) found the planet inclination to be of 87.84° with respect to the plane of the sky. We therefore fixed the stellar inclination to the same value with respect to the line of sight. A uniform prior centred on the value found by Lanza et al. (2009) was set for the stellar period P. Uniform priors were used for all the parameters of the activity features. For each feature, αmax was limited between 0 and 30°. The largest size was found to be sufficient for the modelling.

The analysis of Lanza et al. (2009) showed that bright regions have a minor impact on the light curve of CoRoT-2 compared to dark features. However, these authors obtained their result after removing the transits from the data set. Instead, because our data set includes transits, we decided to model both dark and bright spots. We set a prior for the contrast f between 0.3 and 1.3. These values represent a wider interval than the one defined by the typical sunspot bolometric contrast of 0.67 (Sofia et al. 1982; Lanza et al. 2004) and the bolometric facular contrast 1.115 adopted by Foukal et al. (1991) and Lanza et al. (2004). Latitudes were left free between −90° and 90° in order to allow both non-occulted and occulted spots to be modelled.

For each feature, a uniform prior was set for tmax. Values of tmax external to the time limits of the light curve were allowed in order to induce the fit to exclude possible features in excess. The MCMC could do this by adjusting tmax, , and , to which uniform priors were assigned, too. The contamination factor was fixed to the value found by Gardes et al. (2012; 8.81%). This parameter was fixed because active regions can be considered a contaminant factor (Csizmadia et al. 2013). A uniform prior was set for the flux relative offset and the jitter.

The optimal number of active regions was determined by trials. We increased it until we obtained residuals with a normally distributed dispersion centred on zero and with a width comparable to the photometric data dispersion. Our model needs several tens of active regions to model the entire light curve, in agreement with the results of Silva-Valio et al. (2010) and Silva-Valio & Lanza (2011). This implies hundreds of free parameters, which our computer cluster is not able to handle. Therefore, the light curve was divided into shorter parts which need the modelling of fewer features. The initial and ending time of these segments are indicated in Table A.2. The segments have a duration of ~15−25 days (four to six stellar rotations), and can be fitted with six to nine evolving features. Longer segments tend to produce worse fits. In Fig. 1, the segments are highlighted in colour. It can be seen how the segments are related to the different phases of activity in the light curve. The brightness variations grow in amplitude, reach a maximum, and shrink again. The duration of the segments is consistent with the lifetime of individual spots and active regions found by Lanza et al. (2009), i.e. between 20 and 50 days.

To connect the solutions of consecutive segments, the evolving features with non-zero size at the end of a given segment were kept for the initial values of the features of the next segment. Up to eight features were kept from one segment to the next. For such features, we used the same latitude φ, and constrained the longitude λ around the value of the solution in the previous segment. For this parameter, we chose a normal prior with standard deviation equal to 10°, and constrained c in order to force the feature to remain darker or brighter than the stellar surface. The parameter αmax was again set as a free parameter with a uniform prior, starting from the maximum-likelihood value of the previous segment. The evolution parameters were left free in order to allow for the feature size to grow after a decaying phase and vice versa. Other activity features were added to those kept to connect the segments using the same priors described in Table 2. In this way, the appearance of new features along the light curve were modelled. The total number of features per segment (fixed and free) is between six and nine.

The flux offset was fitted separately in every segment, starting the chains of each segment from the maximum-likelihood value of the previous one. In this way, a possible photometric long-term trend was taken into account.

Ten chains of up to 3 × 105 iterations were run for each segment. Most of the chains of each segment were observed to reach the end of their burn-in phase with this number of iterations. We considered a chain to have reached a steady state if, after its burn-in phase, its likelihood did not show significant increments for a few thousand steps. We then extracted the maximum-likelihood solution among the chains satisfying this criterion.

The probability of finding the global maximum of the likelihood function increases with the number of chains run. However, we note that our method does not allow us to identify the global solution of the starspot inversion problem. With this approach we obtained a fit of the activity features parameters suitable for the fit of the transit parameters. Also, the maximum-likelihood solution allowed us to explore the active regions configuration on the stellar surface. This is discussed in greater detail in the following sections.

Spot-transit fit (fit ).

We carried out a simultaneous fit of the transit parameters and of the activity features both in-transit and out-of-transit. This fit was performed on data set , divided into segments as for fit . This time the transit parameters were set as free parameters, while the parameters of the active regions were fixed to the maximum-likelihood solution of fit for each segment. An exception was made for the size αmax and the longitude λ of the activity features. The parameter αmax affects primarily the planet-to-star radius ratio. Its fit was started from the value of the best solution of fit . The longitude affects the position of a feature in the transit profile, with consequences for the transit duration, and therefore the limb-darkening coefficients, ip, and ρ. We used a normal prior with standard deviation equal to , centred on the best likelihood value.

The priors of the transit parameters are the same as for fit and . Two parameters were set differently. Instead of a/R, KSint uses ρ, for which we used a Gaussian prior centred on the Alonso et al. (2008) value and a standard deviation 0.35 ρ. Normal priors were set for the limb-darkening coefficients, using the results of Alonso et al. (2008).

In this phase, ten chains from 1.5 to 2 × 105 steps were employed to reach convergence. For each segment, the chains were then thinned according to their correlation length and merged into a single one. The results and the uncertainties for each segment were obtained from the respective merged chains, as for fit and .

This fit resulted in slightly scattered transit parameters among the various segments, as discussed in Sect. 4.2. The posterior distributions allowed us to explore the impact of bright spots on the deepest transit and the correlations between activity features and transit parameters. This will be discussed in Sect. 5.

Model-based normalization (fit ).

The fit consists of a standard transit fit where transit normalization is performed via the modelling of the out-of-transit features. The maximum-likelihood solutions obtained from fit were re-computed without the planet. They were merged into a single light curve and used to normalize the transits, obtaining data set (Sect. 3.2). Equation (1) was used, with ni indicating the model and p the flux offset value. For each segment, this parameter was fixed to the maximum-likelihood value obtained with fit . Its uncertainty, obtained with fit , was quadratically added to the standard deviation of the flux. In this way, the uncertainty due to the normalization technique was evaluated.

Once the transits were normalized, a standard transit fit was performed with EBOP. This code was used because of its higher computation speed compared to KSint. The priors are the same as fit and . The contamination was fitted as well because in-transit activity features were no longer modelled.

4. Results

Figure 1 presents the best model of fit plotted over the light curve. The segments are divided by colour. The parameters of the active regions (yielded by fit ) and the transit parameters (fit ) are discussed separately.

4.1. Parameters of the active regions

The maximum-likelihood solutions of fit allowed us to recover some properties of the starspot surface coverage for each segment of the light curve. We do not expect our solutions to be in full agreement with those presented by other authors because our fit 1) includes the fit of activity features inside transits, which provides a strong constraint on their longitudes; 2) uses both dark and bright spots, which obey the same limb-darkening law; and 3) models the evolution of the features. In alternative modellings, instead, the authors drop one of these points.

We computed the effective coverage factor of the stellar surface as a function of the segment of light curve, separately for dark and bright spots. For each segment, this is defined as (2)where the sum is run over all the features i with maximum area αmax and contrast f. As only the maximum-likelihood solutions were used, no error estimate is possible. A bright spot has f> 1, hence yields a negative .

In Fig. 2, the absolute value for dark (in blue) and bright (in red) spots is plotted. It can be noticed that a larger for dark spots is accompanied by a somewhat larger for bright spots. Lanza et al. (2009) found a faculae-to-dark spots surface ratio between 1 and 2.5 in active regions. We chose not to directly compare our results with theirs because, unlike these authors, we used the contrast of the activity features as a free parameter.

Discontinuities of can be observed on the separations between segments (vertical dashed lines). These is due to the lack of constraints on the features evolution parameters tlife, tmax, , and between consecutive segments.

thumbnail Fig. 2

Absolute value of the effective coverage factor as a function of time. The vertical dashed lines mark the separations between segments.

Open with DEXTER

Figure 3 represents the activity features for every segment as a function of their longitude (x-axis) and time (y-axis). No latitudinal information is reported. Filled circles represent dark spots, empty circles represent bright ones. Different colours are used for features in different segments. The sizes of the circles indicate the αmax of the activity features. The relative size compared to the stellar surface was increased for the purpose of illustration. As the maximum size is represented and not its evolution along the light curve, the plot gives an upper-limit coverage of the surface for each segment. Each feature is placed at its respective segment’s midtime for clarity.

Table 3

Transit parameters with their 68.3% credible intervals.

This plot suggests that the activity features are not uniformly distributed along the longitudes. Also, a longitudinal migration towards higher longitudes can be seen for the activity features with λ between 0 and ~50° and those between 200 and 360°. Tentative patterns of migration are shown with diagonal dashed lines. Features between 100 and 200°, instead, do not show the same migration pattern, but keep the rotational period assigned to the star. We argue that assigning different rotation periods to the features, as Fröhlich et al. (2009) did, would allow longer segments of the light curve to be modelled.

The inhomogeneous distribution of the features and their size along the stellar longitudes, shown in Fig. 3, suggest the presence of differently active longitudes. By using a continuous distribution of active regions with fixed contrast or only three dark spots on the stellar surface, respectively, Lanza et al. (2009) and Fröhlich et al. (2009) achieved a similar result. Huber et al. (2010) also confirmed such findings by using a different method. Given our different modelling of the evolution of the features, as well as the degeneracy between their number, their size, and their contrast, we chose not to compare in detail our results with theirs.

thumbnail Fig. 3

Distribution of the activity features as a function of their longitude and time. No latitudinal information is reported. Filled circles represent dark spots, empty circles represent bright ones. Different segments of the light curve are represented by circles of different colour. The relative size of the features compared to the stellar surface was increased for the sake of illustration. Tentative reconstruction of longitudinal migration are represented by diagonal dashed lines.

Open with DEXTER

Most of the features have αmax>. Moreover, the average size of the features does not change among different activity phases (increasing brightness variations, maximum, decrease of the brightness variations; see Sect. 3.3.2 and Fig. 1).

Bright spots are present in every segment. Moreover, in every segment, one or two features (either dark or bright spots) are found to cross the transit chord. As six to nine features are present in each segment, occulted features are a minority. Finally, no preferred values are found for the fitted evolution times tmax, tlife, , and .

4.2. Transit parameters

Table 3 reports the transit parameters of all fits. The results of fit were obtained as the average of all the results obtained on each segment. The values obtained for each segment are reported in Table A.2. For some segments, the results of fit are in poor agreement one to each other. This is likely due to statistical fluctuations introduced by the division of the light curve in segments, as will be discussed in more detail in Sect. 5.2. Because of this, the error bars on the average results of fit , in Table 3, are larger than on the other fits, whose uncertainties are similar to those found by Alonso et al. (2008). For the same reason, the average results of fit are only indicated for comparison to the other results, but are not discussed further. Figure 4 shows the results on each segment. The shaded regions correspond to the results of fit .

thumbnail Fig. 4

From top to bottom: values and uncertainties of kr, ip, ρ, Porb,ua, and ub for fit , indicated as a function of the segment of the light curve are in black. The results from fit with their 68.3% credible intervals are shaded in grey.

Open with DEXTER

Differences were found in the planet-to-star radius ratio kr among all fits carried out on the whole light curve ( and ). This highlights the dependence on the normalization technique which was used. As expected, the parameter kr is larger if a standard normalization is performed (fit , kr = 0.16917 ± 8.5 × 10-4). In fact, in this case, non-occulted spots are neglected and cause an overestimate of this parameter. On the other hand, the Czesla et al. (2009) normalization (fit ) produces a smaller kr (0.16632 ± ,8.0 × 10-4). This occurs because of the division of the transits profiles by the maximum flux value (p), which is assumed to be the least affected by activity. Without a model of activity features, no uncertainty is available for this value. This is the main limitation of this normalization.

thumbnail Fig. 5

Fit on the transits for every segment (from left to right, and top to bottom). The deepest transit, containing an occulted bright spot, is shown in yellow. A transit not affected by spots is in red.

Open with DEXTER

Fit models non-occulted activity features, therefore, it is not affected by the same issues as fits and . It yields a ~1.2σ smaller kr value than that obtained with fit , and 2.3σ larger than obtained with fit .

The kr on segment 5 is lower than those estimated from the other segments, and higher for ip and ρ. The large contribution of bright spots with respect to dark ones, observed in Fig. 2 for this segment, suggests that our correction of the transit profiles for bright spots can be refined.

The stellar densities ρ and the orbital inclinations ip are in agreement between fit and , while they are at less than 2σ agreement between fit and . Except for segments 5 and 8, the results of fit for ip and ρ are consistent with each other.

The orbital periods Porb found with all fits on the whole light curve are between ~1 and 3σ agreement. We note, however, that their measure is scattered among the segments of fit . This can be explained by statistical fluctuations due to number of transits – no more than about ten – in each segment. We exclude the effect of spurious activity-induced TTV such as those recovered by Alonso et al. (2009) in the CoRoT-2 light curve. These authors found a 7.45 d peak in the periodogram of the residuals of the transit midpoints and attributed it to activity features occulted by the planet. The amplitude of the peak they found is 20 s, while the variations in Fig. 4 are about one order of magnitude smaller.

The limb-darkening coefficients ua and ub were found to be compatible for all fits on the whole light curve. In fit , instead, ua results in a poorer agreement among the segments despite the imposed tight normal prior. An explanation can be the difficulty of PASTIS + KSint to recover the limb-darkening coefficients if only a few transits are available as too little information is available to disentangle the effect of activity features on the limb-darkening coefficients. A poorly fitted Porb could contribute to this problem as well. Another cause could be that the limb-darkening coefficients actually change as a function of the varying coverage of the stellar surface by active regions (Csizmadia et al. 2013). This possibility is discussed in Sect. 5.2.

thumbnail Fig. 6

Left: planet (black) and activity features (blue) configuration during the deepest transit, with bright spots allowed. The bright spots are crossed by the planetary disc. Centre: solution with three dark spots (red). Right: deepest transit fitted with the dark-bright spot configuration (blue) and the three-dark spots model (red, shifted). The residuals are shifted for clarity and use the same colour code.

Open with DEXTER

5. Discussion

The results of the modelling were used to explore the impact of the spots on the transit parameters.

5.1. Least distorted transit

According to Czesla et al. (2009), the deepest transits are less affected by occulted dark features than shallower ones. Therefore, their planet-to-star radius ratio should be closer to the true one. They thus interpolated a lower envelope to an average of the deepest transits of CoRoT-2 b. These were found to happen at the moments where the out-of-transit flux reaches the largest level. They fitted only kr and ip, and obtained kr = 0.172 ± 0.001 (i.e. a 3% larger kr than Alonso et al.) and ip = 87.7 ± 0.2°. Their method assumes the dominance of dark spots over faculae, as found by Lanza et al. (2009) by fitting the out-of-transit light curve.

With our approach, the assumption of the prevalence of dark spots inside transits could be checked. Even if we do not model faculae, we are still able to recover the need of features brighter than the stellar surface in the model. Such features would therefore increase the apparent transit depth.

We adopted a similar approach to Czesla et al. (2009), but worked on the deepest transit only in order not to average any activity feature. Figure 5 shows the series of transits with the best solution of fit overplotted. The sixth transit in the first segment, in yellow in Fig. 5, is the deepest.

This transit was isolated from data set . According to fit , this transit is affected by a bright spot, whose position during the transit is shown on the left side of Fig. 6. By fitting a single transit, KSint + PASTIS does not disentangle the bright spot from the transit profile. For this reason, we fixed the configuration of the active regions obtained with fit on the first segment of the light curve. Instead, the transit parameters kr, ip, ρ, and the jitter were set as free parameters. Because of the low number of points in a single transit, the limb-darkening coefficients, the flux offset, and the contamination value were fixed.

The MCMC analysis yielded , ip = 86.02 ± 0.27°, and ρ = 1.093 ± 0.038 ρ. The planet-to-star radius ratio is in 1σ agreement with the Czesla et al. (2009) result. The low values of ip and ρ, instead, have to be attributed to the distorted transit profile.

To check whether a dark spots-only solution can be found by our model, we fixed the planet-to-star radius ratio to the Czesla et al. value, and imposed three dark spots (i.e. with f< 1) for the fit of the deepest transit. A minimum of three spots was considered necessary. Indeed, two occulted dark spots at the borders of a transit can mimic a bright spot at the centre of the transit. A third non-occulted spot is needed to generate a possible out-of-transit flux variation if the other two are not sufficient. The contrast of the spots was fixed to the conservative solar value of 0.67 (Sofia et al. 1982). The latitude of the occulted spots was fixed close to −2.16° in order to lie on the transit chord. Their longitudes were forced to lie in the visible stellar disk to help the fit. The latitude of the non-occulted spot was set to 30°.

The best dark spots-only configuration is plotted in the centre of Fig. 6. This solution and the one with a bright spot are compared on the right side of the figure. Without a bright spot, the distortions of the transit profile are not recovered. The Bayesian Information Criterion (BIC; Schwarz 1978) favours the model with a bright spot over the one with dark spots only: BICfacula – BICdark spots ≃ −4. This corresponds to a Bayes factor ~e-4 ≃ 0.018 between the dark spots-only and the bright spot model.

This suggests that the unperturbed transit profile used by Czesla et al. (2009) may actually be a lower envelope of the transits affected by bright spots or faculae, and may so lead to an overestimate of kr. The transit profiles might be affected by bright spots, even if these spots have a negligible imprint on the out-of-transit part of the light curve. Including the transits in the modelling helps in recovering their presence.

Therefore, the unperturbed transit profile is more likely situated at a lower flux level than at the maximum level, which could indeed be produced by bright spots. Taking this into account, we inspected the modelled configuration of the system as in the previous cases, and we chose the twenty-sixth transit, in red in Fig. 5, as not affected by occulted activity features. We repeated the transit fit and obtained kr = 0.1689 ± 0.0008, i = 87.63 ± 0.53°, and ρ = 1.336 ± 0.062 ρ, in agreement with the results obtained by fit .

Figure 5 also shows that for many transits, the model is not able to correctly reproduce the complex structure of the transit profile. This does not invalidate our previous results, as we used transits which were correctly modelled. However, the plots indicate the presence of short-lived or irregularly shaped features inside transits which are not reproduced by the model. The low level of detail introduced by the re-sampling on data set and the consequent approximated modelling hid the need of additional features for transit modelling.

Silva-Valio & Lanza (2011) needed seven to nine features to model each transit. In our approach, instead, this number of features is used to globally model tens of days of data. In order to get closer to the correct number of features needed for both the in-transit and the out-of-transit data, one possibility is to perform fit on data set from the beginning, i.e. to use full resolution on the data. However, because of the large number of needed features and the large computational weight, we are not yet able to do this. We plan an optimization of the code in order to reduce its execution time and to be able to address the problem in a more efficient way.

5.2. Impact of stellar activity on the transit parameters

In Fig. 7, all the transit parameters found with fit are plotted as a function of the effective coverage factor of the stellar surface, introduced in Sect. 4.1. In each panel, the Spearman rank-order correlation coefficient, followed by its corresponding two-sided p-value in parentheses, is shown. We conservatively adopted a significance level for the p-value of 0.05. Our results therefore show that the hypothesis of non-correlation with is not rejected for all the transit parameters. We conclude that the observed scatter among the transit parameters is mainly due to statistical fluctuations caused by the division of the light curve in segments.

Several studies were presented where, without resorting to joint spot and transit modelling, the transit parameters are found to depend on the level of activity of the star. For example, a lower ρ attributed to starspots was observed for CoRoT-7 (Léger et al. 2009), compared to the derived spectroscopic value. Barros et al. (2014) later showed that this is mainly due to unresolved spot crossing events. Also, Csizmadia et al. (2013) showed that limb-darkening coefficients vary with the fraction of the stellar surface covered by activity features. Our results point to the importance of a joint modelling in order not to incur such effects.

thumbnail Fig. 7

Transit parameters as a function of the effective coverage factor, for all the segments. The values of the Spearman rank-order correlation coefficient are indicated. The corresponding p-values are in parentheses.

Open with DEXTER

5.3. Planet-to-star radius ratio and contamination

In our fit on locally normalized transits (fit ), we used a contamination value of about 8.81%. Alonso et al. (2008) used the same normalization, but used a fixed contamination value of 5.6%. They found kr ≃ 0.1667, while fit yields kr ≃ 0.1692. A difference of 3% in the contamination, therefore, produced an offset in kr of 1.5%. Instead, by keeping the same contamination value and a local (fit ) or model-based normalization (fit , kr ≃ 0.1682) we found a difference of only 0.6%. The local and the Czesla et al. normalization (fit ), finally, yield a 1.7% difference. The choice of the normalization, therefore, may introduce a systematic error in the transit depth just as the choice of the contamination value.

5.4. Planet radius

The combination of fit and allows, in principle, an unbiased measure of the transit parameters. Indeed, it is not affected by the transit normalization, and takes the occulted activity features into account. However, our modelling approach is limited by the need to cut the light curve into segments, which introduces a scatter in the results because of the lower number of points to be fitted.

Fit can be adopted as a good compromise. Although it neglects occulted spots, this fit is not affected by the normalization or by the chop of the light curve into segments. We therefore used kr and ρ obtained from fit to derive the radius of CoRoT-2 b. We used the Geneva stellar evolutionary tracks (Mowlavi et al. 2012) and the updated stellar atmospheric parameters of Torres et al. (2012; Teff=5575 ± 70 K, [Fe/H] = −0.04 ± 0.08).

This yielded RP = 1.475 ± 0.031RJ. If, instead, we adopt the results of a standard normalization (fit ), we obtain RP = 1.485 ± 0.031RJ, i.e. a less than 1% larger radius. These values are about ~1% different from the one found by Alonso et al. (2008; 1.465 ± 0.029RJ), but all of them are compatible. The inflated radius of the planet is confirmed.

5.5. Limitations of the current model

From the presentation of the results and their discussion, we identified some caveats which need to be addressed to improve the quality of our fits.

The first is the fixed number of activity features. An automatic incremental addition of the features could be implemented in order to minimize the residuals. A criterion should be chosen in order to stop the addition of features when residuals are smaller than a given threshold. This method would make the fit of short-lived features inside the transits easier. Second, the circular shape imposed on the features is another limitation of the model, and is linked to the previous one. The possibility of the features to overlap mitigates this problem, but many small features should overlap in order to model an active region with a complex shape.

Fitting λ and αmax might not be the best way to propagate the uncertainties from fit to fit . However, if more parameters of the activity features are left free during fit , the chains do not converge.

We note that all features in our model obey the same limb-darkening law. An additional level of detail could be gained by modelling a limb-angle dependent contrast for bright spots, that is, by including the modelling of real faculae.

The lack of constraints on the features evolution parameters is another point which needs refinements in order to improve the realism of the model. Moreover, the analytic model of KSint can be refined by the addition of differential rotation. Fröhlich et al. (2009) noticed that the information on differential rotation can be recovered only at the expense of some restrictions to the lifetimes of the features. However, if no differential rotation is modelled, two features at the same longitude but at different latitudes have a degenerate contribution to the total flux. Their contribution may therefore be equally modelled by a single, larger feature, which would simplify the convergence of the chains.

Most of these problems cannot be efficiently explored before the computation time required by PASTIS+ KSint is not importantly reduced. In fact, even if an analytic model was employed, some days of calculations were needed for the fit of each segment. This is due to the calculations required to model the overlap of active regions and transits, and to technical details of the code, which need to be optimized.

6. Summary and conclusions

We presented a method for the fit of transit photometry which takes the impact of activity features on the transit parameters into account. This approach is based on an improved version of the analytic code KSint, which models both non-occulted and occulted features and their evolution. The method was applied to the light curve of CoRoT-2 in two ways. The first is based on the simultaneous modelling of activity features and transits in a non-normalized light curve. The second consists in normalizing the transits by using the modelling of the out-of-transit light curve, and then in a standard fit of the so-normalized transits. In particular, we allowed for the presence of bright spots in our model, which was not usually done in previous studies.

The results of our method were compared with other approaches presented in literature. We recovered the total effective coverage factor of the stellar surface all along the light curve. We found different rotational periods of the active regions at different longitudes, as observed by Fröhlich et al. (2009).

We found that the choice of the normalization technique can introduce an offset in the transit depth of CoRoT-2 b as much as the choice of the contamination rate. Our results were then compared to those obtained with a standard local transit normalization, such as the technique used by Alonso et al. (2008). With a local normalization, and given the same contamination value, a 1.2σ increment of the measured planet-to-star radius ratio was found. Instead, the Czesla et al. (2009) normalization, where bright spots or faculae are neglected, yielded a 2.3σ smaller transit depth. The stellar densities, instead, were found to be in agreement whether a model-based or a local normalization was performed, and at less than 2σ agreement if the Czesla et al. normalization was performed. We used the transit parameters obtained by our modelling to calculate the radius of CoRoT-2 b, and confirmed its inflated nature (1.465 ± 0.029RJ).

Our analysis highlights the importance of stellar activity modelling during transit fit. Bright spots are particularly relevant in this respect. Including bright spots allowed us to correctly model the deepest transit profile without the need of increasing the transit depth, unlike previous studies (Czesla et al. 2009; Silva-Valio et al. 2010). Neglecting bright spots leads to the use of the lower envelope of the transits as an estimate of the unperturbed transit profile. We showed that the transit depth might be overestimated this way by almost 3% and worsen the inflated radius issue. Instead, our method shows that the correct transit profile might be closer to the value given by an average run over all the transits. Other cases of planets transiting active stars need to be analysed, however, in order to better constrain this result.

Thanks to the modelling of the time evolution of the features, it was possible to model longer parts of the light curve with respect to previous attempts presented in literature. Nevertheless, it was still necessary to cut the light curve in segments and to model each segment separately. A slight scatter was found among the transit parameters in the various segments. We showed that this has a likely statistical origin, given the small number of transits in each segment. Developments that would allow for a more complex evolution of the activity features and the modelling of longer segments of the light curve would therefore remove the scatter among the transit parameters.

Improvements in the fitting techniques are also needed in order to model light curves affected by stellar activity spanning two or three years of observations, such as those which will be provided by PLATO 2.0 (Rauer et al. 2014).

In order to develop our study in all these directions, we conclude by reaffirming the need of in-depth studies of other benchmark cases.


1

The transit shape is degenerate with respect to the stellar hemisphere that the planet covers. In the formalism used by KSint, inclinations <90° cover the southern hemisphere of the star, and the latitudes of the corresponding occulted activity features are negative.

2

The technical description of the pipeline is available at http://idoc-corotn2-public.ias.u-psud.fr/\jsp/CorotHelp.jsp

3

Hereafter, for simplicity, EBOP.

Acknowledgments

S.C.C.B. acknowledges support by grants 98761 by CNES and the Fundação para a Ciência e a Tecnologia (FCT) through the Investigador FCT Contract No. IF/01312/2014. A.S. is supported by the European Union under a Marie Curie Intra-European Fellowship for Career Development with reference FP7-PEOPLE-2013-IEF, number 627202. Part of this work was supported by Fundação para a Ciência e a Tecnologia, FCT, (ref. UID/FIS/04434/2013 and PTDC/FIS-AST/1526/2014) through national funds and by FEDER through COMPETE2020 (ref. POCI-01-0145-FEDER-007672 and POCI-01-0145-FEDER-016886). We thank S. Csizmadia and N. Espinoza for the fruitful discussions about limb darkening, D. Kipping for his help concerning the use of the macula code (Kipping 2012), R. Alonso for kindly providing us data to test our method, and R. F. Díaz for his explanations about Bayesian methods.

References

Appendix A: Additional tables

Table A.1

Main parameters of the CoRoT-2 system.

Table A.2

Transit parameters with their 68.3% credible intervals for each segment of fit .

Table A.3

Mean and standard deviations of the transit parameters on segments of the light curve with increasing length.

All Tables

Table 1

Data sets prepared for the modellings.

Table 2

Prior distributions used in the combined fit with KSint + PASTIS.

Table 3

Transit parameters with their 68.3% credible intervals.

Table A.1

Main parameters of the CoRoT-2 system.

Table A.2

Transit parameters with their 68.3% credible intervals for each segment of fit .

Table A.3

Mean and standard deviations of the transit parameters on segments of the light curve with increasing length.

All Figures

thumbnail Fig. 1

Model light curve from KSint plotted over the data. The eight segments of the fit are divided by colour. The residuals are shown in the lower panels; the error bars are not shown for clarity. The larger amplitude of the residuals in correspondence with the transits is due to the full resolution kept for the transits. The out-of-transit binning is of 2016 s, and inside the transits 32 s.

Open with DEXTER
In the text
thumbnail Fig. 2

Absolute value of the effective coverage factor as a function of time. The vertical dashed lines mark the separations between segments.

Open with DEXTER
In the text
thumbnail Fig. 3

Distribution of the activity features as a function of their longitude and time. No latitudinal information is reported. Filled circles represent dark spots, empty circles represent bright ones. Different segments of the light curve are represented by circles of different colour. The relative size of the features compared to the stellar surface was increased for the sake of illustration. Tentative reconstruction of longitudinal migration are represented by diagonal dashed lines.

Open with DEXTER
In the text
thumbnail Fig. 4

From top to bottom: values and uncertainties of kr, ip, ρ, Porb,ua, and ub for fit , indicated as a function of the segment of the light curve are in black. The results from fit with their 68.3% credible intervals are shaded in grey.

Open with DEXTER
In the text
thumbnail Fig. 5

Fit on the transits for every segment (from left to right, and top to bottom). The deepest transit, containing an occulted bright spot, is shown in yellow. A transit not affected by spots is in red.

Open with DEXTER
In the text
thumbnail Fig. 6

Left: planet (black) and activity features (blue) configuration during the deepest transit, with bright spots allowed. The bright spots are crossed by the planetary disc. Centre: solution with three dark spots (red). Right: deepest transit fitted with the dark-bright spot configuration (blue) and the three-dark spots model (red, shifted). The residuals are shifted for clarity and use the same colour code.

Open with DEXTER
In the text
thumbnail Fig. 7

Transit parameters as a function of the effective coverage factor, for all the segments. The values of the Spearman rank-order correlation coefficient are indicated. The corresponding p-values are in parentheses.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.