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^*

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 20 s 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, 2014; Smith et al. 2012).

We retrieved TESS Sector 9, 35, and 36 data from from the Mikulski Archive for Space Telescopes (MAST) ³⁷ 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 multitransiting planet search.

GJ 367 was observed with the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph (Mayor et al. 2003), mounted at the ESO-3.6 m telescope of La Silla Observatory in Chile. We collected 295 high-resolution (R ≈ 115,000, λ ∈378–691 nm) spectra between 2020 November 9 and 2022 April 18 (UT), as part of our large observing program 106.21TJ.001 (PI: Gandolfi) to follow-up TESS transiting planets. When added to the 77 HARPS spectra published in Lam et al. (2021), our data includes 371 HARPS spectra.

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, 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 ${{\prime} }_{\mathrm{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 A1 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).

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 subexposures 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.

We derived the spectroscopic parameters using the new coadded 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 $\mathrm{log}{g}_{\star }$, 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 errors—which 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.040}^{+0.035}$ M_⊙ and radius of R_⋆ = 0.411${}_{-0.034}^{+0.032}$ 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.

Table 1. Fundamental Parameters of GJ 367

Parame2ter	Value	Reference
Name	GJ 367
	TOI-731
	TIC 34068865
R.A. (J2000)	09:44:29.15	[1]
Decl. (J2000)	−45:46:44.46	[1]
TESS-band magnitude	8.032 ± 0.007	[2]
V-band magnitude	10.153 ± 0.044	[3]
Parallax (mas)	106.173 ± 0.014	[1]
Distance (pc)	9.413 ± 0.003	[1]
Star mass M* (M⊙)	0.455 ± 0.011	[4]
Star radius R* (R⊙)	0.458 ± 0.013	[4]
Effective temperature Teff (K)	3522 ± 70	[4]
Stellar density ρ* (ρ⊙)	${4.75}_{-0.39}^{+0.44}$	[4]
Metallicity [Fe/H]	−0.01 ± 0.12	[4]
Surface gravity $\mathrm{log}\,{g}_{\star }$	4.776 ± 0.026	[4]
Luminosity L* (L⊙)	${0.0289}_{-0.0027}^{+0.0029}$	[4]
log R ${{\prime} }_{\mathrm{HK}}$	−5.169 ± 0.068	[4]
Spectral type	M1.0 V	[5]

Note. [1] Gaia Collaboration et al. (2021), [2] TESS input catalog (TIC; Stassun et al. 2018, 2019), [3] Paegert et al. (2022), [4] This work, [5] Koen et al. (2010).

Download table as:  ASCIITypeset image

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) also measured a Ca ii H & K chromospheric activity index of log ${R}_{\mathrm{HK}}^{{\prime} }$= −5.214 ± 0.074 from their 77 HARPS spectra. Based on the log ${R}_{\mathrm{HK}}^{{\prime} }$–rotation empirical relationship for M dwarfs from Astudillo-Defru et al. (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}_{\mathrm{HK}}^{{\prime} }$= -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 Section 5.2, and with the period of ${P}_{\mathrm{rot}}={53.67}_{-0.53}^{+0.65}$ d derived by our multidimensional Gaussian process (GP) analysis described in Section 5.3.

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 (Section 3.1). We used about 5000 lines deeper than 20% of the continuum for LSD, reaching an S/N of 7250 per 1 km s⁻¹ 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 〈B_z〉 = − 7.3 ± 3.2 G.

Zoom In Zoom Out Reset image size

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.

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}_{-4.6}^{+3.8}$ 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 18 Gyr.

In this respect, the results presented in our study shows compelling evidence that GJ 367 is not a young star:

• 1.  
  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 days, which translates into a gyrochronological age of ∼4.6–4.8 Gyr (Barnes 2010; Barnes & Kim 2010).
