© The Authors 2024

1 Introduction

To date, 790 exoplanets have been characterized by combining photometric – transits – and spectroscopic – radial velocity (RV) – observations (Christiansen 2022)¹. Among them, only 20 orbit stars of magnitudes of <8 and most are short-period planets (<30 days). In the context of new and future space missions, such as the JWST (Gardner et al. 2006) and ARIEL (2029; Tinetti et al. 2018), exoplanets orbiting bright stars are priority targets for atmospheric characterization.

Most of the known exoplanets are sub-Neptunes² and Super-Earths; that is, planets with a radius of around 2.0–4.0 R_⊕. These planet populations are not present in our Solar System. However, because they lie in the transition regime between rocky planets and gas giants, they can provide strong constraints on planetformation models (Howard et al. 2010b). To conduct statistical studies of these planets at a population level, we require precise knowledge of the physical properties of individual targets. In this context, our goal is to characterize the candidate sub-Neptune planet HD 73344b.

In this paper, we present analyses of new photometric and spectroscopic data for this candidate planet, which orbits the bright F star HD 73344 (V = 6.9 mag) with a period of ≈15 days. This planet was first discovered by Yu et al. (2018) based on six transits observed in the K2 data. While the detection was challenging due to the high activity level of this early-type star, we confirm the detection of this planet by combining K2 data with new Spitzer and TESS photometric data, as well as a set of SOPHIE and HIRES RV observations. In addition, our RV analyses reveal a new sub-Jupiter-mass planet candidate, which is nontransiting and has a period of ≈66 days. We propose a new observing strategy to identify and overcome the different sources of stellar activity that impact the characterization of planetary systems (see e.g., Dumusque et al. 2011; Aigrain et al. 2004; Sulis et al. 2020; Meunier et al. 2023). In particular, we demonstrate the benefits of tracking the star at high cadence for whole nights in order to characterize its short-timescale stellar variability (p-mode oscillations, granulation, supergranulation), which is of very large RV amplitude.

The paper is structured as follows. In Sec. 2, we present the observations. In Sec. 3, we derive the fundamental parameters of the star and characterize various sources of stellar activity. We study the planetary system around HD 73344 in Sec. 4. In Sec. 5, we discuss the stability of the system and the internal composition of the transiting planet. We conclude in Sec. 6.

2 Observations

In this section, we present the various sets of photometric and spectroscopic observations of HD 73344 we used in this study. The main information is summarized in Table 1.

2.1 Photometry

2.1.1 K2

K2 (Howell et al. 2014) observed along the ecliptic a series of 100-square degree zones, each lasting approximately up to ~80 days. The broad bandpass of K2 was ranging from 420 to 900 nm. K2 observed HD 73344 (EPIC 212178066) during campaign 16 (C16), which ran from December 07, 2017 to February 25, 2018, and during campaign 18 (C18), which ran from May 12, 2018 to July 02, 2018. Observations were taken at a long cadence, with an integration time of 30 min. In this work, we used only the data acquired during C16 since the observations are affected by strong systematics during C18. This dataset contains 6 transits of HD 73344 b, originally identified by Yu et al. (2018).

We detrended the C16 light curve using the software EVEREST³ (Luger et al. 2016, 2018). We started by masking the transit events with a window taken as twice the transit duration (T_dur ~ 3.3 h). We then corrected the light curve with a single cotrending basis vector (CBV). We obtained the CBV-corrected detrended flux from which we first removed the 3σ outliers. We then used a second sigma clipping step to remove the remaining outliers that look like “flares” in the dataset. For this second step, we used a median filter of 5-h to smooth the light curve and identified the data points at 1σ above this smoothed light curve. The resulting light curve is shown in Fig. 1 (top), and the individual transits in Appendix A.

2.1.2 TESS

The Transiting Exoplanets Survey Satellite (TESS; Ricker et al. 2015) observed HD 73344 (TIC 175193677) in the red-optical bandpass (600–1100 nm) during sector 45 (November 2021-December 2021) and sector 46 (December 2021-January 2022). The two sectors contain three transits of the planet, which has been identified as TOI 5140.01⁴. The short cadence of these observations (120 s) allows a detailed characterization of the transits. In this work, we used the Simple Aperture Photometry (SAP) flux released by the TESS team on MAST⁵. The resulting light curve (normalized by the median flux) is shown in Fig. 1 (bottom), and the individual transits in Appendix A.

2.1.3 Spitzer

We also observed one transit of HD 73344b with the Spitzer space telescope as part of program 14292 (PI: I. Crossfield; Crossfield et al. 2019). On August 9, 2019, we obtained 3600 × 64 0.1 s subarray frames of HD 73344 with the IRAC2 4.5 μm channel (Fazio et al. 2004), spanning 8.5 h and covering one transit of planet b. The raw and calibrated Spitzer data products are available at the Spitzer Heritage Archive; the analysis is presented in Sec. 4.4.

Fig. 1 Light curves of HD 73344. The transit mid-times of planet b in K2 (top) and TESS (bottom) observations are shown with the dotted vertical lines.

2.2 High-resolution spectroscopy

We carried out the RV follow-up observations of HD 73344 with SOPHIE and HIRES spectrographs over a total time span covering ~715 days. The SOPHIE and HIRES RV are analyzed to get the mass of the transiting planet in Sec. 4.

2.2.1 SOPHIE

We observed HD 73344 with the high-resolution echelle spectrograph SOPHIE (Perruchot et al. 2008) at the Haute-Provence Observatory (OHP, France) as part of the program dedicated to the RV follow-up of K2 planet candidates⁶. The target was observed between 2018-11-02 and 2020-03-01, gathering 345 high-resolution spectra.

The observations were carried out using SOPHIE high resolution (HR) mode (resolving power of λ/Δλ ≈ 75 000 at 550 nm), with simultaneous Fabry-Perot (FP) calibration lamp measurements. The latter enabled us to monitor instrumental drift, ensuring precise and accurate RV measurements. The exposure time was set at 900 s with the classic observational strategy of 3 points per night to average the stellar variability, resulting in a median signal-to-noise ratio (S/N; measured on each points) of 149 per pixel at 550 nm.

Radial velocity calculations were performed using the SOPHIE data reduction system (DRS, Bouchy et al. 2009), employing a G2 mask to extract RVs. To enhance the accuracy of SOPHIE measurements, we implemented the optimized procedure outlined in Heidari (2022) and Heidari et al. (2024). This procedure in particular encompasses: (1) CCD charge transfer inefficiency correction (Bouchy et al. 2013); (2) atmospheric dispersion correction (Modigliani et al. 2019); and (3) RV master constant correction to correct long-term instrumental drifts (Courcol et al. 2015). In addition to the RV observations, using the DRS we also calculated some useful spectroscopic activity indicators such as the Full Width at Half Maximum (FWHM) and the bisector inverse slope (BIS, Queloz et al. 2001). We then calculated the $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ following Noyes et al. (1984) and Boisse et al. (2010), and the Hα index following Boisse et al. (2011). From the raw RV, we removed the 3σ outliers, and the data points with RV uncertainties >5 m s⁻¹ (7 points removed in total). The final RV time series contain 312 data points (hereafter: the “unbinned” dataset), spread over 137 individual nights (used to generate the “binned” dataset). The mean RV uncertainty on all measurements is 2.7 m s⁻¹.

In complement to this long RV campaign, we observed HD 73344 continuously for two consecutive nights to monitor the short timescale stellar variability (dominated by p-mode oscillations, granulation, and supergranulation). The first night (2022-01-11) contains N = 51 data points, taken with an exposure time between τ_exp = [10.4, 16.4] min, during a total of T ~ 9.12 h. The RV shows a significant dispersion, with an RMS of 4.65 m s⁻¹. RV uncertainties on each measurements range from 4.5 to 4.8 m s⁻¹ over the night, and are therefore similar to the observed dispersion. To investigate in more detail this short-term variability, we observed HD 73344 during a second night (2022-01-12) at a shorter temporal cadence (τ_exp = [3.4, 8.3] min, T ~ 9.36 h, N = 152). The RMS of this second RV dataset is 8.72 m s⁻¹, confirming the strong amplitude of the short-term variability. RV uncertainties on each measurements are significantly lower (from 2 to 3.7 m s⁻¹ over the night) compared to the first night due, in particular, to very good atmospheric condition (seeing).

The two sets of observations are shown in Fig. 2. We note a similar pattern across both nights, marked by flux drops at the beginning and end of each night. While these flux drops may indeed stem from instrumental systematics (in particular, we identified a potential issue with the ADC used for observations taken at airmass >1.7), the exact source remains uncertain. When considering only the data points obtained during the middle of the nights (airmass <1.7), we still observe a considerable RV dispersion (exceeding 7 m s⁻¹ for the second night). Similar RV amplitudes are also independently observed in the nightly observations taken with the HIRES/Keck spectrograph. We are therefore confident that the dominant RV variability observed in the high cadence SOPHIE dataset is of stellar origin. The characteristics of this variability are given in Sect. 3.2.2. Based on these two nights of observations, we expect that the classic observing strategy (which consists of observing the target 3 times per night and binning these 3 points) will not be sufficient to significantly reduce the short-term (stellar) variability. This will be confirmed in Sec. 4.2.

Fig. 2 Radial velocity of HD 73344 obtained over two consecutive nights with the SOPHIE spectrograph. Best fitting GP model from their joint analysis is shown in purple. Observations taken at airmass >1.7 are shown in red.

2.2.2 HIRES

We obtained 238 additional RV data points with the HIRES spectrometer (Vogt et al. 1994) installed at the Keck I telescope from 2018-03-17 to 2021-06-03. These observations used the B5 decker, which has a slit width of 0.861″ and gives an effective resolution of 48 000, and HIRES iodine cell. They had typical integration times of 40 s (depending on observing conditions). We followed standard procedures of the California Planet Search for the HIRES observations and reductions (Howard et al. 2010a). The observing strategy was to take 3 sets of 5 consecutive observations over the course of the nights to reduce sensitivity to stellar variations. We grouped these observations into 3 points per night to mimic the sampling of SOPHIE observations. The final, binned time series contains 39 data points, spread over 19 nights. The mean RV uncertainty on all measurements is 1.2 m s⁻¹.

3 Stellar properties

In this section, we first describe how we inferred the fundamental stellar parameters of HD 73344 from the SOPHIE spectra, and the stellar abundances from both the SOPHIE and HIRES spectra. Then, we study the stellar activity signatures in both photometric and spectroscopic data. In particular, we look at the variability modulated with the stellar rotation (spots/faculae), and the variability evolving on short timescales (oscillations, convection).

3.1 Fundamental parameters

