EDP Sciences
Home All issues Volume 595 (November 2016) A&A, 595 (2016) L5 Full HTML
Free Access
Issue
A&A
Volume 595, November 2016
Article Number L5
Number of page(s) 12
Section Letters
DOI https://doi.org/10.1051/0004-6361/201629770
Published online 28 October 2016

© ESO, 2016

1. Introduction

The WASP-47 planetary system is composed of at least four planets. The first planet to be discovered was the transiting hot Jupiter WASP-47 b in a 4.16 days orbit (Hellier et al. 2012), whose mass was determined by radial velocity measurements. Long-term radial velocity follow-up detected a second massive long-period (~570 days) planet, WASP-47 c (Neveu-VanMalle et al. 2016), which is not known to transit the star. Observations from the K2 mission (Howell et al. 2014) that lasted 69 days allowed detecting two additional transiting planets: a Neptune-sized planet on a 9-day orbit (planet d), and a 1.8 RE planet on a 0.8-day orbit (planet e) (Becker et al. 2015). WASP-47 is the only system known to date with a close-in planet inside the orbit of a hot Jupiter and an external massive companion; a benchmark for planetary formation and migration studies.

The transit-timing variations (TTV) detected in the K2 light curve enabled Becker et al. (2015) to determine the mass of planet d, and to set an upper limit to the mass of planet e. In addition, they produced an independent measurement of the mass of planet b that agrees with the mass from radial velocity. Dai et al. (2015) monitored the system with high-precision radial velocity and detected the Doppler signature of planet e. They also obtained a marginal detection of the mass of planet d, and newly determined the mass of planet b. In all cases, the results of Dai et al. (2015) agree with the TTV analysis of Becker et al. (2015). Sanchis-Ojeda et al. (2015) detected the Rossiter-Mclaughlin effect of planet b, obtaining a low projected stellar obliquity that is compatible with an aligned orbit.

In this Letter we perform the first joint analysis of the available photometric light curves and radial velocity measurements, modelling self-consistently the gravitational interactions between the planets. In line with the ideas introduced in Almenara et al. (2015), we continue to explore the ability of this type of modelling to provide absolute measurements of masses, radii, and densities in multi-planetary systems. Additionally, the self-consistent modelling of WASP-47 improves the determination of the system parameters when we add stellar models.

2. Observations

2.1. Light curves

We employed the K2 light curve reduced by Becker et al. (2015), obtained during Campaign 3, with a time-span of 69 days, a time cadence of 1 min, and 350 ppm precision per data-point. This is the data set that contains most of the information on the dynamical interaction between the planets. In addition, we included a high-precision transit of WASP-47 b obtained with EulerCam (Lendl et al. 2012) that was reported in the original discovery paper by Hellier et al. (2012). The EulerCam light curve was observed using the Gunn r filter and has a temporal sampling of ~77 s. We transformed the EulerCam BJDUTC times to BJDTDB (Eastman et al. 2010) to match the time reference of the K2 data. Only the data spanning three transit durations around each transit mid-time were considered for the analysis, after normalisation with a parabola.

No companion was detected with lucky imaging (Wöllert et al. 2015), reducing the probability of unaccounted for additional contamination in the photometric aperture.

2.2. Radial velocities

We used the published data (Hellier et al. 2012; Neveu-VanMalle et al. 2016; Dai et al. 2015) from the CORALIE (Udry et al. 2000) and PFS (Crane et al. 2010) spectrographs. The precision is between 7 and 27 m s-1 for the CORALIE data, and between 2.5 and 4.0 m s-1 for the PFS data. They form a set of 77 radial velocity measurements that we treated as originating in four different instruments, with independent zero-point offsets: the CORALIE measurements presented by Hellier et al. (2012) (CORALIE 2012), those presented by Neveu-VanMalle et al. (2016) and performed before its upgrade in 2014 (CORALIE 2016a), those performed afterwards (CORALIE 2016b), and the PFS measurements from Dai et al. (2015).

