Free Access
Volume 635, March 2020
Article Number A6
Number of page(s) 17
Section Numerical methods and codes
Published online 28 February 2020

© ESO 2020

1. Introduction

A few decades have passed since the discovery of the first exoplanet orbiting the main-sequence (MS) star 51 Pegasi (Mayor & Queloz 1995) and by now ∼4000 exoplanets have been confirmed, ∼75% of which are transiting planets1.

In the initial years the most effective technique for discovering exoplanets was the radial velocity (RV) method (see e.g. Wright 2018), but during the last decade the transit method has become the most prominent, also thanks to past space missions, such as Kepler (Basri et al. 2005; Borucki et al. 2010), its extension K2 (Howell et al. 2014), and CoRoT (Rouan et al. 1999; Auvergne et al. 2009). In the meantime, other space telescopes, such as TESS (Sharma et al. 2018) have just begun detecting exoplanetary transits; CHEOPS (Broeg et al. 2013) has just been launched, whereas PLATO (Rauer et al. 2014) is expected to be launched in 2026.

As is well known, once the light curve (LC) of a planet-hosting star is available, the transit depth gives the squared ratio of planetary to stellar radii , while the orbital period P and the transit shape parameters (impact parameter, depth, and duration) enable the retrieval of the stellar mean density ρ. This is true for circular orbits under the assumption that exoplanetary mass Mp is negligible if compared with stellar mass M; otherwise RV measurements are needed to constrain the orbital eccentricity e and the planetary mass Mp as well. A review of the transit technique is provided in Winn (2010).

From these considerations it is clear that establishing the exoplanetary parameters (in particular, retrieving Rp in case we are just dealing with the transit method) requires not only the availability of transit photometry, but also knowledge of stellar parameters, either M or R, as pointed out by Seager & Mallén-Ornelas (2003). Thus, as a very basic example, LC analyses are usually split into three steps: first the LC-analyser tool recovers ρ, after that ρ together with the stellar effective Teff and metallicity [Fe/H] are employed to infer M and/or R from stellar evolutionary models, and finally the transit-analysis tool is launched again assuming also theoretical stellar parameters as inputs to determine Rp (see e.g. Gillon et al. 2009a).

In this paper we present our custom-developed Markov chain Monte Carlo + isochrones (MCMCI) FORTRAN code, which makes the transit analysis algorithm directly interact with stellar evolutionary models, so that starting from LCs and very basic stellar parameters, it is possible to characterise the whole exoplanetary system directly. Our MCMCI tool is very useful to carry out this kind of analysis and it will be very valuable, especially in coming years when lots of LCs and data from transits are expected to be available as a result of already active and forthcoming missions.

Our unified approach of simultaneously modelling the star and its planets is remarkable even if it is not unique, since it has already been followed, for example by Hartman et al. (2019), who carry out an isochrone-based joint analysis using a differential evolution MCMC procedure (Ter Braak 2006) and PARSEC stellar evolutionary models (Marigo et al. 2017). Also Siverd et al. (2012) simultaneously model the stellar host and its companion, but to break the MR degeneracy they simply consider the mass-radius relation by Torres et al. (2010), as they use the first version of the IDL fitting package EXOFAST (Eastman et al. 2013). We notice that our version of MCMC without the support of stellar evolutionary models also allows us to select a modified version of the Torres et al. (2010)MR law from well-constrained detached binary systems (Gillon et al. 2011); thus we can use this tool during the LC or RV fitting.

Recently, a new version of the EXOFAST package called EXOFASTv2 has been released (Eastman et al. 2019): it facilitates fitting the stellar properties along with the planetary fit with MIST evolutionary tracks (Dotter 2016) or Yonsei-Yale (YY) stellar evolutionary models (Yi et al. 2001), following the same philosophy as our MCMCI. However an optional spectral energy distribution (SED) fitting is also included, which our code does not perform. Moreover, besides jointly dealing with photometric and RV time series similar to our MCMCI, modelling stellar and planetary astrometric signals, modelling Doppler tomography (Collier Cameron et al. 2010), and integrating a mass-radius relation for exoplanets are remarkable features of EXOFASTv2 that are not implemented in MCMCI.

Other codes are available to address both the transit and RV fitting, such as TLCM (Smith et al. 2017), PlanetPack3 (Baluev 2018), juliet (Espinoza et al. 2019), Pyaneti (Barragán et al. 2019), exoplanet2, allesfitter (Günther & Daylan 2019). The TLCM code is written in IDL, PlanetPack3 is written in C++, while the other codes are written in PYTHON (Pyaneti is also composed of PORTRAN routines, which are then wrapped to PYTHON). A notable feature of TLCM is its possibility of simultaneously modelling RV of both components of a binary star, and it incorporates ellipsoidal modulation (Kopal 1959; Morris 1985) and a beaming effects treatment (Rybicki & Lightman 1979; Loeb & Gaudi 2003) in the LCs. The Planet Pack3 tool includes Gaussian processes (Rasmussen & Williams 2005) to model the noise in RV data, but does not allow us to fit multiple planets in the photometry. The juliet, exoplanet, and allesfitter codes account for Gaussian processes in both LC and RV time series using celerite (Foreman-Mackey et al. 2017). Moreover, the core of exoplanet is the implementation of an Hamiltonian MCMC procedure (see e.g. Neal 2012; Betancourt 2017) that promises faster convergence time with respect to traditional MCMC implementation, while alles fitter also enables us to model star spots and stellar flares using aflare (Davenport et al. 2014).

Apart from EXOFASTv2, all the other listed codes do not include a simultaneous analysis of the stellar host through stellar evolutionary models along the LC and RV fitting. Our MCMCI, instead, allows this isochrone-based joint analysis and may also give different constraints for the stellar age (and thus the age of the entire exoplanetary system) using both model-dependent and empirical age indicators. We finally notice that our tool is unique in that it is fully implemented in FORTRAN.

In Sect. 2 an overview of the code is given, in Sect. 3 its application on two different exoplanetary systems is presented, and Sect. 4 reports our conclusions.

2. Code description

The MCMCI is a FORTRAN program that was born from a proper merging of the custom MCMC code widely presented in Gillon et al. (2010, 2012) and the isochrone placement algorithm, which is described in Bonfanti et al. (2015, 2016). In the following subsections we recall key aspects of the two codes according to their most recent updates and we explain how we merged these tools to facilitate their interaction.

2.1. Markov chain Monte Carlo

The MCMC simulation is a stochastic process (see e.g. Holman et al. 2006; Ford 2005) that generates a sequence (chain) of data points (states) starting from an initial state that is then perturbed. The aim of this simulation is to sample the probability distribution function (PDF) of some parameters of interest assuming a given model; our version makes use of the Metropolis-Hastings algorithm (Hastings 1970) and the Gibbs sampler (Casella & George 1992, and references therein).

It is likely that the first states of a chain do not come from the limiting distribution, therefore a common practice in the MCMC approach is to discard these first data points: this process is called burn-in. In this way the effect of initial values on the posterior inference is minimised. Choice of the burn-in length depends upon the initial state and the speed of convergence to the limiting distribution. Establishing the proper burn-in length requires analysis of the output case by case, however experience suggests that setting the burn-in length to 20% represents a conservative compromise. This is the value that has been chosen for all the analysis described in Sect. 3. Anyway, our input form enables us to manually specify the length of the burn-in phase.

2.1.1. Models and input parameters

Our code may deal with any number of LC and RV time series for a complete joint analysis of the transit and dynamics of the exoplanetary system. If only LCs are available, only information recoverable from photometry is retrieved (e.g. the planetary radius Rp), and a determination of planetary mass Mp is not possible unless the RV semi-amplitude K is known from some RV studies. Vice versa, if only RV time series are available, the exoplanetary system are characterised from a dynamic point of view, that is Mp is obtained, but not Rp. Our implementation of the code assumes use of the photometric model by Mandel & Agol (2002) to reproduce the eclipse; and a classical Keplerian model for analysis of the RV signal, in addition to a Rossiter-McLaughlin effect model (Giménez 2006) if RVs are obtained during transit. Our global model may deal with any number of planets, either transiting or not.

The eclipse model is multiplied by a trend model aimed to reproduce all the systematic effects having either astrophysical or instrumental origin, which are responsible for photometric variations beyond those caused by the transit itself. Specifically, through polynomial functions, this baseline is able to model, for example, the effect of inhomogeneous intra-pixel sensitivity of the detector (Knutson et al. 2008; Charbonneau et al. 2008); the so-called ramp effect, according to which the detector gain may increase asymptotically over time, depending on the illumination history of the pixels (Knutson et al. 2008; Charbonneau et al. 2008); the time-dependent photometric modulation of the stellar signal due to rotating spots; and the effect of the variations of the sky background, the airmass, or the width of the point spread function during observation.

The code is also able to face trends in the RV time series in terms of time, cross-correlation function width and bisector parameters (Baranne et al. 1996), and (a typical indicator of stellar magnetic activity; see e.g. Wright et al. 2004, and references therein).

Barycentric Dynamic Time (TDB) is the default time standard that is used by the code in all the analyses; TDB is suitable for precision time monitoring. This time standard is preferable to the Coordinated Universal Time (UTC) because it avoids the drift due to leap seconds; in addition, it refers to the barycentre of the Solar System, thus it corrects for relativistic effect due to the gravitational potential of Earth. A discussion about time standards and timing precision is provided in for example Eastman et al. (2010).

The following set of parameters may be randomly perturbed at each chain step (jump parameters or step parameters).

  • The transit depth dF.

  • The occultation depth dFocc.

  • The impact parameter in the case of circular orbit b′.

  • The eclipse duration W, that is the time between the first and last contact, as defined in Winn (2010).

  • The eclipse timing T0, that is the time of inferior conjunction, when the true anomaly at transit is ft = 90° −ω (with ω argument of the periastron).

  • The orbital period P.

  • The quantities and , where e is the orbital eccentricity and ω is the argument of the periastron.

  • The parameter , where K is the RV semi-amplitude.

  • The quantities and , where v sin i is the stellar rotational velocity along the line of sight and β is the projected angle between the stellar spin axis and orbital axis.

  • The stellar metallicity [Fe/H].

  • The stellar effective temperature Teff or colour index.

  • The stellar radius R.

  • The stellar mass M.

  • The limb-darkening (LD) coefficients.

If more than one planet has to be fitted, all the relevant parameters of each planet are listed in the input form.

The algorithm checks whether the drawn values of the jump parameters are physical. Therefore, for instance, negative values for dF, dFocc, b′, W, T0, P, K2, Teff, R or M are not allowed, any square root arguments must be positive, and the inequalities , must hold since we are dealing with bound orbital systems. If any of the drawn value is not physical, the jump is not accepted, and a copy of the previous state is made.

2.1.2. Comments on jump parameters

