A Radial Velocity Study of the Planetary System of π Mensae: Improved Planet Parameters for π Mensae c and a Third Planet on a 125 Day Orbit

π Men hosts a transiting planet detected by the Transiting Exoplanet Survey Satellite space mission and an outer planet in a 5.7 yr orbit discovered by radial velocity (RV) surveys. We studied this system using new RV measurements taken with the HARPS spectrograph on ESO’s 3.6 m telescope, as well as archival data. We constrain the stellar RV semiamplitude due to the transiting planet, π Men c, as K c = 1.21 ± 0.12 m s−1, resulting in a planet mass of M c = 3.63 ± 0.38 M ⊕. A planet radius of R c = 2.145 ± 0.015 R ⊕ yields a bulk density of ρ c = 2.03 ± 0.22 g cm−3. The precisely determined density of this planet and the brightness of the host star make π Men c an excellent laboratory for internal structure and atmospheric characterization studies. Our HARPS RV measurements also reveal compelling evidence for a third body, π Men d, with a minimum mass M d sin i d = 13.38 ± 1.35 M ⊕ orbiting with a period of P orb,d = 125 days on an eccentric orbit (e d = 0.22). A simple dynamical analysis indicates that the orbit of π Men d is stable on timescales of at least 20 Myr. Given the mutual inclination between the outer gaseous giant and the inner rocky planet and the presence of a third body at 125 days, π Men is an important planetary system for dynamical and formation studies.


INTRODUCTION
The bright (V = 5.65; Table 1) G0 dwarf star π Men has been the target of exoplanet studies for over 20 years.Jones et al. (2002) reported long period (P orb,b ≈ 2100 d) radial velocity (RV) variations that were consistent with the presence of a sub-stellar companion (π Men b) with a minimum mass of ≈ 10 M Jup in a highly eccentric (e b ≈ 0.6) orbit.The star was later found by NASA's Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014) to host a small planet (R c ≈ 2 R ⊕ ) in a 6.27-d orbit (π Men c; Huang et al. 2018;Gandolfi et al. 2018).Due to the brightness of the host star, this planet is a prime candidate for atmospheric characterization studies (García Muñoz et al. 2020, 2021).
Accurate planetary masses and radii are important for exoplanet studies.True planet masses are needed for understanding the architecture of exoplanets and for dynamical studies.Accurate bulk densities are essential for constraining the internal composition of the planet.In the case of π Men c, the measured RV amplitude (and thus mass) was based on archival data taken with the HARPS spectrograph.The observing strategy consisted of sparse observations spread over a relatively long time and thus was not geared to detect small, short period planets.The corresponding error in the mass, ≈ 20 %, made it difficult to distinguish between interior models consisting mostly of water, or having a significant fraction (≈ 50 %) of silicates (Gandolfi et al. 2018).The π Men system has received heightened interest because astrometric studies have determined the orbital inclination of the outer planet.Xuan & Wyatt (2020), De Rosa et al. (2020), and Damasso et al. (2020) combined the long time series of RV measurements for this star along with Hipparcos and Gaia astrometric measurements to determine an orbital inclination of i b ≈ 50 • .Not only does this pin down the true planet mass as M b ≈ 14 M Jup , but more importantly it establishes that the inner and outer planetary orbits are misaligned.π Men joins υ Andromedae (McArthur et al. 2010) and Kepler-108 (Mills & Fabrycky 2017) as stars hosting planets with large mutually inclined orbits.The π Men system is thus important for constraining planet formation theories and studying the dynamical evolution of planetary systems.
In order to measure a more precise mass for π Men c, we included observations of the host star as part of our ESO HARPS large program of spectroscopic follow-up of transiting planet candidates found by TESS.The purpose of these RV measurements was twofold.First, we wished to improve on the error in the planet mass to constrain better compositional models.An accurate mass is also needed to estimate the atmospheric pressure scale height needed for planning observations to detect atmospheric features.Second, we wished to study the architecture of this system by searching for additional planetary companions.This is of especial importance given the recent discovery of the mutual inclination of the orbits between the inner and outer planets.π Men was also intensively observed by the newly commissioned ESPRESSO spectrograph on the ESO's VLT which yielded a K-amplitude of 1.5 ± 0.2 m s −1 for the inner transiting planet π Men c (Damasso et al. 2020).This result offered us a chance to compare the performance of two stateof-the-art spectrographs designed for precise RV work, one a venerable instrument, HARPS (Mayor et al. 2003), mounted on a 3.6-m telescope and in use for almost 20 years, and a more modern one, ESPRESSO, mounted on an 8.2-m telescope (Pepe et al. 2014).The star is bright, so high signal-tonoise ratio (S/N ) data can be obtained on both instruments with relatively short exposure times.In the future, considerable telescope resources will be invested in the RV follow-up of transiting small planets found by the PLATO mission (Rauer et al. 2014), so it is useful to compare the performance of both instruments using contemporaneous observations on the same target.
In this work, we adopted the most recent stellar parameters derived by Damasso et al. (2020).The only exception is the stellar radius, which comes from Csizmadia et al. (2021, submitted).For the sake of completeness, they are listed in Table 1, along with the equatorial coordinates, V-band magnitude, parallax, distance, and proper motion of the star.

