Draftversion July15,2016 PreprinttypesetusingLATEXstyleemulateapjv.12/16/11 DISCOVERY AND VALIDATION OF A HIGH-DENSITY SUB-NEPTUNE FROM THE K2 MISSION N´estor Espinoza1,2, Rafael Brahm1,2, Andr´es Jorda´n1,2, James S. Jenkins3, Felipe Rojas1, Paula Jofr´e4, Thomas Ma¨dler4, Markus Rabus1,5, Julio Chanam´e1,2, Blake Pantoja3, Maritza G. Soto3, Katie M. Morzinski6, Jared R. Males6, Kimberly Ward-Duong7, Laird M. Close6 Draft version July 15, 2016 ABSTRACT 6 WereportthediscoveryofBD+20594b,ahighdensitysub-Neptuneexoplanet,madeusingphotometry 1 from Campaign 4 of the two-wheeled Kepler (K2) mission, ground-based radial velocity follow-up 0 from HARPS and high resolution lucky and adaptive optics imaging obtained using AstraLux and 2 MagAO,respectively. The host star is a bright(V =11.04, K =9.37), slightly metal poor ([Fe/H]= s l −0.15±0.05 dex) solar analogue located at 152.1+9.7 pc from Earth, for which we find a radius of u −7.4 J R∗ =0.928+−00..005450R⊙ and a mass of M∗ =0.961+−00..003229M⊙. A joint analysis of the K2 photometry and HARPSradialvelocitiesrevealthattheplanetisina≈42dayorbitarounditshoststar,hasaradius 4 1 of 2.23+−00..1141R⊕, and a mass of 16.3+−66..01M⊕. Although the data at hand puts the planet in the region of the mass-radius diagram where we could expect planets with a pure rock (i.e. magnesium silicate) ] composition using two-layer models (i.e., between rock/iron and rock/ice compositions), we discuss P more realistic three-layer composition models which can explain the high density of the discovered E exoplanet. The fact that the planet lies in the boundary between “possibly rocky” and “non-rocky” . exoplanets, makes it an interesting planet for future RV follow-up. h p Keywords: kepler, exoplanets - o r 1. INTRODUCTION Howard et al., 2013; Pepe et al., 2013; Grunblatt et al., t s Since the discovery of the first rocky exoplanet 2015). All of these planets have radii smaller than a (term we use here to refer to planets with masses and ∼1.6R⊕, as has been empirically determined. [ Althoughthesampleofrockyplanetsissmall,somein- radii consistent with MgSiO and Fe compositions 3 teresting relationships suggest that some of these rocky 2 following Rogers, 2015), CoRoT-7b (L´eger et al., 2009; v Queloz et al., 2009), effort has been made to find and planets mighthavecommonproperties(Weiss & Marcy, 8 study the formation, composition and evolution of these 2014). Perhaps one of the most interesting rela- 0 systems, since they resemble Earth in many ways. As tions was recently introduced by Dressing et al. (2015) 6 most rocky planets are smaller in size than 1.6R⊕, which, considering the planets with radii and mass 7 measurements measured to better than 20% precision, which correspond to masses of 6M⊕ (Weiss & Marcy, 0 show that the planets follow a common iso-composition 2014; Wolfgang & Lopez, 2015; Rogers, 2015), the . curve on the mass-radius diagram, along with Earth 1 discovery of those type of exoplanets is difficult due and Venus. This relation was recently revised by 0 to the small signals that these radii and masses imply. Zeng, Sasselov & Jacobsen (2016) to be a 74% rock and 6 In fact, in addition to CoRoT-7b, only 9 planets with 26% Fe composition. This suggests that these small, 1 secure masses and radii (i.e., masses and radii with : values more than 3−σ away from zero) in this rocky rocky analogs of Earth might have similar compositions v with small intrinsic scatter. regime exist to date: GJ1132 (Berta-Thompson et al., Xi 2015), Kepler-36b (Carter et al., 2012), K2-3d Here we report what could be a possible interesting addition to the picture of rocky worlds described above: r (Crossfield et al., 2015; Almenara et al., 2015), a Kepler-93b (Dressing et al., 2015), Kepler-10b a 2.23R⊕ exoplanet that falls just where a pure rock (i.e., magnesium silicate) composition is expected in the (Dumusque et al., 2015; Weiss et al., 2016), Kepler-23b mass radius diagram using two-layer models. Although (Ford et al., 2012; Hadden & Lithwick, 2014), Kepler- this does not mean the planet has exactly this composi- 20b (Fressin et al., 2012), Kepler-406b (Marcy et al., tion,itspositiononthediagramdoesmakesitinteresting 2014), and Kepler-78b (Sanchis-Ojeda et al., 2013; due the fact that this has been used in previous works to divide the “non-rocky” and “possibly rocky” planets 1Instituto de Astrof´ısica, Facultad de F´ısica, Pontificia Uni- (Rogers,2015). Thediscoveryismadeinthecontextofa versidadCat´olicadeChile,Av.Vicun˜aMackenna4860,782-0436 Chileanbased effortwhose aimis to follow-upplanetary Macul,Santiago,Chile 2MillenniumInstituteofAstrophysics,Av.Vicun˜aMackenna candidatesselectedusingdatafromthetwo-wheeledKe- 4860, 782-0436Macul,Santiago, Chile pler (K2) mission. K2 has proven to be very effective in 3DepartamentodeAstronom´ıa,UniversidaddeChile,Camino thesearchforexoplanets,enablingaplethoraofnewdis- alObservatorio,CerroCala´n,Santiago, Chile 4InstituteofAstronomy,UniversityofCambridge coveriesof planets of different sizes, which are especially 5MaxPlanckInstituteforAstronomy,Heidelberg,Germany interestingduetothepresenceofseveralbrighthoststars 6Steward Observatory, University of Arizona, 933 N. Cherry in the sample that allow detailed follow-up characteri- Ave,Tucson,AZ85721-0065USA sation (see, e.g., Armstrong et al., 2015; Becker et al., 7School of Earth andSpace Exploration, ArizonaState Uni- 2015; Crossfield et al., 2015; Petigura et al., 2015; versity,Tempe,AZ,85287, USA 2 Espinoza et al. 2016 Sanchis-Ojeda et al., 2015; Vanderburg et al., 2015). on the 1.2m Euler Telescope in La Silla Observatory in The paper is structured as follows. In §2 we present order to obtain rough spectral parameters of the stel- the data, which includes the K2 photometry, archival, lar host, and define whether this was a giant or a dwarf new, adaptive optics (AO) and lucky imagingof the tar- star. Data were reduced and analyzed using the proce- get star, along with high resolution spectra and radial dures decribed in Jord´an et al. (2014). The analysis of velocities obtained with the HARPS spectrograph. §3 the CORALIE spectra gave T =5600 K, log(g)= 4.4 eff presentsa jointanalysisofthe dataand presentsthe de- dex, [Fe/H]=0.0 dex and vsin(i)=2.5 km/s, which re- rived parameters of the planetary system. We discuss vealedthatthestarwasadwarfsolar-typestar. Inaddi- the results in §4 and present our conclusions in §5. tion,nosecondarypeakwasseenonthecross-correlation function indicating no detectable spectroscopic binary. 2. DATA Because of this, the target was promoted to our list of 2.1. K2 Photometry planetary candidates despite the lack of high resolution imaging needed to rule out potential blend events. K2photometryforourtargetwasobtainedbytheKe- pler spacecraft during Campaign 4. This field was ob- 2.3. High precision radial velocities with HARPS served between February and April 2015 and the data was released on September of the same year. We ob- High-precision radial velocities (RVs) were obtained tained the decorrelated versions of all the lightcurves from the HARPS spectrograph mounted on the 3.6m in the campaign which were made publicly available for telescope at La Silla between October and December of download by Vanderburg & Johnson (2014), using the 2015inordertomeasurethereflexmotionofthestardue photometrywiththeoptimalaperture,whichinthecase to the hypothetical planet producing the transit signal. of our target star corresponded to a ≈ 3 pixel radius Theobservationscoveredourpredictednegativeandpos- around the target, or an aperture of ≈ 12′′ radius. We itivequadratures,alongwithepochsinbetween,inorder performed a transit search using a Box Least Squares to probe possible long-term trends in the RVs indicative (BLS, Kov´acs, Zucker & Mazeh, 2002) algorithm. Once of a possible massive companion. 23 spectra were taken aperiodicsignalisdetectedalongwiththebest-fitdepth, intotalwiththesimultaneousThorium-Argonmode;the thetransiteventisflaggedasapontentialplanetarycan- HARPS pipeline (DRS, version 3.8) was used to reduce didateif(1)thedepthisatleast3σ largerthantheaver- these spectra and to obtain the (drift-corrected) radial agenoiselevelofthelightcurve(denotedbyσ)and(2)if velocities,whicharecalculatedviacross-correlationwith there are three or more transit events. Initially, because a G2V mask which is appropiate for the stellar type of of the last requirement, the lightcurve of the target star the host (see §3.1). The typical precision was ∼ 3 m/s was not flagged by our transit searchpipeline. However, for each individual RV measurement. For each spectra, wealsoperformedvisualinspectionofallthelightcurves, thebisectorspan,S-index,andtheintegratedfluxofthe revealing this interesting candidate. In order to dou- H and He I lines were obtained to monitor the activ- α ble check that this was indeed an astrophysical signal ity of the host star and study its influence on the RVs and not a spurious signal arising from the decorrelation (Santos et al.,2010;Jenkins et al.,2011). Themeasured method used to obtain the lightcurve, we also inspected RVs, along with these various calculated activity indica- the detrended lightcurves released by the Kepler team tors, are given in Table 1. Although the times are given using the PDC-MAP algorithm (Stumpe et al., 2012), in UTC, they were convertedto TBD (which is the time andthesamesignalwasobservedattheexactsametimes scale used by Kepler) for our joint analysis, which we as the signals observed in the Vanderburg & Johnson describe in §3.2. (2014)photometry. We werethus confidentthatthe sig- nal is of astrophysical origin and proceeded to analyse 2.4. Archival and New Imaging the light curve. A median filter with a 41 point (∼20.5 hour) window Archival imaging was obtained from the STScI Digi- tized Sky Survey8 at the EPIC coordinates of our tar- was used in order to further filter long-term variations get. DataarefromthePalomarObservatorySkySurvey of this target. The resulting median filter was smoothed (POSS). In Figure 2 we show the best images among usingaGaussianfilterwitha5-pointstandard-deviation, the available archival images in terms of the measured andthis smoothedlightcurvewasusedto normalizethe FWHM. We show images taken at two epochs and with light curve. Using this normalized lightcurve, an initial two filters: one obtained in 1995 using the RG610 filter fit using our transit-fitting pipeline (see below) revealed (red9, 590−715 nm), taken by the POSSII-F and one a P = 41.7d period for this candidate and a lightcurve whose shape resembled that of a planetary transit, with usingthe RG9filter(near-infrarred10,700−970nm)ob- a transit duration consistent with that of a planetary tained in 1996 by the POSSII-N. For reference, we show companion. Using the parameters obtained from this the aperture used to obtain our K2 photometry (black initial fit, we removed outliers from the out-of-transit circle, 12′′) along with circles with 5′′ (white solid line) data,discardinganypointsdeviatingmorethan3-σfrom and 2′′ (white dashed line) radii which are centered on themedianflux. Theresultingnormalizedversionofthis the centroid of our target star, which was obtained by lightcurve is shown on Figure 1. No other significant fitting a 2D Gaussian to the intensity profile. signals were found in the photometry. 8 http://stdatu.stsci.edu/cgi-bin/dss_form 2.2. Reconnaissance spectroscopy 9 Transmission curve available at http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/dss/TransmissionCurves/PO Ahighresolutionspectrumofthistargetwastakenon 10 Transmission curve available at October21stwiththe CORALIEspectrographmounted http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/dss/TransmissionCurves/PO Discovery and Validation of a High-Density sub-Neptune from the K2 Mission 3 1.004 4000 N Normalizedflux 1.0012 02000 ormalizedflux(pp m ) 0.998 −2000 7070 7080 7090 7100 7110 7120 7130 1.00025 250 1 0 N Normalizedflux 00.9.99997955 −−520500 ormalizedflux(ppm ) 0.99925 −750 7070 7080 7090 7100 7110 7120 7130 Time(BJD-2450000) Figure 1. K2photometry(obtained fromVanderburg&Johnson,2014,upperpanel)andlong-termandoutliercorrectedversionofthe photometry(lowerpanel). Thesmooth,long-termvariationobservedintheoriginalphotometrywasremovedbyasmoothedmedian-filter, depictedintheupperpanelbyaredsolidline,whichwasusedforoutlierremoval(seetext). Twocleartransit-likeevents canbeseenon bothversionsofthephotometrycloseto2457070and2457110 BJD(indicatedwithredarrows). Notethattheprecisionobtainedforthis lightcurveis∼55ppm(rms)perpoint. Table 1 RadialvelocitiesobtainedwiththeHARPSspectrographalongwithvariousactivityindicators. BJD RV σRV BIS σBIS SH,K σSH,K Hα σHα HeI σHα (UTC) msec−1 msec−1 msec−1 msec−1 dex dex dex dex dex dex 2457329.63450 −20333.9 4.4 35.0 6.2 0.1748 0.0063 0.10151 0.00013 0.50230 0.00081 2457329.67362 −20340.1 3.6 28.9 5.0 0.1535 0.0050 0.10337 0.00013 0.50279 0.00081 2457329.72375 −20337.9 3.9 34.2 5.6 0.1864 0.0053 0.10367 0.00013 0.50146 0.00081 2457330.80181 −20343.1 2.6 22.6 3.7 0.1483 0.0035 0.10167 0.00013 0.50787 0.00082 2457331.63418 −20342.5 2.4 23.2 3.4 0.1551 0.0029 0.10171 0.00013 0.51233 0.00083 2457331.68695 −20338.4 2.0 11.0 2.8 0.1573 0.0026 0.10209 0.00013 0.50301 0.00081 2457332.64705 −20335.7 2.6 20.2 3.7 0.1549 0.0038 0.10236 0.00013 0.50221 0.00081 2457332.72713 −20338.4 2.0 12.9 2.9 0.1459 0.0025 0.10178 0.00013 0.50273 0.00081 2457336.65528 −20345.3 3.2 20.8 4.5 0.1701 0.0042 0.10397 0.00013 0.49984 0.00081 2457336.73328 −20339.8 3.6 9.5 5.1 0.1547 0.0047 0.10187 0.00013 0.49884 0.00081 2457339.70924 −20343.1 4.9 14.2 6.9 0.1985 0.0068 0.09939 0.00013 0.49982 0.00081 2457339.72063 −20339.2 4.1 31.1 5.8 0.1738 0.0058 0.10496 0.00013 0.50575 0.00082 2457340.69354 −20340.0 3.7 5.9 5.3 0.1833 0.0056 0.10217 0.00013 0.51486 0.00083 2457340.70475 −20336.6 3.1 23.9 4.3 0.1687 0.0044 0.10083 0.00013 0.50038 0.00081 2457341.74523 −20336.6 3.3 18.7 4.7 0.1658 0.0043 0.10536 0.00013 0.50213 0.00081 2457341.75600 −20335.3 3.2 25.4 4.5 0.1491 0.0043 0.10274 0.00013 0.50658 0.00082 2457348.80101 −20330.3 3.1 21.6 4.4 0.1858 0.0056 0.10339 0.00013 0.50277 0.00081 2457360.62435 −20337.2 2.2 14.8 3.2 0.1581 0.0029 0.10354 0.00013 0.50080 0.00081 2457360.63915 −20336.8 2.1 11.0 3.0 0.1585 0.0027 0.10237 0.00013 0.50273 0.00081 2457361.66418 −20337.3 3.1 30.7 4.4 0.1719 0.0042 0.10607 0.00013 0.50020 0.00081 2457361.67814 −20337.0 2.9 32.8 4.1 0.1533 0.0038 0.10231 0.00013 0.49660 0.00080 2457362.66191 −20331.0 2.2 15.4 3.1 0.1517 0.0031 0.10313 0.00013 0.49866 0.00081 2457362.67602 −20335.6 2.2 13.5 3.1 0.1575 0.0032 0.10479 0.00013 0.50121 0.00081 4 Espinoza et al. 2016 1995-09-04(POSSII-F+RG610) 1996-09-11(POSSII-N+RG9) 60 0.25 60 0.25 (a) (b) 40 40 0.2 0.2 20 20 ec) 0.15R ec) 0.15R C(arcs 0 elative C(arcs 0 elative DE 0.1 flu DE 0.1 flu x x ∆−20 ∆−20 0.05 0.05 −40 −40 −60 0 −60 0 −60 −40 −20 0 20 40 60 −60 −40 −20 0 20 40 60 ∆R.A.(arcsec) ∆R.A.(arcsec) Figure 2. Archival imaging for our target at the coordinates given in the EPIC catalog obtained with different versions of POSS: (a) POSSII-F survey, taken with a red filter and (b) POSSII-N survey, taken with an infrared filter. The black circle indicates the aperture used for our K2 data. The white solid circle has a radius of 5′′ and the dashed circle a radius of 2′′ for illustration purposes; these are centeredonthemeasuredcentroidofthetargetstar. Theredcircletotheleftofthetargetstarmarksanobjectwhichis∼8.2magnitudes fainterthanthetargetinR(seetext). Discovery and Validation of a High-Density sub-Neptune from the K2 Mission 5 60 2015-12-27(LCOGT-SBIG-CTIO+R) −2.4 out stars within a 2′′ radius. We follow methods similar to those described in Morzinski et al. (2015) to reduce ourimages,whichwebrieflydescribe here;aPythonim- 40 −2.42 plementation of such methods is available at Github11. First,the imageswerecorrectedby darkcurrentbut not R C(arcsec) 200 −2.44elativeflux flasesantasifirteeivslduitelitdes,ofb(oseepceatiucssaeecltdtihoisnetoflArat.it3osnissnhaonwMdoanrnoztiunonsfkeivineetntrianfllsu.i,xc2pl0ei1xv5ee,ll ∆DE−20 −2.46(log-scale) fmoorradseakr dtpoertoamviiladesdekdebxbapydlaMnpaioxtriezolisnn.sokAfi eftttheiarsl.tehff(e2es0cet1)5c.)orAwreacbstaiuodnsespdiaxirenel applied to each image, we obtain a median image using −2.48 −40 our 32 frames in order to get an estimate of the back- ground flux, which we then subtract from each of the −60 −2.5 individual frames. In order to further correct for differ- −60 −40 −20 0 20 40 60 ences in the sky backgrounds of each image, we apply a ∆R.A.(arcsec) 2D median filter with a 200 pixel (≈ 3′′) window which Figure 3. Modernimagingof ourtarget obtained fromLOCGT takes care of large-scale fluctuations of each image. The fromCTIOusingtheSBIGcamerawiththeRfilteron2015/12/27. Notethat,althoughthecircleshavethesamemeaningsastheones background-subtracted images are then merged by first inFigure2,thescalehereisdifferent. rotating them to the true north (using the astrometric calibrationdescribedinMorzinski et al.,2015)andcom- NewimagingwasobtainedusingtheLasCumbresOb- bined using the centroid of our target star (obtained by servatoryGlobalTelescopeNetwork(LCOGT).Fourim- fitting a 2D Gaussian to the profile) as a common refer- ages were taken using the SBIG camera with the Bessel encepointbetweentheimages. OurresultingAOimage, R filter on UT 2015/12/27from the Cerro Tololo Inter- obtained by combining our 32 images, is shown in Fig- american Observatory (CTIO). Our target star reached ure4. A2DgaussianfittothetargetstargivesaFWHM close-to saturation counts (∼ 47000 counts) in order to of 0′.′2, which we set as our resolution limit. haveenoughphotonstoobservetheclose-bystarspresent The limiting contrasts in our AO observations in the in the POSSimages. Figure 3 shows the resulting image K band were estimated as follows. First, a 2D gaussian s obtained by median-combining our four images, along fittothetargetstarwasmadeandusedtoremoveitfrom with the same circles as those drawn on Figure 2. the image. Although a 2D gaussian does not perform a Giventhatthelargestpotentialsourceoffalse-positive perfect fit at the center, the fit is good enough for the detections in our case comes from blended eclipsing bi- wingsofthePSF,whichisouraim. Then,ateachradial nary systems mimicking a planetary transit event, we distancen×FWHMawayfromthetargetstar,wheren= note that, given that the depth of the observed tran- 1,2,...,15 is an integer, a fake source was injected at 15 sit is ∼ 0.05%, if a blended eclipsing binary system was differentangles. Sourceswithmagnitudedifferencesfrom responsible of the observed depth, then assuming a to- 11 to 0 were injected in 0.1 steps, and a detection was tal eclipse of the primary (which is the worst case sce- definedif3ormorepixelswere5-sigmaabovethemedian nario;all other scenariosshould be easierto detect), the flux level at that position. The results of our injection eclipsedstarwouldhavetobe∼8.23magnitudesfainter and recovery experiments are plotted in Figure 5. thanourtargetstarintheKeplerbandpass. Wecancon- Only one source was detected at ∼ 2′′ from the tar- fidently rule out such a bright star down to a distance get. Theshapeandpositionofthisobjectisinconsistent of 9′′ of the target star with the POSSII and LCOGT withaspecklebutisveryfaint: wemeasureamagnitude images. For reference, the closest star to the left of the difference of ∆K =9.8 with the target and is thus just s targetstar (indicated with a red circle) in Figures 2 and above our contrast level at that position (see Figure 4, 3 is ∼8.2 magnitudes fainter than the target star in the the source is indicated with a grey circle in the upper Rband. Ascanbeseenontheimages,astarthatbright right). A careful assessment of the PSF shape, however, wouldbeevidentinthearchivalPOSSimagesand/oron made it inconsistent with the object having the same our new LCOGT images at distances larger than 9′′. PSF shape as our star. Comparing its PSF with known “ghosts” on the image, on the other hand, revealed that 2.5. Adaptive optics & lucky imaging this source is not of astrophysical but of instrumental Adaptiveoptics(AO)imagingwasobtainedusingMa- origin. gAO+Clio2 instrument mounted at the Magellan Clay Inordertosearchforcompanionsatlargerseparations, telescope in Las Campanas Observatory on December luckyimagingwasobtainedwithAstraLuxSurmounted 6th using the K filter with the full Clio2 1024× 512 ontheNewTechnologyTelescope(NTT)atLaSillaOb- s pixel frames of the narrow camera (f/37.7). The natural servatory(Hippler et al., 2009) on 2015/12/24using the guide star system was used and, because our target is i′ band. Figure 4 shows our final image obtained by relatively bright, it was used as the guide star. 32 im- combiningthebest10%imageswithadrizzlealgorithm. ages with exposure times of 30 sec each were taken in Because the PSF shape obtained for our lucky imaging five different positions of the camera (nodding), all of iscomplexandwealreadyruledoutcompanionsinsidea them at different rotator offset angles. Due to a motor 2′′ radius with Magellan+Clio2, and given that our ob- failure of the instrument, the nodding and rotation pat- jective with lucky imaging was to rule out companions ternswerenotabletocoverthefull16′′×8′′ fieldofview aroundthestar. However,itgaveusenoughdatatorule 11 https://github.com/nespinoza/ao-reduction 6 Espinoza et al. 2016 2015-12-06(MagAO+Ks) 2015-12-24(AstraLuxSur+i’) 0 0 2 −1 4 −1 −2 1 R 2 −2R e e C(arcsec) 0 −−43lativeflux C(arcsec) 0 −3lativeflux ∆DE −5(log-sc ∆DE−2 −4(log-sc −1 ale ale −6) −5) −4 −2 −7 −6 −8 −6 −2 −1 0 1 2 −6 −4 −2 0 2 4 ∆R.A.(arcsec) ∆R.A.(arcsec) Figure 4. (Left)Adaptiveopticsimage(log-scale)obtainedwithMagAO+Clio2on2015-12-06. Theblackdashedcirclehasa2′′ radius forillustrationandcomparisonwithFigure2;thegreycirclemarkafaintsourcefoundonourimage,whichwasaboveourcontrastlimit but we identified as being of instrumental origin (see text). (Right) AstraLux Sur i′-band observations of our candidate on 2015-12-24. Theinnerblackdashedcircleindicates2′′,whiletheouterblacksolidcircleindicates5′′ forcomparisonwithFigure2. 0 0 3.1. Stellar properties In order to obtain the properties of the host star, we 2 2 made use of both photometric andspectroscopic observ- ablesofourtarget. Forthe former,weretrievedB,V,g,r and i photometric magnitudes from the AAVSO Pho- 4 4 tometric All-Sky Survey (APASS, Henden & Munari, Ks ∆ 2014) and J, H and K photometric magnitudes from ∆ 6 6 i 2MASS for our analysis. For the spectroscopic observ- ables, we used the Zonal Atmospherical Stellar Param- eter Estimator (ZASPE, Brahm et al., 2016) algorithm 8 8 using our HARPS spectra as input. ZASPE estimates the atmospheric stellar parameters and vsini from our highresolutionechellespectraviaaleastsquaresmethod 10 10 against a grid of synthetic spectra in the most sensitive 0 1 2 3 4 5 zones of the spectra to changes in the atmospheric pa- Radialdistance(arcsec) rameters. ZASPE obtains reliable errors in the parame- Figure 5. 5-σ contrast curves obtained fromour MagAO+Clio2 ters, as well as the correlationsbetween them by assum- Ksband(blackline)andAstraLuxSuri′-band(redline)observa- ing that the principal source of error is the systematic tionsofourcandidate. mismatch between the data and the optimal synthetic atlargerangulardistances,wedidnotperformPSFsub- spectra, which arises from the imperfect modelling of stractionalgorithmsinordertoobtainthe5−σcontrasts the stellar atmosphere or from poorly determined pa- atthosedistances. Instead,weusedsimpleaperturepho- rameters of the atomic transitions. We used a synthetic tometry in order to estimate the 5−σ contrasts outside the 2′′ radius by performing a procedure similar to that grid provided by Brahm et al. (2016) and the spectral described in W¨ollert et al. (2015). In summary, we esti- region considered for the analysis was from 5000 ˚A to matedthenoiselevelina5×5boxateachradialdistance 6000 ˚A, which includes a large number of atomic transi- at 15 different angles for distances larger than 2′′ from tions and the pressure sensitive Mg Ib lines. The result- the estimated centroid of the image (where the contri- ing atmospheric parameters obtained through this pro- bution ofthe targetstar’PSFto the backgroundlevelis cedure were Teff = 5766± 99 K, log(g) = 4.5± 0.08, low),andcalculatedthemagnitudecontrastbyobtaining [Fe/H] = −0.15±0.05 and vsin(i) = 3.3±0.31 km/s. the flux of the target star using a 5-pixel radius around With these spectroscopic parameters at hand and the it and a 5-pixel radius about the desired distance from photometric properties, we made use of the Dartmouth the star, where 5−σ counts are summed to each pixel Stellar Evolution Database (Dotter et al., 2008) to ob- atthatdistancebeforeperformingtheaperturephotom- tain the radius, mass, age and distance to the host etry. Then, the magnitude contrast at a given distance star using isochrone fitting with the isochrones pack- isobtainedasthe averagevalueobtainedatthedifferent age (Morton et al., 2015). We take into account the angles. Theresulting5−σcontrastsarepresentedinFig- uncertainties in the photometric and spectroscopic ob- ure 5. We study the constrains that our archival, new, servables to estimate the stellar properties, using the AO and lucky imaging put on the false-positive proba- emcee(Foreman-Mackey et al., 2013) implementation of bilities and transit dilutions on the next section. theaffineinvariantMarkovChainMonteCarlo(MCMC) ensemblesamplerproposedinGoodman & Weare(2010) 3. ANALYSIS in order to explore the posterior parameter space. We Discovery and Validation of a High-Density sub-Neptune from the K2 Mission 7 Table 2 StellarparametersofBD+20594. Parameter Value Source IdentifyingInformation EPICID 210848071 EPIC 2MASSID 03343623+2035574 2MASS R.A.(J2000, h:m:s) 03h34m36.23s EPIC DEC(J2000, d:m:s) 20o35′57.23′′ EPIC R.A.p.m. (mas/yr) 36.7±0.7 UCAC4 DECp.m. (mas/yr) −51.8±1.3 UCAC4 Spectroscopic properties Teff (K) 5766±99 ZASPE SpectralType G ZASPE [Fe/H](dex) −0.15±0.05 ZASPE logg∗ (cgs) 4.5±0.08 ZASPE vsin(i)(km/s) 3.3±0.31 ZASPE Photometricproperties Kp (mag) 11.04 EPIC B (mag) 11.728±0.044 APASS V (mag) 11.038±0.047 APASS g′ (mag) 11.352±0.039 APASS r′ (mag) 11.872±0.050 APASS i′ (mag) 10.918±0.540 APASS J (mag) 9.770±0.022 2MASS H (mag) 9.432±0.022 2MASS Ks(mag) 9.368±0.018 2MASS Derivedproperties M∗ (M⊙) 0.961+−00..003229 isochrones+ZASPE R∗ (R⊙) 0.928+−00..005450 isochrones+ZASPE ρ∗ (g/cm3) 1.70+−00..2206 isochrones+ZASPE L∗ (L⊙) 0.88+−00..1152 isochrones+ZASPE Distance(pc) 152.1+9.7 isochrones+ZASPE −7.4 Age(Gyr) 3.34+1.95 isochrones+ZASPE −1.49 Note. Logarithmsgiveninbase10. obtain a radius of R∗ = 0.928+−00..005450R⊙, mass M∗ = (2013). For the quadratic and square-root laws, we use 0.961+−00..003229M⊙, age of 3.3+−11..95 Gyr and a distance to the tdheerttorasnasmfoprlmeathtieonpshydseisccarlliybepdlaiunsiKblieppvainluges(2o0f1t3h)eilnimobr-- host star of 152.1+9.7 pc. The distance to the star was −7.4 darkeningcoefficients. Forthelogarithmiclawweusethe also estimated using the spectroscopic twin method de- transformations described in Espinoza & Jord´an(2016), scribedinJofr´e et al.(2015),whichisindependentofany which presents the sampling of the limb darkening pa- stellar models. The values obtained were 158.3±5.4 pc rameters for the more usual form of the logarithmic law when using 2MASS J band photometry and 160.0±5.7 to allow for easier comparison with theoretical tables (if pc if H band photometry was used instead, where the the geometry of the system is properly taken into ac- stars HIP 1954, HIP 36512, HIP 49728 and HIP 58950 count, see Espinoza & Jord´an, 2015). The code also al- were used as reference for the parallax. Those values lowsthe usertofit the lightcurveassumingeither apure areinverygoodagreementwiththevalueobtainedfrom white-noise model or an underlying flicker (1/f) noise isochronefitting. The stellar parametersofthe hoststar pluswhite-noisemodelusingthewavelet-basedtechnique are sumarized in Table 2. described in Carter & Winn (2009). For the RV mod- 3.2. Joint analysis elling, exonailer assumes Gaussian uncertainties and addsajitterterminquadraturetothem. Thejointanal- We performed a joint analysis of the photometry and ysisisthenperformedusingthe emceeMCMCensemble the radial velocities using the EXOplanet traNsits and sampler (Foreman-Mackey et al., 2013). rAdIalveLocityfittER,exonailer,whichismadepub- For the joint modelling of the dataset presented here, licly available at Github12. For the transit modeling, we tried both eccentric and circular fits. For the radial exonailer makes use of the batman code (Kreidberg, velocities, uninformative priors were set on the semi- 2015), which allows the user to use different limb- amplitude, K, and the RV zero point, µ. The former darkening laws in an easy and efficient way. If chosento was centered on zero, while the latter was centered on be free parameters, the sampling of the limb-darkening the observed mean of the RV dataset. Note that our coefficients is performed in an informative way using priors allow us to explore negative radial velocity am- the triangular sampling technique described in Kipping plitudes, which is intentional as we want to explore the possibility of the RVs being consistent with a flat line 12 http://www.github.com/nespinoza/exonailer 8 Espinoza et al. 2016 (i.e., K = 0). Initially a jitter term was added but we analyze the likelihood function of the radial-velocity was fully consistent with zero, so we fixed it to zero datainorderto comparethe modelsanddecide whichis in our analysis. As for the non-circular solutions, flat preferred by the data. We obtain that both models are priors were set on e and on ω instead of fitting for the indistinguishable, with both the AIC (∆AIC = 2) and Laplace parameters ecos(ω) and esin(ω) because these BIC (∆BIC = 2) values being ∼ 2. We thus choose the imply implicit priors on the parameters that we want to simpler model of those two, which is the circular model, avoid (Anglada-Escud´e et al., 2013). For the lightcurve and report the final parameters using this as our final modelling,weusedtheselectiveresamplingtechniquede- model. scribed in Kipping (2010) in order to account for the 30 The resulting parameters of our fit are tabulated in min cadence of the K2 photometry, which has as a con- Table 3. It is interesting to note that the radial veloc- sequence the smearing of the transit shape. In order to ity semi-amplitude is inconsistent with zero by almost minimize the biases in the retrieved transit parameters 3σ. Moreover, we are confident that those variations we fit for the limb darkening coefficients in our analy- do not arise from activity as all the correlation coeffi- sis (see Espinoza & Jord´an, 2015). In order to decide cients we calculate between our RVs and the different which limb-darkening law to use, we apply the method activity indexes given in Table 1 give correlation coeffi- described in Espinoza & Jord´an (2016) which, through cients which are consistent with 0 at ≈1σ, and all vari- simulations and given the lightcurves properties, aids in ations of the activity indices at the period and time of selecting the best limb-darkening law in terms of both transit-centerfoundforourtargetareconsistentwithflat precision and bias using a mean-squared error (MSE) lines. Interestingly,the radial-velocitysemi-amplitude is approach. In this case, the law that provides the mini- large for a planetary radius of only Rp = 2.23+−00..1141R⊕; mum MSE is the quadratic law, and we use this law in the K = 3.1+1.1 m/s semi-amplitude implies a mass of order to parametrize the limb-darkening effect. In addi- −1.1 tion,theK2photometryisnotgoodenoughtoconstrain Mp = 16.3+−66..01M⊕, which at face value could be con- sistent with a rocky composition, a rare property for a the ingress and egress times because only two transits Neptune-sized exoplanet such as BD+20594b. We cau- wereobservedinlong-cadencemode,whichprovidespoor tion, however, that this interpretation has to be taken phase coverage; this implies that the errors on a/R∗ are with care, as we have poor phase coverage on the “up” rather large. Because of this, we took advantage of the quadrature. We put these values in the context of dis- stellar parameters obtained with our HARPS spectra, coveredexoplanets of similar size in §4. and derived a value for this parameter from them (see Sozzetti et al., 2007) of a/R∗ = 54.83+−23..1196. This value 3.3. Planet scenario validation was used as a prior in our joint analysis in the form of a In order to validate the planet scenario which we have Gaussian prior. We used the largest of the errorbars as implied in the past sub-section, we make use of the for- the standard deviation of the distribution, which is cen- teredonthequotedmedianvalueoftheparameter13. We malism described in Morton (2012) as implemented on tried both fitting a flicker-noisemodel and a white-noise the publicly available vespa14 package. In short, vespa model, but the flicker noise model parameters were con- considers all the false-positive scenarios that might give sistentwithno1/f noisecomponent,sothefitwasfinally rise to the observed periodic dips in the light curve obtainedassumingwhitenoise. 500walkers wereusedto and, using photometric and spectroscopic information of the target star, calculates the false-positive probabil- evolvetheMCMC,andeachoneexploredthe parameter ity (FPP) which is the complement of the probability space in 2000 links, 1500 of which were used as burn-in samples. This gave a total of 500 links sampled from of there being a planet given the observed signal. Be- the posterior per walker, giving a total of 250000 sam- causeourarchivalandmodernimagingpresentedon§2.4 ples fromthe posteriordistribution. Thesesampleswere rule out any companion at distances larger than 9′′ ra- tested to converge both visually and using the Geweke dius, we consider this radius in our search for possible (1992) convergence test. false-positive scenarios using vespa,which considers the Figures 6 and 7 show close-ups to the phased photom- area around the target star in which one might suspect etry and radial velocities, respectively, along with the false-positives could arise. The algorithm calculates the best-fitmodelsforbothcircular(red,solidline)andnon- desired probability as circular (red, dashed line) fits obtained from our joint analysisofthedataset. Thelightcurvefitsforbothmod- 1 FPP= , elsareverysimilar,butintheRVsthedifferencesareevi- 1+f P p dent. Inparticular,theeccentricfitgivesrisetoaslightly smaller semi-amplitude than (yet, consistent with) the where fp is the occurrence rate of the observed planet one obtained with the circular fit. For the eccentric fit, (at the specific observed radius) and P = LTP/LFP, we obtain e = 0.096+0.089, ω = 53+17 degs and a semi- where TP indicates the transiting-planet scenario and amplitude of K = −2.09.0+616.1 m sec−−12.3 For the circular FP the false-positive scenario, and each term is defined −1.0 as L =π L , where π is the prior probability and L is orbit, we find a semi-amplitude of K =3.1+−11..11 m sec−1. the liikelihiooid of the ii-th scenario. For our target, coin- Since the differences on the lightcurves are very small, sideringalltheinformationgatheredandthefactthatno secondary eclipse larger than ≈ 165 ppm (i.e., 3-sigma) 13Performingajointanalysiswithalargeuniformpriorona/R∗ is detected, we obtain a value of P = 4288.79. As for spanningfroma/R∗∈(25,70)givesaposteriorestimateofa/R∗= the occurrence rate of planets like the one observed, we 55.92+5.64 forthisparameter,whichisinexcellentagreementwith −13.11 thisspectroscopicallyderivedvalue. 14 https://github.com/timothydmorton/VESPA Discovery and Validation of a High-Density sub-Neptune from the K2 Mission 9 1.0005 500 N o x 1 0 rm u dfl aliz ze ed mali flu x r No 0.9995 −500 (p p m ) 0.999 −1000 −10 −5 0 5 10 HoursfromMid-Transit Figure 6. Phase-folded photometry (grey points; circles the first transit, triangles the second transit) and best-fit transit lightcurve for thecircular(red,solidline)andeccentric(red,dashedline)fitsforourplanetobtainedfromourjointanalysis. Notethatthedifferencein thelightcurveforbothfitsisverysmall. 2200 Table 3 OrbitalandplanetaryparametersforBD+20594. Parameter Prior PosteriorValue 1100 Lightcurveparameters P (days)................ N(41.68,0.1) 41.6855+0.0030 (m/s)(m/s) T0−2450000(BJDTDB) N(7151.90,0.1) 7151.90−210+−.0000..3001004427 RVRV 00 a/R⋆ ................... N(54.83,3.16) 55.8+−33..33 Rp/R⋆.................. U(0,0.1) 0.02204+−00..0000005587 i(deg).................. U(80,90) 89.55+0.17 −−1100 q1 ...................... U(0,1) 0.38+−−00..021.9614 q2 ...................... U(0,1) 0.52+−00..3320 −−11 −−00..55 PPhh00aassee 00..55 11 σw (ppm) .............. J(50,80) 55.00+−00..7732 RVparameters Figure 7. Phase-foldedHARPSradialvelocities(grey)andbest- fit radial velocity models for both circular (red, solid line) and K (ms−1).............. N(0,100) 3.1+−11..11 eccentric (red,dashed line)fits usingourjointanalysis. Thelight µ(kms−1).............. N(−20.337,0.1) −20.33638+0.00073 bluebands indicateregionsthathavebeenrepeated forbetter vi- −0.00073 sualisationoftheRVcurve. e ....................... — 0(fixed) DerivedParameters consider the rates found by Petigura, Marcy & Howard (2013) for planets between 2−2.83R⊕ with periods be- Mp (M⊕) .............. — 16.3+−66..01 tween 5 and 50 days orbiting solar-type stars, which is Rp (R⊕) ............... — 2.23+−00..1141 b7.i8li%ty,oi.fe.F,PfpP==03.0×781.0T−3h.isGgiivveensuthsaatftahlisse-parlaorbmabpilritoybais- ρp (g/cm3) ............. — 7.89+−33..41 smaller thanthe usual 1%threshold(e.g., Montet et al., loggp (cgs) ............. — 3.50+−00..1241 2015), we consider our planet validated. We note that a(AU) ................. — 0.241+−00..001197 thisFPPisanupperlimitontherealFPPgivenourAO Vesc (km/s) ............ — 30.2+−56..32 and lucky-imaging observations. Both observations rule Teq (K) ................ out an important part of the parameter space for blend- ingscenariosbetween0′.′2and5′′fromthestar,whichare Bondalbedoof0.0 — 546+−1198 the main source of false-positives for our observations. Bondalbedoof0.75 — 386+−1132 Note. Logarithms given inbase10. N(µ,σ) stands for anormal prior 3.4. Transit dilutions with mean µ and standard-deviation σ, U(a,b) stands for a uniform prior with limits a and b and J(a,b) stands for a Jeffrey’s prior with As it will be discussed in the next section, both the thesamelimits. planet radius and mass puts BD+20594 in a very inter- esting part of the mass-radius diagram. Therefore, it changedby a collection of stars inside the aperture used is important to discuss the constraints that our spectro- to obtain the photometry of the target star is given by scopicandnew,AOandluckyimagingobservationspose p1/F , where F is the fractionof the total flux in the % % onpossible backgroundstarsthat mightdilute the tran- aperture added by the star being transited, we estimate sit depth and thus cause us to underestimate the transit that only stars with magnitude differences . 2 are able radius. to change the transit radius by magnitudes similar to Given that the factor by which the planetary radius is the quoted uncertainties in Table 3. We note that such 10 Espinoza et al. 2016 magnitude differences in the Kepler bandpass are ruled to be rocky would result in p ∼ 0. The definition rocky out from 0′.′2 to the aperture radius used to obtain the of“rockyplanet”usedinRogers(2015),whichweadopt photometry for our targetstar: our AO and lucky imag- in this work, is given by those planets spanning com- ingobservationsruleoutcompanionsofsuchmagnitudes positions between 100% rock and 100% Fe. Although from 0′.′2 to 5′′ (see Figure 5). On the other hand, stars this definition is based on simple 2-layer models for the with magnitude differences of that order should be evi- planetarycomposition,andintheoryforagivenpointin dentonour retrievedarchivalandnew images presented the mass-radiusdiagramplanets couldhavedenser com- in §2.4, at least at distances of 5′′ from our target star, positions with a gaseous envelope on top, we use this anduptoandbeyondthe12′′apertureusedtoobtainthe metric anyways in order to compare our newly discov- K2 photometry. Given that the remaining unexplored eredexoplanetin terms ofthe populationofalreadydis- areaonthe skyisverysmall(only0′.′2aroundourtarget covered small planets. This is an important point to star), and that a star of such magnitude should produce make, as p is actually an upper limit on the prob- rocky an evident peak on the cross-correlationfunction on our ability that a planet is indeed rocky. To compute this high resolution spectra which is not seen, we consider value and compare it to the population of exoplanets that our derived transit radius is confidently unaffected with secure masses and radii discovered so far, we use by dilutions of background field stars. the models from Zeng, Sasselov & Jacobsen (2016). To 4. DISCUSSION sample from the posterior distributions given the poste- riorestimatespublishedintheliteratureforthedifferent As mentioned in the previous section, the large mass exoplanets,weusethemethodsdescribedinAppendixA (Mp = 16.3+−66..01M⊕) for the calculated radius (Rp = ofEspinoza & Jord´an(2015)andassumetheseradiiand 2.23+−00..1141R⊕) found for BD+20594b is very interesting. masses are drawn from skew-normal distributions in or- Figure 8 compares BD+20594b with other discovered der to use the asymmetric errorbars published for those exoplanets with radii less than 4R⊕ (∼ Neptune) and parameters, while we use the posterior samples of our massessmallerthan32M⊕ (limits oftheoreticalmodels) MCMC fits described in §3.2 to sample from the poste- as retrieved from exoplanets.eu15 except for the Kepler- riorjointdistributionofmassandradiusofBD+20594b. 10 planets, for which we use the masses obtained by Our results are depicted in Figure 9, where we also in- Weiss et al. (2016), along with 2-layer models obtained dicate the threshold radius found by Rogers (2015) at from Zeng, Sasselov & Jacobsen(2016). As can be seen, whichthereis a significanttransitionbetweenrockyand BD+20594b spans a regime in radius at which most ex- non-rocky exoplanets, with smaller exoplanets having in oplanets have low densities and are composed of large generalrockycompositionsandlargerexoplanetshaving amounts of volatiles (Rogers, 2015). In particular, tak- less dense compositions. ing the mass-radius estimates for BD+20594b at face As evident in Figure 9, BD+20594b is in an interest- value, the best-fit compositionassuming a 2-layermodel ing position in this diagram. The closest exoplanet to for the planet is 100% MgSiO , i.e., a pure rock compo- BD+20594b in this diagram is Kepler-20b, which has a 3 sition, positioning the planet in the boundary of “pos- radius of 1.91+−00..1221R⊕, which is only 2 − σ away from sibly rocky” and “non-rocky” planets. More realistic the “rocky” boundary. BD+20594b, on the other hand, three-layer alternatives, however, can explain the ob- is more than 5 − σ away from it. With a value of served radius and mass of the planet if a rock/Fe core p ∼ 0.43, BD+20594b is the first Neptune-sized rocky has an added volatile envelope, composed either by wa- exoplanet to date with a large (compared to the typi- ter or H/He (see, e.g., the modelling for Kepler-10c in cal Neptune-sized planet) posterior probability of being Weiss et al.,2016). If,forexample,weassumeanEarth- dense enough to be rocky. likeinteriorcompositionfortheplanet(i.e.,74%MgSiO 3 ThelargemassobtainedforBD+20594bimpliesthatif and 26% Fe) and again take the mass and radius esti- theplaneteverhadthechancetoacquireanatmosphere, mates at face value, three-layer models obtained from it should retain it. However, if the planet is indeed ac- Zeng, Sasselov & Jacobsen (2016) give a possible 0.2R⊕ tually primarlycomposedof rock,givenits smallradius, water envelope for the planet (corresponding to 8% in a significant H/He envelope is unlikely in the usual set- mass). This thus gives a maximum radius for a possible tings of planet formation. Calculations using core ac- H/He envelope, which would anyways produce a small cretion theory by Ikoma & Hori (2012), predict that if layerofmuchlessthanapercentinmass;atleastsignif- the mass of rock in the protoplanet is on the order of icantly smaller than the one modelled for Kepler-10c. ∼10M⊕, even for disk dissipation time-scales on the or- Given that the errors on the mass of BD+20594b der of ∼ 10kyr an accretion of a ∼1M⊕ H/He envelope are large enough to be consistent with several compo- should happen. Even in the case of a large opacity of sitions, a careful assessment must be made in order to the protoplanetary disk, a mass of rock similar to the explore its possible rocky nature. To this end, we follow one possible for BD+20594b should imply at least this the approach introduced by Rogers (2015) and compute levelofH/Heaccretion. Giventhebulkcompositionand p , the posterior probability that a planet is suffi- rocky distance of BD+20594bto its parent star, mass loss due ciently dense to be rocky, which is defined as the frac- to X-ray and Extreme UV radiationfrom its parentstar tion of the joint mass-radius posterior distribution that its unlikely. If this indeed is the primary composition of falls between a planet composition consistent with be- this planet, it might be possible that it formed at late ing rocky. A probably rocky planet, then, would have stages in the protoplanetary disk, under conditions sim- procky ∼1,while a planetwitha densitythatis too low ilarto those ontransitiondisks (Lee & Chiang, 2016) or that some external effect removedthe accreted envelope 15 Dataretrievedon23/12/2015