Stepping in dF, that is assuming a uniform prior on dF, creates an implicit bias on towards higher Rp values. If the slope of the Rp prior is flat enough over the range of interest of the posterior, this does not unphysically bias the posterior, however an alternative approach would be to replace the stepping parameter dF with , as proposed for example by Eastman et al. (2013).

The impact parameter in case of circular orbit is given by b′=a cos ip/R, where a is the semi-major axis of planetary orbit and ip is the orbit inclination with respect to the plane of the sky. As pointed out by Gillon et al. (2009a), assuming b′ rather than the actual impact parameter is preferable to minimise the correlation between jump parameters.

Our algorithm can also handle eccentric orbits by enabling the stepping in and . This parametrisation ensures both orthogonality and uniform priors on e and ω, which has been pointed out for the first time by Anderson et al. (2011) and then broadly discussed by Eastman et al. (2013). In case RV data are available, the usual approach is to launch a first run with a fixed e = 0 value, and then launch another run letting e as a free jump parameter and check whether the BIC (Eq. (5), see later) is in favour of an eccentric or circular solution.

Instead, in case of a transit-only fit, treating e is generally problematic. On the one hand, orbits of close-in exoplanets are likely circular. Dawson & Johnson (2018) report that for orbital periods P <  3 days (i.e. orbital semi-major axes a <  0.04 AU, if a Sun-like star is assumed as host), hot Jupiters are consistent with e = 0 orbits, while instead in the 3 <  P <  10 range (in days, i.e. 0.04 AU <  a <  0.09 AU) it is common to observe moderately elliptical orbits. This is consistent with the circularisation process induced by tidal evolution: the shortest the orbital period, the fastest the circularisation timescale. Moreover, the less massive the planet, the fastest the circularisation timescale. Considering an initial semi-major axis a0 = 0.04 AU and an initial orbital eccentricity e0 = 0.1, the orbital circularisation happens after ∼150 Myr, ∼100 Myr, and ∼30 Myr for a hot Jupiter, a hot Neptune, and a super-Earth, all orbiting a Sun-like star (Rodríguez et al. 2010, 2011). As a consequence, less massive planets are expected to exhibit e = 0 orbits at larger host star separation with respect to more massive planets (Pont et al. 2011). Summing up, it is reasonable to assume circular orbits for semi-major axes lower than a few hundredths of astronomical units.

On the other hand, if we are dealing with outer exoplanets, setting e = 0 may cause some biases on the derived parameters of the system. In this case, we could infer the mean stellar density as a result of evolutionary models (say ρ⋆, th its value) and use it together with the transit duration trying to constrain e. To do so, we should set ρ⋆, th as a prior in our MCMCI input form and let P, dF, W, b, e, and ω as jump parameters. The transit-based parameter, which is computed by the code through Eq. (A.4), is translated to ρ through Kepler’s third law and then prior ρ⋆, th drives the ρ value. Thus, favoured values of the transit parameters are those which produce a ρ value more similar to ρ⋆, th. The (e, ω) degeneracy may not lead to a complete convergence of the transit parameters, however, as shown by Dawson & Johnson (2012), it is possible to identify highly eccentric Jupiter-sized planets because the LC alone gives a high lower limit on e. Therefore, if we are studying Jupiter analogues not very close to their hosts (so that an eccentric orbit cannot be excluded a priori), it is worth following the procedure we have just described. Since our algorithm produces many files which list all the values assumed by the jump parameters step by step, it is then possible to produce a plot of e versus ω and compare it with the patterns that are presented in Fig. 2 of Dawson & Johnson (2012). A lower limit emin on the eccentricity can be set if the observed pattern resembles a pattern for which ; that is ρ⋆,th sensibly differs from ρ⋆,circ, which is the mean stellar density that would be derived from the transit fit assuming e = 0. Dawson & Johnson (2012) also note that assuming e = emin implies a transit at periastron (resp. apoastron) if ρ⋆,th <  ρ⋆,circ (resp. ρ⋆,th >  ρ⋆,circ). Thus, launching a second MCMCI run where fixed e = emin and ω = 90° (or 270° according to the case) are assumed may facilitate convergence and reduce possible biases in the transit parameters. Of course, awareness that the claimed e value is just a lower limit and that ω value comes out as a consequence of that must be kept. Instead, if the observed pattern in the e-ω plot is similar to the g(e, ω) = 1 scenario, although any e value could be consistent with the transit observable in principle, Dawson & Johnson (2012) correctly explain that the planetary orbit likely exhibits a low eccentricity. In this case performing a transit-only fit with the assumption e = 0 is reasonable.

In case an e = 0 orbit is assumed, ω (which is undefined) is set equal to 0°: the “periastron” is placed at the ascending node, thus the time of “periastron” is the reference time when RV is maximum if the orbit is circular. Even if exoplanet archives do not have an unambigous convention about the ω value to adopt when e = 0, however we call out that the most popular standard that is adopted in the literature is ω = 90° (see e.g. Eastman et al. 2013; Iglesias-Marzoa et al. 2015; Kreidberg 2015).

The parameter K2 is used as a jump parameter instead of K to minimise the correlation between jump parameters (Gillon et al. 2009a). According to its definition, we note that , that is it is independent of e and P. In particular, for a given stellar mass, assuming a uniform prior on K2 implies assuming a uniform prior on Mp sin ip.

In case of a multi-planetary system, so far only the β value of one single planet is handled by the code. We note that β can be inferred from the Rossiter-McLaughlin effect and it is not so common to observe this effect on multiple exoplanets belonging to the same system, however we leave the update to multiple values of β to a future release of the code.

Stepping in R is optional. This option is used for transits too shallow or noisy to make a precise determination of their impact parameter possible; in this case, a/R and ρ cannot be well constrained either.

Limb-darkening coefficients can be either four (not-linear model, nl) or two (quadratic model, qd). In case the qd model is used, the assumed jump parameters are not the direct LD coefficients u1 and u2, but the linear combinations c1 = 2u1 + u2 and c2 = u1 − 2u2 as in Holman et al. (2006) to minimise the correlation within the uncertainties of LD parameters. For each specified filter, LD coefficients are inferred through proper interpolation in the tables of Claret & Bloemen (2011). Owing to systematics in these tables, we emphasise that simply fixing the LD coefficients at the table values may bias the transit fit or may work backward to bias the log g, Teff or [Fe/H] values from which the LD coefficients are derived. To avoid this scenario, especially in case of high quality photometry, we strongly recommend adding a Bayesian penalty (see later, Eq. (3)) for the LD coefficients with a prior from the tables. Actually, this is the default behaviour of our code and it has been followed for all the analyses described in the paper.

2.1.3. Merit function and its role

The set of input parameters affects both the photometric model reproducing the LC and the RV model reproducing the RV time series (if any); thus at each step we can define


where nLC is the number of available LCs, npi are the number of points of the ith LC, fi, and σfi are the vectors of flux measurements and its errors of the ith LC, while are the fluxes that have been computed according to the eclipse model in use. Adopting the same kind of notation, we can also introduce the analogous referring to the RV time series


Summing together and , we obtain the basic merit function. In addition, the jump parameters specified above – besides the stellar luminosity L, surface gravity log g, density ρ, and projected rotational velocity v sin i – may also be set as priors. Those parameters that have been set as Gaussian priors enter the merit function, as well, in the form of the following Bayesian penalty addendum:


where and are the values with their uncertainties that are attributed to the generic prior parameter, while x comes from the perturbation of at the generic step of the MCMC. In particular, in all the analyses we carried out, we always considered BPs on c1 and c2 (computed from the u1 and u2 LD coefficients). Inference of the initial values for u1 and u2 (with their respective errors) from Claret’s tables follows the same procedure as described in Gillon et al. (2009a). After that, BPc1 and BPc2 control the u1 and u2 floating during the MCMC steps, to obtain a LD solution consistent with theory.

Finally, at each step, it is possible to define the complete merit function, which is given by


where npp is the number of prior parameters.

Once the priors and the jump parameters have been set, the jump towards the new state is accepted if χ2 at that step is lower than the one computed at the previous step. If this does not hold, the new state is accepted with probability , where Δχ2 is the χ2 difference between the two last steps. In all the other cases in which the new state is not accepted, the state is pushed back to the previous set of parameters.

Generation of a new state x from the present state x is controlled by a transition probability function and the optimal choice for this function would be considering the posterior probability distribution. However, the aim of the MCMC is actually to retrieve the posterior distribution itself, which is unknown at the beginning. Thus, as a common choice, we considered a Gaussian distribution centred around x. To make the MCMC efficiently converge, the step size, that is the Gaussian variance, to perform the jump from x to x must be tuned up. This is done by computing the acceptance rate, that is the fraction of accepted states over a window spanning a given number of steps, which we set at 100. As reported by Ford (2005), if x has only one dimension, the optimal acceptance rate is ∼0.44, while if x has many dimensions, the optimal acceptance rate is ∼0.25, as has been proven by Roberts et al. (1997).

We implemented the Gibbs sampling during the burn-in phase, when just a single jump parameter changes at each step; thus, for tuning-up the step size, during the burn-in phase the acceptance rate is set to 40% if the Gibbs sampler is switched on. After the burn-in phase, all the jump parameters are modified at a time, therefore in this phase we monitor the acceptance rate imposing an optimal value equal to 25%.

If several sets of models that may differ in terms of input parameters or baseline choice have to be compared, the adopted criterion to choose the best model is given by the Bayesian information criterion (Schwarz 1978):


where is the smallest χ2 found in the Markov chains, k is the number of free parameters of a given model, and N is the number of data points. The smaller the BIC, the better the model.

At the end, the PDFs of the relevant stellar and planetary parameters are built and meaningful statistics can be retrieved.

2.1.4. Convergence discussion

It is worth launching several chains for a given process and then checking their mutual convergence thanks to the test by Gelman & Rubin (1992), which we briefly present following Ford (2006). If Nc is the number of chains, Lc is the length of each chain, and zic is the ith draw of the generic z parameter at the cth chain, the mean z value within the cth chain is written as


while the global mean of z considering the single mean values computed from each chain is written as


It is then possible to compute W(z) as the average of the variances within each chain, and B(z) as the variance based on the mean values of each chain, that is



An unbiased estimator of the variance of z is then given by the following weighted average between W(z) and B(z):


The Gelman-Rubin statistic is commonly denoted by , the estimator is defined as


and convergence is reached when is close to 1. Different threshold values indicating convergence are adopted in the literature: randomly picking up 100 papers that cite Gelman & Rubin (1992) in 2017, Vats & Knudson (2018) find cut-off values varying from 1.003 to 1.3. Among these papers, 43% of the authors (the relative majority) consider , which is the recommended value by Gelman & Rubin (1992).

