Proxima Centauri b is not a transiting exoplanet

TL;DR

This study used Spitzer Space Telescope data to conclusively rule out transits of Proxima Centauri b at the 4.5 micron wavelength, refining its size constraints and suggesting focus on other observational methods.

Contribution

First infrared transit observations of Proxima Centauri b that set stringent upper limits on its radius and dismissed previous transit claims.

Findings

No transits detected at 4.5 microns at 200 ppm level

Proxima b's radius upper limit is 0.4 Earth radii

Stellar activity is lower in infrared wavelengths

Abstract

We report Spitzer Space Telescope observations during predicted transits of the exoplanet Proxima Centauri b. As the nearest terrestrial habitable-zone planet we will ever discover, any potential transit of Proxima b would place strong constraints on its radius, bulk density, and atmosphere. Subsequent transmission spectroscopy and secondary-eclipse measurements could then probe the atmospheric chemistry, physical processes, and orbit, including a search for biosignatures. However, our photometric results rule out planetary transits at the 200~ppm level at 4.5$~{\mu}m$, yielding a 3$\sigma$ upper radius limit of 0.4~$R_\rm{\oplus}$ (Earth radii). Previous claims of possible transits from optical ground- and space-based photometry were likely correlated noise in the data from Proxima Centauri's frequent flaring. Follow-up observations should focus on planetary radio emission, phase curves, and direct imaging. Our study indicates dramatically reduced stellar activity at near-to-mid infrared wavelengths, compared to the optical. Proxima b is an ideal target for space-based infrared telescopes, if their instruments can be configured to handle Proxima's brightness.

Tables

Table 1. Table 1: Summary of the RV data sets.

Data set	Instrument	No. Observations	Baseline	Cadence
UVES	UVES	77	31/03/2000 $-$ 30/03/2007	Low
pre-2016	HARPS	63	27/05/2004 $-$ 23/03/2013	Low
PRD	HARPS	53	19/01/2016 $-$ 01/04/2016	High
RD	HARPS	62	01/06/2017 $-$ 30/09/2017	High

Table 2. Table 2: Orbital constraints and nuisance parameters for the Proxima b model from the EMPEROR analysis of RV data.

Orbital Parameter	Model Value
Amplitude [ms-1]	${1.32}_{-0.14}^{+0.12}$
Period [d]	${11.1855}_{-0.0014}^{+0.0016}$
Phase [rads]	${3.44}_{-1.79}^{+0.62}$
Longitude [rads]	${4.40}_{-3.47}^{+0.84}$
Eccentricity	${0.08}_{-0.06}^{+0.07}$
$\dot{\gamma }$ [ms-1d-1]	$-{0.0005}_{-0.0002}^{+0.0003}$
${\sigma }_{jit,pre-2016}$ [ms-1]	${1.84}_{-0.11}^{+0.09}$
${\gamma }_{pre-2016}$ [ms-1]	${1.23}_{-0.73}^{+0.67}$
${\varphi }_{pre-2016}$	${0.63}_{-0.11}^{+0.16}$
${\tau }_{pre-2016}$ [d]	${7.57}_{-3.10}^{+1.52}$
${\sigma }_{jit,PRD}$ [ms-1]	${1.44}_{-0.20}^{+0.10}$
${\gamma }_{PRD}$ [ms-1]	${1.98}_{-0.92}^{+1.08}$
${\varphi }_{PRD}$	${0.38}_{-0.11}^{+0.17}$
${\tau }_{PRD}$ [d]	${7.86}_{-4.46}^{+0.97}$
${\sigma }_{jit,RD}$ [ms-1]	${1.14}_{-0.10}^{+0.13}$
${\gamma }_{RD}$ [ms-1]	${2.10}_{-0.90}^{+1.21}$
${\varphi }_{RD}$	${0.50}_{-0.20}^{+0.16}$
${\tau }_{RD}$ [d]	${6.91}_{-2.58}^{+2.33}$
${\sigma }_{jit,UVES}$ [ms-1]	${1.91}_{-0.11}^{+0.09}$
${\gamma }_{UVES}$ [ms-1]	$-{0.20}_{-0.32}^{+0.27}$
${\varphi }_{UVES}$	${0.62}_{-0.17}^{+0.19}$
${\tau }_{UVES}$ [d]	${4.54}_{-1.63}^{+2.35}$

Equations

$$F=F_{s}Tr(t)M(x,y)R(t),$$

$$F=F_s\left(\sum _{i=1}^n{c}_i\widehat{P_i}+Tr(t)+R(t)\right),$$