The stellar spectroscopic parameters (T_eff, log g, microturbulence, [Fe/H]) were estimated using the ARES+MOOG methodology. The methodology is described in detail in Santos et al. (2013); Sousa (2014); Sousa et al. (2021). To consistently measure the equivalent widths (EW) we used the latest version of ARES⁷ (Sousa et al. 2007, 2015). The list of iron lines is the same as the one presented in Sousa et al. (2008). For this we used the combined SOPHIE spectra: we coadded the spectra until reaching S/N~2000, where each individual spectra was corrected in RV to the rest frame prior to be coadded. To find the ionization and excitation equilibrium in this analysis we used a minimization process to converge for the best set of spectroscopic parameters. This process makes use of a grid of ATLAS model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973). We also derived a more accurate trigonometric surface gravity using recent Gaia data following the same procedure as described in Sousa et al. (2021), which provided a consistent value when compared with the spectroscopic surface gravity. In this last process, we also estimated the stellar mass and radius using the calibrations presented in Torres et al. (2010). Furthermore, we determine the Li abundance of this star by performing spectral synthesis also using the code MOOG and ATLAS atmospheres, as well as the above derived stellar parameters. We obtained a value of A(Li) = 2.81 ± 0.05 dex, which is typical of young stars of this T_eff. From this analysis we can also get an estimate of the inclined rotational velocity, after considering the instrumental broadening given by the SOPHIE spectral resolution (R~75 000) and applying the macrotuburbulence velocity empirical calibration from Doyle et al. (2014) dependent on T_eff and logg (V_mac =4.7 km s⁻¹). The measured projected rotational velocity υ sin i_* is 5.3 km s⁻¹. We report the stellar parameters in Table 2.

In addition, we measured the stellar abundances for multiple chemical elements using both SOPHIE and HIRES spectra. For SOPHIE spectra, using the aforementioned stellar atmospheric parameters (we considered the trigonometric surface gravity), we determined the abundances of refractory elements following the classical curve-of-growth analysis method described in Adibekyan et al. (e.g., 2012, 2015); Delgado Mena et al. (e.g., 2017). Similar to the stellar parameter determination, we used ARES to measure the EWs of the spectral lines of these elements, and used a grid of Kurucz model atmospheres along with the radiative transfer code MOOG to convert the EWs into abundances, assuming local thermodynamic equilibrium. Although the EWs of the spectral lines were automatically measured with ARES, for Mg which has only three lines available we performed careful visual inspection of the EWs measurements. Abundances of the volatile elements, C and O, were derived following the method of Delgado Mena et al. (2021); Bertran de Lis et al. (2015) and using the same code and model atmospheres. All the abundance ratios [X/H] are obtained by doing a differential analysis with respect to a high S/N solar (Vesta) spectrum. The final abundances, shown in Appendix B, are typical of a galactic thin-disk star. Moreover, we used the chemical abundances of some elements to derive ages through the so-called chemical clocks (i.e., certain chemical abundance ratios which have a strong correlation for age). We applied the 3D formulas described in Table 10 of Delgado Mena et al. (2019), which also consider the variation in age produced by the effective temperature and iron abundance. The chemical clocks [Y/Mg], [Y/Zn], [Y/Ti], [Y/Si], [Sr/Zn], [Sr/Ti], [Sr/Mg] and [Sr/Si] were used from which we obtain a weighted average age of 2.0 ± 0.2 Gyr. We note that this small uncertainty reflects the high precision of the different chemical clocks for this specific star and is smaller than the true age uncertainty. For HIRES spectra, we measured the stellar abundances following the approach of Polanski et al. (2022, and in prep.,), using the Cannon (Ness 2018), which was designed to be applied to iodine-free spectra from HIRES on Keck I (Rice & Brewer 2020). KeckSpec was trained using a sample of high-quality (S/N > 100) HIRES spectra for which abundances of 15 chemical elements were determined in Brewer et al. (2016). We used an iodine-free spectrum that reached an S/N per pixel of 214. We calculated the α element enhancement and found [α/Fe] values of ~-0.03 dex making HD 73344 chemically consistent with the thin disk. We also report the stellar abundances in Appendix B.

3.2 Stellar activity signatures

3.2.1 Magnetic activity modulated with rotation period

Both the light curves (Fig. 1) and the RV observations (see Sec. 4.3) show strong activity signatures. To analyze the frequency content of this variability, we first performed the generalized Lomb-Scargle periodogram (GLSP; Baluev 2008; Zechmeister & Kürster 2009) of the K2 photometric data, after masking out planet transits (see Fig. 3). The GLS approach is to fit at each frequency a floating mean (a constant) coupled with a periodic term of unknown phase and amplitude. The definition of the GLSP used here is $$P_{\mathrm{G}\mathrm{L}\mathrm{S}}(\nu ):=\frac{{\chi }_0^2-{\chi }^2(\nu )}{{\chi }_0^2},$$(1)

with ${\chi }_0^2$ the residual sum of squares (RSS) resulting by fitting only a constant, and χ²(ν) the RSS by jointly fitting a constant and a sinusoid at frequency ν (see Eq. (4) of Zechmeister & Kürster 2009). The fit is obtained through a weighted least squares problem, with weights provided by the RV uncertainties.

In Fig. 3, we observe two peaks, the highest of which corresponds to the period of ~8. 39 days, close to the value reported by Yu et al. (2018) for the stellar rotation period (P_rot). The same analyses on the individual and combined TESS sectors return slightly longer periods with the highest periodogram peak appearing at ~9. 54 days (sector 45), ~8. 82 days (sector 46), and ~9. 26 days (both sectors, see Fig. 3). The GLSPs of RV data and chromospheric indicators also show a clear peak at 9.1 ± 0.2 days. Finally, the joint analysis of photometric and RV observations with quasi-periodic Gaussian process models finds P_rot = 9.09 ± 0.04 days (see details in Sec. 4.3 and GLSP of the RV data and indicators in Fig. 5). The latter is the value that we reported in Table 2. Combining P_rot with the υ sin i_* ~ 5.3 km s⁻¹ measured from SOPHIE observations⁸, we find an inclination of the rotation axis of i_* ~ 53°. This suggest that the stellar rotation and the orbit of the transiting planet could be misaligned, which would deserve further investigation for the implication on the system history (Huber et al. 2013).

Over the two SOPHIE campaigns, we observe an increase in magnetic activity signatures, with median $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ values decreasing from −4.65 ± 0.03 to −4.56 ± 0.03 (a variation of 2% over an average duration of 224 days). This star is definitely more active than the Sun, as indicated by its $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ value of approximately −4.9 (Brandenburg et al. 2017). Looking in more detail at the temporal variability of the various spectroscopic indicators in Appendix C, we note that these magnetic features remain significant over 3 to 4 P_rot (i.e., 27–36 days). Although it is difficult to identify a precise stellar origin (spot/faculae) of the periodicities observed in the various activity indicators, it is worth noting that the strongest signals at 3 −4 P_rot occur in the CCF area indicator, rather than in the $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ indicator (see Appendix C). As suggested by Costes et al. (2021), this observation may imply that HD 73344 is dominated by faculae.

We then tracked the photometric signatures of the stellar magnetic structures using the photospheric activity proxy S_ph (Mathur et al. 2014a) to place our target into the F-type star population. This global proxy is defined as the standard deviation calculated over subseries of 5 × P_rot in length (Mathur et al. 2014b). The contribution of photon noise (σ_ϕ) is subtracted from this value. For K2 data, a direct relationship between stellar magnitude and photon noise has been derived in Jenkins et al. (2010): it gives σ_ϕ ~ 5 ppm. For the TESS data (both sectors combined), our calculations⁹ also give a value of σ_ϕ ~ 5 ppm. At the end, we found a mean < S_ph >= 1132 ppm in the K2 data, and < S_ph > = 894 ppm in the TESS data. When compared to the sample of 22 F dwarfs studied by Mathur et al. (2014a), our target appears to be much more active. However, we need to keep in mind that the sample is biased because the stars in their sample have detected solar-like oscillations. It is known that strong magnetic activity can lead to smaller mode amplitudes (e.g., García et al. 2010; Chaplin et al. 2011; Mathur et al. 2019) so the F dwarfs sample is mostly constituted of low magnetic activity stars. The recent catalog of rotation periods measured for more than 55 000 Kepler stars (Santos et al. 2019, 2021) gives a better representation of the magnetic activity level of main-sequence solar-like stars. From the Fig. 7 (third row) of Santos et al. (2021), we can see that F-type stars with rotation periods around 8–9 days can have S_ph values similar to the ones obtained for HD 73344. This makes our target less atypical but it is still among the most active stars of the F-dwarf sample. This can be explained by the star’s young age (see Sects. 3.1 and 4.3).

We note that we do not detect “obvious” bright or dark spot crossing events during the planetary transits to study the evolution of such features on the stellar surface. We are therefore most sensitive to the uncrossed magnetic regions that generate the long-term photometric variability observed in Fig. 1.

Fig. 3 Generalized Lomb-Scargle periodograms of the K2 and TESS photometric data between periods of 5 and 40 days. Transit of planet b have been masked. The green vertical lines indicate the period of 9.1 ± 0.2 days found with RV data and chromospheric indicators.

3.2.2 Stellar short-term variability

At short timescales, both photometric and RV observations exhibit correlated noise, likely originating from convective phenomena (granulation, supergranulation) and stellar oscillations (p-modes). Below, we analyzed periodograms of both K2 and TESS photometric data (transits masked), and short-cadence SOPHIE RVs within the period range < 1 day (Fig. 4).

The K2 periodogram (first panel) shows a typical power increase towards low frequency, indicative of granulation signal (Kallinger et al. 2014), but temporal sampling (30 min) hinders precise characterization. The TESS periodogram (middle panel) exhibits a slight power increase, with less significance at low frequency compared to K2, possibly due to the redder wavelength range of TESS observations, reducing granulation amplitude (because the contrast between the rising and falling cells is reduced; e.g., see similar conclusions in Sulis et al. 2023). A 20 s cadence, known to reduce noise and allow the detection of some stellar p-modes signals in TESS data (Huber et al. 2022), could be explored for stellar granulation detection.

The periodogram of SOPHIE RV data (bottom panel) also reveals increased power at low frequency (periods between 50 min and 6.8 h). Note that this power increase does not change when we remove the RV data that are possibly affected by instrumental systematics (see Sect. 2.2.1). The short temporal cadence of this RV dataset (see Table 1) allows to characterize this short-term correlated noise. We then first model the RV periodogram using classical Harvey functions (Harvey 1988) with two components (Kallinger et al. 2014): a WGN (high-frequency region), and a Lorentzian-like function for the granulation noise. Stellar oscillations, not resolved in our observations, were not modeled¹⁰. Best fitting Harvey functions are shown in Fig. 4 to help the visual inspection.

For a more accurate modeling of this short timescale stellar variability, we then fitted various GPs with different covariance matrices. The model based on a square exponential (SE) covariance function best reproduced (in terms of likelihood) the RV data (taken individually or jointly). The GP kernel writes as a decreasing function of the time interval τ = |t_i − t_j|: $$k_{\mathrm{S}\mathrm{E}}(\tau ;{\mathrm{\Phi }}_{\mathrm{S}\mathrm{E}})={\alpha }_{\mathrm{S}\mathrm{E}}^2\exp (-\frac{{\tau }^2}{2{\lambda }_{\mathrm{S}\mathrm{E}}^2}),$$(2)

with the hyperparameters Φ_SE = {α_SE, λ_SE) representing characteristic amplitude and length scale. By jointly fitting this GP to both nights, we found α_SE = 12.8 ± 6.0 m s⁻¹ and λ_SE = 2.4 ± 0.7 h¹¹. The best-fit model is shown in Fig. 2. The RMS value of the residuals is 4.3 m s⁻¹. Asa sanity check, we checked the consistency of the inferred GP parameters when excluding the data points taken at the beginning and at the end of the nights (which are suspected to be of instrumental origin).

