TESS and HARPS reveal two sub-Neptunes around TOI 1062

The Transiting Exoplanet Survey Satellite (\textit{TESS}) mission was designed to perform an all-sky search of planets around bright and nearby stars. Here we report the discovery of two sub-Neptunes orbiting around the TOI 1062 (TIC 299799658), a V=10.25 G9V star observed in the TESS Sectors 1, 13, 27&28. We use precise radial velocity observations from HARPS to confirm and characterize these two planets. TOI 1062b has a radius of 2.265^{+0.095}_{-0.091} Re, a mass of 11.8 +\- 1.4 Me, and an orbital period of 4.115050 +/- 0.000007 days. The second planet is not transiting, has a minimum mass of 7.4 +/- 1.6 Me and is near the 2:1 mean motion resonance with the innermost planet with an orbital period of 8.13^{+0.02}_{-0.01} days. We performed a dynamical analysis to explore the proximity of the system to this resonance, and to attempt at further constraining the orbital parameters. The transiting planet has a mean density of 5.58^{+1.00}_{-0.89} g cm^-3 and an analysis of its internal structure reveals that it is expected to have a small volatile envelope accounting for 0.35% of the mass at maximum. The star's brightness and the proximity of the inner planet to the"radius gap"make it an interesting candidate for transmission spectroscopy, which could further constrain the composition and internal structure of TOI 1062b.


Introduction
The Kepler mission was the first exoplanet mission to perform a large statistical survey of transiting exoplanets (Borucki et al. 2010;Howell et al. 2014), and it impacted the field significantly with the detection of over 2300 exoplanets (see NASA Exoplanet Archive Akeson et al. 2013). Nevertheless, due to the faintness of the stars targeted by Kepler, only a small fraction are suitable for radial velocity (RV) follow-up. The Transiting Exoplanet Survey Satellite (TESS) has been designed to survey 85% of the sky for transiting exoplanets around bright, nearby stars (Ricker et al. 2015). It spent the first year of its mission searching for planets in the Ecliptic Southern Hemisphere. More than 2000 candidates have been detected and over 80 have been confirmed, expanding the number of small planets around cool stars (Guerrero 2020). RV follow-up programs have allowed for a rapid validation of transiting planet candidates and to constrain their masses. Accurate measurements of the planetary mass and radius provide information on the planet's bulk density. However, a full characterization on the planetary composition and internal structure is extremely challenging due to the degenerate nature of this problem, where various compositions and interiors can lead to the same average density. Nevertheless, by studying the population of planets in this mass/size range, it is possible to better understand the dominating processes of planet formation and evolution in a statistical manner (Alibert et al. 2010(Alibert et al. , 2015. A key signature in the planet population is the "hot Neptunian desert", which corresponds to a lack of Neptune-mass planets close to their host stars (Szabó & Kiss 2011). The desert is likely to arise by a combination of photoevaporation and tidal disruption (Beaugé & Nesvorný 2013;Mazeh et al. 2016;Owen & Lai 2018). At small orbital distances, planets are subject to intense radiation from the their host stars, and Neptune-like planets are probably not massive enough to retain their gaseous envelopes. An important advantage of the all-sky nature of TESS resides in its ability to find the brightest examples of stars hosting rare sub-populations of planets, such as the exoplanets lying inside the desert. Indeed, TESS has been able to find the first planets deep inside the desert, for example TOI-849b ), a recently discovered remnant core of a giant planet, and LTT9779b (Jenkins et al. 2020), an ultra-hot Neptune. The NCORES HARPS large program is designed to further study planets which may have undergone substantial envelope loss. Some of the recently discovered planets under this program include the mentioned TOI-849b and the three mini Neptunes TOI-125b, c, d (Nielsen et al. 2020).
Another relevant feature in the exoplanet population is a lack of planets with radii between 1.5R ⊕ and 2R ⊕ (Fulton et al. 2017) known as the "radius valley". This suggests that there is a transition between the super-Earth and sub-Neptune populations (Owen & Jackson 2012). This gap could be a result of stellar irradiation (Lopez & Fortney 2013) or core-powered mass-loss (e.g., Ginzburg et al. 2018) that leads to atmospheric evaporation that strips the planet down to its core or, alternatively, consequence of a gas-poor formation (e.g., Gupta & Schlichting 2019). The observed features of the exoplanet populations are not well understood and more exoplanet discoveries are crucial to reveal the overall picture of planets in the super-Earth and Neptunesize/mass regimes.
In this paper we present the discovery and confirmation of two highly irradiated sub-Neptunes hosted by the TESS Object of Interest (TOI) 1062, a nearby (d= 82 pc) and bright (V mag  The TESS full light curve with the 2-minute cadence data is shown in light-grey, and the same data binned to 10 min in dark-grey. The gaps in the data coverage are due to observation interruptions of the TESS spacecraft to send the data to Earth, that occurs after each TESS orbit of 13.7 days. The dots of the LCOGT are bigger for better visualization. The lower panels shows the phase folded TESS light curves for TOI 1062b (left) and TOI 1062c (right) with the 2 min cadence data in light-grey and binned to 10 minutes in black. We find that TOI 1062c does not transit.
10.2) G9V star (see Table 1 for a full summary of the stellar properties). We also used intensive radial velocity follow-up observations with HARPS to confirm the planetary nature of the transit detected in the TESS data and precisely determine the properties of the planetary system. The paper is organized as follows: in Section 2 we present the data collected on the system to discover and validate the planets, in Section 3 we analyze the host star fundamental parameters, and we characterize and discuss the planets in Section 4. Our conclusions are presented in Section 5. We used the publicly available Presearch Data Conditioning (PDC-SAP) flux time series (e.g. Twicken et al. 2010;Smith et al. 2012) provided by SPOC for the transit modelling, which is corrected for common trends and artifacts, for crowding and for the finite flux fraction in the photometric aperture. Figure 1 shows the 2 min cadence TESS light curve, along with the phase folded curve for TOI 1062b.

