Company for the Ultra-high Density, Ultra-short Period Sub-Earth GJ 367 b: Discovery of Two Additional Low-mass Planets at 11.5 and 34 Days

GJ 367 is a bright (V ≈ 10.2) M1 V star that has been recently found to host a transiting ultra-short period sub-Earth on a 7.7 hr orbit. With the aim of improving the planetary mass and radius and unveiling the inner architecture of the system, we performed an intensive radial velocity follow-up campaign with the HARPS spectrograph—collecting 371 high-precision measurements over a baseline of nearly 3 yr—and combined our Doppler measurements with new TESS observations from sectors 35 and 36. We found that GJ 367 b has a mass of M b = 0.633 ± 0.050 M ⊕ and a radius of R b = 0.699 ± 0.024 R ⊕, corresponding to precisions of 8% and 3.4%, respectively. This implies a planetary bulk density of ρ b = 10.2 ± 1.3 g cm−3, i.e., 85% higher than Earth’s density. We revealed the presence of two additional non-transiting low-mass companions with orbital periods of ∼11.5 and 34 days and minimum masses of Mcsinic = 4.13 ± 0.36 M ⊕ and Mdsinid = 6.03 ± 0.49 M ⊕, respectively, which lie close to the 3:1 mean motion commensurability. GJ 367 b joins the small class of high-density planets, namely the class of super-Mercuries, being the densest ultra-short period small planet known to date. Thanks to our precise mass and radius estimates, we explored the potential internal composition and structure of GJ 367 b, and found that it is expected to have an iron core with a mass fraction of 0.91−0.23+0.07 . How this iron core is formed and how such a high density is reached is still not clear, and we discuss the possible pathways of formation of such a small ultra-dense planet.


INTRODUCTION
Close-in planets with orbital periods of a few days challenge planet formation and evolution theories and play a key role in the architecture of exoplanetary systems (Winn & Fabrycky 2015;Zhu & Dong 2021).
To date, about 132 ultra-short period (USP) planets, namely planets with orbital periods shorter than 1 day (Sahu et al. 2006;Sanchis-Ojeda et al. 2014), have been validated, and only 36 of these were confirmed and have measured radii and masses 1 .USPs are preferred targets for transit and radial velocity (RV) planet search surveys, as the transit probability is higher -it scales as P −2/3 orb -and the Doppler reflex motion is largerit scales as P −1/3 orb .In addition, their orbital period is typically 1 order of magnitude shorter than the rotation period of the star, allowing one to disentangle bona fide planetary signals from stellar activity (Hatzes et al. 2011;Hatzes 2019;Winn et al. 2018).Sanchis-Ojeda et al. (2014) found that the occurrence rate of rocky USP planets seems to depend on the spectral type of the host star, being 0.15 ± 0.05 % for F dwarfs, 0.51 ± 0.07 % for G dwarfs, and 1.10 ± 0.40 % for M dwarfs.In this context, low-mass stars, such as M dwarfs, are particularly suitable to search for closein terrestrial planets.Given the relatively small stellar † elisa.goffo@unito.itelisa@tls-tautenburg.de ‡ NASA Sagan Fellow.
radius and mass, a planet transiting an M dwarf star induces both a deeper transit and a larger RV signal, increasing its detection probability (Cifuentes et al. 2020).
The formation process of USP planets is still not fully understood, and different scenarios have been proposed to explain their short-period orbits: dynamical interactions in multi-planet systems (Schlaufman et al. 2010); low-eccentricity migration due to secular planet-planet interactions (Pu & Lai 2019); high-eccentricity migration due to secular dynamical chaos (Petrovich et al. 2019); tidal orbital decay of USP planets formed in situ (Lee & Chiang 2017); and obliquity-driven tidal migration (Millholland & Spalding 2020).Intensive follow-up observations of systems hosting USP planets can help us to understand the formation and evolution mechanisms of short-period objects and other phenomena related to star−planet interactions (Serrano et al. 2022).Sanchis-Ojeda et al. (2014), Adams et al. (2017), and Winn et al. (2018) found that most USP planets have nearby planetary companions.In multi-planet systems, USP planets show wider-than-usual period ratios with their nearest companion, and appear to have larger mutual inclinations than planets on outer orbits (Rodriguez et al. 2018;Dai et al. 2018).These observations suggest that USP planets experienced a change in their orbital parameters, such as inclination increase and orbital shrinkage, suggesting the presence of long-period planets.
During its primary mission, NASA's Transiting Exoplanet Survey Satellite (TESS ; Ricker et al. 2015) discovered a shallow (∼300 ppm) transit event repeating every ∼7.7 hr (∼0.32 days) and associated to a USP small planet candidate orbiting the bright (V ≈ 10.2), nearby (d ≈ 9.4 pc), M1 V star GJ 367.Lam et al. (2021) recently confirmed GJ 367 b as a bona fide USP sub-Earth with a radius of R b = 0.718 ± 0.054 R ⊕ and a mass of M b = 0.546 ± 0.078 M ⊕ .
As the majority of well-characterized USP systems are consistent with having additional planetary companions (Dai et al. 2021), it is quite realistic to believe that the GJ 367 hosts more than one planet.As part of the RV follow-up program carried out by the KESPRINT consortium2 , we here present the results of an intensive RV campaign conducted with the HARPS spectrograph to refine the mass determination of the transiting USP planet and search for external planetary companions, while probing the architecture of the GJ 367 planetary system.
The paper is organized as follows: we provide a summary of the TESS data and describe our HARPS spectroscopic follow-up in Sect. 2. Stellar fundamental parameters are presented in Sect.3. We report on the RV and transit analysis in Sects.4 and 5, along with the frequency analysis of our HARPS time series.Discussion and conclusions are given in Sects.6 and 7, respectively.
2. OBSERVATIONS 2.1.TESS Photometry TESS observed GJ 367 in Sector 9 as part of its primary mission, from 2019 February 28 to 2019 March 26, with CCD 1 of camera 3 at a cadence of 2 minutes.These observations have been presented in Lam et al. (2021).About 2 yr later, TESS re-observed GJ 367 as part of its extended mission in Sectors 35 and 36, from 2021 February 9 to 2021 April 2, with CCD 1 and 2 of camera 3 at a higher cadence of 20s as well as at 2 minutes.The photometric data were processed by the TESS data processing pipeline developed by the Science Processing Operations Center (SPOC; Jenkins et al. 2016).The SPOC pipeline uses Simple Aperture Photometry (SAP) to generate stellar light curves, where common instrumental systematics are removed via the Presearch Data Conditioning (PDCSAP) algorithm developed for the Kepler space mission (Stumpe et al. 2012;Stumpe et al. 2014;Smith et al. 2012).
We retrieved TESS Sector 9, 35, and 36 data from from the Mikulski Archive for Space Telescopes (MAST)3 and performed our data analyses using the PDCSAP light curve.We ran the Détection Spécialisée de Transits (DST; Cabrera et al. 2012) algorithm to search for additional transit signals and found no significant detection besides the 7.7 h signal associated to GJ 367 b, suggesting that there are no other transiting planets in the system observed in TESS Sector 9, 35, and 36, consistent with the SPOC multi-transiting planet search.
The exposure time varied between 600 and 1200 s, depending on weather conditions and observing schedule constraints, leading to a signal-to-noise ratio (S/N) per pixel at 550 nm ranging between 20 and 90, with a median of ∼55.We used the second fiber of the instrument to simultaneously observe a Fabry-Perot interferometer and trace possible nightly instrumental drifts (Wildi et al. 2010(Wildi et al. , 2011)).The HARPS data were reduced using the Data Reduction Software (DRS; Lovis & Pepe 2007) available at the telescope.The RV measurements, as well as the Hα, Hβ, Hγ, Na D activity indicators, log R ′ HK , the differential line width (DLW), and the chromaticity index (CRX), were extracted using the codes NAIRA (Astudillo-Defru et al. 2017b) and serval (Zechmeister et al. 2017).NAIRA and serval feature template matching algorithms that are suitable to derive precise RVs for M-dwarf stars, when compared to the cross-correlation function technique implemented in the DRS.We tested both the NAIRA and serval RV time series and found no significant difference in the fitted parameters.While we have no reason to prefer one code over the other, we used the RV data extracted with NAIRA for the analyses described in the following sections.
Table 1 lists the HARPS RVs, including those previously reported in Lam et al. (2021), along with the activity indicators and line profile variation diagnostics extracted with NAIRA and serval.Time stamps are given in Barycentric Julian Date in the Barycentric Dynamical Time (BJD TDB ).

