HD 23472: A multi-planetary system with three super-Earths and two potential super-Mercuries

We report the characterisation of a multi-planetary system composed of five exoplanets orbiting the K-dwarf HD~23472 (TOI-174). In addition to the two super-Earths that were previously confirmed, we confirm and characterise three Earth-size planets in the system using ESPRESSO radial velocity observations. The planets of this compact system have periods of $P_d \sim 3.98\,$, $P_e \sim 7.90\,$, $P_f \sim 12.16\,$, $P_b \sim 17.67,\,$ and $P_c \sim 29.80\,$days and radii of $R_d \sim 0.75\,$ , $R_e \sim 0.82,$, $R_f \sim 1.13\,$, $R_b \sim 2.01,\,$ and, $R_c \sim 1.85\,$ \REarth. Because of its small size, its proximity to planet d's transit, and close resonance with planet d, planet e was only recently found. The planetary masses were estimated to be $M_d =0.54\pm0.22$, $M_e =0.76\pm0.30$, $M_f =0.64_{-0.39}^{+0.46}$, $M_b = 8.42_{-0.84}^{+0.83}$, and $M_c = 3.37_{-0.87}^{+0.92}$ \MEarth. These planets are among the lightest planets, with masses measured using the radial velocity method, demonstrating the very high precision of the ESPRESSO spectrograph. We estimated the composition of the system's five planets and found that their gas and water mass fractions increase with stellar distance, suggesting that the system was shaped by irradiation. The high density of the two inner planets ($\rho_d = 7.5_{-3.1}^{+3.9}$ and $\rho_e = 7.5_{-3.0}^{+3.9}\, \mathrm{g.cm^{-3}}$) indicates that they are likely to be super-Mercuries. This is supported by the modelling of the internal structures of the planets, which also suggests that the three outermost planets have significant water or gas content If the existence of two super-Mercuries in the system is confirmed, this system will be the only one known to feature two super-Mercuries, making it an excellent testing bed for theories of super-Mercuries formation. (Abridged)


Introduction
New high-resolution spectrographs, such as the Echelle Spectrograph for Rocky Exoplanet-and Stable Spectroscopic Observations (ESPRESSO) (Pepe et al. 2021), in conjunction with photometric space missions focusing on bright stars, such as K2 (Borucki et al. 2010;Howell et al. 2014), transiting exoplanet survey satellite (TESS) (Ricker et al. 2015), and CHaracterising ExOPlanet Satellite (CHEOPS) (Benz et al. 2021), are pushing the limits of planet characterisation to planets of similar size and mass to Earth. This provides an unprecedented view of the interiors of small exoplanets (e.g. Azevedo Silva et al. 2022;Demangeon et al. 2021;Toledo-Padrón et al. 2020;Lillo-Box et al. 2020). These studies are also pushing the boundaries of small planet characterisations to include cooler planets, thus allowing us to gain a better understanding of how high stellar irradiation shapes a planet's composition.
Multi-planetary systems are especially valuable because they share the same host star and were formed by the same accretion disc (e.g. Ormel et al. 2017;Grimm et al. 2018). The properties of multi-planetary systems have been used to constrain planet formation and evolution models. For example, it has been demonstrated that the sizes of planets in a multi-planetary system are correlated and that there is regular spacing between planets (Lissauer et al. 2011;Ciardi et al. 2013;Weiss et al. 2018), which is referred to as the "peas in a pod" theory. According to a recent study, multi-planetary systems appear to be less similar in mass than in radius (Otegi et al. 2022). However, a larger sample of multi-planetary systems with well-characterised mass and radius is required to confirm this result and uncover additional correlations.
High-energy irradiation received by short period planets causes evaporation of H/He-rich envelopes (e.g. Lammer et al. 2003;Yelle 2004;Owen & Wu 2017;Jin & Mordasini 2018). Atmospheric evaporation has been observed for the low-mass planet GJ 436b (Kulow et al. 2014;Ehrenreich et al. 2015;Lavie et al. 2017). As a result of irradiation close-in planets become denser and smaller (e.g. Lopez et al. 2012;Howe & Burrows 2015). The evaporation theory predicted the existence of a gap in the radius distribution of planets (Owen & Wu 2013). Such a gap was later uncovered using the Kepler sample (Fulton et al. 2017). The California-Kepler Survey revealed that the distribution of planet sizes is bimodal and that there is a radius gap at a planetary radius of " 1.6´1.8 R C (Fulton et al. 2017;Van Eylen et al. 2018). Exoplanets with less than 2 R C are thought to be dry naked cores, whereas exoplanets with more than 2 R C are thought to have icy cores and possibly a gaseous atmosphere (Venturini et al. 2020). Alternatively, the radius gap can also be explained by the theory of core-powered evaporation (e.g. Ginzburg et al. 2018;Gupta & Schlichting 2019). Characterising multi-planetary systems with planets above and below the radius gap can shed light on the mechanism of mass loss that is responsible for these effects.
Some of these naked cores are extremely dense and likely contain an excess of iron, similar to Mercury. Since the discovery of the first super-Mercury, K2-229 b, (Santerne et al. 2018) other exoplanets with an excess of iron have been uncovered: K2-38 b (Toledo-Padrón et al. 2020), K2-106 b (Guenther et al. 2017), Kepler-107 c (Bonomo et al. 2019), Kepler-406 b (Marcy et al. 2014), and HD 137496 b (Azevedo Silva et al. 2022). These planets, with the exception of Kepler-107 c, are all the inner planets of their exoplanetary system. They are all part of multi-planetary systems and have high effective temperatures (ą 1200K). Several theories on planet formation and evolution are under discussion to explain their existence, including a giant impact (Benz et al. 2007), mantle evaporation (Cameron 1985), photophoresis (Wurm et al. 2013), and the possibility of a compressed planetary core (Mocquet et al. 2014).
We chose to characterise HD 23472 (TOI-174) because it serves as an ideal object for gaining further insight into atmosphere evaporation. When we began our observations, the system was reported to have two planets above as well as two planets below the radius gap in the EXOFOP-TESS database (ExoFOP 2019). Furthermore, the planets span a range of stellar incident flux between 8 and 117 of that of the Earth. Hence, HD 23472 evidently stands as a golden target for mass characterisation with ESPRESSO. The new planet that we found below the radius gap makes the system even more interesting. Here, we present the characterisation of the five planets in the HD 23472 system, two of which are likely super-Mercuries. We present our new ESPRESSO radial observations of HD 23472 as well as previous RV and photometric observations in Section 2. We describe our data analysis method in Section 3 and we present the derived system parameters in Section 4. In Section 5.1, we present our analysis of the system's dynamics, which allows us to better constrain the eccentricity of the planets, while in Section 6 we present the analysis of the internal structure of the planets. Finally, our conclusions are summarised in Section 7.