THE RADIAL VELOCITY DATA
Our RV data consists of archival as well as new measurements from our ESO's HARPS large followup program.A total of 77 RV measurements come from the UCLES spectrograph mounted on the 3.9-m telescope of the Anglo Australian Telescope (AAT).These can be found in Butler et al. (2006) or in Gandolfi et al. (2018).R. Wittenmyer kindly provided us with the more recent measurements used by Huang et al. (2018).Archival HARPS data consisting of 145 RV measurements (51 nightly averaged) were also taken from the public archive of the European Southern Observatory (ESO).We also included 275 RV measurements (37 nightly averaged) taken with the ESPRESSO spectrograph (Damasso et al. 2020).
The new data for π Men were taken as part of our ESO observing programs 0101.C-0829, 1102.C-0923, and 106.21TJ.001(PI: Gandolfi) and during technical nights (program IDs 60.A-9700 and 60.A-9709) on HARPS.These consisted of 413 RV measurements 1 spanning September 2018 to December 2020 (hereafter "HARPS-Large" data set).The star is bright, so we typically took several observations per night.Exposure times were 150-300 secs resulting in a median signal-to-noise ratio (S/N) of ∼250 per pixel at 550 nm.If one considers only nightly averages our program resulted in 177 new measurements taken at different epochs.In June 2015 the HARPS fiber bundle was upgraded (Lo Curto et al. 2015).This results in a zero-point offset between the data taken before and after the upgrade.We therefore treated the complete data as four independent sets with different zero point offsets: UCLES, ESPRESSO and HARPS before (HARPS-PRE) and after the fiber upgrade (HARPS-POST, this set also includes HARPS-Large).
We extracted the HARPS spectra using the Data Reduction Software (DRS; Lovis & Pepe 2007).There are a number of reduction pipelines available for the calculation of RVs: The DRS, which uses the cross-correlation method with a digital mask (Pepe et al. 2002), the HARPS Template-Enhanced Radial velocity Re-analysis Application (HARPS-TERRA, Anglada-Escudé & Butler 2012) and the SpEctrum Radial Velocity AnaLyser (SERVAL) pipeline (Zechmeister et al. 2018).For our final analysis we used the RVs calculated with HARPS-TERRA as this produced a final root mean square (rms) scatter that was about 4 % lower than the other two methods.We emphasize, however, that all reduction programs produced consistent orbital parameters that were well within the uncertainties.
Briefly, HARPS-TERRA performs a least squares matching of each observed spectrum to a high signal-to-noise template spectrum produced by co-adding all the observations of the target star after they have all been placed on the same wavelength scale and corrected for the Earth's barycentric motion.The RV for each spectrum is calculated with respect to this master template.
π Men is a high proper motion star.The RV data span 17 years so it is important to remove the secular acceleration.The star has a systemic radial velocity of 10.595 ± 0.0003 km s −1 (as derived from the analysis of the HARPS DRS RVs), a parallax of 54.705 ± 0.067 mas, and proper motion of 311.187 ± 0.127 mas yr −1 and 1048.845± 0.136 mas yr −1 in RA and Declination, respectively (Gaia Collaboration et al. 2018, see also Table 1).This results in a secular acceleration of ∼0.48 m s −1 yr −1 .This was removed for all RVs except for those processed with HARPS-TERRA, which already accounts for the secular acceleration in its pipeline.Table 2 lists the data sets used in our RV analysis.We reiterate that the "HARPS-POST" data contains both the 17 archival HARPS measurements after the upgrade and those from our large program.The "HARPS-Large" is a subset of this which only contains our new 404 useful measurements from the large program.The standard deviation listed is the root mean square (rms) scatter of the data sets after the removal of all periodic signals (see below).Table 3 lists the new RV measurements from our ESO large programs.

PERIODOGRAM ANALYSIS OF THE RV DATA
We first performed a frequency (periodogram) analysis on the RV data in order to confirm that the signal of the transiting planet π Men c is indeed present in the data.Furthermore, identifying all significant signals in the data and modeling these is important for deriving a precise RV amplitude for the transiting planet.To find weak periodic signals we first had to remove the large variations due to the outer planet.Since the UCLES measurements increase the time base of the measurements needed to refine the parameters of the outer planet, we included these in spite of these having poorer RV precision.
An initial orbital solution was performed using the general non-linear least squares fitting program Gaussfit (Jefferys et al. 1987).All orbital parameters were allowed to vary, including the zero-point offset of each individual data set.The orbit parameters listed in Table 4 represent the final ones from our joint fit (see below), which are entirely consistent with the Gaussfit results.The orbital fit for π Men b is shown in Figure 1.The UCLES measurements have an rms scatter more than twice that of the HARPS data and will be excluded from the subsequent analyses.The generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009) of the residual HARPS-Large RV data, after removing the orbital frequency of planet b (f 1 = 4.79 × 10 −4 d −1 ), is shown in the top panel of Figure 2. We only used the HARPS-Large data for this analysis for two reasons.First, it provides a much simpler sampling window.Including the early HARPS-PRE data results in a very complex window with many more alias peaks.Second, the HARPS-PRE data start approximately 10 years earlier and have much sparser sampling.This means that any underlying long term stellar variability, which will be difficult to model, can boost power into an alias frequency thus masking the true one that is present.
The highest peak occurs at a frequency of f 2 = 0.008 d −1 (P = 125 d), although one can see significant power at the orbital frequency of the transiting planet (f 3 = 0.16 d −1 ; P = 6.27 d).Removing f 2 increases the power seen at orbital frequency of the transiting planet (middle panel of Figure 2).The final residuals (lower panel) results in no additional significant peaks.A peak is seen at f = 1.56 × 10 −3 (P = 640 d), but this has low significance with an estimated false alarm probability, FAP ≈ 2 %, which we do not consider significant.
We assessed the statistical significance of the new ∼125 d period via the bootstrap randomization process where the RV values were randomly shuffled keeping the time stamps fixed (Murdoch et al. 1993).In 300,000 realizations of the bootstrap there was no instance where the random data periodogram showed power higher than the real data.This implies a FAP 3.3 × 10 −6 .It is highly unlikely that this peak is due to random noise.5).

THE NATURE OF THE 125-D PERIOD
Before we can adequately model the 125-d signal in our analysis, it is important to establish its nature.If it is planetary in origin it would be an important new member of the π Men system; however, it can also arise from stellar activity.Important criteria for establishing the planetary nature of an RV signal are that it is long-lived and coherent and that it is not present in any activity indicators.