HARPS spectropolarimetric observations
With the aim of measuring the magnetic field of GJ 367, we performed a single circular polarization observation with the HARPSpol polarimeter (Piskunov et al. 2011;Snik et al. 2011) on 2022 November 16 (UT), as part of the our ESO HARPS program 1102.C-0923 (PI: Gandolfi).We used an exposure time of 3600 s split in four T exp = 900 s sub-exposures obtained with different configuration of polarization optics to ensure cancellation of the spurious instrumental signals (see Donati et al. 1997;Bagnulo et al. 2009).The data reduction was carried out with the REDUCE code (Piskunov & Valenti 2002) following the steps described in Rusomarov et al. (2013).The resulting Stokes I (intensity) and Stokes V (circular polarization) spectra cover approximately the same wavelength interval as the usual HARPS observations at a slightly reduced resolving power (R ≈ 110, 000).We also derived a diagnostic null spectrum (e.g.Bagnulo et al. 2009), which is useful for assessing the presence of instrumental artifacts and non-Gaussian noise in the Stokes V spectra.

Photospheric and fundamental parameters
We derived the spectroscopic parameters using the new co-added HARPS spectrum for GJ 367.Following the prescription in Hirano et al. (2018), we first estimate the stellar effective temperature T eff , metallicity [Fe/H], and radius R ⋆ using SpecMatch-Emp (Yee et al. 2017).The code attempts to find a subset of best-matching template spectra from the library to the input spectrum and derives the best empirical values for the above parameters.SpecMatch-Emp returned T eff = 3522 ± 70 K, R ⋆ = 0.452 ± 0.045 R ⊙ , and [Fe/H] = −0.01 ± 0.12.We then used those parameters to estimate the other stellar parameters as well as refine R ⋆ .As described in Hirano et al. (2021), we implemented a Markov Chain Monte Carlo (MCMC) simulation and derived the parameters in a self-consistent manner, making use of the empirical formulae by Mann et al. (2015) and Mann et al. (2019) for the derivations of the stellar mass M ⋆ and radius R ⋆ .We found that GJ 367 has a mass of M ⋆ = 0.455 ± 0.011 M ⊙ and a radius of R ⋆ = 0.458 ± 0.013 R ⊙ , with the latter in very good agreement with the value derived using SpecMatch-Emp.In the MCMC implementation, we also derived the stellar density ρ ⋆ , surface gravity log g ⋆ , and luminosity L ⋆ of the star.Results of our analysis are listed in Table 1.As the quoted uncertainties of the stellar parameters do not account for possible unknown systematic errorswhich in turn might affect the estimates of the planetary parameters -we performed a sanity check and determined the stellar mass and radius using the code BASTA (Aguirre Børsen-Koch et al. 2022).We fitted the derived stellar parameters to the BaSTI isochrones (Hidalgo et al. 2018, science case 4 in Table 1 of the BASTA paper).Starting from the effective temperature, metallicity, and radius, as derived from SpecMatch-Emp, yielded a mass of M ⋆ = 0.435 +0.035 −0.040 M ⊙ and radius of R ⋆ = 0.411 +0.032 −0.034 R ⊙ .The stellar mass is in good agreement with the previous result.The new estimate of the stellar radius is smaller than the value reported in Table 1, but it is still consistent within ∼1.4 σ (where σ is the sum in quadrature of the two nominal uncertainties) with a p-value of ∼15%.Assuming a significance level of 5%, the two radii are consistent, providing evidence that our estimates might not be significantly affected by inaccuracy.

Rotation period
Using archival photometry from the Wide Angle Search for Planets survey (WASP), Lam et al. (2021) found a photometric modulation with a period of 48 ± 2 days.Lam et al. (2021) (2017a), they estimated a stellar rotation period of P rot = 58.0 ± 6.9 days.
We independently derived a Ca II H & K chromospheric activity index of log R ′ HK = -5.169± 0.068 from the 371 HARPS spectra and estimated the rotation period of GJ 367 using the same empirical relationship.We found a rotation period of P rot = 54 ± 6 d, in good agreement with the previous estimated value.We note that our estimate is consistent within 1σ with the value of P rot = 51.30± 0.13 d recovered by our sinusoidal signal analysis described in Sect.5.2, and with the period of P rot = 53.67 +0.65  −0.53 d derived by our multidimensional Gaussian process (GP) analysis described in Sect.5.3.

