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

There are several ways to test whether the 125 day 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 FAP (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.

We first isolated the 125 day signal by removing the contribution of the outer and transiting planets. The blue circles in Figure 3 show the growth of the LS power (P) of the 125 day 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 day RV variations (see below) sampled like the data and with the appropriate random noise added (σ = 1.5 m s⁻¹). 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 day 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 day period is stable and coherent.

Zoom In Zoom Out Reset image size

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 FWHM and bisector inverse slope (BIS) of the cross-correlation function, and 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 1 yr, 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 day⁻¹ (see Section 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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, prewhitening 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 day 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 Figure 6. The S-index and FWHM residuals do not show any significant ²⁸ peaks at 0.008 day⁻¹. The S-index only shows a weak one that is the fifth-highest peak in the frequency range. In the amplitude spectrum it appears to be less than twice the mean height of the surrounding noise peaks, an indication that it is not significant (see Kuschnig et al. 1997).

Zoom In Zoom Out Reset image size

The peak in the BIS residuals at 0.008 day⁻¹ is the highest peak in the frequency range considered (Figure 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.

4.2.1. The Bisector Variations

The bisector is the only activity indicator that shows possible variations at the 125 day 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 day period is planetary in nature, this merits a critical evaluation of the reality of these.

The periodogram of the raw bisector measurements (Figure 4, bottom panel) shows a peak coincident with the 125 day period (f = 0.008 day⁻¹). Without any a priori information, this has an 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 days) 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 unnormalized power, then FAP ∼ e^−z . In this case the FAP is rather low at FAP = 0.13%. This value is consistent with 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 day period results in an amplitude of 0.33 m s⁻¹ or 2.5 times smaller than the rms scatter (σ = 0.83 m s⁻¹) about the fit. In spite of this, the 125 day signal should have been detected at much higher significance owing to the large number of measurements. As a test, we took a 125 day sine function with an amplitude of 0.33 m s⁻¹ and added random noise with the rms scatter of the measurements. The periodogram of the simulated data showed power consistent with an 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 prewhitened analysis provides unconvincing support that the bisector variations are intrinsic to the star. Our analysis of the residuals showed an FAP ∼ 7%. Since we were focused on the 125 day signal, we initially only considered the frequency range out to 0.025 day⁻¹ (Figure 6, bottom panel). If we extend the analysis to higher frequencies, we find that the highest peak in the residuals occurs at 20 days (f = 0.05 day⁻¹). Removing this signal weakens the peak at 125 days further such that its FAP increases to ∼50%.

Third, we examined whether any phase variations of the bisectors can be related to the 125 day signal. Figure 7 shows the orbital fit to the 125 day RV variations. This fit is presented and discussed in detail 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 day RV signal. This strongly indicates that the signal found in the periodogram is due to noise, consistent with the large FAP we derived.

Zoom In Zoom Out Reset image size

Finally, it is worth noting that the case of 51 Peg 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 days 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.

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 day time span of S-index measurements showing sine-like variations with a period of 18.7 ± 0.8 days.

Zoom In Zoom Out Reset image size

Damasso et al. (2020) measured a projected rotational velocity of the star of v sin i = 3.34 ± 0.07 km s⁻¹, which agrees with the value of 3.3 ± 0.5 km s⁻¹ found by Gandolfi et al. (2018). Csizmadia et al. (2021, submitted) measured a stellar radius of R_⋆ = 1.190 ± 0.004 R_⊙, which results in a maximum rotation period of P_rot = 18.0 ± 0.4 days, consistent with the value from the S-index variations and that inferred by Damasso et al. (2020). Interestingly, the 20 day 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.

A preliminary Keplerian fit to the 125 day 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 data and used only the HARPS-POST and ESPRESSO data 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 Tables 4 and 6. They are defined as the median and 68% region of the credible interval of the posterior distributions for each fitted parameter.

Table 6. Orbital Parameters for π Men d

Parameter	Value
Orbital period Porb,d (days)	${124.64}_{-0.52}^{+0.48}$
Time of inferior conjunction T0,d (BJDTDB–2,450,000 days)	${7595.46}_{-6.39}^{+6.90}$
Eccentricity ed	0.220 ± 0.079
Argument of periastron of stellar orbit ω⋆,d (deg)	${323}_{-73}^{+25}$
Radial velocity semiamplitude variation Kd (m s−1)	1.68 ± 0.17
Planetary minimum mass Md × sin id (M⊕)	13.38 ± 1.35

Download table as:  ASCIITypeset image

π Men d has a minimum mass of M_d sini_d = 13.38 ± 1.35 M_⊕ and an orbit that 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 Figure 9, and the time series, focusing on just the time span of the HARPS-Large RVs, is shown in Figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We found that the difference between the Bayesian information criterion of the three-planet (b, c, and d) and two-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 have an rms scatter of σ_HARPS‐Post = 1.40 m s⁻¹, and the ESPRESSO data have σ_ESPRESSO = 0.95 m s⁻¹. Both are about a factor of three larger than the nominal measurement errors, which are typically better than 0.5 m s⁻¹.

For the 125 day RV variations to be due to a planet, it must lie on a stable orbit. A detailed dynamical investigation is beyond the scope of this paper. Instead, we performed a preliminary analysis to assess whether 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 three-planet orbital solutions, while the inclination of π Men b was drawn from i_b = 50° ± 5° according to previous studies (Damasso et al. 2020; De Rosa et al. 2020; Xuan & Wyatt 2020). Since we can only derive the minimum mass ²⁹ (M_d sin i_d) of π Men d, we allowed for an orbital inclination between 20° and 90°.

The left and middle panels of Figure 11 show 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 days 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.

Zoom In Zoom Out Reset image size

Due to limited computational resources, our simulations only covered a time span of 20 Myr, a small fraction of the estimated stellar age of ≈3.3 Gyr (Damasso et al. 2020). Clearly, a numerical simulation is called for, covering a much larger 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 Figure 11 shows that the planet's orbit has a higher chance to be more stable 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 yr were able to show that the actual orbit was nearly circular, e = ${0.07}_{-0.05}^{+0.06}$ (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 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 prewhitening 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 days (π 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⁻¹. An improved approach is to perform a joint fit of all Keplerian orbits that are present. This results in K_c = 1.21 ± 0.12 m s⁻¹, a value in agreement with the result from the prewhitening 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.

Zoom In Zoom Out Reset image size

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 1 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 phase 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 1 day. The only criteria are that the period of the transiting planet is 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 and 125 day periods present in the RV data. Although our inferred rotation period of 18 days 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 2−4 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 χ² 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 π 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⁻¹, a value consistent with the results from other methods.

Zoom In Zoom Out Reset image size

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 meters per second 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⁻¹.

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 and 5 m s⁻¹. The appropriate random noise was added, and to test an extreme case, additional random offsets of several kilometers per second 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 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 days and applied FCO. It yielded K_c = 1.40 ± 0.49 m s⁻¹, 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 (prewhitening, joint fit, and FCO). These have a mean value of 1.22 m s⁻¹ 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.

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 set samples 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⁻¹ is slightly lower than the value of 1.5 ± 0.2 m s⁻¹ found by Damasso et al. (2020). This may be due to authors not removing the underlying 125 day signal, which could not be seen in their data as a result of insufficient sampling. They did note the possible presence of a signal at ≈194 days, 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⁻¹. 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 K_c = 1.16 ± 0.34 m s⁻¹, which compares favorably to our value of 1.25 ± 0.24 m s⁻¹ using the ESPRESSO data. For bright targets where the S/N is high, HARPS should perform as well as ESPRESSO for RV measurements.

The K_c-amplitudes 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⁻¹, 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⁻¹, 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⁻³ 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⁻³. This only highlights that even though we have a formal error of 0.7% in the planet radius, the true uncertainty may be larger.

With a bulk density of 2.03 ± 0.22 g cm⁻³, π 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 by 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.

Zoom In Zoom Out Reset image size

Pi Men c has a relatively large radius, which means that it probably held a significant volatile envelope. It is thus a prime target for atmospheric characterization (Huber et al. 2022). García Muñoz et al. (2020) searched for hydrogen from photodissociation in the atmosphere of π Men c using Lyα in-transit spectroscopy with the Hubble Space Telescope (HST) but found none. One explanation for the nondetection was that the atmosphere was dominated by water or other heavy molecules rather than H₂/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 successfully 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 an H–He envelope from those without such an envelope. The orbital period of π Men c and its radius place it on the boundary between the two classes of planets (Figure 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., Fulton et al. 2017; Owen & Wu 2017; Fulton & 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).

Zoom In Zoom Out Reset image size

Clearly, the addition of just a single point in the period–radius diagram will not provide a breakthrough in understanding the origin of the radius valley. However, exoplanets with precisely determined masses and radii, in particular those in the middle of the gap like π Men c, are needed to shed more light on the cause of the gap.

Our analysis of the RV time series reveals the presence of a 125 day 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 timescales up to 20 Myr, at least over the orbital parameters that we have probed. In order to be certain that this RV signal comes from a planet, it must satisfy three criteria: (1) The signal should be statistically significant with an 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.) that would suggest a stellar origin for the variations. It appears that the 125 day signal satisfies these criteria.

The signal is highly significant with an FAP ≪ 3.3 × 10⁻⁶, 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 increase 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 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 day RV signal whose power behavior in the P–N diagram 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 seem to be no clear variations with the RV period in the activity indicators. The S-index shows periodic variations with a period (P ∼ 500 days) that is much longer than the RV value. The FWHM shows a dominant peak in the periodogram at ∼200 days, while the BIS shows a period of ∼760 days. There is a weak feature in the periodogram of the BIS measurements that is coincident with the RV signal, but as discussed earlier, we do not deem this as significant.

The subset of S-index variations and the inferred value from the radius and rotational velocity of the star indicate a stellar rotational period of approximately 18 days. Clearly, the 125 days is not due to rotational modulation. If an activity cycle is present, it most likely has a period of 500–700 days, 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 day cycle in the F6 star τ Boo. One may speculate that the 125 day RV period is roughly one-fourth of an activity cycle of ∼500 days, and we cannot exclude this with certainty. However, it is puzzling that the third harmonic ³⁰ dominates the RV, yet the "fundamental" period dominates the S-index measurements.

We should stress that although the 125 day 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 roughly a few meters per second. 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 timescales 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 day signal in the RV data.

The planet π Men d may have implications for the formation of the inner 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 protoplanetary disk and then migrated inward if other planets were present. However, one could envision a scenario where π Men c formed first, migrated inward, 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 outer planet, d, remained, but in a highly eccentric orbit. The outermost planet (b) can also perturb the innermost planet (c). Various effects of the giant planet on π Men c were investigated by Xuan & Wyatt (2020) and De Rosa et al. (2020). However, here is not the place to speculate on specific formation scenarios. That is best left for detailed theoretical modeling that can produce the observed configuration of the π Men planetary system. To do that requires knowing the full architecture of the planetary system. 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 relatively high mass of the outer planet and its inclined orbit.

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.

Zoom In Zoom Out Reset image size

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 shorter 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 that 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% higher 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 day signal. It is clear that the precise mass determinations of small planets will come at a very steep price even for bright targets.

Second, π Men seems to be a relatively inactive star, at least at the timescales of the transiting planet's orbital period. It is also a bright star, which results in high-S/N data acquired in short exposure times. In spite of this, the resulting rms scatter for the RV measurements is ≈1.5 m s⁻¹ 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⁻¹, 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 planets, which will severely stress the available telescope resources. Instruments on 2–3 m class telescopes that provide an RV measurement precision of 3–5 m s⁻¹ 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⁻¹. Bringing more 2–3 m class telescopes to the upcoming PLATO follow-up effort could play an important role in the success of the mission.