Ground based photometry with LCOGT
We obtained ground-based time-series follow-up photometry of full transits of TOI 1062.01 on 24th August 2019 in i -band and on 22nd September 2020 in Pan-STARSS z-short band from the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network node at South Africa Astronomical Observatory. We used the TESS Transit Finder which is a customized version of the Tapir software package (Jensen 2013), to schedule the observations. The 4096x4096 LCO SINISTRO cameras have an image scale of 0.389" per pixel, providing a field of view of 26'x26'. The standard LCOGT BANZAI pipeline was used to calibrate the images, and the photometric data were extracted with the AstroImageJ (AIJ) software package ).
The initial i -band observation intentionally saturated the target star to check for a faint nearby eclipsing binary (NEB) that could be contaminating the TESS photometric aperture. To account for possible contamination from the wings of neighboring star PSFs, we searched for NEBs out to 2.5' from the target star. If fully blended in the SPOC aperture, a neighboring star that is fainter than the target star by 8.8 magnitudes in TESS-band could produce the SPOC-reported flux deficit at mid-transit (assuming a 100% eclipse). To account for possible delta-magnitude differences between TESS-band and i -band, we included an extra 0.5 magnitudes fainter (down to TESS-band magnitude 18.3). Our search ruled out NEBs in all 24 neighboring stars that meet our search criteria. The z-short band observation was defocused and exposed appropriately to measure the TOI 1062.01 light curve at high precision. For the data reduction process we used a 10.5"radius aperture and selected four comparison stars with brightness comparable to that of the target star. The inner and outer radii of the background annulus were then set to the corresponding AIJ-calculated values of 18.7" and 28", respectively. This combination was chosen as it minimized the amount of noise in the light curve. A 0.5 ppt transit, with model residuals of 0.3 ppt in 10-minute bins, was detected in the uncontaminated target aperture.

HARPS follow-up
We collected 90 high-resolution spectra of TOI 1062 using the High Accuracy Radial Velocity Searcher (HARPS) spectrograph mounted at the ESO 3.6m telescope of La Silla Observatory, Chile (Mayor et al. 2003a), with the goal of precisely determining the mass of the planet candidate and to search for additional planets in the system. The observations were carried out as part of the NCORES large programme (ID 1102.C-0249, PI: Armstrong) between 31st August 2019 and 17th January 2020. HARPS is a stabilised high-resolution (R 115000) echelle spectrograph which can reach sub-m s −1 RV precision (Mayor et al. 2003b). The instrument was used in high-accuracy mode (HAM), with a 1" fiber on the star and another one to monitor the sky-background. We used exposure times of 30 minutes.
The standard HARPS Data Reduction Software (DRS) was used to reduce the data, using a K0 mask for both the crosscorrelation function (CCF) (Pepe et al. 2002;Baranne et al. 1996) and the colour correction, and reaching a typical signal-to-noise ratio per pixel of 60 and a photon-noise uncertainty of 1.4 m s −1 . For each spectrum the usual activity indicators (S-index, Hα-index, Na-index, Ca-index, log R HK ), the full width half maximum (FWHM), the line bisector and the contrast of the CCF were measured.