Magnetic field
Our spectropolarimetric observation of GJ 367 achieved a median S/N of about 90 over the red HARPS chip.This is insufficient for detecting Zeeman polarization signatures in individual lines even for the most active M dwarfs.To boost the signal, we made use of the least-squares deconvolution procedure (LSD, Donati et al. 1997) as implemented by Kochukhov et al. (2010).The line mask required for LSD was obtained from the VALD database (Ryabchikova et al. 2015) using the atmospheric parameters of GJ 367 and assuming solar abundances (Sect.3.1).We used about 5000 lines deeper than 20% of the continuum for LSD, reaching a S/N of 7250 per 1 km s −1 velocity bin.The resulting Stokes V profile has a shape compatible with a Zeeman polarization signature with an amplitude of ≈0.04% (Figure 1).However, with a false-alarm probability (FAP) = 2.3%, detection of this signal is not statistically significant according to the usual detection criteria employed in high-resolution spectropolarimetry (Donati et al. 1997).The mean longitudinal magnetic field, which represents the disk-averaged line-of-sight component of the global magnetic field, derived from this Stokes V profile is Our spectropolarimetric observation of GJ 367 suggests that this object is not an active M dwarf.Its longitudinal magnetic field was found to be below 10 G, which is much weaker than ≥ 100-700 G longitudinal fields typical of active M dwarfs (Donati et al. 2008;Morin et al. 2008).Considering that the strength and topology of the global magnetic fields of M dwarfs is systematically changing with the stellar mass (Kochukhov 2021), it is appropriate to compare GJ 367 with early-M dwarfs.To this end, the well-known ∼20 Myr old M1V star AU Mic was observed with HARPSpol in the same configuration and with a similar S/N as our GJ 367 observations (Kochukhov & Reiners 2020), yielding consistent detections of polarization signatures and longitudinal fields of up to 50 G.The magnetic activity of GJ 367 is evidently well below that of AU Mic.

Age
Determining the age of M dwarf stars is especially challenging and depends on the methods being utilized.Lam et al. (2021) estimated an isochronal age of 8.0 +3.8  −4.6 Gyr and a gyrochronological age of 4.0 ± 1.1 Gyr for GJ 367.More recently, Brandner et al. (2022) gave two different estimates: (1) by comparing Gaia EDR3 parallax and photometric measurements with theoretical isochrones, they suggested a young age < 60 Myr.However, as pointed out by the authors, this is not in line with the star's Galactic kinematics that exclude membership to any nearby young moving group; (2) by considering the Galactic dynamical evolution, which indicates an age of 1-8 Gyr.
In this respect, the results presented in our study shows compelling evidence that GJ 367 is not a young star: • The time series of our HARPS RVs and activity indicators give a clear detection of a spot-induced rotation modulation with a period of about 52−54 d, which translates into a gyrochronological age of ∼4.6-4.8Gyr (Barnes 2010; Barnes & Kim 2010).
• The HARPS spectra of GJ 367 show no significant lithium absorption line at 6708 Å. Figure 2 displays the co-added HARPS spectrum of AU Mic (Zicher et al. 2022;Klein et al. 2022) in the spectral region around the Li i 6708 Å line.AU Mic is a well-known 22-Myr-old M1 star located in the β Pictoris moving group (Mamajek & Bell 2014).Superimposed with a thick red line is GJ 367's coadded HARPS spectrum, which has been broadened to match the projected rotational velocity of AU Mic (v sin i ⋆ = 7.8 km s −1 ).Since lithium is quickly depleted in young GKM stars, the lack of lithium in the spectrum of GJ 367 suggests an age ≳ 50 Myr (Binks & Jeffries 2014;Binks et al. 2021).
• The low level of magnetic activity inferred by the Ca H & K indicator log R ′ HK (Sect, 3.2) and the weak magnetic field (Sect.3.3) are consistent with an old, inactive M-dwarf scenario (Pace 2013).
• We measured an average Hα equivalent width of EW = 0.0638 ± 0.0014 Å.Using the empirical rela- Superimposed with a thick red line is the co-added HARPS spectrum of GJ 367, which has been rotationally broadened to match the v sin i⋆ = 7.8 km s −1 of AU Mic.
tion that connects the Hα equivalent width with stellar age (Kiman et al. 2021), this translates into an age ≳ 300 Myr.
We therefore conclude that GJ 367 is a rather slowly rotating, old star with a low magnetic activity level, rather than a young M dwarf.This conclusion is consistent with a recent study by Gaidos et al. (2023), who measured an age of 7.95 ± 1.31 Gyr from the M dwarf rotation−age relation.