Our own empirical practice for convergence is considering . This criterion may be easily satisfied in case of one single planet fitting, while a slight tension among jump parameters may arise in case of a multi-planetary system. To check convergence, Ford (2006) also recommends a minimum number of effective independent draws, which can be estimated through


with the estimator .

Since , from (11) we infer


The quantity comes from an average between W(z) and B(z), therefore its value is between W(z) and B(z); in particular, given condition (13), we have


Condition (14) implies


By recovering B(z) from Eq. (10) and substituting it in Eq. (15), we can express as a function of as follows:


where we used the definition of given by (11) in the last passage. By inverting (16), we can finally express as a function of , that is


To give some numbers, if we consider Nc = 5 and Lc = 80 0003, achieving the goal suggested by Ford (2006) () would require . According to our experience, this cut-off appears too demanding in general. We anticipate that this condition is satisfied for all the jump parameters we considered in the analysis of WASP-4 (Sect. 3.2 , star orbited by only one hot Jupiter), but it holds just for a few jump parameters of the HD 219134 system (Sect. 3.1), where the host is orbited by two close-in super Earths (detected both through transit and RV) and by another two, or possibly more, massive planets, whose presence has been inferred from the RV technique. According to all these considerations, we prefer to keep , as a general reference threshold to establish convergence, and we notice that this cut-off condition is stronger than the average practice in the literature.

2.1.5. Noise treatment

As explained in Gillon et al. (2012), to estimate the level of noise correlation in each photometric time series, we can define βw as the ratio of the standard deviation of the residuals (observed minus computed (O–C) values) to the mean photometric error: it represents the under- or over-estimation of the white noise. Furthermore, we can take the effect of red noise into account (Pont et al. 2006) by binning the LC on different timescales, and for each binning we then compute


where Nbin is the mean number of points in each bin, M the number of bins, and σ1 and σN are the standard deviations of the unbinned and binned residuals, respectively. The maximum among the available βr, bin values is set as the reference βr. Then, the photometric errors provided by the observations are multiplied by the correction factor


and eventually, for RVs, a “jitter” noise may be added quadratically to the error budget. After that a new MCMC run is launched with the updated (rescaled) error bars to obtain reliable error bars on the fitted parameters.

2.2. Isochrone placement

The isochrone placement technique aims to infer stellar parameters – above all radius R, mass M, and age t – by interpolating a set of observational input parameters in a grid of stellar evolutionary models; these observational input parameters come from spectroscopy, transit constraints, photometry and, in particular, they may also involve Gaia parallaxes and magnitudes. The grids of theoretical models in use are currently based upon the PARSEC v.1.2S (Chen et al. 2014; Bressan et al. 2012) + COLIBRI PR16 (Marigo et al. 2013; Rosenfield et al. 2016) stellar isochrones and tracks, which are globally presented in Marigo et al. (2017) and available for download through the CMD web interface4.

2.2.1. Tables of isochrones and evolutionary tracks

The code is designed so that the stellar location is compared with that of the isochrones, that is loci of evolutionary states occupied by coeval stars. Each grid of isochrones that has been downloaded from the CMD interface is labelled through the initial metallicity Zini of a star at its birth and contains several isochrones. In particular, we selected an age range satisfying the 6 <  log t <  10.1 (with t in billion years, thereby covering all the plausible stellar age values from pre-MS to post-MS phase) at steps δt = 0.05, which guarantees a fine enough grid considering how slowly stellar parameters evolve and, as a consequence, how large the typical age error bars are. Relevant columns reported by the grids are age, mass, effective temperature, luminosity, surface gravity, and absolute magnitudes in a given photometric system; clearly further theoretical parameters as radius or mean stellar density can be easily computed. Summing up, each Zini grid of values is made of sub-tables, each of which is representative of an isochrone characterised by a specific age value; these are appended sequentially. Each isochrone is made of several rows and the rows represent theoretical stars of the same age characterised by different parameters, most importantly mass.

As described later in Sect. 2.2.4, isochrone placement also takes into account, for example the stellar evolutionary speed and microscopic diffusion phenomenon, for which Z varies with time. To perform these complementary tasks and extract the needed information, it also needs to deal also with tables of evolutionary tracks. Each track table (labelled through Mini, ⋆, Zini) lists the sequence of evolutionary states through which a star having an initial mass Mini, ⋆ and an initial metallicity Zini passes during its life. The rows of the table list the evolution of a specific star passing time in terms of, for example mass, radius, and effective temperature. In particular, evolutionary track tables provide the surface abundance of hydrogen and helium at each evolutionary stage, so that the temporal metallicity evolution of a star can be followed, which is essential to face element diffusion (ED).

Considering that the relation between Z and [Fe/H] is exponential, to guarantee linearly spaced [Fe/H] values, it is necessary to take Z values in geometric progression. In particular, to have isochrone grids differing at most by δ[Fe/H] = 0.05 dex, we downloaded all those grids with Z belonging to a geometric progression with common ratio 1.115 from 0.0005 up to 0.00945 and also all those grids corresponding to Z from 0.010 to 0.070 at steps of 0.001. This implies a [Fe/H] range from −1.48 dex to +0.66 dex.

2.2.2. Input and output parameters

Isochrone placement was originally implemented in MatLAB language as an autonomous and independent routine and then converted into FORTRAN to be ready to interact with the MCMC code described in Sect. 2.1; however it is not based on a MCMC procedure. In other words, isochrone placement is a fitting algorithm that enables us to retrieve t⋆,isoch, M⋆,isoch, and R⋆,isoch in agreement with the isochrones, once a useful set of input parameters is provided.

As a subroutine of the main MCMC, the input parameters of isochrone placement are represented by the relevant set of jump parameters that are generated at each step of the MCMC. The role of M⋆,isoch (and possibly R⋆,isoch) is to drive the jumps in M (and possibly R) as described in Sect. 2.3. At the end, PDFs of M, R and t⋆,isoch are built to extract their respective statistics.

As input parameters, isochrone placement requires the stellar metallicity [Fe/H] and necessarily one parameter among stellar effective temperature Teff and colour index. In addition, it is needed at least one further parameter among distance (observationally available e.g. from parallax) together with apparent magnitude, stellar surface gravity log g, stellar mean density ρ, and stellar radius R. Optionally one or more parameters among v sin i, stellar rotational period Prot, , and the abundance ratio of yttrium to magnesium [Y/Mg] are useful to better constrain the stellar properties.

In principle, colour index can come from any photometric system supposing that that photometric system is available among the evolutionary models; so far, the code is ready to accept Johnson-Cousin or Gaia photometric systems. Stellar surface gravity log g can be available from spectroscopy or via asteroseismology (see e.g. the scaling relation by Chaplin et al. 2014, which is implemented in our code). Stellar mean density ρ can be observationally available from the LC analysis if a star hosts a transiting planet (see e.g. Sozzetti et al. 2007) or via asteroseismic scaling relations such as that in Chaplin et al. (2014) that is implemented in our code. Finally stellar radius R may be retrieved from interferometric observations in case of nearby stars.

2.2.3. Possible removing of isochrone degeneracies

The role of the v sin i/Prot, and [Y/Mg] optional parameters is to give a preliminary age indication, that may help in disentangling degeneracies involving MS and pre-MS isochrones, that is isochrone overlap on the Hertzsprung-Russell diagram (HRD). In the case in which both young and old isochrones are geometrically close to the stellar location even a rough age constraint may exclude unlikely very young ages. From each of these optional parameters, an age estimation can be retrieved through empirical relations. It is enough that a MS age is suggested in order to discard unlikely pre-MS isochrones. If this is the case, from a computational point of view, the following interpolation within isochrones considers only those sub-tables of isochrones not containing extremely young ages.

It is well known that rotation and magnetic activity decrease with age as shown by Skumanich (1972). The physical explanation originates in the presence of a convective envelope, thus any relation linking age with Prot (gyrochronology) or (magnetic activity) is to be applied just to late-type stars. This is the typical research domain of the exoplanetary field, so no particular caution need be taken by the user in this respect.

It is also important to stress that gyrochronological and activity-age relations become less sensible to age variations in mid-to-old stars. In fact, as stated by Soderblom (2010), it is difficult to detect and calibrate the decline of chromospheric activity in stars older than ∼2 Gyr. In addition, Denissenkov (2010) shows that stellar rotational speeds converge towards similar values after a few billion years, regardless of their initial spin rate; consistently Meibom et al. (2015) confirm a well-defined period-age relation up to ∼2.5 Gyr. Taking all these considerations into account, Bonfanti et al. (2016) explain how we use the gyrochronological relation by Barnes (2010) and the -age relation by Mamajek & Hillenbrand (2008) to compute conservative age lower limits, which they denote by τv and τHK, respectively. Just in assessing age lower limits, the code already takes care of the above-mentioned caveats such that the user is not obliged to be aware of them, and we attain the goal of possible discarding of unlikely pre-MS isochrones before interpolation.

Bonfanti et al. (2016) also note that the parameter degeneracies on the HRD under the turn-off (TO) can be broken by using ρ. The mean stellar density is higher in the MS phase than in the pre-MS, thus it enables us to distinguish MS stars from pre-MS stars and to set a conservative age lower limit called τρ.

The updated implementation of isochrone placement now also accepts [Y/Mg] among the input parameters. Nissen (2016) provide the following [Y/Mg]-age relation (t in Gyr; age scatter 0.6 Gyr):


thereby confirming the trend that was first observed by da Silva et al. (2012) for solar-type stars. As explained by Feltzing et al. (2017) and references therein, yttrium is produced through slow neutron capture process (s-process) above all by asymptotic giant branch (AGB) low mass (1−4 M) stars. This implies that the interstellar medium is gradually enriched by Y and, as a consequence, younger stars have higher Y abundances.

Equation (20) holds for solar analogues (Nissen 2016) and also for giant stars of solar metallicity in the helium-core burning phase (Slumstrup et al. 2017). Feltzing et al. (2017) show that the [Y/Mg]-t trend is also [Fe/H]-dependent and it looses its power of predicting ages for [Fe/H] ≃ −0.5 dex, proving that it is only valid for solar metallicity stars. Therefore, if any [Y/Mg] value is provided on input, first of all our code evaluates whether the stellar metallicity belongs to the −0.2 <  [Fe/H] <  0.2 range. If this does not hold, we simply miss the opportunity of removing possible isochronal degeneracies through the [Y/Mg] index, and if neither v sin i, or ρ suggest discarding pre-MS ages, then the isochrone placement scheme would involve the complete set of age values reported by the grids of isochrones. Otherwise, if [Y/Mg] is available and the star has solar metallicity, Eq. (20) is considered with the aim of investigating whether the star is old enough to avoid loading pre-MS isochrones. In particular, the grids used for the following interpolation contain all the isochrones that are not younger than the τY threshold age, which is set as the minimum value between the age computed through Eq. (20) and 2.5 Gyr; as for the other ages derived from empirical relations, τY represents just a conservative age lower limit.