High resolution imaging
To check for the presence of contaminating stars in the TESS photometric aperture, on 14th January 2020 TOI 1062 was observed using the Zorro speckle imager (Scott 2019), mounted on the 8.1m Gemini South telescope in Cerro Pachon, Chile. Zorro uses high speed electron-multiplying CCDs (EMCCDs) to simultaneously acquire data in two bands centered at 562 nm and 832 nm. The data is reduced following the procedures described in Howell et al. (2011), and provides output data including a reconstructed image and robust limits on companion detections. The contrast achieved in the resulting reconstructed image is shown in Fig-J  ure 2. We see that the Zorro speckle images reach a contrast of ∆mag=8.25 at a separation of 1.1" in the 832 nm band, and cover a spatial range of 9.9 to 98.6 au around the star with contrasts between 5 and 8 mag in both bands, showing no close companions to TOI 1062b.

Analysis of the spectrum
The stellar atmospheric parameters (T eff , log g and [Fe/H]) and respective error bars were derived using the methodology described in Sousa (2014); Santos et al. (2013). Briefly, we make use of the equivalent widths (EW) of iron lines, as measured in the combined HARPS spectrum of TOI 1062 using the ARES v2 code 1 , and we assume ionization and excitation equilibrium. This method makes use of a grid of Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973). This approach leads to the following stellar parameters: T eff = 5328 ± 56K, [Fe/H]=0.14 ± 0.04 dex and log g= 4.55 ± 0.09. This is in agreement with the log g derived making use of the Gaia parallax and luminosity (Santos et al. 2004)

Analysis of the spectral energy distribution
As an independent determination of the stellar parameters, we also performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia DR2 parallax (adjusted by +0.08 mas to account for the systematic offset reported by Stassun & Torres 2018), in order to determine 1 The most recent version of ARES code (ARES v2) can be downloaded at http://www.astro.up.pt/∼sousasag/ares an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016); Stassun et al. (2017. We pulled the B T V T magnitudes from Tycho-2, the BVi magnitudes from APASS, the JHK S magnitudes from 2MASS, the W1-W4 magnitudes from WISE, the G, G BP , G RP magnitudes from Gaia, and the NUV magnitude from GALEX. Together, the available photometry spans the full stellar SED over the wavelength range 0.2-22 µm (see Figure 3).
We performed a fit using Kurucz stellar atmosphere models, with the effective temperature (T eff ), metallicity ([Fe/H]), and surface gravity (log g) adopted from the spectroscopic analysis. The only additional free parameter is the extinction (A V ), which we fixed to be zero due to the star's proximity. The resulting fit is very good ( Figure 3) with a reduced χ 2 of 1.8. The reduced χ 2 is improved to 1.2 by excluding the GALEX NUV flux, which exhibits a modest excess; using the empirical relations of Findeisen et al. (2011), the GALEX NUV excess implies an activity level log R HK = −4.9 ± 0.1, consistent with the spectroscopically determined value.
Integrating the (unreddened) model SED gives the bolometric flux at Earth, F bol = 2.429 ± 0.057 × 10 −9 erg s −1 cm −2 . Taking the F bol and T eff together with the Gaia DR2 parallax, gives the stellar radius, R = 0.838 ± 0.020 R . In addition, we can use the R together with the spectroscopic log g to obtain an empirical mass estimate of M = 0.91 ± 0.19M , which is consistent with that obtained via empirical relations of Torres et al. (2010) and a 6% error from the empirical relation itself, M = 0.97 ± 0.06M . These values are in agreement with the ones derived from the spectral analysis in Section 3.1.