Coherence of the 125-d Period
There are several ways to test whether the 125-d period in the RV data is coherent and stable.One can examine the evolution of the power in the standard Lomb-Scargle (LS) periodogram as a function of the number of measurements (Hatzes 2016) or the evolution of the false alarm probability (Trifonov et al. 2018).Alternatively, one can use the Stacked-Bayesian GLS periodogram (Mortier & Collier Cameron 2017) to assess the stability of a signal.However, it can be difficult to assess the stability of a signal merely by looking at the stacked periodogram.The pathology of the sampling often can make a signal appear unstable for a time when in fact it is not, or vice versa.We therefore chose to use the growth of the LS power and compare it to what is expected from a simulated signal with the appropriate noise added.This provides a somewhat more "quantitative" assessment of the stability.5).The veritical dashed line is the frequency of the 125-d period seen in the RV data.
We first isolated the 125-d signal by removing the contribution of the outer and transiting planets.The blue dots in Figure 3 shows the growth of the LS power (P ) of the 125-d signal as a function of the number (N ) of RV measurements (we will refer to this as the P -N relationship) using the HARPS-Large and ESPRESSO data.We then compared this to expectations (red triangles) using an orbital fit to the 125-d RV variations (see below) sampled like the data and with the appropriate random noise added (σ = 1.5 m s −1 ).The error bars indicate the standard deviation in the power using simulated data with 10 different realizations of the noise.The evolution of the power of the 125-d period largely follows the simulated data in that there is a monotonic increase in the power as a function of number of data points, although the slope in the power evolution of the real data is slightly shallower than the simulated data.Overall, the behavior of the P -N relationship seems to indicate that the 125-d period is stable and coherent.

Activity Indicators
Stellar activity can also create periodic RV signals that can mimic a planet and this can be revealed by a periodogram analysis of activity indicators.If an RV signal is absent in all activity indicators then we can be more confident that it is the barycentric motion of the host star due to the presence of a companion.We have three activity indicators at our disposal: The full width at half maximum (FWHM) and bisector inverse slope (BIS) of the cross-correlation function (CCF), as well as the Ca II S-index measurement (see Baliunas et al. 1995).HARPS-TERRA also delivers indices on the Balmer Hα and Na D lines.However, a period analysis of these data showed a dominant period at one year, possibly an indication of telluric contamination.We therefore excluded these from our analysis.
Figure 4 shows the GLS periodograms for the three activity indicators extracted from HARPS-Large data.All have the highest peak at low frequencies that seem to be unrelated to the frequency found in the RV time series.The only exception is the BIS which shows a secondary peak near 0.008 d −1 (see Sect. 4.2.1 for a more detailed discussion of the bisector variations).Table 5 lists the dominant periodic signals found in the activity indicators and the FAP determined via 200,000 realizations of a bootstrap.Figure 5 shows the fit to the time series of each activity indicator using the period of the dominant peak found in the corresponding periodogram.
We investigated whether additional signals could be hidden in the activity indicators.Prewhitening, i.e., the subsequent fitting and removal of dominant frequencies in the periodogram of the time series and its residuals, is a powerful tool for finding weak, underlying signals in data.For instance, pre-whitening of the Hα index measurements for GL 581 was able to isolate variations with the same period as the purported planet GL 581 d (Hatzes 2016), thus supporting the claim of activity as the origin for the 66-d RV period (Robertson et al. 2014).
The GLS periodograms of the residuals of the activity indicators, after removing the appropriate sine fits to the dominant frequency in the data, are shown in Fig. 6.The S-index and FWHM residuals do not show any significant2 peaks at 0.008 d −1 .The S-index only shows a weak one that is the fifth highest peak in the frequency range.In the amplitude spectrum it appears less than twice the mean height of the surrounding noise peaks and indication that it is not significant (see Kuschnig et al. 1997).
The peak in the BIS residuals at 0.008 d −1 is the highest peak in the frequency range considered (Fig. 6).We assessed the significance of this peak via a bootstrap.The FAP that noise can produce a peak exactly coincident at this frequency is about 7 %, which we do not consider statistically significant.Note that a nearby peak has comparable power.

The Bisector Variations
The bisector is the only activity indicator that shows possible variations at the 125-d RV period.Both the raw and residual BIS variations show an isolated and modestly strong peak at the period (frequency) of interest.Since the bisector variations may be the only indicator to cast doubt on a planet hypothesis for the RV period and given the importance of establishing that the 125-d period is planetary in nature, this merits a critical evaluation of the reality of these.
The periodogram of the raw bisector measurements (Fig. 4, lower panel) shows a peak coincident with the 125-d period (f = 0.008 d −1 ).Without any a priori information this has a false alarm probability of FAP ∼ 30 %.However, this FAP represents the probability that noise can create a peak at least this high over a broad frequency range.As pointed out by Scargle (1982), for a known signal (i.e., 125 d) we need to consider the probability that noise can create a peak with the observed power or higher exactly at this frequency.If z represents the un-normalized power, then FAP ∼ e −z .
In this case the FAP is rather low at FAP = 0.13%.This value is consistent the FAP obtained with a bootstrap analysis.In spite of this low FAP we are confident that this signal is not significant for several reasons.
First, the bisector data are quite noisy.A sine fit to these using the 125-d period results in an amplitude of 0.33 m s −1 or 2.5 times smaller than the rms scatter (σ = 0.83 m s −1 ) about the fit.In spite of this, the 125-d signal should have been detected at much higher significance due to the large number of measurements.As a test, we took a 125-d sine function with an amplitude of 0.33 m s −1 and added random noise with the rms scatter of the measurements.The periodogram of the simulated data showed power consistent with a FAP that was more than 1000 times lower than was seen for the real data.A true signal should have shown a much higher power (i.e.lower FAP) than is observed.
Second, the pre-whitened analysis provides unconvincing support that the bisector variations are intrinsic to the star.Our analysis of the residuals showed a false alarm probability of FAP ∼7 %.Since we were focused on the 125-d signal we initially only considered the frequency range out to 0.025 d −1 (Fig. 6, lower panel).If we extend the analysis to higher frequencies we find that the highest peak in the residuals occurs at 20 d (f = 0.05 d −1 ).Removing this signal weakens the peak at 125-d further such that its FAP increases to ∼50 %.
Third, we examined whether any phase variations of the bisectors can be related to the 125-d signal.Figure 7 shows the orbital fit to the 125-d RV variations.This fit is presented and discussed in detail below in Section 5.1.The points show the binned (bin size ≈ 0.05 in phase) residual BIS measurements, after removing the dominant signal, phased to the same period.There is no hint of sinusoidal variations in the BIS that could account for the observed 125-d RV signal.This strongly indicates that the signal found in the periodogram is due to noise, consistent with the large FAP we derived.
Finally, it is worth noting the case of 51 Peg which highlights the pitfalls of relying only on bisector measurements to refute a planet hypothesis.Gray & Hatzes (1997) found bisector variations in 51 Peg that were coincident with the orbital period of the planet.The formal FAP for this signal was 0.3 % and thus seemed to be significant.Subsequently, higher quality bisector measurements were found to be constant (Hatzes et al. 1998).In spite of the low FAP, the previous variations were in fact due to noise.Since the bisector signal at 125 d in π Men is not supported by other activity indicators we conclude that these do not provide convincing enough evidence to refute a planet hypothesis for the RV signal.

