© 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 planets¹.

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 $\mathrm{d}F={\left(\frac{R_{\mathrm{p}}}{R_{\star }}\right)}^2$, 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 M_p is negligible if compared with stellar mass M_⋆; otherwise RV measurements are needed to constrain the orbital eccentricity e and the planetary mass M_p 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 R_p 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 T_eff 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 R_p (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 M_⋆–R_⋆ 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)M_⋆–R_⋆ 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), exoplanet², 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 R_p), and a determination of planetary mass M_p 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 M_p is obtained, but not R_p. 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 $log{R}_{\text{HK}}^{\prime }$ (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 dF_occ.

• 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 T₀, that is the time of inferior conjunction, when the true anomaly at transit is f_t = 90° −ω (with ω argument of the periastron).

• The orbital period P.

• The quantities $\sqrt{e}sin\omega$ and $\sqrt{e}cos\omega$, where e is the orbital eccentricity and ω is the argument of the periastron.

• The parameter $K_2=K\sqrt{1-e^2}{P}^{1/3}$, where K is the RV semi-amplitude.

• The quantities $\sqrt{\mathit{v}sin{i}_{\star }}cos\beta$ and $\sqrt{\mathit{v}sin{i}_{\star }}sin\beta$, 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 T_eff 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, dF_occ, b′, W, T₀, P, K₂, T_eff, R_⋆ or M_⋆ are not allowed, any square root arguments must be positive, and the inequalities $|\sqrt{e}cos\omega |<1$, $|\sqrt{e}sin\omega |<1$ 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 $R_{\mathrm{p}}=\sqrt{\mathrm{d}F}{R}_{\star }$ towards higher R_p values. If the slope of the R_p 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 $\frac{R_{\mathrm{p}}}{R_{\star }}$, as proposed for example by Eastman et al. (2013).

The impact parameter in case of circular orbit is given by b′=a cos i_p/R_⋆, where a is the semi-major axis of planetary orbit and i_p 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 $b=b\prime \frac{1-e^2}{1+esin\omega }$ is preferable to minimise the correlation between jump parameters.

Our algorithm can also handle eccentric orbits by enabling the stepping in $\sqrt{e}sin\omega$ and $\sqrt{e}cos\omega$. 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 a₀ = 0.04 AU and an initial orbital eccentricity e₀ = 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 $\frac{a}{R_{\star }}$ 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 e_min on the eccentricity can be set if the observed pattern resembles a pattern for which $\mathit{g}(e,\omega )={\left(\frac{{\rho }_{\star ,\mathrm{circ}}}{{\rho }_{\star ,\mathrm{th}}}\right)}^{1/3}\ne 1$; 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 = e_min implies a transit at periastron (resp. apoastron) if ρ_⋆,th <  ρ_⋆,circ (resp. ρ_⋆,th >  ρ_⋆,circ). Thus, launching a second MCMCI run where fixed e = e_min 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 K₂ 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 $K_2={(2\pi G)}^{1/3}\frac{M_{\mathrm{p}}sin{i}_{\mathrm{p}}}{M_{\star }^{2/3}}{(1+\frac{M_{\mathrm{p}}}{M_{\star }})}^{-2/3}$, that is it is independent of e and P. In particular, for a given stellar mass, assuming a uniform prior on K₂ implies assuming a uniform prior on M_p sin i_p.

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 u₁ and u₂, but the linear combinations c₁ = 2u₁ + u₂ and c₂ = u₁ − 2u₂ 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, T_eff 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

$$\begin{array}{c}\hfill {\displaystyle {\chi }_{\mathrm{ph}}^2=\sum _{i=1}^{n_{\mathrm{LC}}}\sum _{j=1}^{{\mathrm{np}}_i}{\left(\frac{f_{\mathit{ij}}-{\overline{f}}_{\mathit{ij}}}{{\sigma }_{f_{\mathit{ij}}}}\right)}^2,}\end{array}$$(1)

where n_LC is the number of available LCs, np_i are the number of points of the ith LC, f_i, and σ_fi are the vectors of flux measurements and its errors of the ith LC, while $\overline{f_i}$ 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 ${\chi }_{\text{rv}}^2$ referring to the RV time series

$$\begin{array}{c}\hfill {\displaystyle {\chi }_{\mathrm{rv}}^2=\sum _{i=1}^{n_{\mathrm{RV}}}\sum _{j=1}^{{\mathrm{nr}}_i}{\left(\frac{v_{\mathit{ij}}-{\overline{v}}_{\mathit{ij}}}{{\sigma }_{v_{\mathit{ij}}}}\right)}^2·}\end{array}$$(2)

Summing together ${\chi }_{\text{ph}}^2$ and ${\chi }_{\text{rv}}^2$, 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:

$$\begin{array}{c}\hfill \mathrm{BP}={\left(\frac{x-\overline{x}}{{\sigma }_{\overline{x}}}\right)}^2,\end{array}$$(3)

where $\overline{x}$ and ${\sigma }_{\overline{x}}$ are the values with their uncertainties that are attributed to the generic prior parameter, while x comes from the perturbation of $\overline{x}$ at the generic step of the MCMC. In particular, in all the analyses we carried out, we always considered BPs on c₁ and c₂ (computed from the u₁ and u₂ LD coefficients). Inference of the initial values for u₁ and u₂ (with their respective errors) from Claret’s tables follows the same procedure as described in Gillon et al. (2009a). After that, BP_c1 and BP_c2 control the u₁ and u₂ 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

$$\begin{array}{c}\hfill {\displaystyle {\chi }^2={\chi }_{\mathrm{ph}}^2+{\chi }_{\mathrm{rv}}^2+\sum _{p=1}^{n_{\mathrm{pp}}}{\mathrm{BP}}_{\mathrm{p}},}\end{array}$$(4)

where n_pp 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 χ² 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 $exp(-\frac{1}{2}\mathrm{\Delta }{\chi }^2)$, where Δχ² is the χ² 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):

