NGTS-13b: A hot 4.8 Jupiter-mass planet transiting a subgiant star

NGTS-13 was observed using a single NGTS telescope from UT 2017 December 10 to UT 2018 August 7. During this period 255,709 observations were taken with exposure times of 10 seconds. Image reduction and photometry were performed with the NGTS pipeline (Wheatley et al., 2018) using CASUtools software¹¹1http://casu.ast.cam.ac.uk/surveys-projects/software-release. The lightcurve was detrended and systematic effects were removed using the SysRem algorithm (Tamuz et al., 2005).

The NGTS lightcurves were searched for periodic transit-like signals using an implementation of the box-fitting least-squares (BLS) algorithm (Kovács et al., 2002). We detected a strong peak in the BLS periodogram at 4.12 days and folding the NGTS photometry on this period revealed a clear transit signal.

We binned the NGTS data to 2-minute bins for our analysis as described in Section 3.2. NGTS-13 has some nearby companions (see Figure 1) and the NGTS photometry is expected to suffer from some dilution, which we discuss in Section 3.2.2. The NGTS photometry binned to 30 minutes for display purposes is shown in Figure 2. A full table of reduced NGTS photometric observations is available in a machine-readable format from the online journal.

TESS observed NGTS-13 (TIC 454069765) in Sector 10 from UT 2019 March 28 to UT 2019 April 21 with 30-minute cadence Full Frame Images (FFIs). We created a lightcurve via aperture photometry of the TESS 30-minute cadence FFIs following the method set out in Gill et al. (2020). We determined a threshold for target pixels and for background pixels based on an iterative sigma-clipping to determine the median and standard deviation of the background pixels. To exclude the neighboring T = 12.8 star (TIC 454069772, 57” southeast, see Figure 1), we only include pixels in our final aperture for which neighbouring pixels closer to the center of our target show a higher illumination. For the final pixel mask of NGTS-13, 7 pixels were selected (see Figure 1). To remove the systematic trends from the TESS light curve, we mask out the transits, and then fit a linear spline interpolating across the positions of the transits. Before removing systematic trends, we estimate a signal-to-noise ratio (SNR) as the transit depth relative to the standard deviation of out-of-transit observations and find a single TESS transit SNR$\sim$5.4, which corresponds to SNR$\sim$13.3 when considering all 6 transits in the TESS data (SNR $\propto$ $\sqrt{N_{\mathrm{T}}}$). The TESS photometry is displayed in Figure 2. A full table of reduced TESS photometric observations is available in a machine-readable format from the online journal.

We obtained 13 observations of NGTS-13 with the high resolution echelle spectrograph CORALIE on the Swiss 1.2 m Euler telescope (Queloz et al., 2001) at La Silla Observatory, Chile between UT 2019 December 12 and UT 2020 January 23. CORALIE has a resolution of $R$ $\sim$ 60,000 and is fed by two fibers: a 2 arcsec on-sky science fiber encompassing the star and another fiber that can either connect to a Fabry-Pérot etalon for simultaneous wavelength calibration or on-sky for background subtraction of sky flux. We observed NGTS-13 in the simultaneous Fabry-Pérot wavelength calibration mode using an exposure time of 2700 seconds. The spectra were reduced with the CORALIE standard reduction pipeline and RVs were computed for each epoch by cross-correlating with a binary G2 mask (Pepe et al., 2002). We co-added the 13 CORALIE spectroscopic observations by correcting each epoch to the stellar rest frame with its corresponding barycentric RV to create a single high signal-to-noise spectrum. We list the CORALIE RVs in Table 1. The RVs phased to the 4.12 day period are displayed in Figure 3.

We used wavelet analysis to extract stellar atmospheric parameters from the co-added 13 CORALIE spectroscopic observations of NGTS-13 using the methodology set out in Gill et al. (2018, 2019). With the wavelet analysis we find NGTS-13 to have surface rotational velocity $v\sin i_{*}$ = 6.2 $\pm$ 1.2 km s^-1. We also derived stellar atmospheric parameters using SpecMatch-Emp (Yee et al., 2017). SpecMatch-Emp uses a large library of stars with well-determined parameters to match the input spectra and derive spectral parameters. We use a spectral region that includes the Mg I b triplet (5100 - 5400 $\AA$) to match our spectra. SpecMatch-Emp uses $\chi^{2}$ minimisation and a weighted linear combination of the five best matching spectra in the SpecMatch-Emp library to determine $T_{\rm eff}$, $\log g_{*}$, and $\left[{\rm Fe}/{\rm H}\right]$. We combine the wavelet analysis results and SpecMatch-Emp results to create wide uncertainties that include the uncertainty range of both methods to use as priors for their final derivation as described below. We present these spectroscopically derived parameters from both methods and the adopted combined values in table 2.

