CHARACTERIZATION OF THE K2-19 MULTIPLE-TRANSITING PLANETARY SYSTEM VIA HIGH-DISPERSION SPECTROSCOPY, AO IMAGING, AND TRANSIT TIMING VARIATIONS

In order to improve on estimates of the stellar parameters, we observed K2-19 with the High Dispersion Spectrograph (HDS) on the Subaru 8.2 m telescope on 2015 May 29 (UT). To maximize the S/N as well as to achieve a sufficient spectral resolution, we employed the Image Slicer #2 (Tajitsu et al. 2012, R ∼ 80,000) and the standard I2a setup, simultaneously covering the spectral range between 4950 and 7550 Å with the two CCD chips. The raw spectrum was subjected to the standard IRAF¹⁴ procedures to extract the one-dimensional (1D) spectrum. The wavelength scale was established by observations of the comparison thorium–argon lamp taken during evening and morning twilight. The exposure time was 20 minutes. The 1D spectrum has a S/N of ∼70 per pixel in the vicinity of the sodium D lines.

We observed K2-19 in the H band with the HiCIAO (Tamura et al. 2006) combined with the AO188 (188 element curvature sensor adaptive-optics system: Hayano et al. 2008), mounted on the 8.2 m Subaru Telescape on 2015 May 8 (UT). We used the target star itself as a natural guide star for AO188, and employed an atmospheric dispersion corrector (ADC) to prevent the star from drifting on the detector due to differential refraction between the visible and near-infrared bands. The field of view (FOV) was 20'' × 20'', and the typical AO-corrected seeing was ∼01 on the night of our observations. The observations were conducted in the pupil tracking mode to enable the angular differential imaging (ADI; Marois et al. 2006) technique. We took 35 object frames with an exposure time of 30 s. The total exposure time was 17.5 minutes. As a first step in the reduction, we removed a characteristic stripe bias pattern from each image (Suzuki et al. 2010). Then, bad pixel and flat field corrections were performed. Finally, the image distortions were corrected, using calibration images of the globular cluster M5 that were obtained on the same night.

We observed one transit of K2-19 on 2015 January 28 (UT) with the 1.2 m telescope at the Fred Lawrence Whipple Observatory (FLWO) on Mt. Hopkins, Arizona. We used Keplercam, which is equipped with a 4096 × 4096 pixel CCD with a 231 × 231 FOV. We observed through a Sloan i' filter. The exposure time was 30 s. Debiasing and flat-fielding (using dome flats) were performed using standard IRAF procedures. Aperture photometry was performed with custom routines written in the interactive data language (IDL). We selected the final aperture size of 7 pixels in the 2 × 2 binning configuration, which means the radius of the aperture is about 5''. The sky level per pixel was estimated from the median value in an annulus surrounding the aperture, with a radius that is about twice the aperture radius. The time of each exposure was recorded in UT, and the systematic error of the recorded time with respect to the standard clock was much smaller than the statistical uncertainty for the mid-transit time.

One transit of K2-19b was observed on the night of 2015 February 28 (UT) with the 0.6 m TRAPPIST robotic telescope (TRAnsiting Planets and PlanetesImals Small Telescope), located at ESO La Silla Observatory (Chile). TRAPPIST is equipped with a thermoelectrically cooled 2K × 2K CCD, which has a pixel scale of 065 that translates into a 22' × 22' FOV. For details of TRAPPIST, see Gillon et al. (2011) and Jehin et al. (2011). The transit was observed in an Astrodon Exoplanet (blue-blocking) filter that has a transmittance over 90% from 500 nm to beyond 1000 nm. The exposure time was 10 s. The time of each exposure was recorded in HJD_UTC. During the run, the positions of the stars on the chip were maintained to within a few pixels thanks to a software guiding system that regularly derives an astrometric solution for the most recently acquired image and sends pointing corrections to the mount if needed. After a standard pre-reduction (bias, dark, and flat-field correction), the stellar fluxes were extracted from the images using the IRAF/DAOPHOT aperture photometry software (Stetson 1987). After testing several sets of reduction parameters, we chose the one giving the most precise photometry for stars of similar brightness to the target. Differential photometry of the target star was performed relative to a selected set of reference stars. The set of reference stars was chosen as the one that gives the lowest root-mean-square (rms) transit light curve. It consists of nine stable stars of similar brightness and color to the target.