$$N=\frac{{\left(\sum I(i)\right)}^2}{\sum I(i{)}^2}$$

$$F_o(\lambda )=c[F_t(\lambda )F_i(\lambda +\Delta \lambda )]\ast PSF,$$

$$m(t)=\sum _{i=1}^k\sum _{j=1}^d[K_{i,j}(cos({\omega }_{i,j}+T_{i,j}(t))+e_{i,j}cos{\omega }_{i,j})]+\dot{\gamma }+{\sigma }_{jit,j}+M{A}_j$$

Full text

††pagerange: Proxima Centauri b is not a transiting exoplanet–References††pubyear: 2019

1 Introduction

The search for the nearest small planets has accelerated in recent years with the development of purpose-built instrumentation (e.g., Mayor et al., 2003a; Crane et al., 2006; Cosentino et al., 2012; Pepe et al., 2013; Quirrenbach et al., 2018, amongst others). Some highlights include the multi-planet systems orbiting the nearby stars HD 69830 (Lovis et al., 2006), HD 10180 (Lovis et al., 2011; Tuomi, 2012), HD 40307 (Tuomi et al., 2013), or 61 Virginis (Vogt et al., 2010). The value of these nearby planetary systems significantly increases when the planets are found to transit their host stars, like those orbiting HD 219134 (see for example, Vogt et al., 2015; Motalebi et al., 2015; Gillon et al., 2017a), 55 Cancri e (Winn et al., 2011), or the recent Transiting Exoplanet Survey Satellite discovery of $\Pi$ Mensae c (Huang et al., 2018).

M-dwarf stars are particularly fruitful targets. Their occurrence rate for planets with masses below 10 $M_{\rm{\oplus}}$ (Earth masses) is at least one per star, with a habitable zone (HZ) occurrence rate in the mass range 3–10 $M_{\rm{\oplus}}$ of 0.21${}^{+0.03}_{-0.05}$ planets per star (Tuomi et al., 2014). Within this population are dense multi-planet systems like the seven planet candidates orbiting GJ 667C (Anglada-Escudé et al., 2013) and the four planets around GJ 876 (Rivera et al., 2010; Jenkins et al., 2014). Furthermore, Ribas et al. (2018) detected a small planetary candidate orbiting Barnard’s star.

M dwarfs host some spectacular transiting systems, particularly for the low-mass population, both inside and outside of the HZ. The super-Earth GJ 1214 b transits its host star (Charbonneau et al., 2009), allowing studies of its atmospheric composition (e.g., Bean et al., 2010; Rackham et al., 2017). The star LHS 1140 hosts two transiting planets (Dittmann et al., 2017; Ment et al., 2018). But, the standard-bearer in this class is the TRAPPIST-1 system, with seven small transiting exoplanets, at least three of which are in the HZ (Gillon et al., 2017b).

Proxima Centauri b (Anglada-Escudé et al., 2016, hereafter Proxima b) provides potentially the best possible opportunity for exoplanet characterization. Not only does the planet orbit our nearest stellar neighbour, but it has an equilibrium temperature and minimum mass similar to the Earth, meaning it could be rocky and have liquid surface water. If it transits, characterisation of its atmosphere and surface would be possible. Kipping et al. (2017) used optical photometry from the Microvariability and Oscillations of STars (MOST) Space Telescope. Although they reported a candidate signature matching the planet’s expected properties, they could not confirm the feature. Li et al. (2017) also reported a possible optical transit of Proxima b using data taken with the 30 cm telescope at Las Campanas, but again without confirmation. Using the Bright Star Survey Telescope in Antartica, Liu et al. (2018) report a number of transit-like events that could be ascribed to Proxima b, and assuming large transit timing variations (TTV), they could phase with the event reported in Kipping et al. (2017). On the contrary, Blank et al. (2018) were unable to confirm any of the previously reported transit events using optical data spanning 11 years, albeit with heterogeneous and non-continuous data sets. They do, however, confirm the impact of high stellar activity on optical light curves.

Observations at longer wavelengths mitigate the effects of stellar activity, increasing precision. In $\S$ 2 we discuss an observing campaign with the Spitzer Space Telescope at 4.5 m to search for the Proxima b’s transits. We then present new constraints on the mass provided by the latest radial-velocity (RV) data in $\S$ 3, both for Proxima b and any planets interior to its orbit. Finally, we summarise our findings in $\S$ 4.

2 Photometric Observations and Analysis

