Validation of a Temperate Fourth Planet in the K2-133 Multi-planet System

TL;DR

This paper validates the existence of a fourth temperate planet in the K2-133 system, providing insights into planetary formation and atmospheric evolution, and highlighting the potential for rocky composition in the outer planet.

Contribution

The study confirms a fourth planet in the K2-133 system and discusses its implications for planetary population gaps and atmospheric photo-evaporation theories.

Findings

The fourth planet has a radius of approximately 1.73 Earth radii.

The four planets span the radius gap, supporting atmospheric loss theories.

K2-133 e may be a rocky world on the edge of the temperate zone.

Abstract

We present follow-up observations of the K2-133 multi-planet system. Previously, we announced that K2-133 contained three super-Earths orbiting an M1.5V host star - with tentative evidence of a fourth outer-planet orbiting at the edge of the temperate zone. Here we report on the validation of the presence of the fourth planet, determining a radius of $1.73_{-0.13}^{+0.14}$ R$_{\oplus}$. The four planets span the radius gap of the exoplanet population, meaning further follow-up would be worthwhile to obtain masses and test theories of the origin of the gap. In particular, the trend of increasing planetary radius with decreasing incident flux in the K2-133 system supports the claim that the gap is caused by photo-evaporation of exoplanet atmospheres. Finally, we note that K2-133 e orbits on the edge of the stars temperate zone, and that our radius measurement allows for the possibility that this is a rocky world. Additional mass measurements are required to confirm or refute this scenario.

Figures

Tables

Table 1. Table 1: Refined stellar properties of K2-133 from Gaia astrometry, and their 1 ​ σ 1 𝜎 1\sigma uncertainties.

Property	Value	Unit
Gaia astrometry:
Ref. Epoch	2015.5	Jyr
R.A.	$70.1494$	deg.
Dec	$25.0098$	deg.
${\mu }_{\mathrm{RA}}$	$185.66\pm 0.08$	mas yr-1
${\mu }_{\mathrm{Dec}}$	$-46.36\pm 0.04$	mas yr-1
Radial velocity	$96.9\pm 1.7$	km s-1
Parallax	$13.27\pm 0.04$	mas
Derived properties:
$T_{\mathrm{eff}}$	$3655\pm 80$	K
Distance	$75.2\pm 0.2$	pc
Stellar radius	$0.455\pm 0.022$	R⊙
Stellar mass	$0.461\pm 0.011$	M⊙
Stellar luminosity	$0.0332\pm 0.0013$	L⊙
Space velocity	$107.2\pm 1.4$	km s-1

Table 2. Table 2: Refined planetary parameters of the system with uncertainties given at the 1 ​ σ 1 𝜎 1\sigma confidence level.

Property	Planet b	Planet c	Planet d
Fitted properties:
$T_0$ (BJD$-2454833$)	${2988.3150}_{-0.0015}^{+0.0016}$	$2990.7686\pm 0.0011$	$2993.1709\pm 0.0012$
Period (days)	$3.07133\pm 0.00011$	$4.86784\pm 0.00012$	${11.02454}_{-0.00035}^{+0.00036}$
$R_p/R_s$	$0.0270\pm 0.0008$	$0.0323\pm 0.0009$	${0.0404}_{-0.0010}^{+0.0008}$
Inclination (deg.)	${87.60}_{-0.15}^{+0.18}$	${88.21}_{-0.08}^{+0.11}$	${89.40}_{-0.08}^{+0.10}$
Derived properties:
Radius ($R_{\oplus }$)	${1.340}_{-0.076}^{+0.077}$	${1.603}_{-0.088}^{+0.090}$	${2.003}_{-0.107}^{+0.107}$
Semi-major axis (au)	${0.03194}_{-0.00026}^{+0.00025}$	${0.04341}_{-0.00035}^{+0.00034}$	${0.07487}_{-0.00060}^{+0.00059}$
Impact parameter	${0.630}_{-0.054}^{+0.052}$	${0.640}_{-0.047}^{+0.046}$	${0.370}_{-0.066}^{+0.055}$
Non-linear LDCs	c${}_2{}^{}={3.4}_{-1.9}^{+2.1}$	c${}_3{}^{}=-{5.6}_{-4.4}^{+3.5}$	c${}_4{}^{}={2.9}_{-1.9}^{+2.4}$