We observed one transit of K2-19b with a new multi-color camera, named MuSCAT, installed on the 1.88 m telescope in Okayama Astrophysical Observatory (OAO) on 2015 April 25 (UT). On that night, a full transit was predicted to be observable based on the ephemeris of Foreman-Mackey et al. (2015). However, based on the updated ephemeris of Armstrong et al. (2015) (taking TTVs into account) only a partial transit was predicted to be observable. MuSCAT is equipped with three CCD cameras and is capable of three-color simultaneous imaging in ${g}_{2}^{\prime },$ ${r}_{2}^{\prime },$ and z_s,2 bands of Astrodon Sloan Gen 2 filters (Fukui et al. 2015; Narita et al. 2015). Each CCD camera has a 61 × 61 FOV, with a pixel scale of about 0358. For the ${g}_{2}^{\prime }$ and ${r}_{2}^{\prime }$ band observations, we employed the high-speed readout mode (2 MHz, corresponding to readout time of about 0.58 s and readout noise of ∼10 e⁻). For the z_s,2 band observations we used the low-speed readout mode (100 kHz, corresponding to readout time of about 10 s and readout noise of ∼4 e⁻) because at that time the z_s,2-band CCD was affected by excess readout noise (over 20 e⁻) in the high-speed readout mode.¹⁵ The exposure times were 60 s, 30 s, and 60 s for ${g}_{2}^{\prime },$ ${r}_{2}^{\prime },$ and z_s,2 bands, respectively. Bias subtraction, flat-fielding, and aperture photometry are performed by a customized pipeline (Fukui et al. 2011). Aperture radii are selected as 17, 18, and 19 pixels for ${g}_{2}^{\prime },$ ${r}_{2}^{\prime },$ and z_s,2 bands, respectively. Sky level is measured in the annulus with an inner radius of 45 pixels and an outer radius of 55 pixels. For differential photometry, we select a reference star UCAC4 453-052399 (g = 12.53, r = 12.11, i = 11.96, J = 11.08), which is slightly brighter than the target star (g = 13.36, r = 12.76, i = 12.57, J = 11.60). The aperture radii and the reference star were chosen so that rms of residuals for trial transit fittings are minimized.

For a more comprehensive transit analysis, we also analyze the previously published data from the NITES (Near Infrared Transiting ExoplanetS: McCormac et al. 2014) 0.4 m telescope in the same manner as our own data. The NITES data are identical to those presented by Armstrong et al. (2015). The transit observed by NITES was the same as the one observed by TRAPPIST.

Using the high-resolution spectrum obtained with Subaru/HDS, we extract the atmospheric parameters (the effective temperature T_eff, surface gravity log g, metallicity [Fe/H], and microturbulent velocity ξ) by measuring the equivalent widths of iron i/ii lines (Takeda et al. 2002), which are widely distributed throughout the observed wavelength region. We then translate those atmospheric parameters into estimates of the stellar mass and radius, using the Yonsei–Yale stellar-evolutionary models (Yi et al. 2001). Finally, the projected stellar rotation velocity V sin I_s is estimated by fitting the observed spectrum with the theoretically generated intrinsic stellar absorption line profile, convolved with a kernel taking into account rotational and macroturbulent broadening (Gray 2005) and the instrumental profile (IP) of Subaru/HDS. The theoretical synthetic spectrum was taken from the ATLAS9 model (a plane-parallel stellar atmosphere model in LTE: Kurucz 1993, No. 13). The result of the spectroscopic analysis is summarized in Table 1. See Hirano et al. (2014) for details on how the uncertainties in these spectroscopic parameters are determined. The derived parameters are consistent with those reported by Armstrong et al. (2015), and have a higher precision by about an order of magnitude thanks to higher S/N spectra. We also find that the host star is likely older than ∼8 Gyr, which is far older than the typical timescale for planet formation and migration (∼10⁷ years, see also Section 4.2).

Table 1.  Stellar Parameters Derived from High-dispersion Spectroscopy

Parameter	Value	Errora
Teff (K)	5345	17
log g	4.394	0.050
[Fe/H]	0.07	0.03
Rs (R⊙)	0.914	0.027
Ms (M⊙)	0.902	0.011
V sin Is (km s−1)	0.85	${}_{-0.85}^{+0.95}$
Age (Gyr)	≥8	...

Note.

^aPresented errors are derived from statistical measurement errors in equivalent widths of iron lines. Additional systematic error (say, 40 K in T_eff) may exist. Discussions on such systematic errors are presented in Hirano et al. (2014).

Download table as:  ASCIITypeset image

We employ the locally optimized combination of images (LOCI: Lafrenière et al. 2007) algorithm to maximize the efficiency of the ADI technique and to search for fainter objects in the inner region around K2-19. We estimate our achieved contrast ratio as follows. The final image of a LOCI pipeline is convolved with the point-spread function (PSF) of the unsaturated image. The 1σ limit on the contrast ratio is defined as the ratio of the standard deviation of the pixel values inside an annulus of width 6 pixels to the stellar flux in the unsaturated image. Then the contrast ratio is corrected for a self-subtraction effect, which is estimated from the recovery rate of injected signals. Figure 1 shows a S/N map of a LOCI-reduced image around K2-19 (north is up and east is to the left). There are no other faint sources around K2-19 with S/N exceeding 5, except for a faint source located at the west-northwest edge of the image (∼115 from the center). The faint source has a contrast of ∼4 × 10⁻⁴ relative to the brightness of K2-19 and is located outside Kepler's PSF. Thus the presence of the faint source does not affect the light curves from K2 and other ground-based telescopes. We plot the 5σ contrast limit in Figure 2. Contrast limits of ∼2 × 10⁻³ and ∼10⁻⁴ are achieved at distances of at 04 and 1'' from the host star, respectively. Thus we have not identified any contaminating faint object that can dilute or mimic the transits of K2-19b and K2-19c.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We analyze the light curve for K2-19 that has been produced from the raw pixel data by the ESPRINT collaboration. See Sanchis-Ojeda et al. (2015) for a description of the procedures for extracting the light curve from the K2 pixel data. We fit the K2 light curve by the procedure below and derive the best-fit system parameters along with mid-transit times. First, we separate the entire light curve into 14 segments, each of which spans a transit. Three of the segments involve double transits (see Figure 3), during which two planets simultaneously transit the host star. We note that Armstrong et al. (2015) did not model those mutual transit events. Each segment includes the data within ∼4.5 hr of the predicted ingress/egress times (based on the ephemeris by Armstrong et al. 2015), which is sufficient to span the transit event as well as some time beforehand and afterward. For each segment, we compute the standard deviation of the out-of-transit fluxes, and adopt the standard deviation as an initial estimate of the uncertainty in each flux value.

