PyTranSpot: A tool for multiband light curve modeling of planetary transits and stellar spots

Ines G. Juvan; M. Lendl; P. E. Cubillos; L. Fossati; J. Tregloan-Reed; H. Lammer; E. W. Guenther; A. Hanslmeier

doi:10.1051/0004-6361/201731345

Home

All issues

Volume 610 (February 2018)

A&A, 610 (2018) A15

Full HTML

Free Access

Issue		A&A Volume 610, February 2018


Article Number		A15
Number of page(s)		17
Section		Planets and planetary systems
DOI		https://doi.org/10.1051/0004-6361/201731345
Published online		13 February 2018

A&A 610, A15 (2018)

`PyTranSpot`: A tool for multiband light curve modeling of planetary transits and stellar spots

Ines G. Juvan¹^,2^,3, M. Lendl¹^,4^,5, P. E. Cubillos¹, L. Fossati¹, J. Tregloan-Reed⁶, H. Lammer¹, E. W. Guenther⁷^,8^,⋆ and A. Hanslmeier²

¹ Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042 Graz, Austria
e-mail: Ines.Juvan@oeaw.ac.at
² Institut für Geophysik, Astrophysik und Meteorologie, Karl-Franzens-Universität, Universitätsplatz 5, 8010 Graz, Austria
³ Institut für Astro- und Teilchenphysik, Universität Innsbruck, Technikerstrasse 25, 6020 Innsbruck, Austria
⁴ Observatoire Astronomique de l’Université de Genève, Chemin des Maillettes 51, 1290 Sauverny, Switzerland
⁵ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
⁶ Carl Sagan Center, SETI Institute, Mountain View, CA 94043, USA
⁷ Thüringer Landessternwarte Tautenburg, Sternwarte 5, 07778 Tautenburg, Germany
⁸ Instituto de Astrofísica de Canarias (IAC), vía Láctea s/n, 38205 La Laguna, Tenerife, Spain

Received: 9 June 2017
Accepted: 24 October 2017

Abstract

Several studies have shown that stellar activity features, such as occulted and non-occulted starspots, can affect the measurement of transit parameters biasing studies of transit timing variations and transmission spectra. We present PyTranSpot, which we designed to model multiband transit light curves showing starspot anomalies, inferring both transit and spot parameters. The code follows a pixellation approach to model the star with its corresponding limb darkening, spots, and transiting planet on a two dimensional Cartesian coordinate grid. We combine PyTranSpot with a Markov chain Monte Carlo framework to study and derive exoplanet transmission spectra, which provides statistically robust values for the physical properties and uncertainties of a transiting star-planet system. We validate PyTranSpot’s performance by analyzing eleven synthetic light curves of four different star-planet systems and 20 transit light curves of the well-studied WASP-41b system. We also investigate the impact of starspots on transit parameters and derive wavelength dependent transit depth values for WASP-41b covering a range of 6200−9200 Å, indicating a flat transmission spectrum.

Key words: planetary systems / planets and satellites: individual: WASP-41b / planets and satellites: fundamental parameters / starspots / techniques: photometric

^⋆

Visiting scientist at the IAC.

© ESO, 2018

1. Introduction

To date, over 2600¹ exoplanets have been confirmed and detected by means of the transit method². This method is based on the measurement of a periodic dimming in a stellar light curve, caused by a transiting exoplanet (TEP) passing in front of its host star. By fitting a model to a transit light curve (TLC), one obtains orbital and photometric parameters of the star-planet system, such as the planetary period P_orb, the orbital incination i, and the planet-to-star radius ratio (e.g., Charbonneau et al. 2000; Winn 2010). In combination with radial velocity or transit timing variation (TTV) measurements, it is possible to determine the mass of the planet, hence, the average density. Multi-wavelength transmission spectroscopy measurements can constrain the planet’s atmospheric properties by comparing the wavelength-dependent variations in the planetary radii with theoretical model atmospheres (e.g., Charbonneau et al. 2002; Agol et al. 2005; Fortney et al. 2008; Madhusudhan & Seager 2009; Miller-Ricci et al. 2009).

As already thoroughly discussed (e.g., Sanchis-Ojeda et al. 2011; Ballerini et al. 2012; Oshagh et al. 2013b; Barros et al. 2013; Tregloan-Reed et al. 2013, 2015), the determination of planetary parameters can often be challenging due to the presence of stellar activity features in photometric data sets. Starspots, features which are cooler and thus darker than the surrounding stellar photosphere, can introduce anomalies (“bumps”) in a TLC when they are occulted by a transiting planet (e.g., Silva 2003; Pont et al. 2007; Rabus et al. 2009; Winn et al. 2010b). The improper treatment of starspots in a TLC fitting process can thus lead to an incorrect determination of the depth, duration, and timing of the transit. Additional effects can be introduced by unocculted spots or bright features of stellar activity such as faculae or plages (e.g., Czesla et al. 2009; Désert et al. 2011b; Kipping 2012; Zellem et al. 2017; Rackham et al. 2017). The impact of these features on the light curve depends on their location relative to the planetary path across the stellar disk. Unlike occulted starspots, activity features, located in the non-eclipsed area of the stellar surface, do not cause distinct anomalies in a TLC, but have an impact on the overall level of the light curve (Czesla et al. 2009).

Although spots represent sources of noise in a TLC, they can also be seen as useful features to obtain additional information on the observed star-planet system (Kipping 2012). By modeling spots in TLCs, one can constrain the following properties: the latitudinal stellar rotation period P_⋆, which yields a value for the stellar latitudinal rotational velocity, from which it is then possible to infer the star’s age and activity level (e.g., the gyrochronology relationship Barnes 2007). However, star-planet interactions may influence gyrochronology age estimations (e.g., Ferraz-Mello et al. 2015). Furthermore, one can obtain parameters such as the sky-projected spin-orbit alignment λ, and, together with an estimate of the inclination angle of the stellar rotation axis and the stellar rotational velocity, the true obliquity of the system can be derived. The determination of the true obliquity then also provides clues to the dynamical evolution of the system (e.g., Fabrycky & Winn 2009; Winn et al. 2010a; Nutzman et al. 2011; Sanchis-Ojeda & Winn 2011; Désert et al. 2011a).

So far, several teams have developed routines to model TLCs in the presence of starspots: the spotrod routine (Béky et al. 2014) uses a semi-analytical approach to model the TLC of a spotted star, whereas other authors mainly use numerical methods for their astrophysical models: SOAP-T (Oshagh et al. 2013a), PRISM (Tregloan-Reed et al. 2013, 2015), KSint (Montalto et al. 2014), ellc (Maxted 2016), and StarSim (Herrero et al. 2016). Relatable programs have also been developed within the binary star community, such as the Wilson-Devinney (WD) code (Wilson & Devinney 1971; Wilson 1979, 1990, 2008, 2012, and references therein) and PHEOBE (Prša & Zwitter 2005; Prša et al. 2016).

Motivated by the large number of available multicolor photometric transit observations, we developed PyTranSpot³, a tool which allows for the simultaneous analysis of TLCs in the presence of stellar activity and correlated noise. To perform simultaneous analyses of TLCs, derived from various instruments and in different wavelength bands, we combined PyTranSpot with a Markov chain Monte Carlo (MCMC) framework designed for the determination of exoplanet transmission spectra (Lendl et al. 2017).