ESPRESSO observations
From 20 July 2019 to 9 April 2021, we collected 104 spectra of HD 23472 (K4V, V mag = 9.7) with ESPRESSO (Pepe et al. 2021), each with an exposure time of 900s. The ESPRESSO high resolution echelle spectrograph is mounted on the ESO Paranal Observatory's Very Large Telescope (VLT). ESPRESSO can collect light from any VLT unit telescope or all of them at the same time. The observations were obtained as part of the ESPRESSO Guaranteed Time Observations (programs 1102.C-0744, 1102.C-0958, and 1104.C-0350) which has already al-lowed to constrain the composition of several small planets (Damasso et al. 2020;Toledo-Padrón et al. 2020;Mortier et al. 2020;Sozzetti et al. 2021;Demangeon et al. 2021). ESPRESSO is isolated in order to maintain constant pressure, temperature, and humidity. We used the single UT high resolution mode (HR11, fast-readout) for all observations which has a spectral resolution of R" 140 000 and covers wavelengths from 380 nm to 788nm. We obtained a S/N of 80 per resolution element at 650nm. All measurements were obtained with the Fabry Perrot in Fibre B for simultaneous calibration.This allows the correction of the instrumental drift with a precision better than 10 cm/s (Wildi et al. 2010).
To extract the radial velocities (RVs), we used the version 2.21 of the ESPRESSO pipeline Data-reduction Software (DRS). The RVs are calculated using the DRS by crosscorrelating the spectra with a stellar line mask (Baranne et al. 1996), which in our case was designed for K6 type stars. Then the DRS fits the cross correlating function (CCF) with an inverted Gaussian profile to obtain the centre of the Gaussian (RV measurement), the full width at half maximum (FWHM), and the amplitude (contrast of the CCF). The RVs' uncertainties are computed using the Bouchy et al. (2001) technique. Moreover, the DRS also computes other activity indicators: the BIS , the depth of the H α line (Suárez Mascareño et al. 2015), the depth of the Sodium doublet (NaD, Díaz et al. 2007), and the S-index (Lovis et al. 2011;Noyes et al. 1984). The median uncertainty of the RV measurements is 0.38 m/s, and their peak-to-peak amplitude is 16.33 m/s. The generalised Lomb Scargle Periodogram (GLS, Zechmeister & Kürster 2009) of the RVs, activity indicators (FWHM, BIS, S-index, NaD, H α ), contrast, and the associated window function are shown in Figure 1. We also show the time series of the RV and activity indicators. The two strongest peaks in the RVs correspond to the known planet b and c that have an orbital period of 17 days and 29.9 days, respectively. There is not a strong evidence for these periodicities in the indicators. The third strongest peak in the RVs " 40 days, on the other hand, is most likely due to stellar activity, as it is also present in all of the indicators except the NaD.

Previous RV observations
Following the public release of transiting candidates in the EXOFOP-TESS database (ExoFOP 2019), Trifonov et al. (2019) analysed 14 public HARPS RV measurements of HD 23472 to constrain the mass of the outermost planets (P b " 17.7 days and P c " 29.8 days). The HARPS RVs allowed them to set an upper limit on the RV semi-amplitude of both planets: K b " 5.33`0 .67 4.2 m s´1 and K c " 4.29`0 .26 3.4 m s´1. Assuming a stellar mass of M˚" 0.75˘0.04 M @ they placed an upper limit on the absolute mass of the planets to be m b " 17.9`1 .4 14 M C and m c " 17.18`1 .1 14 M C . Their calculations suggested a possible 5:3 mean motion resonance (MMR) between the two planets that would result in an oscillation of their period ratio (P c {P b ). We analysed the 14 HARPS RVs and 7 publicly available CORALIE RVs alongside our ESPRESSO measurements, but due to their higher uncertainties, low number and sparseness, they do not improve our results. Hence, they were excluded from the final analysis. The HARPS and CORALIE measurements have a median uncertainty of 1.7 m/s and 6.0 m/s, respectively, while the peak-to-peak amplitude of the variation of the measurements is 19.18 m/s and 33 m/s, respectively. S. C. C. Barros et al.: HD 23472: A multi-planetary system with three super-Earths and two potential super-Mercuries  Later, HD 23472 was observed with the Carnegie Planet Finder Spectrograph (PFS) mounted on the Magellan II telescope (Teske et al. 2021). They obtained 64 observations, allowing them to constrain the masses of the two outermost planets. The median of the uncertainty in the RV measurements is 0.63 m/s and the RV show a peak to peak amplitude of 12.10 m/s. Using their first method, which they refer as the 'juliet' method (without a correction for the activity), they found that the RV semi-amplitude for planet b is K b " 3.39˘0.27 m s´1, which corresponds to m b " 11.7˘1.3 M C , and for planet c it is K c " 1.18˘0.30 m s´1, which corresponds to m c " 4.85˘1.28 M C . This is based on an assumed stellar mass of M˚" 0.780˘0.089 M @ . They obtained slightly different results for the RV semi-amplitude, K b " 2.99˘0.39 m s´1and K c " 1.39˘0.48 m s´1using their second method, called 'radvel' method. This second method included the correction of stellar activity using a Gaussian Process (GP). We included the PFS measurements in our final analysis because they help to constrain the planetary RV signatures.