In November 2016, we observed Proxima Centauri for over 48 hours in the 4.5 m band of the InfraRed Array Camera (IRAC, Fazio et al., 2004) on the Spitzer Space Telescope (Werner et al., 2004), with target reacquisition roughly every 16 hours. We centered the observation on the predicted transit time of 2,457,708.02 $\pm$ 0.33 BJD, calculated from the original orbital solution for Proxima b published in Anglada-Escudé et al. (2016). The subarray mode frame time of 0.02 s resulted in $\sim 600,000$ individual frames (Spitzer Proposal ID 13155, PI: James Jenkins). Due to on-board data storage limits, there are short gaps between successive sets of 64 subarray images. The IRAC heater was off for the duration of the stare.

We used Basic Calibrated Data frames from Spitzer pipeline version S19.2.0. We performed twice-iterated 4$\sigma$ bad pixel rejection at every pixel position within each 64-frame set of subarray images to mask cosmic ray hits, and combined these with masks supplied by Spitzer.

Two groups within our team, RC and JH at UCF and NT, RR, and JJ at U. de Chile, analyzed all the data with completely independent codes, obtaining closely similar results. The UCF group used its Photometry for Orbits, Eclipses, and Transits pipeline (POET; Stevenson et al., 2012; Cubillos et al., 2014), while the Chilean group wrote a new code, in consultation with the UCF group, but not sharing code in either direction. The codes performed centering, aperture photometry, and light-curve modeling.

We considered Gaussian, center-of-light, and least-asymmetry (Lust et al., 2014) centering, as well as fixed and (at UCF only) variable-aperture (Lewis et al., 2013) photometry. We selected the optimal centering and photometry method by minimizing the standard deviation of normalized residuals (SDNR) and the binned-$\sigma$ $\chi^{2}$ ($\chi^{2}_{\rm bin}$, Deming et al., 2015) of the model. This second metric looks for a broad-bandwidth solution by comparing a curve of SDNR vs. bin size to the expected inverse square root. We find that $\chi^{2}_{\rm bin}$ more sucessfully selects against correlated noise, so we present the results using that selection criterion. The raw photometry and the position of the target on the detector relative to pixel center are shown in Figure 1.

To remove IRAC’s intrapixel sensitivity variations, we applied both BiLinearly Interpolated Subpixel Sensitivity mapping (BLISS; Stevenson et al., 2012) and Pixel-Level Decorrelation (PLD; Deming et al., 2015), using independent codes we each developed. In brief, BLISS iteratively computes a subpixel-resolution sensitivity grid from the light curve. We account for other effects (transit features, non-flat baselines, etc.) in other model components, all of which fit simultaneously. The flux is then:

[TABLE]

where $F_{s}$ is the stellar flux, $Tr$ is a transit model (e.g., Mandel & Agol, 2002; Rappaport et al., 2014), $M$ is the subpixel sensitivity grid, and $R$ is the non-flat baseline (typically linear or quadratic). PLD corrects the same effect by noting that motion of the target on the detector will be correlated with individual pixel flux values, and models the light curve as a weighted sum of the brightest pixels, after normalisation:

[TABLE]

where $i$ denotes each of $n$ pixels, $c_{i}$ are the weights, and $\hat{P_{i}}$ are normalized pixel values. Both methods have been used extensively to correct Spitzer photometry (e.g., Stevenson et al., 2012; Blecic et al., 2014; Cubillos et al., 2014; Deming et al., 2015; Buhler et al., 2016; Wong et al., 2016).

We fit each model to every combination of centering method, photometry method, and aperture size. We used fixed apertures with 1.5–4.0 pixel radii in 0.25 pixel increments, and variable apertures with radii from $\sqrt{N}$ to $\sqrt{N}$+2.0 pixels in 0.25 pixel increments. $N$ is the ”noise pixel” parameter (Lewis et al., 2013; Spitzer IRAC handbook), defined as

[TABLE]

where $I(i)$ is the intensity of pixel $i$, and all pixels within the centering aperture are considered. We used a 17x17 pixel box, centered on the pixel containing Proxima Centauri, for centering. We take the combination of centering and photometry that result in the lowest $\chi^{2}_{\rm bin}$, Gaussian centering with a fixed 2.0 pixel radius aperture, as the best.