In the following sections of this study, we employ these values as priors to model the short timescale stellar signal. However, when applying model (2) to the RV data collected during the long RV observation campaign (in Sec. 4), we anticipate identifying a signal of diminished amplitude. This expectation stems from the longer exposure time (t_exp > 15 min), which attenuates stellar signals, including oscillations (Chaplin et al. 2019).

Fig. 4 Generalized Lomb-Scargle periodograms of K2, TESS, and SOPHIE RV data for periods of less than one day. The axes are presented in a log-scale format. The yellow and red dashed curves represent best Harvey-function fits to the periodograms, helping for a visual representation of both white Gaussian noise (WGN) and short-term stellar variability.

4 Characterization of the planetary system

In this section, we carried out a three-step analysis of our data. First, we modeled the transits of planet b in the K2 and TESS photometric data to confirm this planet candidate (Sec. 4.1). We then used the planet’s ephemeris inferred from this first analysis as priors for the analysis of the SOPHIE RV to evaluate the best strategy to mitigate the stellar activity noise (Sec. 4.2). We finally analyzed the photometric and RV data (SOPHIE+HIRES) jointly (Sec. 4.3). The final adopted parameters result from this last analysis. They are reported in Table D.1.

4.1 K2 and TESS transit analyses

We began by jointly analyzing the K2 and TESS observations, containing a total of nine transits of planet b. First, using the Box Least-Square algorithm (Kovács et al. 2016), we performed a transit search analysis but did not detect any transit signatures other than those attributed to planet b. In the following, we isolated the transit events to save computational time and avoid complex modeling of the stellar activity signals. We took out all data at 3 and 5 times the transit duration from the transit center in K2 and TESS light curves, respectively (K2 data are strongly affected by instrumental systematics). The individual transits curves are shown in Appendix A.

We used the Planet Analysis and Small Transit Investigation Software (PASTIS; Díaz et al. 2014) to characterize the nine transits of planet b. To account for the different temporal sampling of the K2 (29.6 min) and TESS (2 min) observations, the software oversamples the transit model at the 2-min rate and then calculates the likelihood over the original rate of the input observations.

The spectral energy distribution (SED) was computed using the BT-SETTL stellar atmosphere models (Allard et al. 2012). The host star was modeled using the Dartmouth stellar-evolution tracks (Dotter et al. 2008). The priors on the stellar parameters (T_eff, [Fe/H], density ρ_*) were set to follow Gaussian distributions parameterized by the values given in Table 2. We used a quadratic law to model the stellar limb darkening for each passband, with parameters (u_a; u_b) interpolated from the Claret & Bloemen (2011)’s table. These interpolations are done for each iteration of the stellar parameters.

Regarding the planet parameters, we used Gaussian priors on the ephemeris from Yu et al. (2018), with uncertainties on P_b and T_0,b enlarged by a factor 100. For the eccentricity e_b, we used a truncated zero-mean Gaussian distribution with a dispersion of 0.083 following the recommendation from Van Eylen et al. (2019). For the other parameters (inclination i_b, radius ratio R_b/R_*, argument of periapsis ω_b), we used uniform priors.

As the light curves analyzed here are restricted to observations taken in the vicinity of the transits of planet b, we modeled the variability around each of the nine transits by GPs with a SE-type convariance function (see Eq. (2)), with uniform priors on their hyperparameters.

A total of 40 Markov chains of 500 000 samples were run. Convergence of each chain was ensured by a Kolmogorov-Smirnov test, and the converged chains were then merged after removing half of the samples as a burn-in phase.

The inferred parameters of planet b are reported in Table 3. The joint analysis indicates a transit depth of ~2.2%, which corresponds to a mini-Neptune size planet with R_b ~ 2.8 R_⊕. The best-fitting model is shown in Appendix A. We reconfirm that the residuals show no signature of spot-crossing events during TESS transits.

By performing a linear propagation of the transit ephemeris from the analysis of the K2 observations alone, we found those estimated from the TESS observations analyzed individually (see Table 3, values compatible within 1σ). We therefore measure no transit timing variations induced by a nontransiting nearby exoplanet.

4.2 SOPHIE RV analysis

As a preliminary study based on the SOPHIE time series, we followed the approach of observing the star three times a night and grouping these data points together within each night to mitigate short timescale stellar variability (Dumusque et al. 2011). The resulting binned dataset consists of N = 137 nights of observations, spanning a time coverage of 485 days across two SOPHIE campaigns.

In Fig. 5, we show the GLSP of these observations. The GLSP for RV (top panel) reveals a prominent peak at ~66 days. This period does not align with a harmonic of the stellar rotation period, which is identified as the second most dominant peak at around 9 days (excluding the structure of peaks around one day induced by the time sampling). Importantly, none of the five activity indicators exhibit the presence of a signal with a period ~66 days, as illustrated in the bottom panels. However, they do exhibit strong periodicities at intermediate periods, specifically in the period range of 3–4 P_rot. This is discussed in Appendix C.1. These findings, coupled with supplementary analyses detailed in Appendix C.2, strongly support the existence of a nontransiting planetary candidate¹², called planet c below, with period P_c ~ 66 days. No significant signal is observed in the RV periodogram at the period of the transiting exoplanet (P_b ~ 15.6 days), despite the planet being a mini-Neptune orbiting at close distance to the star. This is a direct consequence of the high level of stellar activity masking the planet’s signature. Detailed modeling of the stellar signal is required to detect a peak at this period (see text below).

We then ran the PASTIS software (Díaz et al. 2014) on this dataset to extract the minimum mass estimates of the two planets. We constrained the ephemeris for the transiting planet based on the photometric data analyses (see Sec. 4.1), incorporating Gaussian priors on P_b and T_0,b, along with truncated Gaussian priors on the eccentricity e_b. Uniform priors were applied to the argument of periapsis ω_b, RV semiamplitude K_b, and the five Keplerian parameters of planet c (K_c, P_c, T_p,c, e_c, ω_c). We employed a quasi-periodic (QP) GP model to capture the stellar variability induced by rotational modulation and evolution of the magnetic regions on the stellar surface (Haywood et al. 2014; Aigrain et al. 2012; Stock et al. 2023). This kernel is defined as: $$k_{\mathrm{Q}\mathrm{P}}(\tau ;{\mathrm{\Phi }}_{\mathrm{Q}\mathrm{P}})={\alpha }_{\mathrm{Q}\mathrm{P}}^2\exp (-\frac{{\tau }^2}{2{\lambda }_{1,\mathrm{Q}\mathrm{P}}^2}-\frac{2}{{\lambda }_{2,\mathrm{Q}\mathrm{P}}^2}{\sin }^2\left[\frac{\pi \tau }{{\mathrm{P}}_{\mathrm{r}\mathrm{o}\mathrm{t}}}\right]),$$(3)

where Φ_QP = {α_QP, λ_1,QP, λ_2,QP, P_rot} represents the set of hyperparameters corresponding to the characteristic amplitude, decoherence timescales, harmonic complexity (or roughness of the signal), and rotation period. The stellar rotation period was constrained by a Gaussian prior based on results from Sec. 3.1. Summary of all the priors used in this work is provided in Table D.1. We used 40 Markov chains of 500000 samples.

The posteriors of the best fitting parameters are shown in Fig. 6. For planet b, we found K_b = 1.5 ± 0.9 m s⁻¹, P_b = 15.537 ± 0.049 days, and e_b = 0.058 ± 0.05 among the inferred set of parameters. For planet c, we found K_c = 15.0 ± 1.8 m s⁻¹, P_c = 66.46 ± 0.44 days, and e_c = 0.07 ± 0.06. The RV jitter is estimated to be σ_SOPHIE = 6.3 ± 1.3 m s⁻¹, far above the initial RV uncertainties (see Table 1). The GP model converged towards a period of P_rot = 9.16 ± 0.17 days, in agreement with Sec. 3.2. It converged towards a characteristic amplitude α_QP = 11.6 ± 2 m s⁻¹, and a decoherence timescale ${\lambda }_{1,\mathrm{Q}\mathrm{P}}={22.1}_{-5.9}^{+11.6}$ days. We note a very wide tail of the posterior distributions of the stellar activity model parameters. This indicates they are not well constrained by RV observations taken at the rate of one point per night.

The left panels of Fig. 7 show the GLSP of RV observations (top), iteratively subtracted (from top to bottom) by the best-fitting models for planet c, stellar activity, and planet b. The final RMS¹³ of the RV residuals is around 4.4 m s⁻¹. In the GLSP of these residuals, we observe strong peaks at short periods. Without correcting the RV time series from these short-term noises, the peak at P_b is not prominent¹⁴ (see third row). This stellar signal remains significant over periods longer than a day, affects the RV characterization of the transiting planet, and needs to be corrected for Meunier et al. (2015).

To this end, we analyzed the SOPHIE observations without grouping the 3 data points by night. Based on the detailed analysis of the two full nights taken at high cadence rates (see Sect. 3.2.2), we modeled short-timescale stellar variability with an SE covariance function (see Eq. (2)). We parameterized Gaussian priors on these two hyperparameters with the values deduced in Sect. 3.2.2. All other priors were kept identical to the analysis of the binned SOPHIE observation.

First, we find that the inferred Keplerian parameters for the two planets are consistent with the analysis of the binned RV data (see Fig. 6). This means that adding a second GP to model stellar activity did not degrade the inferred planetary signals. We also observed narrower posterior distribution of the stellar activity parameters, indicating that the model is more well constrained than previously.

Second, when we compare the GLSPs of the unbinned (right panels of Fig. 7) and binned (left panels) SOPHIE RV dataset, we see the planet b has now the largest peak (third row). Moreover, we see no strong residual periodic component in the GLSP of the RV residuals. This leads us to conclude that the main contribution of short-term correlated variability has been well constrained by the second GP noise model.

Third, we note an RMS of the data residuals of 2.1 m s⁻¹. This RMS corresponds now to the inferred RV jitter σ_SOPHIE = 2.2 ± 0.4 m s⁻¹, and is also in agreement with the initial RV uncertainties (see Table 1).

Fig. 5 Generalized Lomb-Scargle periodograms for RV and activity indicators. GLSP are calculated on the basis of binned SOPHIE data, where we have eliminated long-term variation using a two-degree polynomial fit. From top to bottom, we show the GLSP for the RV, FWHM, the Area of the Gaussian fit to the CCF, Bisector, $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$, and Hα lines. The red, green, and blue vertical lines indicate the stellar rotation period, the orbital period of the transiting planet (b), and the orbital period of the candidate planet (c), respectively.

Fig. 6 Normalized posterior distribution of the parameters fitted to the RV data. The distributions resulting from the analyses of the SOPHIE binned and unbinned observations are shown in red and black, respectively. The distributions resulting from the joint analysis combining the photometric (K2+TESS) and RV (SOPHIE+HIRES) observations are shown in yellow. Top: Five Keplerian parameters of planet b (T0,b, Pb, Kb, eb, ωb) and V0. Middle: Five Keplerian parameters of the candidate planet (T0,c, Pc, Kc, ec, ωc), and the RV jitter (σSOPHIE). Bottom: GP hyperparameters of the stellar magnetic activity model (Prot, αQP, λ1,QP, λ2,QP), and the short-term stellar noise model (αSE, λSE). For the latter, the values resulting from the analysis of the two SOPHIE nights of observation (Sect. 3.2.2) are shown for comparison (blue).