• 2.  
  The HARPS spectra of GJ 367 show no significant lithium absorption line at 6708 Å. Figure 2 displays the coadded 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⁻¹). 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).
• 3.  
  The low level of magnetic activity inferred by the Ca H & K indicator log ${R}_{\mathrm{HK}}^{{\prime} }$ (Section 3.2) and the weak magnetic field (Section 3.3) are consistent with an old, inactive M-dwarf scenario (Pace 2013).
• 4.  
  We measured an average Hα equivalent width of EW = 0.0638 ± 0.0014 Å. Using the empirical relation that connects the Hα equivalent width with stellar age (Kiman et al. 2021), this translates into an age ≳300 Myr.

Zoom In Zoom Out Reset image size

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.

We used the floating chunk offset (FCO) method to determine the semiamplitude 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 longer-period 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 long-period 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 source software suite pyaneti (Barragán et al. 2019, 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 nonzero 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., $\sqrt{e}\sin {\omega }_{\star }$ and $\sqrt{e}\cos {\omega }_{\star }$). 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 Section 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 best-fitting 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 best-fitting transit model. We found an RV semiamplitude variation of K_b = 1.003 ± 0.078 m s⁻¹, 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⁻³ (12.7% precision). The eccentricity of the USP planet (${e}_{{\rm{b}}}={0.06}_{-0.04}^{+0.07}$) is consistent with zero, as expected given the short tidal evolution timescale and the age of the system. Assuming a circular orbit, our fit gives an RV semiamplitude of K_b = 1.001 ± 0.077 m s⁻¹, in excellent agreement with the values listed in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. GJ 367 b Parameters from the Joint FCO and Transit Modeling with pyaneti

GJ 367 b	Prior	Derived Value
Model parameters
Orbital period Porb, b [days]	${ \mathcal U }$[0.3219221,0.3219229]	0.3219225 ± 0.0000002
Transit epoch T0,b [BJDTDB − 2,450,000]	${ \mathcal U }$[8544.13235,8544.14035]	8544.13635 ± 0.00040
Planet-to-star radius ratio Rb/R⋆	${ \mathcal U }$[0.001,0.025]	0.01399 ± 0.00028
Impact parameter bb	${ \mathcal U }$[0,1]	${0.584}_{-0.037}^{+0.034}$
$\sqrt{{e}_{{\rm{b}}}}\sin {\omega }_{\star ,{\rm{b}}}$	${ \mathcal U }$[-1.0,1.0]	${0.16}_{-0.22}^{+0.17}$
$\sqrt{{e}_{{\rm{b}}}}\cos {\omega }_{\star ,{\rm{b}}}$	${ \mathcal U }$[-1.0,1.0]	0.04 ± 0.14
Radial velocity semiamplitude variation Kb [m s−1]	${ \mathcal U }$[0,50]	1.003 ± 0.078
Derived parameters
Planet mass Mb [M⊕]	⋯	0.633 ± 0.050
Planet radius Rb [R⊕]	⋯	0.699 ± 0.024
Planet mean density ρb [g cm−3]	⋯	10.2 ± 1.3
Semimajor axis of the planetary orbit ab [au]	⋯	0.00709 ± 0.00027
Orbit eccentricity eb	⋯	${0.06}_{-0.04}^{+0.07}$
Argument of periastron of stellar orbit ω⋆,b [deg]	⋯	${66}_{-108}^{+41}$
Orbit inclination ib [deg]	⋯	${79.89}_{-0.85}^{+0.87}$
Transit duration τ14,b [hr]	⋯	0.629 ± 0.008
Equilibrium temperature Teq,b [K] a	⋯	1365 ± 32
Received irradiance Fb [F⊕]	⋯	${579}_{-52}^{+57}$
Additional model parameters
Stellar density ρ⋆ [g cm−3]	${ \mathcal N }$[6.68,0.59]	6.76 ± 0.59
Parameterized limb-darkening coefficient q1	${ \mathcal N }$[0.3766,0.1000]	0.343 ± 0.095
Parameterized limb-darkening coefficient q2	${ \mathcal N }$[0.1596,0.1000]	${0.163}_{-0.088}^{+0.096}$
Radial velocity jitter term σRV,HARPS [m s−1]	${ \mathcal J }$[0,100]	0.43 ± 0.08
TESS jitter term σTESS	${ \mathcal J }$[0,100]	0.00004 ± 0.00003