Zoom In Zoom Out Reset image size

All the light curve segments are fitted simultaneously with the common system parameters along with mid-transit time(s) for each segment. We compute the posterior distributions for those parameters based on the Markov Chain Monte Carlo (MCMC) algorithm, assuming that the likelihood is proportional to exp(−χ²/2) where

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}\;\displaystyle \frac{{({f}_{\mathrm{LC},\mathrm{obs}}^{(i)}-{f}_{\mathrm{LC},\mathrm{model}}^{(i)})}^{2}}{{\sigma }_{\mathrm{LC}}^{(i)2}},\end{eqnarray} \tag{ 1 }$$

where ${f}_{\mathrm{LC},\mathrm{obs}}^{(i)}$ and ${\sigma }_{\mathrm{LC}}^{(i)}$ are the ith observed K2 flux and its error, respectively. We employ the transit model by Ohta et al. (2009) integrated over the 29.4 minute averaging time of the Kepler data. We neglect the impact of any possible planet–planet eclipses (Hirano et al. 2012) on the light curve shape, due to the rather sparse sampling of the K2 data. Finally, to obtain the model flux ${f}_{\mathrm{LC},\mathrm{model}}^{(i)},$ the integrated light curve model for each segment is multiplied by a second-order polynomial function of time, representing the longer-timescale flux variations at the time of the transit. Therefore, the adjustable parameters in our model are (for each planet) the scaled semimajor axis a/R_s, the transit impact parameter b, and the planet-to-star radius ratio R_p/R_s, the quadratic limb-darkening parameters u₁ and u₂, the mid-transit time(s) T_c for each segment, and the coefficients of the second-order polynomials. The sparse sampling of the K2 data relative to the timescale of ingress and egress prohibits us from constraining the orbital eccentricity e of either planet based only on the transit light curve (Dawson & Johnson 2012), so we simply assume e = 0 for both planets.

Following the procedure described in Hirano et al. (2015), we compute the posterior distributions for the parameters listed above; we first optimize Equation (1) using Powell's conjugate direction method, allowing all the relevant parameters to vary. At this point, we compute the level of time-correlated noise β (so-called red noise: Pont et al. 2006; Winn et al. 2008) for each segment, and inflate the errors in each segment by β. We then fix the coefficients of the polynomial functions to the optimized values and perform the MCMC computation starting from the initial best-fit values for the other parameters. The step size for each parameter is iteratively optimized such that the final acceptance ratio over the whole chain is between 15% and 35%. We extend the chains to 10⁷ links, given the relatively large number of free parameters. We employ the median, and 15.87 and 84.13 percentiles of the marginalized posterior distribution of each parameter to convey the representative value and its ±1σ errors. As a check on convergence, we repeat the analysis after changing the initial input values for the system parameters, and find no dependence on the initial values. We report the values and uncertainties of the basic transit parameters in Table 2, and we report the mid-transit times as well as the β factor for each segment in Table 3. Combining the R_p/R_s values with the R_s value listed in Table 1, we estimate the planetary radii of K2-19b and K2-19c as R_p,b = 7.34 ± 0.27 R_⊕ and R_p,c = 4.37 ± 0.22 R_⊕, respectively. The planetary radii are also presented in Table 2.

Table 2.  Best Fitting Parameters and Errors for K2 Transit Light Curves

Parameter	Value	Error
a/Rs,b	19.35	${}_{-1.45}^{+0.56}$
a/Rs,c	24.09	${}_{-3.01}^{+1.72}$
bb	0.233	${}_{-0.163}^{+0.213}$
bc	0.367	${}_{-0.246}^{+0.260}$
Rp/Rs,b	0.0737	${}_{-0.0011}^{+0.0016}$
Rp/Rs,c	0.0439	${}_{-0.0009}^{+0.0018}$
T14,b (day)	0.1362	${}_{-0.0014}^{+0.0020}$
T14,c (day)	0.1536	${}_{-0.0027}^{+0.0042}$
${u}_{1}+{u}_{2}$	0.69	${}_{-0.14}^{+0.18}$
${u}_{1}-{u}_{2}$	0.06	${}_{-0.38}^{+0.35}$
Rp,b (R⊕)	7.34	±0.27
Rp,c (R⊕)	4.37	±0.22

Note. Derived planetary radii are also presented.

Download table as:  ASCIITypeset image