The Rotational Period
In the full time series we found no significant peaks that may be indicative of stellar rotation.This is expected given that rotational modulation from activity may not be coherent over the long time span of our data.However, our time series did have long stretches when closely spaced measurements were made spanning 10-20 days.We searched for evidence of rotational modulation in these subsets and found only one instance, at the end of the S-index time series, which may show variations due to rotational modulation.Figure 8 shows a 16-d time span of S-index measurements showing sine-like variations with a period of 18.7 ± 0.8 d.Damasso et al. (2020) measured a projected rotational velocity of the star of v sin i = 3.34 ± 0.07 km s −1 which agrees with the value of 3.3 ± 0.5 km s −1 found by Gandolfi et al. (2018).Csizmadia et al. (2021, submitted) measured a stellar radius of as R = 1.190 ± 0.004 R , which results in a maximum rotation period of P rot = 18.0 ± 0.4 d, consistent with the value from the S-index variations and that inferred by Damasso et al. (2020).Interestingly, the 20-d period found in the bisector measurements is close to this value.Given the low significance of the bisector signal we are uncertain that it is truly related to stellar rotation.

KEPLERIAN MOTION FOR THE 125-D PERIOD
The most likely explanation for the 125-d period is that it stems from the presence of a third companion that we refer to as π Men d.Here we determine the orbital parameters and investigate the stability of the system.

Orbital Solution
A preliminary Keplerian fit to the 125-d period using the residuals after removing the contributions of the inner and outer planets indicates an eccentric orbit (e ∼ 0.2).We performed a joint fit for all three Keplerian signals using the code pyaneti (Barragán et al. 2019), which employs a Bayesian approach combined with Markov chain Monte Carlo sampling to estimate the planetary parameters.We derived the best-fitting orbital solution for π Men b from the UCLES, HARPS-PRE, HARPS-POST, and ESPRESSO date-sets, and used only the HARPS-POST and ESPRESSO date-sets to determine the parameters of π Men c and d.We imposed uniform uninformative priors for all the fitted parameters, except for the orbital period and time of first transit of π Men c, for which we adopted Gaussian priors from Damasso et al. (2020).A circular orbit for the inner transiting planet was assumed, but we fitted the eccentricity for π Men b and d.We accounted for the RV offsets between the different instruments/setups and fitted for jitter terms to account for both instrumental noise not included in the nominal uncertainties and RV variation induced by stellar activity.A parameter space with 500 Markov chains was explored to generate a posterior distribution of 250 000 independent points for each model parameter.The inferred parameters are given in Table 4 and  Table 6.They are defined as the median and 68 % region of the credible interval of the posterior distributions for each fitted parameter.
π Men d has a minimum mass of M d sini d = 13.38 ± 1.35 M ⊕ whose orbit is modestly eccentric (e = 0.220 ± 0.079).Interestingly, the angle of periastron passage, ω is comparable to that of the outer planet.The phase curve of the orbit is shown in Fig. 9 and the time series, focusing on just the time span of the HARPS-Large RVs, is shown in Fig. 10.
We found that the difference between the Bayesian information criterion of the 3-planet (b, c, and d) and 2-planet (b and c) models is ∆BIC = −73, providing very strong evidence for the existence of a third Doppler signal in the data.After removing the contribution of the orbital motions of all planetary signals the HARPS data has an rms scatter of σ HARPS−Post = 1.40 m s −1 and ESPRESSO, σ ESPRESSO = 0.95 m s −1 .Both are about a factor of three larger than the nominal measurement errors which are typically better than 0.5 m s −1 .

Orbital Stability
For the 125-d RV variations to be due to a planet, it must lie on a stable orbit.It is beyond the scope of this paper for a detailed dynamical investigation.Instead, we performed a preliminary analysis to assess if a stable orbit is at least feasible.For this dynamical study we employed Rebound (Rein & Liu 2012) and drew samples from the posterior distributions of our 3-planet orbital solutions while the inclination of π Men b was drawn from i b = 50 ± 5 • according to previous studies previous The left and center panels of Figure 11 shows the stability probability as a function of the eccentricity, period, and mass of planet d for our simulation.For the period range of 110-130 d the orbit of planet d is stable, although there are isolated regions where it is unstable.The planet must have a mass < 20 M ⊕ for stability, For our orbital solution (Table 6) this implies an orbital inclination i d > 40 • .The right panel of Figure 11 shows the stability probability in the mass versus mutual inclination plane.In general orbits are more likely to be unstable for high mutual inclinations and a low mass planet, or high mass even for relatively low mutual inclinations.
Due to limited computational resources our simulations only covered a time span of 20 Myrs, a small fraction of the estimated stellar age of ≈ 3.3 Gyrs (Damasso et al. 2020).Clearly, a numerical simulation is called for covering a much large fraction of the stellar life.For such a simulation covering a longer time span it would be important to obtain more accurate orbital parameters in order to restrict the parameter space of the simulations.
The left panel of Fig. 11 shows that the planet's orbit has a higher chance of stability if it lies on a more circular orbit, e < 0.3 and our measured eccentricity of 0.22 lies below this.Further RV measurements may reveal that the orbital eccentricity may in fact be lower.For instance, Eri b is an exoplanet where the K-amplitude is also comparable to the rms scatter of the RV measurements as is the case here.The discovery paper initially reported a high eccentricity (e = 0.6) for this planet (Hatzes et al. 2000).However, additional measurements spanning 30 years were able to show that the actual orbit was nearly circular, e = 0.07 +0.06 −0.05 (Mawet et al. 2019).The initial high eccentricity most likely stemmed from the low K-amplitude that was comparable to the measurement error and the influence of activity on the RV variations from the orbital motion.