$$\begin{array}{c}\hfill \mathrm{BIC}={\widehat{\chi }}^2+klogN,\end{array}$$(5)

where ${\widehat{\chi }}^2$ is the smallest χ² 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 N_c is the number of chains, L_c is the length of each chain, and z_ic is the ith draw of the generic z parameter at the cth chain, the mean z value within the cth chain is written as

$$\begin{array}{c}\hfill {\displaystyle {\overline{z}}_{·c}=\frac{1}{L_{\mathrm{c}}}\sum _{i=1}^{L_{\mathrm{c}}}{z}_{\mathit{ic}},}\end{array}$$(6)

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

$$\begin{array}{c}\hfill {\displaystyle {\overline{z}}_{··}=\frac{1}{N_{\mathrm{c}}}\sum _{c=1}^{N_{\mathrm{c}}}{\overline{z}}_{·c}.}\end{array}$$(7)

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 ${\overline{z}}_{·c}$ mean values of each chain, that is

$$\begin{array}{cc}& {\displaystyle W(z)=\frac{1}{N_{\mathrm{c}}}\sum _{c=1}^{N_{\mathrm{c}}}\frac{1}{L_{\mathrm{c}}-1}\sum _{i=1}^{L_{\mathrm{c}}}{(z_{\mathit{ic}}-{\overline{z}}_{·c})}^2}\hfill \end{array}$$(8)

$$\begin{array}{cc}& {\displaystyle B(z)=\frac{L_{\mathrm{c}}}{N_{\mathrm{c}}-1}\sum _{c=1}^{N_{\mathrm{c}}}{({\overline{z}}_{·c}-{\overline{z}}_{··})}^2.}\hfill \end{array}$$(9)

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

$$\begin{array}{c}\hfill {\widehat{\mathrm{var}}}^{+}(z)=\frac{L_{\mathrm{c}}-1}{L_{\mathrm{c}}}W(z)+\frac{1}{L_{\mathrm{c}}}B(z).\end{array}$$(10)

The Gelman-Rubin statistic is commonly denoted by $\widehat{R}$, the estimator is defined as