Table 3. Table 3: Parameters of planet e with uncertainties given at the 1 ​ σ 1 𝜎 1\sigma confidence level. ‡ Computed with supersampling factor = 1. † Assuming an albedo of 0.3.

Property	Value
Fitted properties:
T0 (BKJD)	${3004.8657}_{-0.0023}^{+0.0022}$
Period (days)	${26.5841}_{-0.0017}^{+0.0018}$
Rp/R⋆	${0.0348}_{-0.0019}^{+0.0023}$
Inclination (deg.)	${89.164}_{-0.010}^{+0.012}$
Derived properties:
Radius (R⊕)	${1.73}_{-0.13}^{+0.14}$
Semi-major axis (R⋆)	${63.6}_{-3.0}^{+3.3}$
Semi-major axis (AU)	$0.1346\pm 0.0011$
Impact parameter	${0.928}_{-0.045}^{+0.049}$
Transit depth (ppm)‡	${987}_{-81}^{+86}$
Duration (hours)‡	${1.44}_{-0.06}^{+0.09}$
Incident flux (S⊕)	${1.80}_{-0.23}^{+0.25}$
Equilibrium temp. (K)†	$296\pm 10$

Equations

$$ESI=1-\sqrt{\frac{1}{2}\left({\left(\frac{S-1}{S+1}\right)}^2+{\left(\frac{R-1}{R+1}\right)}^2\right)}$$

Full text

††pubyear: 2018††pagerange: Validation of a Temperate Fourth Planet in the K2-133 Multi-planet System–Validation of a Temperate Fourth Planet in the K2-133 Multi-planet System

1 Introduction

M-dwarf stars provide the best opportunity to study potentially habitable planets in the near-future. They are the most common stars in the Galaxy and have been estimated to host a rate of 0.45 small ($1-2R_{\earth}$) temperate planets per star, using the recent Venus and early Mars boundary criteria (Dressing & Charbonneau, 2015). Due to their small stellar radii, transit signals of planets orbiting M-dwarfs are larger than for those around earlier-type stars. In addition, the orbital periods of planets in the temperate zones of M-dwarfs are very short (typically days to a few weeks), making their detection and characterisation more accessible. Similarly, the lower masses of these stars allow more accurate mass measurements of their planets. The greater transit signals also improve the studies of the atmospheres of these worlds via techniques such as transmission spectroscopy, potentially allowing searches for biosignatures that are currently extremely challenging for Earth-like planets in Earth-like orbits around solar-like stars (Irwin et al., 2014).

A multi-planet system was discovered around K2-133 – an early M-star, using data from campaign 13 of the Kepler-K2 mission (Wells et al., 2018). The system consists of three confirmed super-Earths increasing in size with distance from the star – possibly caused by photo-evaporation of their atmospheres. A candidate fourth super-Earth was also stated in the discovery paper, on a longer period orbit. In this work, through a careful re-analysis of the K2 data, we statistically validate this planet candidate with a period of 26.6 days and a radius of $1.73_{-0.13}^{+0.14}$ R*⊕*, which orbits this quiet M-dwarf close to the optimistic habitable zone. The radius of this planet means it is likely to have a gaseous envelope (Rogers, 2015), however we cannot rule out a rocky composition with the current uncertainty. Future photometric observations with higher cadence are needed to further constrain the radius.

Since the discovery, new data has been acquired allowing a more detailed characterisation of the system. In Section 2 we present AO imaging, NIR spectra and Gaia DR2 data which further constrain the properties of the host star and therefore also the planets. In Section 3, we give the results from a statistical test of the probability that the candidate planet is real, and provide its properties. The planet is discussed in detail in Section 4, considering its potential for habitability and scope for future observations. The planetary system as a whole is also discussed in the context of the known exoplanet radius gap.

2 Methods

2.1 Keck AO imaging