THE K-AMPLITUDE OF π MEN C
Having determined the presence of all stable periodic signals in the RV time series we now focus on measuring a precise mass for the transiting planet π Men c.This hinges on the radial velocity K-amplitude, K c .

Pre-whitening and Joint Fitting
The simplest procedure is to fit sequentially a Keplerian orbit to each periodic signal, remove it, and fit the next periodic signal in a so-called pre-whitening procedure.We did this using only the HARPS-POST and ESPRESSO RVs, as these have the lowest rms scatter in our data-set.We first fitted the orbit of the outer planet (π Men b), removed this and then fitted the orbit of the planet at 125-d (π Men d).A fit was then made for planet c on the residual RVs.This procedure results in K c = 1.30 ± 0.13 m s −1 .An improved approach is to perform a joint fit an all Keplerian orbits that are present.This results in K c = 1.21 ± 0.12 m s −1 , a value in agreement with the result from the pre-whitening procedure to well within the nominal uncertainty.We adopt this as our best solution.
Figure 12 shows the phase-folded Doppler reflex motion induced by π Men c on the star, following the subtraction of the RV signals of the other two planets.The zero phase of the RV curve corresponds to the time of first transit (see Table 6).

Floating Chunk Offset Method
As an independent approach to determining the K-amplitude we applied the so-called floating chunk offset (FCO) method.This technique was developed in order to extract the Doppler reflex motion induced by ultra-short period (USP) planets, i.e., planets with orbital periods shorter than one day (Hatzes et al. 2010).The basis is that if one takes several measurements in one night the dominant variations come from the stellar orbital motion induced by the short-period planet.If the RV variations from other sources (rotation, long period planets, systematic errors, etc.) are much longer than the short-period planet, then the nightly variations stem predominantly from orbital motion of the USP planet and all other phenomena merely add a constant value to the RV for that night.One fits a Keplerian orbit keeping the period and phased fixed, but varying the K-amplitude and the zero-point offset velocity for each night (i.e.chunk).
The FCO method can also be applied to planets with orbital periods longer than one day.The only criteria are that period of the transiting planet be smaller than periods from other sources and the observing cadence is reasonably high.This is the case for π Men c whose orbital period is much less than the 2088-d and 125-d periods present in the RV data.Although our inferred rotation period of 18 d and its harmonics could have an influence, there is no evidence for these in the periodogram of the RV time series.Our observing strategy for π Men was such that RV measurements were taken on several consecutive nights.The advantage of the FCO method is that we do not have to remove the large orbital motion of the outer planet as it is naturally filtered out in the analysis.We divided the RV data from the HARPS large program into time chunks with consecutive measurements spanning two to four days.Figure 13 shows the FCO periodogram produced by fitting the RV chunk data using a trial period, finding the best amplitude and plotting the reduced χ 2 as a function of input period.The best fit was obtained for the period of the transiting planet; the FCO method can detect the signal of the π Men c.We note that the raw RV data with the large orbital motion of the outer planet as well as the respective instrumental offsets were present in the chunks, so it is relatively insensitive to the details in how the other signals are removed.This resulted in K c = 1.29 ± 0.22 m s −1 , a value consistent with the results from other methods.
A more refined value to the FCO result can be obtained by first subtracting the large orbital RV variations due to the outer planet (π Men b).The high amplitude and eccentricity can still result in a large, short-term variation, especially since our measurements caught the maximum in the orbital curve where the RV changes by many m s −1 over several days.This can introduce a systematic error in the K-amplitude determination.To remove the variations of π Men b we simply subtracted the orbital solution from Table 4 from the raw data, which retained all the zero point offsets of the individual data sets.For this final step we included the ESPRESSO measurements in the analysis.This resulted in K c = 1.16 ± 0.13 m s −1 .
As a test to ensure that FCO recovers known input signals we created synthetic simulated data consisting of the Keplerian orbits of planets c and d.These were sampled and divided into chunks in the same manner as the data.The input K-amplitude of π Men c was varied between 0 -5 m s −1 .The appropriate random noise was added, and to test an extreme case, additional random offsets of several km s −1 were added to each chunk.Over the full range of the considered K-amplitudes the FCO method was able to recover the input amplitude.
One advantage of the FCO method is that, unlike for other methods, it can produce good results on sparse data where you do not have a good knowledge of timescales of other signals in the data.As a test, we took a subset of only 27 measurements taken over a time span of 3400 d and applied FCO.It yielded K c = 1.40 ± 0.49 m s −1 , a ∼3 σ result that is consistent to within the error of our nominal value.Table 7 lists the K c -amplitude determinations using the different methods (pre-whitening, joint fit, and FCO).These have a mean value of 1.22 m s −1 and agree well within their nominal uncertainties, which indicates that the K c amplitude is robust and insensitive to the method one uses to extract it.