POET finds the best-fitting model using least squares. Since we find that the Spitzer pipeline tends to overestimate uncertainties, we rescale our photometric uncertainties such that the reduced $\chi^{2}$ of the best fit is 1. For fits with a BLISS map, we set the $x$ and $y$ widths of the subpixel grid equal to the root-mean-square of the point-to-point variation in the $x$ and $y$ positions found from centering. We also require that each subpixel bin contain at least 4 frames.

We then explore the parameter space using Multi-Core Markov-Chain Monte Carlo (MC3, Cubillos et al., 2017), a Markov-Chain Monte Carlo (MCMC) wrapper, to determine accurate parameter uncertainties. Our Markov chains use DEMCzs, or “snooker”, a form of differential evolution Markov Chain, to explore the parameter space efficiently (ter Braak, 2006; ter Braak & Vrugt, 2008). We run sufficent iterations for all parameters to pass the Gelman & Rubin convergence test within 1% of unity (Gelman & Rubin, 1992).

With this analysis, only a single asymmetric transit-like feature appears, towards the end of the time series (see Fig. 2 upper panel). At $\sim$0.3% max depth below the continuum, it is smaller than the 0.5% transit depth that we predict for Proxima b using the parameters determined from the RV modeling effort. However, if we consider variable-aperture photometry (which is not preferred by our noise-minimization metrics, as this method significantly increases the white noise in the light curve), the feature disappears completely (see Fig. 2 lower panel). The asymmetric transit-like feature corresponds to a telescope vibration that smears the point-spread function, frequently associated with the ”noise pixel” parameter. The strength of the vibrational effect depends upon the noise-minimization metric (i.e., the choice of centering and photometry techniques) as well as the decorrelation model (BLISS, PLD, etc.). Further study of this effect, including detection and mitigation techniques, will appear in a forthcoming paper (Challener et al. 2019, in prep).

Proxima flares over 60 times per day (Davenport et al., 2016), giving an optical light curve stability at the 0.5–1% level. We find an SDNR of 7527 and 9300 ppm for the fixed and variable aperture cases, respectively, in this infrared filter. When binned over a typical $\sim 2$ hour transit, these SDNR drop to 170 and 222 ppm. Taking these as uncertainties and the stellar radius to be 0.154 $R_{\odot}$, we rule out transiting objects with radii $>$0.43 $R_{\rm{\oplus}}$ at the 3$\sigma$ level of confidence, using the more-conservative variable-aperture photometry. Previously detected features in optical light curves for this star (e.g., Kipping et al., 2017; Li et al., 2017; Liu et al., 2018) are not due to Proxima b. They may be residual correlated noise from the stellar activity.

3 Spectroscopic Observations and Analysis

3.0.1 The Ultraviolet and Visual Echelle Spectrograph (UVES)

The RV data from the Ultraviolet and Visual Echelle Spectrograph (UVES) and High Accuracy Radial velocity Planet Searcher (HARPS; see Anglada-Escudé et al., 2016 for details) that were used to discover Proxima b, along with new HARPS data observed as part of the Red Dots111https://reddots.space/ campaign, allowed further confirmation of the existence of Proxima b, along with improved upper limits on the mass of any additional body orbiting Proxima with an orbit interior to that of Proxima b. The 77 observations from UVES span a baseline of over seven years, with RVs acquired from Julian Date 2,451,634.731 to 2,454,189.714 at signal-to-noise ratios over 100 at 5500Å, necessary to measure optical RVs at the 1 ms*-1* level, and generally taken at low observing cadence (for instance, no measurements were acquired on successive nights). The resolving power of the UVES spectra was $R=\lambda/\Delta\lambda$ = 100,000 – 120,000, where $\lambda$ is wavelength, through application of image slicer #3, which redistributes the light from the $1^{\prime\prime}$ opening along the 0.3*′′* slit. An iodine cell placed in the optical light path before entrance into the echelle spectrograph had an operational temperature of 70∘ C. It imprinted a dense forest of molecular iodine lines on the stellar spectra between $\sim$5000$\--$6200Å. More details on the observational strategy for the UVES spectra can be found in Kürster et al. (2003), Endl & Kürster (2008), Zechmeister et al. (2009), and Anglada-Escudé et al. (2016).

The treatment of the spectra follow the classical reduction steps for such data (see Baranne et al., 1996; Jenkins et al., 2017), including debiasing, cosmic ray removal, echelle order location, flatfielding, scattered-light removal, spectral extraction, spectral deblazing, and wavelength calibration. For UVES, the resulting spectrum, when compared against a previously measured Fourier transform spectrometer (FTS) iodine spectrum at much higher resolution, enables modeling of the iodine spectrum. However, this requires a prior template observation of the star without the iodine cell and without the spectrograph’s point-spread function included. The mathematical form of the process is