4.3 Combined photometry and RV analyses

Despite the robustness of the inferred results coming from the separate analyses of the K2+TESS light curve and the SOPHIE RV, we performed a joint orbital analysis of the photometric and spectroscopic dataset, along with stellar evolution tracks to refine the parameters and derive self-consistent uncertainties in the model parameters (taking into account the underlying correlations between some of them). We also added to this combined analysis the HIRES RV. Both the SOPHIE and HIRES observations are taken at the rate of three data points per night (the raw HIRES data have been binned into three data points per night, as described in Sect. 2.2.2).

We used again the PASTIS software, with the same setting as in Sects. 4.1 and 4.2: planet b in transit, planet c not transiting, two GPs to model stellar activity signals in the RV, and nine distinct GPs models for the photometric and instrumental variability observed around the transits of planet b. For the stellar parameters (T_eff, logg, [Fe/H]), we used normal priors centered on the values derived from spectral analysis (see Table 2). For distance to Earth (d), we used a normal prior centered on the Gaia DR3 value (see Table 2), and for stellar extinction E(B-V) we used a uniform prior. Priors on each parameter are listed in Table D.1. In total the fitting procedure involved 75 free parameters. The main inferred planet (Keplerian), stellar (GPs), and instrumental parameters are reported in Table D.1.

Overall, the results are in agreement with those derived with the individual analyses of the photometric and RV data.

Concerning the stellar fundamental parameters, the results of this joint analysis are compatible with the results from spectral analysis (Sec. 4.3) withing the 1σ errorbars. Based on isochrones, PASTIS derived an age for the host star of 1.3 ± 0.3 Gyrs, which is not in exact agreement with the age derived from chemical clock relationships (see Sec. 3.1). This could also be compared with the age estimated with gyrochronological relationships, as proposed by Mamajek & Hillenbrand (2008); Angus et al. (2019); or Mathur et al. (2023). A visual examination of Figs. 1 and 5 in Mathur et al. (2023) (which relate P_rot and S_ph to stellar age) confirms that the age of HD 73344 should be between 1 and 2 Gyrs. The age-activity relationship (using $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$) described in Eq. (3) of Mamajek & Hillenbrand (2008) gives an age of ~1.15 Gyrs. Despite the lack of consensus on the stellar age derived by these different techniques (and in the absence of a proper asteroseismology study), all these age estimates nevertheless indicate that HD 73344 is certainly a young star, which is consistent with the high level of activity discussed in Sec. 3.

The orbital period of planet b is refined for the RV analysis thanks to the photometric data, and it gives an RV semiamplitude for the planet smaller than with the analysis based on the RV data alone¹⁵. We find a marginal RV signal for the planet b ($K_{\mathrm{b}}={0.667}_{-0.426}^{+0.559}\mathrm{m}{\mathrm{s}}^{-1}$), despite the fact it is a mini-Neptune ($R_{\mathrm{b}}={2.884}_{-0.072}^{+0.082}{R}_{\oplus }$) planet at short orbit (P_b = 15.611 ± 0.00003). The marginal detection of planet b indicates the need of injection tests to secure the estimate of the planet mass, in presence of strong stellar variability noise (Meunier et al. 2023). If we consider the 3σ uncertainties on K_b, we find a signal of <2.34 m s⁻¹, which is compatible with the RV jitter attributed to both SOPHIE and HIRES instruments. This corresponds to a planet with a mass $M_{\mathrm{b}}={2.983}_{-1.90}^{+2.50}$ M_⊕ (or M_b < 10.48 M_⊕ at 3σ). However, this leaves us with a minor constraint on the bulk planet density of ρ_b < 2.45 g cm⁻³ at 3σ. According to the planet distance to its host star, we estimate the equilibrium temperature of the planet to be around $T_{\mathrm{e}\mathrm{q},\mathrm{b}}={910}_{-7}^{+9}$ K (assuming zero albedo) and $T_{\text{lock },\mathrm{b}}={1066}_{-12}^{+15}$ K if the planet is tidally locked to its star (assuming homogeneous redistribution of heat in the atmosphere; Cowan & Agol 2011).

For planet c, we find an RV signal with an amplitude of K_c = 16.1 ± 1.8 m s⁻¹. The planet is found on a nearly circular orbit (e_c = 0. 061 ± 0. 02). This corresponds to a planet with a minimum mass of $M_{\mathrm{c}}\sin i_{\mathrm{c}}={116.3}_{-13.0}^{+12.8}{M}_{\oplus }\sim 0.37\pm 0.04M_{\mathrm{J}}$.

For the stellar activity parameters, we find values fully consistent with the previous analyses (see Fig. 6). The stellar rotation period is P_rot = 9. 09 ± 0. 04 days. The signatures of the stellar activity sources that are modulated with the stellar rotation have an amplitude around 11.8 m s⁻¹, which is large compared to the signal of the transiting planet. It evolved over long timescales (λ_1,QP ~ 19 days), close to the ~15 days orbital period of the transiting planet. The short-term stellar variability signals is correlated over longer timescales than solarlike stars (λ_SE = 4.0 ± 0.4 h), and generate a significant RV noise (α_SE = 4.8 ± 0.5 m s⁻¹), also above the RV signal of the transiting planet. As we anticipated in Sect. 3.2.2, the amplitude of the granulation signal is smaller than the one derived from the two nights of SOPHIE observations. This is explained by the longer exposure time used during the long observation campaign.

The final RV observations are shown in Fig. 8, and the phased folded transit light curves in Fig. 9. The RMS of their residuals are 1.79 m s⁻¹ and 264 ppm, respectively.

Fig. 7 Generalized Lomb-Scargle periodograms of the SOPHIE binned (left) and unbinned (right) data. From top to bottom: Raw data, then subtraction of the best-fitting Keplerian model for the planet candidate, the stellar activity model(s), and planet b. For the stellar activity model(s), we considered the mean of the predictive GP distribution resulting from the best fits. The stellar rotation period, and the orbital period of planets b and c are indicated by the red, green and blue vertical lines, respectively. In each panel, the red dots indicate the highest periodogram peak taken over periods >1 day. We note that the planetary and activity signals are jointly estimated from the RV data, and then iteratively removed in the successive panels.

Fig. 8 Radial velocity of HD 73344 resulting from the joint analysis of the photometric and RV data. Panel a: RV of SOPHIE (black) and HIRES (yellow) observations and best fitting model (purple). Panel b: RV residuals. Panel c: RV phased at the period of planet b (planet c and activity models subtracted). Panel d: RV phased at the period of planet c (planet b and activity models subtracted). Best-fitting models for planets b and c are shown in red.

Fig. 9 Phased light curve of HD 73344b resulting from the joint analysis of the photometric and RV data. K2 data are shown in black, TESS data in blue. These light curves have been corrected by the best-fitting noise models. The best-fitting transit model is shown in red. Yellow dots represent the 120min binned light curves (K2 and TESS combined).

4.4 Spitzer data analysis

We analyzed the Spitzer photometry of HD 73344 using the POET¹⁶ (Photometry for Orbits, Eclipses, and Transits) code (Stevenson et al. 2012; Cubillos et al. 2013), following the same analysis approach as used by Crossfield et al. (2020). We identified a 2.5 pix aperture with 0.01 pix resolution on the pixel map as giving the optimal precision. Because the available K2 and TESS light curves give tighter constraints, we held all transit parameters fixed except the time of transit and R_b/R_*. The quadratic limb-darkening coefficients were constrained by Gaussian priors, set to the mean and standard deviation of all the values tabulated by Claret et al. (2013) for model grid parameters closest to the stellar parameters of HD 73344 − that is, u_a = 0.0106 + 0.0670 and u_b = 0.229 + 0.142. The final detrended light curve and best-fit transit model are shown in Fig. 10; the median and standard deviation on the derived parameters are T_0,b;Spitzer = 2458704.7289 ± 0.0014 BJD and R_b;Spitzer/R_* = 0.0203 ± 0.0017. The transit depth observed by Spitzer in the NIR being fully compatible with the optical one measured by K2 and TESS, it confirms the achromatic property of this planetary transit (Fressin et al. 2012; Désert et al. 2015) which nature is now considered as validated.

5 Discussion

In this section, we first discuss the stability of the two-planet system (Sect. 5.1). We then study how we can constrain the inclination of the candidate planet from the transit probability of planet b (Sect. 5.2). We conclude with a general discussion of the composition of the transiting planet’s interior, bearing in mind that our knowledge of this composition is drastically limited by the impact of stellar activity on the measurement of the planet’s mass (Sect. 5.3).

5.1 Stability of the two-planet system

We investigated whether it is possible to refine the orbital parameters of planets HD 73344 b and c using constraints on their orbital stability (see e.g., Stalport et al. 2022), and perhaps constrain the true mass of planet c. To this end, we used a similar approach to what has been done to other two-planet systems by Correia et al. (2005, 2009), Laskar & Correia (2009), or Couetdic et al. (2010).

The stability of a dynamical system can be quantified by running a frequency analysis on the output of a numerical integration (Laskar 1988, 1993, 2003). In our case, we are interested in the short-term stability of the planets¹⁷, and so it is enough to analyze their mean longitude λ, which is related to their orbital motion. For a given planet, we define the stability coefficient $$\delta =\frac{|n_2-n_1|}{n_0},$$(4)

where n₁ and n₂ are the ‘mean’ mean motions (i.e., the linear part of λ) obtained by running a frequency analysis on the first and second half of a numerical integration, and n₀ is a reference value, taken here to be the mean motion of the planet for its best-fit parameters. For a stable system, δ should be 0 up to the numerical accuracy of the analysis, while values close to 1 or above denote strongly unstable systems. As planet b is much less massive than planet c, it is much more sensitive to chaos, so we focus here on its stability coefficient δ_b.

We integrated the system with the numerical scheme SABA(10,6,4) of Blanes et al. (2013) implemented in the REBOUND package (Rein et al. 2019), and performed the frequency analysis with the dedicated function of TRIP (Gastineau & Laskar 2011). We set the duration of our integrations to 200 yr, which represents about 4700 orbits of planet b and 1100 orbits of planet c. On such a short duration, tidal dissipation and general relativistic precession can be neglected. We checked that the oblateness of the star (expected to be J₂ ≈ 10⁻⁷; see Batygin & Adams 2013; Spalding & Millholland 2020) also produces negligible orbital perturbations compared to planet-planet interactions. We first mapped the parameter space around the best-fit solution, by varying the semimajor axis and eccentricity of planet c on a regular grid while all other parameters were set to their nominal values. As planet c is only detected in radial velocity, there is a degeneracy between its mass M_c and inclination i_c with respect to the sky plane; hence, we repeated the same experiment for several values of i_c and modified M_c accordingly. In these simulations, we chose equal longitudes of node in the sky plane for the two planets, such that their mutual orbital inclination is simply i_b − i_c, where i_b ≈ 88° (see Table D.1).

