Issue 
A&A
Volume 635, March 2020



Article Number  A6  
Number of page(s)  17  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201936326  
Published online  28 February 2020 
MCMCI: A code to fully characterise an exoplanetary system^{⋆}
^{1}
Space Sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, 19C Allée du 6 Août, 4000 Liège, Belgium
email: a.bonfanti@uliege.be
^{2}
Astrobiology Research Unit, Université de Liège, Allée du 6 Août 19, 4000 Liège, Belgium
Received:
16
July
2019
Accepted:
29
December
2019
Context. Useful information can be retrieved by analysing the transit light curve of a planethosting star or induced radial velocity oscillations. However, inferring the physical parameters of the planet, such as mass, size, and semimajor axis, requires preliminary knowledge of some parameters of the host star, especially its mass or radius, which are generally inferred through theoretical evolutionary models.
Aims. We seek to present and test a whole algorithm devoted to the complete characterisation of an exoplanetary system thanks to the global analysis of photometric or radial velocity time series combined with observational stellar parameters derived either from spectroscopy or photometry.
Methods. We developed an integrated tool called MCMCI. This tool combines the Markov chain Monte Carlo (MCMC) approach of analysing photometric or radial velocity time series with a proper interpolation within stellar evolutionary isochrones and tracks, known as isochrone placement, to be performed at each chain step, to retrieve stellar theoretical parameters such as age, mass, and radius.
Results. We tested the MCMCI on the HD 219134 multiplanetary system hosting two transiting rocky super Earths and on WASP4, which hosts a bloated hot Jupiter. Even considering different input approaches, a final convergence was reached within the code, we found good agreement with the results already stated in the literature and we obtained more precise output parameters, especially concerning planetary masses.
Conclusions. The MCMCI tool offers the opportunity to perform an integrated analysis of an exoplanetary system without splitting it into the preliminary stellar characterisation through theoretical models. Rather this approach favours a close interaction between light curve analysis and isochrones, so that the parameters recovered at each step of the MCMC enter as inputs for purposes of isochrone placement.
Key words: planets and satellites: fundamental parameters / stars: fundamental parameters
Data and code are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/635/A6
© ESO 2020
1. Introduction
A few decades have passed since the discovery of the first exoplanet orbiting the mainsequence (MS) star 51 Pegasi (Mayor & Queloz 1995) and by now ∼4000 exoplanets have been confirmed, ∼75% of which are transiting planets^{1}.
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 planethosting star is available, the transit depth gives the squared ratio of planetary to stellar radii , while the orbital period P and the transit shape parameters (impact parameter, depth, and duration) enable the retrieval of the stellar mean density ρ_{⋆}. This is true for circular orbits under the assumption that exoplanetary mass 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énOrnelas (2003). Thus, as a very basic example, LC analyses are usually split into three steps: first the LCanalyser 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 transitanalysis 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 customdeveloped 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 isochronebased 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 massradius 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 wellconstrained 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 YonseiYale (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 massradius 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^{2}, 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 (ForemanMackey 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 isochronebased joint analysis and may also give different constraints for the stellar age (and thus the age of the entire exoplanetary system) using both modeldependent 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 MetropolisHastings 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 burnin. In this way the effect of initial values on the posterior inference is minimised. Choice of the burnin length depends upon the initial state and the speed of convergence to the limiting distribution. Establishing the proper burnin length requires analysis of the output case by case, however experience suggests that setting the burnin 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 burnin 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 semiamplitude 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 RossiterMcLaughlin 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 intrapixel sensitivity of the detector (Knutson et al. 2008; Charbonneau et al. 2008); the socalled 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 timedependent 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, crosscorrelation function width and bisector parameters (Baranne et al. 1996), and (a typical indicator of stellar magnetic activity; see e.g. Wright et al. 2004, and references therein).
Barycentric Dynamic Time (TDB) is the default time standard that is used by the code in all the analyses; TDB is suitable for precision time monitoring. This time standard is preferable to the Coordinated Universal Time (UTC) because it avoids the drift due to leap seconds; in addition, it refers to the barycentre of the Solar System, thus it corrects for relativistic effect due to the gravitational potential of Earth. A discussion about time standards and timing precision is provided in for example Eastman et al. (2010).
The following set of parameters may be randomly perturbed at each chain step (jump parameters or step parameters).

The transit depth dF.

The occultation depth 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_{0}, 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 and , where e is the orbital eccentricity and ω is the argument of the periastron.

The parameter , where K is the RV semiamplitude.

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

The stellar metallicity [Fe/H]_{⋆}.

The stellar effective temperature T_{eff} or colour index.

The stellar radius R_{⋆}.

The stellar mass M_{⋆}.

The limbdarkening (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_{0}, P, K_{2}, T_{eff}, R_{⋆} or M_{⋆} are not allowed, any square root arguments must be positive, and the inequalities , must hold since we are dealing with bound orbital systems. If any of the drawn value is not physical, the jump is not accepted, and a copy of the previous state is made.
2.1.2. Comments on jump parameters
Stepping in dF, that is assuming a uniform prior on dF, creates an implicit bias on towards higher 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 , 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 semimajor 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 is preferable to minimise the correlation between jump parameters.
Our algorithm can also handle eccentric orbits by enabling the stepping in and . This parametrisation ensures both orthogonality and uniform priors on e and ω, which has been pointed out for the first time by Anderson et al. (2011) and then broadly discussed by Eastman et al. (2013). In case RV data are available, the usual approach is to launch a first run with a fixed e = 0 value, and then launch another run letting e as a free jump parameter and check whether the BIC (Eq. (5), see later) is in favour of an eccentric or circular solution.
Instead, in case of a transitonly fit, treating e is generally problematic. On the one hand, orbits of closein exoplanets are likely circular. Dawson & Johnson (2018) report that for orbital periods P < 3 days (i.e. orbital semimajor axes a < 0.04 AU, if a Sunlike 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 semimajor axis a_{0} = 0.04 AU and an initial orbital eccentricity e_{0} = 0.1, the orbital circularisation happens after ∼150 Myr, ∼100 Myr, and ∼30 Myr for a hot Jupiter, a hot Neptune, and a superEarth, all orbiting a Sunlike 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 semimajor 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 transitbased 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 Jupitersized 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 ; 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 transitonly 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; IglesiasMarzoa et al. 2015; Kreidberg 2015).
The parameter K_{2} is used as a jump parameter instead of K to minimise the correlation between jump parameters (Gillon et al. 2009a). According to its definition, we note that , that is it is independent of e and P. In particular, for a given stellar mass, assuming a uniform prior on K_{2} implies assuming a uniform prior on M_{p} sin i_{p}.
In case of a multiplanetary system, so far only the β value of one single planet is handled by the code. We note that β can be inferred from the RossiterMcLaughlin 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.
Limbdarkening coefficients can be either four (notlinear 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_{1} and u_{2}, but the linear combinations c_{1} = 2u_{1} + u_{2} and c_{2} = u_{1} − 2u_{2} 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
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 are the fluxes that have been computed according to the eclipse model in use. Adopting the same kind of notation, we can also introduce the analogous referring to the RV time series
Summing together and , we obtain the basic merit function. In addition, the jump parameters specified above – besides the stellar luminosity L_{⋆}, surface gravity log g, density ρ_{⋆}, and projected rotational velocity v sin i_{⋆} – may also be set as priors. Those parameters that have been set as Gaussian priors enter the merit function, as well, in the form of the following Bayesian penalty addendum:
where and are the values with their uncertainties that are attributed to the generic prior parameter, while x comes from the perturbation of at the generic step of the MCMC. In particular, in all the analyses we carried out, we always considered BPs on c_{1} and c_{2} (computed from the u_{1} and u_{2} LD coefficients). Inference of the initial values for u_{1} and u_{2} (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_{1} and u_{2} floating during the MCMC steps, to obtain a LD solution consistent with theory.
Finally, at each step, it is possible to define the complete merit function, which is given by
where 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 χ^{2} at that step is lower than the one computed at the previous step. If this does not hold, the new state is accepted with probability , where Δχ^{2} is the χ^{2} difference between the two last steps. In all the other cases in which the new state is not accepted, the state is pushed back to the previous set of parameters.
Generation of a new state x′ from the present state x is controlled by a transition probability function and the optimal choice for this function would be considering the posterior probability distribution. However, the aim of the MCMC is actually to retrieve the posterior distribution itself, which is unknown at the beginning. Thus, as a common choice, we considered a Gaussian distribution centred around x. To make the MCMC efficiently converge, the step size, that is the Gaussian variance, to perform the jump from x to x′ must be tuned up. This is done by computing the acceptance rate, that is the fraction of accepted states over a window spanning a given number of steps, which we set at 100. As reported by Ford (2005), if x has only one dimension, the optimal acceptance rate is ∼0.44, while if x has many dimensions, the optimal acceptance rate is ∼0.25, as has been proven by Roberts et al. (1997).
We implemented the Gibbs sampling during the burnin phase, when just a single jump parameter changes at each step; thus, for tuningup the step size, during the burnin phase the acceptance rate is set to 40% if the Gibbs sampler is switched on. After the burnin phase, all the jump parameters are modified at a time, therefore in this phase we monitor the acceptance rate imposing an optimal value equal to 25%.
If several sets of models that may differ in terms of input parameters or baseline choice have to be compared, the adopted criterion to choose the best model is given by the Bayesian information criterion (Schwarz 1978):
where is the smallest χ^{2} found in the Markov chains, k is the number of free parameters of a given model, and N is the number of data points. The smaller the BIC, the better the model.
At the end, the PDFs of the relevant stellar and planetary parameters are built and meaningful statistics can be retrieved.
2.1.4. Convergence discussion
It is worth launching several chains for a given process and then checking their mutual convergence thanks to the test by Gelman & Rubin (1992), which we briefly present following Ford (2006). If 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
while the global mean of z considering the single mean values computed from each chain is written as
It is then possible to compute W(z) as the average of the variances within each chain, and B(z) as the variance based on the mean values of each chain, that is
An unbiased estimator of the variance of z is then given by the following weighted average between W(z) and B(z):
The GelmanRubin statistic is commonly denoted by , the estimator is defined as
and convergence is reached when is close to 1. Different threshold values indicating convergence are adopted in the literature: randomly picking up 100 papers that cite Gelman & Rubin (1992) in 2017, Vats & Knudson (2018) find cutoff values varying from 1.003 to 1.3. Among these papers, 43% of the authors (the relative majority) consider , which is the recommended value by Gelman & Rubin (1992).
Our own empirical practice for convergence is considering . This criterion may be easily satisfied in case of one single planet fitting, while a slight tension among jump parameters may arise in case of a multiplanetary system. To check convergence, Ford (2006) also recommends a minimum number of effective independent draws, which can be estimated through
with the estimator .
Since , from (11) we infer
The quantity comes from an average between W(z) and B(z), therefore its value is between W(z) and B(z); in particular, given condition (13), we have
Condition (14) implies
By recovering B(z) from Eq. (10) and substituting it in Eq. (15), we can express as a function of as follows:
where we used the definition of given by (11) in the last passage. By inverting (16), we can finally express as a function of , that is
To give some numbers, if we consider N_{c} = 5 and L_{c} = 80 000^{3}, achieving the goal suggested by Ford (2006) () would require . According to our experience, this cutoff appears too demanding in general. We anticipate that this condition is satisfied for all the jump parameters we considered in the analysis of WASP4 (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 closein super Earths (detected both through transit and RV) and by another two, or possibly more, massive planets, whose presence has been inferred from the RV technique. According to all these considerations, we prefer to keep , as a general reference threshold to establish convergence, and we notice that this cutoff 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 overestimation of the white noise. Furthermore, we can take the effect of red noise into account (Pont et al. 2006) by binning the LC on different timescales, and for each binning we then compute
where N_{bin} is the mean number of points in each bin, M the number of bins, and σ_{1} and σ_{N} are the standard deviations of the unbinned and binned residuals, respectively. The maximum among the available β_{r, bin} values is set as the reference β_{r}. Then, the photometric errors provided by the observations are multiplied by the correction factor
and eventually, for RVs, a “jitter” noise may be added quadratically to the error budget. After that a new MCMC run is launched with the updated (rescaled) error bars to obtain reliable error bars on the fitted parameters.
2.2. Isochrone placement
The isochrone placement technique aims to infer stellar parameters – above all radius R_{⋆}, mass M_{⋆}, and age t_{⋆} – by interpolating a set of observational input parameters in a grid of stellar evolutionary models; these observational input parameters come from spectroscopy, transit constraints, photometry and, in particular, they may also involve Gaia parallaxes and magnitudes. The grids of theoretical models in use are currently based upon the PARSEC v.1.2S (Chen et al. 2014; Bressan et al. 2012) + COLIBRI PR16 (Marigo et al. 2013; Rosenfield et al. 2016) stellar isochrones and tracks, which are globally presented in Marigo et al. (2017) and available for download through the CMD web interface^{4}.
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 preMS to postMS 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 subtables, 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}, , 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 JohnsonCousin 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}, and [Y/Mg] optional parameters is to give a preliminary age indication, that may help in disentangling degeneracies involving MS and preMS isochrones, that is isochrone overlap on the HertzsprungRussell 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 preMS isochrones. If this is the case, from a computational point of view, the following interpolation within isochrones considers only those subtables 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 (magnetic activity) is to be applied just to latetype 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 activityage relations become less sensible to age variations in midtoold 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 welldefined periodage relation up to ∼2.5 Gyr. Taking all these considerations into account, Bonfanti et al. (2016) explain how we use the gyrochronological relation by Barnes (2010) and the age relation by Mamajek & Hillenbrand (2008) to compute conservative age lower limits, which they denote by τ_{v} and τ_{HK}, respectively. Just in assessing age lower limits, the code already takes care of the abovementioned caveats such that the user is not obliged to be aware of them, and we attain the goal of possible discarding of unlikely preMS isochrones before interpolation.
Bonfanti et al. (2016) also note that the parameter degeneracies on the HRD under the turnoff (TO) can be broken by using ρ_{⋆}. The mean stellar density is higher in the MS phase than in the preMS, thus it enables us to distinguish MS stars from preMS stars and to set a conservative age lower limit called τ_{ρ}.
The updated implementation of isochrone placement now also accepts [Y/Mg] among the input parameters. Nissen (2016) provide the following [Y/Mg]age relation (t in Gyr; age scatter 0.6 Gyr):
thereby confirming the trend that was first observed by da Silva et al. (2012) for solartype stars. As explained by Feltzing et al. (2017) and references therein, yttrium is produced through slow neutron capture process (sprocess) 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 heliumcore burning phase (Slumstrup et al. 2017). Feltzing et al. (2017) show that the [Y/Mg]t_{⋆} trend is also [Fe/H]dependent and it looses its power of predicting ages for [Fe/H] ≃ −0.5 dex, proving that it is only valid for solar metallicity stars. Therefore, if any [Y/Mg] value is provided on input, first of all our code evaluates whether the stellar metallicity belongs to the −0.2 < [Fe/H] < 0.2 range. If this does not hold, we simply miss the opportunity of removing possible isochronal degeneracies through the [Y/Mg] index, and if neither v sin i_{⋆}, or ρ_{⋆} suggest discarding preMS 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 preMS 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
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 bidimensional Gaussian window function
where y_{⋆} may indicate either log L_{⋆}, log ρ_{⋆}, or log g, depending on the specific HRD we are working in.
Fourth, we define the weight
to be attributed to the ith row of Is’, where y_{⋆} has the same meaning as at the previous point, while may represent, if directly available, the input M_{⋆} and log g (being ). The weight takes into account both the geometrical distance between theoretical and observational parameters and the evolutionary speed v_{evo} of a star in that particular region of the HRD, that is
where (log T_{eff, 1}, y_{1}) and (log T_{eff, 2}, y_{2}) are two consecutive points with ages t_{1} and t_{2} 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 “theoryobservations” 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
where is the number of points (rows) of the grid Is’.
At this point the ED phenomenon, also known as microscopic or atomic diffusion (Burgers 1969; Chapman & Cowling 1970; Chaboyer et al. 2001), enters the picture of the algorithm. Element diffusion is caused by the interaction among various chemical species in stars with a convective envelope, which leads to a progressive surface depletion of elements heavier than hydrogen, which sink downwards.
In particular, it turns out that very low mass stars up to 0.5 M_{⊙} exhibit constant Z during their entire life; this is also the case for stars more massive than ∼2 M_{⊙}. In the intermediate mass range, instead, Z diminishes during the MS phase. After that, as the convective envelope deepens, the consequent mixing of elements makes the star approximately as metallic as it was at its birth (see e.g. Bonfanti et al. 2015, Fig. 3).
As already said, the isochrone grids are identified through 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 th curve such that , 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 WASP4 (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^{5} 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^{5} 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 submatrices 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}, is computed through Eq. (A.4). From Kepler’s third law, is translated into ρ_{⋆,step}, and R_{⋆,step} is derived from M_{⋆,step} and ρ_{⋆,step}. In addition, ρ_{⋆,step} enters isochrone placement (together with [Fe/H]_{step} and T_{eff, step}), to derive the theoretical M_{⋆,isoch} at that step, according to Eq. (25). Finally, the Bayesian penalty
(where ΔM_{⋆} is the uncertainty of M_{⋆,isoch}, that is computed according to Eq. (26)) is added to the merit function given by Eq. (4). In this way, M_{⋆,step} is driven by the theoretical mass value M_{⋆,isoch}, and thus R_{⋆,step} is computed from M_{⋆,step} (constrained by evolutionary models) and ρ_{⋆,step} (constrained by the transit);
In the second case scenario the transit is very shallow, such that b is not directly inferrable from the LC. Besides [Fe/H], 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:
to be added to the merit function expressed by Eq. (4). The quantity ρ_{⋆,step} is computed consistently with M_{⋆,step} and R_{⋆,step}, whose values are driven by the theoretical M_{⋆,isoch} and R_{⋆,isoch} coming out at each step. Thanks to ρ_{⋆,step}, is computed through Kepler’s third law and finally W_{step} is expressed as a function of by inverting Eq. (A.4), such that transit parameters are constrained a posteriori by evolutionary models. We note that if the stellar radius has been also set as a prior as a result of the BP given by Eq. (3), its value also feeds back into ρ_{⋆} and into the corresponding transit parameters in general.
From M_{⋆,step} and R_{⋆,step}, all the relevant planetary parameters can be computed, especially 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 i77820 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 switchedon) took ∼30 h. One single run for analysing Wasp4 system (14 LC time series for a total amount of 4010 data points, one orbiting planet, element diffusion switchedon) took ∼20 h. If we considered a very basic fit of one single planet orbiting a star with just 1 available LC, the expected runtime 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 runtime. 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; ForemanMackey et al. 2013) technique. We are currently working on the implementation of the AIMCMC 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 runtime speed by a factor 10. The full FORTRAN implementation of the AIMCMC 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 planethosting 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 WASP4, 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), SanchisOjeda 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 closeby 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.
HD 219134: input stellar parameters considered in this work.
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 Rapproach; 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 colourmagnitude diagram (left panel) and on the log T_{eff}–M_{V} 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 supersolar metallicity Z = 0.070 → [Fe/H] = 0.66. 
Fig. 2.
Values of T_{eff} found in the literature for HD 219134. Values on the xaxis 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 HARPSN 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 detrending 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 multiplanetary system as HD 219134.
For all four planets, T_{0}, P, and K_{2} 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}, , 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}, , 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_{1} and c_{2} 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 GelmanRubin test), that are summarized in Tables 2 and 3. In particular, the transiting HD 219134 b and c are confirmed to be two superEarth rocky planets. Their corresponding LCs with the superimposed transit models as inferred from the Lapproach 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. 
HD 219134: output stellar parameters derived considering L approach or R approach.
Planets hosted by HD 219134.
A comparison involving radii and masses of HD 219134 and its planets as derived by various authors is presented in Table 4. Gillon et al. (2017) analyse the system assuming stellar 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.
Comparison involving masses and radii derived by different authors for HD 219134 system.
Finally, we comment on the results we obtained with our two different and independent L and R approaches with a particular reference to our derived R_{⋆} values. Input V magnitude and Gaia parallax are totally independent of (and possibly inconsistent with) the interferometric measure of the radius. Actually, the comparison between the L and R approach aims to investigate the mutual consistency of input parameters by looking at the output results.
In Fig. 4 the red histogram is the output R_{⋆} PDF, which has been inferred considering ([Fe/H]_{step}, B − V_{step}, L_{⋆, step}) as inputs for isochrone fitting and assuming a prior on the parallaxbased 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.
WASP4: regression of M_{⋆} vs. (left panel), and of M_{p} vs. (right panel). A strong correlation is expected because ρ_{⋆} is constrained from the transit, and (M_{p}, R_{p}) are related to (M_{⋆}, R_{⋆}), transit parameters and the RV semiamplitude. 
It is important to stress that the red histogram is the isochronalconsistent 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 isochronalconsistent 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 isochronebased 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 parallaxbased 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. WASP4
WASP4 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 WASP4 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_{0}, 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_{1} and c_{2} 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 SanchisOjeda 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 semiamplitude 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 starplanet 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 GelmanRubin 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 WASP4 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 WASP4 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^{6} 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 versus that should be expected from error propagation; similarly, we have versus .
WASP4 planetary system.
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 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 semiamplitude; 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 (resp. ) values. Linear correlations are clear and the coefficients of determination are (left panel) and (right panel).
If a strict linear correlation held (r^{2} = 1), no dispersion in and in would be registered, thus both their uncertainties would be zero. The equation represents the fraction of root mean square that is not explained by the linear regression, thus it is statistically expected that error bar extensions coming from simple error propagation (where the involved quantities are considered independent) may be reduced up to a factor σ_{ne} or so. Statistically, uncertainties may go down to the following estimated levels:
It turned out that the statistically expected uncertainties (indicated by the subscript “exp”) are in agreement with the actual derived ones (indicated by the subscript “out”).
Table 6 reports a comparison involving stellar and planetary masses and radii, as derived by different authors. Our results have the advantage of retrieving M_{⋆} and R_{⋆} owing to the strict interaction of the MCMC algorithm with the stellar evolutionary models, which yields to very precise M_{p} and R_{p} values, if compared to what is stated in the literature.
WASP4: masses and radii of both the star and the planet derived by different authors.
3.3. Comments on the uncertainties of isochronal parameters
Output parameters that are computed from isochrones (i.e. t_{⋆}, M_{⋆} and R_{⋆}) are given as if stellar evolutionary models are perfect. Therefore their respective errors are likely to be underestimated and should be considered as internal. One way to attribute more realistic uncertainties to isochronal parameters is by comparing our results with another set of evolutionary models. If the target stars and the interpolating algorithm are the same, differences arising in the output are to be attributed to the different input physics of the evolutionary models, such as the equations of state, adopted solar mixture, initial chemical composition of stars, nuclear reaction rates, opacities, overshooting treatment, mixinglength 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 WASP4 falls inside. Among the several differences in terms of stellar input physics between PARSEC and CLES models, one is related to the helium content Y. In PARSEC models, Y is assumed to increase with Z, according to
where Y_{p} = 0.2485 (Komatsu et al. 2011) is the primordial helium abundance, while the heliumtometal enrichment ratio is derived from solar calibration. Instead, is not fixed in CLES, such that Y is not a function of Z, and actually for any given value of Z, four possible values of Y are available, from Y = 0.25 to Y = 0.28 at steps of 0.01.
To evaluate the impact of stellar model systematics, we decided to analyse WASP4 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 WASP4. 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 WASP4, 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 WASP4, 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.
NASA Exoplanet Archive: https://exoplanetarchive.ipac.caltech.edu/docs/counts_detail.html
It is just a reference value that considers the effective length of a chain, where 20% burnin length is subtracted to an original chain made of 100 000 steps. Actually Eq. (17) is almost L_{c} independent for L_{c} of order of thousands, which is a common practice in MCMC applications.
The MCMCI code is available at https://github.com/Bonfanti88/MCMCI
Constraining stellar density from transit with a relative uncertainty lower than 2% as in this case is not so common, and it has been possible because of the prominent transit depth and high quality LCs. According to NASA Exoplanet Archive, the median uncertainty for transitinferred ρ_{⋆} is δρ_{⋆, t} ∼ 9%, and only ∼5% of the systems have δρ_{⋆, t} < 2%.
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.SFNRS. The research leading to these results has received funding from the ARC grant for Concerted Research Actions, financed by the WalloniaBrussels 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. CollierCameron for fruitful discussions that help improving the paper as well.
References
 Anderson, D. R., Collier Cameron, A., Hellier, C., et al. 2011, ApJ, 726, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Angus, R., Morton, T. D., ForemanMackey, D., et al. 2019, AJ, 158, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baluev, R. V. 2018, Astron. Comput., 25, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [Google Scholar]
 Barnes, S. A. 2010, ApJ, 722, 222 [NASA ADS] [CrossRef] [Google Scholar]
 Barragán, O., Gandolfi, D., & Antoniciello, G. 2019, MNRAS, 482, 1017 [NASA ADS] [CrossRef] [Google Scholar]
 Basri, G., Borucki, W. J., & Koch, D. 2005, New Astron. Rev., 49, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Beerer, I. M., Knutson, H. A., Burrows, A., et al. 2011, ApJ, 727, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Betancourt, M. 2017, ArXiv eprints [arXiv:1701.02434] [Google Scholar]
 Bonfanti, A., Ortolani, S., Piotto, G., & Nascimbeni, V. 2015, A&A, 575, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonfanti, A., Ortolani, S., & Nascimbeni, V. 2016, A&A, 585, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bouma, L. G., Winn, J. N., Baxter, C., et al. 2019, AJ, 157, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, ApJ, 757, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, Eur. Phys. J. Web Conf., 47, 03005 [CrossRef] [EDP Sciences] [Google Scholar]
 Buldgen, G., Reese, D. R., & Dupret, M. A. 2016, A&A, 585, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Burgers, J. 1969, Flow Equations for Composites Gases (New York: Academic Press) [Google Scholar]
 Caffau, E., Ludwig, H.G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Casella, G., & George, E. I. 1992, Am. Stat., 46, 167 [Google Scholar]
 Chaboyer, B., Fenton, W. H., Nelan, J. E., Patnaude, D. J., & Simon, F. E. 2001, ApJ, 562, 521 [NASA ADS] [CrossRef] [Google Scholar]
 Chaplin, W. J., Basu, S., Huber, D., et al. 2014, ApJS, 210, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Chapman, S., & Cowling, T. G. 1970, The Mathematical Theory of Nonuniform Gases. An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases, 3rd edn. (Cambridge: Cambridge University Press) [Google Scholar]
 Charbonneau, D., Knutson, H. A., Barman, T., et al. 2008, ApJ, 686, 1341 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A., & Bloemen, S. 2011, A&A, 529, A75 [Google Scholar]
 Collier Cameron, A., Bruce, V. A., Miller, G. R. M., Triaud, A. H. M. J., & Queloz, D. 2010, MNRAS, 403, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Cosentino, R., Lovis, C., Pepe, F., et al. 2012, in Groundbased and Airborne Instrumentation for Astronomy IV, Proc. SPIE, 8446, 84461V [Google Scholar]
 da Silva, R., Porto de Mello, G. F., Milone, A. C., et al. 2012, A&A, 542, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Davenport, J. R. A., Hawley, S. L., Hebb, L., et al. 2014, ApJ, 797, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Dawson, R. I., & Johnson, J. A. 2012, ApJ, 756, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Dawson, R. I., & Johnson, J. A. 2018, ARA&A, 56, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Denissenkov, P. A. 2010, ApJ, 719, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Dotter, A. 2016, ApJS, 222, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Dotter, A., Conroy, C., Cargile, P., & Asplund, M. 2017, ApJ, 840, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Eastman, J., Siverd, R., & Scott Gaudi, B. 2010, PASP, 122, 935 [NASA ADS] [CrossRef] [Google Scholar]
 Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Eastman, J. D., Rodriguez, J. E., Agol, E., et al. 2019, PASP, submitted [arXiv:1907.09480] [Google Scholar]
 Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 90, 2262 [Google Scholar]
 Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Feltzing, S., Howes, L. M., McMillan, P. J., & Stonkutė, E. 2017, MNRAS, 465, L109 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B. 2005, AJ, 129, 1706 [Google Scholar]
 Ford, E. B. 2006, ApJ, 642, 505 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
 ForemanMackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [Google Scholar]
 Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457 [Google Scholar]
 Gillon, M., Demory, B.O., Triaud, A. H. M. J., et al. 2009a, A&A, 506, 359 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gillon, M., Smalley, B., Hebb, L., et al. 2009b, A&A, 496, 259 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gillon, M., Lanotte, A. A., Barman, T., et al. 2010, A&A, 511, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gillon, M., Doyle, A. P., Lendl, M., et al. 2011, A&A, 533, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Fortney, J. J., et al. 2012, A&A, 542, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gillon, M., Demory, B.O., Van Grootel, V., et al. 2017, Nat. Astron., 1, 0056 [NASA ADS] [CrossRef] [Google Scholar]
 Giménez, A. 2006, ApJ, 650, 408 [NASA ADS] [CrossRef] [Google Scholar]
 Günther, M. N., & Daylan, T. 2019, Astrophysics Source Code Library [record ascl:1903.003] [Google Scholar]
 Hartman, J. D., Bakos, G. Á., Bayliss, D., et al. 2019, AJ, 157, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Hastings, W. K. 1970, Biometrika, 57, 97 [Google Scholar]
 Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Holman, M. J., Winn, J. N., Latham, D. W., et al. 2006, ApJ, 652, 1715 [NASA ADS] [CrossRef] [Google Scholar]
 Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398 [NASA ADS] [CrossRef] [Google Scholar]
 IglesiasMarzoa, R., LópezMorales, M., & Morales, M. J. A. 2015, PASP, 127, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Jørgensen, B. R., & Lindegren, L. 2005, A&A, 436, 127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., & Megeath, S. T. 2008, ApJ, 673, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Kopal, Z. 1959, Close Binary Systems (London: Chapman & Hall) [Google Scholar]
 Kreidberg, L. 2015, PASP, 127, 1161 [NASA ADS] [CrossRef] [Google Scholar]
 Lebreton, Y., & Goupil, M. J. 2014, A&A, 569, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Loeb, A., & Gaudi, B. S. 2003, ApJ, 588, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Madhusudhan, N., & Winn, J. N. 2009, ApJ, 693, 784 [NASA ADS] [CrossRef] [Google Scholar]
 Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [NASA ADS] [CrossRef] [Google Scholar]
 Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [Google Scholar]
 Meibom, S., Barnes, S. A., Platais, I., et al. 2015, Nature, 517, 589 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Mishenina, T. V., Soubiran, C., Bienaymé, O., et al. 2008, A&A, 489, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morris, S. L. 1985, ApJ, 295, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Motalebi, F., Udry, S., Gillon, M., et al. 2015, A&A, 584, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Murray, C., & Correia, A. 2011, in Exoplanets, ed. S. Seager (Tucson, AZ: University of Arizona Press), 526, 15 [NASA ADS] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Neal, R. M. 2012, Handbook of Markov Chain Monte Carlo [arXiv:1206.1901] [Google Scholar]
 Nikolov, N., Henning, T., Koppenhoefer, J., et al. 2012, A&A, 539, A159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nissen, P. E. 2016, A&A, 593, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pont, F., & Eyer, L. 2004, MNRAS, 351, 487 [NASA ADS] [CrossRef] [Google Scholar]
 Pont, F., Zucker, S., & Queloz, D. 2006, MNRAS, 373, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Pont, F., Husnoo, N., Mazeh, T., & Fabrycky, D. 2011, MNRAS, 414, 1278 [NASA ADS] [CrossRef] [Google Scholar]
 Prugniel, P., Vauglin, I., & Koleva, M. 2011, A&A, 531, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ramírez, I., Allende Prieto, C., & Lambert, D. L. 2013, ApJ, 764, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Rasmussen, C. E., & Williams, C. K. I. 2005, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) (The MIT Press) [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Roberts, G. O., Gelman, A., & Gilks, W. R. 1997, Ann. Appl. Probab., 7, 110 [CrossRef] [Google Scholar]
 Rodríguez, A., & FerrazMello, S. 2010, in EAS Pub. Ser., eds. K. Gożdziewski, A. Niedzielski, & J. Schneider, 42, 411 [CrossRef] [Google Scholar]
 Rodríguez, A., FerrazMello, S., Michtchenko, T. A., Beaugé, C., & Miloni, O. 2011, MNRAS, 415, 2349 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenfield, P., Marigo, P., Girardi, L., et al. 2016, ApJ, 822, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Rouan, D., Baglin, A., Barge, P., et al. 1999, Phys. Chem. Earth C, 24, 567 [Google Scholar]
 Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley) [Google Scholar]
 SanchisOjeda, R., Winn, J. N., Holman, M. J., et al. 2011, ApJ, 733, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
 Scuflaire, R., Théado, S., Montalbán, J., et al. 2008, Ap&SS, 316, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Seager, S., & MallénOrnelas, G. 2003, ApJ, 585, 1038 [NASA ADS] [CrossRef] [Google Scholar]
 Sharma, S., Stello, D., Buder, S., et al. 2018, MNRAS, 473, 2004 [NASA ADS] [CrossRef] [Google Scholar]
 Siverd, R. J., Beatty, T. G., Pepper, J., et al. 2012, ApJ, 761, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Skumanich, A. 1972, ApJ, 171, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Slumstrup, D., Grundahl, F., Brogaard, K., et al. 2017, A&A, 604, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smith, A. M. S., Gandolfi, D., Barragán, O., et al. 2017, MNRAS, 464, 2708 [NASA ADS] [CrossRef] [Google Scholar]
 Soderblom, D. R. 2010, ARA&A, 48, 581 [NASA ADS] [CrossRef] [Google Scholar]
 Soubiran, C., Bienayme, O., Mishenina, T. V., & Kovtyukh, V. V. 2008, VizieR Online Data Catalog: J/A+A/480/91 [Google Scholar]
 Southworth, J., Hinse, T. C., Burgdorf, M. J., et al. 2009, MNRAS, 399, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Sozzetti, A., Torres, G., Charbonneau, D., et al. 2007, ApJ, 664, 1190 [NASA ADS] [CrossRef] [Google Scholar]
 Takeda, G., Ford, E. B., Sills, A., et al. 2007, ApJS, 168, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Ter Braak, C. J. F. 2006, Stat. Comput., 16, 239 [Google Scholar]
 Torres, G., Andersen, J., & Giménez, A. 2010, A&ARv, 18, 67 [Google Scholar]
 Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 van Leeuwen, F. 2007, A&A, 474, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vats, D., & Knudson, C. 2018, ArXiv eprints [arXiv:1812.09384] [Google Scholar]
 Vogt, S. S., Burt, J., Meschiari, S., et al. 2015, ApJ, 814, 12 [NASA ADS] [CrossRef] [Google Scholar]
 White, T. R., Huber, D., Mann, A. W., et al. 2018, MNRAS, 477, 4403 [NASA ADS] [CrossRef] [Google Scholar]
 Wilson, D. M., Gillon, M., Hellier, C., et al. 2008, ApJ, 675, L113 [NASA ADS] [CrossRef] [Google Scholar]
 Winn, J. N. 2010, ArXiv eprints [arXiv:1001.2010] [Google Scholar]
 Winn, J. N., Holman, M. J., Carter, J. A., et al. 2009, AJ, 137, 3826 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, J. T. 2018, Radial Velocities as an Exoplanet Discovery Method, 4 [Google Scholar]
 Wright, J. T., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2004, ApJS, 152, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Yi, S., Demarque, P., Kim, Y.C., et al. 2001, ApJS, 136, 417 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Relation between transit duration W and the scaled semimajor axis