High resolution imaging can be used to rule out background eclipsing binary scenarios which can mimic transit signals, and also to identify blended sources which cause us to underestimate transit depths. K2-133 was observed by PI Crossfield on the night of August 20 2017, using the NIRC2 adaptive optics imager on the 10m Keck-II telescope. The data were obtained using the narrow Br-$\gamma$ filter ($\lambda_{\mathrm{eff}}=2.157\mu m$, $\Delta\lambda=0.043\mu m$), combining 9 dithered images with integration times of 10.3 seconds each. More details of the setup used can be found in Ciardi et al. (2018), Crossfield et al. (2016) and Livingston et al. (2018). The data are publicly available at the ExoFOP website111https://exofop.ipac.caltech.edu/k2/edit_target.php?id=247887989.

Fig. 1 shows the final image after flat-fielding, sky subtraction and then co-adding the dithered images. The $5\sigma$ sensitivity curve is shown, resulting from injection and recovery tests of fake sources at different distances from the main source. K2-133 is clearly seen in centre and no extra sources are detected above the contrast curve, i.e. to a contrast of 7.7 mag at $0.5^{\prime\prime}$ and 8.4 mag further than $1^{\prime\prime}$.

2.2 Spectroscopic follow-up

We obtained zJ and HK NIR spectra of K2-133 using LIRIS (Manchado et al., 1998) at the 4.2m William Herschel Telescope. We also obtained spectra of HD 27267, a nearby rapidly-rotating A0V star for use as a telluric standard. The observations were taken in service mode on October 09 2017 (proposal SW2017b04). K2-133 and HD 27267 are located $5.4^{\circ}$ apart on the sky, resulting in a ca. 0.03 difference in airmass between the targets at the time of observation. We used the ABBA nod pointing pattern with exposures of 25 and 12 seconds for K2-133 for the zJ and HK grisms, respectively. Four exposures were taken for each grism and pointing position. The spectra were reduced and combined using the lirisdr iraf package222https://github.com/jaacostap/lirisdr, and wavelength calibrated using spectra of Argon arc lamps. The K2-133 spectra were then corrected for tellurics using the extracted A0V spectra and the idl routines of Vacca et al. (2003).

While bulk tellurics were corrected, some features were left in the resulting spectrum, possibly due to slightly differing airmass between the target and comparison star or atmospheric changes between exposures. We therefore masked regions of leftover telluric features in the corrected spectrum of K2-133 ($1.11-1.16$, $1.34-1.48$ and $1.80-1.92$ $\mu m$) to not interfere with our interpretations. The corrected spectrum can be seen in Fig. 2, split into three parts bounded by two of the masked telluric regions.

The corrected K2-133 spectrum was compared to template M-dwarf spectra from the IRTF Spectral Library333http://irtfweb.ifa.hawaii.edu/~spex/IRTF_Spectral_Library/index.html (Cushing et al., 2005; Rayner et al., 2009). The template spectra were first binned to the same wavelength using the pysynphot Python package (STScI development Team, 2013). Then each template was shifted and scaled to best match the K2-133 spectrum, and finally a $\chi^{2}$ test was conducted between them. The result of the $\chi^{2}$ tests can be seen in Fig. 3, showing increased agreement of the data with templates of spectral types M2/M1. This is close to the M1V spectral type found in Wells et al. (2018) and is consistent with the closer distance from Gaia (see Section 2.3). The IRTF template spectrum with the lowest $\chi^{2}$ was HD 42581 – an M1V star which is over-plotted on Fig. 2. An almost equally good $\chi^{2}$ was given by the M2V HD 95735 template.

2.3 Gaia DR2

Properties of the host star are further constrained by a precise parallax measurement from Gaia DR2 (Gaia Collaboration et al., 2018), providing a distance of $75.2\pm 0.2$ pc (Bailer-Jones et al., 2018). This is consistent with the distance of $78\pm 7$ pc reported in Wells et al. (2018), but far more precise. The slightly closer distance to K2-133 supports the analysis of the NIR spectrum, suggesting the star is closer to an M2 spectral type than previously reported.

We used the new distance to determine revised values of the stellar radius, mass and luminosity. We first calculated the absolute K-band magnitude using the apparent 2MASS K magnitude of $10.279\pm 0.018$ (Skrutskie et al., 2006) and the Gaia-derived distance. We then found the bolometric correction, BCK, using Equation 10 from Mann et al. (2015) (note the erratum Mann et al. 2016) with the apparent J and SDSS r magnitudes. We derived the bolometric luminosity directly from the absolute K mag and BCK, assuming Solar metalicity and no K-band extinction. The effective temperature was computed using the relations of Mann et al. (2015), and the stellar mass from Mann et al. (2019). The stellar radius was computed from the derived bolometric luminosity and effective temperature. The updated stellar values are given in Table 1, with uncertainties given at $1\sigma$ confidence level from an MCMC analysis.