This paper is structured as follows: Sect. 2 introduces PyTranSpot. Section 3 presents the validation of the code using four synthetic star-planet systems. Section 4 describes the analysis and results of 20 WASP-41b TLCs. Section 5 presents the summary and conclusions of our study.

2. Modeling transits and starspots with `PyTranSpot`

2.1. Astrophysical model and geometry

PyTranSpot is based on a pixellation approach, similar to the one used in the PRISM code (Tregloan-Reed et al. 2013, 2015). Within this approach, the pixels are defined as squares on a two-dimensional grid, on which the stellar sphere (and its respective limb darkening), the transit cord, and spots are projected (see Fig. 1). The pixel size of the star is determined through dividing the a-priori defined planetary pixel radius by the planet-to-star radius ratio. We assume dark and bright features of stellar activity to be homogeneous and circular over the stellar surface, and they deform to ellipses as they approach the stellar limb. From solar observations we know that sunspots can also appear as a complex group of multiple spots with differing sizes. However, the quality of currently obtained TLCs is in general not high enough to detect such fine structures. Furthermore, we assume that the stellar rotation period is much longer than the transit.

Fig. 1

Projection of the stellar sphere with its respective limb darkening, the location of the spot, and the transit cord on the two-dimensional grid in Cartesian coordinates. The colorbar on the right hand side indicates the intensity of the stellar flux. Stellar activity features on the stellar sphere are assumed to be homogeneous and circular.

We calculate the transit light cuve model using: $\begin{matrix} \begin{matrix} \begin{matrix} \end{matrix} Δ F = \frac{F_{out} - F_{transit}}{F_{out}}, \end{matrix} \end{matrix}$ $\begin{eqnarray} \begin{gathered} \Delta F = \frac{F_{\rm{out}}-F_{\rm{transit}}}{F_{\rm{out}}}, \label{eq:flux} \end{gathered} \end{eqnarray}$ (1)where ΔF represents the flux measurement for every timestamp (in and out of transit). The total out-of-transit flux F_out also takes into account stellar limb darkening and apparent starspots. From this, we subtract F_transit, which describes the fraction of flux on the stellar sphere occulted by the transiting planet. To derive a normalized light curve, we divide F_out−F_transit through F_out. It is also possible to multiply ΔF with a selected photometric baseline model.

Figure 2 shows the geometry used within our model. The center of the stellar sphere is located at the origin of the three-dimensional spherical coordinate system. PyTranSpot does not take into account a fractional area correction, as used within the WD code (Wilson & Devinney 1971; Wilson 1979, 1990, 2008, 2012). We argue that this effect is negligible, as the resulting loss of accuracy is much smaller than the noise currently present in observations. However, to obtain precise photometric transit and spot parameters, we recommend the use of a planetary pixel radius between 15 and 50 pixels (Tregloan-Reed et al. 2015). On the stellar sphere, every point is described by the two angles (longitude θ, co-latitude φ) and the distance to the stellar center (stellar radius r_s in pixels). The longitude θ varies between − 90° and 90°, with the center of the stellar disk corresponding to a value of 0°. The co-latitude φ ranges from 0° to 180°, with the stellar equator set at 90°. Figure 3 illustrates the projection of a spot on to the stellar sphere, as seen from a two-dimensional perspective. The observer is assumed to lie far along the z-axis. To determine the pixels on the stellar sphere, which correspond to the starspot, we implement the following boundary condition: If the angle Δσ between the pixel on the sphere and the spotcenter is greater than the angular radius of the spot α, then this pixel does not belong to the spot. The values for Δσ are derived by using the spherical law of cosines: $\cos (Δ σ) = \cos (φ_{spot}) \cdot \cos (φ) + \sin (φ_{spot}) \cdot \sin (φ) \cdot \cos (Δ θ),$ $\begin{equation} \cos(\Delta\sigma) = \cos(\phi_{\rm{spot}}) \cdot \cos(\phi) + \sin(\phi_{\rm{spot}}) \cdot \sin(\phi) \cdot \cos(\Delta\theta), \label{eq:slaw} \end{equation}$ (2)where φ_spot and φ are the co-latitudes of the spotcenter and the surrounding pixels, respectively. The value Δθ represents the absolute difference in longitude between the spotcenter and the pixel center.

Fig. 2

Geometry of the stellar sphere and spot feature. The origin of the coordinate system is the center of the star and spot parameters are defined using spherical coordinates. The gray circle on the stellar surface represents a spot at the distance r_s (here: stellar radius in pixels) from the origin. The longitude θ (blue) can have values between − 90° and 90°, whereas 0° represents the center of the stellar disk. The co-latitude φ (purple) is defined between 0° to 180°, with the equator set at 90°, as seen from an observer lying far along the z-axis.

Fig. 3

Two-dimensional cut through the stellar sphere, at the position of the spotcenter in the direction perpendicular to the line of sight. The green dashed line represents the spot feature as it would be seen in three dimensions on the stellar sphere. The size of the spot is determined by the angle α (blue). The minor arc d (red) between the spotcenter and any pixel location on the spot can be found through the angle Δσ (red) times the stellar radius in pixels r_s. Values of Δσ are derived from Eq. (2). Every point on the stellar sphere, corresponding to the spot feature, needs to fulfil the condition: Δσ ≤ α.

2.2. Model parameters

PyTranSpot calculates the transit and spot model by using the following photometric and orbital parameters:

phase offset from the orbital phase = 0.0, which indicates thetransit midpoint;
planet-to-star radius ratio r_P/r_S;
orbital inclination i in degrees;
semi-major axis a in units of the stellar radius a/R_S;
planetary orbital period P_orb in days;
orbital eccentricity e;
argument of periastron ω in degrees;
linear u₁ and quadratic u₂ limb darkening coefficients;
coefficients of the baseline model functions.

To model one or more spots on the stellar sphere, each starspot is characterized by the following set of parameters:

spot longitude θ in degrees ( $\hbox{$-90 \degr \leqslant \theta \leqslant 90 \degr $}$ );
spot co-latitude φ in degrees (0° to 180°);
spot size α in degrees;
spot contrast ρ_spot (ρ_spot ∈ [0,1], where 1 equals the surrounding stellar photosphere. This limitation corresponds to the modeling of (un)occulted dark spots. When modeling bright features (e.g., faculae or plages), the contrast can take values of ρ_spot> 1.0. However, defining the contrast of, for example, a plage region is more complex as it depends on its location on the stellar sphere and the stellar type (e.g., Beeck et al. 2013a,b; Thaler & Spruit 2014)).

2.3. Treatment of stellar limb darkening

PyTranSpot employs the quadratic limb darkening law (LDL; Kopal 1950), which is the most commonly used LDL in TLC analyses. The quadratic LDL describes the specific intensity of a star I(μ) on the surface as $I (μ) / I (0) = 1 - u_{1} (1 - μ) - u_{2} (1 - μ)^{2},$ $\begin{equation} I(\mu)/I(0) = 1 - u_{1}(1-\mu) - u_{2}(1-\mu)^2, \label{eq:LD} \end{equation}$ (3)where I(0) defines the intensity at the center of the stellar disk, μ is the cosine of the angle between the line of sight of the observer and the unit normal to the stellar surface, and u₁ and u₂ are the quadratic limb darkening coefficients (LDCs).