[TABLE]

where $F_{o}$ is the observed spectrum, $F_{t}$ is the observed stellar template without iodine, $F_{i}$ is the FTS iodine spectrum, $\Delta\lambda$ is the subsequent shift in wavelength between the iodine model and the FTS observation, and $c$ is a normalisation constant. Forward modeling of the iodine spectrum, the stellar spectrum (Proxima), and the instrumental response can yield precise measurements of the star’s RV signature. The final UVES internal precision using this technique is $\sim$0.6 m/s.

3.0.2 High Accuracy Radial velocity Planet Searcher (HARPS)

The HARPS data cover Julian Dates from 2,457,406.870 to 2,458,027.479, giving rise to 115 measurements after cleaning outliers. Contrary to the UVES observing cadence, the HARPS data were acquired at very high observing frequency, with multiple observations being taken on individual nights, and many covering successive nights.

HARPS uses a different observing methodology, but there are many similarities in the processing. For initial reduction we used the HARPS Data Reduction Software (Mayor et al., 2003b) to perform the steps outlined above for UVES, extracting calibrated spectra from the raw images. To calculate the RVs, we used the HARPS-TERRA package (Anglada-Escudé & Butler, 2012). Without a gas cell in front of the spectrograph, the calculation can be performed directly by comparing the observed spectra. TERRA first shifts and combines all high S/N observed spectra. Their mean shift is the reference velocity for the rest of the observations. Precise velocities come from a Gaussian fit to the minimum of the cross-correlation function between the template and individual observations. We correct this for internal wavelength drift, measured by a ThAr gas lamp spectrum taken at the same time as the observations. Internal velocity precision for Proxima is $\sim$0.9 m/s.

3.1 Radial Velocity Constraints

In order to analyse the latest RV data for this work, which had the aim of confirming the existence of Proxima b and searching for additional companions on orbits interior to that of the planet, we employed the Exoplanet Mcmc Parallel tEmpering Radial VelOcity fitteR (EMPEROR) code (Peña Rojas & Jenkins, 2018). The algorithm uses MCMC to explore the posterior parameter space, along with Bayesian statistics to determine if any signal exists. In this work we employed EMPEROR with a first-order moving-average (MA), correlated-noise model, to smooth out the high-frequency noise that tends to dominate RV measurements. No linear correlation terms were included, therefore we did not model any impact from stellar activity that is tracked by measured indices drawn from the stellar spectra themselves (for example, see Díaz et al., 2018), beyond the MA model, since when included they were mostly found to be statistically similar to zero. The model ($m(t)$) we employ as function of time for a given planet ($k$) and dataset ($d$) is described by

[TABLE]

where $K$ is the semiamplitude of the planet model, $e$ is the eccentricity, $\omega$ is longitude of periastron, $T$ corresponds to the time of periastron passage, $\dot{\gamma}$ is the linear trend added to the model, $\sigma_{jit}$ is the excess jitter noise, and $MA$ is the moving average model.

The full timeseries was modeled as four separate datasets simultaneously, those data coming from four separate programs, three using HARPS and one using UVES (see Table 1 for a brief summary). The pre-2016 HARPS data, HARPS Pale Red Dot (PRD), and UVES data were discussed in (Anglada-Escudé et al., 2016) and the HARPS Red Dots (RD) data, which was an extension and expansion of the PRD program and followed the observing procedure set out there, can be found at the website222https://reddots.space/. The simultaneous modeling contained five Keplerian parameters to model the planet, along with independent data offsets, MA coefficients, and excess noise (jitter) parameters, and finally a linear trend (Equation 5). The MA modeling finds all four datasets have correlation coefficients that are statistically significantly different from zero at around the 3$\sigma$ level of confidence (see Table 2). Therefore, both the high cadence and low cadence data require a correlated noise model to extract the most RV information from the negative effects of the noise. It also appears that the jitter level slightly decreases between the PRD timeseries and the RD data, possibly due to a decrease in the activity state of Proxima. The pre-2016 and UVES jitter values are significantly larger, likely due to the lower precision of these datasets compared to the post-HARPS upgrade observations. No signficant linear trend was found in the full timeseries, however this does not rule out longer period companions as we split the data up into individual runs and we also nightly binned any data that were observed on the same night. Such data handling would ultimately disfavour long period signals in the data, which is fine for these purposes since we were looking to constrain any planets with orbital periods less than that of Proxima b.