We verified these stellar properties by comparing them to stellar models. We utilised the isochrones Python package (Morton, 2015a), using the MIST stellar isochrones (Choi et al., 2016; Dotter, 2016), and the 2MASS JHK magnitudes and Gaia parallax as inputs. This resulted in stellar parameters of T${}_{\mathrm{eff}}=3770_{-75}^{+104}$ K, R${}_{\star}=0.440_{-0.013}^{+0.010}$ R*☉, M${}_{\star}=0.461\pm 0.014$ M☉* and L${}_{\star}=0.0351_{-0.0015}^{+0.0019}$ L*☉*. These are in good agreement with the values from the empirical relations, in particular, both the mass and radius agree better than $1\sigma$.

2.4 Transit fitting

We have computed refined properties of the transiting planets, using the new stellar properties with the batman transit modelling code (Kreidberg, 2015; Mandel & Agol, 2002) and emcee (Foreman-Mackey et al., 2013; Goodman & Weare, 2010). We also fit for the limb-darkening coefficients using the three-parameter law and the Kipping (2016) sampling method. We note that the limb-darkening coefficients are not well defined, but they are important when estimating uncertainty in the inclinations of the planets. Since the original discovery, the K2 team has begun an effort to reduce each campaign with a uniform version of the reduction pipeline. We therefore use the reprocessed lightcurve in our analysis, with systematics removed using a modified version of the Lightkurve (Vinícius et al., 2018) implementation of the k2sff method (Vanderburg & Johnson, 2014).

3 Validation of planet e

With our refined stellar properties and lightcurve, we revisited the planetary candidate transit signal of K2-133 e. We note the transit is now detected at a ${\sim}13\sigma$ significance, leading us to evaluate the probability that it is caused by a fourth planet in the system and is not due to random or correlated noise. In order to validate the candidate, we require a false-positive probability (FPP) of less than one per cent, which is typically used to validate transit events.

We first note that K2-133 was located on Module 15, Output 2 (Channel 50) during observations. This channel has no reported detector irregularities from the K2 engineering test (Van Cleve & Caldwell, 2016). We also inspected the background counts (pixels outside the aperture) for K2-133 and two nearby stars – EPIC 247894072 and EPIC 247876252. We found no disagreements between background levels of the three stars and no signs of variability or trends, such as rolling band effects, which could affect the detected transit events.

3.1 Detrending

Our first test of the planet candidate was to use multiple detrending methods and see whether the transit signal was apparent in each detrended lightcurve. We utilised the k2sc algorithm (Aigrain et al., 2016), the everest high-level science product (Luger et al., 2018) and our SFF lightcurve from Section 2.4. We note that the everest lightcurve is from detrending the original (non-reprocessed) K2 data. We find the same transit feature in all three lightcurves with good agreement between each method as seen in Fig. 4, indicating the transit signal is not due to a detrending artefact. For the rest of our analysis we use the SFF lightcurve because it has the best 6.5-hour combined differential photometric precision (CDPP) of 145 ppm and the least visible correlated noise. We further removed data points corresponding to an asteroid passing through the aperture at time (BJD$-2454833$) ca. 3018 days and any exposures with bad Kepler quality flags.

3.2 Injection test

Our next goal was to test the sensitivity of our detection pipeline. To do this, we injected planetary transits into the full, non-detrended lightcurve with all four planets removed. The injected transits covered the parameter space of: 1-40 day orbital periods and transit depths of 400-1400 ppm. For reference, the candidate has an orbital period of 26.6 days and a depth of 1000 ppm. We used the batman transit model (Kreidberg, 2015) with a super-sampling factor of 15, varying the parameters of orbital phase, orbital period, planetary radius and orbital inclination. The semi-major axis was computed using Kepler’s third law. The eccentricity was fixed to zero, as this parameter does not largely effect K2 lightcurves due to the time sampling of 30 minutes. The stellar radius, mass and quadratic limb darkening coefficients ($\mu_{1}=0.5079,\mu_{2}=0.2239$ from Sing 2010) were also kept fixed for each injected model. Each injected planet had a minimum of three transits present in the lightcurve.