If more than one empirical age is available, among τv, τHK, τρ, and τY, the maximum age indicator is set as the threshold value, such that only evolutionary grids containing higher ages are loaded and used in all the following computation processes to neglect unlikely young ages. We finally observe that it is possible to compute stellar ages from isochrones only if the observational input parameters vary appreciably with time. This does not hold for very low MS stars, which are characterised by very slow evolution. For instance, according to PARSEC evolutionary tracks, a 0.3 M MS star exhibits a Teff variation ∼0.1% over the whole age of the universe, which is lower than any reasonable input Teff error bar.

2.2.4. Interpolation scheme

Loading of isochrones is based upon stellar metallicity; for convenience we note the adopted relation among Z and [Fe/H] (see Bonfanti et al. 2015), that is


such that [Fe/H] = 0 for Z = Z = 0.01524 consistent with the PARSEC calibration of the Sun (Bressan et al. 2012) taken from Caffau et al. (2011). According to (21), the code initially loads the grid of isochrones corresponding to [Fe/H].

If only Teff and ρ (or log g) are available as inputs, we work in the modified HRD log Teff– log ρ (or log Teff–log g). In all the other combinations of input parameters, it is possible to work in the ordinary HRD, even in the case in which only photometric measurements are available. In fact, if the colour index is among the input parameters, a proper interpolation in the isochrone grid is made to retrieve the corresponding input Teff. Similarly, if both the apparent magnitude and distance are available as inputs, the absolute magnitude is computed and then an interpolation in the grid gives the input stellar luminosity.

If both the gravity proxies log g and ρ are available, a rough value for the input stellar mass may be retrieved. Then a preliminary search for at least an evolutionary state that could be contemporarily compatible with the effective temperature, mass, and the two gravity proxies (within their error bars) is performed. In case this preliminary consistency check fails, the input gravity proxy with the largest error bars is discarded to avoid an insolvable tension between mutually dependent input parameters and evolutionary models, which would produce a not satisfactory interpolation.

The star is placed in the HRD and then we take the following steps.

First, further metallic grids of isochrones are loaded compatibly with the [Fe/H] error bars, thereby giving a huge table of isochrones within which we make interpolations (hereafter this table is called Is).

Second, for each isochrone in each grid, the theoretical evolutionary state that is the closest to the star location in the HRD is considered to then build a new theoretical grid of values (hereafter Is’), which is customised according to the stellar location. Theoretical values of these custom points are inferred through linear interpolation considering the two nearest tabulated points. In this way, all the following computations disregard the specific descrete grid of points that were originally tabulated in theoretical models.

Third, isochrone points entering the Is’ customised grid (say i the row index of the grid) are smoothed through a bi-dimensional Gaussian window function


where y may indicate either log L, log ρ, or log g, depending on the specific HRD we are working in.

Fourth, we define the weight


to be attributed to the ith row of Is’, where y has the same meaning as at the previous point, while may represent, if directly available, the input M and log g (being ). The weight takes into account both the geometrical distance between theoretical and observational parameters and the evolutionary speed vevo of a star in that particular region of the HRD, that is


where (log Teff, 1, y1) and (log Teff, 2, y2) are two consecutive points with ages t1 and t2 along a track displaying the evolutionary path of a theoretical star having the same mass and metallicity of the star to be analysed. The quantity vref is a normalisation factor corresponding to the minimum evolutionary speed that has been detected along the whole track.

The less vevo, the greater the probability of finding a star in that evolutionary stage. According to Eq. (23), the weight of a theoretical point increases as the distance “theory-observations” decreases and vevo decreases.

Finally, the generic X parameter and its uncertainty ΔX to be retrieved according to stellar evolutionary models are then computed through



where is the number of points (rows) of the grid Is’.

At this point the ED phenomenon, also known as microscopic or atomic diffusion (Burgers 1969; Chapman & Cowling 1970; Chaboyer et al. 2001), enters the picture of the algorithm. Element diffusion is caused by the interaction among various chemical species in stars with a convective envelope, which leads to a progressive surface depletion of elements heavier than hydrogen, which sink downwards.

In particular, it turns out that very low mass stars up to 0.5 M exhibit constant Z during their entire life; this is also the case for stars more massive than ∼2 M. In the intermediate mass range, instead, Z diminishes during the MS phase. After that, as the convective envelope deepens, the consequent mixing of elements makes the star approximately as metallic as it was at its birth (see e.g. Bonfanti et al. 2015, Fig. 3).

As already said, the isochrone grids are identified through Zini, while we know only the present day metallicity of a star Z. As shown by Bonfanti et al. (2015, Fig. 4), neglecting the ED effect may result in older age derivations, especially for intermediate age stars. Also Dotter et al. (2017) observed that when initial metallicity is equated to the current metallicity, the resulting isochronal ages of stars can be systematically overestimated by up to ∼20%. Assessing a wrong evolutionary state results in bad mass and radius estimation as well. For this reason, starting from evolutionary tracks (which tabulate Z values along time), we built up a piecewise third degree polynomial interpolation of the curves Zk, l = Zk, l(t), for any given stellar mass and Zini (represented through the subscripts l and k, respectively). In this way, after an iterative process, it is possible to infer what the Zini of a given star was, starting from its present Z and a first raw guess of its age and mass. In fact, at the first iteration, the isochrone grid identified by the present day metallicity of the star is employed to derive a first guess for stellar age (t1, ⋆) and mass (M1, ⋆). After computing Zk, M1, ⋆(t = t1, ⋆) for all the available curves characterised by different Zini (identified by the subscript k), we identify the th curve such that , that is the first guess for the initial metallicity of the star Z1, ini,⋆.

By loading the grid of isochrones labelled Z1, ini,⋆ and repeating all the algorithms described above for retrieving theoretical stellar parameters, new output values for stellar age and mass are obtained. The entire ED procedure is repeated iteratively until a convergence in the Zini, ⋆ values is reached.

According to what we have just described, our isochrone placement evaluates many model stars at each chain step and it emerges with the stellar properties that best match the rest of the parameters. Since an optimisation methodology within each step may lead to underestimated uncertainties, we decided to compare the age uncertainties inferred from the t⋆,isoch PDFs of HD 219134 and WASP-4 (see our analysis and results in Sect. 3) to that stated in the literature. From the review by Soderblom (2010), who evaluates the works by Pont & Eyer (2004), Jørgensen & Lindegren (2005), and Takeda et al. (2007) in assessing isochronal ages, it turns out that typical isochronal age uncertainties are in the range ∼20%–50% (before considering systematics on evolutionary models). Age uncertainties we provide for our test stars are ∼36%, ∼19%, and ∼34%, thus they are consistent with what specialists in the field usually report. We also note that our ages do not purely come from isochrones, and in fact activity index and gyrochronology are used to better constrain them. In this regard, Angus et al. (2019) point out that the combination of isochrone fitting and gyrochronology improves age precision, thus age uncertainties that are lower than those coming from the isochronal analysis alone are expected and our claimed formal uncertainties do not appear underestimated. Therefore, our specific implementation of the optimisation algorithm does not seem to lead to underestimated uncertainties, and we also stress the importance of adding systematic uncertainties (see Sect. 3.3) to our formal uncertainties. However, an optimisation algorithm may be a potential source of underestimated uncertainties and letting the MCMC sample the isochrones directly (like in Eastman et al. 2019) represents a valid alternative to our implementation scheme.

2.3. MCMCI

If in its original implementation the MCMC required a specification of a type of priors or relation involving M and R to fully characterise the exoplanetary system, the current MCMCI5 input form enables us to select the use of theoretical models. Isochrone placement, which is called as a subroutine by the main MCMC, is in charge of deriving M and R (besides the stellar age t) as a consequence of the stellar input parameters to be specified, as explained in Sect. 2.2.2.

The step parameters always adopted on the stellar side are comprised of [Fe/H], Teff (or the colour index), and M, where [Fe/H] and Teff (or the colour index) are also set as priors. Since [Fe/H] is always present among the prior parameters that are randomly perturbed at each step of the chains and several chains of order ∼105 steps are expected to be launched, it would be very time consuming to load the specific metallic grids of isochrones and tracks that are expected to vary from step to step each time. Therefore we decided to load a reasonable amount of evolutionary models since the beginning of the algorithm and to store these in very big matrices. By keeping track of initial and ending indices identifying the various models in the matrices, at each step only the needed sub-matrices are actually used in computations.

Two basic situations are possible in which we specifically draw attention to what is essential to retrieve M and R thanks to evolutionary models.

In the first case scenario the transit is very deep and well defined, such that its impact parameter b can be precisely measured from the LC. We step in [Fe/H], Teff, M, and W. From Wstep, is computed through Eq. (A.4). From Kepler’s third law, is translated into ρ⋆,step, and R⋆,step is derived from M⋆,step and ρ⋆,step. In addition, ρ⋆,step enters isochrone placement (together with [Fe/H]step and Teff, step), to derive the theoretical M⋆,isoch at that step, according to Eq. (25). Finally, the Bayesian penalty


(where ΔM is the uncertainty of M⋆,isoch, that is computed according to Eq. (26)) is added to the merit function given by Eq. (4). In this way, M⋆,step is driven by the theoretical mass value M⋆,isoch, and thus R⋆,step is computed from M⋆,step (constrained by evolutionary models) and ρ⋆,step (constrained by the transit);

In the second case scenario the transit is very shallow, such that b is not directly inferrable from the LC. Besides [Fe/H], Teff and M, R is a also jump parameter in this case and may be also set as prior. At each step, input parameters for isochrone placement are [Fe/H]step, Teff, step, and at least one indirect piece of information involving the stellar radius (e.g. the stellar luminosity through apparent magnitude and parallax, the spectroscopic log g, or an interferometric measurement of the radius). Then, according to Eqs. (25) and (26), isochrone placement gives M⋆,isoch ± ΔM and R⋆,isoch ± ΔR as output values, which enter the following Bayesian penalties:


to be added to the merit function expressed by Eq. (4). The quantity ρ⋆,step is computed consistently with M⋆,step and R⋆,step, whose values are driven by the theoretical M⋆,isoch and R⋆,isoch coming out at each step. Thanks to ρ⋆,step, is computed through Kepler’s third law and finally Wstep is expressed as a function of by inverting Eq. (A.4), such that transit parameters are constrained a posteriori by evolutionary models. We note that if the stellar radius has been also set as a prior as a result of the BP given by Eq. (3), its value also feeds back into ρ and into the corresponding transit parameters in general.

From M⋆,step and R⋆,step, all the relevant planetary parameters can be computed, especially Rp, step (from dFstep, if any LC is available) and Mp, step (if any RV time series or, at least, K is available). At the end of all the chains, PDFs of the parameters of the exoplanetary system are built upon the step values.