For the transit analysis our code computes the Cartesian coordinates of the Keplerian orbit and then it feeds the Mandel & Agol (2002) model with the planet to star skyprojected distance and the size ratio. Those Cartesian coordinates are computed using notably ; relying on transit observables, comes from P, W, b, and dF (besides e and ω if the orbit is eccentric).
Following Winn (2010), establishing the total transit duration W, that is the temporal interval between the first and fourth contact requires us to compute the integral
where r(f) is the ellipse equation (i.e. the starplanet distance) as a function of the true anomaly f
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 starplanet distance in the sky plane, that is
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– relation that has been implemented in our algorithm is written as
where
Equation (A.4) is a revised version of Eq. (8) by Seager & MallénOrnelas (2003), where we added the eccentricitydependent factors (A.5) and (A.6), as suggested by Winn (2010), and we considered the following approximations:
We evaluated the impact of our approximations, also thanks to some reference data of interest that we inferred from the NASA Exoplanet archive. It turned out that values which are computed through (A.4) differ from those derived from the actual formula of Seager & MallénOrnelas (2003) by less than 1% if , which holds for more than 97% of the exoplanets that have been confirmed so far. It is even more reassuring that the median value of produces negligible variations on (∼0.05%), while in the rare cases in which , we may register variations on of order of a few percent, which are still not worrying in general since the median relative uncertainty on is ∼10%.
Without paying a high price, the introduced approximations enable us to easily express W as a function of . This is useful in case the transit is very shallow and it is preferable not to rely on it to infer W. As also explained in Sect. 2.3, in this scenario ρ_{⋆} is established from evolutionary models, then is computed from ρ_{⋆} thanks to Kepler’s third law, and finally also W is available by easily inverting Eq. (A.4).
All Tables
HD 219134: output stellar parameters derived considering L approach or R approach.
Comparison involving masses and radii derived by different authors for HD 219134 system.
WASP4: masses and radii of both the star and the planet derived by different authors.
All Figures
Fig. 1.
HD 219134 location on the colourmagnitude diagram (left panel) and on the log T_{eff}–M_{V} 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 supersolar metallicity Z = 0.070 → [Fe/H] = 0.66. 

In the text 
Fig. 2.
Values of T_{eff} found in the literature for HD 219134. Values on the xaxis 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.
WASP4: regression of M_{⋆} vs. (left panel), and of M_{p} vs. (right panel). A strong correlation is expected because ρ_{⋆} is constrained from the transit, and (M_{p}, R_{p}) are related to (M_{⋆}, R_{⋆}), transit parameters and the RV semiamplitude. 

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

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.