Rotational period and age
The stellar age can be derived from empirical rotation-age relations of Mamajek & Hillenbrand (2008). The HARPS calibration from the CCF-FWHM and the color index B-V gives v sin i= 2.13 ± 0.5 km/s, then together with R , we estimate the stellar rotation period P rot / sin i = 20.5 ± 2.9 d. In addition, Fig 4 shows that the peak of the periodogram of the TESS light-curves of Sectors 27 and 28 is at 21.82.4 days. The measured rotational modulation is nearly identical in the other TESS sectors. Another estimate of the expected rotation can be inferred from the average value of the Ca ii H & K chromospheric activity indicator for TOI 1062, which is log R HK = −4.80 ± 0.05. Using the empirical relations in Mamajek & Hillenbrand (2008) we get a rotation period of 21.9 +3 −2 days, in agreement with the result derived from the HARPS calibration and the modulation of the RAW lightcurves. Another way of constraining the rotational period of the star consists of fitting the RVs with time-dependent gaussian processes (GP). We use the celerite package (e.g., Foreman-Mackey et al. 2017) with a quasi-periodic kernel, and obtain a rotational period of 22.4 +4.6 −6.7 days, which is in agreement with the previous estimations.
Using the rotational period obtained from the TESS data and the relations of Mamajek & Hillenbrand (2008), the age estimate is τ =2.5±0.3 Gyr. Note that, because the v sin i is a projected rotational velocity, it is a lower limit on the true rotational velocity, so P rot / sin i is an upper limit on the true P rot , which implies that the age estimate above provides an upper limit to the stellar age.  We also attempted to derive the stellar age by making use of the so-called chemical clocks (e.g. Delgado Mena et al. 2019), however we obtain a value much higher than with other methods described below. The reason is that the low T eff of TOI 1062 lies at the limit of applicability of such empirical relations and thus we refrained to use the age provided by this method.

Stellar abundances
Stellar abundances of refractory elements were derived using the classical curve-of-growth analysis method assuming local thermodynamic equilibrium (e.g. Adibekyan et al. 2012Adibekyan et al. , 2015Delgado Mena et al. 2017). Abundances of the volatile elements, C and O, were derived following the method of Delgado Mena et al. (2010); Bertran de Lis et al. (2015). Since the two spectral lines of oxygen are usually weak and the 6300.3Å line is blended with Ni and CN lines, the EWs of these lines were manually measured with the task splot in IRAF. All the [X/H] ratios are obtained by doing a differential analysis with respect to a high Table 2: Prior parameter distribution of the global fit with juliet. U(a, b) indicates a uniform distribution between a and b; J(a, b) a Jeffrey or log-uniform distribution between a and b; and N(a, b) a normal distribution with mean a and stardard deviation b.

Parameter
Prior distribution Instrumental Parameters: Planetary parameters: S/N solar (Vesta) spectrum from HARPS. The obtained values for this star are normal considering its metallicity except for the lower than expected value of oxygen. We find a high dispersion between both oxygen indicators that has to be taken with caution due to the difficulty of measuring reliable oxygen abundances for cool metal rich stars. The stellar abundances of the elements are presented in Table 1.