After a first MCMCI run when CFs referring to each LC are computed through Eq. (19), the photometric error bars are properly rescaled to obtain reliable error bars on the fitted parameters. If everything else is the same, the MCMCI is launched a second time to retrieve reliable PDFs of both stellar and planetary parameters.

Algorithm tests and all the subsequent analysis that is presented in Sect. 3 have been run on a DELL Latitude laptop (1 TB SSD; 16 GB RAM; Intel core i7-7820 HQ @ 2.90 GHz; FSB = 3.68 GHz). One single run for analysing HD 219134 system (4 LC time series for a total amount of 3322 data points, 1 RV time series made of 663 data points, four orbiting planets, element diffusion switched-on) took ∼30 h. One single run for analysing Wasp-4 system (14 LC time series for a total amount of 4010 data points, one orbiting planet, element diffusion switched-on) took ∼20 h. If we considered a very basic fit of one single planet orbiting a star with just 1 available LC, the expected run-time is seven to eight hours if ED is switched on.

By switching off ED in performing the isochrone placement technique, we save a factor of 2 in run-time. However, as we already stressed, neglecting ED may imply an age overestimation up to ∼20%. Instead, the impact on the isochronal M and R can be estimated up to a few percent.

We are aware of other implementations of the MCMC algorithm, such as the differential evolution (DE; Ter Braak 2006) or the affine invariant (AI; Foreman-Mackey et al. 2013) technique. We are currently working on the implementation of the AI-MCMC procedure, which will also contain Gaussian processes (Rasmussen & Williams 2005) for modelling correlated noise, and according to our first estimations, this technique will increase the run-time speed by a factor 10. The full FORTRAN implementation of the AI-MCMC requires likely one year, since it is necessary to revise the architecture of our current MCMC code. Thus, we prefer to share our current MCMCI version now and then considering the ongoing update as part of a future work.

3. Analysis and results

To test our algorithm we selected two planet-hosting stars that are representative of two different test cases. One of these stars is HD 219134, hosting several low mass planets, which has already been studied by, for example Motalebi et al. (2015), Vogt et al. (2015), and Gillon et al. (2017). The other is WASP-4, which is known to host a heavily bloated hot Jupiter owing to transit analyses already carried out by Wilson et al. (2008), Gillon et al. (2009b), Winn et al. (2009), Southworth et al. (2009), Sanchis-Ojeda et al. (2011), and Nikolov et al. (2012).

3.1. HD 219134

HD 219134 (also known as GJ 892 or HIP 114622) is a bright close-by K3 V star. Its relevant input parameters used during our analyses are listed in Table 1. Our first attempt to analyse this system took the tabulated Teff among input parameters, but no convergence was reached in theoretical models.

Table 1.

HD 219134: input stellar parameters considered in this work.

Actually, as shown in Fig. 1 (right panel), the Teff value provided by Boyajian et al. (2012), which has been also considered, for example by Motalebi et al. (2015) and Gillon et al. (2017), is completely inconsistent with theoretical models, regardless of the metallicity. This value was empirically computed by Boyajian et al. (2012), considering the interferometric measurement of the stellar diameter (which yields R = 0.778 R) and an empirical estimate of the bolometric flux by fitting photometry to spectral templates. On the other hand we observed that the B − V colour index (which is more straightforward to obtain) has a reasonable value, if compared with stellar models (Fig. 1, left panel). In addition, interpolation inside theoretical models suggests corresponding Teff = 4836 K or 4902 K (depending on the L- or R-approach; see later), which agree with other estimates provided in the literature as shown in Fig. 2. For these reasons, we decided to start from B − V, rather than Teff.

thumbnail Fig. 1.

HD 219134 location on the colour-magnitude diagram (left panel) and on the log TeffMV diagram (right panel) together with some reference isochrones. Red stands for 1 Gyr models, while blue stands for 10 Gyr models; solid lines are representative of the nominal stellar metallicity (Z = 0.020 → [Fe/H] = 0.11); dashed lines refer to extremely super-solar metallicity Z = 0.070 → [Fe/H] = 0.66.

thumbnail Fig. 2.

Values of Teff found in the literature for HD 219134. Values on the x-axis correspond to the following references: (1) Valenti & Fischer (2005); (2) Mishenina et al. (2008); (3) Soubiran et al. (2008); (4) Prugniel et al. (2011); (5) Ramírez et al. (2013); (6) Motalebi et al. (2015), EW approach; (7) Motalebi et al. (2015), SPC approach; (8) Our work, L approach; and (9) our work, R approach.

The transit analysis was carried out because of the four LCs observed through the Spitzer/IRAC detector (Fazio et al. 2004) and the RV time series gathered by the HARPS-N spectrograph (Cosentino et al. 2012) already considered in Gillon et al. (2017). Both the nominal model made of one star and four planets and the choice of the baseline functions were taken from Gillon et al. (2017). We recall that two LCs contain the transit of planet b, while the other two the transit of planet c. Planets d and f are not transiting and their presence was deduced from the RV analysis.

We notice that there are claims of further orbiting planets in the literature. Motalebi et al. (2015) talk about four planets containing the transit of b, c, d, plus one long period planet that orbits the star with P = 1842 days, but without considering planet f of our analysis. Instead, Vogt et al. (2015) claim the presence of six planets; besides b, c, d, and f which we considered, they also find planet g with P = 94.2 days and planet h with period P = 2247 days. The high polynomial degree (order 4) that enters the temporal de-trending of the RV time series suggests the likely presence of other (long period) planets besides the four we selected. On the other hand, there is not full agreement in the literature about the total number and period of the further planets, and deeply investigating and characterising the presence of other planets is beyond the scope of this paper; instead, we want to show the behaviour of our MCMCI in front of a challenging multi-planetary system as HD 219134.

For all four planets, T0, P, and K2 were assumed as jump parameters. Also the eccentricity e has been let free to vary, except for planet b, for which a fixed e = 0 value was set, according to Gillon et al. (2017), who did not register any improvement in the BIC by let e vary. Consistently, Gillon et al. (2017) state that the closiness of HD 219134 b to its star implies a circularisation timescale of 80 Myr even assuming as tidal quality factor the maximum value that is derived for terrestrial planets and satellites of the Solar System (Murray & Dermott 1999). The same computation for planet c yields a timescale of 2.5 Gyr, therefore we preferred letting e free to jump for this and the other planets. In case of the two transiting planets b and c as well as dF and b′ were set as jump parameters. The transits result to be too shallow to reliably constrain W. We preferred not to set this a priori, but rather to count on observational stellar parameters and on theoretical evolutionary models to characterise the star completely so that W was later inferred.

On the stellar side, two different sets of input parameters were considered. The first set considers [Fe/H], B − V, Prot, , V magnitude and parallactic distance d. Interpolation within isochrones enabled us to recover input stellar luminosity L from V and d; Bayesian penalties (BPs) involving [Fe/H], B − V and L were added to the merit function by setting the corresponding parameters as priors. Hereafter the analysis starting from this input set is called Luminosity prior approach (L approach).

The second set considers [Fe/H], B − V, Prot, , input stellar radius Rprior from interferometry. In this case BPs entering the merit function involved [Fe/H], B − V and Rprior. Hereafter the analysis starting from this input set is called Radius prior approach (R approach).

The MCMCI was launched considering five chains of 100 000 steps each, for both the L and R approach. In both cases, a qd model for LD was assumed and c1 and c2 were properly computed to build their corresponding BPs to be added to the merit function. Furthermore, also the sum BPM + BPM given by (28) was added to the merit function.

After a first run, the MCMCI was launched again by applying the CFs recovered at the first run through Eq. (19), to retrieve all the relevant stellar and planetary parameters (convergence checked thanks to the Gelman-Rubin test), that are summarized in Tables 2 and 3. In particular, the transiting HD 219134 b and c are confirmed to be two super-Earth rocky planets. Their corresponding LCs with the super-imposed transit models as inferred from the L-approach are displayed in Fig. 3.

thumbnail Fig. 3.

Superimposed transiting models on LCs, as derived from the L approach; the R approach gives analogous results. Two top panels: transit of HD 219134 b, two bottom panels: HD 219134 c.

Table 2.

HD 219134: output stellar parameters derived considering L approach or R approach.

Table 3.

Planets hosted by HD 219134.

A comparison involving radii and masses of HD 219134 and its planets as derived by various authors is presented in Table 4. Gillon et al. (2017) analyse the system assuming stellar Rprior as prior similarly to our R approach and find that very good agreement holds between these two estimations. From Table 4 we also notice that our uncertainties affecting Mp and Rp are generally lower than those reported in the literature, especially in the case of planetary masses. This suggests that the integration of stellar theoretical models with the MCMC algorithm is able to reduce the PDF width of the output parameters.

Table 4.

Comparison involving masses and radii derived by different authors for HD 219134 system.

Finally, we comment on the results we obtained with our two different and independent L and R approaches with a particular reference to our derived R values. Input V magnitude and Gaia parallax are totally independent of (and possibly inconsistent with) the interferometric measure of the radius. Actually, the comparison between the L and R approach aims to investigate the mutual consistency of input parameters by looking at the output results.

In Fig. 4 the red histogram is the output R PDF, which has been inferred considering ([Fe/H]step, B − Vstep, L⋆, step) as inputs for isochrone fitting and assuming a prior on the parallax-based stellar luminosity Lprior. Instead, the blue histogram is the output R PDF that has been inferred considering ([Fe/H]step, B − Vstep, R⋆, step) as inputs for isochrone fitting, and assuming a prior (black Gaussian) on the interferometric radius Rprior.

thumbnail Fig. 4.

Red histogram: posterior R PDF, as derived from the L approach. Blue histrogram: posterior R PDF, as from the R approach. Black Gaussian: interferometric prior on stellar radius that has been imposed in the R approach. The black Gaussian prior essentially overrides the isochronal constraint (which is looser) and drives the output R value. See text for a complete discussion.

thumbnail Fig. 5.

WASP-4: regression of M vs. (left panel), and of Mp vs. (right panel). A strong correlation is expected because ρ is constrained from the transit, and (Mp, Rp) are related to (M, R), transit parameters and the RV semi-amplitude.

It is important to stress that the red histogram is the isochronal-consistent R PDF when L⋆,step enters the fitting algorithm. On the other hand, if R⋆,step is assumed to be an input for isochrone placement instead of L⋆,step, the isochronal-consistent R PDF that is produced in output is expected to be consistent with R⋆,step, and thus possibly shifted with respect to the red histogram if the input radius and luminosity are not consistent for the given Teff and [Fe/H]. Actually, if R⋆step is assumed on input, the isochrone-based Bayesian penalty BPR alone (Eq. (28)) would drive to an output radius Rout, isoch ≈ 0.767 ± 0.020 R: if compared with the L approach, the amplified uncertainty suggests that isochrones better agree with the luminosity input values, rather than the radius input values. If we now consider that the R approach is also characterised by the presence of the additional black Gaussian prior whose σ = 0.005 R, it turns out that its correspondent BP (Eq. (3)) plays the most prominent role in determining the output R PDF. The blue histogram enters this scenario consistently, as it is just slightly shifted towards lower values with respect to the black Gaussian.