FREQUENCY ANALYSIS OF THE HARPS TIME SERIES
In order to search for the Doppler reflex motion induced by the USP planet GJ 367 b and unveil the presence of potential additional signals associated with other orbiting companions and/or stellar activity, we performed a frequency analysis of the HARPS RV measurements and activity indicators.For this analysis we did not include the HARPS measurements4 presented in Lam et al. (2021) and used only the Doppler data collected between 2020 November 9 and 2022 April 18 (UT) as part of our HARPS large program, to avoid spurious peaks introduced by the poor sampling of the existing old data set, and to avoid having to account for the RV offset caused by the refurbishment of the instrument.
Figure 3 shows the generalized Lomb−Scargle (GLS; Zechmeister & Kürster 2009) periodograms of the HARPS RVs and activity indicators in two fre-quency ranges, i.e. 0.000 − 0.130 d −1 (left panels) and 3.075 − 3.125 d −1 (right panels), with the former including the frequencies at which we expect to see activity signals at the rotation period of the star, and the latter encompassing the orbital frequency of the transiting planet GJ 367 b.The horizontal dashed lines mark the GLS powers at the 0.1%, 1%, and 10% FAP.The FAP was estimated following the bootstrap method described in Murdoch et al. (1993), i.e., by computing the GLS periodogram of 10 6 mock time series obtained by randomly shuffling the measurements and their uncertainties, while keeping the time stamps fixed.In this work we assumed a peak to be significant if its FAP < 0.1%.
The GLS periodogram of the HARPS RVs (Fig. 3, upper panel) shows its most significant peak at f 1 = 0.086 day −1 , corresponding to a period of about 11.5 days.This peak is not detected in the activity indicators, providing evidence that the 11.5 day signal is caused by a second planet orbiting the host star, hereafter referred to as GJ 367 c.
We used the pre-whitening technique (Hatzes et al. 2010) to identify additional significant signals and successively remove them from the RV time series.We employed the code pyaneti (Barragán et al. 2019;Barragán et al. 2022) to subtract the 11.5 day signal from the HARPS RVs assuming a circular model, adopting uniform priors centered around the period and phase derived from the periodogram analysis, while allowing the systemic velocity and RV semi-amplitude to uniformly vary over a wide range.
The periodogram of the RV residuals following the subtraction of the signal at f 1 (Fig. 3, second panel) shows its most significant peak at f 2 = 0.019 day −1 (52.2 days).Iterating the pre-whitening procedure and removing the signal at f 2 , we found a significant peak at f 3 = 0.029 day −1 (34 days).This peak has no significant counterpart in the activity indicators, suggesting it is associated to a third planet orbiting the star, hereafter referred as GJ 367 d (Fig. 3, third panel).The periodograms of the CRX, DLW, Hα, Hβ, Hγ, log R ′ HK , and Na D show significant peaks in the range 48−54 days (Fig. 3, lower panels), i.e., close to the rotation period of the star, suggesting that the peak at f 2 = 0.019 day −1 (52.2 days) seen in the RV residuals is very likely associated to the presence of active regions appearing and disappearing on the visible stellar hemisphere as the star rotates about its axis.
After removing the signal at f 3 , we found significant power at f 4 = 0.009 day −1 (∼115 days; Fig. 3, fourth panel).The periodograms of the activity indicators show also significant power around f 4 , providing evidence this signal is associated to stellar activity.As we will discuss in Sect.5.3, we interpreted the power at f 4 as the evolution timescale of active regions.Finally, we found that the RV residuals also show a significant fifth peak at f 5 = f b = 3.106 day −1 (0.322 days), the orbital frequency of the USP planet GJ 367 b, further confirming the planetary nature of the transit signal identified in the TESS data and announced by Lam et al. (2021).
Finally, we computed the GLS periodogram of the HARPS RVs including all of the available data, i.e., the data acquired before and after the refurbishment of the instrument.Adding the old data points increases the baseline of our observations and, consequently, the frequency resolution.However, the resulting periodogram is "jagged", owing to the presence of aliases with very small frequency spacing, making it more difficult to identify the true peaks.

DATA ANALYSIS
We modeled the TESS transit light curves and HARPS RV measurements using three different approaches.The methods differ in the way the Doppler data are fitted, as described in the following subsections.

Floating chunk offset method
We used the floating chunk offset (FCO) method to determine the semi-amplitude K b of the Doppler reflex motion induced by the USP planet GJ 367 b.Pioneered by Hatzes et al. (2011) for the mass determination of the USP planet CoRoT-7 b, the FCO method relies on the reasonable assumption that, within a single night, the RV variation of the star is mainly induced by the orbital motion of the USP planet rather than stellar rotation, magnetic cycles, or orbiting companions on longerperiod orbits.As the RV component due to long-period phenomena can thus be treated as constant within a given night, introducing nightly offsets filters out any long-term RV variations, allowing one to disentangle the reflex motion of the USP planet from additional longperiod Doppler signals.
With an orbital period of only 7.7 hr, the USP planet GJ 367 b is suitable to the FCO method (see, e.g., Gandolfi et al. 2017;Barragán et al. 2018).GJ 367 is accessible for up to 8 hr at an airmass < 1.5 (i.e, altitude > 40 • ) from La Silla Observatory, allowing one to cover one full orbit in one single night by acquiring multiple HARPS spectra per night.Within the nightly visibility window of GJ 367, the phase of the long-term signals does not change significantly, the variation being 0.029, 0.010, and 0.006 for the 11.5, 34, and 52 day signals, respectively.
We simultaneously modeled the TESS transit light curves and HARPS RV measurements using the open Frequency (d 1 ) Power 0.32d 11.51d 34.0d  source software suite pyaneti (Barragán et al. 2019;Barragán et al. 2022).The code utilizes a Bayesian approach in combination with MCMC sampling to infer the parameters of planetary systems.The photometric data included in the analysis are subsets of the TESS light curve.We selected 2.5 hr of TESS data points centered around each transit and detrended each light curve segment by fitting a second-order polynomial to the out-of-transit data.Following Hatzes et al. (2011), we divided the HARPS RVs into subsets ("chunks") of nightly measurements and analyzed only those chunks containing at least two RVs per observing night, leading to a total of 96 chunks.
The RV model includes one Keplerian orbit for the transiting planet GJ 367 b and 96 nightly offsets.We fitted for a non zero eccentricity adopting the parameterization proposed by Anderson et al. (2010) for the eccentricity e and the argument of periastron of the stellar orbit ω ⋆ (i.e., √ e sin ω ⋆ and √ e cos ω ⋆ ).We fitted for a photometric and an RV jitter term to account for any instrumental noise not included in the nominal TESS and HARPS uncertainties.We used the limb-darkened quadratic model by Mandel & Agol (2002) for the transit light curve.We adopted Gaussian priors for the limb darkening coefficients, using the values derived by Claret (2017) for the TESS passband, and we imposed conservative error bars of 0.1 on both the linear and the quadratic limb-darkening term.As the shallow transit light curve of GJ 367 b poorly constrains the scaled semimajor axis (a b /R ⋆ ), we sampled for the stellar density ρ ⋆ using a Gaussian prior on the star's mass and radius as derived in Sect.3, and recovered the scaled semimajor axis of the planet using the orbital period and Kepler's third law of planetary motion (see, e.g.Winn 2010).We adopted uniform priors over a wide range for all of the remaining model parameters.We ran 500 independent Markov chains.The posterior distributions were created using the last 5000 iterations of the converged chains with a thin factor of 10, leading to a distribution of 250,000 data points for each model parameter.The chains were initialized at random values within the priors ranges.This ensured a homogeneous sampling of the parameter space.We followed the same procedure and convergence test as described in Barragán et al. (2019).The final estimates and their 1σ uncertainties were taken as the median and 68% of the credible interval of the posterior distributions.Table 2 reports prior ranges and posterior values of the fitted and derived system parameters.Figure 5 displays the phase-folded RV curve with our HARPS data, along with the bestfitting Keplerian model.Different colors refer to different nights.Figure 4 shows the phase-folded transit light curve of GJ 367 b, along with the TESS data and bestfitting transit model.We found an RV semi-amplitude variation of K b = 1.003 ± 0.078 m s −1 , which translates into a planetary mass of M b = 0.633 ± 0.050 M ⊕ (7.9 % precision).The depth of the transit light curve implies a radius of R b = 0.699 ± 0.024 R ⊕ (3.4% precision) for GJ 367 b.When combined together, the planetary mass and radius yield a mean density of ρ b = 10.2 ± 1.3 g cm −3 (12.7% precision).The eccentricity of the USP planet (e b = 0.06 +0.07 −0.04 ) is consistent with zero, as expected given the short tidal evolution time-scale and the age of the system.Assuming a circular orbit, our fit gives an RV semi-amplitude of K b = 1.001 ± 0.077 m s −1 , in excellent agreement with the values listed in Table 2.
We note that 27 of the HARPS RV measurements were taken during transits of GJ 367 b.Assuming the star is seen equator-on (i ⋆ = 90°), its rotation period of P rot ≈ 51−54 days (Sect.3.2) implies an equatorial rotational velocity of v rot ≈ 0.45 km s −1 .Using the equations in Triaud (2018), we estimated the semi-amplitude of the Rossiter-McLaughlin effect to be ≈ 0.05 m s −1 , which is too small to cause any detectable effect given our RV uncertainties.