3. Analysis

The light curve and radial velocity model accounts for the gravitational interactions between all components of the system. The model is described in detail in Almenara et al. (2015). Briefly, the projected velocity of the star is obtained as a direct output of the the N-body integrator. The light-curve model is obtained using the analytical formulae of Mandel & Agol (2002) with a quadratic limb-darkening law (Manduca et al. 1977). The projected centre-to-centre distance between the star and the planets is obtained from the N-body integration. In this analysis, we used REBOUND (Rein & Liu 2012) with the WHFast integrator (Rein & Tamayo 2015). The integration time step was set to 0.001 days, which results in a maximum error of 8  cm s-1 for the radial velocity model. The combined error in the photometric model is smaller than 6 ppm per observed data point because of the discrete N-body time step and the oversampling factor of 10 (Kipping 2010).

The N-body integrator was initialized with the positions and velocities of the system bodies (parametrised by osculating astrocentric asteroidal orbital elements, see Table A.1), at a given reference time, tref = 2456979.5 BJDTDB. The reference time tref was chosen immediately before the first K2 transit, which is close to the moment containing most of the dynamical information. This means that we used two N-body integrations initialized at tref, one that was integrated forward in time until the last PFS data point in 2015, and a second one that was integrated backwards in time to the first observation of CORALIE in 2010.

The physical parameters of the model are the bodies’ masses and radii, the orbital parameters at tref, and the limb-darkening coefficients for each photometric instrument. To minimise correlations, we parametrised the model using the combinations listed in Table A.1. Additionally, we fitted a normalisation factor for each light curve, an offset value for each radial velocity data set (we assumed zero systemic velocity), and a multiplicative factor to account for additional noise not included in the reported uncertainties (the jitter parameter) for each data set. Although we did not use stellar models, our model implicitly assumes spherical shapes for the star and planets. Under this assumption the model does not depend on the absolute value of the longitude of the ascending node (Ω). Therefore, we fixed Ωb at tref to 180°. Finally, there is also a reflection symmetry that allows choosing the hemisphere on which one of the planets transits. We therefore constrained the orbital inclination of planet b, ib< 90° (Fig. 1).

thumbnail Fig. 1

Schematic view of the WASP-47 system. The sizes of the star and the planets are to scale, and the distances from the surface of the star are in logarithmic scale. The vertical positions of the filled circles representing the planets correspond to maximum a posteriori estimates of the impact parameters, and the shaded areas (offset in the x direction for clarity) represent the corresponding posterior distributions. The distribution of the planet c is much larger and is truncated in this figure; its size is unknown.
Open with DEXTER

We coupled the dynamical model with the emcee algorithm (Goodman & Weare 2010; Foreman-Mackey et al. 2013) to sample from the posterior distribution of the parameters. The orbital inclination of planet d exhibits a bimodal distribution with well-separated modes, each corresponding to the planet transiting a different stellar hemisphere. As emcee does not sample well-separated multi-modal distributions efficiently (Foreman-Mackey et al. 2013), we ran two different emcee runs, one for each mode, and combined the results as explained below. The dynamical interactions are different between these orbital configurations, and it is therefore important to consider them separately, as they could in principle be distinguished. We note that most TTV analyses assume coplanar orbits, as was the case in Becker et al. (2015). Concerning the remaining planets, their inclination distributions are either very wide and do not exhibit a bimodal distribution (planet c), or the modes are not separate (planet e). The MCMC algorithm should therefore be able to sample from them correctly. As explained above, the inclination of planet b was limited to ib< 90°, preventing its bi-modality.