We make this choice because, unlike three- or four-parameter LDLs, the quadratic two-parameter law preserves the curvature in the intensity of the star, without trying to model bumps due to starspots in the light curve (Csizmadia et al. 2013). Also, using a two-parameter law reduces the number of free parameters in the model (Kipping 2012), which is especially important when analyzing large data sets that require a great number of MCMC jump parameters. To guarantee physical values of the quadratic LDCs u₁ and u₂, we implement the following conditions proposed by Kipping (2013) within our astrophysical model: $\begin{matrix} \begin{matrix} \begin{matrix} \end{matrix} u_{1} + u_{2} < 1, \begin{matrix} \end{matrix} \\ \begin{matrix} \end{matrix} u_{1} > 0, \\ \begin{matrix} \end{matrix} u_{1} + 2 u_{2} > 0. \end{matrix} \end{matrix}$ $\begin{eqnarray} \begin{aligned} u_{1} + u_{2} & < 1, \\ u_{1} & > 0, \\ u_{1} + 2u_{2} & > 0 . \label{eq:kipp1} \end{aligned} \end{eqnarray}$ (4)These conditions make sure that the intensity profile remains positive everywhere, and guarantees a monotonically decreasing intensity profile from the center of the stellar sphere to the limb.

2.4. Determination of transmission spectra

We use PyTranSpot in combination with the MCMC framework developed by Lendl et al. (2017). This MCMC framework employs the statistical package MCcubed (see Sect. 3.2.1) and allows for the simultaneous analysis of multiband transit observations yielding filter dependent planet-to-star radius ratios r_P/r_S. When calculating a transmission spectrum using multiband TLCs, the light curves share the same model parameters (see Sect. 2.2). Exceptions are the filter dependent LDCs, the baseline model coefficients, and the spot parameters. These parameters are fitted individually. To derive the multiband planet-to-star radius ratios, we fit a wavelength dependent offset to a reference r_P/r_S value (Lendl et al. 2017).

3. Code validation using synthetic light curves

3.1. Synthetic data sets

As a first approach to test the performance of PyTranSpot, we create four synthetic star-planet systems yielding a total of eleven light curves. The first system SYNTH-1 (Table 1) consists of a hot Jupiter planet with a solar-like host star. We derive one TLC without stellar activity features (SYNTH-1a). The star-planet system SYNTH-2 (Table 2) hosts a Saturn-size planet orbiting an active, solar-like star. The generated light curve (SYNTH-2a) shows a starspot crossing at the limb of the star. The third system SYNTH-3 is similar to SYNTH-2, but with a shorter orbital period, and we create five synthetic observations of SYNTH-3. The five light curves of this system are referred to as SYNTH-3a, -3b, -3c, -3d, and -3e. We further create one of the light curves (SYNTH-3e) with an anomaly due to an occulted starspot. Table 3 shows the system properties and the spot parameters that refer to the SYNTH-3e light curve. System SYNTH-4 (Table 4) also describes a hot Jupiter planet orbiting a solar-like star. For this system, we simulate a simultaneous observation of the same transit event, measured in the Johnson B, V, and Cousins R, I filters. We further assume that the corresponding r_P/r_S value has no wavelength dependence. The resulting light curves SYNTH-4a to SYNTH-4d show the same occulted starspot. We calculate the wavelength-dependent spot contrasts using Eq. (1) of Silva (2003), using a blackbody approximation, a solar-like effective temperature for the host star of T_eff = 5772 K, and a spot temperature of T_spot = 4772 K. We generate all simulated TLCs using PyTranSpot and add Gaussian noise to the calculated flux.

Table 1

Photometric properties of the synthetic light curve SYNTH-1a.

Table 2

Photometric properties of the synthetic light curve SYNTH-2a.

Table 3

Photometric properties of the synthetic system SYNTH-3.

Table 4

Photometric properties of the synthetic system SYNTH-4.

3.2. Light curve analysis

We analyze the synthetic light curves of the systems SYNTH-1 and SYNTH-2 individually, whereas the TLCs SYNTH-3a to SYNTH-3e, and SYNTH-4a to SYNTH-4d are analyzed simultaneously. For each MCMC analysis, we run ten chains with a total of 600 000 samples. For the simultaneous analysis of SYNTH-4a to SYNTH-4d, we use 23 chains with a total of 2 000 000 samples due to the larger number of free parameters. The MCMC jump parameters are those listed in Sect. 2.2, except for the planetary orbital period, the eccentricity, and the argument of periastron, which are fixed. We also do not consider photometric baseline models during these analyses. Furthermore, we analyze the spotted TLCs (SYNTH-2a and SYNTH-3e) assuming a spot-free transit model to investigate the impact of starspots on transit parameters. The simultaneous multiband observations of SYNTH-4 are used to study the ability of PyTranSpot to reproduce the transit and spot parameters while fitting for the filter-dependent LDCs and spot contrasts.

We recalculate errorbars of each data set using uncorrelated (white) and correlated (red) noise factors (Winn et al. 2008; Gillon et al. 2010). This guarantees that uncertainties are not being underestimated in the course of the analysis. We perform the analysis for each system multiple times (at least ten repetitions for SYNTH-1 and SYNTH-2, and three for SYNTH-3 and SYNTH-4) to make sure that the obtained results are consistent and thus robust.

3.2.1. Statistical package

To carry out the statistical analysis, we use the open-source package Multi-Core MCMC (MCcubed; Cubillos et al. 2017). MCcubed⁴ is a Python/C code that provides statistically-robust model optimization via Levenberg-Marquardt minimization and credible-region estimation via MCMC sampling. MCcubed assesses the goodness-of-fit between the model and data through χ² statistics, considering user-defined flat or Gaussian priors. To sample the parameter space, we choose the Snooker Differential-Evolution algorithm (ter Braak & Vrugt 2008), which automatically adjusts the scale and orientation of the proposal distribution. The code checks for MCMC convergence through the Gelman & Rubin (1992) statistics.

Fig. 4

Transit light curve (blue dots) of the synthetic star-planet system SYNTH-1, with the derived best-fit model (red line) and the original light curve model (green line). The obtained photometric parameters can be found in Table B.1. The (O−C) residuals and the difference between the original and best-fit light curve models are presented in Fig. A.1.

3.3. Results

3.3.1. Transit and spot parameters