K-amplitude from Individual Data Sets
It is instructive to compare the K-amplitude of the transiting planet π Men c derived from individual and combinations of the data sets.In particular, single instrument sets will have the same zero point offset and systematic errors.We can see how consistent these solutions and their errors are to the nominal values.
For this exercise we first used the complete data set to fit and remove the orbital motion of the two outer planets.This was done because each data sets sample different parts and a smaller fraction of the longer period orbits and we want a fair comparison.Using subsets of the data to fit all planetary orbits will have a strong influence on K c .The remaining residuals contained just the orbital variations of the inner planet.
Table 8 shows the resulting K c amplitudes for the individual data sets.All have consistent values with each other and to our nominal value.The HARPS-PRE data, however, give a larger amplitude and error compared to the ESPRESSO result in spite of having a comparable number of measurements.It is not known whether this is due to a higher intrinsic stellar jitter during this time, or to the different quality of the optical fiber used for scrambling.Note that the ESPRESSO amplitude of K c = 1.25 ± 0.24 m s −1 is slightly lower than the value of 1.5 ± 0.2 m s −1 found by Damasso et al. (2020).This may be due to authors not removing the underlying 125-d signal, which could not be seen in their data due to insufficient sampling.They did note the possible presence of a signal at ≈ 194 d, but this was not significant and they chose not to include it in the modeling.
Adding the ESPRESSO data to the HARPS-POST data does not alter the K-amplitude significantly.In spite of the larger scatter of the HARPS-PRE data these also have little influence on the result when combining all the data sets.Although not listed in the table, we also performed a solution including the lower quality UCLES data and this resulted in K c = 1.24 ± 0.12 m s −1 .Clearly, the solution is driven by the higher quality HARPS and ESPRESSO measurements.
It is of interest to make a direct comparison between the performance of HARPS and ESPRESSO for comparable data taken at roughly the same time.To do this we took a subset of our HARPS measurements so that they had the same number as those from ESPRESSO (37 nightly averages) and roughly covered the same time span.This resulted in a K c = 1.16 ± 0.34 m s −1 , which compares favorably to our value of 1.25 ± 0.24 m s −1 using the ESPRESSO data.For bright targets where the S/N ratio is high, HARPS should perform as well as ESPRESSO for RV measurements.

The mass of π Men c
The K c -amplitude from the various methods and data sets (Tables 7 and 8, excluding the UCLES and HARPS-PRE RV measurements) are all consistent and tightly clustered around a mean value of 1.21 m s −1 , indicating a well constrained K c -amplitude and thus planet mass.Since all the RV measurements were taken over a relatively short time span with excellent sampling we take the K c -amplitude determined using a joint fit to the HARPS-Post + ESPRESSO data sets, namely K c = 1.21 ± 0.12 m s −1 , as our adopted value.Using the stellar mass of M = 1.07 ± 0.04 M from Damasso et al. (2020) results in a planet mass, M c = 3.63 ± 0.38 M ⊕ .Damasso et al. (2020) derived a ratio of the planet to star radii of R c /R = 0.0165 ± 0.0001 and a stellar radius of R = 1.17 ± 0.02 R .Recently Csizmadia et al. (2021;submitted) showed that the color indices of a star can be used to check stellar parameters.For π Men they determined a stellar radius of R = 1.190 ± 0.004 R .This results in a planet radius of R c = 2.145 ± 0.015 R ⊕ and a bulk density of ρ c = 2.03 ± 0.22 g cm −3 for π Men c.
We note that Huang et al. (2018) and Gandolfi et al. (2018) derived a planet radius of R c = 2.21 ± 0.04 R ⊕ and R c = 2.23 ± 0.04 R ⊕ , respectively.This results in a lower planet density of ρ c = 1.8 ± 0.2 g cm −3 .This only highlights that even though we have formal error of 0.7 % in the planet radius, the true uncertainty may be larger.

DISCUSSION
Our investigation of the RV variations of π Men has produced two main results -an accurate and precise mass for the transiting planet π Men c and the discovery of a third planet in the system.

π Men c in the Mass-Radius Diagram
With a bulk density of 2.03 ± 0.22 g cm −3 , π Men c appears to have a rather low density compared to most exoplanets with its mass.This is highlighted by Figure 14 which shows the location of π Men c in the mass-radius diagram for exoplanets with well determined radii and masses (better to 15 %).π Men c lies closest to the internal structure model for pure water.We stress that these are simple planet composition models for a planet without an atmosphere and, as noted Lopez & Fortney (2014), an atmospheric envelope may only account for a few percent of the planet mass, but it can have a huge impact on the measured planet radius.
With its relatively large radius π Men c means it probably held a significant volatile envelope and thus is a prime target for atmospheric characterization (Huber et al. 2021).García Muñoz et al. (2020) searched for hydrogen from photodissociation in the atmosphere of pi Men c using Ly-α intransit spectroscopy with the Hubble Space Telescope (HST) but found none.One explanation for the non-detection was that the atmosphere was dominated by water or other heavy molecules rather than H 2 /He.Subsequently, García Muñoz et al. ( 2021), again using HST, were able to detect C II ions at the 3.6-σ level which seems to support this scenario.The planet may still be in transition into a bare rocky planet.More progress on the characterization of the atmosphere of π-Men c will surely come with investigations using the successfuly launched James Webb Space Telescope.
The planet lies near the middle of the so-called radius valley, a gap in the radius distribution of small planets around 1.75 -2.00 R ⊕ that separates planets with masses sufficient to maintain a H-He envelope from those without such an envelope.The orbital period of π Men c and its radius places it on the boundary between the two classes of planets (Fig. 15) It is still an open question as to whether this gap results from planet formation or evolution.One hypothesis is that it is caused by atmospheric erosion of short-period planets due to photoevaporation from the close proximity to the host star (e.g.Owen & Wu 2017;Fulton et al. 2017Fulton & Petigura 2018).Alternatively, the gap may simply result from the planet formation process itself.Ginzburg et al. (2018) showed that planet formation with a core-powered mass loss mechanism could account for the radius distribution of planets without invoking photoevaporation.This scenario is able to match the radius valley, location, shape and slope (Gupta & Schlichting 2019).
Clearly, the addition of just a single point in the period-radius diagram will not provide the breakthrough in understanding the origin of radius valley.However, exoplanets with precisely determined masses and radii, in particular those in the middle of the gap like π Men c, are needed shed more light on the cause of the gap.