Sinusoidal activity signal modeling
Using the FCO method, we cannot determine the semi-amplitude of the Doppler signals induced by the two outer planets and by stellar activity (Sect.4).In the analysis described in this section, we treated the RV signals associated with stellar activity as coherent sinusoidal signals.Once again, we used the code pyaneti and performed an MCMC analysis similar to the one described in Sect.5.1.The RV model includes three Keplerians, to account for the Doppler reflex motion induced by the three planets GJ 367 b, c, and d, and two additional sine functions, to account for the activity-induced signals at the rotation period of the star (∼52 days) and at the evolution timescale of active regions (∼115 days), as described in Sect. 4. We used Gaussian priors for the orbital period and time of first transit of GJ 367 b as derived in Sect.5.1, and uniform wide priors for all of the remaining parameters.We fitted for a RV jitter term, as well as for a non zero eccentricity both for the USP planet, and for the two outer companions, following the e-ω ⋆ parameterization proposed by Anderson et al. (2010).Details of the fitted parameters and prior ranges are given in Table 3.We used 500 independent Markov chains initialized randomly inside the prior ranges.Once all chains converged, we used the last 5000 iterations and saved the chain states every 10 iterations.This approach generated a posterior distribution of 250,000 points for each fitted parameter.
The RV semi-amplitude variations induced by the three planets are K b = 1.10 ± 0.14 m s −1 , K c = 2.01 ± 0.15 m s −1 , and K d = 1.98 ± 0.15 m s −1 , which imply planetary masses and minimum masses of M b = 0.699 ± 0.083 M ⊕ , M c sin i c = 4.08 ± 0.30 M ⊕ , and M d sin i d = 5.93 ± 0.45 M ⊕ for GJ 367 b, c, and d, respectively, whereas the RV semi-amplitudes induced by stellar activity signals at 51.3 and 138 days are of K ⋆,Rot = 2.52 ± 0.13 m s −1 and K ⋆,Evol = 1.25 ± 0.69 m s −1 .The RV time series along with the best-fitting model are shown in Fig. 6.The phase-folded RV curves for each signal are displayed in Fig. 7.

Multi-dimensional Gaussian process approach
We also followed a multidimensional Gaussian process (GP) approach to account for the stellar signals in our RV time series (Rajpaul et al. 2015).This approach models the RVs along with time series of activity indicators assuming that the same GP, a function G(t), can describe them both.The function G(t) represents the projected area of the visible stellar disk that is covered by active regions at a given time.For our best GP analysis, we selected the activity indicator that shows the strongest signal in the periodograms, i.e., the DLW, and modeled the RVs alongside this activity index.We created a two-dimensional GP model via (1) The amplitudes A RV , B RV , and A DLW are free parameters, which relate the individual time series to G(t).The RV data are modeled as a function of G(t) and its time derivative Ġ(t), since they depend both on the fraction of the stellar surface covered by active regions, and on how these regions evolve and migrate on the disk.The DLW, which measures the width of the spectral lines, has been proven to be a good tracer of the fraction of the surface covered by active regions, and is thus expected to be solely proportional to G(t) (Zicher et al. 2022).For our covariance matrix, we used a quasi-periodic kernel and its derivatives (Barragán et al. 2022;Rajpaul et al. 2015).The parameter P GP in Eq. 3 is the characteristic period of the GP, which is interpreted as the stellar rotation period; λ p is the inverse of the harmonic complexity, which is associated with the distribution of active regions on the stellar surface (Aigrain et al. 2015); λ e is the long-term evolution timescale, i.e., the active region lifetime on the stellar surface.We performed the multidimensional GP regression using pyaneti (described in Barragán et al. 2022), adding three Keplerians to account for the Doppler reflex motion of the three planets, as described in Sect.5.2.We performed an MCMC analysis setting informative Gaussian priors based on the transit ephemeris of the innermost planet GJ 367 b, as derived in Sect.5.1, and uniform priors for all the remaining parameters.We also used uniform priors to sample for the multidimensional GP hyper-parameters.We included a jitter term to the diagonal of the covariance for each time series.
We performed our fit with 500 Markov chains to sample the parameter space.The posterior distributions were created using the last 5000 iterations of the converged chains with a thin factor of 10, leading to a distri-bution of 250,000 data points for each fitted parameter.The chains were initialized at random values within the priors ranges.This ensured a homogeneous sampling of the parameter space.
Priors and results are listed are listed in Table 4. Planets GJ 367 b, c, and d are significantly detected in the HARPS RV time series with Doppler semi-amplitudes of K b = 0.86 ± 0.15 m s −1 , K c = 1.99 ± 0.17 m s −1 , and K d = 2.03 ± 0.16 m s −1 , respectively.
These imply planetary masses and minimum masses of M b = 0.546 ± 0.093 M ⊕ , M c sin i c = 4.13 ± 0.36 M ⊕ , and M d sin i d = 6.03 ± 0.49 M ⊕ .The resulting GP hyperparameters are P GP = 53.67 +0.65  −0.53 days, λ p = 0.44 ± 0.05, and λ e = 114 ± 19 days.The characteristic period P GP is in agreement with the stellar rotation period, as discussed in Sect.3.2, as well as the long-term evolution timescale λ e , which is in agreement with the long-period signal found in the analyses of Sects. 4 and 5.2.
Figure 8 shows the RV and DLW time series, along with the inferred models, whereas Fig. 9 displays the phase-folded RV curves of the three planets and the bestfitting models.