EMPEROR provides excellent constraints on the orbital characteristics of Proxima b, particularly refining the orbital period of 11.1855 days to a level better than $\pm$2.3 minutes (see Table 2). The eccentricity is also better constrained, with a 3$\sigma$ upper limit of 0.29, a movement towards zero of 0.06 compared with the value published in (Anglada-Escudé et al., 2016). Further limits on the eccentricity are warranted, since lower values of eccentricity require the planet’s orbit to be tidally locked to the star, providing additional constraints on the habitability of Proxima b (see for example, Ribas et al., 2016). These results highlight that the latest HARPS data are in excellent agreement with the previous data (see Fig. 3) and strongly confirm the existence of Proxima b. We find that at orbital periods shorter than that of Proxima b, the radial-velocity precision we can reach is $\sim$0.5 m/s, mainly coming from the high-cadence datasets, placing an upper limit of $\sim$0.5 $M_{\rm{\oplus}}$  on any possible inner Proxima c. The semiamplitude of the signal we find here is also in excellent agreement with that already found for Proxima b, with a difference of only 0.06 ms*-1*  between the model published in Anglada-Escudé et al. (2016) and the model found here; although the uncertainties we find are almost half those found previously.

4 Summary

We have addressed the recent claims of transit-like events in optical photometry arising from the habitable-zone terrestrial planet Proxima b. We observed the system with the Spitzer Space Telescope for $\sim$48 hours at 4.5 m. The observations covered the 99% probability window predicted for the transit using the published RV model in Anglada-Escudé et al. (2016). The limits on this window were drawn from the posterior density distribution of the model, assuming 99% uncertainty limits on the model parameters like period and eccentricity.

Our observations and BLISS analysis allowed us to reach an unbinned photometric precision of 7500 ppm, with a 2 hr (rough transit duration) binned precision of 200 ppm. No transit-like event could be attributed to the passage of Proxima b in front of its star. The previously witnessed transit-like events may result from residual correlated noise arising from the star’s complex and frequent flaring and activity patterns. Our photometric precision places a 3$\sigma$ upper limit on the size of a transiting Proxima b of 0.4 $R_{\rm{\oplus}}$. This corresponds to an implausible minimum density of $\sim$112 g cm*-3*.

We performed a short radial-velocity experiment to search for additional small planets interior to the orbit of Proxima b, whilst constraining better Proxima b’s orbital characteristics. Beyond the data published in Anglada-Escudé et al. (2016), we also included newly observed HARPS data from the Red Dots program. After fitting for the orbit of Proxima b, the residuals reveal no inner planet down to the 0.5 m/s level, which relates to planets with minimum masses of 0.5 $M_{\rm{\oplus}}$. The orbital period and eccentricity of Proxima b’s orbital solution were also better constrained in this process, with a precision in period of better than $\pm$30 s found, and a 3$\sigma$ upper limit on the eccentricity of 0.29.

Finally, we did witness a transit-like event at the 0.3% depth level and with an asymmetric morphology. However, we found we could remove the feature completely from the time series by using variable-radius photometry apertures. A study of this and similar features in additional Spitzer data of Proxima and beyond, as well as detection and treatment methods, will be published in a future paper (Challener et al. 2019, in prep).

Acknowledgments

We thank the Spitzer Science Center staff for making these observations possible. This work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. We also thank the anonymous referee for their efficient and detailed review. The authors acknowledge support from the following: CATA-Basal/Chile PB06 Conicyt and Fondecyt/Chile project #1161218 (JSJ). CONICYT Chile through CONICYT-PFCHA/Doctorado Nacional/2017-21171752 (JP). Spanish MINECO programs AYA2016-79245-C03-03-P, ESP2017-87676-C05-02-R (ER), ESP2016-80435-C2-2-R (EP) and through the ”Centre of Excellence Severo Ochoa” award SEV-2017-0709 (PJA, CRL and ER). STFC Consolidated Grant ST/P000592/1 (GAE). NASA Planetary Atmospheres Program grant NNX12AI69G, NASA Astrophysics Data Analysis Program grant NNX13AF38G. Spanish Ministry of Science, Innovation and Universities and the Fondo Europeo de Desarrollo Regional (FEDER) through grant ESP2016-80435-C2-1-R (IR). We thank contributors to SciPy, Matplotlib, and the Python Programming Language; the free and open-source community; and the NASA Astrophysics Data System for software and services. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programmes 096.C-0082, 191.C-0505, and 099.C-0880.