Note. ${ \mathcal U }[a,b]$ refers to uniform priors between a and b; ${ \mathcal N }[a,b]$ refers to Gaussian priors with mean a and standard deviation b; ${ \mathcal 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.

^a Assuming zero albedo and uniform redistribution of heat.

Download table as:  ASCIITypeset image

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 (Section 3.2) implies an equatorial rotational velocity of v_rot ≈ 0.45 km s⁻¹. Using the equations in Triaud (2018), we estimated the semiamplitude of the Rossiter-McLaughlin effect to be ≈0.05 m s⁻¹, which is too small to cause any detectable effect given our RV uncertainties.

Using the FCO method, we cannot determine the semiamplitude of the Doppler signals induced by the two outer planets and by stellar activity (Section 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 Section 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 Section 4. We used Gaussian priors for the orbital period and time of first transit of GJ 367 b as derived in Section 5.1, and uniform wide priors for all of the remaining parameters. We fitted for an RV jitter term, as well as for a nonzero 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,000250000 points for each fitted parameter.

Table 3. System Parameters as Derived Modeling the Stellar Signals with Two Sine Functions

Parameter	Prior	Derived Value
GJ 367 b
Model parameters
Orbital period Porb, b [days]	${ \mathcal N }$[0.3219225, 0.0000002]	0.3219225 ± 0.0000002
Transit epoch T0,b [BJDTDB-2,450,000]	${ \mathcal N }$[8544.1364,0.0004]	8544.13632 ± 0.00040
$\sqrt{{e}_{{\rm{b}}}}\sin {\omega }_{\star ,{\rm{b}}}$	${ \mathcal U }$[-1.0,1.0]	−0.23${}_{-0.23}^{+0.30}$
$\sqrt{{e}_{{\rm{b}}}}\cos {\omega }_{\star ,{\rm{b}}}$	${ \mathcal U }$[-1.0,1.0]	−0.07 ± 0.13
Radial velocity semiamplitude variation Kb [m s−1]	${ \mathcal U }$[0.00, 0.05]	1.10 ± 0.14
Derived parameters
Planet mass Mb [M⊕] a	⋯	0.699 ± 0.083
Orbit eccentricity eb	⋯	${0.10}_{-0.07}^{+0.14}$
Argument of periastron of stellar orbit ω⋆,b [deg]	⋯	${251}_{-102}^{+23}$
GJ 367 c
Model parameters
Orbital period Porb, c [days]	${ \mathcal U }$[11.4858,11.5858]	11.543 ± 0.005
Time of inferior conjunction T0,c [BJDTDB-2,450,000]	${ \mathcal U }$[9152.6591,9154.6591]	9153.46 ± 0.21
$\sqrt{{e}_{{\rm{c}}}}\sin {\omega }_{\star ,{\rm{c}}}$	${ \mathcal U }$[-1,1]	0.38 ${}_{-0.13}^{+0.10}$
$\sqrt{{e}_{{\rm{c}}}}\cos {\omega }_{\star ,{\rm{c}}}$	${ \mathcal U }$[-1,1]	0.27 ${}_{-0.14}^{+0.11}$
Radial velocity semiamplitude variation Kc [m s−1]	${ \mathcal U }$[0.00, 0.05]	2.01 ± 0.15
Derived parameters
Planet minimum mass Mc $\sin {i}_{{\rm{c}}}$ [M⊕]	⋯	4.08 ± 0.30
Orbit eccentricity ec	⋯	0.23 ± 0.07
Argument of periastron of stellar orbit ω⋆,c [deg]	⋯	55 ± 18
GJ 367 d
Model parameters
Orbital period Porb, d [days]	${ \mathcal U }$[34.0016,34.6016]	34.39 ± 0.06
Time of inferior conjunction T0,d [BJDTDB-2,450,000]	${ \mathcal U }$[9179.2710,9183.2710]	9180.90 ${}_{-0.81}^{+0.70}$
$\sqrt{{e}_{{\rm{d}}}}\cos {\omega }_{\star ,{\rm{d}}}$	${ \mathcal U }$[-1,1]	−0.10${}_{-0.18}^{+0.20}$
$\sqrt{{e}_{{\rm{d}}}}\cos {\omega }_{\star ,{\rm{d}}}$	${ \mathcal U }$[-1,1]	${0.16}_{-0.20}^{+0.16}$
Radial velocity semiamplitude variation Kd [m s−1]	${ \mathcal U }$[0.00, 0.05]	1.98 ± 0.15
Derived parameters
Planet minimum mass Md $\sin {i}_{{\rm{d}}}$ [M⊕]	⋯	5.93 ± 0.45
Orbit eccentricity ed	⋯	${0.08}_{-0.05}^{+0.07}$
Argument of periastron of stellar orbit ω⋆,d [deg]	⋯	${277}_{-242}^{+58}$
Stellar activity-induced RV signal
Rotation period P⋆,Rot [days]	${ \mathcal U }$[50.0903,52.0903]	51.30 ± 0.13
Rotation RV semiamplitude K⋆,Rot [m s−1]	${ \mathcal U }$[0.00, 0.05]	2.52 ± 0.13
Active region evolution period P⋆,Evol [days]	${ \mathcal U }$[103.1797,163.1797]	138 ± 2
Active region evolution RV semiamplitude K⋆,Evol [m s−1]	${ \mathcal U }$[0.00, 0.05]	1.25 ± 0.14
Additional model parameters
Systemic velocity γHARPS [m s−1]	${ \mathcal U }$[47.806,48.025]	47.91674 ± 0.00013
Radial velocity jitter term σRV,HARPS [m s−1]	${ \mathcal J }$[0,100]	1.59 ± 0.07

Note. ${ \mathcal U }[a,b]$ refers to uniform priors between a and b; ${ \mathcal N }[a,b]$ refers to Gaussian priors with mean a and standard deviation b; ${ \mathcal 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.

^a Assuming an orbital inclination of i_b = 79.89${}_{-0.85}^{+0.87}$°, from the modeling of TESS transit light curves (Section 5.1).

Download table as:  ASCIITypeset image

The RV semiamplitude variations induced by the three planets are K_b = 1.10 ± 0.14 m s⁻¹, K_c = 2.01 ± 0.15 m s⁻¹, and K_d = 1.98 ± 0.15 m s⁻¹, which imply planetary masses and minimum masses of M_b = 0.699 ± 0.083 M_⊕, ${M}_{{\rm{c}}}\sin {i}_{{\rm{c}}}$ = 4.08 ± 0.30 M_⊕, and ${M}_{{\rm{d}}}\sin {i}_{{\rm{d}}}$ = 5.93 ± 0.45 M_⊕ for GJ 367 b, c, and d, respectively, whereas the RV semiamplitudes induced by stellar activity signals at 51.3 and 138 days are of K_⋆,Rot = 2.52 ± 0.13 m s⁻¹ and K_⋆,Evol = 1.25 ± 0.69 m s⁻¹. The RV time series along with the best-fitting model are shown in Figure 6. The phase-folded RV curves for each signal are displayed in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also followed a multidimensional 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

$$\begin{eqnarray}&&\mathrm{RV}={A}_{\mathrm{RV}}G(t)+{B}_{\mathrm{RV}}\dot{G}(t),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\mathrm{DLW}={A}_{\mathrm{DLW}}G(t).\end{eqnarray} \tag{ 2 }$$

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 $\dot{G}(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 quasiperiodic kernel

$$\begin{eqnarray}&&\gamma ({t}_{i},{t}_{j})=\exp \left[-\displaystyle \frac{{\sin }^{2}\left[\pi ({t}_{i}-{t}_{j})/{P}_{\mathrm{GP}}\right]}{2{\lambda }_{{\rm{p}}}^{2}}-\displaystyle \frac{{\left({t}_{i}-{t}_{j}\right)}^{2}}{2{\lambda }_{{\rm{e}}}^{2}}\right]\end{eqnarray} \tag{ 3 }$$

and its derivatives (Barragán et al. 2022; Rajpaul et al. 2015). The parameter P_GP in Equation (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 Section 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 Section 5.1, and uniform priors for all of 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 distribution of 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 semiamplitudes of K_b = 0.86 ± 0.15 m s⁻¹, K_c = 1.99 ± 0.17 m s⁻¹, and K_d = 2.03 ± 0.16 m s⁻¹, respectively. These imply planetary masses and minimum masses of M_b = 0.546 ± 0.093 M_⊕, ${M}_{{\rm{c}}}\sin {i}_{{\rm{c}}}$ = 4.13 ± 0.36 M_⊕, and ${M}_{{\rm{d}}}\sin {i}_{{\rm{d}}}$ = 6.03 ± 0.49 M_⊕. The resulting GP hyper-parameters are P_GP = 53.67${}_{-0.53}^{+0.65}$ 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 Section 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 Sections 4 and 5.2.

Table 4. System Parameters as Derived Modeling the Stellar Signals with a GP

Parameter	Prior	Derived Value
GJ 367 b
Model parameters
Orbital period Porb,b [days]	${ \mathcal N }$[0.3219225, 0.0000002]	0.3219225 ± 0.0000002
Transit epoch T0,b [BJDTDB − 2,450,000]	${ \mathcal N }$[8544.1364, 0.0004]	8544.13631 ± 0.00038
$\sqrt{{e}_{b}}\sin {\omega }_{\star ,{\rm{b}}}$	${ \mathcal U }$[−1.0,1.0]	${0.13}_{-0.32}^{+0.29}$
$\sqrt{{e}_{b}}\cos {\omega }_{\star ,{\rm{b}}}$	${ \mathcal U }$[−1.0,1.0]	−0.05 ± 0.21
Radial velocity semiamplitude variation Kb [m s−1]	${ \mathcal U }$[0.00, 0.05]	0.86 ± 0.15
Derived parameters
Planet mass Mb [M⊕] a	⋯	0.546 ± 0.093
Orbit eccentricity eb	⋯	${0.10}_{-0.07}^{+0.13}$
Argument of periastron of stellar orbit ω⋆,b [deg]	⋯	${71}_{-173}^{+60}$
GJ 367 c
Model parameters
Orbital period Porb,c [days]	${ \mathcal U }$[11.4, 11.6]	11.5301 ± 0.0078
Time of inferior conjunction T0,c [BJDTDB − 2,450,000]	${ \mathcal U }$[9153.0, 9155.0]	9153.84 ± 0.30
$\sqrt{{e}_{{\rm{c}}}}\sin {\omega }_{\star ,{\rm{c}}}$	${ \mathcal U }$[−1.0,1.0]	−0.11${}_{-0.20}^{+0.23}$
$\sqrt{{e}_{{\rm{c}}}}\cos {\omega }_{\star ,{\rm{c}}}$	${ \mathcal U }$[−1.0,1.0]	${0.14}_{-0.15}^{+0.19}$
Radial velocity semiamplitude variation Kc [m s−1]	${ \mathcal U }$[0.00, 0.05]	1.99 ± 0.17
Derived parameters
Planet minimum mass ${M}_{{\rm{c}}}\,\sin {i}_{{\rm{c}}}[{M}_{\oplus }]$	⋯	4.13 ± 0.36
Orbit eccentricity ec	⋯	0.09 ± 0.07
Argument of periastron of stellar orbit ω⋆,c [deg]	⋯	−34${}_{-54}^{+74}$
GJ 367 d
Model parameters
Orbital period Porb,d [days]	${ \mathcal U }$[30.0,40.0]	34.369 ± 0.073
Time of inferior conjunction T0,d [BJDTDB − 2,450,000]	${ \mathcal U }$[9173.0, 9180.0]	9181.82 ± 1.10
$\sqrt{{e}_{{\rm{d}}}}\sin {\omega }_{\star ,{\rm{d}}}$	${ \mathcal U }$[−1.0,1.0]	−0.09 ± 0.19
$\sqrt{{e}_{{\rm{d}}}}\cos {\omega }_{\star ,{\rm{d}}}$	${ \mathcal U }$[−1.0,1.0]	−0.30${}_{-0.13}^{+0.20}$
Radial velocity semiamplitude variation Kd [m s−1]	${ \mathcal U }$[0.00, 0.05]	2.03 ± 0.16
Derived parameters
Planet minimum mass ${M}_{{\rm{d}}}\,\sin {i}_{{\rm{d}}}[{M}_{\oplus }]$	⋯	6.03 ± 0.49
Orbit eccentricity ed	⋯	0.14 ± 0.09
Argument of periastron of stellar orbit ω⋆,d [deg]	⋯	−126${}_{-38}^{+287}$
Additional model parameters
Characteristic period of the GP PGP [days]	${ \mathcal U }$[35.0,65.0]	${53.67}_{-0.53}^{+0.65}$
Inverse of the harmonic complexity λp	${ \mathcal U }$[0.1,3.0]	0.44 ± 0.05
Long-term evolution timescale λe [days]	${ \mathcal U }$[1,300]	114 ± 19
Amplitude ARV [km s−1]	${ \mathcal U }$[0.0,0.005]	0.0016 ± 0.0004
Amplitude BRV [km s−1]	${ \mathcal U }$[−0.05,0.05]	${0.009}_{-0.0019}^{+0.0025}$
Amplitude ADLW	${ \mathcal U }$[0.0,30.0]	${6.31}_{-1.02}^{+1.42}$
Systemic velocity γHARPS [km s−1]	${ \mathcal U }$[47.40,48.42]	47.9168 ± 0.0005
Offset DLW [103 m2 s−2 ]	${ \mathcal U }$[−16.90,14.11]	−0.11${}_{-1.98}^{+2.02}$
Radial velocity jitter term σHARPS [m s−1]	${ \mathcal J }$[0,100]	1.83 ± 0.08
DLW jitter term σDLW [m2 s−2]	${ \mathcal J }$[0,1000]	${378}_{-335}^{+145}$

Note. ${ \mathcal U }[a,b]$ refers to uniform priors between a and b; ${ \mathcal N }[a,b]$ refers to Gaussian priors with mean a and standard deviation b; ${ \mathcal 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.

^a Assuming an orbital inclination of i_b = 79.89${}_{-0.85}^{+0.87}$°, from the modeling of TESS transit light curves (Section 5.1).

Download table as:  ASCIITypeset image

Figure 8 shows the RV and DLW time series, along with the inferred models, whereas Figure 9 displays the phase-folded RV curves of the three planets and the best-fitting models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 semiamplitude 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. Since we do not have measured values for the stellar [Si/H] and [Mg/H], we left them unconstrained and sampled the 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}_{-0.83}^{+1.10}$ 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.23}^{+0.07}$ (median with 5th and 95th percentiles) to be very high.

Zoom In Zoom Out Reset image size

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, 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 as large a 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.