We present in Figs. 4–7 the synthetic light curves with the resulting best-fit models. Figures A.1−A.6 show the (O−C) residuals and a visualization of the differences between the input- and best-fit light curve models. Tables B.1−B.4 give the derived transit and spot properties. A comparison of the obtained parameters of systems SYNTH-1, SYNTH-2, and SYNTH-3 with their original system parameters (Tables 1, 2, and 3) shows that we can recover the majority of the input values consistently within one sigma. Only the LDCs of SYNTH-1a slightly differ from the original input values, but the differences are smaller than 1.3σ. We also find that the derived results of the repeated MCMC analysis for each system yield consistent values. The transit parameters, spot locations, filter-dependent LDCs, and spot contrasts of our multiband SYNTH-4 TLCs could also be reproduced within 1.2σ. An exception is the derived spot contrast for light curve SYNTH-4a (B filter), which differs by 2σ from the original value. However, using the derived spot contrasts, a stellar effective temperature of T_eff = 5772 K, and Eq. (1) of Silva (2003), we derive an average spot temperature of T_spot = $495 6_{-175}^{+ 245}$ $\hbox{$4956^{+245}_{-175}$}$ K. This value agrees within 1.1σ with the original spot temperature of T_spot = 4772 K. The ability of our code to reproduce also the wavelength-dependent spot contrasts is an important result as simultaneous multiband observations of starspots can help to constrain starspot temperatures. Without such simultaneous measurements, the spot temperature remains strongly correlated with its radius (e.g., Tregloan-Reed et al. 2013; Mancini et al. 2014).

Fig. 5

Synthetic light curve SYNTH-2a (blue dots), with the derived best-fit model (red line) and the original light curve model (green line). This synthetic light curve shows a starspot crossing at the limb of the star around phase −0.01. The obtained photometric properties and spot parameters can be found in Tables B.1 and B.3. The (O−C) residuals and the difference between the original and best-fit TLC models are shown in Fig. A.2.

Fig. 6

Synthetic light curves SYNTH-3a to SYNTH-3e, with the derived best-fit models (red line) and the original light curve model (green line). Light curve SYNTH-3e is showing a starspot anomaly located close to the limb of the star around phase −0.005. The obtained photometric properties and spot parameters can be found in Tables B.1 and B.3. The (O−C) residuals and the difference between the original and best-fit TLC models are shown in Figs. A.3 and A.4, respectively.

Fig. 7

Synthetic light curves SYNTH-4a to SYNTH-4d, with the derived best-fit models (red line) and the original light curve model (green line). The TLCs all show the same transit event and occulted starspot around phase −0.002. The obtained photometric and spot parameters can be found in Tables B.2 and B.4. The (O−C) residuals and the difference between the original and best-fit TLC models are shown in Figs. A.5 and A.6, respectively.

We also perform a test to calculate the theoretically expected errorbar on the planet-to-star radius ratio to compare it with the one derived from our analysis. For this, we calculate the expected uncertainty on r_P/r_S of SYNTH-1. We assume that our data is only affected by white noise and obtain a theoretical uncertainty σ_{r_P/r_S} of 0.00052. We find that our derived uncertainty is more than approximately three times larger than the theoretical one. This is a reasonable result as we assume the transit shape to be a simple trapezoid when calculating the expected uncertainty, but in practice the TLC model is more complex (i.e., it has more free parameters). In addition, to quantify the effect of a lower pixel resolution on the derived parameters and uncertainties, we re-analyze the TLCs SYNTH-1 and SYNTH-2. We perform identical analyzes, but vary the planetary pixel radius r_P = 10,15,20 from the original one r_P = 50. We find that for all cases, the derived values using a lower pixel resolution differ by ≤ 0.1σ from the results using r_P = 50. We also calculate the differences between the low resolution and the r_P = 50 light curve models and derive rms values between 8−20 ppm, which are much lower than the noise in SYNTH-1 (700 ppm) and SYNTH-2 (400 ppm). This agrees with Tregloan-Reed et al. (2015), who also find that a higher pixel resolution does not significantly increase the numerical resolution of the resulting parameters and uncertainties. However, the use of r_P = 50 is recommended, as a smaller planet pixel radius does affect the smoothness of the resulting best-fit model.

3.3.2. Impact of starspots on transit parameters

We compare the transit parameters, which we derive from the analyses of the spotted TLCs SYNTH-2a and SYNTH-3e (assuming both a spot and spot-free model), and find the following: The majority of the transit parameters agree well within 1.4σ. We also identify that the phase offset value, which we derive from SYNTH-2a’s spot-free model, differs by 2.5σ. Not taking into account the spot feature, which is located at the limb of SYNTH-2, seems to affect the determination of the transit duration and hence, the transit midtime. The effect of starspots on the measured transit duration and timing confirms findings by various authors (e.g., Silva-Válio 2010; Sanchis-Ojeda et al. 2011; Oshagh et al. 2013b).

4. WASP-41b: a broadband 6200–9200 Å transmission spectrum in the presence of starspots

To further test the performance of PyTranSpot, we use data of the well-studied WASP-41 system. WASP-41 is one of the targets of the Wide Angle Search for Planets project (WASP⁵; Pollacco et al. 2006), and is a V = 11.6 G8V star, which is known to show magnetic activity and rotational modulation on a period of 18.41 ± 0.05 days (Maxted et al. 2011). The system’s transiting hot Jupiter, WASP-41b, has a measured planetary mass and radius of 0.94M_Jup and 1.06R_Jup, respectively. The study of Southworth et al. (2016) on WASP-41b states that some of the TLCs show anomalies in brightness due to occulted spots. Table 5 presents the derived spot parameters of Southworth et al. (2016). From modeling two spot features, and assuming that these are caused by the same spot, the authors determined the rotation period of the host star to be P_⋆ = 18.6 ± 1.5 days, and a sky-projected orbital obliquity of λ = 6 ± 11°. Since the host star is magnetically active showing TLCs with and without starspots, WASP-41b represents an ideal object to further test our routine. Our aim is to reproduce the results of Southworth et al. (2016), to investigate the effect of starspots on the transit parameters, and to derive a broadband transmission spectrum in the range of 6200−9200 Å for WASP-41b. To accomplish this, we use archival data in different wavelength bands together with yet unpublished TLCs of WASP-41b.

Table 5

Photometric and starspot properties of the WASP-41 system, as taken from Southworth et al. (2016).

Table 6

Observing log of EulerCam observations of WASP-41.

Table 7

WASP-41b observations analyzed in this work.

4.1. Data

4.1.1. EulerCam observations

We observed a total of nine transits of WASP-41 between January 2011 and April 2012 with EulerCam, the CCD imager installed at the 1.2 m Euler telescope at ESO La Silla, Chile. From these nine unpublished TLCs, three show evidence of occulted stellar spots. The observations were carried out using a r′-Gunn filter and the telescope was slightly defocused to improve efficiency and point spread function (PSF) sampling. Table 6 gives a summary of the individual observations. We reduced the data using aperture photometry and tested a range of apertures and reference stars, selecting those which produce the minimal residual scatter of the fitted TLC. Lendl et al. (2012) give details on the instrument and the data reduction procedures.

4.1.2. DFOSC I observations

For our study, we also consider a set of four TLCs from Southworth et al. (2016). The authors observed four transits of WASP-41b with the Danish Faint Object Spectrograph and Camera (DFOSC) instrument, which is installed on the 1.54 m Danish Telescope at ESO La Silla, Chile. The object was observed using a Bessell I filter. Southworth et al. (2016) describe the observations and data reduction of these data. The DFOSC data are of specific interest for this paper as two of the four light curves show occulted starspots. The authors also observed two additional transits of WASP-41b using the 84 cm telescope at Observatorio Cerro Amazones in Antofagasta. Due to their lower quality, these two TLCs are not included in our study.

4.1.3. TRAPPIST and DFOSC R observations