We used non-informative uniform priors for all 48 emcee jump parameters. We ran emcee with 250 walkers, from a starting point constructed as follows: first, we modelled the radial velocity alone using Keplerian orbits with zero eccentricity for all planets except for planet c; second, we included the dynamical interactions between the planets, allowed for eccentric orbits, and also modelled the light-curve data. The transit times from Becker et al. (2015) were included as an ancillary data set to reduce the algorithm burn-in phase. Then, the final emcee state was used as starting point for the final modelling, for which the transit times were removed. This procedure reduces the burn-in phase in the final photodynamical modelling. For the last step we ran 215000 steps of the emcee algorithm and considered the first 190000 steps as an additional burn-in period, which was removed for the final inference. The results of the two different orbital configurations described above were combined assuming equal probability for each mode.

4. Results

Table A.1 lists the maximum a posteriori (MAP) estimate, the median, the 68% confidence interval (CI), and the 95% highest density interval (HDI) of the inferred system parameters. Figures A.1 and A.2 show the radial velocity measurements and the transit light curves, together with their respective models, Fig. A.3 show the posterior of the TTV. The posterior distributions coming from the two orbital configurations for planet d are very similar (Fig. A.4). The analysis seems largely insensitive to the hemisphere on which planet d transits, which justifies combining samples from both modes assuming equal probability for each.

The bulk densities are determined with a precision between 1.5% for the stellar host and 38% for 1.8-RE planet e. The precision in the stellar density is similar to what can be obtained through asteroseismology (Huber et al. 2013). The density of planet d is 1.8 ± 0.16 gcm-3, which represents a precision of 9% and constitutes an improvement of almost a factor of 10 with respect to the analysis of Dai et al. (2015). The radii of the star and of all transiting planets are determined with a precision of 22%. As all planet-to-star radius ratios are very precisely constrained by the transit observations, the knowledge on the absolute radii is dominated by the determination of the stellar radius. Concerning the masses, their posterior distributions are skewed, and it is therefore not trivial to quote a precision value. Using the half-width of the 68.3%-HDI and its central value, we obtain a precision of between 50% and 60% for the star and all planets, except for planet b, whose absolute mass is known to a precision of around 40%. We note that this precision is a factor of three lower than the one obtained for the radii. This is surely linked to the fact that as a result of the nature of the gravitational force, the best constrained parameters are the bulk densities of the objects, which involves the third power of the radii but depends only linearly on the masses. We note that our model enables inferring all these physical parameters without resorting to theoretical stellar models, and their accuracy depends, therefore, only on the validity of the model assumptions.

Concerning the orbital eccentricities, Becker et al. (2015) obtained upper limits using long-term dynamical simulations and concluded that the orbits would be unstable for eccentricities above 0.06. This result was later used as a prior for their TTV analysis. Dai et al. (2015) assumed circular orbits for the three interior planets. Here, we obtain stringent eccentricity upper limits directly from the analysis of the data, without requiring additional assumptions about the long-term stability of the system. The eccentricity 99% upper confidence limit for planets e, b, and d are 0.17, 0.014, and 0.033, respectively.

The dynamical interactions between the planets allow constraining the true mass of the non-transiting planet c. The inclination is constrained to be higher than 18° at 99%, with a corresponding mass upper limit of 4.7 MJ. On the other hand, the longitude of the ascending node for this planet is completely unconstrained. The model permits even retrograde orbits, which would not seem to produce a detectable effect on the timescale of the observations. Long-term stability arguments should probably distinguish between prograde and retrograde orbits.

The limb-darkening parameters are precisely constrained thanks to the presence of a giant transiting planet, which produces high signal-to-noise ratio transits. This produces precisely determined radius ratios, which are usually covariate with the limb-darkening coefficients. This is particularly important for the smaller planets in the system, whose transits do not provide strong information on the stellar flux distribution across the disk.

thumbnail Fig. 2

Mass-radius, radius-density, and mass-density diagrams for known exoplanets with parameters determined with a precision <50% (Han et al. 2014). The WASP-47 planets are shown as 68.3% coloured contours for the model-free photodynamical analysis and as dots with 1σ error bars for the analysis using the Dartmouth model (planet b in green, planet c in orange, planet d in blue, and planet e in red). Planet c, whose radius and density are not known, is represented at the top of the mass-radius and mass-density diagrams. The solid lines represent theoretical models for different composition (Zeng & Sasselov 2013).
Open with DEXTER