π Men d
Our analysis of the RV time series reveals the presence of a 125-d periodic signal in the data very likely due to a third planet in the system.Orbital solutions yield a minimum planet mass of M d sin i d = 13.38 ± 1.35 M ⊕ in an eccentric orbit (e = 0.220 ± 0.079).A preliminary dynamical study indicates that its orbit is stable on time scales up to 20 Myrs, at least over the orbital parameters that we have probed.In order to be certain this RV signal comes from a planet it must satisfy three criteria: 1) The signal should be statistical significant with a FAP < 0.1%.2) The signal should be long-lived and coherent with no change in the period, amplitude and phase.3) There should be no variations with the RV with any indicators (photometry, Ca II, etc.), which would suggest a stellar origin for the variations.It appears that the 125-d signal satisfies these criteria.
The signal is highly significant with a FAP 3.3 × 10 −6 , which makes it unlikely to arise from noise.It also appears to be stable and coherent.Adding measurements causes the statistical significance as shown by the P -N , to increases in the expected manner given the signal, sampling, and noise level.There may be a slight concern that the slope of the P -N behavior is shallower than expected from simulations, but this may be due to the noise not following a strictly a Gaussian distribution.We stress, however, that such a stable growth in P -N should not be used as a final confirmation of the planetary nature of an RV signal.A case in point is the proposed planet around α Tau.Hatzes et al. (2015) found evidence for a 629-d RV signal whose power behavior in the P -N diagran followed the expected growth for a stable signal.In spite of this, additional RV measurements seemed to contradict the planet hypothesis (Reichert et al. 2019).
There seems to be no clear variations with the RV period in the activity indicators.The S-index shows periodic variations with a period (P ∼ 500 d) which is much longer than RV value.The FWHM shows a dominant peak in the periodogram at ∼200 d, while the BIS shows a period of ∼ 760 d.There is a weak feature in the periodogram of the BIS measurements is coincident with the RV signal, but as discussed earlier, we do not deem this as significant.
The subset of S-index variations as well as the inferred value from the radius and rotational velocity of the star indicate a stellar rotational period of approximately 18 d.Clearly, the 125-d is not due to rotational modulation.If an activity cycle is present it most likely has a period of 500-700 d, as manifested in the activity indicators.G-type stars are expected to have activity cycles with timescales of years to decades.However, shorter period activity cycles are not unprecedented for other late-type stars.Schmitt & Mittag (2017) found evidence for a ∼120-d cycle in the F6 star τ Boo.One may speculate that the 125-d RV period is roughly one-fourth of an activity cycle of ∼ 500 d and we cannot exclude this with certainty.However, it is puzzling that the third harmonic4 dominates the RV, yet the "fundamental" period dominates the S-index measurements.
We should stress that although the 125-d signal seems to satisfy our criteria for planet confirmation, all of these criteria are necessary, but not sufficient conditions for planet confirmation.That is to say that an RV planet candidate must satisfy these criteria, but that is still no guarantee that the signal stems from the orbital motion of a planet.This is especially true for weak signals with RV amplitudes of a few ∼ m s −1 .Conceivably, a stellar activity cycle or its harmonics could be seen in the RV data without a strong presence in classic activity indicators.We are just starting to understand the influence of stellar activity on precise RV measurements and the time scales involved, so surprises could be in store.However, all the best available evidence at hand suggests that the presence of a third planet in the system is the most likely explanation for the 125-d signal in the RV data.
The planet π Men d may have implications for the formation of the innner planet.De Rosa et al. (2020) argued that the presence of nearby planets to π Men c would favor in-situ formation.At first glance, it would seem difficult for π Men c to have formed in the outer region of the proto-planetary disk and then migrated inwards if other planets were present.However, one could envision a scenario where π Men c formed first, migrated inwards and then π Men d formed at a later time.Dissipation of the disk would then halt the migration of π Men d.
Alternatively, the inner planet could have formed via tidal migration.An interaction of two planets (say c and d) would have scattered planet c into the inner regions where its orbit would have been circularized via tidal effects.The planet d remained, but in a highly eccentric orbit.The outermost giant planet can also perturb the inner planets.Various effects of the giant planet on π Men c were investigated by Xuan &Wyatt (2020) andDe Rosa et al. (2020).However, here is not the place to speculate on specific formation scenerios.That is best left for detailed theoretical modeling which can produce the observed configuration of the π Men planetary system.To do that requires knowing the full architecture of the planetary system and π Men b is a key component.It is therefore important to derive more accurate orbital elements for π Men d.Our simple dynamical analysis of the planetary system to π Men indicates that the orbit of planet d should be stable.Clearly, a more detailed dynamical analysis is required, which is best left for a dedicated investigation.The π Men planetary system offers us a very interesting system for such studies, especially given the relative high mass of the outer planet and its inclined orbit.

Lessons learned for the RV follow-up of small transiting planets
As a closing remark, our study of π Men offers us important lessons for the follow-up of small transiting planets in the era of TESS, and in the near future, PLATO.First, π Men c has one of the most precise mass determinations for a small planet, with an error of about 10 %.This precision required approximately 200 (nightly averaged) measurements taken with superb instruments on a bright star.If one wants to increase the precision on the mass measurement one naturally requires more measurements, but the number of these may not always follow the expectations from white noise.
Figure 16 shows the percent error in the mass of π Men c as a function of the number of RV measurements, N , using the homogeneous data set provided by the HARPS-Large RVs.For white noise one expects an error proportional to N −0.5 , but in reality it is slightly worse, being ∝ N −0.36 .At the beginning our RV measurements performed as expected, but at about N = 70 the error in the mass measurement increased.For subsequent measurements, the error does follow the expected behavior σ ∝ N −0.5 , but from a higher starting point.There must be a "red" noise component, most likely attributable to a variable contribution of stellar jitter.One cannot always rely on a few more RV measurements to improve the error on the K-amplitude.
In the case of π Men c, improving the mass measurement to better than 5 % would require more than 1000 RV measurements if one were to adhere to the same observing strategy and sampling of our study.This is 50 % more than the number estimated based on the trend shown at the start of the measurements.Clearly, if one is only interested in deriving the mass of the transiting planet it is more effective to obtain a large number of RV measurements for π Men over a much short time span, preferably over a single orbital period.Even then, a considerable number of observations will be required to find and remove additional signals due to stellar rotation or other planets which also can have a large influence on the K-amplitude of the transiting planet.This is demonstrated by the ESPRESSO data.Damasso et al. (2020) found a K c -amplitude 20% lower than our value for the same data set, most likely by not including the Keplerian motion of the third planet.Considerably more RV measurements were required to be certain of the presence of the 125-d signal.It is clear that the precise mass determinations of small planets will come with at a very steep price even for bright targets.
Second, π Men seems to be a relatively inactive star, at least at the time scales of the transiting planet's orbital period.It is also a bright star which results in high signal-to-noise ratio data acquired in short exposure times.In spite of this, the resulting rms scatter for the RV measurements is ≈1.5 m s −1 even after removing all periodic signals.This is a factor of 2-3 larger than the measurement errors, even when using the premier instruments in the world for RV measurements.Regardless of the effort taken to minimize instrumental and photon noise, the noise floor will be set by the intrinsic variability of the star.We will be hard pressed to find stars that are "RV quieter" than about 1 m s −1 and most stars will have an intrinsic jitter higher than this.
For bright targets with high stellar jitter, telescope size will not necessarily matter.The RV confirmation of small transiting planets will require an inordinate amount of telescope resources.π Men has shown us that for bright targets with high intrinsic jitter, using a larger aperture telescope does not always result in higher precision measurements.PLATO is expected to produce a large number of small transiting planet which will severely stress the available telescope resources.Instruments on 2-3m class telescopes that provide an RV measurement precision of 3-5 m s −1 can clearly play an important role in the PLATO era.The lack in measurement precision can be compensated with the shear number of measurements that can be invested on one target.A simple simulation shows that with ∼150 RV measurements over nine consecutive nights one can recover the K-amplitude of π Men c (4σ result) with a modest measurement precision of 3 m s −1 .Bringing more 2-3m class telescopes to the upcoming PLATO follow-up effort could play an important role in the success of the mission.