To detect planets, we utilised the Box-fitting Least Squares (BLS) algorithm of (Kovács et al., 2002) which searches for periodic box-shaped dips in a light curve. A number of trial periods are tested and each return a signal detection efficiency (SDE) signal-to-noise statistic. To recover an injected planet, we required the BLS SDE to be above 8 and the orbital period to be identified within 1% of the injected value. This is very similar to the tests performed in Rizzuto et al. (2017). We ran 64,000 injection tests which are summarised in Fig. 5, showing the recovered percentage as a function of the injected transit depth and orbital period. There is a clear trend of increased recovery rate for shorter periods and larger transit depths, which greatly impact the BLS algorithm. Planets b, c and e are marked on the figure and reside in regions of 70-90%. Planet d has a transit depth of 2000 ppm and orbital period of 11 days, placing it well into the 90+% region. We therefore conclude that the K2 data is of high enough quality that we should be able to detect transit signals of planet e.

3.3 Bootstrap tests

After deciding our detection pipeline was sufficiently sensitive to the properties of the planetary candidate, we next tested the probability that the transit signal could be caused by random or correlated noise. To do so, we first computed the SDE of the candidate transit signal in our refined lightcurve. We used the PyBLS444https://github.com/hpparvi/PyBLS implementation of the BLS algorithm using a period range of $1-40$ days, a transit duration to period ratio range of $0.001-0.2$, 50,000 frequency bins and no binning of the phase folded lightcurves. We then subtracted the running median from the resulting power spectrum using a window size of 1500 points, resulting in a peak SDE value of $10.7$. We then removed the transit signal by subtracting our best-fit transit model, resulting in a lightcurve composed of a combination of red and white noise. The highest BLS SDE of this lightcurve was $4.8$, far below any detection limit. This lightcurve was used for all false positive tests.

A transit signal can be mimicked by both white and red noise, so we first checked how the level of correlated noise compared to the Gaussian white noise. We used the mc3 Python package (Cubillos et al., 2017) to compute the RMS of our data and the RMS expected from only white noise for a range of bin sizes. The result is plotted in Fig. 6, where the red noise is roughly 40% above the white noise. I.e. our lightcurve is white noise dominated but contains a non-negligible amount of red noise as well. Consequently, we decided to test the candidate transit signal against red and white noise separately.

To test against white noise, we implemented a bootstrap test. For each test in the bootstrap, we randomly sampled the flux value at each cadence from the whole light curve, while keeping the times and flux errors the same. We repeated this 21,000 times, for each re-sample searching for periodic signals using the BLS algorithm with the same parameters as for the candidate. The resulting BLS SDE distribution is shown in Fig. 7, over-plotted with a kernel density estimation (KDE). We calculated the cumulative probability of the KDE at the SDE value of the candidate (10.7) and then the FPP was found by subtracting this value from unity, giving a FPP of 0.017%.

To test against correlated noise, we repeated this process with a different sampling method, this time using 50,000 re-samples. We opted for the bootstrapping method of Efron & Tibshirani (1991), used in the analysis of Roche tomograms by Watson & Dhillon (2001). In short, we randomly selected data points from the entire light curve, placing them at their original positions in the new light curve. The flux uncertainty values were then reduced by the square-root of the number of times the data point was selected. I.e. if a data point was selected $n$ times the associated error value was reduced by $\sqrt{n}$; if a point was not selected at all then it was removed from the lightcurve. This method should show up correlated noise as periodic signals, given enough samples. The resulting distribution is again shown in Fig. 7, over-plotted with a kernel density estimation (KDE). The FPP was found to be $2.7\times 10^{-6}\%$.

3.4 Astrophysical false positives

We also tested the transit signal against eclipsing binary scenarios, using the VESPA Python package (Morton, 2015b, 2012). VESPA computes likelihoods of the signal being caused by eclipsing binaries (EBs), hierarchical eclipsing binaries (HEBs) and background eclipsing binaries (BEBs), as well as planetary transits. Scenarios are also considered at twice the candidate orbital period. VESPA then excludes systems that are inconsistent with the input data.