Figure 11 shows the result obtained for i_c = 88° and 10°. We checked that integrating over a longer duration (e.g., 2000 yr) does not substantially alter our maps; this shows that 200 yr is long enough here for the frequency analysis to give a pertinent result. We rule out a strong mean-motion resonance between the two planets. The closest large resonance visible in Fig. 11 is the 4:1 mean-motion resonance (vertical structure on the left of the best-fit location), but it is more than 3-σ away from the most probable parameters of the planets. For the coplanar configuration, for which planet c has its minimum mass, the best-fit solution lies in a very stable region (Fig. 11a). If the sky-plane inclination of planet c is small, however, its corresponding large mass and large mutual inclination with respect to planet b are a source of instability (Fig. 11b). We must therefore quantify the emergence of this instability.

Using the posterior distribution of the system’s parameters obtained from the joint analysis in Sect. 4.3, we computed the histogram of log₁₀ δ_b for various inclinations of planet c from i_c = 88° (coplanar case) to i_c = 1^o (almost perpendicular case). Examples of the histograms obtained can be found in Appendix E. For inclination values i_c larger than about 30°, the posterior distribution of log₁₀ δ_b has a single peak located below −4; this means that the whole sample is stable (similarly to the black and dark blue regions in Fig. 11). For 5° ≲ i_c ≲ 30°, an unstable subsample appears as a second peak located above −4. As we decrease i_c, and therefore increase the mass of planet c, this unstable subsample grows. For i_c ≲ 5°, there is again a single peak in the distribution, but located above −2; this means that the whole sample is now unstable (similarly to the yellow and red regions in Fig. 11).

From this analysis, we deduce that the dynamical stability of the system would be able to constrain the planets’ parameters only if i_c ≲ 30°. However, such a small inclination i_c would correspond to a very large mutual inclination between the two planets, which we consider unlikely. First, the statistical distribution of multiplanetary systems shows that planets having small eccentricities tend to have small mutual inclinations, and vice versa (Xie et al. 2016). This can be understood by the statistical equipartition of angular momentum deficit as a result of chaotic diffusion (see Laskar & Petit 2017). Second, a large mutual inclination would result in a fast precession of the orbital plane of planet b in and out of transiting configuration, which would reduce its transit probability (see e.g., Becker & Adams 2016). In the next section, we use this last property to put more stringent constraints on the unknown parameters.

Fig. 10 Spitzer observations of HD 73344. Top: raw Spitzer light curve (black points) of HD 73344b and best-fit models with (blue) and without (black) a transit fit. Middle: detrended Spitzer data and best-fit transit light curve. Bottom: residuals to the fit.

Fig. 11 Stability of the HD 73344 system as a function of the semimajor axis and eccentricity of planet c. The cross shows the best-fit values and their 1-σ uncertainty interval. The colour scale is chosen so that dark blue denotes stable orbits and red is highly unstable. The inclination chosen for planet c is labelled; its corresponding mass is Mc = 0.36 MJ for panel a and Mc = 2.05 MJ for panel b.

Fig. 12 Transit probability of planet HD 73344 b computed as the fraction of time its impact parameter is smaller than 1. The map is drawn as a function of the unknown longitude of node and inclination of planet c with respect to the sky plane. Other parameters are set to their best-fit values (see Table D.1). The transit probability (colour scale) is obtained from a 50 kyr numerical integration (see text). Points are coloured white if a planet is ejected before the end of the simulation.

5.2 Transit probability of planet b

During the orbital evolution of the system due to mutual planetary perturbations, the fraction of time the orbit of planet b passes in front of the star (as observed today) gives an indication of the likelihood for the considered parameters. Assuming that the system does not contain additional unseen planets, the only parameters that are unconstrained by transit and RV data are the inclination of planet c (linked to its mass through M_c sin i_c) and the longitudes of node of the two planets in the sky plane. As the choice of origin for measuring the longitudes is arbitrary, only the difference Ω_c − Ω_b actually matters, which reduces the unknown parameters to only two.

Figure 12 shows the transit probability of planet b as a function of the two unknown parameters. For each pixel of the figure, a numerical integration is performed over 50 kyr and the fraction of time steps the orbit of planet b passes in front of the star is recorded¹⁸. We used the integration scheme SABA(10,6,4) of Blanes et al. (2013), with the inclusion of the general relativistic precession implemented in the same way as Saha & Tremaine (1994).

The circular features in Fig. 12 roughly correspond to curves of constant mutual inclination between the two orbits. As obtained in Sect. 5.1, a large mutual inclination is ruled out because it would make the system unstable. Figure 12 shows that the transit probability of planet b sharply peaks at 100% in two very small regions. These regions correspond to near coplanarity between the orbits of the two planets, which are either prograde (Ω_c − Ω_b ≈ 0°) or retrograde (Ω_c − Ω_b ≈ 180°) between each other. We point out that observational data cannot tell whether the inclination value of planet b is i_b or 180° − i_b. In the latter case, the corresponding transit probability map is obtained from Fig. 12 by the transformation i_c → 180° − i_c and Ω_c → 180° + Ω_c.

A zoom-in view of the regions of highest transit probability can be seen in Fig. 13. As planet c is not observed to transit the star, its inclination today is necessarily smaller than about 89.1° (or larger than 180° minus this value). Then, if we require the transit probability of planet b to be $\mathcal{P}$ > 0.9, Fig. 13 gives that the inclination of planet c must be i_c ≳ 87.9° (or smaller than 180° minus this value). Likewise, we obtain i_c ≳ 87.4° for $\mathcal{P}$ > 0.5, and i_c ≳ 80.0° for $\mathcal{P}$ ≳ 0.3. The constraint obtained here is therefore much more stringent than the mere dynamical stability of the system (see Sect. 5.1).

To estimate the uncertainty on the probability values depicted in Fig. 13, we generated a similar map where, for each pixel, we propagated the full posterior distribution of the MCMC fit to the data (see Sect. 4.3) instead of just the best-fit solution. This was made possible by reducing the resolution of the map to 15 × 15 and by propagating the trajectories using the Lagrange-Laplace theory (see e.g., Murray & Dermott 1999). The statistics obtained for each pixel show that the 1-σ uncertainty on the probability value is everywhere smaller than 0.01. This uncertainty translates into a 2-σ range of less than 0.1° on the values of i_c cited above for $\mathcal{P}$ > 0.9 and $\mathcal{P}$ > 0.5, and a 2-σ range of about 2° on the value cited for $\mathcal{P}$ > 0.3. These very small uncertainties mainly come from the tight constraint that we have on the inclination of planet b (see Table D.1).

The above analysis translates into a likelihood estimate for the inclination and mass of planet c; however, we still do not have a direct measurement of their values. A way to break this degeneracy would be to detect the nodal precession of planet b through a variation of its impact parameter b_b over time (see e.g., Judkovsky et al. 2022). We examined the behavior b_b in our simulations, and conclude that for an uncertainty of about 0.02 on the measurement of b_b (see Table D.1), we would need observations spanning at least 10 yr in order to detect a substantial noncoplanarity of several degrees between the orbits of the two planets. Conversely, no variation in b_b should be detectable if the system lies in the most likely region of Fig. 13, which corresponds to near coplanarity between the two orbits.

Fig. 13 Same as Fig. 12, but showing enlarged views. In the hatched band, planet c would be observed to transit the star; this region is therefore excluded.

5.3 Internal composition of planet b

This section proposes a preliminary investigation of the internal composition of planet b. However, it is crucial to bring to the reader’s attention the significant degeneracy observed among the various parameters (e.g., core mass fraction and envelope size) of planet internal composition models. This degeneracy is particularly critical when dealing with planets with large uncertainties on their fundamental parameters (mass, radius), as it is the case here with the large uncertainty in the planet mass.

In the mass-radius diagram, planet b stands out the rocky super-Earth population and appears as a likely member of the sub-Neptune population (see Fig. 14). At first, we compare the position of planet b in the mass-radius diagram with simple mass-radius relationships. Given the mean density of the planet and its strong irradiation (T_eq,b = 910 ± 7 K assuming a zero albedo, see Sect. 4.3), we considered the mass-radius relations inferred from water-rich composition models (Mousis et al. 2020; Aguichine et al. 2021; Acuña et al. 2022). These models assume a water-dominated atmosphere either in vapour or supercritical state, on top of a high-pressure water layer or a mantle. For interior analysis, we computed the mass-radius relationships by employing Eq. (29) from Aguichine et al. (2021) for different H₂O mass fractions at an equilibrium temperature of T_eq = 900 K. We assumed a core mass fraction of 30%, which aligns with the characteristics of an Earth-like interior. Results are shown with the colored lines in Fig. 14. This preliminary comparison shows that HD 73344b is compatible with a very high water-mass fraction of at least 80%. A more quantitative analysis also confirms that the water content is at least 75% (see Appendix F).

However, such a high value is well above the most water-rich bodies in the solar system, such as the icy moons. This suggests that a water-dominated atmosphere is not inflated enough to account for the low density of the planet. We performed then a similar analysis but assuming that the envelope is made of gas of solar composition (i.e., dominated by H₂ and He), on top of a core of Earth-like composition. To do so, we used the tabulated mass-radius relations from the interior structure model of Lopez & Fortney (2014). Intermediate mass-radius relations were obtained by linear interpolation on the grid provided in Lopez & Fortney (2014). In this case, we found that HD 73344b properties are compatible with an envelope that accounts for 2–3% of the mass of the planet (see gray lines in Fig. 14 and Appendix F).

If we consider the scenario described above as realistic, it implies that HD 73344b lies in a parameter space where extreme atmospheric escape of hydrogen is expected (Fossati et al. 2017; Zahnle & Catling 2017; Rogers et al. 2023). The restricted Jeans escape parameter can be computed as Λ = GM_pm_H/(k_BT_eqR_p) (Fossati et al. 2017), where m_H is the mass of the proton. Studies of Owen & Wu (2016); Cubillos et al. (2017); Fossati et al. (2017) conclude that planets with Λ < 10–15 are hydrodynamically unstable and their atmospheres must experience extreme atmospheric escape (e.g., Parker wind or boil-off). The analysis of Vivien et al. (2022) showed that pure water atmospheres can remain stable down to values of Λ as small as ~0.5. For $\mathrm{H}\mathrm{D}73344\mathrm{b}\mathrm{\Lambda }={8.5}_{-5.6}^{+6.0}$ meaning that a hydrogen dominated atmosphere is unlikely. To account for the planet’s unusually large radius, the atmosphere and envelope would have to be a mixture of hydrogen and heavier volatile elements. In this case, an atmosphere with a higher mean molecular weight would result in a reduced escape rate, potentially explaining how this planet has retained its volatile envelope. This underlines the need for spectroscopic characterization to break the degeneracy on the planet’s composition. It also shows that HD 73344b is a promising target for testing the escape of H₂ and He. We calculated HD 73344b’s Transmission Spectroscopy Metric (TSM; Kempton et al. 2018) while propagating all parameter uncertainties and found $\mathrm{T}\mathrm{S}\mathrm{M}={260}_{-30}^{+640}$. Despite the large upper error bar (corresponding to lower planet masses), this metric suggests that HD 73344b could be a particularly promising target for transmission spectroscopy. However, transmission spectra depends on the scale height of the atmosphere, which depends on the planet mass. It is therefore necessary to improve the measurement of the planet mass (see discussions in Almenara et al. 2022 for similar conclusions) to reach, at the very least, the 50% level of accuracy on the mass required to produce reliable atmospheric retrievals (Di Maio et al. 2023). Moreover, we expect the transmission spectra to be contaminated by stellar activity signatures, while this may complicate the interpretation of such observations it could also bring very useful information on the chromaticity of the stellar signal.