Our procedure for modeling the ground-based transit light curves follows Narita et al. (2013). First, the time stamps of the photometric data are placed onto the BJD_TDB system using the code by Eastman et al. (2010). We check all the light curves by eye and eliminate obvious outliers. We then fit the transit light curves simultaneously, using an analytic transit light curve model and various choices to describe the more gradual out-of-transit (baseline) flux variations.

We adopt baseline model functions F_oot described as follows:

$$\begin{eqnarray*}{F}_{{\rm{oot}}} & = & {k}_{0}\times {10}^{-0.4{\rm{\Delta }}{m}_{{\rm{cor}}}},\\ {\rm{\Delta }}{m}_{\mathrm{cor}} & = & \displaystyle \sum \;{k}_{i}{X}_{i},\end{eqnarray*}$$

where k₀ is the normalization factor, $\{{\boldsymbol{X}}\}$ are observed variables, and $\{{\boldsymbol{k}}\}$ are the coefficients (Fukui et al. 2013). In order to select the most appropriate baseline models for the light curves, we adopt the Bayesian Information Criteria (BIC: Schwarz 1978). The BIC value is given by $\mathrm{BIC}\equiv {\chi }^{2}+k\mathrm{ln}N,$ where k is the number of free parameters and N is the number of data points. For the variables $\{{\boldsymbol{X}}\}$, we test various combinations of time (t), airmass (z), the relative centroid positions in x (dx) and y (dy), and the sky background counts (s). In addition, we also allowed an additional factor k_flip to account for a change of the normalization factor after the meridian flip for the TRAPPIST light curve. Based on minimum BIC, we adopt k₀, t, and z to be fitting parameters for FLWO and MuSCAT; k₀, z, and k_flip for TRAPPIST; and k₀ and t for NITES.

For the transit light curve model, we employ a customized code (Narita et al. 2007) that uses the analytic formula by Ohta et al. (2009). The formula is equivalent to that of Mandel & Agol (2002) when using the quadratic limb-darkening law, $I(\mu )=1-{u}_{1}(1-\mu )-{u}_{2}{(1-\mu )}^{2},$ where I is the intensity and μ is the cosine of the angle between the line of sight and the line from the position of the stellar surface to the stellar center. We refer to the tables of quadratic limb-darkening parameters by Claret et al. (2013), and compute allowed ${u}_{1}+{u}_{2}$ and ${u}_{1}-{u}_{2}$ values for T_eff = 5300 or 5400 K and log g = 4.0 or 4.5, based on the stellar parameters presented in Table 1. We adopt uniform priors for ${u}_{1}+{u}_{2}$ and ${u}_{1}-{u}_{2}$ as follows: [0.59, 0.61] and [0.24, 0.34] for FLWO i', [0.77, 0.85] and [0.67, 0.82] for MuSCAT ${g}_{2}^{\prime },$ [0.68, 0.70] and [0.35, 0.46] for MuSCAT ${r}_{2}^{\prime },$ [0.52, 0.55] and [0.17, 0.26] for MuSCAT z_s,2, and [0.67, 0.70] and [0.38, 0.50] for TRAPPIST and NITES. We assume an orbital period P = 7.921 days and a reference epoch for the transits T_c,0 = 2457082.6895, the values reported by Armstrong et al. (2015) based on the NITES data. We note that these assumptions have little impact on the resultant transit parameters, because we allow T_c to be a free parameter, and the uncertainty in P is negligible. To constrain the mid-transit times precisely even without complete coverage of the entire transit event, we place an a priori constraint on the total transit duration, T₁₄ = 0.1365 ± 0.0017 day, based on our K2 analysis. We allow R_p/R_s, a/R_s, and the orbital inclination i to be free parameters.

First we optimize free parameters for each light curve, using the AMOEBA algorithm (Press et al. 1992). The penalty function is also given by Equation (1), where ${f}_{{\rm{LC,model}}}^{(i)}$ is a combination of the baseline model and the analytic transit formula mentioned above. Then if the reduced χ² is larger than unity, we rescale the photometric errors of the data such that the reduced χ² for each light curve becomes unity. We also estimate the level of time-correlated noise for each light curve, and further rescale the errors by multiplying by the β factors listed in Table 4. Finally, we use the MCMC method (Narita et al. 2013) to evaluate values and uncertainties of the free parameters. For the FLWO, TRAPPIST, and NITES light curves, we fit each light curve independently. We create five chains of 10⁶ points for each light curve, and discard the first 10⁵ points from each chain (the "burn-in" phase). The jump sizes are adjusted such that the acceptance ratios are 20%–30%. For the OAO/MuSCAT light curves, we fit the three-band data simultaneously, requiring consistency in the parameters T_c, i, and a/R_s. We create five chains of 3 × 10⁶ points for MuSCAT light curves, and discard the first 3 × 10⁵ points from each chain. The acceptance ratios are set to about 25%. Table 4 lists the median values and ±1σ uncertainties, which are defined by the 15.87 and 84.13 percentile levels of the merged posterior distributions. The baseline-corrected transit light curves are plotted in Figure 4 (FLWO), Figure 5 (TRAPPIST and NITES), and Figure 6 (OAO).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Best Fitting Parameters and Uncertainties for Ground-based Transit Light Curves

	FLWO		TRAPPIST		NITES