DISCUSSION
The three techniques used to determine the mass of GJ 367 b give results that are consistent to within ∼1σ.We adopted the results from the FCO method, which gives a planetary mass of M b = 0.633 ± 0.050 M ⊕ (7.9% precision).Our result differs by about ∼1σ from the mass of M b = 0.546 ± 0.078 M ⊕ reported by Lam et al. (2021), which was also derived using the FCO method applied on 20 HARPS data chunks that do not entirely cover the orbital phase of the USP planet, as opposed to our 96 chunks.We found that GJ 367 b has a radius of R b = 0.699 ± 0.024 R ⊕ (3.4% precision), consistent with the value of R b = 0.718 ± 0.054 R ⊕ from Lam et al. ( 2021), but more precise, thanks to the additional TESS photometry and increased cadence.Lam et al. (2021) reported a density of ρ b = 8.1 ± 2.2 g cm −3 .The higher mass but similar radius measured in this work make the USP planet denser, with a ultra-high density of ρ b = 10.2 ± 1.3 g cm −3 .
GJ 367 b belongs to the handful of small USP planets (R p < 2 R ⊕ , M p < 10 M ⊕ ) whose masses and radii are known with a precision better than 20%.Figure 10 shows the mass-radius diagram for small USP planets along with the theoretical composition models for rocky worlds (Zeng et al. 2016).GJ 367 b is the smallest and densest USP planet known to date.The position of the planet on the mass-radius diagram sug-  The other configuration for the outer planets to transit is if they are mutually inclined with planet b (i c ̸ = i b ; i d ̸ = i b ).This is not an unusual architecture for systems with USP planets (Dai et al. 2018).With minimum masses of 4.13 M ⊕ and 6.03 M ⊕ , GJ 367 c and d are expected to have a minimum radius of ∼1.5 R ⊕ given the maximum collisional stripping limit (Marcus et al. 2009(Marcus et al. , 2010)).The inclination of their orbits should be larger than ∼86.5°to produce non-grazing transits.If the two planets had radii of 1.5 R ⊕ , the transit depths would be ∼900 ppm.For 2.5 R ⊕ , the non grazing transit depth is expected to be ∼2500 ppm.Similarly, transit depth of a 4 R ⊕ planet is ∼ 6400 ppm.The rms of the TESS light curve of GJ 367 is approximately 500 ppm.The expected transit times of planets c and d also fall well within the baseline of the TESS data.Therefore, GJ 367 c and d would be easily detected if they produced non-grazing transits.
The presence of two additional planetary companions to GJ 367 b is in line with the tendency of USP planets to belong to multi-planet systems (Dai et al. 2021).While the origin of USP planets is debated, it is likely that the presence of outer companions could be responsible for the migration processes that carried USP planets to their current positions (Pu & Lai 2019;Petrovich et al. 2019), which would lie inside the magnetospheric cavity of a typical protoplanetary disk.In these models, planet migration occurs after the dissipation of the protoplanetary disk, by a combination of eccentricity forcing from the planetary companions and tidal dissipation in the innermost planet's interior, although the precise dynamical forcing mechanisms are debated.Serrano et al. (2022) showed, for TOI-500, that the loweccentricity secular forcing proposed by Pu & Lai (2019) can explain the migration of that system's USP planet from a formation radius of ∼0.02 au to its observed location.For GJ 367 b, we tested this migration scenario with similar initial conditions (USP planet starting at 0.02 au, initial eccentricities of all planets set to 0.2), and found only modest migration of planet b, from 0.02 au to ∼0.01 au, short of the observed 0.007 au.This would support an alternative migration history, such as by high-eccentricity secular chaos (Petrovich et al. 2019).However, we note that GJ 367 c and GJ 367 d lie close to a 3:1 mean motion commensurability, having a period ratio of 2.98, and that the dynamics may be affected by the 3:1 resonance, which can provide additional eccentricity forcing in the system.A deeper analysis would be needed to draw definite conclusions on the formation, migration, and dynamics of the system.The GJ 367 system is thus an excellent target for studying planetary system formation and evolution scenarios.