Neveu-VanMalle et al. (2016) present eight transits of the WASP-41 system from which we adopt the five data sets obtained with the TRAnsiting Planets and PlanetesImals Telescope-South (TRAPPIST; Gillon et al. 2011; Jehin et al. 2011) in the I + z filter, and two TLCs observed with the DFOSC instrument in the Bessell R filter (see Table 7). We note that only one of the DFOSC R light curves covers the full transit. In addition, we decide not to include the TLC observed with the Faulkes Telescope South (FTS) telescope located at Siding Spring Observatory, as the observation was affected by poor weather conditions.

4.2. Light curve analysis

We model the WASP-41b TLCs individually, as well as simultaneously, using PyTranSpot within the MCMC transmission spectroscopy framework (Lendl et al. 2017). For all individual and simultaneous fitting processes, we additionally analyze the spotted TLCs assuming a spot-free model. From the simultaneous analyses, we further derive a transmission spectrum in the range of 6200−9200 Å for WASP-41b. We remark that we infer only one R-band planet-to-star radius ratio for the transmission spectrum, combining the EulerCam r′ and DFOSC R TLCs.

Within the simultaneous MCMC analysis, the TLCs share the same transit parameters (see Sect. 2.2), except for the bandpass-dependent LDCs. In addition, the spot parameters and baseline coefficients are analyzed separately for each light curve. The wavelength-dependent r_P/r_S values are derived through fitting an offset to a reference planet-to-star radius ratio. Following Gillon et al. (2010), the coefficients describing the baseline models are calculated for every MCMC step by applying least-square minimization (Lendl et al. 2017). As discussed in Southworth et al. (2016), we also fix the orbital inclination to the value i = 88.7°, which restricts the strong correlation between the planet’s orbital inclination and the spot latitude. Furthermore, we fix the planetary orbital period to P_orb = 3.05 days, as well as the eccentricity e and the argument of periastron ω, which are both set to zero. We use the quadratic LDCs, which we inferred from JKTLD⁶ (Southworth 2008), as starting values for our MCMC analysis. Errorbars for each light curve are rescaled as discussed in Sect. 3.2. Each run consists of 10−20 parallel MCMC chains with a total of up to 1 600 000 samples. The final sample size depends on the number of light curves to be analyzed.

We visually inspect all TLCs and select a general baseline model of a quadratic polynomial in time to correct for time-dependent modulations. We also assume that some of the EulerCam observations must have suffered from coordinate drifts of the telescope, hence, we test the application of an additional quadratic polynomial in the telescope drift. The TRAPPIST light curves seemed to have experienced difficulties with the autofocus, as discussed in Neveu-VanMalle et al. (2016), resulting in significant variations of the full-width-half-maximum values (FWHM) for each image. Therefore, we also consider multiplying our light curve models with second-order polynomials with respect to the FWHM values. However, we only choose more complex models over our general (minimal) baseline model, if the derived Bayes factor (e.g., Schwarz 1978) implies a higher probability. The final baseline model for each light curve is presented in Table 7.

We find anomalies due to occulted starspots in the TLCs obtained with the DFOSC instrument (2015/05/13 and 2015/05/17), in the EulerCam (2011/04/02, 2011/05/15, and 2011/01/31), and TRAPPIST (2011/04/02) observations. The EulerCam (2011/04/02) and the TRAPPIST (2011/04/02) measurements observed the same transit event, hence, they show the same spot. To obtain accurate values for the spot location, size, and contrast, we also fit these data sets separately, and compare the results to the values that we derive from the simultaneous analysis.

Fig. 8

EulerCam light curves (blue dots) and best-fit models (red line) with the corresponding residuals shown in Fig. 9. The TLCs are presented with their respective date of observation. We note that the light curves with observing dates 2011/01/31, 2011/04/02, and 2011/05/15 were modeled with an occulted starspot feature. Results obtained from the simultaneous analysis are given in Tables B.5 and B.6, and the derived spot parameters from the simultaneous and individual analysis are presented in Tables B.7 and B.8.

Fig. 9

Residuals of the EulerCam light curves. The light curves with observing dates 2011/01/31, 2011/04/02, and 2011/05/15 were modeled with an occulted starspot feature around phase −0.002.

Fig. 10

DFOSC I light curves (blue dots) and best-fit models (red line). The corresponding residuals are shown in Fig. 11. The TLCs are listed with their respective date of observation, and the data with observing date 2015/05/13 and 2015/05/17 show one and two occulted spots, respectively. Results obtained from the simultaneous analysis are given in Tables B.5 and B.6 and the derived spot parameters from the simultaneous and individual analysis are presented in Tables B.7 and B.8.

Fig. 11

Residuals of the DFOSC I light curves, listed with their respective date of observation. The transit light curves with observing date 2015/05/13 and 2015/05/17 show one (phase −0.01) and two occulted spots (phases −0.005 and 0.005), respectively.

4.3. Results

4.3.1. WASP-41b system parameters

Figures 8−15 present the WASP-41b TLCs with their best-fit models and residuals from the individual analyses. Tables B.5 and B.6 then give the results from the simultaneous fitting processes for the two discussed cases (assuming a spot and a spot-free model). We decide to thoroughly present the simultaneous fitting results only, as the results of the individual TLCs agree within one sigma with the results from the combined fit. The only exceptions are the EulerCam TLC (2011/05/12) and the DFOSC R TLC (2013/04/23), which show discrepancies in the transit midtime, planet-to-star radius ratio, relative semi-major axis, and limb darkening coefficients to an extent of 4σ. This is likely the result of the lack of data points in the first (or second) half of the transit, which directly affects the accurate determination of these transit parameters.

Fig. 12

TRAPPIST light curves (blue dots) with the derived best-fit models (red line). The corresponding residuals are shown in Fig. 13. The TLCs are presented with their corresponding observing date. The TLC with the observing date 2011/04/02 shows a starspot at phase −0.005. Results, which we obtained from the simultaneous analysis, are given in Tables B.5 and B.6 and the derived spot parameters from the simultaneous and individual analysis are presented in Tables B.7 and B.8.

Fig. 13

Residuals of the TRAPPIST light curves (blue dots), given with their corresponding observing date. The TLC with the observing date 2011/04/02 shows a starspot at phase −0.005.

Fig. 14

DFOSC R light curves (blue dots) and obtained best-fit models (red line) with the corresponding residuals shown in Fig. 15. On the right hand side, we indicated the corresponding observing dates. We note that only one of the light curves shows a full transit event. Results obtained from the simultaneous analysis are given in Tables B.5 and B.6.

Fig. 15

Residuals of the DFOSC R light curves (blue dots) with the corresponding observing dates indicated on the right hand side.

We compare all our results, from the individual as well as the simultaneous analyses, with those presented in Southworth et al. (2016) and find that they agree within one sigma (with the exception of the EulerCam TLC (2011/05/12) and DFOSC R (2013/04/23)). We remark also that in the course of the simultaneous analysis, we identify that one of the TLCs observed with TRAPPIST (2011/03/21) leads to an excess value of the transit depth of r_P/r_S ~ 0.144, which differs by almost 2σ from the r_P/r_S results of the remaining TRAPPIST TLCs. We believe that this must be a result of the TRAPPIST autofocus issues (Neveu-VanMalle et al. 2016), as this TLC is affected the most by variations of the FWHM. Therefore, we exclude this light curve from the analysis.