Fig. 14 Mass–radius diagram of confirmed exoplanets with mass and radius precision better than 3σ and 10σ, respectively. Planetary parameters (black dots) are taken from the NASA Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu), and updated to October 2023. HD 73344b is represented by the red dot. The red arrow represents the upper 10σ limit on its mass. Composition models calculated from Aguichine et al. (2021) with a core mass fraction of 30% and various H2O mass fraction (from 10% to 80%) are represented by coloured lines (see text). Composition models calculated from Lopez & Fortney (2014) for a 1 Gyr-old system, Teq,b ~ 900 K, and an envelope mass fraction fenv = [1, 2, 3, 4, 5] % are shown in gray. For comparison, the relation for an Earth-like composition (32.5% core, 67.5% mantle) from Brugger et al. (2017) is shown by the green dashed curve, and some solar system planets (Venus, the Earth, Uranus and Neptune) are displayed.

6 Conclusions

We observed the bright star HD 73344 (V_mag = 6.9) with SOPHIE and HIRES in order to confirm the transiting planet with a period of P_b ~ 15.6 days, which was initially a candidate in the K2 data (Yu et al. 2018). This planet was also confirmed by contemporary TESS (TOI 5140.01) and Spitzer observations. Our main results are listed below.

• Analysis of the spectroscopic SOPHIE and HIRES spectra made it possible to refine the parameters of the host star and measure the stellar abundances.

• In the RV data, a candidate planet with $P_{\mathrm{c}}={66.456}_{-0.250}^{+0.100}$ days and a minimum mass of M_c sin i_c = 0.37 ± 0.04 M_J is detected. A two-planet system in this configuration is dynamically stable if i_c ≳ 30° (which translates into M_c ≲ 0.7 M_J). Moreover, imposing that planet b transits the star (as observed) more than 50% of the time requires a very near coplanarity between the two planetary orbits, such that 87.4° ≲ i_c ≲ 89.1°. This inclination interval translates into a very tight range for the true mass Mc –much tighter than the observational uncertainty on M_c sin i_c.

• In the RV data, the variability of the host star completely masks the RV signal of the transiting planet. Both spectroscopic and photometric data show that the star is indeed particularly active. The rotation period (P_rot ~ 9 days) is close to the orbital period of the transiting planet. The coherence of the activity signal extends to 3–4 Prot and is particularly strong in the area of the Gaussian fit to the CCF indicator (see Appendix C). According to Costes et al. (2021), this may indicate that HD 73344 is faculae-dominated.

• The prevailing approach in RV observations today involves averaging a limited number of data points gathered over a single night to mitigate the influence of short-term stellar variability, such as p-modes, granulation, and supergranulation. However, this strategy falls short when the goal is to detect planetary signals of just a few meters per second around evolved stars. In our efforts to better understand and characterize the short-term variability exhibited by the star HD 73344, we conducted observations spanning two entire nights using SOPHIE. In doing so, we identified a signal with an amplitude of 12.8 ± 6.0 m s⁻¹, and we observed variation in this signal over a coherence time of 2.4 ± 0.7 h. These findings provide us with robust constraints for modeling these sources of noise when employing nonbinned RV data.

• Tests based on SOPHIE RV data reveal that not binning the data (1) provides a more accurate estimation of stellar activity, and (2) yields planetary parameters that are consistent with the binned case, indicating that binning does not enhance precision in planetary parameter determination. The planetary signal is clearly evident in the periodograms when using unbinned data (unlike the binned case), as the high-frequency noise has been effectively modeled. We note, however, that without the use of priors on P_b (coming from photometry), the planet RV signal at P_b is not detected.

• The joint analysis of photometric and spectroscopic data using a model featuring two planets and two Gaussian processes (one for capturing the effects of rotationally modulated magnetic activity and another for the variability occurring on short timescales) allows a more comprehensive understanding of both the planetary system and the overall activity of the host star. For the transiting planet, we infer a radius of $R_{\mathrm{b}}={2.884}_{-0.072}^{+0.082}{R}_{\oplus }$ and amass of $M_{\mathrm{b}}={2.983}_{-1.905}^{+2.500}{M}_{\oplus }$ (marginal detection). This gives an average density of ${\rho }_{\mathrm{b}}={0.681}_{-0.438}^{+0.590}\mathrm{g}/{\mathrm{c}\mathrm{m}}^3$, which is consistent with the density expected for gaseous planets.

• Our initial assessment suggests the presence of an atmosphere enriched in volatile gases, such as hydrogen and helium. However, due to the significant uncertainty regarding the mass of planet b, we refrained from conducting an in depth analysis of its internal structure, including a detailed examination of its atmospheric composition (e.g., employing models from Acuña et al. 2022). This also underscores the critical importance of obtaining a precise estimate of the planet’s mass in order to reveal its true nature.

As a perspective of this study, we note that alternative data analysis techniques based on multidimensional Gaussian processes could improve the RV detection of the transiting planet (see e.g., Rajpaul et al. 2015; Barragán et al. 2021; Hara & Delisle 2023). Not yet implemented in the PASTIS software, such analyses were beyond the scope of the present study. Another perspective is to carry out injection tests to secure the uncertainty on the inferred planet mass (see Meunier et al. 2023), and to optimize the observational strategy for future observations of this system.

Finally, HD 73344 is a very bright star, and planet b is a sub-Neptune planet with an ideal orbital period for future observations by JWST and/or ARIEL (according to the TSM metric). However, the high activity level of the host star may complicate the interpretation of transmission spectra (Rackham et al. 2023; Rackham & de Wit 2023), and refining the planet mass should be considered first. On the other hand, if observed with ARIEL, this target could serve as a benchmark for testing stellar activity diagnostic tools and correction methods for transmission spectra (Cracchiolo et al. 2021; Thompson et al. 2024).

Acknowledgements

The authors thank the anonymous referee for her/his helpful comments that improved the quality of the paper. We warmly thank the OHP and Keck staff for their support on the observations. In particular, we thank Dr. L. Bouma, Dr. E. Petigura, B. Famaey, and Dr. J.-B. Salomon for the time they devoted to the observations of HD 73344 with the HIRES and SOPHIE spectrographs. M.S. thanks Aurélie Astoul and Yubo Su for valuable discussions about the stellar structure and star-planet interactions. We also thank Dr. F. Thevenin for the helpful discussion about Gaia’s observations of this system. This paper includes data collected by the K2 and TESS missions, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for these missions is provided by NASA’s Science Mission directorate. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. This work is based in part on observations made with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. S.S., N.M., and D.M acknowledge support from the Programme National de Planétologie (PNP), and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU. N.C.S, S.S., E.D.M. and V.A. acknowledge funding by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work was supported by FCT – Fundação para a Ciência e Tecnologia through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UIDB/04434/2020; UIDP/04434/2020. E.D.M. acknowledges the support from FCT through Stimulus FCT contract 2021.01294. V.A. was supported by FCT through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UIDB/04434/2020; UIDP/04434/2020; 2022.06962.PTDC. This material is based upon work supported by NASA’S Interdisciplinary Consortia for Astrobiology Research (NNH19ZDA001N-ICAR) under award number 19-ICAR19_2-0041. S.M. acknowledges support by the Spanish Ministry of Science and Innovation with the Ramon y Cajal fellowship number RYC-2015-17697, the grant number PID2019-107187GB-I00, the grant no. PID2019-107061GB-C66, and through AEI under the Severo Ochoa Centres of Excellence Programme 2020-2023 (CEX2019-000920-S). J.C. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 757561). S.B. gratefully acknowledges support from MIT International Science and Technology Initiatives (MISTI) and the MIT-France program. S.D. is funded by the UK Science and Technology Facilities Council (grant number ST/V004735/1). The project leading to this publication has received funding from the Excellence Initiative of Aix-Marseille Université–A*Midex, a French “Investissements d’Avenir program” AMX-21-IET-018. This research holds as part of the project FACOM (ANR-22-CE49-0005-01_ACT) and has benefited from a funding provided by l’Agence Nationale de la Recherche (ANR) under the Generic Call for Proposals 2022. We acknowledge funding from the French ANR under contract number ANR18CE310019 (SPlaSH). This work is supported in part by the French National Research Agency in the framework of the Investissements d’Avenir program (ANR-15-IDEX-02), through the funding of the “Origin of Life” project of the Grenoble-Alpes University. O.V. acknowledges funding from the ANR project ‘EXACT’ (ANR-21-CE49-0008-01), from the Programme de Planétologie (PNP) and from the Centre National d’Études Spatiales (CNES). E.W. acknowledges support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, G.A. no. 772293 http://www.asterochronometry.eu).

Appendix A Individual transits of HD 73344b seen by K2 and TESS

Figure A.1 illustrates the individual transits of HD 73344b observed by K2 (six transits) and TESS (three transits). The corresponding residuals are displayed in the right panels. The RMS values of the residuals are (from top to bottom): [80, 62, 87, 90, 17, 40] ppm for K2 transits, and [229, 227, 341] ppm for TESS transits.

The K2 observations exhibit pronounced instrumental systematics that distort the transit shapes, particularly noticeable near the transit bottoms. The details of these systematics are poorly captured by the GP noise models described in Sec. 4.1. However, these systematics remain of relatively low in amplitude (residuals <200 ppm if we zoom into the right panels), which is attributed to their averaging effect over the 30-min K2 integration time. In contrast, TESS observations offer high-quality data with precise transit events characterization, owing to its 120-s sampling rate.

The K2 and TESS data present notable signals of stellar variability, underlining the importance of incorporating them into joint modeling with planetary transits. However, K2 observations are also particularly sensitive to instrumental systematics. This complicates the modeling that also needs to take into account the different bandwidths and integration times of K2/TESS observations. Given our aim to conduct a comprehensive analysis integrating photometric and RV data, we opted for a pragmatic approach in Sec. 4.3: i.e., rather than using overly complex noise models for the photometric part, we used simpler GP noise models for each individual transits. These models effectively capture the global variability evolving around individual transits. Nonetheless, there is scope for further refinement in future analyses.

Appendix B Stellar abundances from HIRES spectra

The stellar abundances of HD 73344, extracted from the SOPHIE and HIRES spectra (with KeckSpec), are reported in Table B.1. They were extracted using two different codes, which are compared here. For the chemical elements extracted from the two spectra, all are in agreement at 1σ, with the exception of oxygen. However, this is easily explained by the sensitivity of oxygen abundances to different line indicators.

From the stellar abundances derived from the SOPHIE spectra, we also estimated the stellar age from 3D chemical clock formulas based on T_eff and [Fe/H] (see Table 10 of Delgado Mena et al. 2019). The results are presented in Table B.2. These ages are consistent with the value given in Table D.1. The weighted average age of 2.0 ± 0.2 Gyr is however twice the value extracted from isochrones analysis with PASTIS (see Sec. 4.3 and discussion therein).