Parameter	Value	Error	Value	Error	Value	Error
epoch	30		34		34
Tc	2457051.00413	${}_{-0.00225}^{+0.00218}$	2457082.69550	${}_{-0.00116}^{+0.00140}$	2457082.69398	${}_{-0.00435}^{+0.00342}$
Rp/Rs	0.0633	${}_{-0.0045}^{+0.0042}$	0.0682	${}_{-0.0046}^{+0.0043}$	0.0645	${}_{-0.0109}^{+0.0082}$
i (deg)	90.00	${}_{-0.83}^{+0.80}$	90.00	±0.46	90.00	±0.72
a/Rs	19.22	${}_{-0.56}^{+0.44}$	19.58	${}_{-0.34}^{+0.30}$	19.35	${}_{-1.12}^{+0.43}$
k0	0.9899	±0.0025	1.0012	±0.0017	0.9982	${}_{-0.0013}^{+0.0011}$
kt	−0.0077	±0.0024	...		−0.0043	${}_{-0.0035}^{+0.0031}$
kz	0.00199	${}_{-0.00079}^{+0.00075}$	−0.00020	${}_{-0.00035}^{+0.00036}$	...
kflip	...		−0.00030	${}_{-0.00031}^{+0.00033}$	...
β	1.258		1.000		1.000
	MuSCAT (g${}_{2}^{\prime }$)		MuSCAT (r${}_{2}^{\prime }$)		MuSCAT (zs,2)
Parameter	Value	Error	Value	Error	Value	Error
epoch	...		41		...
Tc	...		2457138.15047	${}_{-0.00176}^{+0.00145}$	...
Rp/Rs	0.0753	${}_{-0.0070}^{+0.0067}$	0.0789	${}_{-0.0041}^{+0.0043}$	0.0730	${}_{-0.0061}^{+0.0058}$
i (deg)	...		88.94	${}_{-0.95}^{+0.73}$	...
a/Rs	...		18.80	${}_{-2.10}^{+0.94}$	...
k0	0.9943	${}_{-0.0047}^{+0.0046}$	0.9945	±0.0028	0.9913	±0.0040
kt	−0.0243	${}_{-0.0047}^{+0.0046}$	−0.0107	±0.0027	−0.0095	±0.0041
kz	0.0014	±0.0015	0.00188	±0.00089	0.0030	±0.0013
β	1.059		1.000		1.069

Download table as:  ASCIITypeset image

Since K2-19b and K2-19c are close to the 3:2 MMR, it is not surprising that large TTVs have been observed in this system. To explain the observed mid-transit times, we adopt the analytic formulae for TTV by Deck & Agol (2015) and Nesvorný & Vokrouhlický (2014), who considered the synodic chopping effect. Successful modeling of the synodic chopping effect enables us to estimate the mass of a perturbing body. In particular the resulting mass determination is not subject to the mass-eccentricity degeneracy that was discussed by Lithwick et al. (2012), who derived analytic formulae for a pair of planets near a first-order MMR. Although the chopping analytic formulae are only applicable to systems with nearly circular and coplanar orbits, those requirements are likely fulfilled in this case, given the presence of mutual transits and the requirement for long-term dynamical stability.

The chopping formulae predict the TTVs of the inner and outer planets as follows:

$$\begin{eqnarray}&&\delta {t}_{{\rm{b}}}=\displaystyle \sum _{j=1}^{\infty }\;\displaystyle \frac{{P}_{{\rm{b}}}}{2\pi }{\mu }_{{\rm{c}}}{f}_{{\rm{b}}}^{(j)}(\alpha )\mathrm{sin}{\psi }_{j},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\delta {t}_{{\rm{c}}}=\displaystyle \sum _{j=1}^{\infty }\;\displaystyle \frac{{P}_{{\rm{c}}}}{2\pi }{\mu }_{{\rm{b}}}{f}_{{\rm{c}}}^{(j)}(\alpha )\mathrm{sin}{\psi }_{j},\end{eqnarray} \tag{ 3 }$$

where μ_b,c are the planet-to-star mass ratios (M_b,c/M_s) and