4.3.2. Spot modeling results

At first, we verify if the results derived from the individual and simultaneous analyses are consistent. We find that all obtained spot locations, sizes, and contrasts agree within one sigma, which shows that our simultaneous fitting process can also produce reliable spot parameters. We also investigate the EulerCam light curves for possible reappearing spot features. We do not detect reappearing starspots as the observations were taken too far apart when compared to the estimated stellar rotation period of 18.6 ± 1.5 days (Southworth et al. 2016) and from rotational modulation, 18.41 ± 0.05 days (Maxted et al. 2011).

When we compare our individually and simultaneously derived spot parameters of the DFOSC I light curves (Tables B.7 and B.8) with the results published by Southworth et al. (2016), we find that they agree well within the uncertainties. We also conclude that the spot longitudes are in general well determined, while we find large uncertainties for the spot latitudes. Southworth et al. (2016) mention that they encountered the same issue and further argue that this is to be expected for systems in which the transit cord is close to the center of the stellar sphere (as is the case for WASP-41b).

One of the TRAPPIST TLCs (2011/04/02) shows the same spot feature as the EulerCam observation obtained during the same night. We remark that these observations were performed in different filters and, furthermore, the TRAPPIST measurements suffered from problems affecting the telescope autofocus, which resulted in large systematics. Nevertheless, the obtained spot parameters (spot longitude, co-latitude, contrast, and size) agree within 1.5σ. We also use the derived spot contrasts from the simultaneous multiband measurements of the same starspot (Table B.7) together with T_eff from Table 5 to calculate the spot temperatures using Eq. (1) of Silva (2003), assuming a black body approximation. We find that the obtained spot temperatures $T_{spot, TRAPPIST} = 5296_{-245}^{+ 119}$ $\hbox{$T_{\rm{spot, TRAPPIST}} = {5296}^{+119}_{-245}$}$ K and $T_{spot, Euler Cam} = 5220_{-126}^{+ 58}$ $\hbox{$T_{\rm{spot, {\it Euler}Cam}} = {5220}^{+58}_{-126}$}$ K agree well within their uncertainties. The indicated temperature difference (between the stellar photosphere and the spot) of about 300 K is consistent with literature values obtained for other main-sequence stars, as illustrated in Fig. 8 of Mancini et al. (2017).

4.3.3. WASP-41b transmission spectrum

In Fig. 16 we show the wavelength dependent planet-to-star radius ratios, which we derive from the simultaneous analyses assuming a spot and a spot-free model. Our results are compatible with a flat transmission spectrum. We compare our data to a model transmission spectrum computed with the Python Radiative-transfer in a Bayesian framework package (Pyrat Bay, Cubillos et al., in prep.), which is based on Blecic (2016) and Cubillos (2016). The model assumes a solar-abundance atmosphere in thermo-chemical equilibrium (Blecic et al. 2016), hydrostatic equilibrium, and an isothermal temperature profile (at the WASP-41b equilibrium temperature) for the system parameters given in Table 5. The radiative-transfer calculation considers opacities from Na and K (Burrows et al. 2000), H₂O and CO₂ (Rothman et al. 2010), collision-induced absorption from H₂-H₂ (Borysow et al. 2001; Borysow 2002) and H₂-He (Borysow et al. 1988, 1989; Borysow & Frommhold 1989), and H₂ Rayleigh scattering (Lecavelier Des Etangs et al. 2008). While a flat transmission spectrum can be interpreted as a sign of aerosols (clouds or hazes) in the planet’s atmosphere, our data do not have the necessary resolution to distinguish a cloudy from a cloud-free atmosphere.

Fig. 16

Transmission spectrum obtained from the simultaneous light curve analyses for two cases: the blue triangles show the radius ratios with their respective uncertainties from the analysis of all light curves assuming a spot-free model. The red circles present the radius ratios and uncertainties obtained through the analysis of all TLCs, taking starspots into account. We remark that for both scenarios, we excluded the TRAPPIST (2011/03/21) light curve from the analysis.

4.3.4. Impact of starspots on transit parameters

To investigate the impact of neglecting starspots in the analysis, we also model the spotted TLCs assuming they are spot-free. We further use these results to verify (by means of Bayesian Information Criterion (BIC) comparison) that models, which take into account starspots, have a higher probability than spot-free models. We then compare the spot and spot-free model results and find that the transit parameters, derived from the simultaneous analyses and, for most cases, from the individual analyses, agree well within 1.4σ. The transit depth values derived in the Gunn r′ (including Bessell R), the Bessell I, and the I + z band, assuming a spot-free model, agree well with those inferred when fitting for starspots. However, we find the following discrepancies for the individually analyzed EulerCam TLCs (2011/01/31, 2011/04/02) and the DFOSC I observation (2015/05/13):

The results from the spot-free analysis of the EulerCam TLC(2011/04/02) show differences in phase offset,a/R_S and limb darkening coefficients by <3σ, <2σ, and <5σ, respectively. We find similar discrepancies for the DFOSC I (2015/05/13) TLC affecting also the a/R_S value and the limb darkening coefficients.
The EulerCam (2011/01/31) TLC only shows differences in the limb darkening coefficients to an extent of <2σ.

The impact of stellar spots on the determination of LDCs, the relative semi-major axis a/R_S, and the transit timing and duration is a known issue, and has been investigated by several authors (e.g., Silva-Válio 2010; Ballerini et al. 2012; Tregloan-Reed et al. 2013; Oshagh et al. 2013b; Csizmadia et al. 2013). Especially starspots, which are located at the limb of the star, can bias the transit ingress (or egress) time of the planet. This directly affects the measurement of the total transit duration and hence, the relative semi-major axis a/R_S. It is thus important to account for spot features in TLC fitting processes to derive precise transit parameters, and when studying TTVs and transmission spectra.

5. Summary and conclusions

We have presented PyTranSpot, a routine designed to model TLCs with stellar activity features. PyTranSpot uses a pixellation approach to model the transiting planet, stellar limb darkening, and starspots on the stellar surface. We have merged PyTranSpot with the MCMC framework developed by Lendl et al. (2017) to derive accurate exoplanet transmission spectra in the presence of starspots and correlated instrumental noise.

We validated our routine by analyzing eleven synthetic light curves of four different star-planet configurations, including six synthetic light curves, which show anomalies due to occulted starspots. By comparing our derived results with the original system parameters, we found that PyTranSpot could reproduce the properties of our synthetic systems.