$$\begin{array}{c}\hfill \widehat{R}(z)=\sqrt{\frac{{\widehat{\mathrm{var}}}^{+}(z)}{W(z)}}\end{array}$$(11)

and convergence is reached when $\widehat{R}$ 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 $\widehat{R}$ cut-off values varying from 1.003 to 1.3. Among these papers, 43% of the authors (the relative majority) consider $\widehat{R}=1.1$, which is the recommended value by Gelman & Rubin (1992).

Our own empirical practice for convergence is considering $\widehat{R}\lesssim 1.05$. 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

$$\begin{array}{c}\hfill \widehat{T}(z)=N_c{L}_{\mathrm{c}}min(\frac{{\widehat{\mathrm{var}}}^{+}(z)}{B(z)},1)\end{array}$$(12)

with the estimator $\widehat{T}(z)\ge 1000$.

Since $\widehat{R}(z)>1$, from (11) we infer

$$\begin{array}{c}\hfill W(z)<{\widehat{\mathrm{var}}}^{+}(z).\end{array}$$(13)

The quantity ${\widehat{\mathrm{var}}}^{+}(z)$ 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

$$\begin{array}{c}\hfill W(z)<{\widehat{\mathrm{var}}}^{+}(z)<B(z)\Rightarrow {\widehat{\mathrm{var}}}^{+}(z)<B(z).\end{array}$$(14)

Condition (14) implies

$$\begin{array}{c}\hfill min(\frac{{\widehat{\mathrm{var}}}^{+}(z)}{B(z)},1)=\frac{{\widehat{\mathrm{var}}}^{+}(z)}{B(z)}\Rightarrow \widehat{T}(z)=N_c{L}_{\mathrm{c}}\frac{{\widehat{\mathrm{var}}}^{+}(z)}{B(z)}·\end{array}$$(15)

By recovering B(z) from Eq. (10) and substituting it in Eq. (15), we can express $\widehat{T}(z)$ as a function of $\widehat{R}(z)$ as follows:

$$\begin{array}{cc}\hfill \widehat{T}(z)& =N_c{L}_{\mathrm{c}}\frac{{\widehat{\mathrm{var}}}^{+}(z)}{L_{\mathrm{c}}{\widehat{\mathrm{var}}}^{+}(z)-(L_{\mathrm{c}}-1)W(z)}\hfill \\ \hfill & =N_c{L}_{\mathrm{c}}\frac{1}{L_{\mathrm{c}}-(L_{\mathrm{c}}-1)\frac{W(z)}{{\widehat{\mathrm{var}}}^{+}(z)}}\hfill \\ \hfill & =N_c{L}_{\mathrm{c}}\frac{1}{L_{\mathrm{c}}-\frac{(L_{\mathrm{c}}-1)}{{\widehat{R}}^2(z)}},\hfill \end{array}$$(16)

where we used the definition of $\widehat{R}(z)$ given by (11) in the last passage. By inverting (16), we can finally express $\widehat{R}(z)$ as a function of $\widehat{T}(z)$, that is

$$\begin{array}{c}\hfill \widehat{R}(z)=\sqrt{\frac{L_{\mathrm{c}}-1}{L_{\mathrm{c}}(1-\frac{N_{\mathrm{c}}}{\widehat{T}(z)})}}·\end{array}$$(17)

To give some numbers, if we consider N_c = 5 and L_c = 80 000³, achieving the goal suggested by Ford (2006) ($\widehat{T}(z)\ge 1000$) would require $\widehat{R}(z)\le 1.0025$. According to our experience, this $\widehat{R}(z)$ 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 $\widehat{R}(z)\lesssim 1.05$, 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

$$\begin{array}{c}\hfill {\beta }_{r,\mathrm{bin}}=\frac{{\sigma }_N}{{\sigma }_1}\sqrt{\frac{N_{\mathrm{bin}}(M-1)}{M}},\end{array}$$(18)

where N_bin is the mean number of points in each bin, M the number of bins, and σ₁ 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