The inconsistence of the derived R values (red versus blue histogram in Fig. 4) is the consequence of the inconsistence between the parallax-based stellar luminosity and the interferometric radius, for the given photometry and metallicity. Reasons for this tension may be found in the input photometry and [Fe/H], surface gravity that is implied by stellar position in the HRD (which all influence Teff and bolometric correction, upon which L is inferred), intrinsic limits of evolutionary models, and the interferometric technique for retrieving stellar radius, which may suffer some systematics (see e.g. White et al. 2018).

3.2. WASP-4

WASP-4 is a V = 12.5 mag G7V star orbited by a transiting hot Jupiter. We carried out the transit analysis with two LCs by Winn et al. (2009) (Sloan z′ band exposures taken at Las Campanas Observatory in Chile using the MagIC camera) and 12LCs derived from three transits simultaneously observed in the four Sloan g′, r′, i′, z′ bands by Nikolov et al. (2012) at La Silla Observatory. We are aware of the availability of other LCs since other studies of WASP-4 have been already done, as specified at the beginning of Sect. 3. Carrying out an analysis considering all the available material in the literature is beyond the scope of the paper, which instead tests the MCMCI and its ability to contemporaneously deal with several LCs taken in different bands.

By launching multiple small runs of the MCMCI (chains of 10 000 steps each) and monitoring the BIC value after varying the baseline, we realised that the photometric time series do not need a baseline function more complex than a simple scalar (normalisation). We considered dF, b′, W, T0, and P as jump parameters and, in addition, we set Teff and [Fe/H] as priors. Thus BPs on Teff and [Fe/H] were considered, besides those on the c1 and c2 LD coefficients (qd model assumed). Moreover, BPM from Eq. (27) was also added to the merit function. The orbit was assumed to be circular, following the indications by Sanchis-Ojeda et al. (2011); these authors report that no eccentricity has been revealed in the study of several RV time series (Wilson et al. 2008; Madhusudhan & Winn 2009; Pont et al. 2011) and occulation data (Beerer et al. 2011). We complemented this information by specifying also v sin i to improve the convergence within theoretical models and the RV semi-amplitude K to get an output value for the planetary mass.

Bouma et al. (2019) has stressed that the sequence of transit recently observed by TESS happened earlier than expected. The authors speculate that this timing variation may be caused by tidal orbital decay or apsidal precession, while it seems unlikely to be the presence of a third body that could perturb the star-planet system. Anyway, we decided to launch a run in which transit timing variations (TTVs; Holman & Murray 2005) were also allowed in our model, but the BIC we obtained was higher with respect to the case in which TTVs were not considered. Therefore we avoided including TTVs in our analysis.

As in the case of HD 219134, first the MCMCI code was launched to estimate the CFs (Eq. (19)) and then it was launched a second time after applying the CFs. Five chains of 100 000 steps each were considered at each run and the convergence of the derived parameters was checked using the Gelman-Rubin test. As a reference, the two LCs by Winn et al. (2009) and the three LCs observed in the g′ band by Nikolov et al. (2012) are shown in Fig. 6, together with the superimposed transit models.

thumbnail Fig. 6.

Light curves of WASP-4 b. First two panels show the data provided by Winn et al. (2009), who observed the transit in the Sloan z′ band. The other three panels show the observation by Nikolov et al. (2012) in the Sloan g′ band. The red thick line represents the transit model as inferred from our analysis.

Input and derived parameters of the WASP-4 system are listed in Table 5. We note that the relative uncertainty on ρ is lower than those affecting M and R. In fact, unlike the case of HD 219134, in this work ρ is constrained from the transit6 and then used to derive (M, R), which depend on ρ itself and on the evolutionary stellar models. Various considerations hold for log g and ρp: they are computed from the pairs (M, R) and (Mp, Rp), respectively, thus according to error propagation, we could expect their relative uncertainties to be higher than those affecting M, R, Mp, and Rp. But this is not the case since we have a derived output uncertainty versus that should be expected from error propagation; similarly, we have versus .

Table 5.

WASP-4 planetary system.

However, uncertainties in agreement with standard error propagation are expected only if the pairs (M, R) and (Mp, Rp) are uncorrelated. Actually, a correlation between M and R exists because, besides isochrones, we have a strong constraint on thanks to the transit. In addition, Mp (resp. Rp) differ from M (resp. R) by factors that are constrained from the transit and the RV semi-amplitude; thus a correlation between Mp and Rp is also expected, although it is a bit more diluted. What has just been described is shown in Fig. 5, where the M (resp. Mp) values coming from all the chain steps are plotted versus the corresponding (resp. ) values. Linear correlations are clear and the coefficients of determination are (left panel) and (right panel).

If a strict linear correlation held (r2 = 1), no dispersion in and in would be registered, thus both their uncertainties would be zero. The equation represents the fraction of root mean square that is not explained by the linear regression, thus it is statistically expected that error bar extensions coming from simple error propagation (where the involved quantities are considered independent) may be reduced up to a factor σne or so. Statistically, uncertainties may go down to the following estimated levels:



It turned out that the statistically expected uncertainties (indicated by the subscript “exp”) are in agreement with the actual derived ones (indicated by the subscript “out”).

Table 6 reports a comparison involving stellar and planetary masses and radii, as derived by different authors. Our results have the advantage of retrieving M and R owing to the strict interaction of the MCMC algorithm with the stellar evolutionary models, which yields to very precise Mp and Rp values, if compared to what is stated in the literature.

Table 6.

WASP-4: masses and radii of both the star and the planet derived by different authors.

3.3. Comments on the uncertainties of isochronal parameters

Output parameters that are computed from isochrones (i.e. t, M and R) are given as if stellar evolutionary models are perfect. Therefore their respective errors are likely to be underestimated and should be considered as internal. One way to attribute more realistic uncertainties to isochronal parameters is by comparing our results with another set of evolutionary models. If the target stars and the interpolating algorithm are the same, differences arising in the output are to be attributed to the different input physics of the evolutionary models, such as the equations of state, adopted solar mixture, initial chemical composition of stars, nuclear reaction rates, opacities, overshooting treatment, mixing-length parameter, and atmospheric models. Setting up all these ingredients depends on our knowledge of stellar evolution, different choices produce different models, and consequent variations in the output parameter estimations are a measure of the uncertainty of evolutionary models.

Our isochrone placement as autonomous routine is also able to interact with the grids of tracks and isochrones produced by CLES (Code Liégeois d’Évolution Stellaire; Scuflaire et al. 2008). These models span just a limited range in mass (0.90 ≤ M ≤ 1.30M) and metallicity (0.008 ≤ Z ≤ 0.018) so far: while HD 219134 is outside the mass range, luckily WASP-4 falls inside. Among the several differences in terms of stellar input physics between PARSEC and CLES models, one is related to the helium content Y. In PARSEC models, Y is assumed to increase with Z, according to


where Yp = 0.2485 (Komatsu et al. 2011) is the primordial helium abundance, while the helium-to-metal enrichment ratio is derived from solar calibration. Instead, is not fixed in CLES, such that Y is not a function of Z, and actually for any given value of Z, four possible values of Y are available, from Y = 0.25 to Y = 0.28 at steps of 0.01.

To evaluate the impact of stellar model systematics, we decided to analyse WASP-4 through our autonomous isochrone placement routine. Input parameters were Teff, [Fe/H], ρ as derived from transit, and v sin i. Interpolations were performed considering both PARSEC models and the four sets of CLES grids (one per each Y value). The helium content cannot be generally constrained unless high quality asteroseismic oscillation frequencies are available (see e.g. Lebreton & Goupil 2014; Buldgen et al. 2016).

According to Eq. (31), after considering element diffusion, PARSEC models assign an initial helium content Y = 0.28 to WASP-4. By comparing the output t, M, and R values derived by the two models at the same Y = 0.28 content, the relative differences (CLES versus PARSEC) in age, mass, and radius are −19%, +1.3% and +0.5%, respectively. If we also consider the effect of helium content on theoretical models, output variations on age, mass, and radius rise to −23%, +5.8% and +1.9%, respectively.

From this test case, we deduce that for having realistic error bars on the isochronal output parameters, conservative reference systematics to be added to the internal uncertainties are of the order of ∼20%, ∼6%, and ∼2% for age, mass, and radius, respectively. A complete study about isochrone precision is beyond the scope of the paper; properly comparing two different evolutionary models would require a wide stellar sample, possibly spanning different locations in the HRD, while the just mentioned percentage values are provided to give an initial idea about isochrone precision.

4. Conclusions

We presented our MCMCI code that was born by merging the MCMC code (see e.g. Gillon et al. 2010, 2012) and the isochrone placement algorithm developed by Bonfanti et al. (2015, 2016). As a reference, we took the opportunity to recall the working details of the two codes separately considering their most recent implementations.

Summing up, as the name suggests, the MCMC code is based upon Markov chain Monte Carlo statistics and it aims to retrieve the main parameters of an exoplanetary system once that photometric or RV time series are provided and stellar observational parameters are available. Given the indirect nature of the transit and RV methods, M and R have to be known to derive the mass and radius of the planet, which is possible as a consequence of stellar evolutionary models. Starting from this necessity, we integrated isochrone placement into the MCMC, thereby building a new MCMCI code, which at each step of the chain “asks” theoretical models to provide M and R, given the available observational data.

The main advantage is dealing with a powerful tool that retrieves PDFs of all the planetary and stellar parameters at a time without the need to split the analysis for separately recovering theoretical stellar parameters. In particular, the derived PDFs of M and R are a strict consequence of the values that are obtained at each step of the chains from the perturbation of the observational parameters. In addition, our code is also able to establish the stellar age t, considering the constraints coming both from stellar evolutionary models and several empirical age relations.

The MCMCI was tested by selecting two already studied planetary systems that may be considered representative of two different reference cases.

First we considered HD 219134, which hosts several planets (four according to the model we followed), two of which are transiting super Earths. This is the typical situation of shallow transits for which the transit impact parameter b cannot be directly inferred from the LCs and knowledge of R is fundamental. We followed two different approaches, both considering an interferometric measurement of R and also recovering it once Teff was calibrated from B − V and L was deduced from V magnitude and distance.

Second we considered WASP-4, that is known to host a hot Jupiter. The transit features are evident such that ρ can be inferred from the LCs and it enters the set of input parameters for inferring M and R from theoretical models.