Appendix C Temporal evolution of stellar activity and origin of the ~66 days signal

In this appendix, we first analyze the temporal evolution of stellar activity over the two years of observations with the SOPHIE spectrograph. Then, we discuss the origin of the ~66 days signal spotted in Sec. 4.2. Finally, we look at the impact of stellar activity on the planet derived parameters.

C.1 Temporal evolution of stellar activity

The temporal evolution of stellar activity in the RV observations is here compared with that of the chromospheric indicators discussed in Sec. 4.2. We used the binned RV observations and the GLSP computed as in Eq. (1). These observations were obtained over two long campaigns of 144 and 115 days respectively, separated by a loss of observations during 224 days. First campaign (C1) and second campaign (C2) of observations contain 78 and 61 individual nights, respectively. On each time series, we removed all signals generated by long-term changes in stellar activity with a second-order polynomial, following the recommendations given by Costes et al. (2021).

Fig. A.1 Light curves of HD 73344 showing the individual transits of planet b as seen with K2 (6 transits; top left panels) and TESS (3 transits; bottom left panels). The best-fitting models are shown in red. The corresponding residuals are shown in the right panels.

In Fig. C.1, we show the RV data (first line) followed by the chromospheric indicator data (lines 2 to 7). The last row represents a linear trend introduced to track the influence of gaps in the frequency domain. The last three columns show the periodograms of each series calculated for the two campaigns separately and jointly (Sec. 4.3). The stellar rotation period, its first harmonic at half the rotation period, and the orbital periods of the planets are highlighted by large colored vertical bars to facilitate visual inspection of the results.

During the first campaign, we observe significant peaks at the star’s rotation period in all indicators. However, these peaks lose significance during the second campaign, while we see that other structures appear both at short and long periods (in $\log {\mathrm{R}}_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ for instance). This suggests a strong temporal evolution of the star’s activity between the two campaigns and potentially the presence of several active regions at different latitude of the stellar surface during the second campaign.

We also observe a non-negligible activity signal close to the transit period of the planet, particularly in the Area indicator (in both campaigns). During the second campaign, a prominent peak at P_b is also observed in the $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ and Hα activity indicators. For all activity indicators, however, the peak at P_b disappears when the two campaigns are analyzed together. This shows the temporal decoherence of stellar activity at this particular period, which may have helped to disentangle the planet’s coherent signal at P_b from the activity signal in our previous analyses (Sec. 4).

Close to the candidate planet’s period P_c, weak periodicities also appear in some activity indicators and in the linear trend signal. However, all such peaks are absent in the joint analysis of the two campaigns.

In the joint analysis, we observe a set of peaks in some indicators at periods around ~3 to 4 P_rot. These peaks are particularly visible in the Area indicator. It is not easy to interpret the evolving signatures of these activity indicators in periodograms and associate them with a clearly identified stellar origin (spot/faculae). Placing HD 73344 in Radick et al. (2018) correlation plot between luminosity and $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ (see their Fig. 14, with B-V ~0.547 for HD 73344), we see that HD 73344 is located in the transition regime between spot- and faculaedominated stars. However, spots lifetime rarely last more than two stellar rotations, in comparison with faculae that have a much longer lifetime (e.g., from a few months up to a year in the case of the Sun; Collier Cameron et al. 2019). In addition, F-type stars are known to be predominantly faculae-dominated (see e.g., Costes et al. 2021) due to their shallow convective envelope. Thus, we suggest that the peaks seen around ~3 to 4 P_rot in the periodogram of some activity indicators may come from active faculae regions.

C.2 Origin of the RV signal at ~66 days period

In Sec. 4 we spotted a periodicity at ≈66 days. Indeed, there are several aspects that point towards a sustained periodicity at ≈66 days. First, the peak appears in the GLSP of the SOPHIE RV data of each campaign, and jointly (see Fig. C.1). The phase of the periodic signal is found to be consistent for the two individual campaigns (73° ± 6° for campaign 1 and 55° ± 26° for campaign 2¹⁹). This coherence is not found, for instance, for the ≈9 days peak, which is present in the Area indicator for individual campaigns, but not jointly.

Second, we see that the linear trend signal presents a peak at ≈66 days (at least in individual campaigns), so that the ≈66 d signal in the RV could simply be the trace of a remaining trend in the data. To check whether this is correct, we extended the GLSP to the case where the model includes not only a constant but also a trend (that is, the GLSP P(ν) now computes the reduction in the RSS when jointly fitting constant plus trend plus sinusoid at frequency ν, with respect to the RSS when fitting only constant plus trend). As a result, the ≈66 d signal remained²⁰.

Third, since the activity indicators capture (some of the) frequency contents about activity, we investigated whether these signals could be used to capture the 66 d periodicity of the RV data. To do so, we computed another modified (‘activity aware’) GLSP, where the model now includes all the activity indicators in addition to the constant and the linear trend of the previous case. Again, the result showed a dominant peak at ≈66 days.

By construction, activity indicators can trace activity signals but not planetary signals (see e.g., Queloz et al. 2001; Boisse et al. 2011) and the results above support the fact that the 66 d peak is not caused by activity. If this signal is indeed caused by stellar activity, we cannot explain why it is invisible in the activity indicators.

Another possibility, besides a planetary signature, is that this periodic signal is caused by an unmodeled instrumental noise. To our knowledge, no instrumental variation at this timescale is known, and a coherent signal with an amplitude of K ≳>~ 16 m/s coming from an instrumental systematics would have been easily spotted by the SOPHIE technical team. As an additional check, we analyzed the RV data taken with the HIRES spectrograph alone to exclude possible instrumental systematics from SOPHIE. We indeed observe that the ≈66 d peak is also present in the GLSP of HIRES data (third largest peak, not shown here), but many other high peaks also exist, which prevents us from drawing a very clear conclusion on the interpretation of these data. We note that the RV from the telluric lines alone do not show any ≈66 d periodicity.

In summary, there is a periodic signal at ≈66 d in the RV observations. We cannot exclude the hypothesis of a nonplan-etary origin, but none of the dataset we analyzed favor this hypothesis. We assume then that the peak spotted in the GLSP at P_c ~ 66 days period is the signature of a non transiting planet.

As a final note, we comment on the ‘significance’ of this peak. Any statistical significance estimation algorithm relies on a noise model. In our case, the main question is whether the 66 d is caused by the stellar activity noise or not. If we test this hypothesis with a noise model free from this component, and investigate the probability that such a large peak occurs with this noise (that is, the p-value of this peak), then this peak is so large that it is declared highly significant by standard procedures. For instance, we obtain a p-value less than 10⁻⁴ for a classical bootstrap procedure based on the permutation of the residuals. This procedure is often used but provides reliable estimates only if the noise is white. The more elaborated 3SD procedure (Sulis et al. 2022) based on an activity noise model composed with the two estimated GPs also leads to a low p-value (8 × 10⁻⁸). However, considering these p-values as a definitive support for the significance of this peak would be somewhat adventurous for such an active star, as these values are strongly conditioned to the noise model, which is not well constrained.

Fig. C.1 Temporal evolution of the stellar activity. From top to bottom: RV data, FWHM, the Area of the Gaussian fit to the CCF, Bisector, $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ Hα lines, and linear trend. From left to right: Time series, GLSP for first campaign, second campaign and both campaigns. In all time series, we have eliminated long-term variation by a two-degree polynomial fit. The large red, green and blue vertical lines indicate the stellar rotation period (and half its period), the orbital period of the transiting planet (b), and the orbital period of the candidate planet (c), respectively.

C.3 Impact of the temporal evolution of stellar activity on planet-derived parameters

As a final test, we evaluated the impact of the temporal evolution of stellar activity on planet-derived parameters. This assessment involved comparing outcomes derived from SOPHIE RV data collected during the first campaign (C1), the second campaign (C2), and the combined dataset (C1+C2; as detailed in Sec. 4.2). Our benchmark was established through the joint analysis integrating both photometric (K2, TESS) and nonbinned RV data (SOPHIE, HIRES), as outlined in Sec. 4.3.

The distribution of pertinent posterior probabilities resulting from these diverse analyses is shown in Fig. C.2, incorporating both binned and unbinned RV data.

The RV signature of planet c is not well constrained when examining SOPHIE’s individual RV campaigns, as their duration closely aligns with the planet orbital period P_c.

We see that, during the second SOPHIE observing campaign (represented by cyan and pink posteriors), the stellar activity signal exhibited a notably higher amplitude (see GP amplitude and Sect. C.1). As a result, this has been accompanied by larger uncertainties in the planetary parameters (see the enlarged posterior tails in Fig. C.2). We also note a smaller difference between the parameters derived during the C1 and C2 campaigns when using unbinned RV data (2 GP) compared to binned RV data (1 GP).

We see that a better parameter convergence is consistently observed in both campaign when utilizing the nonbinned dataset (comparing blue to magenta; or cyan to pink).

Ultimately, we present the RMS of the RV residuals associated with the various analyses in Table C.1. Notably, the application of the two-GP noise models on the unbinned RV dataset considerably improves the quality of our noise modeling and leads to similar residual RMS for both campaigns; while larger differences are observed on the RMS of the binned RV dataset between C1 and C2.

Fig. C.2 Normalized posteriors distribution of the parameters fitted to the RV data. The distributions resulting from the analyses of the SOPHIE binned observations are shown in magenta (C1), pink (C2), and red (C1+C2). The distributions corresponding to the analyses of the unbinned data are shown in blue (C1), cyan (C2), and black (C1+C2). The distributions resulting from the joint analysis combining the photometric (K2+TESS) and nonbinned RV (SOPHIE+HIRES) observations are shown in yellow (benchmark). Top: the five Keplerian parameters of planet b (T0,b, Pb, Kb, eb, ωb), and V0. Middle: the five Keplerian parameters of the candidate planet (T0,c, Pc, Kc, ec, ωc), and the RV jitter (σSOPHIE). Bottom: the GP hyperparameters of the stellar magnetic activity model (Prot, αQP, λ1,QP, λ2,QP), and the short-term stellar noise model (αSE,/λSE).

Appendix D List of the main physical parameters of the HD 73344 planetary system

Appendix E Posterior distribution of the stability coefficient

In Sec. 5.1, we quantify the dynamical stability of the HD 73344 system by using the stability coefficient δ_b). As noted by Stalport et al. (2022), an instructive insight of the data can be obtained by computing the distribution of δ_b that corresponds to the posterior distribution of the parameters. To this aim, we use the MCMC posterior sample coming from a joint fit (radial velocity and transit data; see Sec. 4.3), made of 10⁵ realizations of the system²¹. The histogram of δ_b drawn from these realizations is presented in Fig. E.1 for different values of the inclination of planet c. Histograms for values 30° < i_c < 88° are not shown; they all present a single peak at log₁₀ δ_b ≈ −5 or below. When we decrease i_c below 30°, however, Fig. E.1 shows that the distribution gradually transitions from a single stable population (for i_c ≳ 30°) to a single highly unstable population (for i_c ≲ 5°). In between, the population contains both stable, metastable, and unstable subsamples, visible as the multiple bumps in the histogram.