4.1. Stellar models

To improve on the precision of our results, at the expense of accuracy, we added the constraints coming from the theoretical stellar models. To do this, we kept only the posterior samples that are compatible with the Dartmouth models (Dotter et al. 2008), given the bulk density from our analysis and the stellar atmospheric parameters from Mortier et al. (2013): Teff = 5576 ± 68 K, [Fe/H] = 0.36 ± 0.05. This significantly narrowed the posterior distributions of the stellar mass and radius, and, therefore, those of the planetary parameters. Curiously, the MAP estimate of the stellar mass agrees with the value obtained using stellar models. In Table A.1 we report in the last column the parameters for which the precision was improved significantly with respect to our model-free analysis.

In addition to improving the precision of the masses, radii, and semi-major axes for all the planets, the inclusion of theoretical stellar models also narrowed the posterior distributions of other parameters. Radial velocities set stringent constraints on the properties of planet b, but are much less efficient in determining the mass of planet d, for which the dynamical interaction with the remaining planets is a major asset. As a consequence, including stellar models improved the determination of the radial velocity amplitude induced by planet d, which depends on the ratio , but did not affect the determination of the mass ratio Md/M⋆ significantly. Conversely, the mass ratio Mb/M⋆ was improved by the inclusion of the models but not the radial velocity amplitude. As a corollary, given that the mean stellar density and the radius ratios were not affected by the inclusion of the theoretical models, the bulk density of planet b and the surface gravity of planet d were improved as well1. This allows distinguishing whether the mass detection for a given planet is dominated by the radial velocity or timing measurements.

Concerning the stellar age, the Dartmouth models provided an age of of 7.1 ± 1.5 Gyr. Using the projected velocity  km s-1 from Sanchis-Ojeda et al. (2015) and the radius from the Dartmouth tracks and assuming the stellar rotational axis lies in the plane of the sky, i⋆ = 90°, we derived a rotational period of Prot = 31.2 ± 3.5 days. Coupling this with the mass determination, we derived a gyrochronological age of 6.5 Gyr (Barnes 2010; Barnes & Kim 2010; Meibom et al. 2015). Both age determinations agree well, but they disagree with a previous gyrochronological estimation by Brown (2014) and Hellier et al. (2012), who used a slightly higher value of v sin i⋆ = 3.0 ± 0.6.

The next three transit windows of planet c are 2457741 ± 22, 2458327 ± 30, and 2458914 ± 35 BJD. The expected centre transit duration for a Jupiter-sized planet is 14.0 ± 1.6 h. If the transit windows are reduced with follow-up radial velocity observations, this could be a target for the upcoming CHEOPS mission (Broeg et al. 2013).

5. Discussion

Our analysis allowed characterising the planets and the star of the WASP-47 system without resorting to theoretical stellar models at all, except for the implicit assumption of spherical bodies. The self-consistent modelling is mandatory to optimally exploit the data of this system and led to an improved determination for the planetary bulk densities for some of the planets with respect to previous studies. In addition, these measurements are probably more accurate as well.

In Fig. 2 we placed the WASP-47 planets in mass-radius, mass-density, and radius-density diagrams, together with the planets reported in the literature (Han et al. 2014). Our analysis confirms that the structure of planet d must be dominated by volatiles, and the rocky character of planet e. Although Rogers (2015) determined that most planets with radii larger that 1.6 RE are not dense enough to be rocky, the model predicts a non-zero fraction of rocky planets as large as WASP-47 e.