Figure 1 .
Figure 1.(Top) The RV time series from UCLES (brown squares), HARPS-PRE (red dots) and HARPS-POST (blue triangles).For clarity we do not plot the ESPRESSO measurements as these are contemporaneous with the HARPS-POST data.The curve is the Keplerian orbital solution for π Men b. (Bottom) The residuals after subtracting the orbit of π Men b.

Figure 2 .
Figure 2. (Top) The GLS periodogram of the HARPS-Large RV measurements after removing the motion of the outer planet (orbital frequency, f 1 = 4.79 × 10 −4 d −1 ).(Middle) The GLS periodogram of the RV residuals after removing the contribution of the dominant peak at f 2 = 0.008 d −1 .(Bottom) The GLS periodogram of the final RV residuals after removing the orbital frequency of the transiting planet, f 3 = 0.16 d −1 .

Figure 3 .
Figure3.Growth of the Lomb-Scargle power of the 125-d period attributed to a third body in the system as a function of the number of RV measurements for the real data (HARPS-Large and ESPRESSO, blue dots) and simulations (red triangles).

Figure 4 .
Figure 4.The GLS periodograms of the activity indicators extracted from the HARPS-Large data-set.From top to bottom: S-index, FWHM, and BIS.The vertical dashed line marks the orbital frequency of the 125-d detected in the RVs.

Figure 5 .
Figure5.The time series of the S-index (top), FWHM (middle) and BIS (bottom) measurements extracted from the HARPS-Large data-set.The red curves are sine fits using the period of the dominant peak in each periodogram (Table5).

Figure 6 .
Figure 6.GLS periodogram of the activity indicators after removal of the dominant signal in each (Table5).The veritical dashed line is the frequency of the 125-d period seen in the RV data.

Figure 7 .
Figure 7.The phase-binned BIS measurements using the ≈ 125-d period found in the RVs.The red curve represents an orbital solution to the 125-d period (see Section 5.1).

Figure 8 .
Figure 8. S-index measurements over 16 consecutive nights.The curve represents the sine fit with a period of 18.7 d.

Figure 9 .
Figure 9. RV variations and Keplerian orbit for π Men d phased to the orbital period of ∼125 d.The contributions from the inner and outer planets have been removed.Blue circles and red diamonds are nightly binned HARPS-POST and ESPRESSO RV measurements, respectively.The gray error bars include jitter.The curve represents the orbital solution.The zero phase of the RV curve corresponds to the time of inferior conjunction (see Table6).

Figure 10 .
Figure 10.The RV versus time from HARPS large program (blue circles) and ESPRESSO (red diamonds) nightly binned measurements for π Men d (125-d period) after removing the contribution of the inner and outer planets.The error bars include the jitter term.The curve represents the orbital solution.

Figure 11 .
Figure 11.Stability maps for π Men d.Yellow and green regions indicate regions with a high probability of stability whereas dark regions are unstable.(Left) Stability regions in the eccentricity -planet mass plane.(Right) Stability regions in the of orbital period -planet mass plane.

Figure 12 .
Figure 12.Phase-folded RV variations and Keplerian orbit for π Men c.The contributions from planets b and d have been removed.Blue circles and red diamonds are nightly binned HARPS-POST and ESPRESSO RV measurements, respectively.The gray error bars include jitter.The curve represents the orbital solution.The zero phase of the RV curve corresponds to the time of first transit (see Table6).

Figure 13 .
Figure 13.The FCO periodogram of the raw RV data.The red dashed vertical line marks the orbital period of the transiting planet.

Figure 14 .
Figure 14.The location of π Men c in the mass-radius diagram for exoplanets with well determined radii and masses (better to 15 %).Composition models from Zeng et al. (2016) are displayed with different lines and colors.

Figure 15 .
Figure 15.The location of the radius valley from Van Eylen et al. (2018).The blue dashed line marks the hyperplane of maximum separation.The location of π Men c is marked by the red square.

Figure 16 .
Figure 16.The error in the mass determination as a function of the number of measurements.The solid line represents a fit with error ∝ N −0.36 .The red dashed lines represent a nominal error ∝ N −0.50 expected for white noise.The top red dashed line has its origin at the N = 70 data point.

Table 1 .
Equatorial coordinates, proper motion, parallax, distance, V -band magnitude, and fundamental parameters of π Men.

Table 2 .
The RV Data Sets.

Table 3 .
RV measurements and activity indicators from the ESO HARPS Large Programs (full table in electronic version) BJD TDB RV σ RV BIS FWHM S-index σ S−index (d) m s −1 m s −1 m s −1 km s −1

Table 5 .
Periods in Activity Indicators and Their False Alarm Probabilities.

Table 8 .
K c -Amplitude from Data Sets