The Kepler-19 system: a thick-envelope super-Earth with two Neptune-mass companions characterized using Radial Velocities and Transit Timing Variations
Luca Malavolta, Luca Borsato, Valentina Granata, Giampaolo Piotto, Eric Lopez, Andrew Vanderburg, Pedro Figueira, Annelies Mortier, Valerio Nascimbeni, Laura Affer, Aldo S. Bonomo, Francois Bouchy, Lars A. Buchhave, David Charbonneau, Andrew Collier Cameron, Rosario Cosentino, Courtney D. Dressing, Xavier Dumusque, Aldo F. M. Fiorenzano, Avet Harutyunyan, Raphaëlle D. Haywood, John Asher Johnson, David W. Latham, Mercedes Lopez-Morales, Christophe Lovis, Michel Mayor, Giusi Micela, Emilio Molinari, Fatemeh Motalebi, Francesco Pepe, David F. Phillips, Don Pollacco, Didier Queloz, Ken Rice, Dimitar Sasselov, Damien Ségransan, Alessandro Sozzetti, Stéphane Udry, Chris Watson
DD raft version S eptember
11, 2018
Preprint typeset using L A TEX style AASTeX6 v. 1.0
THE KEPLER-19 SYSTEM: A THICK-ENVELOPE SUPER-EARTH WITH TWO NEPTUNE-MASS COMPANIONSCHARACTERIZED USING RADIAL VELOCITIES AND TRANSIT TIMING VARIATIONS L uca M alavolta , L uca B orsato , V alentina G ranata , G iampaolo P iotto , E ric L opez , A ndrew V anderburg , P edro F igueira ,A nnelies M ortier , V alerio N ascimbeni , L aura A ffer , A ldo S. B onomo , F rancois B ouchy , L ars A. B uchhave , D avid C harbonneau ,A ndrew C ollier C ameron , R osario C osentino , C ourtney D. D ressing , X avier D umusque , A ldo F. M. F iorenzano , A vet H arutyunyan , R apha ¨ elle D. H aywood , J ohn A sher J ohnson , D avid W. L atham , M ercedes L opez -M orales , C hristophe L ovis , M ichel M ayor , G iusi M icela , E milio M olinari , F atemeh M otalebi , F rancesco P epe , D avid F. P hillips , D on P ollacco , D idier Q ueloz ,K en R ice , D imitar S asselov , D amien S´ egransan , A lessandro S ozzetti , S t ´ ephane U dry , C hris W atson Dipartimento di Fisica e Astronomia “Galileo Galilei”, Universita’di Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy; [email protected] INAF - Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy SUPA, Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH93HJ, UK Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, Massachusetts 02138, USA Instituto de Astrof´ısica e Ciˆencias do Espac¸o, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal Centre for Exoplanet Science, SUPA, School of Physics and Astronomy, University of St. Andrews, St. Andrews KY16 9SS, UK INAF - Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90124 Palermo, Italy INAF - Osservatorio Astrofisico di Torino, via Osservatorio 20, 10025 Pino Torinese, Italy Observatoire Astronomique de l’Universit´e de Gen`eve, 51 ch. des Maillettes, 1290 Versoix, Switzerland Centre for Star and Planet Formation, Natural History Museum of Denmark & Niels Bohr Institute, University of Copenhagen, Øster Voldgade 5-7, DK-1350Copenhagen K, Denmark INAF - Fundaci´on Galileo Galilei, Rambla Jos´e Ana Fernandez P´erez 7, 38712 Bre˜na Baja, Spain Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA NASA Sagan Fellow INAF - IASF Milano, via Bassini 15, 20133, Milano, Italy Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK Cavendish Laboratory, J J Thomson Avenue, Cambridge CB3 0HE, UK Astrophysics Research Centre, School of Mathematics and Physics, Queens University Belfast, Belfast BT7 1NN, UK
ABSTRACTWe report a detailed characterization of the
Kepler -19 system. This star was previously known to host atransiting planet with a period of 9.29 days, a radius of 2.2 R ⊕ and an upper limit on the mass of 20 M ⊕ . Thepresence of a second, non-transiting planet was inferred from the transit time variations (TTVs) of Kepler -19bover 8 quarters of
Kepler photometry, although neither mass nor period could be determined. By combiningnew TTVs measurements from all the
Kepler quarters and 91 high-precision radial velocities obtained with theHARPS-N spectrograph, we measured through dynamical simulations a mass of 8 . ± . ⊕ for Kepler -19b.From the same data, assuming system coplanarity, we determined an orbital period of 28.7 days and a massof 13 . ± . ⊕ for Kepler -19c and discovered a Neptune-like planet with a mass of 20 . ± . ⊕ on a63 days orbit. By comparing dynamical simulations with non-interacting Keplerian orbits, we concluded thatneglecting interactions between planets may lead to systematic errors that could hamper the precision in theorbital parameters when the dataset spans several years.With a density of 4 . ± .
87 g cm − (0 . ± . ρ ⊕ ) Kepler -19b belongs to the group of planets with a rockycore and a significant fraction of volatiles, in opposition to low-density planets characterized by transit-timevariations only and the increasing number of rocky planets with Earth-like density.
Kepler -19 joins the smallnumber of systems that reconcile transit timing variation and radial velocity measurements. INTRODUCTIONAfter the discovery of thousands of planets with radiismaller than 2.7 R ⊕ from the NASA Kepler mission (Boruckiet al. 2011; Coughlin et al. 2016), a consistent e ff ort has been devoted to understanding the formation scenario and chemi-cal composition of such planets (e. g., Weiss & Marcy 2014,Dressing et al. 2015 and Wolfgang & Lopez 2015). To dis-tinguish between a rocky composition and the presence ofa thick envelope of water or volatile elements, the radius a r X i v : . [ a s t r o - ph . E P ] M a r M alavolta et al .derived from the transit depth must be coupled with a pre-cise mass determination (better than 20%), either from ra-dial velocity (RV) measurements or transit timing variations(TTVs). Planets that have been characterized with such levelof precision appear to fall into two populations, a first onefollowing an Earth-like composition and a second one withplanets larger than 2 Earth radii requiring a significant frac-tion of volatiles (e. g., Rogers 2015, Gettel et al. 2016 andL´opez-Morales et al. 2016). Recently, improved mass andradius determinations of known planets have uncovered theexistence of super-Earths which fall between these two popu-lations, such as 55 Cancri e (Demory et al. 2016) and Kepler -20b (Buchhave et al. 2016).We need more planets in that mass regime with precisemass and radius measurements to understand the undergoingphysics. For this reason we carried out a RV follow-up of
Kepler -19b (hereafter K . . ± .
048 R ⊕ , orbiting a relatively bright( V = . , K = .
9) solar-type star ( T e ff = ±
60 K,log g = . ± .
10, [Fe / H] = − . ± . Kepler light curve,
Spitzer observations to verify the achro-maticity of the transit and Keck-HIRES high-resolution spec-troscopy to rule out the presence of massive, non-planetaryperturbers. From
BLENDER analysis (Torres et al. 2011) theprobability of a false-positive scenario was constrained toless than 1 . × − . RVs measured on high-resolution spec-tra were consistent with a mass of 1.5 M ⊕ (0.5 m s − ) andactivity-induced RV jitter of 4 m s − , but ultimately theylacked the required precision for a robust determination ofthe planetary mass, and only an upper limit of 20.3 M ⊕ wasset. The existence of an additional planet, Kepler -19c (here-after K (cid:46)
160 days and mass (cid:46) jup , anda further confirmation of the planetary nature of K
19b wereinferred by B11 from the presence of TTVs on 8 quarters of K
19b light curve.In this paper we couple high-precision RV measurementsobtained with HARPS-N with updated measurements of tran-sit times ( T ) encompassing all 17 quarters of Kepler data todetermine the orbital parameters of K K
19c and a third,previously unknown planet in the system,
Kepler -19d (here-after K ff ectof di ff erent mutual inclinations on the goodness of the fit.We confirm that only the inner planet is seen transiting thehost star. A comparison between RV obtained from dynam-ical simulations and when assuming non-interacting planets is performed. We conclude describing the role of K
19b inunderstanding the bulk densities of small planets. RADIAL VELOCITIESWe collected 101 spectra using HARPS-N at the Telesco-pio Nazionale Galileo, in La Palma. Observations spannedover two years, from June 2012 to November 2014, over-lapping the
Kepler observations during the first year. Ev-ery observation consisted of a 30 minute exposure, with amedian signal-to-noise ratio of 37 at 550 nm, correspond-ing to a RV nominal error of 2.8 m s − . Given the faint-ness of the target, observations were gathered with the ob-jAB setting, i. e., the second fiber (fiber B) was observingthe sky instead of acquiring a simultaneous thorium-argon(ThAr) lamp spectrum. Several observations demonstratedthat the stability of the instrument over 24 hours is within 1m s − (e. g., Cosentino et al. 2014), thus the precision of themeasurements was dictated largely by photon noise. Datawere reduced using the standard Data Reduction Software(DRS) using a G8 flux template (the closest available oneto the spectral type of the target) to correct for variations inthe flux distribution as a function of the wavelength, and aG2 binary mask was used to perform the cross-correlation(Baranne et al. 1996; Pepe et al. 2002). The resulting RVdata with their formal 1- σ uncertainties, the Full Width HalfMaximum (FWHM) of the cross-correlation function (CCF)and its contrast (i. e., the depth normalized to the continuum),the bisector inverse span (BIS), and the log R (cid:48) HK activity in-dex are listed in Table 1.2.1. E ff ect of Moon Illumination A simple procedure was adopted to check the influenceof the moon illumination on the science fiber (labeled asfiber A). First, the cross-correlation function of the sky spec-trum acquired with fiber B, CCF B , was recomputed usingthe same flux correction coe ffi cients as for the target (CCF A )for that specific acquisition. Then, CCF B was subtracted tothe corresponding CCF A and radial velocities were computedagain using the script from Figueira et al. (2013), which usesthe same algorithm implemented in the DRS. For 10 obser-vations the di ff erence between the sky-corrected RV and theDRS RV was greater than twice the photon noise, so we re-jected those observations, and used the remaining 91 RVsfrom the DRS in the following analysis. A flag has been in-cluded in Table 1 to identify the rejected observations.While the rejected observations have in common a frac-tion of the illuminated Moon greater than 0.9 and a barycen-tric RV correction within 15 km s − from the absolute RV ofthe target star, not all the observations that met this criteriumwere a ff ected by the sky contamination, suggesting that other When using the objAB setting, CCF B is computed by the DRS withoutflux correction by default. he Kepler -19 system Table 1 . HARPS-N radial velocities and ancillary measurements of Kepler-19. Epoch of the observations, RV with associated noise, FWHMand contrast of the CCF, inverse bisector span,log R (cid:48) HK activity index with the associated error, and a flag indicating if the data has beencontaminated by the Moon (1) or not (0), see Section 2.1.BJD UTC RV σ RV FWHM Contrast BIS log R (cid:48) HK σ logR (cid:48) HK Moon flag[d] [m s − ] [m s − ] [km s − ] [m s − ] [dex] [dex]2456100.608 -10608.31 2.31 6.745 48.83 -41.90 -4.975 0.039 02456100.629 -10610.99 2.09 6.735 48.91 -40.73 -5.039 0.031 02456101.606 -10613.07 3.56 6.732 48.70 -55.78 -4.961 0.060 0... ... ... ... ... ... ... ... .N ote —Table 1 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content. unidentified factors can determine whether or not contamina-tion is negligible. An in-depth analysis of the outcome of theobservations is then advised when measuring RVs for faintstars. KEPLER
PHOTOMETRY
Kepler -19 was initially observed in long-cadence (LC)mode during the first 0-2 quarters, and then in short-cadence(SC) mode from quarter 3 until the end of the mission in 2013(at quarter 17). At the time of publishing, B11 had at theirdisposal only the first 8 quarters. In an e ff ort to make useof any additional information coming from Kepler photome-try, we redetermined the transit times for all the quarters nowavailable, in both LC and SC light curves. Quarters alreadyanalyzed by B11 were examined as well, to validate our T determination and to provide a homogeneous set of measure-ments.Transit identification was performed by propagating thelinear ephemeris of B11, with the inclusion of 3 hours ofpre-ingress and post-egress around the expected transit time.For each transit time, we firstly detrended the transit lightcurve with a polynomial between 1st and 10th degree, withthe best-fit degree chosen according to the Bayesian infor-mation criterion (BIC). Then, we determined a new T guesswith an automatic selection among di ff erent search algo-rithms fitting the Mandel & Agol (2002) transit model, im-plemented in PyTransit (Parviainen 2015), fixing all otherparameters to the literature value. Finally, the T s were re-fined using JKTEBOP program (Southworth et al. 2004) andthe associated errors were determined with a classical boot-strap approach.Transit times from LC and SC light curves were matchedtogether, keeping the SC measurements when available.Transit times are reported in Table 2. A comparison withB11 measurements of the observed minus predicted time of Levenberg-Marquardt (Mor´e et al. 1980), Nelder-Mead (Nelder &Mead 1965; Wright 1996), COBYLA (Powell 1994) as implemented in scipy.optimize , and the A ffi ne Invariant Markov Chain Monte Carlo(MCMC) Ensemble sampler implemented in emcee (Foreman-Mackey et al.2013) available at https://github.com/dfm/emcee Available at https://github.com/hpparvi/PyTransit transit ( O − C ), using their linear ephemeris for both datasets,is shown in Figure 1. The scatter of the residuals, well withinthe error bars, shows that the methodologies are perfectlyconsistent, i. e., that we are limited by photon noise, datasampling and / or unknown systematics rather than the exactprocedure followed to measure the transit times. Ballard et al. (2011)This work O - C [ m i n ] −15−10−505101520 r e s [ m i n ] −10−5051015 transit number (N) Figure 1 . The di ff erence between the observed and the predicted(from the linear ephemeris) times of transit for K Kepler quarters. In the lower plot the di ff er-ence between the two measurements is shown for the data pointsin common, with the error bars obtained by summing in quadraturethe errors from the two estimates. The small scatter of the residualswith respect to the size of the error bars demonstrate that we are notinfluenced by the exact methodology used to measure the T s. Due to an error in
Kepler archiving system, at the time ofB11 publication the time stamps of all the
Kepler light curveswere reported in Coordinated Universal Time system (UTC)instead of the Barycentric Dynamical Time system (TDB) . http://archive.stsci.edu/kepler/timing_error.html M alavolta et al . Table 2 . Transit times of
Kepler -19 from Q0-Q17.Transit Number T [BJD UTC ] σ T [d]0 2454959.7074 0.00141 2454968.9935 0.00232 2454978.2801 0.0020... ... ...N ote —Table 2 is published in its entirety in the machine-readableformat. A portion is shown here for guidance regarding its formand content. While this error was not a ff ecting the internal consistency ofB11 analysis, it must be taken into account when comparingtime-series with timing accuracy better than a few minutes.We corrected for this error before comparing B11 data withour new T measurements. PHYSICAL PARAMETERSAtmospheric stellar parameters of
Kepler -19 were deter-mined in B11 using Spectroscopy Made Easy (
SME , Valenti& Fischer 2005). Since this method may su ff er of correla-tion between derived parameters (Torres et al. 2012), andhaving at our disposal several high-resolution spectra fromHARPS-N, we decided to carry out an independent deter-mination with an alternative approach, i. e., equivalent widthmeasurements of individual spectral lines instead of fittingof the whole spectrum. We used all the spectra free from skycontamination to obtain a coadded spectrum with an averageS / N of 350.Stellar atmospheric parameters were determined using theclassical line-of-growth approach. For this purpose we usedthe 2014 version of the line analysis and synthetic code
MOOG (Sneden 1973), which works in the assumption of lo-cal thermodynamic equilibrium (LTE), and the ATLAS9 gridof stellar model atmosphere from Castelli & Kurucz (2004)with the new opacity distribution functions and no convec-tive overshooting. Equivalent Width measurements were car-ried out with the code
ARESv2 (Sousa et al. 2015) coupledwith the updated linelist of Malavolta et al. (2016), wherethe oscillator strength of the atomic lines have been modifiedto correctly take into account the chemical abundances fromAsplund et al. (2009).Temperature and microturbulent velocity were determinedby minimizing the trend of iron abundances from individuallines with respect to excitation potential and reduced equiva-lent width respectively, while the gravity log g was adjustedby imposing the same average abundance from neutral andionized iron lines. For a detailed description of the proce-dure for the atmospheric parameters and associated errors werefer the reader to Dumusque et al. (2014). The derived at- Available at Available at
Table 3 . Astrophysical parameters of the star.Parameter B11 This work T e ff [K] 5541 ±
60 5544 ± g . ± .
10 4 . ± . ξ t [km s − ] - 0 . ± . / H] − . ± . − . ± . (cid:63) [M (cid:12) ] 0 . ± .
040 -R (cid:63) [R (cid:12) ] 0 . ± .
018 -Age (Gyr) 1 . ± . (cid:48) HK − . ± . − . ± . mospheric parameters are summarized in Table 3.Our stellar atmospheric parameters agree within the un-certainties with the ones determined by B11, including thesurface gravity which is usually the parameter most di ffi cultto derive, and there is only a di ff erence of 3 K in T e ff despitethe use of two complementary approaches and independentdatasets. For this reason we adopted their determination formass and radius of the star and physical radius of K
19b basedon light curve analysis. STELLAR ACTIVITYRecently a considerable e ff ort has been devoted to analyzethe e ff ect of stellar activity on RV measurements as well as T determination (Mazeh et al. 2015; Ioannidis et al. 2016).The HARPS-N DRS automatically delivers several diagnos-tics for activity such as the Full Width Half Maximum of theCCF, the bisector inverse span and the log R (cid:48) HK index, whileseveral other chromospheric indexes such as the H α index(Gomes da Silva et al. 2011; Robertson et al. 2013) can be de-termined from the spectra themselves . In Figure 2 the analy-sis of BIS and log R (cid:48) HK are reported as representative of all theindexes. For each index we checked the presence of any cor-relation with time, either by visual inspection (upper panelsof Figure 2) and with the Generalized Lomb-Scargle (GLS)periodogram (Zechmeister & K¨urster 2009), where the 1%and 0.1% false alarm probability (FAP) have been computedwith a bootstrap approach (middle panels of Figure 2). Fi-nally, the presence of any correlation with RVs was verifiedby calculating the Spearmans rank correlation coe ffi cient ρ ,the slope of the linear fit m with its error and the p-valueusing the weighted least-square regression (lower panels ofFigure 2). We omitted the FWHM from the analysis since afew changes of the spectrograph focus during the first yearof observations modified the instrumental profile, hence themeasured FWHM, without however a ff ecting the measuredRVs. The code to retrieve the activity indexes is available at https://github.com/LucaMalavolta/ StatsModels available at http://statsmodels.sourceforge.net/ he Kepler -19 system l og R ' M W -5.20-5.10-5.00-4.90-4.80 BJD-2450000 [d]6200 6400 6600 6800 7000 log R' HK G L S po w e r l og R ' HK -5.20-5.10-5.00-4.90-4.80 RV [Km s -1 ]−10.635 −10.625 −10.615 −10.605 B I S [ m s - ] −60−50−40−30−20−10 BJD-2450000 [d]6200 6400 6600 6800 7000 Bisector Inverse Span G L S po w e r B I S [ m s - ] −60−50−40−30−20−10 RV [Km s -1 ]−10.635 −10.625 −10.615 −10.605 Figure 2 . The bisector inverse span (panels on the left side) andthe log R (cid:48) HK index (panels on the right side) are shown as exampleof the analysis conducted on CCF asymmetry and activity indexes.Upper panels: indexes as a function of time, the seasonal medianswith the first and third quartiles indicated in red. Middle panels:GLS periodograms of the indexes, the rotational period of the staris indicated with a red vertical line. The 1% and 0.1% FAP levels aredisplayed as dashed and dotted horizontal lines, respectively. Lowerpanels: indicators as a function of RV. The best fit is represented bythe dashed red line. The absence of significant peaks in the periodogram of theindexes under analysis around the expected rotational periodof 32 days (from B11 following Noyes et al. 1984), 34 ± ± (cid:48) HK , followingNoyes et al. 1984 and Mamajek & Hillenbrand 2008 respec-tively), and the lack of statistically significant correlationsbetween the indexes and RVs, confirmed the low activitylevel of the star already deduced by B11 and consistent withlog R (cid:48) HK = − . ± . .We also searched for stellar variability on the most re-cent Kepler Pre-search Data Conditioning (PDC) light curve,which presents several improvements in correcting instru-mental trends (and thus is better suited to search for activitymodulation) with respect to the light curve available at thetime of B11 publication. For most of the time the star is pho-tometrically quiet, while in some parts of the light curve aclear signal, likely due to stellar activity, is detected. We ap-plied the autocorrelation function technique over these por-tions of light curve and we estimated a rotational period of30 days. This signal is characterized by a short time scale ofdecay (a few rotational periods) and a rapid loss in coherence. Note however that this value could be a ff ected by interstellar mediumabsorption and be higher than measured (Fossati et al. 2017). We comment on the impact on RV in section 7. TTV ANALYSISDynamical analysis of the system was performed with
TRADES (TRAnsits and Dynamics of Exoplanetary Sys-tems, Borsato et al. 2014), a N-body integrator with the capa-bility of fitting RVs and T s simultaneously to determine theorbital parameters of the system through χ minimization.Since only the innermost planet is transiting, it is extremelydi ffi cult to constrain the orbital parameters of all the planetsin the system by TTV alone. In fact, attempts to fit the T swith a two-planet model resulted in a strong degeneracy be-tween the mass and the period of the second planet, i. e., O-Cdiagrams with similar shape could be produced by jointly in-creasing the mass and the period of K T s however can still give us upper limits onthe mass and period of the non-transiting planet if we com-pare the outcome of the dynamical simulations with the max-imum mass compatible with the observed semi-amplitude ofthe RVs.We proceeded as follows. We used TRADES to perform afit of the T s by assuming a two-planet model and fixing themass of K
19b to a grid of values between 2.5 and 20 M ⊕ anda spacing of 2.5 M ⊕ , with its period already measured fromthe Kepler observations. To each point of this grid we as-signed several values for the mass of K
19c randomly selectedbetween 5 and 250 M ⊕ . For each combination of K
19b and K
19c masses, the other orbital parameters of the system wereleft free to vary and their best-fitting values were determinedby χ minimization through the Levenberg-Marquardt algo-rithm. We discarded those solutions with at least one planethaving eccentricity greater than 0.3, assuming that interact-ing planets meeting this condition are likely to be unstable.In Figure 3 we show the results obtained for the mass of K
19c as a function of its period. The expected RV semi-amplitude of K
19c as a function of mass and period are su-perimposed. We note that many of the orbital configurationsreported in the plot could be unstable, since dynamical sta-bility was not checked yet at this stage (stability analysis isintroduced in Section 8). We can attempt to estimate a lowerlimit to the periods and masses of the non-transiting planetsaccording to the observed semi-amplitude of the RVs, by tak-ing advantage of the correlation between the mass and periodof K − ,so if we keep into account the additional signal of K
19b (afew m s − in the case of a Neptune-like density), K
19c ampli-tude should necessarily lie below the K =
10 m s − line. Thisfact suggests that a short period ( (cid:46)
50 days) for K
19c shouldbe expected, while nothing can be said regarding K
19b fromTTV alone. We remind the reader that this analysis cannot be TRADES is available at https://github.com/lucaborsato/trades M alavolta et al . K19b Mass [M e ] K c =10 ms -1 K c =20 ms -1 K c =30 ms -1 M a ss K c [ M e ] Figure 3 . The mass of K
19c versus its period, determined by fitting K T s for a grid of masses of the two planets. The three linesrepresent the expected RV semi-amplitude of K
19c only as a func-tion of its mass and period. A linear fit to the data is marked with adash-dot blue line. considered conclusive due to the reduced number of pointsand their scatter around the best linear fit. RV ANALYSISWhile only one planet is transiting
Kepler -19, the presenceof at least one additional planet was inferred by the presenceof TTVs. To reveal such planets, we first performed a fre-quentist analysis on the RV dataset by computing the GLSperiodogram and iterating over the residuals until no signifi-cant periodicity was present. This analysis revealed two sig-nals at 28.6 days and 62.3 days with false alarm probabilitylower than 1% , while the signal at 9.3 days corresponding to K
19b was barely detected (Figure 4).The signal at 28 days is very close to the second order 3:1mean motion resonance (MMR), which is among the possi-ble configurations listed by B11 as the cause of the TTVs of K (cid:39) − ), we can compare it with the probability of observ-ing K
19c near another MMR resonance. In order to do so,we scrambled once again the RV observations and calculatedthe fraction of periodograms which had a stronger peak withrespect to the untouched dataset, in a frequency range of 5%around the interior and exterior 1:2, 2:3, 3:4 first order reso-nances, the 1:3 and 3:5 second order resonances, and the 1:4and 1:5 higher order resonances. The FAP of the signal at 28days computed in this way decreases to 8 × − , i. e., it isunlikely that the signal at 28 days is a spurious signal causedby a planet in another MMR configuration. New CCDOld CCD R V - . [ m s - ] −20−15−10−5051015 BJD-2450000 [d]6000 6200 6400 6600 6800 7000 G L S po w e r Figure 4 . Upper panel: the RV dataset. The vertical line dividesdata taken before and after the substitution of HARPS-N CCD.Lower panels: GLS periodograms computed starting from the ini-tial dataset and iterating over the residuals. The 1% and 0.1% falsealarm probabilities were estimated following a bootstrap approachfor each periodogram. Red lines mark the best-fit period at 28.6days, 9.29 days and 62.3 days.
The signal at (cid:39) . Kepler light curve, which however have the characteristic of rapidlydecaying and reappearing later at di ff erent phase and inten-sity. We checked if the RV signal has the same properties byperforming a jacknife analysis by splitting the RV dataset intwo parts (at BJD = . We obtainedK = . ± . − , φ = . ± . = . ± . − , φ = . ± . / N is increasing with the squareof the number of observations, as expected for a coherent(planetary) signal, although the plot is heavily a ff ected bythe poor sampling and the signal of K
19b near the 3:1 MMR,as we verified by performing the same analysis on syntheticdatasets with the same temporal sampling. The absence ofRV modulation due to stellar activity is supported by theanalysis performed in Section 5, where no periodicity or cor-relation with RVs is seen for any of the activity indexes.The semi-amplitudes of these signals however cannot be We imposed as a prior the period P = . ± .
02 day to compensatefor the scarcity and poor sampling of the split datasets he Kepler -19 system Table 4 . Orbital parameters for a 3-planet model for the
Kepler -19system obatined from RV only, except for the period of K
19c whichis constrained by the
Kepler light curve. The reference time for theorbital elements is T ref = . K K K . ± − . ± .
24 63 . ± .
3K [m s − ] 2 . ± . . ± . . ± . φ [deg] 194 . ± . ±
43 174 ± ⊕ ] 7 . ± . . ± . . ± . . ± − . ± .
27 62 . ± .
3K [m s − ] 2 . ± . . ± . . ± . φ [deg] 198 ± ±
29 172 ± √ e cos ω . ± . . ± . . + . − . √ e sin ω . ± . . ± . − . . − . Mass [M ⊕ ] 8 . ± . . ± . . ± . e ≤ . ≤ . ≤ . ω [deg] unconstrained unconstrained − + − inferred by the frequentist analysis alone, since an o ff set inRV between the data taken before and after September 2012may exist due to the failure of the first CCD of HARPS-N(see Bonomo et al. 2014). The value of this o ff set cannotbe determined a priori, even when observations of RV stan-dard stars are available, since it may depends primarily onthe spectral type of the star (i. e., the two CCD may have adi ff erent e ffi ciency as a function of wavelength). Introducingan RV o ff set as a free parameter is a common procedure fornon-overlapping datasets, e. g., see Benatti et al. (2016).We then performed a tentative fit using the Markov-ChainMonte Carlo (MCMC) code PyORBIT (Malavolta et al.2016), allowing for two di ff erent systemic velocity γ of thestar and using a 3-planet model to fit the data. We attemptedtwo di ff erent fits with both circular and Keplerian orbits. Forplanet c and d we set uniform priors in the range [10 ,
50] daysand [50 ,
90] days respectively, while the other parameterswere left to vary within their physically meaningful range(e. g., positive-definite RV semi-amplitude). The results areshown in Table 4.Eccentricity and argument of pericenter of each of the threeplanets are poorly constrained by the RVs alone, thus a ff ect-ing the precision on mass of the presumed planets. A com-bined analysis of RVs and TTV is then required to unambigu-ously detect and characterize the planets in the system. COMBINED RV AND TTV ANALYSISThe amplitude of dynamical perturbations between planetsis very sensitive to the eccentricity and angular parameters Available at https://github.com/LucaMalavolta/PyORBIT (i. e., argument of pericenter ω and mean anomaly at refer-ence time M ) of the planetary orbits, which however can beonly poorly constrained by the RVs, especially for planets innearly circular orbits. For this reason we expect an overallimprovement on the precision of the orbital parameters bysimultaneously fitting RVs and TTVs.Dynamical simulations are extremely time-consuming,and we have to use all the information at our disposal to re-duce the extension of the parameter space. We used the re-sults from the RV fit in Section 7 to put a constraint to therange of period, mass and orbital phase of each planet.We followed an iterative approach to avoid being trappedin a local minimum of the χ . We started 10 separate runs of TRADES with the Particle Swarm Optimization algorithm andloose priors on the periods ( ±
10 days around the expectedperiod of each planet from the RV analysis), masses (between0 and 40 M ⊕ ) and eccentricities ( e ≤ . TRADES again on a range of parameters which was half the size than inthe previous run and centered of the outcome of the previousrun with the lowest χ among those that satisfied the stabilityrequirement. Convergence was considered achieved when allthe runs resulted in similar parameters (within 5% from themean) and similar χ (10% from the mean).Following Gladman (1993), during the numerical integra-tion we checked the stability criterion for each pair of plan-ets: ∆ = √ R H ( i , j ), where ∆ = a j − a i is the semi-majoraxis di ff erence between the j th and i th planet, and R H ( i , j ) isthe mutual Hill radius between planet i and j . At the endof each TRADES fit we performed a N-body integration with
SymBA (Duncan et al. 1998) and checked the stability of theresult with the Frequency Map Analysis tool (FMA, Laskaret al. 1992; Laskar 1993a,b) with the prescriptions of Marzariet al. (2002). A system is considered stable if the coe ffi cientof di ff usion is lower than 10 − , in an unstable or chaotic stateotherwise.We used the value measured by B11 for the inclination i of K Kepler multi-planetsystems (Fabrycky et al. 2014) and the additional informa-tion coming from systems characterized with high-precisionRVs (Figueira et al. 2012) to impose coplanarity with K Ω was fixed to zero for all the planets. This assumptionallowed us to drastically reduce computational time. The or-bital period and mass of the planet, the eccentricity e , theargument of periastron ω and the mean anomaly at the refer-ence epoch M were left as free parameters.Di ff erently from (Borsato et al. 2014), we fitted √ e cos ω and √ e sin ω instead of e cos ω and e sin ω , and the mean lon- M alavolta et al .gitude at the reference epoch λ = ω + M + Ω (where Ω isthe longitude of the ascending node, fixed to zero for beingunconstrained by the data) instead of M . We used scaledstellar and planetary masses as commonly done in TTV anal-ysis (e. g., Nesvorn´y et al. 2012). A RV o ff set between thedata taken before and after September 2012 was included asa free parameter to take into account the change of the CCD(see Section 7).We used the solution obtained with the global explorationof the parameter space as a starting point for the Bayesiananalysis. For this purpose we expanded TRADES functionali-ties with the emcee package (Foreman-Mackey et al. 2013),an a ffi ne invariant Markov chain Monte Carlo (MCMC) en-semble sampler, and made it available to the community inthe TRADES repository. Following Feigelson & Babu (2012),we calculated the log-likelihood ln L from the χ using Equa-tion 1, where do f is the degree of the freedom of the problem.ln L = − ln (2 π ) dof2 − (cid:80) ln σ − χ ff erent scenarios. The first model as-sumed that the signal at (cid:39)
28 days (the first one being de-tected in the RV periodogram) is the only planet in the sys-tem other than K (cid:39)
28 days isdue to activity (without any e ff ects on the T s) and that thesystem consists of two planets at (cid:39) . (cid:39)
62 days. Thethird model assumed that all the RV signals have planetaryorigin. The results are presented in the following subsections.We note that a planet in a strong MMR with K
19b could stillproduce the TTVs while having a RV semi-amplitude belowour detection sensitivity. Following B11, the TTVs of K ⊕ respectively. The signal of the perturb-ing planet would nave an RV semi-amplitude on the orderof [1, 0.6, 0.3] m s − , i. e., beyond the reach of modern ve-locimeters given the magnitude of our target. This scenariohowever requires an activity origin for the 28 days signals,which is not supported by the analysis of Section 5 and thecoherent nature over years of the RV signal in opposition tothe short time-scale decay of the spots observed in Kepler photometry, as shown in Section 7. For these reasons we canregard this scenario as unlikely.8.1.
Two-planet model
We tested the two-planet model following the steps de-scribed in Section 8, with the only additional constraint of K
19c having a period lower than 40 days. We ran the MCMCsampler for 50000 steps, with a number of chains twice thedimensionality of the problem. We checked the convergenceof the chains using the Gelman-Rubin statistic ( ˆ R < . χ tothe median of its distribution and with each parameter withinits confidence interval. Confidence intervals are computedby taking the 15.87th and 84.14th percentiles of the distribu-tions, and are reported in the table with respect to the selectedsample. 8.2. Two-planet and stellar activity model
The presence of the activity signal can potentially a ff ectthe analysis, but TRADES is not equipped to deal with non-planetary signals in RV datasets. For this reason we decidedto remove such signal from the RV time series before start-ing the global exploration of the parameter space, using asdataset the residuals of the first iteration of the GLS analysisin Section 7. The global solution was used as starting pointfor a MCMC analysis with
PyORBIT , this time using the orig-inal RV dataset and the activity signal modeled with a Kep-lerian curve (as similarly done in Pepe et al. 2013). RVs andT s calculations regarding the two interacting planets wereperformed by calling the dynamical integrator of TRADESthrough a FORTRAN90 wrapper. As additional test-cases, werun TRADES using the two-planet model and no activitymodeling on both the original RV dataset (imposing a lowerlimit of 40 days for the period of K GLS-cleaned dataset. Inall cases we repeated the analysis without imposing any con-straint on eccentricity, since the eccentricities of both planetswere moving towards the upper boundary imposed in Sec-tion 8. We followed the same methodology described in Sec-tion 8.1 to run the MCMC and extract the results reported inTable.The eccentricity of K
19b is extremely high in all the caseswe considered. The most likely explanation is that high ec-centricities for both planets are required to produce the sameT s while keeping their masses within the boundaries set bythe RVs. Following Burke (2008) we computed the ratio τ of the transit duration of an eccentric orbit with respect to acircular orbit, and we obtained τ = . ± .
04 for the orig-inal RV dataset, τ = . ± .
06 for the
GLS-cleaned case,and 0 . ± .
10 for the planets + activity case. All these valuesare well below the value of τ = . Kepler light curve. We note that the period-mass combina-tions we obtain for the outer planet fall in an empty regionof Figure 3, which however was obtained by selecting thosesolutions with e < . Three-planet model
Bayesian analysis for the three-planet model were per-formed using
TRADES combined with emcee (see Section 8).From a preliminary analysis we noticed that the chains were he Kepler -19 system ff ected by poor mixing, with a Gelman-Rubin statistic ˆ R (cid:39) . χ . After the first 150000 iter-ations we did not see any variation of the posterior distribu-tions while increasing the length of the chains, making usconfident of the robustness of our result. We finally built theposterior probability by drawing 40000 independent samplesfrom the chains, after removing the burn-in part and applyinga thinning factor equal to their auto-correlation time ( (cid:39) solution of the orbital parameters of the planets. InFigure 6 we show the solution overplotted on T s and RVs,with their respective residuals.In Table 6 we compare the outcome of the di ff erent modelsunder examination. The three-planet solution is the favoriteone according to the BIC ( having a ∆ BIC >
10 with re-spect to all the other models under analysis, Kass & Raftery1995), and the considerations at the end of Section 8.2 furtherstrengthen this result. Our three-planet solution has a total re-duced χ r = .
25 ( χ = .
7, BIC = . .The standard deviation of the residual RVs is 2.9 m s − , con-sistent with the average error of the measurements, and whileseveral peaks can be found in the periodogram, none of themreaches the 1% false-alarm probability threshold (Figure 7).Similarly, the standard deviation of the TTV of 2.2 minutesis consistent with the average error of 2.1 minutes.For the three planets we determined their masses with aprecision of σ M b = . ⊕ (error of 19% on planetarymass) for K σ M c = . ⊕ (error of 21%) for K
19c and σ M d = . ⊕ (error of 17%) for K ff ectof orbital inclinations, assuming as representative distribu-tions for the latter a normal distribution with i c = i d = . ◦ and dispersion σ i = ◦ for K
19c and σ i = . ◦ for K K
19c and K
19d should have atransit depth greater than 0.5 mmag (Winn 2010), a visual in-spection of the
Kepler light-curve confirmed that the planetsare not transiting.Our solution agrees with the solution in Table 4, except forthe mass of K
19c , with a 2- σ di ff erence between the RV-only These two values do not correspond to the individual χ red of RVs andTTVs, since the sum of residuals of each dataset is divided by the total num-ber of data points. and RV + TTV determinations. The origin of this discrepancyis likely due to the uncertainty in the orbital phase and eccen-tricity from the RV-only analysis, with one planet absorbingthe signal of the other, possibly coupled with the e ff ect de-scribed in Section 10.The RV o ff set between the old and new CCD is (cid:39) − ,i. e., within the RV noise. We do not expect any influence ofthe RV o ff set on the outcome of the analysis, since the threeplanets have periods fully contained within the time span ofboth old and new CCD datasets. MUTUAL INCLINATIONSIn Section 8 we assumed co-planarity with K
19b (the onlyplanet with known inclination) to derive the orbital parame-ters of the additional planets in the system. To check if thisassumption was still valid after determining the orbital pe-riod of the additional planets, we determined the upper limiton inclination for which a transit is visible for a given planet,by inverting Equation 7 of Winn (2010) with the assumptionfor the impact parameter b tra =
1. To properly take into ac-count variation with time of e and ω induced by dynamicalinteractions, we selected 1000 ( e , ω ) pairs for each planet byrandomly sampling in time the integration our best solutionover the Kepler observational time span. Error in the stel-lar radius was included by generating random samples froma normal distribution with mean R (cid:63) and standard deviation σ R (cid:63) as reported in Table 3. We obtained i min = . ± . ◦ for K
19c and i min = . ± . ◦ for K i b = .
94 (from B11), then the coplanarity assump-tion cannot hold for K i c and i d (from 60 to 120degrees, with step of 0 . ◦ for i c and 0 . ◦ for i d ), with theremaining orbital parameters fixed to our best solution, anddetermining the reduced χ as in Section 8. As shown inFigure 8, the reduced χ increases rapidly with K
19c goingfarther from coplanarity with K ≈ . | i b − i c | (cid:39) ◦ , while K
19d can span a larger interval in in-clinations without a ff ecting the outcome (although increasing | i b − i d | would a ff ect negatively the stability of the system).Assuming i c = .
72 (grazing scenario) the outcome of thefit is nearly the same χ = .
26, which is very close to thevalue obtained when assuming coplanarity ( χ = .
25, seeSection 8). It is likely that the system is very close to orbitalalignment, as observed for the majority of transiting multi-planet systems, e. g., Figueira et al. 2012; Fabrycky et al.2014; Ballard & Johnson 2016; Becker & Adams 2016.
DYNAMICAL VERSUS NON-INTERACTINGORBITSDynamical simulations include by definition the e ff ects ofgravitational interactions between planets. This is di ff erentfrom the usual approach followed in the exoplanet literature,where a series of non-interacting Keplerian orbits are used to0 M alavolta et al . Table 5 . Orbital parameters of the planets in the
Kepler -19 system obatined from TTV + RV MCMC analysis using di ff erent assumptions forthe number of planets and the stellar activity. The reference time for the orbital elements is T ref = . ⊕ ] √ e cos ω √ e sin ω λ [deg] e ω [deg] M [deg]Two-planet model, P c <
40 days K
19b 9 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . K
19c 28 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . Two-planet model, P c >
40 days K
19b 9 . + . − . . + . − . − . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . K
19c 63 . + . − . . + . − . − . + . − . − . + . − . . + . − . . + . − . − . + . − . . + . − . Two-planet model,
GLS-cleaned
RV dataset K
19b 9 . + . − . . + . − . − . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . K
19c 63 . + . − . . + . − . − . + . − . − . + . − . . + . − . . + . − . − . + . − . . + . − . Two-planet and activity model K
19b 9 . + . − . . + . − . − . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . Act 28 . + . − . (2 . + . − . ) a . + . − . . + . − . . + . − . . + . . . + . − . . + . − . K
19c 63 . + . − . . + . − . . + . − . − . + . − . . + . − . . + . − . − . + . − . . + . − . Three-planet model K
19b 9 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . K
19c 28 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . K
19d 62 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . a Semi-amplitude of the signal in m s − Table 6 . Statistical indexes for the di ff erent models under exams. The number of parameters in the fit, the degree of freedom ( do f ), the χ andits reduced value, the log-likelihood ln L and the Bayesian Information Criterion (BIC) are reported.Model N. parameters do f χ χ ln L BIC2 planets P c <
40 days 12 223 323 . + . − . . + . − . . + . − . . + . − . c >
40 days 12 223 324 . + . − . . + . − . . + . − . . + . − . GLS-cleaned
RVs 12 223 318 . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . . + . − . derive the planet parameters in multiple system. While theassumption of negligible interactions between planets mayhold in most of the cases, in the presence of TTVs we knowthat such interactions are happening. It is then worthy to an-alyze the di ff erences in the RVs between the two approaches,i. e., dynamical vs. non-interacting Keplerian orbits. In or-der to do so, we have simulated the expected RVs of the Kepler -19 system at the observational epochs of our dataset,using the planetary parameters in Table 5 and assuming non-interacting orbits, and then we have subtracted the outcometo the dynamically-derived RVs from the same orbital param-eters. The results are shown in Fig. 9. The RV
TTV − RV Kep residuals show a peak-to-peak variation of 0.30 m s − with aprominent periodicity at 29.3 days, i. e., very close to the or-bital period of K ff erence between dynamical and Keplerian RVs di-verges while moving further from the reference time of theorbital parameters T ref . The reason resides in the small varia- tions of the orbital parameters with time caused by planet in-teractions, which are implicitly taken into account in dynam-ical integration but ignored in the Keplerian approach. Forthe Kepler -19 system, after two years of observations from T ref the peak-to-peak di ff erence has already reached (cid:39) .
6m s − , i. e., a value that can negatively a ff ect the mass deter-mination of the planets. Assuming non-interacting planetswhen modeling the RVs can thus mislead the determinationof the orbital parameters, e. g. in the case of dataset spanningseveral years. DISCUSSIONThe transiting planet
Kepler -19b was previously validatedby Ballard et al. (2011). A period of 9.23 days and an upperlimit of 20 M ⊕ were determined from the analysis of the lightcurve and high-resolution spectroscopy. From the presenceof TTV in 8 Kepler quarters they deduced the presence ofa second, non-transiting planet with period (cid:46)
160 days andmass (cid:46) jup . In this paper we presented the first precise he Kepler -19 system . . . . M b M ? (cid:20) M ⊕ M (cid:12) (cid:21) λ d [ ◦ ] -0.3717-0.18470.00230.18920.3762 √ e s i n ω d -0.4593-0.2514-0.04340.16460.3725 √ e c o s ω d P d [ d a y s ] M d M ? (cid:20) M ⊕ M (cid:12) (cid:21) λ c [ ◦ ] -0.02920.06120.15170.24210.3326 √ e s i n ω c √ e c o s ω c P c [ d a y s ] M c M ? (cid:20) M ⊕ M (cid:12) (cid:21) λ b [ ◦ ] √ e s i n ω b -0.03030.02740.08500.14260.2002 √ e c o s ω b P b [ d a y s ] M b M ? (cid:20) M ⊕ M (cid:12) (cid:21) . . . . P b [days] P b [days] - . . . . √ e cos ω b √ e cos ω b . . . . √ e sin ω b √ e sin ω b . . . . λ b [ ◦ ] λ b [ ◦ ] . . . . M c M ? (cid:20) M ⊕ M (cid:12) (cid:21) M c M ? (cid:20) M ⊕ M (cid:12) (cid:21) . . . . P c [days] P c [days] . . . . √ e cos ω c √ e cos ω c - . . . . √ e sin ω c √ e sin ω c . . . . λ c [ ◦ ] λ c [ ◦ ] . . . . M d M ? (cid:20) M ⊕ M (cid:12) (cid:21) M d M ? (cid:20) M ⊕ M (cid:12) (cid:21) . . . . P d [days] P d [days] - . - . . . √ e cos ω d √ e cos ω d - . - . . . √ e sin ω d √ e sin ω d λ d [ ◦ ] Figure 5 . Posterior distributions of the fitted parameters for the three-planet model for the
Kepler -19 system. Dashed blue lines identifythe reference solution listed in Table 5. mass measurement for K
19b ( 8 . ± . ⊕ ) and the charac-terization of two non-transiting Neptune-mass companions, K
19c ( P = . ± .
01 days, M = . ± . ⊕ ) and K P = . ± . M = . ± . ⊕ ), obtained by simul-taneously modeling TTVs and RVs through dynamical inte-gration. We excluded stellar activity as a possible origin forthe RV signal at P (cid:39)
28 days using all the activity diagnosticsat our disposal, including the latest reduction of
Kepler lightcurve. Nevertheless we performed standard model selectionbetween the three-planet model and two-planet models, ei- ther with the hypothesis of the outer planet causing the TTVsand di ff erent assumptions for the activity signal, and with thehypothesis of the perturber having period (cid:39)
28 days and noouter planets. In all cases the three-planet model resulted thefavorite one with high degree of confidence. A planet in astrong MMR with K
19b could still produce the TTVs whilehaving a RV semi-amplitude below our detection sensitivity,but only at the condition of assuming that the 28 signal is dueto stellar activity, which is not supported by our data.With a period ratio P c / P b = .
09, the system is very close2 M alavolta et al . O - C [ m ] -10-50510 simulations observations r e s [ m ] -10-50510 BJD
TDB - 2450000.05000 5250 5500 5750 6000 6250 R V [ m / s ] -15-10-50510 r e s [ m / s ] -10-50510 BJD TDB - 2450000.06200 6400 6600 6800 7000
Figure 6 . Top panels: Observed-Calculated (O-C) transit timesand residuals for K Kepler -19system. Red empty circles represent the predicted valuesfrom the solution obtained with
TRADES . The radial velocity curveof the
Kepler -19 system obtained from the dynamical integration isshown in gray. to a 3:1 mean motion resonance. The sinusoidal shape of theTTV is then induced by the 3:1 MMR of the inner planets,with a modulation caused by the outer planet. We performeda comparison between dynamical RVs with those calculatedassuming non-interacting planets, using the orbital parame-ters of the
Kepler -19 system, and showed that, for this spe-cific system, the di ff erence for a dataset spanning severalyears is at the limit of detection with the state-of-the-art in-struments used for planet search and characterization.Our new determination of K
19b mass is in disagreementwith the most-likely value of 1.6 M ⊕ (semi-amplitude of 0.5m s − ) obtained by B11 while attempting a fit with only 8Keck-HIRES RVs and assuming only one planet in the sys-tem. As a consequence of this assumption, they derived alikely RV jitter of (cid:39) − , which is not confirmed by ouranalysis. In general, RV analysis performed on a small num-ber of measurements should always be handled with extremecare, since the results could be a ff ected by additional not-transiting planets with little influence on the TTV of the tran-siting planets but with significant RV semi-amplitude, suchas K K
19b falls in the region of super-Earths with rocky coresand a significant fraction of volatiles or H / He gas, in op-position to the low-density planets characterized by TTVonly (Weiss & Marcy 2014; Jontof-Hutter et al. 2016), andwell separated from the group of rocky planets with radiismaller than 2 R ⊕ , as can be seen in Figure 10. If we as-sume that the planet composition is a mixture of H / He enve- lope with solar composition atop a rocky core with Earth-likerock / iron abundances, we can estimate the internal structureof K19b. By employing theoretical models from Lopez &Fortney (2014), and assuming an age of 2 Gyr, an envelopeof 0 . ± .
3% of the total mass is required to explain the ob-served mass and radius of K (cid:39) K
19c and K
19d to haveaccreted significantly larger envelopes than K ff ected by photo-evaporation due to the larger masses ofthese two planets. Using Equation 22 from Lee & Chiang(2015), we can get a rough estimate for the size of these en-velopes and find that K
19c should have an H / He envelope of (cid:39) .
5% of its total mass and a radius of (cid:39) . ⊕ , while K − ⊕ .Our results confirm that TTV and RV techniques can con-verge to planetary densities similar to the ones obtained withRV dataset only, when enough data is involved from bothsides. This result support the analysis of Ste ff en (2016),where the discrepancy between planetary density obtainedfrom TTV and RV noted by Weiss & Marcy (2014) is ex-plained as a selection e ff ect rather than an intrinsic problemof one of the two techniques. Previous notable examples ofagreement between RV and TTV masses are represented bythe WASP-47 (Becker et al. 2015; Weiss et al. 2016) and the Kepler -18 systems (Cochran et al. 2011), and the indepen-dent RV confirmation by Dai et al. (2016) of the TTV-derivedmasses of the planets in the K2-19 system (Barros et al.2015). It should be noted however that Nespral et al. (2016)derived an higher RV mass for K2-19b, while the mass ob-tained by combining the two datasets resides halfway the twoextreme values. In the case of the KOI-94 system, RV (Weisset al. 2013) and TTV (Masuda et al. 2013) planetary massesagree except for planet d , where the TTV mass measurementis half the mass obtained by high-precision RVs. Up to now,most of the targets characterized with TTVs are too faint for aprecise, independent mass characterization with RVs. Futurespace-borne missions such as TESS (Transiting ExoplanetSurvey Satellite, Ricker et al. 2014) and PLATO (PLAne-tary Transits ans stellar Oscillations, Rauer et al. 2014) willfinally shed light on this problem by providing a large num-ber of targets bright enough for mass measurement with bothtechniques independently.The HARPS-N project was funded by the Prodex Pro- he Kepler -19 system Figure 7 . Orbital solution and RV residuals for K
19b (upper-left panels), K
19c (lower-left panels) and K
19d (upper-right panels), phased onthe period of the corresponding planet after removing the RV contribution from the other planets. These plots have been obtained using non-interacting Keplerian orbits, and are intended for illustrative purpose only. In the lower-right panels, both the RV residuals after subtracting thedynamical solution from
TRADES and their periodogram show no evidence for additional signals that are statistically significant. . . . . . . . χ i d [ d e g ] i c [deg]
84 86 88 90 92 94 96
Figure 8 . Distribution of χ red as a function of K
19c and K
19d in-clinations, with the other orbital parameters fixed to our solution(Table 5). Contour lines for several values of the χ red are shown forreference. Values of inclination for which one of the planets wouldtransit are shadowed in gray. gram of the Swiss Space O ffi ce (SSO), the Harvard- Univer-sity Origin of Life Initiative (HUOLI), the Scottish Universi-ties Physics Alliance (SUPA), the University of Geneva, theSmithsonian Astrophysical Observatory (SAO), and the Ital-ian National Astrophysical Institute (INAF), University of St.Andrews, Queen’s University Belfast and University of Ed-inburgh. R V D yn - R V K e p [ m / s ] −0.4−0.200.20.4 BJD-2450000 G L S po w e r Period [d] R V D yn - R V K e p [ m / s ] -0.10-0.050.000.050.100.150.20 Phase −0.2 0 0.2 0.4 0.6 0.8 1
Figure 9 . Analysis of the RV di ff erence ( ∆ RV) obtained by sub-tracting to the RV from dynamical simulation (RV
Dyn ) the RV fromnon-interacting Keplerian orbits (RV
Kep ), assuming planets with thesame orbital parameters in both cases. In the upper panel, ∆ RV atthe HARPS-N epochs and for a 5-year time span are shown as blackdots and a black line respectively. The reference time of the orbitalparameters T is marked with a dashed red line, to emphasize theincrease in ∆ RV when moving further from T . In the middle panel,the periodogram of such di ff erence reveals the presence of a peri-odicity at 29.3 d (red line). As comparison, the periods of the threeplanets in the system are highlighted in blue. The ∆ RV phased atsuch period is shown in the lower panel, with the sinusoidal modelmarked with a red line. alavolta et al . F p [F e ]10 100 1000 Kepler-19b σ ρ < 20%20% < σ ρ < 40%40% < σ ρ < 60%σ ρ > 60%1% H/He envelopePure WaterPure RockyEarth comp.Collisional strip R a d i u s [ R e ] e ]2 4 6 8 10 12 14 Figure 10 . Mass-radius for transiting exoplanets with measuredmasses b . Planets are color-coded according to the incident flux onthe planet F p , relative to the solar constant F ⊕ . Line thickness re-flects the precision on density measurements. Blue error bars iden-tify planets with TTV detection. The dashed lines are theoreticalmass-radius curves for di ff erent internal compositions, while Earth-like composition and a mixture of Earth-like core plus H / He en-velope of 1% in mass are represented by solid red and blue linesrespectively (Lopez et al. 2012). The shaded gray region marks themaximum value of iron content predicted by collisional stripping(Marcus et al. 2010). K
19b is highlighted with a black dot. a Automatically retrieved in February 2017 from the Nasa ExoplanetArchive, http://exoplanetarchive.ipac.caltech.edu b Automatically retrieved in February 2017 from the Nasa ExoplanetArchive, http://exoplanetarchive.ipac.caltech.edu
The research leading to these results received fundingfrom the European Union Seventh Framework Programme(FP7 / / Jet Propulsion Laboratoryfunded by NASA through the Sagan Fellowship Program ex-ecuted by the NASA Exoplanet Science Institute.AV is supported by the NSF Graduate Research Fellow-ship, Grant No. DGE 1144152.This publication was made possible by a grant from theJohn Templeton Foundation. The opinions expressed in thispublication are those of the authors and do not necessarilyreflect the views of the John Templeton Foundation. Thismaterial is based upon work supported by NASA under grantNo. NNX15AC90G issued through the Exoplanets ResearchProgram.PF acknowledges support by Fundac¸ ˜ao para a Ciˆencia ea Tecnologia (FCT) through Investigador FCT contract ofreference IF / / / FSE (EC) by FEDERfunding through the program “Programa Operacional de Fac-tores de Competitividade - COMPETE”, and further sup-port in the form of an exploratory project of referenceIF / / / CT0001.XD is grateful to the Society in Science − Branco Weiss Fel-lowship for its financial support.REFERENCES
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47,481Ballard, S., & Johnson, J. A. 2016, ApJ, 816, 66Ballard, S., Fabrycky, D., Fressin, F., et al. 2011, ApJ, 743, 200Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373Barros, S. C. C., Almenara, J. M., Demangeon, O., et al. 2015, MNRAS,454, 4267Becker, J. C., & Adams, F. C. 2016, MNRAS, 455, 2980Becker, J. C., Vanderburg, A., Adams, F. C., Rappaport, S. A., &Schwengeler, H. M. 2015, ApJL, 812, L18Benatti, S., Desidera, S., Damasso, M., et al. 2016, ArXiv e-prints,arXiv:1611.09873Bonomo, A. S., Sozzetti, A., Lovis, C., et al. 2014, A&A, 572, A2Borsato, L., Marzari, F., Nascimbeni, V., et al. 2014, A&A, 571, A38Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19Buchhave, L. A., Dressing, C. D., Dumusque, X., et al. 2016, AJ, 152, 160Burke, C. J. 2008, ApJ, 679, 1566Castelli, F., & Kurucz, R. L. 2004, arXiv, astro-phCochran, W. D., Fabrycky, D. C., Torres, G., et al. 2011, ApJS, 197, 7Cosentino, R., Lovis, C., Pepe, F., et al. 2014, in Society of Photo-OpticalInstrumentation Engineers (SPIE) Conference Series, Vol. 9147, Societyof Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 8Coughlin, J. L., Mullally, F., Thompson, S. E., et al. 2016, ApJS, 224, 12 Dai, F., Winn, J. N., Albrecht, S., et al. 2016, ApJ, 823, 115Demory, B.-O., Gillon, M., de Wit, J., et al. 2016, Nature, 532, 207Dressing, C. D., Charbonneau, D., Dumusque, X., et al. 2015, ApJ, 800,135Dumusque, X., Bonomo, A. S., Haywood, R. D., et al. 2014, ApJ, 789, 154Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146Feigelson, E. D., & Babu, G. J. 2012, Modern Statistical Methods forAstronomyFigueira, P., Santos, N. C., Pepe, F., Lovis, C., & Nardetto, N. 2013, A&A,557, A93Figueira, P., Marmier, M., Bou´e, G., et al. 2012, A&A, 541, A139Ford, E. B. 2006, ApJ, 642, 505Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP,125, 306Fossati, L., Marcelja, S. E., Staab, D., et al. 2017, ArXiv e-prints,arXiv:1702.02883Gelman, A., & Rubin, D. B. 1992, Statistical Science, 7, 457Gettel, S., Charbonneau, D., Dressing, C. D., et al. 2016, ApJ, 816, 95Gladman, B. 1993, Icarus, 106, 247Gomes da Silva, J., Santos, N. C., Bonfils, X., et al. 2011, A&A, 534, A30Ioannidis, P., Huber, K. F., & Schmitt, J. H. M. M. 2016, A&A, 585, A72Jontof-Hutter, D., Ford, E. B., Rowe, J. F., et al. 2016, ApJ, 820, 39 he Kepler -19 system Kass, R. E., & Raftery, A. E. 1995, Journal of the American StatisticalAssociation, 90, 773Laskar, J. 1993a, Physica D: Nonlinear Phenomena, 67, 257—. 1993b, Celestial Mechanics and Dynamical Astronomy, 56, 191Laskar, J., Froeschl, C., & Celletti, A. 1992, Physica D: NonlinearPhenomena, 56, 253Lee, E. J., & Chiang, E. 2015, ApJ, 811, 41Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59L´opez-Morales, M., Haywood, R. D., Coughlin, J. L., et al. 2016, AJ, 152,204Malavolta, L., Nascimbeni, V., Piotto, G., et al. 2016, A&A, 588, A118Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264Mandel, K., & Agol, E. 2002, ApJL, 580, L171Marcus, R. A., Sasselov, D., Hernquist, L., & Stewart, S. T. 2010, ApJL,712, L73Marzari, F., Tricarico, P., & Scholl, H. 2002, ApJ, 579, 905Masuda, K., Hirano, T., Taruya, A., Nagasawa, M., & Suto, Y. 2013, ApJ,778, 185Mazeh, T., Holczer, T., & Shporer, A. 2015, ApJ, 800, 142Mor´e, J. J., Garbow, B. S., & Hillstrom, K. E. 1980, User Guide forMINPACK-1, Report ANL-80-74Mortier, A., & Collier Cameron, A. 2017, ArXiv e-prints,arXiv:1702.03885Nelder, J. A., & Mead, R. 1965, The Computer Journal, 7, 308Nespral, D., Gandolfi, D., Deeg, H. J., et al. 2016, ArXiv e-prints,arXiv:1604.01265Nesvorn´y, D., Kipping, D. M., Buchhave, L. A., et al. 2012, Science, 336,1133Noyes, R. W., Weiss, N. O., & Vaughan, A. H. 1984, ApJ, 287, 769Parviainen, H. 2015, MNRAS, 450, 3233Pepe, F., Mayor, M., Galland, F., et al. 2002, A&A, 388, 632 Pepe, F., Cameron, A. C., Latham, D. W., et al. 2013, Nature, 503, 377Powell, M. J. D. 1994, A Direct Search Optimization Method That Modelsthe Objective and Constraint Functions by Linear Interpolation, ed.S. Gomez & J.-P. Hennart (Dordrecht: Springer Netherlands), 51–67Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38,249Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Proc. SPIE, Vol.9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, andMillimeter Wave, 914320Robertson, P., Endl, M., Cochran, W. D., & Dodson-Robinson, S. E. 2013,ApJ, 764, 3Rogers, L. A. 2015, ApJ, 801, 41Sneden, C. 1973, ApJ, 184, 839Sousa, S. G., Santos, N. C., Adibekyan, V., Delgado-Mena, E., & Israelian,G. 2015, A&A, 577, A67Southworth, J., Maxted, P. F. L., & Smalley, B. 2004, MNRAS, 351, 1277Ste ff en, J. H. 2016, MNRAS, 457, 4384Torres, G., Fischer, D. A., Sozzetti, A., et al. 2012, ApJ, 757, 161Torres, G., Fressin, F., Batalha, N. M., et al. 2011, ApJ, 727, 24Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141Weiss, L. M., & Marcy, G. W. 2014, ApJL, 783, L6Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14Weiss, L. M., Deck, K., Sinuko ff , E., et al. 2016, ArXiv e-prints,arXiv:1612.04856Winn, J. N. 2010, ArXiv e-prints, arXiv:1001.2010Wolfgang, A., & Lopez, E. 2015, ApJ, 806, 183Wright, M. H. 1996, in Pitman Research Notes in Mathematics, Vol. 344,Numerical Analysis 1995 (Proceedings of the 1995 Dundee BiennialConference in Numerical Analysis), ed. D. F. Gri ffiffi