TESS observations
TESS observed HD 23472 (TIC 425997655, TOI-174) in five sectors (1,2,3,4,11) at a 2 minute cadence and four sectors at a 2 minute-and-20-second cadence (29,30,31,34). The noncontinuous observations span is " 900 days. We downloaded light curves computed by the TESS pipeline (Jenkins et al. 2016) from the Mikulski Archive for Space Telescopes (MAST) 1 . We used the pre-search data-conditioned simple aperture photometry PDCSAP light curve Stumpe et al. 2012Stumpe et al. , 2014, which corrects for the systematics of the light curves by removing trends that are common to all stars in the same CCD. We removed points with a quality flag other than zero, as well as points that deviated more than 5 σ from a smooth version of the light curve. The individual sector light curves were normalised separately before being combined to produce the final light curve.
Four planet candidates in HD 23472 were detected by the TESS Science Processing Operations Center (SPOC) using a wavelet-based adaptive matched filter (Jenkins 2002;Jenkins et al. 2010Jenkins et al. , 2020. The data validation reports (Twicken et al. 2018;Li et al. 2019) were reviewed by the TESS Science Office (TSO) and issued alerts in May 2019 (Guerrero et al. 2021). The four candidates have orbital periods of P b " 17.7 days, P c " 29.8 days, P d " 3.98 days, and P f " 12.2 days as well as radii R b " 1.9 R C , R c " 2.1 R C , R d " 0.8 R C , and R f " 1.2 R C . A subsequent search of sectors 1-34 by the SPOC revealed a 5th planetary signature at a period of 7.908 days which was alerted by the TSO on 21 October 2021. The TESS light curve shows no clear rotational modulation variability, and the GLS of the light curve shows no peak at the 40 days activity signal seen in the RVs, neither half nor double of this period.

Stellar parameters
We combined the individual S1D ESPRESSO spectra after correcting for their radial velocities. The combined spectrum was then used to derived the stellar atmospheric parameters (T eff , log g, microturbulence, [Fe/H]) and the respective uncertainties using ARES+MOOG, following the same methodology described in Sousa et al. (2021); Sousa (2014); Santos et al. (2013). The analysis starts with the measurement of the equivalent widths (EW) of iron lines using the ARES code 2 (Sousa et al. 2007(Sousa et al. , 2015. We used a minimisation process to find ionisation and excitation equilibrium and converge to the best set of spectroscopic parameters. This process makes use of a grid of Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973). The line-list used for this analysis was taken from Tsantaki et al. (2013), which is more reliable for stars with effective temperature below 5200 K. The values derived for the temperature, log g, [Fe/H], and microturbulence are given in Table 1. Following the same methodology as described in Sousa et al. (2021), we used the distance derived from the GAIA eDR3 parallaxes (Gaia Collaboration et al. 2021) and estimated the trigonometric surface gravity to be 4.53˘0.06 dex and we adopted this value.
The stellar abundances of Mg and Si were derived using the classical curve-of-growth analysis method (Griffin & Griffin 1967;Adibekyan et al. 2012) and assuming local thermodynamic equilibrium. We used the same tools and models as for the stellar parameter determination. Although the EWs of the spectral lines were automatically measured with ARES, we performed a careful visual inspection of these measurements to ensure that systematic noise sources, such as cosmic rays, do not affect the results. For the derivation of the abundance values, we followed the methods described in, for instance, Adibekyan et al. (2012Adibekyan et al. ( , 2015. The mean log R 1 HK was derived by first extracting the S Ca II index from the co-added spectra with ACTIN 3 (Gomes da Silva et al. 2018), converting it to the S MW scale using the calibration for ESPRESSO in pyrhk 4 and calibrated to R 1 HK via the methodology described in Gomes da Silva et al. (2021).
To determine the mass and radius of HD 23472, we used the Bayesian tool PARAM (da Silva et al. 2006;Rodrigues et al. 2014Rodrigues et al. , 2017, in which a set of observed quantities (namely, T eff , [Fe/H], luminosity) are matched to a well-sampled grid of stellar evolutionary tracks. The stellar luminosity is inferred from the 2MASS K s band (Cutri et al. 2003), corrected by distance (Gaia parallax) and the bolometric correction. The latter is estimated by YBC (Chen et al. 2019) that interpolate the observed effective temperature, log g, and metallicity within a series of spectra libraries. Although extinction was considered in the calculation (estimated by STILISM Lallement et al. 2014;Capitanio et al. 2017), due to the proximity of the star, its impact on the luminosity is negligible with respect to overall statistical error (0.4% against 15%). As in Rodrigues et al. (2014), the grid of stellar evolutionary tracks and isochrones is from the PARSEC 5 code (Bressan et al. 2012).
For comparison, we computed the stellar parameters using the SteParSyn code 6 (Tabernero et al. 2022). The code im-plements the spectral synthesis method with an MCMC sampler to retrieve the stellar atmospheric parameters. We employed a grid of synthetic spectra computed with the Turbospectrum (Plez 2012) code alongside MARCS stellar atmospheric models (Gustafsson et al. 2008) and the atomic and molecular data of the Gaia-ESO line list (Heiter et al. 2021). We employed a set of Fe i,ii lines that are well suited for the analysis of FGKM stars (Tabernero et al. 2022). In addition, STEPARSYN allowed us to compute the following stellar atmospheric parameters: T eff " 4825˘120 K, log g " 4.70˘0.15 dex, [Fe/H] "´0.13˘0.1 dex, including systematic errors. Using these stellar parameters, we used PARAM to derive a stellar mass of 0.72˘0.04 M d and a stellar radius of 0.710 .03 R d . These values are within 1 σ of the values derived with ARES+MOOG. We adopted the values of ARES+MOOG, which is the method used for most of the planetary systems analysed by the ESPRESSO GTO.
Assuming the relation between log R 1 HK and the stellar rotation period derived by (Suárez Mascareño et al. 2015) we estimated the stellar rotation to be P rot " 44˘7days. From the v sin i˚and the stellar radius, we can estimated the rotation period of the star to be 25.1˘3.0 sini days. Hence, to match the rotation period, we would need a stellar inclination of 40`9 .4 6.8 degrees.

Light curve analysis: A new planet
We performed our own transit search for additional transiting planets in the system. We used a method similar to that of Barros et al. (2016), which was further optimised for multi-planetary systems search. We began by applying a spline filter with breakpoints every 0.5 days to correct long-term variability in each sector. Then we performed a transit search using a box least squares algorithm (BLS) (Kovács et al. 2002). We analysed the phasefolded light curve at the period corresponding to the highest peak of the BLS periodogram to assess whether it is a good transit candidate. In any case, the duration and epoch provided by the algorithm were used to remove all points within the possible transit and close to it -two durations before and after the midtransit time. Due to the relatively high uncertainty in the period and epoch derived by the BLS algorithm, a cut wider than the transit duration is required. Depending on the complexity of the system or the noise in the light curve, this process was repeated several times. For HD 23472, the first transit signal found by our procedure is the 17.7 day period planet (planet b), the second is the 29.8 day period planet (planet c), and the third is the 12.2 day period planet (planet f). However, the fourth detection is a new transiting planet candidate with a period of 7.9 days (planet e) and the fifth is the 3.98 day period planet (planet d). The orbital period of the new planet is nearly the double of the period of the previously known candidate at 3.98 days, hence, we conducted several tests to confirm that there are indeed two planets. We discovered that when only the first four sectors of TESS are used, the fourth planet detected is planet d, and the new candidate is not detected. In our case, this is due to the fact that planet e transits just before planet d, and if planet d is discovered first, all transits of planet e are removed from the light curve when we remove transits of planet d, given our wide window cut. However, if planet e is discovered first, some of the transits of planet d remain. This could have been the reason planet e had not been previously detected. All transiting planet candidates pass the standard false positive tests (Barros et al. 2016). There is: no significant difference between odd and even transits; no detectable  (Høg et al. 2000) and K, H and J from (Cutri et al. 2003) 1 trigonometric surface gravity derived using the GAIA eDR3 parallaxes (Gaia Collaboration et al. 2021) secondary eclipse; no ellipsoidal modulation at the planet period; and the transit duration is consistent with the stellar density derived by spectroscopy. We also used the open-source Python package (FULMAR) (Rodrigues et al. 2022) which uses the Transit Least-Squares (TLS) algorithm (Hippke & Heller 2019) to retrieve and analyse the light curve and independently confirm the detection of planet e. At the time of this writing, a new candidate was presented the EXOFOP-TESS database (Ex-oFOP 2019) (P 05 " 7.9 days), supporting our new detection.
For further analysis of the transit photometry, we removed points in the light curve that are more than 1.5 transit durations away from each of the planet transits. Each transit was normalised by a linear trend computed from the out-of-transit flux close to each transit. Due to the uncertainty of the BLS-derived period and epoch, we used the ephemerides derived from the first iteration's multi-transit fit to recut the light curve and renormalise the data for our final analysis.

Radial velocity and light curve modelling
The light curve and the RVs were analysed simultaneously using the LISA code (Demangeon et al. 2018(Demangeon et al. , 2021, which uses the RadVel python package (Fulton et al. 2018) to model the RV observations and a modified version of the batman transit model 7 (Kreidberg 2015) to model the transits. The system is param- 7 The modified version of batman is available at https://github.com/odemangeon/batman. The modification prevents an error for very eccentric orbits. eterised by the systemic velocity (v 0 ), stellar density (ρ ‹ ), two limb darkening parameters for the TESS bandpass corresponding to the quadratic limb darkening law, and for each planet, the semi-amplitude of the RV signal (K), the planetary period (P), the mid-transit time (T 0 ), and the products of the planetary eccentricity by the cosine and sine of the stellar argument of periastron e cos ω, e sin ω, the planet-to-star radius ratio (r p {R ‹ ), and the cosine of the orbital inclination (cos i). We included an offset between the two RV data sets (∆RV PFS{ESPRESSO ) and for each data set we included one additive jitter parameter (σ RV,ESPRESSO , σ RV,PFS , σ T ES S . Moreover, we modelled the stellar activity seen in the RVs using a Gaussian process with a quasi periodic Kernel of the form: where A RV is the amplitude of the activity signal, τ decay is the decay timescale, P rot is the period of the activity signal usually related to the stellar rotation period (e.g., Barros et al. 2020), and γ is the periodic coherence scale (e.g. Grunblatt et al. 2015). The GP was implemented with the Python package george (Ambikasaran et al. 2015). We used uniform priors for the parameters: systemic velocity, RV offset, jitter, planet-to-star radius ratio, semi-amplitude of the RV signal, the hyper-parameters of the GP, and the impact parameter. We included a prior in the impact parameter in-stead of the inclination to ensure the planets are transiting and to help convergence. We used Gaussian priors for the stellar density given in Table 1, the planetary period, and the mid-transit time derived from the BLS analysis and the limb darkening parameters derived with the LDTK code (Parviainen & Aigrain 2015;Husser et al. 2013) for the TESS bandpass ( u 1 " 0.472˘0.056 and u 2 " 0.186˘0.055). We used a Jeffrey prior for the eccentricity and constrained the eccentricity to be smaller than 0.15. Numerical simulations (Section 5.1) show that for the system to be stable, the eccentricity of all planets must be less than 0.1. If we do not include an informative prior, the eccentricity is poorly constrained to be less than 0.3, resulting in unstable solutions in the majority of the explored parameter space. Therefore, we imposed a strong prior in the eccentricity to ensure that the stable region of the parameter space was correctly sampled.
LISA employs the Bayesian inference framework (e.g. Gregory 2005) for parameter inference by maximising the posterior probability density function. It explores the parameter space with the affine-invariant Markov-chain Monte-Carlo ensemble sampler implemented in emcee (Goodman & Weare 2010;Foreman-Mackey et al. 2013). The number of walkers is set to 2.5 times the number of fitted parameters by default. To speed up convergence, the starting parameters are randomly drawn from the prior and a pre-minimisation is performed using the Nelder-Mead simplex algorithm (Nelder & Mead 1965) implemented in the Python package scipy.optimize. The chains were checked for convergence using Geweke test (Geweke 1992) and the burning-in part of the chain was removed before merging the chains. The final clean chain was composed of 5 310 000 values. The best value for each parameter was derived from the median of the posterior distribution and the uncertainties from its 68% confidence interval. More details on the LISA fitting procedure is available in Demangeon et al. (2018Demangeon et al. ( , 2021. Figure 2 shows the best model for the five planets overploted on the RV observations corrected for stellar activity with the fitted GP model, while Figure 3 shows the best transit model overploted on the transit light curves. The five planets' transits are clearly evident, but the RV signature is only clear for the two largest and longest period planets, b and c, which have already been confirmed. Table A.1 displays the best-fit model parameters. The system has three inner planets that are similar in size to the Earth but have a lower mass. The two innermost ones are smaller and more massive, so they have higher densities than the three outermost ones. The density of the middle planet f is lower than that of the Earth. The two outermost planets have the largest radii, " 2 R C , with planet b having more than twice the mass of planet c and hence being denser. As the figures show, we obtained a good precision, that is, better than 8.7%, for the radius measurement of all five planets in the system. Our precision, however, is lower for the masses. The two outermost planets are detected at more than 4σ, planets d and e have a hint of detection at 2.7 σ, and planet f's semi-amplitude of the RV is measured at 1.9σ. As a result, for the three smaller inner planets, particularly planet f, we advise caution in interpreting the planet's composition. The relative precisions in the derived semi-amplitude is 9%, 25%, 37%, 38%, and 55% for planets b, c, d, e, and f, respectively. Comparing our results with previous estimates of the RV semi-amplitude for planet b and planet c derived using the PFS data (Teske et al. 2021), we find that our results are in better agreement with the second (RadVel) method, which includes a GP to correct stellar activity than with the first (juliet) method, which, for this target, did not include the correction of stellar activity. As can be seen in Figure A.1, the start of the ESPRESSO observations coincides with the end of the PFS observations, and there is no substantial difference in the amplitude of the stellar variability between the earlier season of PFS observations and the later season of ESPRESSO observations. Furthermore, a GLS of the PFS observations reveals a clear peak at the stellar rotation period of the star " 40.5 days. Therefore, the new ESPRESSO observations allow to improve the precision and the accuracy of the mass measurements of the outermost super-Earths in the HD 23472, as well as to measure the mass and density of planets d and e and derive a strong upper limit on the mass of planet f.

Exploring the correction of the stellar activity
The ESPRESSO DRS pipeline, as mentioned in Section 2.1, provides time series for several activity indicators (FWHM, BIS, Sindex, NaD, H α , contrast). These activity indicators assess the distortions and depth of stellar lines, which serve as proxies for stellar activity. They may also indicate the presence of blending caused by another star or an eclipsing binary. Because they are not affected by the presence of planets, they are ideal for distinguishing the periodicities caused by stellar activity. Figure 1 compares the time series and GLS periodograms of the RVs and activity indicators. The activity indicators show several statistically significant peaks at long periods, but none at the planet's periods. The third strongest peak in the RVs corresponds to significant peaks in all activity indicators, a clear sign of stellar activity. We identify this as the rotation period of the star, which is approximately 40 days. This is in agreement with the values of the rotation period of the star estimated in Section 3.1 from log R 1 HK . As mentioned above, the TESS light curve shows no clear signs of rotational modulation and no periodicity at 40 days. Hence, the light curve was not used to constrain the stellar activity signal.
Using the information contained in the activity indicators can assist in correcting for stellar activity. This was accomplished by modelling the activity indicators with a GP using the same Kernel as the RVs (equation 1). The RVs and one activity indicator are fitted simultaneously. The GPs modelling the RV and the activity indicator are independent but the value of their hyperparameters are assumed to be the same except for the amplitude. For the activity indicators, we assumed that the mean function is a constant value that we also fit. We repeated this procedure using the activity indicators FWHM, H α , and contrast, individually, and we also performed the fit using the two activity indicators, FWHM and H α , simultaneously. We found that there was no significant improvement in the activity modelling when using the activity indicators. The fitted hyper-parameters with and without activity indicators have similar values and uncertainties. We also found no improvements in the RVs' residuals. This is most likely due to two factors. First, the RVs are well sampled and their precision allows for a good constraint on the activity model's hyper-parameters. This is supported by the well constrained hyper-parameters' posteriors shown in Figure A.3. Second, the power spectra of the activity indicators and the RVs differ slightly, which could result in different best fit hyper-parameter values. For example, there may be periodicities in the activity indicators that are not present in the RVs, which could affect the hyper-parameter values. This could be related to the stellar activity of HD 23472 being plage-dominated rather than spot-dominated. Stars with spot-dominated stellar activity have been found to have RVs that are better correlated with the activity indicators. The best model is shown in green. The RVs were corrected for systemic velocity, offset between the two instruments, and stellar activity with the fitted GP model. For clarity, we also show the binned RVs and, in addition, the uncertainties of the PFS data are not displayed, but they are slightly larger than the ESPRESSO ones, as mentioned in Section 2. The residuals relative to the best-fit model are shown in the bottom panel.
A better understanding of the correlation and usefulness of RV activity indicators is of extreme importance to obtain accurate masses of exoplanets. Therefore, we decided to test the bestfit values of individual activity indicators. We performed a simultaneous fit of RV, FWHM, H α , and contrast. In this case, we did not simultaneously fit the transit light curve, but instead we included Gaussian priors for the period and the transit epoch of the five planets. For the activity indicators, we used a constant mean function whose value we fit as before as well a GP to model the covariance matrix using the same RV model as presented above. In this case the only hyper-parameter shared between the activity models for the four data sets is the rotation period. Although the amplitude of the activity models for the activity indicators is expected to be different, the values of the decay timescale (τ decay ) and the periodic coherence scale (γ) are typically considered to be the same (Demangeon et al. 2021;Suárez Mascareño et al. 2020). The values we derived for the decay timescale and the coherence scale are given in Table 2. The rotation period value was found to be P rot " 39.9`1 0.91 days, which is the same as the value for the period derived if we fit only the RVs (see Table A.1 ).
We found that the derived values for the decay timescale and the coherence scale from the FWHM, H α , and contrast are significantly different from the values derived from the RVs, in particular for the coherence scale. The coherence timescale derived for the activity indicators is more similar to each other than any of them is to the RVs. The derived value of the decay timescale for the FWHM is the closest to the derived value for the RVs agreeing within 1σ. This difference in the hyper-parameters of the activity models for the activity indicators is probably the cause for the lack of improvement in including the activity indicators in the global RV fit.
With our new insight into how the hyper-parameters vary for different activity indicators we tested the correction of stellar activity including the FWHM in our systems fit of the RV and the photometry. However, in this new test, we assume only P rot and τ decay to be the same for both the RV and the FWHM activity model and A RV and γ are independent for the RV and the FWHM activity model. The derived values for rotation period is similar to our previous results with a similar error P rot " 39.5`0 .92 0.82 days and the derived the decay timescale smaller and with half of the uncertainty of the previous result. We also confirmed that the values of coherence scale for the RV and FWHM are significantly different. These improvements, however, did not affect the derived planetary parameters. All the planetary parameters are within 0.1σ of those previously derived and have similar errors. There was no reduction of the jitter derived for the RVs or the rms of the residuals. Given that there was no improvement, we conclude that the added complexity of the model is not justified. Further insights could be obtained by probing whether this  behaviour is shared by other stars. A better understanding of the correlation and utility of RV activity indicators is essential to obtain accurate exoplanet masses and more research in this area is required. With a log R 1 HK "´5 dex, TOI-174 is a chromospherically inactive K dwarf, typical of old main sequence (MS) stars (Henry et al. 1996). Inactive K dwarfs have higher activity variability than inactive FG dwarf stars (Gomes da Silva et al. 2021) but have lower induced RV variability (Lovis et al. 2011). The mea-sured value of the amplitude of the RV variation for HD 23472 ( A RV " 2m s´1) is in agreement with the typical RV variability (ď 2.5m{s) of MS stars with similar mass and activity levels (Luhn et al. 2020). The measured rotation of "40 days, obtained independently from the activity-rotation relation and activity time series, is at the upper envelope of the rotation periods distribution of stars with similar effective temperatures and is typical of an old K dwarfs with ages of ě4.5 Gyrs (McQuillan et al. 2014;Angus et al. 2020).

ESPRESSO performance
Because ESPRESSO is a new instrument, it is interesting to compare the results of our final adopted model combining the ESPRESSO and PFS data with the results of the analysis using only ESPRESSO data. Table A.2 displays the most relevant parameters for the ESPRESSO + TESS fit with the same priors as in the final model (Table A.1).
All of the results are found to be within 1σ of the combined ESPRESSO and PFS results. However, the amplitudes of the RV variation measured using only the ESPRESSO have a slightly lower level of precision, which is expected due to the decrease in the number of RV measurements. When the PFS data are in- cluded, the precision increases by 12´9%, 35´25%, 563 7%, 57´38%, and 56´55% for planets b, c, d, e, and f, respectively. The mean of the residuals after removing the GP is 0.796 m/s for the PFS data and 0.555 m/s for the ESPRESSO data and the jitter parameter for the PFS data is also slightly larger than for the ESPRESSO data as can be seen in Table A.1.

Dynamical stability and refinement of eccentricities
In orbital period ascending order, the period ratios of the HD 23472 planets are as follows: P e {P d " 1.99, P f {P e " 1.54, P b {P f " 1.45, and P c {P b " 1.68. These pairs are somewhat close to mean motion resonances (MMRs) where P out {P in « pk`qq{k. The three innermost pairs are closer to first-order MMRs (q=1, k=1 or 2), while the outermost pair is closer to a second-order MMRs (q=2, k=3). In case of multi-planetary systems close to several two-planet MMRs, similar super-periods amongst successive pairs P " 1{ppk`qq{P out´k {P in q  Izidoro et al. 2017). In particular, the inner pair displays a period ratio less than the exact commensurability (P e {P d ă 2). This is not what we would expect for a pair of planets that are initially inside a first order MMR and evolve smoothly through tides (typically effective for orbital periods less than 10 days): for such a pair, we would expect tidal evolution to gradually increase the period ratio toward P e {P d ą 2. The period ratio is less than the exact commensurability, which could be due to a disruption early in the system's history, possibly near the end of the protoplanetary disc phase or before tides had a chance to come into play.
We first performed a short-term stability analysis of the posterior samples in the case where the eccentricity of the planets is left free. For this analysis, we used the frequency analysis stability index, which is based on the diffusion of the system's main . Blue is short-term stable while red is short-term unstable. For each map, the parameters of the rest of the system are set to the nominal value given in Table 1, except for the eccentricities that are set to zero. The white dashed lines show the .16-.84 quantiles of the posterior of the mass of each planet. Only the .84 quantile is shown for planets d and f. frequencies and is defined as (Laskar 1990(Laskar , 1993: where n p1q and n p2q are the proper mean motion of the planets, computed over the first half and the second half of the integration. We found that a majority the posterior samples in the case where the eccentricity of the planets is left free is unstable in the short term. Given the system's compactness, this might be due to the relatively large eccentricity values for the free-eccentricity solution. Because the eccentricities and some of the masses are highly uncertain, we explore these parameters to gain a better understanding of their impact on the system's stability. Figure 5 displays stability maps varying the mass and eccentricity of each planet. The orbital parameters and masses of the other planets are set to their median value in each case, with the exception of eccentricities, which are set to zero. The colour code corresponds to the maximum of the stability index given in Equation 2 for all planets, implying that the stability increases from red to blue (for more details, see Section 4.1 of Petit et al. 2018). The maps show unstable structures that are hard to interpret in the masseccentricity plane, but are typically associated with the presence of nearby MMRs (e.g. Leleu et al. 2021), the precise position and boundaries of which also depend on the masses and eccentricities (see Henrard & Lemaitre 1983). A significant fraction of the 1σ mass interval is stable for all planets, pointing to eccentricities as the main parameters indicative of a system's stability. Another common trend is the appearance of instabilities for eccentricities typically greater than " 0.1.
We performed a second MCMC, exploring the eccentricities in the [0,0.15] range for all planets, because setting upper limits on each eccentricity while forcing the other eccentricities to be equal to zero is not satisfactory. We then cut the posterior in two parts using the method presented in Stalport et al. (2022): we set the stability index threshold to -4 (orange in Fig. 5), removing, from the posterior, all trajectories that are short-term unstable (37% of the trajectories were discarded as such). When examining the new, short-term stable posterior, the 0.16 and 0.84 quantiles of most parameters remain roughly unchanged, with the exception of some eccentricities. The estimated eccentrici-ties, when accounting for the stability of the system are given in Table 3. The posterior of the eccentricities of HD 23472 b and HD 23472 f shifted towards a lower value, which is required to ensure the system's stability.
What is noteworthy, according to our final result, the three inner pairs are outside of MMR. The outer pair, on the other hand, is quite close to the 5:3 MMR, and additional monitoring could confirm its resonant state. The transit timing variations (TTVs) between resonant pairs are strongly dependent on eccentricity, which is not well constrained in our case. Assuming zero eccentricity, we estimate the peak-to-peak amplitude of TTVs for planet b to be 4 minutes and for planet c to be 6 minutes. We used the same method as Barros et al. (2022) to calculate the observed TTVs from the TESS light curve and found no significant TTVs. Our estimated errors, however, have the same amplitude as the expected TTVs, 6 minutes for planet b and 7 minutes for planet c. Therefore, we cannot rule out the existence of TTVs. Photodynamical modelling the system (e.g. Barros et al. 2015) could improve the precision of the measured TTVs and help constrain the system parameters, but this approach is beyond the scope of this paper.

Architecture
A five planet system is an excellent test of the 'peas in the pod' trends reported by Weiss et al. (2018). We investigated whether HD 23472 follows the main trends presented in Weiss et al. (2018), listed below in italic. Table 4 gives the values of the "peas in the pod" metrics. Our main conclusions are as follows: The size of adjacent planets is similar or the outer planet is larger than the inner planet: The planets in HD 23472 have similar radii and the planetary radius increases outwards from the star for the first four planets, but decreases slightly for the planet with the longest period. The orbital spacing between planets is similar: This also occurs in HD 23472, although the orbital spacing increases slightly as one moves away from the star. Smaller planets are more packed than larger planets: This is in contrast to HD 23472, where the larger outer planets are more densely packed than the inner ones (i.e. they have a smaller separation in units of mutual hill radius). The temperature difference correlates with the planet size ratios: This occurs in the HD 23472 system as well.
Our high errors in terms of mass prevent us from testing whether the planets are more similar in mass than in radius, as found by Otegi et al. (2022). Additional RV observations to improve mass precision would be extremely helpful in testing this hypothesis. HD 23472 follows the general trends of "peas in the pod" but it appears to be breaking the trends for the longest period planet. As a result, probing longer period planets in this system would be very beneficial.

Planetary composition
Our analysis enabled us to achieve relative precisions in the densities of the five planets of 18%, 31%, 46%, 46%, and 59% for planets b, c, d, e, and f, respectively. Although this result is remarkable for such low-mass planets, a higher level of precision would be necessary to uniquely characterise the composition of the planets especially the three smaller ones and we advise taking care in the interpretation of the results. Figure 6 displays the position of the five planets in the mass-radius diagram as well as the compositional models of Zeng et al. (2016) and the ra-dius gap (Fulton et al. 2017). The three inner planets are below the radius gap while the two outermost planets are above it. This could be caused by either irradiation shaping the system or corepowered evaporation. The two inner planets, d and e, are very dense, indicating that they likely contain more iron than Earth according to the models of Zeng et al. (2016). Planets b, c, and f appear to have a significant amount of water or gas in their composition, with planet c having the most. Moreover, the two inner planets are in a sparsely populated region of the parameter space near L 98-59 b, whose mass was recently measured with ESPRESSO (Demangeon et al. 2021) and Trappist-1 h (Gillon et al. 2017), whose mass was measured with transit timing variations. Finding and characterising planets in this region of the parameter space would greatly improve our understanding of the composition of small exoplanets.   Fig. 6: HD 23472 planets in the context of other known transiting planets with measured mass and radius precision better than 50 %. The exoplanets data have been extracted from the NASA exoplanet archive. Planets whose masses were determined using the RV technique are represented as circles, whereas planets whose masses were determined using transit timing variations are represented as squares. The intensity of the incident flux is indicated by the colour of the points.The planetary bulk density's relative precision is proportional to the transparency of the error bars. The planets of the solar system are illustrated as blue stars. We also show the mass-radius models of Zeng et al. (2016) as dashed lines. The radius gap (Fulton et al. 2017) is shown as a shaded horizontal blue line and the maximum collision stripping of the mantle region is shown in grey. This graph was created using the mass radius diagram code 8 . .
To constrain the internal structure of the planets, we performed a Bayesian analysis following the method of Dorn et al. (2015Dorn et al. ( , 2017. The same method was used to analyse other sys- Mode between 10 and 20 ∆p f, bq " 10.39˘0.34 ∆pe, f q " 22.3˘3.0 ∆pd, eq " 36.9˘3.9 T eq,d´Teq,e " 185˘32 K T eq,i´Teq,i`1 positively correlated with R i`1 {R i T eq,e´Teq, f " 96˘28 K T eq, f´Teq,b " 79˘21 K T eq,b´Teq,c " 80˘18 K Notes. -∆pi, jq is the separation in mutual Hill radii (see Eq. 5 in W18) When the notation x˘y is used in the column " W18 distribution", x is the median of the observed distribution, and y is its standard deviation. tems, including, L 98-59 b (Demangeon et al. 2021), TOI-178 (Leleu et al. 2021)m and Nu2 Lupi (Delrez et al. 2021). The model assumes the planet has four layers: an iron and sulphur inner core, a silicate mantle ( made of Si, Mg and Fe), a water layer, and a gas layer (made of H and He). We used an improved equation of state for the water layer  and an improved equation of state of the iron core (Hakim et al. 2018) for our analysis.
The model includes two parts: a forward model that computes the planetary radius as a function of the internal structure parameters and a Bayesian analysis that computes the posterior distribution of the internal structure parameters needed to fit the observed radii, masses, equilibrium temperatures, and stellar parameters. Since absolute planetary masses and radii are measured relative to the same host star, they are correlated. As a result, rather than fitting the absolute masses and radii, we fit the radius ratio and radial velocity semi-amplitudes of all planets in the system simultaneously. The input parameters of our model are the stellar mass, stellar radius, stellar effective temperature, stellar age, stellar chemical abundances of Fe, Mg, and Si, and the planetary radial velocity semi-amplitudes, planetary radius ratio, and orbital periods. The fitted parameters are the mass fractions of core, mantle, water layer, and gas layer. Except for the mass of the gas layer, which is assumed to follow a uniform-inlog prior, we used uniform priors for these parameters. We note that the results we obtain are influenced by these priors to some extent. Assuming a uniform prior for the gas mass, for example, would result in a large mass fraction of gas and a smaller mass fraction of water. Finally, the mass fractions of the different layers add up to one, and water mass fractions greater than 50% are excluded Marboeuf et al. 2014).
We ran a first set of models, assuming, as in other similar studies that the planets share the same Si/Mg/Fe molar ratio equal to the stellar one (see, however, Adibekyan et al. (2021)). In this case, it turns out to be impossible to fit the mass and radius of the two innermost planets (d and e), whose densities are too high to be reproduced by a core-mantle structure with stellar abundances. We then relaxed the assumption of stellar composi-tion for these two planets, letting the Si/Fe and Mg/Fe ratio free, but assuming the Si/Mg ratio to be stellar. We note that this last assumption comes from the idea that in the case of the formation of a Mercury-like planet by a giant collision that would have removed part of the mantle, it is not expected that Si and Mg would fractionate 9 . All other assumption regarding the forward model are kept the same. In this case, it is possible to reproduce the five planets in an accurate way.
The gas fraction, water fraction, and core mass fraction for the system are shown in Figure 7 (see also Table 5). These examples show a clear increase of the water mass fraction and the gas mass fraction with decreasing planetary effective temperature. Our Bayesian analysis shows that the two inner planets have a small water mass fraction and a negligible gas mass fraction. The three outermost planets, particularly planets b and c, are likely to have much more water. Similarly, the three outermost planets could have a non-negligible mass of gas, up to a few percent of an Earth mass for the two outermost ones. The two innermost planets, on the other hand, are very likely to be devoid of gas. This could be explained by the loss of volatiles and gas due to irradiation on these planets. It could also be explained if the two inner planets formed inside the snow line and are dry, while the three outer planets formed beyond the snow line and are water worlds (Venturini et al. 2020;Luque & Pallé 2022) . The two innermost planets have a large iron core (on the order of 30 to 60 % in mass) and have a Si/Fe and Mg/Fe ratio smaller than the stellar one, whereas the three outer planets have smaller cores (less than 20 % in mass). This dichotomy in the structure of the planets (the two innermost ones versus the three outermost ones) is remarkable and puzzling in terms of identifying their formation processes (for more, see the next section). , and core mass fraction (bottom right panel) for the five planets orbiting HD 23472. The green triangle represents the distribution mean, the orange line represents the distribution median, the red star represents the distribution mode, and the box represents the 25 and 75 percentiles. The opacity of the black line (outside the boxes) is proportional to the posterior distribution (normalised to the mode). Table 5: Internal structure parameter for the five planets. 'Core' is the core mass fraction, 'Mantle' the mantle mass fraction, 'Water' the water mass fraction, 'Gas' the log (base 10) of the mass of gas (in Earth masses), 'Si/Fe' the Si/Fe bulk molar fraction and 'Mg/Fe/ the Mg/Fe bulk molar fraction. The uncertainties correspond to the 5 and 95 % percentile, while the central value corresponds to the median.

Two likely super-Mercuries in HD 23472
Recently, Adibekyan et al. (2021) reported that all super-Mercuries are formed in proto-planetary disks with enhanced iron abundance when compared to Mg and Si. Following Adibekyan et al. (2021) and using the stellar abundances, we estimated the iron-to-silicate-mass fraction of the proto-planetary disc of HD 23472. We found that it is high (32.3 +/-3.9%), despite the star being relatively metal-poor. We noticed that this occurs because the abundances of Si and Mg in the star are similar to the abundance of Fe. This is uncommon since most stars with low iron content are slightly enhanced in Si and Mg relative to iron due to galactic chemical evolution (Adibekyan et al. 2012) . We also derived the planetary density normalised to the density of an Earth-like composition to check whether the planets of HD 23472 follow the correlation reported by Adibekyan et al. (2021) between iron mass fraction of the star and the normalised density of super-Earths. The normalisation accounts for the increase in density with mass for planets with the same composition. Figure 8 compares the scaled planet density with the iron mass fraction of the star for planets HD 23472 d, HD 23472 e, and HD 23472 f within the context of known small planets relation proposed by Adibekyan et al. (2021). Planets b and c have significant water-or gas-rich envelopes (or a combination of both) and are above the radius gap; hence, they are not expected to follow the correlation. Planet f is below the radius gap but it might not follow the correlation if it has a significant gas rich envelope (see Fig. 7). The two inner planets are in the region of the super-Mercuries, which have higher densities than what would be expected from the host star composition. However, within the relative large mass uncertainties, they are also compatible with the Super-Earth population. A better level of precision for the planetary mass is needed to confirm that both planets are super-Mercuries. If the presence of two super-Mercuries in HD 23472 is confirmed, it would make this an excellent test-bed for theories of super-Mercuries formation and evolution, and it may hold the key to solving their mystery. In the context of the giant impact theory, it has been shown that forming super-Mercuries would require a series of strong giant impacts (Scora et al. 2020). This implies that the formation of super-Mercuries via giant impact is rare, and the presence of two in the same system is extremely unlikely. The gap in planetary density between super-Mercuries and super-Earths reported by Adibekyan et al. (2021) also calls into question the hypothesis of giant impact as its stochastic nature predicts a continuous distribution of densities.
The theory of mantle evaporation was proposed to explain the high density of Mercury, whose high dayside temperature would be sufficient to cause mantle evaporation into an atmosphere of silicate vapour (Cameron 1985;Perez-Becker & Chiang 2013). However, a very high rate of evaporation is required, which is not supported by theory. Furthermore, this mechanism should be applicable to all super-Earths with similar equilibrium temperatures, which has not been observed .
Photophoresis, that is, the depletion of silicates at the inneredge of the proto-planetary disc, could also be responsible for the formation of super-Mercuries. Other mechanisms that also change the conditions of the inner disc, such as, rocklines (Aguichine et al. 2020), magnetic erosion (Hubbard 2014), and magnetic boost (Kruss & Wurm 2018), have also been proposed to form super-Mercuries and might also be able two form two super-Mercuries in the same system.
Finally, another possibility is that super-Mercuries are planetary cores that have 'recently' lost their envelopes due to a significant mass-loss event and are still compressed (Mocquet et al. 2014). Characterising these cores would provide us with unique information about the planet's interiors.

Conclusions
We obtained RV observations with the ESPRESSO spectrograph mounted in the VLT to measure the masses of the five planets around HD 23472. Combining our new observations with previous PFS RV observations and TESS photometry, we estimated the composition of the five planets in the system. We found slightly smaller masses for the two exoplanets that were previously confirmed, namely, planets b and c. We also constrained the mass of the other three smaller inner planets in the system. The two outermost planets are approximately twice the size of the Earth, and the mass of planet b, M b " 8.42`0 .83 0.84 M C , is more than twice the mass of planet c, M c " 3.37`0 .92 0.87 M C ; hence, it is much denser and has a lower gas and water content than planet c. The middle planet (planet f) is slightly larger than the Earth and the lightest of the planets, weighing approximately half as much as the Earth. The semi-amplitude of the RV signature of this planet was detected at 1.9σ (55% relative precision) with an upper limit on the mass of 1.5 M C at 95% confidence level. Hence, more observations are required to confirm its low density and large water and gas mass fraction. The two inner planets are both smaller and lighter than Earth. The semi-amplitude of the RV signature of the two inner planets was detected at higher significance (2.7σ) but more observations are also advisable to confirm their derived high density. We show that their high density and properties match those of super-Mercuries previously discovered. Using Bayesian internal structure model we find that the two inner planets have much higher iron core than the outer planets and they have a Si/Fe and Mg/Fe ratio smaller than the stellar one consistent with super-Mercuries composition. Further RV observations to improve the precision of the planetary masses would greatly improve our understanding of this unique system. If they confirm that the composition of the two inner planets, this would be the first time two super-Mercuries have been discovered in the same system, making it a golden target for further characterisation.
Atmosphere observations may shed light on the formation of the two super-Mercuries as well as the system architecture. We computed the transmission spectrum metric (TSM, Kempton et al. 2018) to access the observability with the James Webb space telescope. Planets b, c, d, e, and f have TSMs of 36, 59, 7, 5, and 14, respectively. It is only planet f that has a TSM greater than the threshold proposed by Kempton et al. (2018). The TSM values for the five planets of HD 23472 system are shown in Figure 9, set in the context of known and well-characterised exoplanets. Planets b and c are good targets for transmission spectroscopy due to their low effective temperature T eq ă 550 K. As expected, the two likely super-Mercuries are difficult targets for transmission spectroscopy, but they may be better targets for emission spectroscopy using, for example, the future spectrograph ANDES@ELT. Furthermore, HD 23472 is the brightest star with super-Mercuries (H mag = 7.3), making this the best system for studying a potential atmosphere around Super-Mercuries. For an Earth-mass planet orbiting HD 23472, the optimistic limits of the habitable zone correspond to periods between 98 days and 355 days. If the near-resonant chain continues to longer orbital periods and the system contains two longer period planets (g and h), planet h may be on the edge of the optimistic habitable zone. As a result, HD23472 is shown to be an excellent candidate for the search for habitable Earths. We also encourage monitoring studies of the transit times of the two outermost planets to assess whether they are in MMR. . Fig. 9: Transmission spectrum metric as a function of planetary radius for the five planets in the HD 23472 and the well characterised small planet population. We show all planets from the exoplanet archive with precision in mass and radius better than 50%Ṫhe planets with mass derived by RVs are shown as circles and the planets with mass derived by TTVs are shown as squares.  Table A.3. Figure A.1 depicts the RV time series, as well as the best-fit five-planet model and the GP activity model. The relative high RV variability due to stellar activity is clearly evident. We also show a zoom of the plot during the overlap of ESPRESSO and PFS observations, demonstrating good agreement between the two instruments.
In Figure A.2, we compare the GLS of the RV data with the GLS of the five Keplerian model, the GP model, and the residuals. With the exception of planet c, which is modelled by the Keplerian model, the GP aptly models all the periodicities longer than half the rotation period of the star (P rot " 40 days). The residuals show low power periodicities at shorter time scales, which could be due to short timescale stellar variability such as granulation. Figure A.3 shows the corner plots of the GP model's hyperparameters. The posteriors of the GP hyper-parameter demonstrate that the RV aptly data constrain hyper-parameters.  Notes. ‚ indicates that the parameter is a main or jumping parameter for the mcmc explorations  Up0, 400q Notes. Upa; bq is a uniform distribution between a and b; Jpa; bq is a Jeffreys distribution between a and b; Npa; bq is a normal distribution with mean a and standard deviation b.
S. C. C. Barros et al.: HD 23472: A multi-planetary system with three super-Earths and two potential super-Mercuries