We derive both planet and stellar parameters using EXOFASTv2 (Eastman et al., 2013; Eastman, 2017; Eastman et al., 2019), which can globally fit all photometry and RV data. A full description of EXOFASTv2 is given in Eastman et al. (2019), which can fit any number of transits and RV sources while exploring the vast parameter space through a differential evolution Markov Chain coupled with a Metropolis-Hastings Monte Carlo sampler. A built-in Gelman-Rubin statistic (Gelman & Rubin, 1992; Gelman et al., 2003; Ford, 2006) is used to check the convergence of the chains.

We use Modules for Experiments in Stellar Astrophysics (MESA) Isochrones & Stellar Tracks (MIST; Choi et al., 2016; Dotter, 2016) isochrones and the spectral energy distribution (SED) within the EXOFASTv2 fit to determine host star parameters. For the SED fit within EXOFASTv2 (Stassun & Torres, 2016) we use photometry from APASS DR9 BV (Henden et al., 2016), 2MASS JHK (Skrutskie et al., 2006), SDSS DR12 gri (Alam et al., 2015), ALL-WISE W1, W2 and W3 (Wright et al., 2010), and GALEX NUV (Bianchi et al., 2011), which are presented in Table 3. The SED fit within EXOFASTv2 puts sytematic floors on the broadband photometry errors (Stassun & Torres, 2016). The EXOFASTv2 SED fit is presented in Figure 4. We use our spectroscopically derived $T_{\rm eff}$, $\log g_{*}$, and $\left[{\rm Fe}/{\rm H}\right]$ values as priors in the global model as well as a Gaussian prior for parallax from the value and uncertainty in $Gaia$ DR2. We put an upper limit on the extinction based on the maximum in line-of-sight value from the Schlafly & Finkbeiner (2011) galaxy dust map. All other fitted and derived parameters from our EXOFASTv2 model have conservative physical boundaries that are detailed in Table 3 of Eastman et al. (2019), which also gives a thorough explanation of each parameter.

Due to the long exposure time of the TESS full frame exposures we must account for smearing of the light curve. Therefore, within EXOFASTv2 we specify an exposure time of 30 minutes and average over 10 data points to integrate a model over the exposure time equivalent to a midpoint Riemann sum (Eastman et al., 2019). We assign the TESS FFI photometry to the built-in EXOFASTv2 ”TESS” band and the NGTS photometry to the ”R” band for the computation of the limb darkening coefficients. We found the EXOFASTv2 ”R” band to have the most similar bandpass to NGTS in the EXOFASTv2 framework, which only allows filters defined in Claret & Bloemen (2011). We tested the use of this filter by running a fit without TESS data and NGTS data assigned to the ”TESS” filter; we obtained limb darkening parameters within uncertainties to those presented in Table 4. We also ran a fit including TESS data assigned to the ”TESS” filter and the NGTS data assigned to the ”Kepler” filter; we again found limb darkening parameters within uncertainties to those presented. We account for dilution of the photometry lightcurves, which we describe in Section 3.2.2. For our EXOFASVTv2 global model and Markov chain Monte Carlo (MCMC) fit we used 56 walkers (2 $\times$ n${}_{\text{parameters}}$), or chains, and allowed the fit to run until convergence (132776 steps).

Previous studies have suggested that 4 $\mathrm{M_{\rm Jup}}$ could be an approximate border for a distinct massive giant planet population compared to lower-mass giant planets. Santos et al. (2017) find that massive giant planets tend to orbit stars with $\left[{\rm Fe}/{\rm H}\right]$ distributions statistically similar to field stars, unlike lower-mass giant planets which occur more frequently around metal-rich stars. This could be evidence for two different formation mechanisms where lower-mass giants may form via core accretion (Pollack et al., 1996) and higher mass giants may form via disk instability (Boss, 1997), as formation via disk instability is not as metallicity dependent as core-accretion (Boss, 2002; Cai et al., 2006).

Additionally, Schlaufman (2018) find that substellar companions with masses $\lesssim$ 4 $\mathrm{M_{\rm Jup}}$ preferentially orbit metal rich stars (a property associated with core accretion), while companions with masses $\gtrsim$ 10 $\mathrm{M_{\rm Jup}}$ do no share this property. Maldonado et al. (2019) also find that massive giant planets tend to have host stars with lower metallicities than lower-mass giant planets and conclude that the core-accretion planet formation mechanism achieves its maximum efficiency for planets with masses in the range 0.2 – 2 $\mathrm{M_{\rm Jup}}$. Goda & Matsuo (2019) conclude that the mean metallicity of stars hosting companions between 4 and 25 $\mathrm{M_{\rm Jup}}$ is lower than that of the sample of companions with masses 0.3 to 4 $\mathrm{M_{\rm Jup}}$, and the mean metallicity of stars with companions more massive than 25 $\mathrm{M_{\rm Jup}}$ is much lower than those of the other two sub-samples. Goda & Matsuo (2019) use 25 $\mathrm{M_{\rm Jup}}$ as the lower mass limit of brown dwarfs (instead of the deuterium burning limit) because they found the mean metallicity of G-type host stars with objects lighter than 25 $\mathrm{M_{\rm Jup}}$ to be much higher than those with objects more massive than 25 $\mathrm{M_{\rm Jup}}$; they suggest 25 $\mathrm{M_{\rm Jup}}$ is the upper mass limit for core-accreted planets. However, other previous studies suggest lower mass brown dwarfs (13 - 42.5 $\mathrm{M_{\rm Jup}}$) form via disk instability, e.g., Ma & Ge (2014).