Signal identification
In order to search for periodicities in the RV data, we computed the l 1 periodogram (e.g. Hara et al. 2017Hara et al. , 2020. This tool is designed to find periodicities in time-series data corresponding to exoplanets detected using radial velocity data. It can be used similarly to a typical Lomb-Scargle periodogram or its variants (e.g. Baluev 2008), but it is based on a different principle.
The l 1 periodogram can analyze the radial velocity without the need to estimate the frequency iteratively, using the theory of compressed sensing adapted for handling correlated noise. Instead of fitting a sinusoidal function at each frequency as the Lomb-Scargle-type periodograms, the l 1 periodograms searches directly for a representation of the input signal as a sum of sinusoids. As a result, all the frequencies of the signal are searched simultaneously. This relies on the so-called "sparse recovery" tools, which are designed to find a representation of an input signal as a linear combination of sine functions.
The l 1 periodogram outputs a figure which has a similar aspect as a standard Lomb-Scargle-type periodogram, but with fewer peaks due to aliasing. We use the l 1 periodogram of the HARPS data on a grid of frequencies between 0.01 and 0.3 cycles per day. In order to measure the significance of a possible signal the False Alarm Probability (FAP) is calculated, which encodes the probability of measuring a peak of a given height conditioned on the assumption that the data consists of Gaussian noise with no periodic component. Typically a signal at the 1% FAP level is considered suggestive, and 0.1 % FAP level statistically significant. Figure 5 shows two l 1 periodograms of the HARPS RVs. The top panel corresponds to the undetrended RVs. In this case we detect three significant signals at 4.11d (FAP 1.10 −4 %), 7.98 days (FAP 7.10 −6 %) and 26.6 days (FAP 3.10 −5 %). Although the coefficient amplitude of the second and third peak are relatively close, their FAP can differ significantly since it is estimated in a sequential manner, i.e. the highest peak is computed against zero, the second peak against the first and the third peak against the second. Figure 6 clearly shows the effects of activity in the FWHM at the end of the time series. It suggests that the star was more active at that moment, probably due to one spot crossing event during one single rotational period, standing for only 20 days (from BJD 58830 to 58850). However, we find that signal at 26.6 days disappears when we detrend the RVs using GP modelling with SPLEAF (e.g. Delisle et al. 2020) (see bottom panel of Figure 5). The GP is trained on the logR HK indicator and then it is used to detrend the RVs using a linear fit. This method gives two significant peaks at at 4.11d with FAP 6.10 −5 % and 8.00 days with FAP 4.10 −4 %. Another way to detrend the stellar activity is to fit the correlations between activity indicators as the full width at half maximum (FWHM), the bisector or the S-index, and subtract them in order to correct for any long-term effects resulting from a magnetic cycle. We do so dividing the data into different time regimes selected in order to maximize the correlation of the RV data with the activity indicators in each bin. This approach also identifies two peaks at 4.11d with FAP 7.10 −8 % and 8.01 days with FAP 3.10 −4 %. Finally, we also find the same two significant signals when we do not use the data TBJD 58730 to TBJD 58750, which corresponds to period in which the star was more active. The signal at 26.6 days only appears when using the full undetrended RV sample which is clearly affected by stellar activity and it is within 1.6-sigma of the rotation period derived from the chromospheric activity indicator. We therefore conclude that only the peaks at ∼4.11 days and 8 days are caused by planets. In addition, despite the 8 day signal is nearly double the ∼4.11 day, in all the cases presented we find the 8 day signal after removing the ∼4.11 day signal. We also evaluate the evidence of the second planet by comparing the difference in Bayesian Information Criterion (∆BIC) between the one-eccentric-planet and the two-eccentric-planets models. We get a value of ∆BIC=11, which suggests a strong evidence.

Global fit with juliet
We determine the planetary parameters using the publicly available software juliet (e.g., Espinoza et al. 2019). juliet allows us to jointly fit the TESS and LCOGT photometry and the HARPS radial velocities with GP. The software is built on several publicly available tools for the modeling of the photometric data (using the batman package, Kreidberg 2015), the radial velocities (using the radvel package, Fulton et al. 2018), and also to incorporate GP (via the celerite package, Foreman-Mackey et al. 2017). The parameter space is explored using nested sampling, with the MultiNest algorithm (e.g., Feroz & Hobson 2008) in its Python implementation, PyMultiNest (e.g., Buchner et al. 2014). In addition, juliet computes the Bayesian evidence using dynesty (e.g., Speagle 2020), a Python package to estimate Bayesian posteriors and evidences using Dynamic Nested Sampling. In short, nested sampling algorithms work as follows. The algorithm samples some number of live points randomly from the prior distribution, and the likelihood is evaluated at each of these points. At each iteration the point with the lowest likelihood is replaced by a new point sampled, keeping the number of live points constant. This process is continued until Bayesian evidence reaches some specified value. The number of live points used has to be large enough to adequately sample the parameter space.
The transit model fits the stellar density ρ together with the planetary and jitter parameters. For the stellar density we use the value obtained in Section 3, and the priors of the orbital parameters of the inner planet are taken from ExoFOP. We use the quadratic limb darkening coefficients (q 1 , q 2 ) introduced by Kipping (2013) for the photometric data, since it was shown to be appropriate for space-based missions (Espinoza & Jordán 2015). In addition, instead of fitting the planet-to-star radius ratio and the impact parameter of the orbit, use the parametrization introduced in Espinoza (2018) and fit the parameters r 1 and r 2 to make sure the full exploration of physically plausible values in the (p,b) space. We parametrize the eccentricity and the argument of periastron with √ e sin ω and √ e cos ω, always ensuring that e ≤ 1. We use the celerite aproximate Matern multiplied by exponential kernel to account for the stellar activity. The fit obtained with this kernel is nearly identical to the one obtained with other kernels as the quasi-periodic or the exp-sine-squared kernel.
We fit both planets with juliet and find no hint of transit for TOI 1062c, as shown in Figure 1. The final median planetary parameters determined by the juliet fit are listed in Table 2. We find that TOI 1062b has an orbital period of 4.114 days, a radius of 2.264 +0.095 −0.091 R ⊕ , a mass of 10.17 +0.83 −0.81 M ⊕ and a density of 4.86 +0.85 −0.74 g cm −3 , slightly below the Earth's mean density. Nevertheless, as discussed in detail in Section 5, its composition and internal structure are unlikely to be similar to that of the Earth. We find an eccentricity of 0.18 +0.08 −0.07 for TOI-1062b and 0.14 +0.09 −0.07 for TOI-1062c. Figure 8 shows the eccentricity against the orbital period for the population of observed exoplanets from the NASA Exoplanet Archive, and we see that TOI-1062b is slightly higher than usual for its orbital period. However, it may not be significant due to the bias towards larger eccentricities for nearly circular planets (Lucy & Sweeney 1971). TOI 1062 c instead is found to have an orbital period of 7.988 ± 0.04 days and a minimum mass of 9.35 ± 1.23 M ⊕ . TOI 1062 c is in near the 2:1 motion resonance with its inner companion.

Internal structure
In order to characterize the internal structure of TOI 1062 b, we model its interior considering a pure-iron core, a silicate mantle, a pure-water layer, and a H-He atmosphere. The equations of state (EOSs) used for the iron core are taken from Hakim et al. (2018), the EOS of the silicate-mantle is calculated with PERPLE_X from Connolly (2009)  . We assume an envelope luminosity of L=10 22.52 erg s −1 (equal to Neptune's luminosity). The thickness of the planetary layers were set by defining their masses and solving the structure equations. To obtain the transit radius, we follow Guillot (2010) and evaluate the location where the chord optical depth τ ch is 2/3. We do not use stellar abundances as an additional constraint since it is not clear whether they are a good proxy for the planetary bulk abundances (Wang et al. 2019;Plotnykov & Valencia 2020) and it has been shown that the stellar abundances are not always useful for constraining the internal composition (Otegi et al. 2020b). Figure 9 shows M-R curves tracing the compositions of Earth-like (with a CMF=0.33) and pure-water subjected to a For reference, we also show exoplanets with accurate and reliable mass and radius determinations (Otegi et al. 2020a, accessible on the Data & Analysis Center for Exoplanet DACE 2 ). TOI 1062 b sits above the Earth-like curve and below the pure-water curve, suggesting that it contains a small amount of volatile materials of less than 1% in mass. Figure 9 also displays the insolation flux relative to Earth against radii for the known exoplanets extracted from the NASA Exoplanet Archive, which shows the separate populations of super-Earths and mini-Neptunes. We see that TOI 1062 b sits in the mini-Neptune regime, close to the radius valley. In Otegi et al. (2020a) we identified two distinct populations of volatile-rich and "rocky" (referring to exoplanets expected to have low amounts of volatiles) separated by the water-line. Using these results, we see that TOI 1062b lies in the "rocky" exoplanet regime defined in Otegi et al. (2020a) even if its radius is above the radius valley, suggesting that it is mostly composed of refractory materials by mass.
We use a generalized Bayesian inference analysis using a Nested Sampling scheme (Buchner et al. 2014) to quantify the degeneracy between various interior parameters and produce posterior probability distributions. Fig. 10 shows ternary diagrams of the inferred composition of TOI 1062 b. The ternary diagram shows the degeneracy associated with the determination of the composition of exoplanets with measured mass and radius. We Article number, page 9 of 14   find a median H-He mass fraction of 0.1%, which corresponds to a lower-bound since enriched H-He atmospheres are more compressed and, therefore, increase the planetary H-He mass fraction. Indeed, formation models suggest that sub-Neptunes are likely formed by envelope enrichment (Venturini & Helled 2017). We also find that TOI 1062 b is expected to have a very significant iron core and silicate mantle, accounting for nearly 40% of the planetary mass and thicknesses of 1R ⊕ and 0.5R ⊕ respectively. The water layer has an estimated relative mass fraction of 15%. Nevertheless, the degeneracy between the core, silicate mantle, and water layer in this M-R regime is particularly high (Otegi et al. 2020b), and it does not allow accurate estimates of the masses of these constituents. Interior models cannot distinguish between water and H-He as the source of low density material, so we also run a 3-layer model without the H 2 O envelope. Table 4 lists the inferred mass fractions of the core, mantle,water-layer and H-He envelope for both the 4-layer model and the water-free one. In this case we find that the planet is 0.35% H-He, 26% iron and 73% rock by mass, setting maximum limits for the atmospheric and rock mass since any water added would decrease these mass fractions.

Dynamical analysis
As a general note of this section, the studies presented here were done in the hypothesis of a co-planar planetary system, i.e. the orbits of planets b and c evolve in the same plane. This is consistent with the observations, since TOI 1062c would not transit if it had the same orbital inclination as TOI 1062b.

Resonance
As can be seen from Table 2, the period ratio of the planets is P c /P b = 1.941. The planet pair thus lies close to the 2:1 mean-motion resonance (MMR). How close is the system to the MMR? We explore the structure of the parameter space of TOI-1062 near that resonance, in order to investigate its dynamical state.
We therefore design a two dimensional section of the parameter space defined by the period ratio on one axis and the eccentricity of the outer planet e c on the other. Our section has a resolution of 101x101, i.e., we explore the dynamics of 10201 initial configurations defined by a unique set of (P c /P b ,e c ). All the other orbital parameters were initially fixed at the best values reported in Table 2. We numerically compute the future evolution of each configuration over 20 kyr. This is performed with the adaptive time-step high order N-body integrator IAS15, which is available from the python package REBOUND (Rein & Liu 2012;Rein & Spiegel 2015). The perturbative effect of general relativity described in Anderson et al. (1975) was included via the library REBOUNDx (Tamayo et al. 2019). From these numerical simulations, the level of chaos of each configuration was evaluated with the NAFF fast chaos indicator (Laskar et al. 1992;Laskar 1993). The result is presented in Fig. 11.
The NAFF computes precisely, for each planet, the average mean motion n = 2π P over the two halves of the integration and compares these two estimations. Due to the secular constancy of the semi-major axis of regular orbits, this difference should be small in non-chaotic orbits. The higher is the drift in the average mean-motion, the more chaotic is the orbit. In this work, we took as the NAFF of the system the maximal value of this drift over the planetary orbits, in logarithmic scale: NAFF = max i log 10 ∆n i n 0,i , where the subscript i refers to the planet b or c, ∆n i is the drift of the average mean motion over the two halves of the integration, and n 0,i is the initial mean motion of planet i. The color code in Fig. 11 depicts the NAFF defined here above: the bluer, the more regular is the system. We also distinguish the strongly chaotic configurations from the ones that did not finish the integration because of either a close-encounter or an escape of a body (white boxes). We finally explain our choice of 20 kyr for the total integration time. Over this time span, the planetary orbits in the TOI-1062 system are expected to cover several secular cycles during which the semi-major axes oscillate. Covering several of these cycles allows to properly average the secular variations and isolate the chaotic diffusion. We verified numerically that several secular cycles are made up over an integration.
In this map, the 2:1 mean motion resonance appears clearly as the orange band in the middle of the plot. For the current system's parameters and the resolution of our chaoticity map, this resonance seems therefore chaotic. It is important to stress that this picture is highly dependent on the system's parameters. For instance in that case, the arguments of periastron ω b and ω c are close to the alignement. The opposite configuration where the orbits are anti-aligned would show a drastically different picture, with a different strength and apparent stability of the 2:1 MMR. The two vertical lines depict the 1σ window of P c . With the currently estimated orbital parameters and planet masses, the TOI 1062 system certainly lies outside of the 2:1 MMR. Despite the influence that a revision of the parameters may have on the resonance, it seems very unlikely that this conclusion changes.
Finally, we stress that the uncertainty on the eccentricity of the inner planet e b is somewhat large. Modifying this parameter will directly impact the strength of the 2:1 MMR, as the resonance width is expected to increase with the eccentricity. The same note applies for the planetary masses.

Stability constraints
While the chaoticity map presented above (Fig. 11) is informative about the maximal eccentricity of the outer planet allowed  Table 2). The eccentricity of the outer planet e c and the period ratio P c /P b are explored on a grid of 101x101 different configurations. A color is assigned to each box after a short numerical integration, according to its level of chaos using the NAFF indicator (see text). by orbital stability, this estimate should be considered with caution. Indeed, this map explores the parameter space on two dimensions, without taking into account potential correlations with the other parameters. As already discussed, if the best fit estimate of the other parameters changes, the entire map might change as well. In order to give proper stability constraints on the orbital parameters and masses of the planets, we also perform a global stability analysis from the posterior distribution of the joint photometry + RV analysis. The only fixed parameter is the inclination of the outer planet, taken equal to the inclination of the inner planet (i.e. we explored the co-planar case, which is compatible with the outer planet non transiting). We selected a sample of 10000 solutions from that posterior. Each of these is numerically integrated over 20 kyr using the same integrator set-up described in the previous section. The level of chaos is again computed with the NAFF indicator. Over the 10000 configurations, 8455 did survive the entire 20 kyr simulation. The others were thus unstable (escape or close encounter of two bodies). Imposing a more strict stability threshold by further removing all solutions with NAFF > −5 leaves us with a stable posterior of 7224 configurations.
In any case, no stringent constraint can be added on the orbital parameters and masses of the planets. With the sorting in NAFF as defined above, only a slight cut in the eccentricities is observed. The new median estimates and 1σ confidence intervals are e b = 0.162 [0.098,0.228] and e c = 0. 129 [0.065,0.205]. In particular, the global stability analysis does not allow us to further constrain the planetary masses (again, assuming co-planar orbits).

Atmospheric characterization
TOI 1062b poses an interesting target for atmospheric characterisation given its equilibrium temperature is ∼1000 K. We expect a high metallicity atmosphere for a planet with such a small radius and mass, as a result of the strong stellar irradiation received which would result in a significant loss of the H 2 /He envelope. However, such a target is an excellent laboratory to study carbon chemistry at high temperature given its high expected metallicity. Beyond metallicities of 300× solar, the primary carbon bearing species in the atmosphere transitions from being CH 4 /CO to CO 2 in thermochemical equilibrium (e.g. Moses et al. 2013). This planet is thus an ideal target to explore the CO, CH 4 , CO 2 chemical stability boundary which occurs at ∼800 K in photospheric conditions at such high metallicities. The mean molecular weight for such an atmosphere at ∼ 300× metallicity is expected to be ∼4 g/mol, and hence significant features are still present in transmission spectra compared to a solar metallicity atmosphere (at ∼ 2.35 g/mol). Using the transmission spectroscopy metric (TSM) defined in Kempton et al. (2018), we determine the that TOI 1062b has a TSM value of 31.8 assuming a mean molecular mass of 2.3 g/mol and 18.3 assuming 4 g/mol. The lower TSM value for the higher mean molecular weight is due to the reduced atmospheric scale height of a heavier atmosphere. This value is comparable to targets such as Trappist-1f and several simulated TESS targets from Sullivan et al. (2015).
On the other hand, TOI 1062b may have completely lost all of its H 2 /He envelope which would result in an ultra-high metallicity secondary atmosphere. At metallicities 3000× solar, the mean molecular weight of the heavier species such as CO 2 and H 2 O now completely dominate due to the lack of significant H 2 or He. Therefore, the mean molecular weight is likely to be > 18 g/mol, and thus the scale height of spectral features in the transmission spectra are significantly reduced by a factor of 7× over a solar metallicity atmosphere. This reduces the TSM value for such an atmosphere to 5, making constraints on the abundances difficult. However, observations of the terminator may still be able to place a lower limit on the molecular weight and thus the metallicity. These would also provide a direct contrast to planets such as 55 Cancri e which also indicate a high mean molecular weight atmosphere (Jindal et al. 2020). In addition, secondary eclipse spectroscopy for such a cool target will be challenging, but has been achieved for GJ436b (Stevenson et al. 2010). However, for the case of TOI 1062b, the host star is a G-type and thus will result in a very weak planet/star flux contrast in emission spectroscopy.

Conclusions
We present the discovery of two new planets from the TESS mission in the TOI 1062 system. The analysis is based on 2-min cadence TESS observations from 4 sectors, ground-based followup from LCOGT and RV data from the HARPS spectrograph. High-resolution imaging from Zorro speckle imager rules out the presence of nearby companions and potential nearby eclipsing binaries. We find that the host star is rotating with a period of 21.9 days, as well as the existence of a stellar spot whose life-time is similar to the stellar rotational period.
TOI 1062b is expected to be a mini-Neptune with a period of 4.11 days, radius of 2.26±0.1 R ⊕ , and mass of 10.2±0.8 M ⊕ . Internal structure models indicate that TOI 1062b is expected to be composed of significant iron core and silicate mantle, accounting for nearly 40% of the planetary mass each, and that it is expected to have a small volatile envelope of 0.35% of the mass at maximum. TOI 1062c, which is not transiting, is found in the HARPS RV data. We find that its minimum mass is inferred to be 9.35±1.23 M ⊕ , and that it is close to a 2:1 motion resonance with its inner companion, with a period of 8.12 days. The position of TOI 1062b with respect to the radius valley and its high equilibrium temperature of 1000K make it a interesting candidate for atmospheric characterization. The strong stellar irradiation may result in a significant loss of the H-He envelope leaving a high metallicity atmosphere that can be an excellent laboratory to study carbon chemistry at high temperature.