$$\begin{array}{c}\hfill \mathrm{CF}={\beta }_w·{\beta }_r\end{array}$$(19)

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 interface⁴.

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 Z_ini 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 Z_ini 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 M_ini, ⋆, Z_ini) lists the sequence of evolutionary states through which a star having an initial mass M_ini, ⋆ and an initial metallicity Z_ini 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 T_eff 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 P_rot, $log{R}_{\text{HK}}^{\prime }$, 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_⋆/P_rot, $log{R}_{\text{HK}}^{\prime }$ 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 P_rot (gyrochronology) or $log{R}_{\text{HK}}^{\prime }$ (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 $log{R}_{\text{HK}}^{\prime }$-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):

$$\begin{array}{c}\hfill [\mathrm{Y}/\mathrm{Mg}]=0.170-0.0371t,\end{array}$$(20)

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_⋆, $log{R}_{\text{HK}}^{\prime }$ 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 T_eff variation ∼0.1% over the whole age of the universe, which is lower than any reasonable input T_eff 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

$$\begin{array}{c}\hfill Z={10}^{[\mathrm{Fe}/\mathrm{H}]-1.817}\end{array}$$(21)

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 T_eff and ρ_⋆ (or log g) are available as inputs, we work in the modified HRD log T_eff– log ρ (or log T_eff–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 T_eff. 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

$$\begin{array}{c}\hfill \mathcal{G}(y_i,T_{\mathrm{eff},i})=\frac{e^{-\frac{1}{2}[{\left(\frac{log{T}_{\mathrm{eff},i}-log{T}_{\mathrm{eff}}}{\mathrm{\Delta }log{T}_{\mathrm{eff}}}\right)}^2+{\left(\frac{log{y}_{\star ,i}-log{y}_{\star }}{\mathrm{\Delta }log{y}_{\star }}\right)}^2]}}{2\pi \mathrm{\Delta }log{T}_{\mathrm{eff}}\mathrm{\Delta }log{y}_{\star }},\end{array}$$(22)

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

$$\begin{array}{cc}\hfill p_i=& [{\left(\frac{log{T}_{\mathrm{eff}}-log{T}_{\mathrm{eff},i}}{\mathrm{\Delta }log{T}_{\mathrm{eff}}}\right)}^2+{\left(\frac{y_{\star }-y_{\star ,i}}{\mathrm{\Delta }{y}_{\star }}\right)}^2+\sum _j{\left(\frac{{\stackrel{\sim }{y}}_{j\star }-{\stackrel{\sim }{y}}_{j\star ,i}}{\mathrm{\Delta }{\stackrel{\sim }{y}}_{j\star }}\right)}^2\hfill \\ \hfill & {+{log}^2\left(\frac{v_{\mathrm{ref}}}{v_{\mathrm{evo}}}\right)]}^{-1}\hfill \end{array}$$(23)

to be attributed to the ith row of Is’, where y_⋆ has the same meaning as at the previous point, while ${\stackrel{\sim }{\mathit{y}}}_{j\star }$ may represent, if directly available, the input M_⋆ and log g (being ${\stackrel{\sim }{\mathit{y}}}_{j\star }\ne {\mathit{y}}_{\star }\forall j$). The weight takes into account both the geometrical distance between theoretical and observational parameters and the evolutionary speed v_evo of a star in that particular region of the HRD, that is

$$\begin{array}{c}\hfill v_{\mathrm{evo}}=\sqrt{{\left(\frac{y_2-y_1}{t_2-t_1}\right)}^2+{\left(\frac{log{T}_{\mathrm{eff},2}-log{T}_{\mathrm{eff},1}}{t_2-t_1}\right)}^2},\end{array}$$(24)

where (log T_eff, 1, y₁) and (log T_eff, 2, y₂) are two consecutive points with ages t₁ and t₂ 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 v_ref is a normalisation factor corresponding to the minimum evolutionary speed that has been detected along the whole track.

The less v_evo, 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 v_evo decreases.

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

$$\begin{array}{cc}& X_{\star }=\frac{{\displaystyle \sum _{i=1}^{\overline{n}}{X}_i\mathcal{G}(L_i,T_{\mathrm{eff},i})p_i}}{{\displaystyle \sum _{i=1}^{\overline{n}}\mathcal{G}(L_i,T_{\mathrm{eff},i})p_i}}·\hfill \end{array}$$(25)

$$\begin{array}{cc}& \mathrm{\Delta }{X}_{\star }=\sqrt{\frac{{\displaystyle \sum _{i=1}^{\overline{n}}{(X_i-X_{\star })}^2\mathcal{G}(L_i,T_{\mathrm{eff},i})p_i}}{{\displaystyle \sum _{i=1}^{\overline{n}}\mathcal{G}(L_i,T_{\mathrm{eff},i})p_i}}},\hfill \end{array}$$(26)

where $\overline{n}$ 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 Z_ini, 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 Z_k, l = Z_k, l(t), for any given stellar mass and Z_ini (represented through the subscripts l and k, respectively). In this way, after an iterative process, it is possible to infer what the Z_ini 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 (t_1, ⋆) and mass (M_1, ⋆). After computing Z_{k, M1, ⋆}(t = t_1, ⋆) for all the available curves characterised by different Z_ini (identified by the subscript k), we identify the $\widehat{k}$th curve such that $Z_{\widehat{k},M_{1,\star }}(t=t_{1,\star })=Z_{\star }$, that is the first guess for the initial metallicity of the star Z_1, ini,⋆.

By loading the grid of isochrones labelled Z_1, 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 Z_ini, ⋆ 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 MCMCI⁵ 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], T_eff (or the colour index), and M_⋆, where [Fe/H] and T_eff (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 ∼10⁵ 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], T_eff, M_⋆, and W. From W_step, ${\frac{a}{R_{\star }}}_{\mathrm{step}}$ is computed through Eq. (A.4). From Kepler’s third law, ${\frac{a}{R_{\star }}}_{\mathrm{step}}$ 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 T_eff, step), to derive the theoretical M_⋆,isoch at that step, according to Eq. (25). Finally, the Bayesian penalty