Adibekyan (2019) analyzes the known giant planet population and his results do not support that 4 $\mathrm{M_{\rm Jup}}$ is a transition point between two separate formation channels, but his results do suggest that high mass planets can form through different mechanisms depending on their initial environment. With a metallicity of $\left[{\rm Fe}/{\rm H}\right]$ = 0.25 $\pm$ 0.17 NGTS-13b’s host star is metal-rich and does not provide additional evidence that more massive-giant planets prefer a different formation mechanism than core accretion. We put NGTS-13b in this context in Figure 6 which shows the metallicity of giant planets from 0.5 - 13 $\mathrm{M_{\rm Jup}}$.

We put NGTS-13b in terms of its density, $\rho_{P}$ = 4.02 $\pm$ 0.55 g cm^-3, in Figure 7 by displaying it in the mass-radius diagram with other giant planets and brown dwarfs as well as density as a function of mass. We color symbols in the mass-radius diagram based on equilibrium temperature to show the effect high equilibrium temperatures can have on radii. We can see in Figure 7 the larger maximum radii for planets with masses lower than $\sim$4 $\mathrm{M_{\rm Jup}}$, but this is explained by the redder coloring of these symbols indicating higher equilibrium temperatures (resulting in more inflation). NGTS-13b has a typical radius compared to planets with similar masses suggesting that even with its high equilibrium temperature it is not significantly affected by bloating (inflated radius due to a planet’s atmosphere being heated from stellar irradiation and expanding), which is expected for massive planets (Sestovic et al., 2018). The downward trend in radii (and increasing trend in density) with larger mass in Figure 7 highlights the importance of the larger surface gravity in more massive objects.

NGTS-13 has an effective temperature of $T_{\rm eff}$ = 5819 $\pm$ 73 K, suggesting it is a G2 type star (Pecaut & Mamajek, 2013). However, the rather large stellar radius of $R_{*}$ = 1.79 $\pm$ 0.06 $\mathrm{R_{\odot}}$ indicates NGTS-13 has evolved slightly as this radius is larger than expected for a main sequence star of this temperature, and its surface gravity $\log g_{*}$ = 4.04 $\pm$ 0.05 suggests the star is a subgiant ($\log g_{*}$ $<$ 4.1). Our MIST isochrone analysis (see Section 3.2) finds an equivalent evolutionary point (EEP) value of 413${}^{+47}_{-18}$, which still puts NGTS-13 in the main sequence evolutionary phase, as the terminal age main sequence EEP does not start until EEP = 454, see Choi et al. (2016), Dotter (2016), and the MIST documentation⁴⁴4http://waps.cfa.harvard.edu/MIST/README_tables.pdf. Additionally, it is less likely in terms of evolutionary timescales to observe a star near or past the turnoff point than to observe it in the middle of the main sequence, therefore NGTS-13 is likely still fusing Hydrogen in its core but has likely begun the transition from the main sequence to the red giant branch placing it in the subgiant branch.

Lillo-Box et al. (2016) found a relative lack of hot Jupiters around giant and subgiant stars, suggesting that the close-in massive planets around main-sequence stars are engulfed by the star as it evolves. However, other studies do not find a lack of hot Jupiters around evolved hosts; additionally NGTS-13 is likely still burning hydrogen in its core and has not evolved long enough to engulf planets near the orbital distance of NGTS-13b. Grunblatt et al. (2019) found that low-luminosity red giant branch evolved stars have a tentatively higher population of close-in giant ($R_{P}$ $>$ 1 $\mathrm{R_{\rm Jup}}$) planets than main-sequence stars suggesting that stellar evolution does not affect close-in giant planet occurrence significantly until stars are substantially in the red giant branch with radii of 5-6 $\mathrm{R_{\odot}}$. Zhou et al. (2019) studied the occurrence rate of hot Jupiters as a function of stellar mass and did not find any statistically consistent trends with stellar mass and that hot Jupiters are just as abundant around main-sequence A stars as they are around F and G stars.