Very good agreement was found by comparing our results to those reported in the literature, which suggests that our tool represents a powerful and integrated solution to analyse an exoplanetary system fully giving refined and reliable stellar and planetary parameters at the same time. In particular, planetary mass values are more precise if compared with previous determinations available in the literature.

Aside from the specific exoplanetary context, we also took the opportunity to speculate about isochrone precision in this specific work. By employing PARSEC and CLES stellar evolutionary models for studying WASP-4, we quantify the systematics affecting isochronal age (δt ∼ 20%), mass (δM ∼ 6%), and radius (δR ∼ 2%). Although we are aware that these percentage values should be considered as just indicative estimates of the isochrone precision because they are inferred from a single test case, these values draw attention to the fact that any isochronal parameter suffers intrinsic uncertainties, which are related to our current ability to model stars.


It is just a reference value that considers the effective length of a chain, where 20% burn-in length is subtracted to an original chain made of 100 000 steps. Actually Eq. (17) is almost Lc independent for Lc of order of thousands, which is a common practice in MCMC applications.


The MCMCI code is available at


Constraining stellar density from transit with a relative uncertainty lower than 2% as in this case is not so common, and it has been possible because of the prominent transit depth and high quality LCs. According to NASA Exoplanet Archive, the median uncertainty for transit-inferred ρ is δρ⋆, t ∼ 9%, and only ∼5% of the systems have δρ⋆, t <  2%.


This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. M. Gillon is Senior Research Associate at F.R.S-FNRS. The research leading to these results has received funding from the ARC grant for Concerted Research Actions, financed by the Wallonia-Brussels Federation. A.B. is a postodoctoral researcher funded by the ESA PRODEX programme for CHEOPS (PEA 4000113509). The authors thank the referee, Jason Eastman; all his suggestions and comments lead to a remarkable improvement of the paper. Moreover, the authors thank A. Noëls, V. Van Grootel and A. Collier-Cameron for fruitful discussions that help improving the paper as well.


  1. Anderson, D. R., Collier Cameron, A., Hellier, C., et al. 2011, ApJ, 726, L19 [NASA ADS] [CrossRef] [Google Scholar]
  2. Angus, R., Morton, T. D., Foreman-Mackey, D., et al. 2019, AJ, 158, 173 [NASA ADS] [CrossRef] [Google Scholar]
  3. Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Baluev, R. V. 2018, Astron. Comput., 25, 221 [NASA ADS] [CrossRef] [Google Scholar]
  5. Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [Google Scholar]
  6. Barnes, S. A. 2010, ApJ, 722, 222 [NASA ADS] [CrossRef] [Google Scholar]
  7. Barragán, O., Gandolfi, D., & Antoniciello, G. 2019, MNRAS, 482, 1017 [NASA ADS] [CrossRef] [Google Scholar]
  8. Basri, G., Borucki, W. J., & Koch, D. 2005, New Astron. Rev., 49, 478 [NASA ADS] [CrossRef] [Google Scholar]
  9. Beerer, I. M., Knutson, H. A., Burrows, A., et al. 2011, ApJ, 727, 23 [NASA ADS] [CrossRef] [Google Scholar]
  10. Betancourt, M. 2017, ArXiv e-prints [arXiv:1701.02434] [Google Scholar]
  11. Bonfanti, A., Ortolani, S., Piotto, G., & Nascimbeni, V. 2015, A&A, 575, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bonfanti, A., Ortolani, S., & Nascimbeni, V. 2016, A&A, 585, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  14. Bouma, L. G., Winn, J. N., Baxter, C., et al. 2019, AJ, 157, 217 [NASA ADS] [CrossRef] [Google Scholar]
  15. Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, ApJ, 757, 112 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 [NASA ADS] [CrossRef] [Google Scholar]
  17. Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, Eur. Phys. J. Web Conf., 47, 03005 [CrossRef] [EDP Sciences] [Google Scholar]
  18. Buldgen, G., Reese, D. R., & Dupret, M. A. 2016, A&A, 585, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Burgers, J. 1969, Flow Equations for Composites Gases (New York: Academic Press) [Google Scholar]
  20. Caffau, E., Ludwig, H.-G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255 [NASA ADS] [CrossRef] [Google Scholar]
  21. Casella, G., & George, E. I. 1992, Am. Stat., 46, 167 [Google Scholar]
  22. Chaboyer, B., Fenton, W. H., Nelan, J. E., Patnaude, D. J., & Simon, F. E. 2001, ApJ, 562, 521 [NASA ADS] [CrossRef] [Google Scholar]
  23. Chaplin, W. J., Basu, S., Huber, D., et al. 2014, ApJS, 210, 1 [NASA ADS] [CrossRef] [Google Scholar]
  24. Chapman, S., & Cowling, T. G. 1970, The Mathematical Theory of Non-uniform Gases. An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, 3rd edn. (Cambridge: Cambridge University Press) [Google Scholar]
  25. Charbonneau, D., Knutson, H. A., Barman, T., et al. 2008, ApJ, 686, 1341 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  26. Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525 [NASA ADS] [CrossRef] [Google Scholar]
  27. Claret, A., & Bloemen, S. 2011, A&A, 529, A75 [Google Scholar]
  28. Collier Cameron, A., Bruce, V. A., Miller, G. R. M., Triaud, A. H. M. J., & Queloz, D. 2010, MNRAS, 403, 151 [NASA ADS] [CrossRef] [Google Scholar]
  29. Cosentino, R., Lovis, C., Pepe, F., et al. 2012, in Ground-based and Airborne Instrumentation for Astronomy IV, Proc. SPIE, 8446, 84461V [Google Scholar]
  30. da Silva, R., Porto de Mello, G. F., Milone, A. C., et al. 2012, A&A, 542, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Davenport, J. R. A., Hawley, S. L., Hebb, L., et al. 2014, ApJ, 797, 122 [NASA ADS] [CrossRef] [Google Scholar]
  32. Dawson, R. I., & Johnson, J. A. 2012, ApJ, 756, 122 [NASA ADS] [CrossRef] [Google Scholar]
  33. Dawson, R. I., & Johnson, J. A. 2018, ARA&A, 56, 175 [NASA ADS] [CrossRef] [Google Scholar]
  34. Denissenkov, P. A. 2010, ApJ, 719, 28 [NASA ADS] [CrossRef] [Google Scholar]
  35. Dotter, A. 2016, ApJS, 222, 8 [NASA ADS] [CrossRef] [Google Scholar]
  36. Dotter, A., Conroy, C., Cargile, P., & Asplund, M. 2017, ApJ, 840, 99 [NASA ADS] [CrossRef] [Google Scholar]
  37. Eastman, J., Siverd, R., & Scott Gaudi, B. 2010, PASP, 122, 935 [NASA ADS] [CrossRef] [Google Scholar]
  38. Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83 [NASA ADS] [CrossRef] [Google Scholar]
  39. Eastman, J. D., Rodriguez, J. E., Agol, E., et al. 2019, PASP, submitted [arXiv:1907.09480] [Google Scholar]
  40. Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 90, 2262 [Google Scholar]
  41. Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10 [NASA ADS] [CrossRef] [Google Scholar]
  42. Feltzing, S., Howes, L. M., McMillan, P. J., & Stonkutė, E. 2017, MNRAS, 465, L109 [NASA ADS] [CrossRef] [Google Scholar]
  43. Ford, E. B. 2005, AJ, 129, 1706 [Google Scholar]
  44. Ford, E. B. 2006, ApJ, 642, 505 [NASA ADS] [CrossRef] [Google Scholar]
  45. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
  46. Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220 [NASA ADS] [CrossRef] [Google Scholar]
  47. Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [Google Scholar]
  48. Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457 [Google Scholar]
  49. Gillon, M., Demory, B.-O., Triaud, A. H. M. J., et al. 2009a, A&A, 506, 359 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Gillon, M., Smalley, B., Hebb, L., et al. 2009b, A&A, 496, 259 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Gillon, M., Lanotte, A. A., Barman, T., et al. 2010, A&A, 511, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Gillon, M., Doyle, A. P., Lendl, M., et al. 2011, A&A, 533, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Gillon, M., Triaud, A. H. M. J., Fortney, J. J., et al. 2012, A&A, 542, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Gillon, M., Demory, B.-O., Van Grootel, V., et al. 2017, Nat. Astron., 1, 0056 [NASA ADS] [CrossRef] [Google Scholar]
  55. Giménez, A. 2006, ApJ, 650, 408 [NASA ADS] [CrossRef] [Google Scholar]
  56. Günther, M. N., & Daylan, T. 2019, Astrophysics Source Code Library [record ascl:1903.003] [Google Scholar]
  57. Hartman, J. D., Bakos, G. Á., Bayliss, D., et al. 2019, AJ, 157, 55 [NASA ADS] [CrossRef] [Google Scholar]
  58. Hastings, W. K. 1970, Biometrika, 57, 97 [Google Scholar]
  59. Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  60. Holman, M. J., Winn, J. N., Latham, D. W., et al. 2006, ApJ, 652, 1715 [NASA ADS] [CrossRef] [Google Scholar]
  61. Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398 [NASA ADS] [CrossRef] [Google Scholar]
  62. Iglesias-Marzoa, R., López-Morales, M., & Morales, M. J. A. 2015, PASP, 127, 567 [NASA ADS] [CrossRef] [Google Scholar]
  63. Jørgensen, B. R., & Lindegren, L. 2005, A&A, 436, 127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., & Megeath, S. T. 2008, ApJ, 673, 526 [NASA ADS] [CrossRef] [Google Scholar]
  65. Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [NASA ADS] [CrossRef] [Google Scholar]
  66. Kopal, Z. 1959, Close Binary Systems (London: Chapman & Hall) [Google Scholar]
  67. Kreidberg, L. 2015, PASP, 127, 1161 [NASA ADS] [CrossRef] [Google Scholar]
  68. Lebreton, Y., & Goupil, M. J. 2014, A&A, 569, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Loeb, A., & Gaudi, B. S. 2003, ApJ, 588, L117 [NASA ADS] [CrossRef] [Google Scholar]
  70. Madhusudhan, N., & Winn, J. N. 2009, ApJ, 693, 784 [NASA ADS] [CrossRef] [Google Scholar]
  71. Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [NASA ADS] [CrossRef] [Google Scholar]
  72. Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [NASA ADS] [CrossRef] [Google Scholar]
  73. Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488 [NASA ADS] [CrossRef] [Google Scholar]
  74. Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77 [NASA ADS] [CrossRef] [Google Scholar]
  75. Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [Google Scholar]
  76. Meibom, S., Barnes, S. A., Platais, I., et al. 2015, Nature, 517, 589 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  77. Mishenina, T. V., Soubiran, C., Bienaymé, O., et al. 2008, A&A, 489, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Morris, S. L. 1985, ApJ, 295, 143 [NASA ADS] [CrossRef] [Google Scholar]
  79. Motalebi, F., Udry, S., Gillon, M., et al. 2015, A&A, 584, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Murray, C., & Correia, A. 2011, in Exoplanets, ed. S. Seager (Tucson, AZ: University of Arizona Press), 526, 15 [NASA ADS] [Google Scholar]
  81. Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
  82. Neal, R. M. 2012, Handbook of Markov Chain Monte Carlo [arXiv:1206.1901] [Google Scholar]
  83. Nikolov, N., Henning, T., Koppenhoefer, J., et al. 2012, A&A, 539, A159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Nissen, P. E. 2016, A&A, 593, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Pont, F., & Eyer, L. 2004, MNRAS, 351, 487 [NASA ADS] [CrossRef] [Google Scholar]
  86. Pont, F., Zucker, S., & Queloz, D. 2006, MNRAS, 373, 231 [NASA ADS] [CrossRef] [Google Scholar]
  87. Pont, F., Husnoo, N., Mazeh, T., & Fabrycky, D. 2011, MNRAS, 414, 1278 [NASA ADS] [CrossRef] [Google Scholar]
  88. Prugniel, P., Vauglin, I., & Koleva, M. 2011, A&A, 531, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Ramírez, I., Allende Prieto, C., & Lambert, D. L. 2013, ApJ, 764, 78 [NASA ADS] [CrossRef] [Google Scholar]
  90. Rasmussen, C. E., & Williams, C. K. I. 2005, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) (The MIT Press) [Google Scholar]
  91. Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [NASA ADS] [CrossRef] [Google Scholar]
  92. Roberts, G. O., Gelman, A., & Gilks, W. R. 1997, Ann. Appl. Probab., 7, 110 [CrossRef] [Google Scholar]
  93. Rodríguez, A., & Ferraz-Mello, S. 2010, in EAS Pub. Ser., eds. K. Gożdziewski, A. Niedzielski, & J. Schneider, 42, 411 [CrossRef] [Google Scholar]
  94. Rodríguez, A., Ferraz-Mello, S., Michtchenko, T. A., Beaugé, C., & Miloni, O. 2011, MNRAS, 415, 2349 [NASA ADS] [CrossRef] [Google Scholar]
  95. Rosenfield, P., Marigo, P., Girardi, L., et al. 2016, ApJ, 822, 73 [NASA ADS] [CrossRef] [Google Scholar]
  96. Rouan, D., Baglin, A., Barge, P., et al. 1999, Phys. Chem. Earth C, 24, 567 [Google Scholar]
  97. Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley) [Google Scholar]
  98. Sanchis-Ojeda, R., Winn, J. N., Holman, M. J., et al. 2011, ApJ, 733, 127 [NASA ADS] [CrossRef] [Google Scholar]
  99. Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
  100. Scuflaire, R., Théado, S., Montalbán, J., et al. 2008, Ap&SS, 316, 83 [NASA ADS] [CrossRef] [Google Scholar]
  101. Seager, S., & Mallén-Ornelas, G. 2003, ApJ, 585, 1038 [NASA ADS] [CrossRef] [Google Scholar]
  102. Sharma, S., Stello, D., Buder, S., et al. 2018, MNRAS, 473, 2004 [NASA ADS] [CrossRef] [Google Scholar]
  103. Siverd, R. J., Beatty, T. G., Pepper, J., et al. 2012, ApJ, 761, 123 [NASA ADS] [CrossRef] [Google Scholar]
  104. Skumanich, A. 1972, ApJ, 171, 565 [NASA ADS] [CrossRef] [Google Scholar]
  105. Slumstrup, D., Grundahl, F., Brogaard, K., et al. 2017, A&A, 604, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Smith, A. M. S., Gandolfi, D., Barragán, O., et al. 2017, MNRAS, 464, 2708 [NASA ADS] [CrossRef] [Google Scholar]
  107. Soderblom, D. R. 2010, ARA&A, 48, 581 [NASA ADS] [CrossRef] [Google Scholar]
  108. Soubiran, C., Bienayme, O., Mishenina, T. V., & Kovtyukh, V. V. 2008, VizieR Online Data Catalog: J/A+A/480/91 [Google Scholar]
  109. Southworth, J., Hinse, T. C., Burgdorf, M. J., et al. 2009, MNRAS, 399, 287 [NASA ADS] [CrossRef] [Google Scholar]
  110. Sozzetti, A., Torres, G., Charbonneau, D., et al. 2007, ApJ, 664, 1190 [NASA ADS] [CrossRef] [Google Scholar]
  111. Takeda, G., Ford, E. B., Sills, A., et al. 2007, ApJS, 168, 297 [NASA ADS] [CrossRef] [Google Scholar]
  112. Ter Braak, C. J. F. 2006, Stat. Comput., 16, 239 [Google Scholar]
  113. Torres, G., Andersen, J., & Giménez, A. 2010, A&ARv, 18, 67 [Google Scholar]
  114. Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  115. van Leeuwen, F. 2007, A&A, 474, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Vats, D., & Knudson, C. 2018, ArXiv e-prints [arXiv:1812.09384] [Google Scholar]
  117. Vogt, S. S., Burt, J., Meschiari, S., et al. 2015, ApJ, 814, 12 [NASA ADS] [CrossRef] [Google Scholar]
  118. White, T. R., Huber, D., Mann, A. W., et al. 2018, MNRAS, 477, 4403 [NASA ADS] [CrossRef] [Google Scholar]
  119. Wilson, D. M., Gillon, M., Hellier, C., et al. 2008, ApJ, 675, L113 [NASA ADS] [CrossRef] [Google Scholar]
  120. Winn, J. N. 2010, ArXiv e-prints [arXiv:1001.2010] [Google Scholar]
  121. Winn, J. N., Holman, M. J., Carter, J. A., et al. 2009, AJ, 137, 3826 [NASA ADS] [CrossRef] [Google Scholar]
  122. Wright, J. T. 2018, Radial Velocities as an Exoplanet Discovery Method, 4 [Google Scholar]
  123. Wright, J. T., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2004, ApJS, 152, 261 [NASA ADS] [CrossRef] [Google Scholar]
  124. Yi, S., Demarque, P., Kim, Y.-C., et al. 2001, ApJS, 136, 417 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Relation between transit duration W and the scaled semi-major axis