Internal structure of GJ 367 b
Given the high precision of both the derived mass and radius, we used a Bayesian analysis to infer the internal structure of GJ 367 b.We followed the method described in Leleu et al. (2021), which is based on Dorn et al. (2017).Prior to presenting the results of our modeling, we provide the reader with a brief overview over the used model.
Modeling the interior of an exoplanet is a degenerate problem: there is a multitude of different compositions that could explain the observed planetary density.This is why a Bayesian modeling approach is used, with the goal of finding posterior distributions for the internal structure parameters.We assumed a planet that is made up of four fully distinct layers: an inner iron core, a silicate mantle, a water layer, and a gas layer made up of hydrogen and helium.In our forward model, this atmospheric layer is modeled separately from the rest of the planet following Lopez & Fortney (2014) and it is assumed that it does not influence the inner layers.While the presence of a gas and water layer is not expected given the high equilibrium temperature of GJ 367 b, we still included them in the initial model setup, as this is the most general way to model the planet and give us all possible compositions that could lead to the observed planetary density.As input parameters, we used the planetary and stellar parameters listed in Tables 1  and 2, i.e. the transit depth, RV semi-amplitude and period of the planet and the mass, radius, effective temperature, and metallicity of the star.In addition, we chose a prior distribution of 5 ± 5 Gyr for the age of the star.
For our Bayesian analysis, we chose a prior that is uniform in log for the gas mass fraction of the planet.For the mass fractions of the inner iron core, the silicate mantle, and the water layer (all with respect to the solid part of the planet without the H/He gas layer), our chosen prior is a distribution that is uniform on the simplex on which these three mass fractions add up to one.Additionally, we added an upper limit of 0.5 for the water mass fraction in the planet, as a planet with a pure water composition is not physical (Thiabaud et al. 2014;Marboeuf et al. 2014).Note that choosing very different priors influences the results of the internal structure analysis.
There are multiple studies that show a correlation between the composition of planets and their host star (e.g.Thiabaud et al. 2015;Adibekyan et al. 2021), but it is not clear if this correlation is 1:1 or has a different slope.As a first step, we assumed the composition of the planet to match that of its host star exactly.  .Mass-radius diagram of small USP planets (P < 1 days, R < 2 R⊕, M < 10 M⊕) with masses and radii known with a precision better than 20%, as retrieved from the Transiting Extrasolar Planet Catalogue (Southworth 2011).From bottom to top, the solid and dashed curves are theoretical models from Zeng et al. (2016).We highlighted GJ 367 b with a red dot and the previous position on the diagram with a black dot (Lam et al. 2021).
stellar parameters from [Si/H] = [Mg/H] = 0 +1 −1 .However, our analysis showed that the observed mass and radius values of GJ 367 b are not compatible with a 1:1 relationship between the Si/Mg/Fe ratios of the planet and of these sampled synthetic host stars, as it is not possible to reach such a high density under these constraints.The same is true when assuming the steeper slope between stellar and planetary Fe/(Si+Mg) ratios found by Adibekyan et al. (2021).We then repeated our analysis allowing the iron to silicon and iron to magnesium ratios in the planet to be up to a factor 1000 higher than in the sampled stars.
The results of this second analysis are summarized in Figure 11.Indeed, we can see that the Fe/Si ratio in the planet (and therefore also the Fe/Mg ratio) was a factor of 10 to the power of 2.01 +1.10 −0.83 higher than in the sampled host star.We further expect GJ 367 b to host no H/He envelope and no significant water layer.Conversely, we expect the mass fraction of the inner iron core with 0.91 +0.07 −0.23 (median with 5th and 95th percentiles) to be very high.

Formation and evolution of the ultra-high density sub-Earth GJ 367 b
It is not clear how a low-mass high-density planet like GJ 367 b forms.Possible pathways may include the formation out of material significantly more iron rich than thought to be normally present in protoplanetry disks.Although it is not clear if disks with such a large relative iron content specifically near the inner edge (where most of the material might be obtained from) exist (Dullemond & Monnier 2010; Aguichine et al. 2020;Adibekyan et al. 2021) Another possibility is the formation of a larger planet with a lower bulk density.Subsequently the planet differentiates into a denser core and less dense outer layers.These outer layers are then removed.This may be accomplished through two processes.
(i) A first process is collisional stripping.The preferred removal of outer material in giant collisions (e.g., Marcus et al. 2009;Reinhardt et al. 2022) might have increased the bulk density.Indeed, this is what might have lead to the high iron fraction and therefore high density of Mercury for its small size (Benz et al. 1988(Benz et al. , 2007)), as a Mercury with a chondritic composition would have a lower density.The amount of maximum stripping is governed by the mass, impact velocity, and impact parameter of the impactor (Marcus et al. 2010;Leinhardt & Stewart 2012).Preferable removal of outer layers requires the right mass ratio (close to unity), right impact parameter (close to grazing), and right relative velocity.There is also evidence that, at least in some systems, densities have been altered by collisions (Kepler-107; Bonomo et al. 2019) A problem to effectively remove mass might be the removal of collision debris and the avoidance of a re-accretion of debris material onto the planet, as re-accretion would leave the bulk density largely unchanged.However this might not be such a large problem as originally thought (Spalding & Adams 2020).Our measurement of the bulk density of GJ 367 b suggests that collisional stripping has to be remarkably effective in removing non-iron material from the planet if it is the only process at work.
(ii) A second process that might have played a role in shaping GJ 367 b after core formation is removal of outer layers of material facilitated by the enormous stellar radiation to which this planet is subjected.At the equilibrium temperature of 1365 ± 32 K, material might sublimate, be uplifted, and transported away from the surface.
Of course, all of the above discussed processes could have contributed to create the nearly pure ball of iron, known as GJ 367 b.

CONCLUSIONS
We report refined mass and radius determinations for GJ 367 b, the ultra-high density, USP sub-Earth planet transiting the M-dwarf star GJ 367 recently discovered by Lam et al. (2021).We used new TESS observations from sectors 35 and 36 combined with 371 Doppler measurements collected as part of an intensive RV followup campaign with the HARPS spectrograph.We derived a precise planetary mass of M b = 0.633 ± 0.050 M ⊕ (7.9 % precision) and a radius of R b = 0.699 ± 0.024 R ⊕ (3.4% precision), resulting in an ultra-high density of ρ b = 10.2 ± 1.3 g cm −3 (∼13%).According to our internal structure analysis, GJ 367 b is predominantly composed of iron with an iron core mass fraction of 0.91 +0.07 −0.23 , which accounts for the aforementioned planetary density.In addition, our HARPS RV follow-up observations, which span a period of nearly ∼3 yr, allowed us to discover two additional non-transiting small companions with orbital periods of ∼11.5 and 34.4 days, and minimum masses of M c sin i c = 4.13 ± 0.36 M ⊕ and M d sin i d = 6.03 ± 0.49 M ⊕ .GJ 367 joins the small group of well-characterized multiplanetary systems hosting a USP planet, with the inner planet GJ 367 b being the densest and smallest USP planet known to date.This unique multiplanetary system hosting this ultra-high density, USP sub-Earth is an extraordinary target to further investigate the formation and migration scenarios of USP systems.Note-The entire RV data set is available in its entirety in machine-readable form.Only a portion of this table is shown here to demonstrate its form and content.
also measured a Ca II H & K chromospheric activity index of log R ′ HK = −5.214± 0.074 from their 77 HARPS spectra.Based on the log R ′ HK −rotation empirical relationship for M-dwarfs from Astudillo-Defru et al.