$$\begin{array}{c}\hfill {\mathrm{BP}}_{\mathrm{M}}={\left(\frac{M_{\star ,\mathrm{step}}-M_{\star ,\mathrm{isoch}}}{\mathrm{\Delta }{M}_{\star }}\right)}^2\end{array}$$(27)

(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], T_eff 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, T_eff, 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:

$$\begin{array}{c}\hfill \mathrm{B}{P}_M+\mathrm{B}{P}_R={\left(\frac{M_{\star ,\mathrm{step}}-M_{\star ,\mathrm{isoch}}}{\mathrm{\Delta }{M}_{\star }}\right)}^2+{\left(\frac{R_{\star ,\mathrm{step}}-R_{\star ,\mathrm{isoch}}}{\mathrm{\Delta }{R}_{\star }}\right)}^2\end{array}$$(28)

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, ${\frac{a}{R_{\star }}}_{\mathrm{step}}$ is computed through Kepler’s third law and finally W_step is expressed as a function of ${\frac{a}{R_{\star }}}_{\mathrm{step}}$ 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 R_p, step (from dF_step, if any LC is available) and M_p, 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 T_eff among input parameters, but no convergence was reached in theoretical models.

Actually, as shown in Fig. 1 (right panel), the T_eff 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 T_eff = 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 T_eff.

Fig. 1. HD 219134 location on the colour-magnitude diagram (left panel) and on the log Teff–MV 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.

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, T₀, P, and K₂ 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, P_rot, $log{R}_{\text{HK}}^{\prime }$, 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, P_rot, $log{R}_{\text{HK}}^{\prime }$, input stellar radius R_prior from interferometry. In this case BPs entering the merit function involved [Fe/H], B − V and R_prior. 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 c₁ and c₂ were properly computed to build their corresponding BPs to be added to the merit function. Furthermore, also the sum BP_M + BP_M 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.

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.

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 R_prior 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 M_p and R_p 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.

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 − V_step, L_⋆, step) as inputs for isochrone fitting and assuming a prior on the parallax-based stellar luminosity L_prior. Instead, the blue histogram is the output R_⋆ PDF that has been inferred considering ([Fe/H]_step, B − V_step, R_⋆, step) as inputs for isochrone fitting, and assuming a prior (black Gaussian) on the interferometric radius R_prior.

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.