We further performed a multi-wavelength analysis of 20 TLCs of the well-studied system WASP-41b using archival and as yet unpublished data. From this data set, seven TLCs were affected by starspot occultations. We analyzed the light curves simultaneously for two cases (analysis of all TLC assuming a spot and a spot-free model) as well as individually. We found that our derived results agree well within one sigma with the values given in the literature. In our study, assuming a spot-free model for the spotted TLCs did not seem to have a significant impact on the determination of the transit depth. However, we identified that not taking into account stellar spots in the (individual) TLC analyses affected measurements of the LDCs, the relative semi-major axis a/R_S, and the transit midtime. These results confirm findings and predictions from various authors (e.g., Silva-Válio 2010; Sanchis-Ojeda et al. 2011; Ballerini et al. 2012; Tregloan-Reed et al. 2013; Oshagh et al. 2013b; Csizmadia et al. 2013; Zellem et al. 2017). For each simultaneous analysis, we additionally obtained three wavelength-dependent r_P/r_S values for WASP-41b covering a range of about 6200−9200 Å. We did not observe any significant variation of the transit depth with wavelength, however, our data are fully compatible both with a cloud-free and cloudy (i.e., flat) transmission spectrum. From the simultaneous multiband observation of the same starspot in the TLCs EulerCam and TRAPPIST (2011/04/02), we further derived a temperature difference between the stellar photosphere and the starspot of about 300 K, which is consistent with literature values for similar dwarf stars (e.g., Mancini et al. 2017).

We conclude that the outcome of our light curve analyses illustrates the importance of accounting for stellar activity features in TLCs for the correct interpretation of exoplanet transit parameters, TTVs, and transmission spectra. Furthermore, having simultaneous multiband observations of occulted starspots can help constrain a starspot’s temperature, disentangling the correlation between spot size and temperature.

¹

http://exoplanet.eu/catalog/ (Schneider et al. 2011).

²

A more detailed description of the transit method can be found in Winn (2010).

³

Researchers interested in PyTranSpot should contact the author.

⁴

https://github.com/pcubillos/MCcubed

⁵

http://wasp-planet.net

⁶

http://www.astro.keele.ac.uk/jkt/codes/jktld.html

Acknowledgments

We acknowledge the Austrian Forschungsförderungsgesellschaft FFG projects “RASEN” P847963. We also thank the anonymous referee for helping to improve the original manuscript. We thank contributors to Numpy (van der Walt et al. 2011), SciPy (Jones et al. 2001), Matplotlib (Hunter 2007), the Python Programming Language, and contributors to the free and open-source community. I.J. thanks M. Oshagh for fruitful discussions and J. Southworth for providing the corrected spot parameter. H. Lammer and M. Lendl also acknowledge support from the FWF project P25256-N27 “Characterizing Stellar and Exoplanetary Environments via Modelling of Lyman-α Transit Observations of Hot Jupiters”.