Figure 1 .
Figure 1.Least-squares deconvolved (LSD) Stokes I, V and null profiles of GJ 367.The polarization profiles are shifted vertically and expanded by a factor of 100 relative to the intensity profile.The vertical dashed lines indicate velocity interval adopted for the longitudinal field measurement.

Figure 2 .
Figure 2.Co-added HARPS spectrum of AU Mic (black line) in the spectral region encompassing the Li i 6708 Å absorption line.Superimposed with a thick red line is the co-added HARPS spectrum of GJ 367, which has been rotationally broadened to match the v sin i⋆ = 7.8 km s −1 of AU Mic.

Figure 3 .
Figure3.Generalized Lomb-Scargle (GLS) periodograms of the HARPS RV measurements (upper panel ); RV residuals after subtracting the f1 signal at 11.5 days (second panel ), the f2 signal related to stellar activity at 52.2 days (third panel ), the f3 signal at 34 days (fourth panel ), the f4 and f5=f b signals at 115 days and 0.322 days (fifth panel ).Also shown are periodograms of the activity indicators (remaining panels).The 10%, 1% and 0.1% false alarm probabilities (FAPs) estimated using the bootstrap method are shown with horizontal blue lines.The red vertical lines mark the orbital frequencies of the transiting planet GJ 367 b (f5 = f b = 3.106 day −1 ), and of the additional Doppler signals we found in the HARPS data, which are associated to the presence of two additional orbiting planets (f1 = fc = 0.086 day −1 and f3 = f d = 0.029 day −1 ).The shaded yellow bands indicate the rotation period of the star centered around f2, and the long-period stellar signal f4.

Figure 4 .
Figure 4. TESS transit light curve of GJ 367 b and bestfitting model folded at the orbital period of the planet.

Figure 5 .
Figure 5. HARPS RVs of GJ 367 phase-folded at the orbital period of the USP planet and best-fitting model as derived using the FCO method.The different colors refer to the 96 different nightly chunks, which include at least two measurements per night.

−0. 22 √Figure 6 .
Figure 6.HARPS RV time series of GJ 367 along with the best-fitting five-signal model (three planets + stellar rotation + long-period stellar signal).Data are shown as blue filled circles with their nominal uncertainties.The vertical gray lines mark the error bars including the RV jitter.

Figure 7 .
Figure 7. Phase-folded RV curves of GJ 367 b (a), GJ 367 c (b), GJ 367 d (c), stellar rotation (d), and long-period stellar signal (e).Data are shown as blue filled circles with their nominal uncertainties.The vertical gray lines mark the error bars including the RV jitter.
[m s −1 ] J [0,100] 1.59 ± 0.07 Note-U[a, b] refers to uniform priors between a and b; N [a, b] refers to Gaussian priors with mean a and standard deviation b; J [a, b] refers to Jeffrey's priors between a and b.Inferred parameters and uncertainties are defined as the median and the 68.3% credible interval of their posterior distributions. (* ) Assuming an orbital inclination of i b = 79.89+0.87 −0.85 °, from the modeling of TESS transit light curves (Sect.5.1).

Table 4 .Figure 8 .
Figure 8. RV and differential line width (DLW) time series with best-fitting models from the multi-GP (solid black lines) and the 1σ and 2σ credible intervals of the corresponding GP models (light gray shaded areas).The upper panel shows the RV data with the full model in black, whilst the planetary signal (blue), and stellar (red) inferred models are shown with a vertical offset for clarity.The lower panel shows the DLW along with the stellar inferred model.Data are shown with filled circles, with their nominal error bars and semitransparent error bar extensions accounting for the inferred jitter term.

Figure 9 .
Figure 9. Phase-folded RVs curve of GJ 367 b (a), GJ 367 c (b), and GJ 367 d (c).Data are shown as filled green circles with the error bars and the light-green error bars accounting for the inferred RV jitter term.gests that its composition is dominated by iron.Taking into account its mean density, GJ 367 b leads the class of super-Mercuries, namely, extremely dense planets containing an excess of iron, analogous to Mercury: K2-229 b (Santerne et al. 2018), 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), HD 137496 b (Azevedo Silva et al. 2022), HD 23472 b (Barros et al. 2022), and TOI-1075 b (Essack et al. 2022).The two non-transiting planets GJ 367 c and GJ 367 d have orbital periods of ∼11.5 d and 34.4 days, respectively, and minimum masses of M c sin i c = 4.13 ± 0.36 M ⊕ and M d sin i d = 6.03 ± 0.49 M ⊕ , as derived adopting the multidimensional GP approach to model stellar activity (Sect.5.3).We note that our minimum mass determinations are in very good agreement with those of M c sin i c = 4.08 ± 0.30 M ⊕ and M d sin i d = 5.93 ± 0.45 M ⊕ we determined modeling stellar activity with two sinusoidal components (Sect.5.2).

Figure 10
Figure10.Mass-radius diagram of small USP planets (P < 1 days, R < 2 R⊕, M < 10 M⊕) with masses and radii known with a precision better than 20%, as retrieved from the Transiting Extrasolar Planet Catalogue(Southworth 2011).From bottom to top, the solid and dashed curves are theoretical models fromZeng et al. (2016).We highlighted GJ 367 b with a red dot and the previous position on the diagram with a black dot(Lam et al. 2021).

Figure 11 .
Figure 11.Results of a Bayesian analysis of the internal structure of GJ 367 b.The depicted internal structure parameters are: the mass fractions of the inner iron core and a possible water layer with respect to the solid planet without a possible H/He layer), the molar fractions of Si and Mg within the silicate mantle layer, the molar fraction of iron in the inner core, the logarithm with base 10 of the mass of the H/He layer in Earth masses, and the shift in the Fe/Si ratio of the planet with respect to the Fe/Si ratio in the host star, again in a logarithmic scale.The titles in each column show the median of the posterior distributions and the 5th and 95th percentiles.

Table 2 .
GJ 367 b parameters from the joint FCO and transit modeling with pyaneti.
Since we do not have measured values for the stellar [Si/H] and [Mg/H], we left them unconstrained and sampled the

Table 1 .
Radial velocities and spectral activity indicators.