Fig. 5. WASP-4: regression of M⋆ vs. $R_{\star }^2$ (left panel), and of Mp vs. $R_{\text{p}}^3$ (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 T_eff and [Fe/H]. Actually, if R_⋆step is assumed on input, the isochrone-based Bayesian penalty BP_R alone (Eq. (28)) would drive to an output radius R_out, 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 T_eff 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, T₀, and P as jump parameters and, in addition, we set T_eff and [Fe/H] as priors. Thus BPs on T_eff and [Fe/H] were considered, besides those on the c₁ and c₂ LD coefficients (qd model assumed). Moreover, BP_M 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.

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 transit⁶ 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 (M_p, R_p), respectively, thus according to error propagation, we could expect their relative uncertainties to be higher than those affecting M_⋆, R_⋆, M_p, and R_p. But this is not the case since we have a derived output uncertainty ${\frac{\mathrm{\Delta }{g}_{\star }}{g_{\star }}|}_{\mathrm{out}}=0.020$ versus ${\frac{\mathrm{\Delta }{g}_{\star }}{g_{\star }}|}_{\mathrm{ep}}=0.063$ that should be expected from error propagation; similarly, we have ${\frac{\mathrm{\Delta }{\rho }_{\mathrm{p}}}{{\rho }_{\mathrm{p}}}|}_{\mathrm{out}}=0.026$ versus ${\frac{\mathrm{\Delta }{\rho }_{\mathrm{p}}}{{\rho }_{\mathrm{p}}}|}_{\mathrm{ep}}=0.065$.

However, uncertainties in agreement with standard error propagation are expected only if the pairs (M_⋆, R_⋆) and (M_p, R_p) are uncorrelated. Actually, a correlation between M_⋆ and R_⋆ exists because, besides isochrones, we have a strong constraint on ${\rho }_{\star }\sim \frac{M_{\star }}{R_{\star }^3}$ thanks to the transit. In addition, M_p (resp. R_p) differ from M_⋆ (resp. R_⋆) by factors that are constrained from the transit and the RV semi-amplitude; thus a correlation between M_p and R_p is also expected, although it is a bit more diluted. What has just been described is shown in Fig. 5, where the M_⋆ (resp. M_p) values coming from all the chain steps are plotted versus the corresponding $R_{\star }^2$ (resp. $R_{\text{p}}^3$) values. Linear correlations are clear and the coefficients of determination are $r_{\star }^2=0.91$ (left panel) and $r_{\text{p}}^2=0.87$ (right panel).

If a strict linear correlation held (r² = 1), no dispersion in $log{g}_{\star }\sim \frac{M_{\star }}{R_{\star }^2}$ and in ${\rho }_{\mathrm{p}}\sim \frac{M_{\mathrm{p}}}{R_{\mathrm{p}}^3}$ would be registered, thus both their uncertainties would be zero. The equation ${\sigma }_{\mathrm{ne}}=\sqrt{1-r^2}$ 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:

$$\begin{array}{cc}& {\frac{\mathrm{\Delta }{g}_{\star }}{g_{\star }}|}_{\exp }={\frac{\mathrm{\Delta }{g}_{\star }}{g_{\star }}|}_{\mathrm{ep}}·\sqrt{1-r_{\star }^2}=0.019\approx {\frac{\mathrm{\Delta }{g}_{\star }}{g_{\star }}|}_{\mathrm{out}}\hfill \end{array}$$(29)

$$\begin{array}{cc}& {\frac{\mathrm{\Delta }{\rho }_{\mathrm{p}}}{{\rho }_{\mathrm{p}}}|}_{\exp }={\frac{\mathrm{\Delta }{\rho }_{\mathrm{p}}}{{\rho }_{\mathrm{p}}}|}_{\mathrm{ep}}·\sqrt{1-r_{\mathrm{p}}^2}=0.024\approx {\frac{\mathrm{\Delta }{\rho }_{\mathrm{p}}}{{\rho }_{\mathrm{p}}}|}_{\mathrm{out}}·\hfill \end{array}$$(30)

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 M_p and R_p values, if compared to what is stated in the literature.

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

$$\begin{array}{c}\hfill Y=Y_{\mathrm{p}}+\frac{\mathrm{\Delta }Y}{\mathrm{\Delta }Z}Z,\end{array}$$(31)

where Y_p = 0.2485 (Komatsu et al. 2011) is the primordial helium abundance, while the helium-to-metal enrichment ratio $\frac{\mathrm{\Delta }Y}{\mathrm{\Delta }Z}=1.78$ is derived from solar calibration. Instead, $\frac{\mathrm{\Delta }Y}{\mathrm{\Delta }Z}$ 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 T_eff, [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 T_eff 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.

Acknowledgments

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.

References

Appendix A: Relation between transit duration W and the scaled semi-major axis $\frac{a}{R_{\star }}$

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 $\frac{R_{\mathrm{p}}}{R_{\star }}$ size ratio. Those Cartesian coordinates are computed using notably $\frac{a}{R_{\star }}$; relying on transit observables, $\frac{a}{R_{\star }}$ 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

$$\begin{array}{c}\hfill W=\frac{P}{2\pi \sqrt{1-e^2}}{\displaystyle {\int }_{f_{\mathrm{I}}}^{f_{\mathrm{IV}}}\frac{r^2(f)}{a^2}\mathrm{d}f,}\end{array}$$(A.1)

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

$$\begin{array}{c}\hfill r=\frac{a(1-e^2)}{1+ecosf},\end{array}$$(A.2)

while f_I and f_IV 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

$$\begin{array}{c}\hfill r_{\mathrm{sky}}=\frac{a(1-e^2)}{1+ecosf}\sqrt{1-{sin}^2(\omega +f)sin{i}_{\mathrm{p}}}.\end{array}$$(A.3)

By imposing r_sky = R_⋆ + R_p and then solving Eq. (A.3) in terms of f, f_I, and f_IV 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–$\frac{a}{R_{\star }}$ relation that has been implemented in our algorithm is written as

$$\begin{array}{c}\hfill \frac{a}{R_{\star }}=\frac{P}{\pi W}{\stackrel{\sim }{e}}_1\sqrt{{(1+\sqrt{\mathrm{d}F})}^2-{(b^{\prime }{\stackrel{\sim }{e}}_2)}^2},\end{array}$$(A.4)

where

$$\begin{array}{cc}\hfill {\stackrel{\sim }{e}}_1& =\frac{\sqrt{1-e^2}}{1+esin\omega }\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill {\stackrel{\sim }{e}}_2& =\frac{1-e^2}{1+esin\omega }·\hfill \end{array}$$(A.6)

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:

$$\begin{array}{c}\hfill sin\frac{\pi W}{P}\approx \frac{\pi W}{P}{cos}^2\frac{\pi W}{P}\approx 1.\end{array}$$(A.7)

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 $\frac{a}{R_{\star }}$ 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 $\frac{W}{P}<0.07$, 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 $\frac{W}{P}=0.015$ produces negligible variations on $\frac{a}{R_{\star }}$ (∼0.05%), while in the rare cases in which $\frac{W}{P}>0.07$, we may register variations on $\frac{a}{R_{\star }}$ of order of a few percent, which are still not worrying in general since the median relative uncertainty on $\frac{a}{R_{\star }}$ is ∼10%.

Without paying a high price, the introduced approximations enable us to easily express W as a function of $\frac{a}{R_{\star }}$. 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 $\frac{a}{R_{\star }}$ is computed from ρ_⋆ thanks to Kepler’s third law, and finally also W is available by easily inverting Eq. (A.4).

All Tables

All Figures

Fig. 1. HD 219134 location on the colour-magnitude diagram (left panel) and on the log Teff–MV 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

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

	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

	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

	Fig. 5. WASP-4: regression of M⋆ vs. $R_{\star }^2$ (left panel), and of Mp vs. $R_{\text{p}}^3$ (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

	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