The $M_{*}$ = 1.30${}^{+0.11}_{-0.18}$ $\mathrm{M_{\odot}}$ stellar mass of NGTS-13 is larger than a typical G-type main-sequence star, which suggests it evolved to a cooler temperature and is a ”retired” F star. A relatively larger host star mass is consistent with previous studies that found massive planets to be more common around massive stars (Santos et al., 2017; Maldonado et al., 2019). Santos et al. (2017) find that more massive giant planets are more common around more massive hosts, which are often more evolved; they suggest the two distinct populations of giant planets could be related to evolved stars not showing a clear metallicity-giant planet correlation.

We display host star mass by color in Figure 7. Only one host star has a mass below 0.9 $\mathrm{M_{\odot}}$ in our sample of 17 well-defined massive giant planets. Grunblatt et al. (2018) found that close-in giant planets around evolved stars tend to have more eccentric orbits than those around main sequence stars, but they focus on giant stars and orbital periods longer than 4.5 days. We explore NGTS-13b’s eccentricity further in Section 4.4.

NGTS-13b has a slightly eccentric, $e$ = 0.086 $\pm$ 0.034, orbit with a short period of 4.12 days. Analyzing the posterior distribution of the eccentricity displayed in Figure 8 we see that it is not centered around 0. We performed the test of Lucy & Sweeney (1971), and find the statistical significance of the eccentric fit to be P$(e>0)$ = 0.9591, which just passes the 5% significance level suggested by Lucy & Sweeney (1971).

To further test the eccentric model we fit a Keplerian to the 13 RVs with priors on the time of conjunction and period and compare fits with and without forcing $\sqrt{e}\cos\omega_{*}$ and $\sqrt{e}\sin\omega_{*}$ to 0. We also allow an RV jitter term to vary, which corresponds to $k$ = 7 free parameters for the eccentric fit and $k$ = 5 for the circular fit. We determine the Bayesian Information Criterion (BIC), where given the free parameters $k$, number of measurements $n$, and maximized value of the likelihood function of the model $\mathscr{L}$:

When picking from several models the one with the lowest BIC is preferred. We find a $\Delta$BIC = 7.17 showing the eccentricity model is moderately favored over the circularized model. For this test we found an RV jitter of 31 m s^-1 for the circular fit but 0 m s^-1 for the eccentric fit as no additional RV uncertainty is needed for the eccentric model.

We note that for our final EXOFASTv2 analysis presented the model has a RV jitter variance $\sigma_{J}^{2}$ = $-800^{+2400}_{-1100}$ which essentially reduces the error size of the RVs when evaluating the model. As noted in Equation 13 of Hara et al. (2019), the uncertainty in the eccentricity is approximately proportional to the error of the RV measurement; therefore, this added RV jitter variance may cause an underestimation of our eccentricity uncertainties. Hara et al. (2019) show that the eccentricity estimate can be affected by an undetected correlated signal. Given the amplitude of the RV signal of NGTS-13b, we study the possibility that an undetected outer planetary companion would cause a nonzero eccentricity. We test this by adding a linear drift to our global EXOFASTv2 model. When allowing for an unconstrained linear RV drift we find a slope of only 1.1${}^{+1.2}_{-1.4}$ m s^-1 day^-1 and an eccentricity $e$ = 0.084${}^{+0.023}_{-0.021}$, similar to our presented model with no RV drift. The similar eccentricity suggests that NGTS-13b’s eccentricity measurement is not caused by the model attempting to account for an undetected planet.

The slight eccentricity of NGTS-13b is not uncommon for massive planets (see Figure 9). Hansen (2010) found that hot giant planets more massive than 3 $\mathrm{M_{\rm Jup}}$ have the upper envelope of their eccentricity distribution shifted to lower periods compared to less massive hot Jupiters. Also Dawson & Johnson (2018) found that most hot Jupiters with periods $<$ 3 days have circular orbits, but some hot-Jupiters in the 3-10 day orbital period range occupy moderately eccentric orbits. Assuming a general tidal circularization effect for close-in planets from equation 3 of Adams & Laughlin (2006), NGTS-13b has a circularization timescale of 0.5 Gyr for a tidal quality factor $Q_{p}$ = 10⁵, 5.0 Gyr for $Q_{p}$ = 10⁶, and 50.3 Gyr for $Q_{p}$ = 10⁷. Given the 4.23${}^{+2.65}_{-1.59}$ Gyr age of the system, the current eccentricity may be an indication that NGTS-13b underwent dynamical interactions with other components in the system during its migration history, e.g., high-eccentricity tidal migration (see Section 3.1 of Dawson & Johnson (2018)). We put NGTS-13b’s eccentricity into context with other giant planets and transiting brown dwarfs in Figure 9.