The SDSS, 2MASS and Kepler magnitudes alongside the Gaia parallax were used as inputs for the stellar populations. The fitted orbital period and planet-to-star radius ratio, transit light curve, Keck contrast curve and stellar RA and Dec were used as constraints to compute the likelihood of each scenario. We used a simulation sample size of 100,000 for each population and ran 10 bootstrap resamples. We obtained mean probabilities of 0.08, 2.28, 0.00 and 97.65 per cent for the EB, HEB, BEB and planetary scenarios, respectively. This equates to a FPP of 2.35 per cent. Lissauer et al. (2012) have shown that false-positives are less-likely in multi-planet systems, and the false-positive probability of a transit signal in a $2+$ planetary system is decreased by a factor of approximately 50. Applying this to our candidate reduces the FPP to ca. 0.05%.

With both noise and astrophysical false-positive probabilities well below one per cent, we confirm the planet candidate, now named K2-133 e.

3.5 Properties

We evaluated properties for the newly validated planet by running batman/emcee with wide priors on the transit parameters. The phase-folded K2 data and MCMC-fitted transit model are shown in Fig. 8, and the fitted transit parameters are given in Table 3. All four planet transits were fit simultaneously, allowing the limb-darkening coefficients to vary. The MCMC samples of fitted parameters are shown in Fig. 9. Of interest are the planetary radius of 1.7 R*⊕* and the orbital distance from the host star, which are discussed in terms of potential habitability in Section 4.1. We note the lowest flux point has a rather large effect on the fitted parameters – masking this point leads to a more flat-bottomed transit with a $\chi^{2}$ value of approximately 1. However, we do not exclude any data points in our analysis.

4 Discussion

4.1 Temperate Zone

Planet e orbits the host star at a distance of 0.14 AU (64 stellar radii), meaning it receives roughly 1.8 times the stellar flux that Earth does. The $2\sigma$ confidence of the derived radius means the planet could be equally likely to be rocky or gaseous (1.6 $R_{\earth}$) to very likely to be gaseous (2.0 $R_{\earth}$) (Rogers, 2015). Therefore the first test to the planets habitability would be to constrain the mass, which will be challenging (see Section 4.3).

Continuing, assuming a rocky planet, we evaluated the stars “habitable” zone boundaries using the parametric equations of Kopparapu et al. (2014). We find that the planet orbits just interior to the optimistic habitable zone set by the “recent Venus” criterion. However, the boundaries predicted by different models vary considerably depending on assumptions and the inclusion of clouds.

[TABLE]

We calculated the Earth Similarity Index555http://phl.upr.edu/projects/earth-similarity-index-esi (ESI; Schulze-Makuch et al. 2011) for planet e using Equation (1), where $S$ is the incident stellar flux and $R$ is the planetary radius, both relative to Earth. We find a ESI score of 0.72, which places the planet in the top-25 planets most similar to Earth in the Optimistic Sample of the Planetary Habitability Laboratory666http://phl.upr.edu/projects/habitable-exoplanets-catalog.

Due to the tight orbit and likely >1 Gyr age of the system, K2-133 e is expected to be tidally locked to the host star (Grießmeier et al., 2009) which may reduce the prospect for life. The K2 light curve shows no evidence of stellar flares – which are capable of stripping atmospheres. The radii of the planets, however, do span the radius gap hinting that their atmospheres may have been eroded in the past and would therefore be uninhabitable. This is discussed more in Section 4.2.

4.2 Radius gap

A radius gap has been identified in the distribution of known exoplanets between 1.5 and 2 $R_{\earth}$ (Fulton et al., 2017; Fulton & Petigura, 2018). The slope of this gap as a function of orbital period has been found to be negative, which is consistent with the gap caused by photo-evaporation rather than being inherent in planet formation (Van Eylen et al., 2018). This would mean planets are formed with gaseous envelopes and then are either stripped to their rocky cores by high stellar X-ray/UV irradiation (shorter periods), or retain their atmospheres in a lower radiation environment (longer periods) while the star becomes less active as it ages (Lopez & Fortney, 2013; Owen & Wu, 2013). The core mass of a planet is also important for retaining an atmosphere, as heavier planets will suffer lower mass-loss rates.