References

Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, MNRAS, 359, 567 [NASA ADS] [CrossRef] [Google Scholar]
Ballerini, P., Micela, G., Lanza, A. F., & Pagano, I. 2012, A&A, 539, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Barnes, S. A. 2007, ApJ, 669, 1167 [NASA ADS] [CrossRef] [Google Scholar]
Barros, S. C. C., Boué, G., Gibson, N. P., et al. 2013, MNRAS, 430, 3032 [NASA ADS] [CrossRef] [Google Scholar]
Beeck, B., Cameron, R. H., Reiners, A., & Schüssler, M. 2013a, A&A, 558, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Beeck, B., Cameron, R. H., Reiners, A., & Schüssler, M. 2013b, A&A, 558, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Béky, B., Kipping, D. M., & Holman, M. J. 2014, MNRAS, 442, 3686 [NASA ADS] [CrossRef] [Google Scholar]
Blecic, J. 2016, ArXiv e-prints [arXiv:1604.02692] [Google Scholar]
Blecic, J., Harrington, J., & Bowman, M. O. 2016, ApJS, 225, 4 [NASA ADS] [CrossRef] [Google Scholar]
Borysow, A. 2002, A&A, 390, 779 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Borysow, A., & Frommhold, L. 1989, ApJ, 341, 549 [NASA ADS] [CrossRef] [Google Scholar]
Borysow, J., Frommhold, L., & Birnbaum, G. 1988, ApJ, 326, 509 [NASA ADS] [CrossRef] [Google Scholar]
Borysow, A., Frommhold, L., & Moraldi, M. 1989, ApJ, 336, 495 [NASA ADS] [CrossRef] [Google Scholar]
Borysow, A., Jorgensen, U. G., & Fu, Y. 2001, J. Quant. Spec. Radiat. Transf., 68, 235 [NASA ADS] [CrossRef] [Google Scholar]
Burrows, A., Marley, M. S., & Sharp, C. M. 2000, ApJ, 531, 438 [NASA ADS] [CrossRef] [Google Scholar]
Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M. 2000, ApJ, 529, L45 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377 [NASA ADS] [CrossRef] [Google Scholar]
Csizmadia, S., Pasternacki, T., Dreyer, C., et al. 2013, A&A, 549, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cubillos, P. E. 2016, ArXiv e-prints [arXiv:1604.01320] [Google Scholar]
Cubillos, P., Harrington, J., Loredo, T. J., et al. 2017, AJ, 153, 3 [NASA ADS] [CrossRef] [Google Scholar]
Czesla, S., Huber, K. F., Wolter, U., Schröter, S., & Schmitt, J. H. M. M. 2009, A&A, 505, 1277 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Désert, J.-M., Charbonneau, D., Demory, B.-O., et al. 2011a, ApJS, 197, 14 [NASA ADS] [CrossRef] [Google Scholar]
Désert, J.-M., Sing, D., Vidal-Madjar, A., et al. 2011b, A&A, 526, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fabrycky, D. C., & Winn, J. N. 2009, ApJ, 696, 1230 [NASA ADS] [CrossRef] [Google Scholar]
Ferraz-Mello, S., Tadeu dos Santos, M., Folonier, H., et al. 2015, ApJ, 807, 78 [NASA ADS] [CrossRef] [Google Scholar]
Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419 [NASA ADS] [CrossRef] [Google Scholar]
Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457 [Google Scholar]
Gillon, M., Lanotte, A. A., Barman, T., et al. 2010, A&A, 511, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gillon, M., Jehin, E., Magain, P., et al. 2011, in Eur. Phys. J. Web Conf., 11, 06002 [Google Scholar]
Herrero, E., Ribas, I., Jordi, C., et al. 2016, A&A, 586, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
Jehin, E., Gillon, M., Queloz, D., et al. 2011, The Messenger, 145, 2 [NASA ADS] [Google Scholar]
Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, http://www.scipy.org/ [Google Scholar]
Kipping, D. M. 2012, MNRAS, 427, 2487 [NASA ADS] [CrossRef] [Google Scholar]
Kipping, D. M. 2013, MNRAS, 435, 2152 [Google Scholar]
Kopal, Z. 1950, Harvard College Observatory Circular, 454, 1 [Google Scholar]
Lecavelier DesEtangs, A., Pont, F., Vidal-Madjar, A., & Sing, D. 2008, A&A, 481, L83 [NASA ADS] [CrossRef] [Google Scholar]
Lendl, M., Anderson, D. R., Collier-Cameron, A., et al. 2012, A&A, 544, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lendl, M., Cubillos, P. E., Hagelberg, J., et al. 2017, A&A, 606, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Madhusudhan, N., & Seager, S. 2009, ApJ, 707, 24 [NASA ADS] [CrossRef] [Google Scholar]
Mancini, L., Southworth, J., Ciceri, S., et al. 2014, MNRAS, 443, 2391 [NASA ADS] [CrossRef] [Google Scholar]
Mancini, L., Southworth, J., Raia, G., et al. 2017, MNRAS, 465, 843 [NASA ADS] [CrossRef] [Google Scholar]
Maxted, P. F. L. 2016, A&A, 591, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Maxted, P. F. L., Anderson, D. R., Collier Cameron, A., et al. 2011, PASP, 123, 547 [NASA ADS] [CrossRef] [Google Scholar]
Miller-Ricci, E., Seager, S., & Sasselov, D. 2009, ApJ, 690, 1056 [NASA ADS] [CrossRef] [Google Scholar]
Montalto, M., Boué, G., Oshagh, M., et al. 2014, MNRAS, 444, 1721 [NASA ADS] [CrossRef] [Google Scholar]
Neveu-VanMalle, M., Queloz, D., Anderson, D. R., et al. 2016, A&A, 586, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Nutzman, P. A., Fabrycky, D. C., & Fortney, J. J. 2011, ApJ, 740, L10 [NASA ADS] [CrossRef] [Google Scholar]
Oshagh, M., Boisse, I., Boué, G., et al. 2013a, A&A, 549, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Oshagh, M., Santos, N. C., Boisse, I., et al. 2013b, A&A, 556, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pollacco, D. L., Skillen, I., Collier Cameron, A., et al. 2006, PASP, 118, 1407 [NASA ADS] [CrossRef] [Google Scholar]
Pont, F., Gilliland, R. L., Moutou, C., et al. 2007, A&A, 476, 1347 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
Prša, A., & Zwitter, T. 2005, ApJ, 628, 426 [NASA ADS] [CrossRef] [Google Scholar]
Prša, A., Conroy, K. E., Horvat, M., et al. 2016, ApJS, 227, 29 [NASA ADS] [CrossRef] [Google Scholar]
Rabus, M., Alonso, R., Belmonte, J. A., et al. 2009, A&A, 494, 391 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Rackham, B., Espinoza, N., Apai, D., et al. 2017, ApJ, 834, 151 [NASA ADS] [CrossRef] [Google Scholar]
Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spectr. Rad. Transf., 111, 2139 [NASA ADS] [CrossRef] [Google Scholar]
Sanchis-Ojeda, R., & Winn, J. N. 2011, ApJ, 743, 61 [NASA ADS] [CrossRef] [Google Scholar]
Sanchis-Ojeda, R., Winn, J. N., Holman, M. J., et al. 2011, ApJ, 733, 127 [NASA ADS] [CrossRef] [Google Scholar]
Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schwarz, G. 1978, Ann. Stat., 6, 461 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Silva, A. V. R. 2003, ApJ, 585, L147 [NASA ADS] [CrossRef] [Google Scholar]
Silva-Válio, A. 2010, in Solar and Stellar Variability: Impact on Earth and Planets, eds. A. G. Kosovichev, A. H. Andrei, & J.-P. Rozelot, IAU Symp., 264, 440 [Google Scholar]
Southworth, J. 2008, MNRAS, 386, 1644 [NASA ADS] [CrossRef] [Google Scholar]
Southworth, J., Tregloan-Reed, J., Andersen, M. I., et al. 2016, MNRAS, 457, 4205 [NASA ADS] [CrossRef] [Google Scholar]
terBraak, C. J. F., & Vrugt, J. A. 2008, Stat. Comput., 18, 435 [CrossRef] [Google Scholar]
Thaler, I., & Spruit, H. C. 2014, A&A, 566, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Tregloan-Reed, J., Southworth, J., & Tappert, C. 2013, MNRAS, 428, 3671 [NASA ADS] [CrossRef] [Google Scholar]
Tregloan-Reed, J., Southworth, J., Burgdorf, M., et al. 2015, MNRAS, 450, 1760 [NASA ADS] [CrossRef] [Google Scholar]
van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [CrossRef] [Google Scholar]
Wilson, R. E. 1979, ApJ, 234, 1054 [NASA ADS] [CrossRef] [Google Scholar]
Wilson, R. E. 1990, ApJ, 356, 613 [NASA ADS] [CrossRef] [Google Scholar]
Wilson, R. E. 2008, ApJ, 672, 575 [NASA ADS] [CrossRef] [Google Scholar]
Wilson, R. E. 2012, AJ, 144, 73 [NASA ADS] [CrossRef] [Google Scholar]
Wilson, R. E., & Devinney, E. J. 1971, ApJ, 166, 605 [NASA ADS] [CrossRef] [Google Scholar]
Winn, J. N. 2010, ArXiv e-prints [arXiv:1001.2010] [Google Scholar]
Winn, J. N., Holman, M. J., Torres, G., et al. 2008, ApJ, 683, 1076 [NASA ADS] [CrossRef] [Google Scholar]
Winn, J. N., Fabrycky, D., Albrecht, S., & Johnson, J. A. 2010a, ApJ, 718, L145 [NASA ADS] [CrossRef] [Google Scholar]
Winn, J. N., Johnson, J. A., Howard, A. W., et al. 2010b, ApJ, 723, L223 [NASA ADS] [CrossRef] [Google Scholar]
Zellem, R. T., Swain, M. R., Roudier, G., et al. 2017, ApJ, 844, 27 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Synthetic light curve residuals

Fig. A.1

(O−C) residuals (a) and the difference between original and best-fit light curve model (b) of SYNTH-1. The obtained photometric parameters can be found in Table B.1.

Fig. A.2

(O−C) residuals (a) and the difference between original and best-fit light curve model (b) of system SYNTH-2. The SYNTH-2 light curve shows an occulted starspot at the limb of the star around phase –0.01. The obtained photometric and spot parameters can be found in Tables B.1 and B.3.

Fig. A.3

(O−C) residuals of the synthetic light curves SYNTH-3a to SYNTH-3e. SYNTH-3e shows a starspot anomaly located close to the limb of the star around phase –0.005. The obtained photometric and spot parameters can be found in Tables B.1 and B.3.

Fig. A.4

Differences between the original and best-fit light curve models of SYNTH-3a, -3b, -3c, and -3e.

Fig. A.5

(O−C) residuals of the light curves SYNTH-4a to SYNTH-4d. All TLCs show the same transit event and occulted starspot around phase –0.002. The obtained photometric and spot parameters can be found in Tables B.2 and B.4.

Fig. A.6

Differences between the original and best-fit light curve models of SYNTH-4a, -4b, -4c, and -4d.

Appendix B: Derived transit and spot parameters

Table B.1

Derived best-fit values resulting from the individual (SYNTH-1a and SYNTH-2a) and simultaneous (SYNTH-3a to SYNTH-3e) analysis.

Table B.2

Derived best-fit transit parameters resulting from the simultaneous analysis of SYNTH-4a to SYNTH-4d.

Table B.3

Derived spot parameters for the synthetic light curves SYNTH-2a, SYNTH-3e.

Table B.4

Derived spot parameters for the synthetic light curves SYNTH-4a to SYNTH-4d.

Table B.5

Derived best-fit parameters of the simultaneous multiband WASP-41b TLC analysis, assuming a spot-free model for all TLCs.

Table B.6

Derived best-fit parameters of the simultaneous multiband transit and starspot LC analysis of WASP-41b.

Table B.7

Derived spot parameters for the simultaneously analyzed WASP-41b light curves.

Table B.8

Derived spot parameters for the individually analyzed WASP-41b light curves.

All Tables