$$\begin{eqnarray}&&\left({f}_{{\rm{b}}}^{(j)},(\alpha ),=,-\alpha ,(,j({\beta }^{2}+3){b}_{1/2}^{j}(\alpha )+2\alpha \beta \displaystyle \frac{\partial }{\partial \alpha }{b}_{1/2}^{j}(\alpha )\right.\\ &&-\left.\alpha {\delta }_{j,1}({\beta }^{2}+2\beta +3),\phantom{\displaystyle \frac{\partial }{\partial \alpha }},\phantom{\rule{-1.6em}{0ex}}\right)/\left({\beta }^{2}({\beta }^{2}-1)\right),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\left({f}_{{\rm{c}}}^{(j)},(\alpha ),=,(,j({\kappa }^{2}+3){b}_{1/2}^{j}(\alpha )+2\kappa \left(\alpha \displaystyle \frac{\partial }{\partial \alpha }{b}_{1/2}^{j}(\alpha )\right.\right.\\ &&+\left.\left.{b}_{1/2}^{j}(\alpha ),\phantom{\displaystyle \frac{\partial }{\partial \alpha }},\phantom{\rule{-1.6em}{0ex}}\right)-{\alpha }^{-2}{\delta }_{j,1}({\kappa }^{2}-2\kappa +3),\phantom{\displaystyle \frac{\partial }{\partial \alpha }},\phantom{\rule{-1.6em}{0ex}}\right)/\\ &&\times \left({\kappa }^{2}({\kappa }^{2}-1)\right),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{b}_{1/2}^{j}(\alpha )=\displaystyle \frac{1}{\pi }{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{\mathrm{cos}(j\theta )}{\sqrt{1-2\alpha \mathrm{cos}\theta +{\alpha }^{2}}}d\theta .\end{eqnarray} \tag{ 6 }$$

Note that P_b,c and μ_b,c are defined as non-negative variables. The parameters α, β, and κ are defined as α = a_b/a_c, $\beta =j({n}_{{\rm{b}}}-{n}_{{\rm{c}}})/{n}_{{\rm{b}}},$ and $\kappa =j({n}_{{\rm{b}}}-{n}_{{\rm{c}}})/{n}_{{\rm{c}}},$ where a_b,c and n_b,c are the semimajor axes and the mean motions of the planets. Those parameters can be expressed by the periods of the planets: α ≃ (P_b/P_c)^2/3 (using Kepler's third law and neglecting the planet masses in comparison to the stellar mass), $\beta =j(1-{P}_{{\rm{b}}}/{P}_{{\rm{c}}}),$ and $\kappa =j({P}_{{\rm{c}}}/{P}_{{\rm{b}}}-1).$ Note that δ_j,1 represents Kronecker's delta, which is 1 when j = 1 and 0 for j > 1. When the time of conjunction (t_conj) is used as the origin of the time scale, the phase ψ_j can be expressed simply as

$$\begin{eqnarray}&&{\psi }_{j}=2\pi {tj}(1/{P}_{{\rm{b}}}-1/{P}_{{\rm{c}}})\equiv 2\pi {tj}/{P}_{{\rm{syn}}},\end{eqnarray} \tag{ 7 }$$

where P_syn is the synodic period. For the current system, we derive t_conj = 2456852.9344 ± 0.0022 BJD_TDB from one of the mutual transit events, assuming that both planets are orbiting in the same direction. Consequently, the observed mid-transit times (T_c) for the planets can be modeled by six free parameters T_c(0)_b,c, P_b,c, and μ_b,c as

$$\begin{eqnarray}&&{T}_{{\rm{c}}}{({E}_{{\rm{b}}})}_{{\rm{b}}}={T}_{{\rm{c}}}{(0)}_{{\rm{b}}}+{P}_{{\rm{b}}}\times {E}_{{\rm{b}}}+\delta {t}_{{\rm{b}}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{T}_{{\rm{c}}}{({E}_{{\rm{c}}})}_{{\rm{c}}}={T}_{{\rm{c}}}{(0)}_{{\rm{c}}}+{P}_{{\rm{c}}}\times {E}_{{\rm{c}}}+\delta {t}_{{\rm{c}}},\end{eqnarray} \tag{ 9 }$$

where E_b,c are the transit epochs of the planets with origins at the first K2 transits.

We first try to find an optimal parameter set for the six parameters (T_c(0)_b, T_c(0)_c, P_b, P_c, μ_b and μ_c) by minimizing

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}\;\displaystyle \frac{{({T}_{{\rm{c,obs}}}^{(i)}-{T}_{{\rm{c,model}}}^{(i)})}^{2}}{{\sigma }_{{T}_{{\rm{c}}}}^{(i)2}},\end{eqnarray} \tag{ 10 }$$

using the AMOEBA algorithm. In evaluating Equations (2) and (3), we truncate the series at j = 7 to reduce the computational cost. We model the observed the mid-transit times for both planets listed in Tables 3 and 4 (21 in total) using the equations above. However, we find that the AMOEBA algorithm does not converge. Note that we have tested the effects of including higher orders of j in the calculation, but the result remains unchanged. Apparently there are multiple local minima of χ² in the parameter space, which are almost equally favored. This result is quite understandable: we have not observed any additional transits of planet c since the end of the K2 campaign, and hence the TTVs of planet c are poorly constrained.

For this reason, we decide to fit only the mid-transit times of the inner planet (14 in total) using the corresponding four parameters (${T}_{{\rm{c}}}{(0)}_{{\rm{b}}},$P_b, P_c, and μ_c). This time the AMOEBA algorithm converges. Table 5 summarizes the optimal parameters and Figure 7 plots the optimal model explaining the observed TTVs of the inner planet. We note that we find a somewhat large reduced χ², namely χ²/ν ∼ 4.4, where ν = 10 is the number of degrees of freedom. This may indicate that the uncertainties in T_c, particularly those based on K2 data, have been underestimated by a factor of about 2.1. To account for possible systematic errors, the uncertainties given in Table 5 have been enlarged by a factor of $\sqrt{{\chi }^{2}/\nu }$ from the original 1σ errors estimated by the usual criterion Δχ² = 1.0.

Zoom In Zoom Out Reset image size

Table 5.  Planetary Parameters Derived from the TTV Analysis

	Minimum χ2		Second Minimum χ2
Parameter	Value	Errora	Value	Errora
Pb (day)	7.920994	0.000071	7.921122	0.000176
Pc (day)	12.0028	0.0092	11.7748	0.0142
μc	0.0000713	0.0000057	0.0000672	0.0000087
Tc(0)b	2456813.37624	0.00050	2456813.37466	0.00109
χ2/ν	43.98/10	...	49.38/10	...
Mc	21.4	1.9	20.2	2.7
Pc/Pb	1.51531	0.00117	1.48651	0.00180
${f}_{{\rm{b}}}^{(1)}(\alpha )$	11.91	...	13.51	...
${f}_{{\rm{b}}}^{(2)}(\alpha )$	20.09	...	21.19	...
${f}_{{\rm{b}}}^{(3)}(\alpha )$	−139.71	...	178.84	...
${f}_{{\rm{b}}}^{(4)}(\alpha )$	−4.10	...	−5.56	...
${f}_{{\rm{b}}}^{(5)}(\alpha )$	−1.22	...	−1.57	...
${f}_{{\rm{b}}}^{(6)}(\alpha )$	−0.50	...	−0.64	...
Psyn	23.29234	...	24.20293	...
α	0.75799	...	0.76776	...

Notes. Coefficients for the synodic chopping formulae are also presented.

^aPresented errors are inflated by $\sqrt{{\chi }^{2}/\nu }$ to account for possible systematic errors.

Download table as:  ASCIITypeset image

We then check on the existence of possible local minima of χ² by fixing P_c close to the optimal value and optimizing the other parameters using the AMOEBA algorithm. We repeat this exercise for values of P_c ranging from 11.5 to 12.5. Figure 8 shows the variation in χ² and the optimal parameters with varying P_c. We find that the second-minimum χ² is located at P_c ∼ 11.775 days with Δχ² = 5.4 from the minimum χ². Figure 9 presents a TTV model for the second-minimum χ². The behavior of the second-minimum χ² model is almost identical to that of the minimum χ² model in the observing period. Thus we cannot exclude this solution statistically at this point. Therefore, Table 5 also gives the parameters and errors based on the second-minimum χ² solution. We also note that the third-minimum χ² is located at P_c ∼ 12.17 days with χ² ∼ 67.50, which is statistically less favorable than the above two solutions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As seen in Figure 8, the orbital period for the inner planet P_b is robustly determined to be about 7.921 days. The minimum and second-minimum χ² are thus located slightly above (P_c/P_b = 1.51531 ± 0.00117) and below (P_c/P_b = 1.48651 ± 0.00180) the exact 3:2 commensurability (P_c ∼ 11.882 days). The orbital period of the outer planet is not permitted to lie exactly on the commensurability, where larger TTVs are predicted and the synodic chopping formulae are less reliable. This is why the optimization results in the mass of the planet c being driven to zero when P_c is held fixed near the 3:2 commensurability. The planet-to-star mass ratio for the outer planet μ_c is estimated as 0.0000713 ± 0.0000027 for the minimum χ² solution and 0.0000672 ± 0.0000039 for the second-minimum. According to the mass of the host star presented in Table 1, these values correspond to a mass for planet c of 21.4 ± 1.9 M_⊕ and 20.2 ± 2.7 M_⊕ for the minimum and second-minimum solutions, respectively. Thus the mass of the outer planet is well constrained to be ∼20 M_⊕. We note that the error in t_conj has less impact on the above results. Finally, we present coefficients of the synodic chopping formulae in Table 5 so as to facilitate transit predictions for the inner planet in future years.

Very recently, Barros et al. (2015) reported a photodynamical modeling of TTVs and radial velocities of the K2-19 system. They derived the mass of the planet c as ${15.9}_{-2.8}^{+7.9}{M}_{\oplus },$ which is consistent with our result. We note that the derived orbital periods of both planets are also consistent when taking account of a difference in the definitions¹⁶ of P_b and P_c.

K2-19c has a relatively low mean density of 1.43 g cm⁻³. It is therefore likely to have a gaseous atmosphere on top of any solid component (see also Figure 10). We explore models for the interior structure of K2-19c (21.4 M_⊕, 4.37 R_⊕ at 0.1024 AU) based on the condition of hydrostatic equilibrium, using equations of state (EOSs) for four constituents: SCvH EOS for H/He (Saumon et al. 1995); the ab initio EOS (French et al. 2009) and the SESAME EOS 7150 (Lyon & Johnson 1992) for water; and the SESAME EOSs 7100 and 2140 for rocky material and iron, respectively. We calculate a pressure–temperature profile of an irradiated planet's atmosphere in radiative equilibrium, using the analytical formulae of Parmentier & Guillot (2014) and gas opacities derived by Freedman et al. (2008). We assume that both the incoming and outgoing radiation fields are isotropic, the bond albedo of the planet is zero, and the incoming flux is averaged over the dayside hemisphere. The atmospheric mass fraction of K2-19c is estimated to be ∼7% for an Earth-like core model (32.5% iron core, 67.5% silicate mantle) and ∼1.3% for an icy core model (ice:rock = 2.7:1), where we assumed that K2-19c has a H/He atmosphere with the stellar metallicity of [Fe/H] = 0.07. We find that K2-19c likely has a H/He atmosphere of ≲10 wt% despite the fact that its total mass exceeds the critical core mass for gas accretion within the protoplanetary disk, which is thought to be ∼10 M_⊕ (e.g., Pollack et al. 1996). This suggests that the atmosphere of K2-19c might have been removed via giant impacts and/or a mass loss driven by stellar irradiation.

Zoom In Zoom Out Reset image size

The mass of K2-19b remains undetermined. Here, we discuss the possible mass range of K2-19b (7.34 R_⊕ at 0.0823 AU) from the standpoint of the formation and evolution of close-in planets. Under the extreme assumption that K2-19b has no gaseous atmosphere and is purely solid, its observed radius would require its mass to be >100 M_⊕. The formation of such a massive solid core without any gaseous envelope seems unlikely. Such a massive core would be theoretically expected to initiate the (in situ) accretion of the ambient disk gas (e.g., Ikoma & Hori 2012). A core of this mass has been inferred for HD 149026b (Sato et al. 2005), but in that case the core is accompanied by a substantial gaseous atmosphere. Consequently, for K2-19b, it seems most plausible to consider models that possess an atmosphere. However, we note that K2-19b, orbiting at only ∼0.08 AU, may have undergone atmospheric loss due to stellar irradiation. As K2-19 is relatively old (≥8 Gyr), a Neptune-mass planet could have lost a substantial amount of its atmosphere, whereas a Jupiter-mass planet would likely have retained its original atmosphere (e.g., Owen & Wu 2013).

With these considerations in mind, we consider two models: K2-19b is either a Neptune-like planet (≲20 wt% of a H/He atmosphere) or a gas giant. We first suppose that K2-19b has a H/He atmosphere of 1 wt%, 10 wt%, or 20 wt% on top of a solid core. The atmosphere is taken to have the same metallicity as the host star, and the core is taken to have an Earth-like composition. Under these assumptions, the mass of K2-19b is expected to be ∼5.1 M_⊕, 18 M_⊕, or 37 M_⊕, respectively. Second, we suppose that K2-19b has a mass of 1 M_Jup. In this case the models indicate that the atmospheric mass fraction is ∼35%. Accordingly, K2-19b can be a sub-Neptune-mass planet with a tenuous atmosphere, a super-Neptune planet, or a close-in gas giant with an extremely massive solid core. Recently, Barros et al. (2015) derived the mass of K2-19b as 44 ± 12 M_⊕ from photodynamical modeling. Their result supports the idea that K2-19b may be surrounded by a thick H/He atmosphere of ∼20 wt%.

The K2-19 system is the first discovery of close-in Neptune-sized planets near a 3:2 MMR. A slow convergent migration allows low-mass planets to be locked into a first-order commensurability (e.g., Papaloizou & Szuszkiewicz 2005). In fact, Figure 11 shows that all of the other pairs of planets close to the 3:2 commensurability found by Kepler are smaller than 3 R_⊕. Also, planets with P ≲ 10 days near a 2:1 MMR have radii smaller than about 3 R_⊕. The peculiar pair close to the 3:2 MMR around K2-19 is suggestive of inward transport of Neptune-sized planets through a convergent migration. We find there is a lack of close-in gas giants among both 3:2 MMR and 2:1 MMR systems.¹⁷ Besides, 3:2 MMR systems seem to be less common than 2:1 MMR ones. Dynamical instability of more closely packed 3:2 MMR systems may be responsible for the feature. In addition, there is no gas giant in a 3:2 MMR so far discovered within P = 50 days. For 2:1 MMR systems, gas giants such as Kepler-9b, 9c (Holman et al. 2010), and Kepler-30c (Fabrycky et al. 2012) form the 2:1 commensurability beyond P ∼ 10 days. The distinct habitat may imply that large planets with radius of ≳3 R_⊕ most likely fall into a 2:1 MMR rather than a 3:2 MMR due to the strength of planetary migration.

Zoom In Zoom Out Reset image size

An outer planet formed in a 2:1 MMR or 3:2 MMR tends to be larger than an inner one (see Figure 11). This may imply that an outer planet has caught up with an inner one to form a first-order MMR system in a protoplanetary disk. The inner planet closer to the central star loses its atmosphere more readily via stellar X-ray and UV (XUV) irradiation. The strength of the stellar XUV flux that a pair of planets in a p:p + 1 MMR receives is related as $\frac{{F}^{{\rm{out}}}}{{F}^{{\rm{in}}}}\propto {\left(\frac{p}{p+1}\right)}^{4/3},$ where p is the commensurability integer and F is the stellar XUV flux. The subtle but non-negligible flux difference may lead to the prevalence of larger outer planets in a 2:1 MMR or a 3:2 MMR. However, the radius of K2-19b is about 1.68 times as large as that of the Neptune-sized K2-19c. This means that the K2-19 system has a unique 3:2 MMR configuration, unlike other close-in MMR systems with smaller planets. The determination of the mass of K2-19b should play a crucial role in understanding what causes their contrast in radius, and why the close-in super-Neptune and Neptune pair was trapped into a 3:2 MMR rather than a 2:1 MMR, for example a difference in migration speed (Ogihara & Kobayashi 2013).