Fig. E.1 Posterior distribution of the stability coefficient δb coming from the joint fit (radial velocity and transit data). Different values for the inclination of planet c in the sky plane are assumed (see labels). The corresponding mass of planet c ranges from 0.7 MJ for panel a to 20 MJ for panel d.

Appendix F Inferring the bulk H₂O and H₂-He content using interior structure models

To infer the bulk composition of HD 73344b, we used the open source tool SMINT²² (Piaulet et al. 2021). This tool performs a MCMC retrieval on grids of interior structure models existing in the literature, to infer the composition of a planet based on its physical properties. We considered two possible compositions for the interior: i) an Earth-like core with a H₂-He envelope of solar metallicity (Lopez & Fortney 2014), and ii) a refractory core with a variable core mass fraction and a pure H₂O envelope and atmosphere on top (Aguichine et al. 2021). These compositions represent end-member cases between an envelope that would form with a Sun-like composition, and a H₂-He free envelope. Our goal was to determine the range of possible bulk volatile contents in HD 73344b.

For the H₂-He case, the MCMC takes as input Gaussian priors on planet mass, planet radius, incident stellar flux, and age, and produces as output the posterior for the envelope mass fraction f_env that best matches the radius. For the pure H₂O case, the inputs are planet mass, planet radius, and equilibrium temperature, and the outputs are the water mass fraction $f_{{\mathrm{H}}_2\mathrm{O}}$, and the composition of the core $f_{\text{core }}^{\mathrm{\prime }}$. We found that the properties of HD 73344b are compatible with an interior where f_env = 2.5 ± 0.3% or $f_{{\mathrm{H}}_2\mathrm{O}}={86}_{-10}^{+7}\mathrm{\%}$. The posteriors on all parameters for the two cases are shown in Figure F.1. We tested the case of an envelope with a metallicity 50 times solar (also from Lopez & Fortney 2014), and found a smaller value f_env,50 = 2.1 ± 0.4%. This is likely due to the fact that higher metallicity planets have greater atmospheric opacity, and therefore cool down (and contract) slower than solar metallicity planets. We also tested the pure water case where we fixed composition of the core to the Earth value $f_{\text{core }}^{\mathrm{\prime }}$ = 0.325. The results were extremely similar, owing to the fact that the core represents only ~20% of the planet mass (and, consequently, radius), so that its composition has a marginal impact on the total planet radius. We noticed that in all our cases, the posterior on the mass is centered at a mass of ~4 M_⊕ instead of the measured ~3 M_⊕. This is very likely due to the Gaussian prior on the mass, which would allow masses of < 1 M_⊕ at 1-σ (limit on the validity range of interior models), and even negative masses at 2-σ. Such values are discarded from the fit, favoring the higher-end distribution of masses and pushing the mean mass to a slightly greater value. This implies that the volatile content is slightly underestimated. Therefore, from this analysis we conclude that the volatile content of HD 73344b would be > 75% if it was water, and 2 − 3% if it was gas of solar composition.

Fig. F.1 Results of the SMINT tool on parameters of HD 73344b. Top panel corresponds to the H2-He envelope case (Lopez & Fortney 2014). Bottom panel corresponds to the pure H2O envelope case (Aguichine et al. 2021).

References

All Tables

All Figures

	Fig. 1 Light curves of HD 73344. The transit mid-times of planet b in K2 (top) and TESS (bottom) observations are shown with the dotted vertical lines.
In the text

	Fig. 2 Radial velocity of HD 73344 obtained over two consecutive nights with the SOPHIE spectrograph. Best fitting GP model from their joint analysis is shown in purple. Observations taken at airmass >1.7 are shown in red.
In the text

	Fig. 3 Generalized Lomb-Scargle periodograms of the K2 and TESS photometric data between periods of 5 and 40 days. Transit of planet b have been masked. The green vertical lines indicate the period of 9.1 ± 0.2 days found with RV data and chromospheric indicators.
In the text

	Fig. 4 Generalized Lomb-Scargle periodograms of K2, TESS, and SOPHIE RV data for periods of less than one day. The axes are presented in a log-scale format. The yellow and red dashed curves represent best Harvey-function fits to the periodograms, helping for a visual representation of both white Gaussian noise (WGN) and short-term stellar variability.
In the text

Fig. 5 Generalized Lomb-Scargle periodograms for RV and activity indicators. GLSP are calculated on the basis of binned SOPHIE data, where we have eliminated long-term variation using a two-degree polynomial fit. From top to bottom, we show the GLSP for the RV, FWHM, the Area of the Gaussian fit to the CCF, Bisector, $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$, and Hα lines. The red, green, and blue vertical lines indicate the stellar rotation period, the orbital period of the transiting planet (b), and the orbital period of the candidate planet (c), respectively.

In the text

Fig. 6 Normalized posterior distribution of the parameters fitted to the RV data. The distributions resulting from the analyses of the SOPHIE binned and unbinned observations are shown in red and black, respectively. The distributions resulting from the joint analysis combining the photometric (K2+TESS) and RV (SOPHIE+HIRES) observations are shown in yellow. Top: Five Keplerian parameters of planet b (T0,b, Pb, Kb, eb, ωb) and V0. Middle: Five Keplerian parameters of the candidate planet (T0,c, Pc, Kc, ec, ωc), and the RV jitter (σSOPHIE). Bottom: GP hyperparameters of the stellar magnetic activity model (Prot, αQP, λ1,QP, λ2,QP), and the short-term stellar noise model (αSE, λSE). For the latter, the values resulting from the analysis of the two SOPHIE nights of observation (Sect. 3.2.2) are shown for comparison (blue).

In the text

Fig. 7 Generalized Lomb-Scargle periodograms of the SOPHIE binned (left) and unbinned (right) data. From top to bottom: Raw data, then subtraction of the best-fitting Keplerian model for the planet candidate, the stellar activity model(s), and planet b. For the stellar activity model(s), we considered the mean of the predictive GP distribution resulting from the best fits. The stellar rotation period, and the orbital period of planets b and c are indicated by the red, green and blue vertical lines, respectively. In each panel, the red dots indicate the highest periodogram peak taken over periods >1 day. We note that the planetary and activity signals are jointly estimated from the RV data, and then iteratively removed in the successive panels.

In the text

Fig. 8 Radial velocity of HD 73344 resulting from the joint analysis of the photometric and RV data. Panel a: RV of SOPHIE (black) and HIRES (yellow) observations and best fitting model (purple). Panel b: RV residuals. Panel c: RV phased at the period of planet b (planet c and activity models subtracted). Panel d: RV phased at the period of planet c (planet b and activity models subtracted). Best-fitting models for planets b and c are shown in red.

In the text

	Fig. 9 Phased light curve of HD 73344b resulting from the joint analysis of the photometric and RV data. K2 data are shown in black, TESS data in blue. These light curves have been corrected by the best-fitting noise models. The best-fitting transit model is shown in red. Yellow dots represent the 120min binned light curves (K2 and TESS combined).
In the text

	Fig. 10 Spitzer observations of HD 73344. Top: raw Spitzer light curve (black points) of HD 73344b and best-fit models with (blue) and without (black) a transit fit. Middle: detrended Spitzer data and best-fit transit light curve. Bottom: residuals to the fit.
In the text

Fig. 11 Stability of the HD 73344 system as a function of the semimajor axis and eccentricity of planet c. The cross shows the best-fit values and their 1-σ uncertainty interval. The colour scale is chosen so that dark blue denotes stable orbits and red is highly unstable. The inclination chosen for planet c is labelled; its corresponding mass is Mc = 0.36 MJ for panel a and Mc = 2.05 MJ for panel b.

In the text

Fig. 12 Transit probability of planet HD 73344 b computed as the fraction of time its impact parameter is smaller than 1. The map is drawn as a function of the unknown longitude of node and inclination of planet c with respect to the sky plane. Other parameters are set to their best-fit values (see Table D.1). The transit probability (colour scale) is obtained from a 50 kyr numerical integration (see text). Points are coloured white if a planet is ejected before the end of the simulation.

In the text

	Fig. 13 Same as Fig. 12, but showing enlarged views. In the hatched band, planet c would be observed to transit the star; this region is therefore excluded.
In the text

Fig. 14 Mass–radius diagram of confirmed exoplanets with mass and radius precision better than 3σ and 10σ, respectively. Planetary parameters (black dots) are taken from the NASA Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu), and updated to October 2023. HD 73344b is represented by the red dot. The red arrow represents the upper 10σ limit on its mass. Composition models calculated from Aguichine et al. (2021) with a core mass fraction of 30% and various H2O mass fraction (from 10% to 80%) are represented by coloured lines (see text). Composition models calculated from Lopez & Fortney (2014) for a 1 Gyr-old system, Teq,b ~ 900 K, and an envelope mass fraction fenv = [1, 2, 3, 4, 5] % are shown in gray. For comparison, the relation for an Earth-like composition (32.5% core, 67.5% mantle) from Brugger et al. (2017) is shown by the green dashed curve, and some solar system planets (Venus, the Earth, Uranus and Neptune) are displayed.

In the text

	Fig. A.1 Light curves of HD 73344 showing the individual transits of planet b as seen with K2 (6 transits; top left panels) and TESS (3 transits; bottom left panels). The best-fitting models are shown in red. The corresponding residuals are shown in the right panels.
In the text

Fig. C.1 Temporal evolution of the stellar activity. From top to bottom: RV data, FWHM, the Area of the Gaussian fit to the CCF, Bisector, $\log R_{\mathrm{H}\mathrm{K}}^{\mathrm{\prime }}$ Hα lines, and linear trend. From left to right: Time series, GLSP for first campaign, second campaign and both campaigns. In all time series, we have eliminated long-term variation by a two-degree polynomial fit. The large red, green and blue vertical lines indicate the stellar rotation period (and half its period), the orbital period of the transiting planet (b), and the orbital period of the candidate planet (c), respectively.

In the text

Fig. C.2 Normalized posteriors distribution of the parameters fitted to the RV data. The distributions resulting from the analyses of the SOPHIE binned observations are shown in magenta (C1), pink (C2), and red (C1+C2). The distributions corresponding to the analyses of the unbinned data are shown in blue (C1), cyan (C2), and black (C1+C2). The distributions resulting from the joint analysis combining the photometric (K2+TESS) and nonbinned RV (SOPHIE+HIRES) observations are shown in yellow (benchmark). Top: the five Keplerian parameters of planet b (T0,b, Pb, Kb, eb, ωb), and V0. Middle: the five Keplerian parameters of the candidate planet (T0,c, Pc, Kc, ec, ωc), and the RV jitter (σSOPHIE). Bottom: the GP hyperparameters of the stellar magnetic activity model (Prot, αQP, λ1,QP, λ2,QP), and the short-term stellar noise model (αSE,/λSE).

In the text

	Fig. E.1 Posterior distribution of the stability coefficient δb coming from the joint fit (radial velocity and transit data). Different values for the inclination of planet c in the sky plane are assumed (see labels). The corresponding mass of planet c ranges from 0.7 MJ for panel a to 20 MJ for panel d.
In the text

	Fig. F.1 Results of the SMINT tool on parameters of HD 73344b. Top panel corresponds to the H2-He envelope case (Lopez & Fortney 2014). Bottom panel corresponds to the pure H2O envelope case (Aguichine et al. 2021).
In the text