The habitable zone (HZ) of WASP-47 extends between 1.06 au and 1.8 au (Kopparapu et al. 2013). The MAP estimate of the semi-major axis of the orbit of WASP-47 c is 1.37 au, but owing to the orbital eccentricity, the planet makes excursions outside of the HZ (Fig. A.5). However, Williams & Pollard (2002) argued that long-term climate stability is dictated by the mean incident flux throughout the orbit rather than by the time spent in the HZ. The effective incident flux on WASP-47 c is 0.640 ± 0.052 times the flux received by Earth, which would place it in the middle of the HZ (see Díaz et al. 2016; Kopparapu et al. 2013, Fig. 8). Therefore, hypothetical rocky satellites orbiting WASP-47 c have good prospects for being habitable.

Including theoretical stellar models improves the precision for some parameters, in particular masses and radii, but does not change our conclusions about the planet compositions significantly. This is probably because WASP-47 is a solar-type star for which models are expected to perform correctly. On the other hand, these determinations are probably less accurate, and might exhibit systematic errors larger than the measured precision.

The current precision of the photometric measurements of extrasolar planets already calls for the use of realistic models including the interactions between the planets. In the coming years, TESS (Ricker et al. 2014) and PLATO (Rauer et al. 2014) will continue to provide such high-quality light curves, but the brightness of the typical target will also permit obtaining high-precision radial velocity measurements, further increasing the need for realistic self-consistent modelling of the data.

1

The bulk density is ρp = ρ⋆(Mp/M⋆)(Rp/R⋆)-3, while the surface gravity is .

Acknowledgments

We thank C. Hellier, who kindly sent us the EulerCam transit data, R. Mardling for discussions about dynamics, T. Fenouillet and Y. Revaz for their assistance with the computing clusters used in this work, and L. Kreidberg for her Mandel & Agol code. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. Simulations in this paper made use of the REBOUND code which can be downloaded freely at http://github.com/hannorein/rebound. Part of these simulations have been run on the Regor cluster kindly provided by the Observatoire de Genève. J.M.A. and X.B. acknowledges funding from the European Research Council under the ERC Grant Agreement No. 337591-ExTrA. This work has been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF).

References

Appendix A: Additional figures and tables

Table A.1

Inferred system parameters.

thumbnail Fig. A.1

Photometric data of planets b) (top, left), d) (top, right), and e) (bottom). Each transit is centred relative to a linear ephemeris. Small dots represents the photometric data, whereas circles represent the data binned. Each panel is labelled with the epoch, with zero the first transit after tref. The black curve is the median value of the distribution of oversampled models corresponding to 1000 random MCMC steps, the different shades of grey represent the 68.3, 95.5, and 99.7% confidence intervals. In the lower part of each panel the residuals after subtracting the MAP model to the observed data are shown.
Open with DEXTER

thumbnail Fig. A.2

Idem Fig. A.1 for the radial velocity as a function of time. From top to bottom: CORALIE 2012, CORALIE 2016a, CORALIE 2016b, and PFS.
Open with DEXTER

thumbnail Fig. A.3

Posterior TTV of planets e), b), and d) (from top to bottom). Black dots with coloured error bars are the posterior TTV from the analysis (median and 68.3% confidence interval). The grey open circles are the Becker et al. (2015) TTV, measured individually from each observed transit. The gain in the mean transit times precision is a factor 9 for planets b) and e), and 3 for planet d).
Open with DEXTER

thumbnail Fig. A.4

Two-parameter joint posterior distributions for the most relevant MCMC model parameters. The 39.3, 86.5, and 98.9% two-variable joint confidence regions (in the case of a Gaussian posterior, these regions project on to the one-dimensional 1, 2, and 3σ intervals) are denoted by three different grey levels. The histogram of each parameter is shown at the top of each column, except for the parameter on the last line, which is shown at the end of the line. Units are the same as in Table A.1.
Open with DEXTER

thumbnail Fig. A.5