The inner three planets in the K2-133 system span the radius gap, with radii of 1.3, 1.6 and 2.0 $R_{\earth}$. The trend of increasing radius with decreasing incident flux hints that the system may have been sculpted by photo-evaporation when the star was young and more active. While the radius of the fourth planet is not well constrained ($2\sigma$ confidence $1.7\pm 0.3$ R*⊕*), the radius-period trend of the inner planets indicates that it may have been sufficiently far from the star to keep its primordial atmosphere. However, another explanation could be that the planets have different core masses.

K2-133 is comparable to the GJ-9827 system – three super-Earths orbiting a late K star with radii 1.6, 1.2 and 2.0 $R_{\earth}$ (Niraula et al., 2017), which has been studied in great detail with precise radial velocity measurements (Teske et al., 2018; Prieto-Arranz et al., 2018; Rice et al., 2019). A similar follow-up study of K2-133 would be interesting as a comparison to the GJ 9827 system. The induced RV shifts should be similar (${\sim}{2}{+}$ m s*-1*), however the star is much fainter (V = 14).

4.3 Follow-up prospects

The system is challenging for follow-up observations, due to the planets producing shallow transits and weak RV shifts. The star is reasonably bright towards the IR (K = 10) but fainter in the optical (V = 14), making ground observations particularly difficult. Ideally follow-up observations would provide information on the planetary masses and atmospheres, plus a stronger constraint on the radius of planet e or possibly identify additional planets in the system.

It is possible to infer masses from transit timing variations (TTVs) (Grimm et al., 2018; Jontof-Hutter et al., 2016; Nesvorný et al., 2013). This involves precisely measuring the central transit times of individual events which can vary by seconds to hours due to gravitational interactions between the planets. Using the mass-radius relation of Weiss & Marcy (2014), we estimate the mid-transit times of the K2-133 planets to vary by 5-40 seconds using the TTVFaster code (Agol & Deck, 2016), however these amplitudes could be larger due to heavier planets or the presence of more planets in the system. We do not identify any significant TTVs in the K2 data, however due to the cadence we are limited to a $1\sigma$ transit timing precision of 2.2+ minutes.

K2-133 (TIC 150096001) lies close to the ecliptic plane with an ecliptic latitude of $2.8^{\circ}$. This means it will not be observed in TESS Cycle 2 which will view a minimum latitude of $6^{\circ}$. However, TESS may still observe the ecliptic in its extended mission. K2-133 has an estimated TESS mag of 12.3 in the TESS Input Catalog (Stassun et al., 2018), which corresponds to a 1-hour noise level of roughly 750ppm (Barclay, 2017). Individual transits should therefore be detected at ca. 1.5-3 sigma, with three of the planets transiting more than once. New transit data could also be acquired with the forthcoming PLATO mission (Rauer et al., 2014), which is planned to launch in 2026. K2-133 may be observed in one of the stop-and-stare pointings, which will be defined 2 years before launch (Rauer & Heras, 2018). Unfortunately K2-133 is too faint for observations with CHEOPS, which will have an approximate limiting V-band magnitude of 12 (Broeg et al., 2013).

5 Conclusions

We have validated a fourth planet in the K2-133 system. This planet has a radius of $1.73_{-0.13}^{+0.14}$ R*⊕* and orbits on the inner edge of the stars optimistic habitable zone. We have also presented updated properties of the host star and three already validated planets. The system warrants further study to test theories of the cause of the observed radius gap in the exoplanet population, such as photo-evaporation. Planetary mass measurements are needed for this, which will be a challenge to obtain due to the faintness of the host star.

Acknowledgements

We thank the reviewer for helpful comments to improve the manuscript. We also thank Dr. Ian Crossfield and Dr. David Ciardi for allowing the use of their AO imaging observations and analysis. This research made use of the WHT and its service programme, which are operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. This research made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. RW acknowledges funding from the Northern Ireland Department for Education. KP acknowledges funding from the Leibniz Association through its 2019 Leibniz Competition Program. CAW acknowledges support from the STFC grant ST/P000312/1.