For the transit analysis our code computes the Cartesian coordinates of the Keplerian orbit and then it feeds the Mandel & Agol (2002) model with the planet to star sky-projected distance and the size ratio. Those Cartesian coordinates are computed using notably ; relying on transit observables, comes from P, W, b, and dF (besides e and ω if the orbit is eccentric).

Following Winn (2010), establishing the total transit duration W, that is the temporal interval between the first and fourth contact requires us to compute the integral


where r(f) is the ellipse equation (i.e. the star-planet distance) as a function of the true anomaly f


while fI and fIV are the true anomalies at the moment of first and last contact, respectively. By projecting r onto the plane of the sky from Eqs. (53) and (54) by Murray & Correia (2011), we obtain the star-planet distance in the sky plane, that is


By imposing rsky = R + Rp and then solving Eq. (A.3) in terms of f, fI, and fIV are recovered, such that W can be finally computed from (A.1).

It is clear that all this procedure leads to lengthy algebra and to numerical resolution of equations. Nonetheless, there are very useful approximations which relate W to transit observables in a more straightforward way, and the W relation that has been implemented in our algorithm is written as





Equation (A.4) is a revised version of Eq. (8) by Seager & Mallén-Ornelas (2003), where we added the eccentricity-dependent factors (A.5) and (A.6), as suggested by Winn (2010), and we considered the following approximations:


We evaluated the impact of our approximations, also thanks to some reference data of interest that we inferred from the NASA Exoplanet archive. It turned out that values which are computed through (A.4) differ from those derived from the actual formula of Seager & Mallén-Ornelas (2003) by less than 1% if , which holds for more than 97% of the exoplanets that have been confirmed so far. It is even more reassuring that the median value of produces negligible variations on (∼0.05%), while in the rare cases in which , we may register variations on of order of a few percent, which are still not worrying in general since the median relative uncertainty on is ∼10%.

Without paying a high price, the introduced approximations enable us to easily express W as a function of . This is useful in case the transit is very shallow and it is preferable not to rely on it to infer W. As also explained in Sect. 2.3, in this scenario ρ is established from evolutionary models, then is computed from ρ thanks to Kepler’s third law, and finally also W is available by easily inverting Eq. (A.4).

All Tables

Table 1.

HD 219134: input stellar parameters considered in this work.

Table 2.

HD 219134: output stellar parameters derived considering L approach or R approach.

Table 3.

Planets hosted by HD 219134.

Table 4.

Comparison involving masses and radii derived by different authors for HD 219134 system.

Table 5.

WASP-4 planetary system.

Table 6.

WASP-4: masses and radii of both the star and the planet derived by different authors.

All Figures

thumbnail Fig. 1.

HD 219134 location on the colour-magnitude diagram (left panel) and on the log TeffMV diagram (right panel) together with some reference isochrones. Red stands for 1 Gyr models, while blue stands for 10 Gyr models; solid lines are representative of the nominal stellar metallicity (Z = 0.020 → [Fe/H] = 0.11); dashed lines refer to extremely super-solar metallicity Z = 0.070 → [Fe/H] = 0.66.

In the text
thumbnail Fig. 2.

Values of Teff found in the literature for HD 219134. Values on the x-axis correspond to the following references: (1) Valenti & Fischer (2005); (2) Mishenina et al. (2008); (3) Soubiran et al. (2008); (4) Prugniel et al. (2011); (5) Ramírez et al. (2013); (6) Motalebi et al. (2015), EW approach; (7) Motalebi et al. (2015), SPC approach; (8) Our work, L approach; and (9) our work, R approach.

In the text
thumbnail Fig. 3.

Superimposed transiting models on LCs, as derived from the L approach; the R approach gives analogous results. Two top panels: transit of HD 219134 b, two bottom panels: HD 219134 c.

In the text
thumbnail Fig. 4.

Red histogram: posterior R PDF, as derived from the L approach. Blue histrogram: posterior R PDF, as from the R approach. Black Gaussian: interferometric prior on stellar radius that has been imposed in the R approach. The black Gaussian prior essentially overrides the isochronal constraint (which is looser) and drives the output R value. See text for a complete discussion.

In the text
thumbnail Fig. 5.

WASP-4: regression of M vs. (left panel), and of Mp vs. (right panel). A strong correlation is expected because ρ is constrained from the transit, and (Mp, Rp) are related to (M, R), transit parameters and the RV semi-amplitude.

In the text
thumbnail Fig. 6.

Light curves of WASP-4 b. First two panels show the data provided by Winn et al. (2009), who observed the transit in the Sloan z′ band. The other three panels show the observation by Nikolov et al. (2012) in the Sloan g′ band. The red thick line represents the transit model as inferred from our analysis.

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.