Schematic view of the orbital plane of WASP-47 c. The host star is represented by the orange circle at the centre. The maximum a posteriori orbit is indicated by the empty black points, which are equally spaced in time over the orbit. The orbital movement is counter-clockwise. The angles are measured counter-clockwise from the negative x-axis. The semi-major axis of the orbit is shown as a thin grey line, and the concentric circles are labelled in astronomical units. The black thick arrow points towards the observer. The filled green area is the habitable zone comprised between the runaway greenhouse limit and the maximum greenhouse limit, according to the model of Kopparapu et al. (2013). The red area corresponds to the increased habitable zone if the outer edge is defined by the early-Mars limit.
Open with DEXTER

All Tables

Table A.1

Inferred system parameters.

In the text

All Figures

thumbnail Fig. 1

Schematic view of the WASP-47 system. The sizes of the star and the planets are to scale, and the distances from the surface of the star are in logarithmic scale. The vertical positions of the filled circles representing the planets correspond to maximum a posteriori estimates of the impact parameters, and the shaded areas (offset in the x direction for clarity) represent the corresponding posterior distributions. The distribution of the planet c is much larger and is truncated in this figure; its size is unknown.
Open with DEXTER
In the text
thumbnail Fig. 2

Mass-radius, radius-density, and mass-density diagrams for known exoplanets with parameters determined with a precision <50% (Han et al. 2014). The WASP-47 planets are shown as 68.3% coloured contours for the model-free photodynamical analysis and as dots with 1σ error bars for the analysis using the Dartmouth model (planet b in green, planet c in orange, planet d in blue, and planet e in red). Planet c, whose radius and density are not known, is represented at the top of the mass-radius and mass-density diagrams. The solid lines represent theoretical models for different composition (Zeng & Sasselov 2013).
Open with DEXTER
In the text
thumbnail Fig. A.1

Photometric data of planets b) (top, left), d) (top, right), and e) (bottom). Each transit is centred relative to a linear ephemeris. Small dots represents the photometric data, whereas circles represent the data binned. Each panel is labelled with the epoch, with zero the first transit after tref. The black curve is the median value of the distribution of oversampled models corresponding to 1000 random MCMC steps, the different shades of grey represent the 68.3, 95.5, and 99.7% confidence intervals. In the lower part of each panel the residuals after subtracting the MAP model to the observed data are shown.
Open with DEXTER
In the text
thumbnail Fig. A.2

Idem Fig. A.1 for the radial velocity as a function of time. From top to bottom: CORALIE 2012, CORALIE 2016a, CORALIE 2016b, and PFS.
Open with DEXTER
In the text
thumbnail Fig. A.3

Posterior TTV of planets e), b), and d) (from top to bottom). Black dots with coloured error bars are the posterior TTV from the analysis (median and 68.3% confidence interval). The grey open circles are the Becker et al. (2015) TTV, measured individually from each observed transit. The gain in the mean transit times precision is a factor 9 for planets b) and e), and 3 for planet d).
Open with DEXTER
In the text
thumbnail Fig. A.4

Two-parameter joint posterior distributions for the most relevant MCMC model parameters. The 39.3, 86.5, and 98.9% two-variable joint confidence regions (in the case of a Gaussian posterior, these regions project on to the one-dimensional 1, 2, and 3σ intervals) are denoted by three different grey levels. The histogram of each parameter is shown at the top of each column, except for the parameter on the last line, which is shown at the end of the line. Units are the same as in Table A.1.
Open with DEXTER
In the text
thumbnail Fig. A.5

Schematic view of the orbital plane of WASP-47 c. The host star is represented by the orange circle at the centre. The maximum a posteriori orbit is indicated by the empty black points, which are equally spaced in time over the orbit. The orbital movement is counter-clockwise. The angles are measured counter-clockwise from the negative x-axis. The semi-major axis of the orbit is shown as a thin grey line, and the concentric circles are labelled in astronomical units. The black thick arrow points towards the observer. The filled green area is the habitable zone comprised between the runaway greenhouse limit and the maximum greenhouse limit, according to the model of Kopparapu et al. (2013). The red area corresponds to the increased habitable zone if the outer edge is defined by the early-Mars limit.
Open with DEXTER
In the text

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

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

Initial download of the metrics may take a while.