© The Authors 2025

1 Introduction

One of the greatest achievements in astrophysics in recent times is the detection of over 5000 planets orbiting stars beyond our Solar System. Since the first detection of exoplanets (Latham et al. 1989; Mayor & Queloz 1995), the vast majority of planets have been discovered using the radial velocity (RV) and transit techniques. Of particular interest is the characterisation of the atmospheres and internal composition of rocky-Earth mass planets orbiting in the liquid water habitable zones of their host stars. In the Habitable Worlds Catalog (HWC)¹, however, there are only 23 detections of potential rocky (0.5 R_⊙ < R_p ≤ 1.6 R_⊙ or 0.1 M_⊙ < M_p,min ≤ 3 M_⊙) planets that orbit in their host star’s conservative habitable zone (Kopparapu et al. 2014). Among these, only the planets in the TRAPPIST-1 system have been robustly determined to be rocky from accurate mass and radius measurements. Nine out of the 23 potentially habitable worlds were detected by the RV method, and the other 14 by the transit method.

The lack of robustly determined planetary mass values is primarily because the RV technique can only provide a lower limit on a planet’s mass unless the inclination is known from other measurements, such as via planetary transits. However, the probability that a planet transits is low in general. A planet in the habitable zone of a late M dwarf such as GJ 1061 has a transit probability of about 1%. Moreover, the majority of known transiting systems are at larger distances to our Sun, with a relatively faint host star apparent brightnesses, making characterisation of their atmospheres more challenging. The transit planets in the HWC for example have a mean distance of 105 pc while the RV planets have a mean distance of only 3.4 pc. The latter are therefore prime targets for direct spectroscopic investigations of planet atmospheres in the era of extremely large telescopes with instruments such as the ArmazoNes high Dispersion Echelle Spectrograph (Palle et al. 2023) or the Planetary Camera Spectrograph (Kasper et al. 2021) and with space missions such as the Habitable Worlds Observatory (National Academies of Sciences, Engineering and Medicine 2023) or the complementary Large Interferometer For Exoplanets mission (Quanz et al. 2022a,b).

In multi-planet systems, the strength of the gravitational planet-planet interactions depends on the masses and mutual inclinations of the planets, and opens up the possibility to obtaining more precise planetary parameters. For dynamically compact systems, even small deviations from Keplerian orbits can become detectable within just a few years for close orbiting planets, especially with the advent of ultra-high precision spectrographs such as ESPRESSO (Pepe et al. 2021). An early example of the application of planet-planet interactions as a probe of the planetary masses is GJ 876. It is among the first detected exoplanets (Marcy et al. 1998; Delfosse et al. 1998), and the dynamical interaction between the two Jupiter-mass planets in a 2:1 mean motion resonance was first modelled by Rivera et al. (2005) and later refined by Rivera et al. (2010) and Trifonov et al. (2018). Another interesting case is the GJ 581 system. The first three planets were detected by Bonfils et al. (2005) and Udry et al. (2007). Vogt et al. (2010) announced GJ 581 as a possible six-planet system, but stellar activity has been identified as the likely cause of three of the signals (Baluev 2013; Robertson et al. 2014; Hatzes 2016). The planet-planet interaction of the confirmed three planets in the system were recently used to infer the orbital inclinations and hence the masses from RV measurements spanning over about 7500 d (von Stauffenberg et al. 2024).

Additionally, from the point of view of planetary system formation and evolution, planets in multiple systems are very interesting, because comparative planetology allows constraining formation and evolution models. An example is the TRAPPIST-1 system, where precise mass and radius determinations of all planets indicate the formation out of a slightly metal-poor proto-planetary disk (Agol et al. 2021). Other examples include systems such as HD 15337, which contains two planets (Gandolfi et al. 2019) with a different composition on both sides of the radius gap (Dempsey & Nelson 2018), and both GJ 876 and TOI 216, which host planetary systems that exhibit mean motion resonances (Dempsey & Nelson 2018; Nesvorný et al. 2022). A total of 326 multi-planet systems with at least three planets are currently known (NASA Exoplanet Data from August 27, 2024). For those that do not transit, modelling of the planet-planet interaction using RV data can reveal their masses, given a sufficiently long observational time base.

Finally, the planets orbiting stars closest to our Solar System are the optimal targets for follow-up study of atmospheric chemistry. The mid M-dwarf GJ 1061 with a mass of 0.12 M_⊙ (Table 1), is of particular interest: (i) The system harbours three potentially rocky planets, two of them in the conservative habitable zone, detected using the RV technique: GJ 1061 b: P=3.204 ± 0.001 d, M sin i = 1.37 ± 0.16 M_⊙; GJ 1061 c: P=6.689 ± 0.005 d, M sin i = 1.74 ± 0.23 M_⊙; and GJ 1061 d: P=13.03 ± 0.03 d, M sin i = 1.64 ± 0.24 M_⊙ (Dreizler et al. 2020), (ii) it is located at a distance of only 3.6 pc making it the second closest compact multi-planetary system to our Solar System with at least three planets. The radii of the planets are unknown, as the planets do not transit (Lovos et al. 2022). The planetary system is dynamically packed, so that planet-planet interactions may be detectable using intensive RV monitoring, and we have sufficient data taken over a long baseline to provide a robust solution.

In this paper, we investigate the gravitational planet-planet interactions of the three planets orbiting the mid-M dwarf GJ 1061. In Section 2 we describe the observations and data processing, and in Section 3 we specify the methods and computer codes used in this work. In Section 4 we present our results and interpret them.

2 Observations and data

2.1 Spectroscopic data

For the measurement of the RV variation, we used data that are publicly available on the ESO Science Archive² from the HARPS³ instrument (Mayor et al. 2003) at the 3.6 m telescope of La Silla observatory, Chile. Table 2 lists the five campaigns that comprise the 184 observations of the full dataset. A total of seven spectra were recorded prior to the fibre change of the instrument in June 2015 and 177 post fibre change but before the 2020 warm-up. We treated these subsets of HARPS data as distinct datasets in this analysis.

The spectroscopic observations were reduced for both subsets separately using the serval algorithm, (Zechmeister et al. 2018) which uses an iterative template matching process to first obtain RVs relative to the template and then recalculate the template by co-adding the RV corrected spectra. It begins with the highest signal-to-noise (S/N) spectrum as the template, and is programmed to stop after the first co-added template is created and matched. Three spectra were discarded during this process by serval for their low S/N, leaving seven pre- and 174 post-fibre change spectra. By default, serval further discards the first ten HARPS echelle orders due to their lower S/N. Due to the late spectral type of GJ 1061 and correspondingly low flux in the blue-most orders, we decided to restrict this further to orders 20 to 72.

Outliers were removed by applying a conservative 5σ clipping, leaving 172 individual RV measurements. Two RV datasets were created: (1) averaged in the case of multiple observations per night, resulting in a final total of 142 data points, and (2) an unbinned dataset sensitive to any ultra-short period signals.

The HARPS data were complemented with 26 spectra from ESPRESSO⁴ (Pepe et al. 2021). We also used serval to retrieve the RV measurements, again skipping the bluest orders and using orders 16 up to 169. A 5σ clipping had no effect, the nightly averaging reduced the number of measurements to 25.

Along with the RV measurements, serval determines spectroscopic activity indicators including the Hα equivalent width, the chromatic index, and the differential line width (dLW). More details are described in Zechmeister et al. (2018).

2.2 Photometric data

GJ 1061 has been observed with MEarth⁵, a photometric monitoring project to search for transiting planets around ~3000 nearby mid-to-late M dwarfs, each using eight 40 cm-telescopes in the northern and southern hemisphere (Berta et al. 2012). GJ 1061 was observed by Telescopes 13 and 11 of MEarth-South at Cerro Tololo Inter-American Observatory (CTIO) in Chile from June 1, 2017 to February 27, 2022 and from October 20, 2016 to March 10, 2020, respectively. The 18 948 and 53 510 epochs were filtered with a 5σ clipping and then binned to 663 and 276 nightly averages.

There are also photometric observations from the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017; Christy et al. 2023). We retrieved the light curve from their online server⁶, which does not include proper-motion correction. We therefore corrected the coordinates for each year of observation (from 2014 to 2024), which were specified in each query, and then downloaded the data and combined them into a single (corrected) light curve. The observations were made with the V and g filters using different cameras.

As with MEarth, we also binned the ASAS-SN data nightly (1 d) and performed a 5σ clipping. After this, the cameras bE, bm, and bi (all g filters) had only one data point, so we discarded them. The bn camera (V filter) was also discarded as it had only 44 points with a time span of 50 d, covering the same time series as the bf camera (also V filter). Both the MEarth and ASAS-SN time series are shown in Fig. B.1. We note that the MEarth-11 data follows a long-term trend, and that the MEarth-13 dataset has a long observation gap (~700 d).

We also retrieved ‘cyan’ (c-band, 5200å) and ‘orange’ (o-band, 6600å) forced photometry from the Asteroid Terrestrialimpact Last Alert System (ATLAS, Tonry et al. 2018; Heinze et al. 2018). Due to the high proper motion (0.83″ yr⁻¹) of GJ 1061 compared to the 1.86″ pixel scale of ATLAS, we performed the retrieval in 15 six-month intervals. GJ 1061 is saturated (median magnitude of 11.2) in o-band but marginally not in c-band (13.1) and we performed a Lomb-Scargle periodogram analysis (Scargle 1982) on the 444 points in the c-band light curve with a median precision of 0.002 mags (Fig. 1). The periodogram (Fig. 2) contains a marginally significant peak at a ~125 d signal, which might be the rotational signal tentatively observed in the RV analysis. An overview of the photometric data is listed in Table B.2.

Fig. 1 ATLAS c-band light curve of GJ 1061. Some error bars are smaller than the points.

Fig. 2 Lomb-Scargle periodogram of the ATLAS c-band photometry of GJ 1061. The green-line is the p = 0.01 false-alarm threshold calculated using, and the blue curve is the maximum of the ensemble of 100 periograms of the randomly scrambled photometry. The red dotted lines mark the lunar synoptic period and its 1:2 and 3:2 harmonics.

3 Methods

The measured RV values were analysed using a multi-planetary model, employing either Keplerian or Newtonian orbits. In the Keplerian case, the radial velocity is derived analytically, while in the Newtonian case, it is obtained through direct N-body integration using rebound with the IAS15 integrator (Rein & Liu 2012; Rein & Spiegel 2015). The main difference is that in the Keplerian models the orbit of each planet is calculated as if it were the only planet in the system, whereas in the Newtonian models the mutual gravitational interactions between all bodies are explicitly taken into account. Small perturbations due to these interactions manifest as apsidal precession. Over time, this causes a Keplerian orbit to exhibit a phase shift relative to the observed RV orbit. The strength of the perturbation depends on the planetary masses, while the RV amplitude is determined by the planetary mass and orbital inclination. Due to the perturbations, the orbital elements are not constant but evolve with time. Energy and angular momentum exchange lead to variations in the orbital period and eccentricity. Therefore, RV measurements over a sufficiently long time baseline allow for the determination of planetary masses and inclinations, constrained by both the RV amplitudes and the phase drift of the orbital motions.

The Keplerian model is parameterised with the semiamplitude K, the orbital period P, the eccentricity e, the longitude of periastron ω, and the mean longitude λ. While the inclination i and the longitude of the ascending node Ω remain undetermined in the Keplerian model, both can in principle be determined in the case of dynamical orbits from the strength of the planet-planet gravitational interaction. For Newtonian orbits, we replace K by the planetary mass m and add i and Ω as additional free parameters. It should be noted that the dynamical model can only provide Ω relative to for example that of the innermost planet, where Ω = 0 is assumed. In the scenario where planet-planet interactions are weak, a simplification assuming co-planar orbits might be appropriate. In that case, only a common inclination is determined and Ω = 0 is set for all planets.

As discussed in Liebing et al. (2024), the parametrisation of e and ω has an impact on the posteriors of other planetary parameters. We therefore used either e and ω directly or alternatively $h^{\mathrm{\prime }}=\sqrt{e}\sin \omega$ and $k^{\mathrm{\prime }}=\sqrt{e}\cos \omega$ to check the impact of the parametrisation on the results. Since the planetary system of GJ 1061 is compact, we expect small eccentricities in order to keep it dynamically stable.

For the Newtonian orbits, we checked the critical two-planet separation in units of the mutual Hill radius (Equation (5), Obertas et al. 2017). We rejected parameter combinations violating this critical Hill-radius separation criterion to avoid computational time being spent on physically unlikely parameters.

The planetary model of varying complexity is complemented with Gaussian Process (GP) regression (Sect. 3.1) to model possible rotational modulated contributions from stellar activity. The posterior distributions of the parameters and hyperparameters are determined with a Nested Sampling algorithm (Sect. 3.2) which also provides the Bayesian evidence which is used to identify the best model (Sect. 4.2). Finally, the parameters of the best model are checked for orbital stability (Sect. 3.3).

3.1 Gaussian process regression

Gaussian process regression is a powerful non-parametric Bayesian method for modelling complex functions. GPs provide a flexible way to model data by assuming that the underlying function is a sample from a Gaussian distribution. The python package (celerite, Foreman-Mackey et al. 2017) is an efficient implementation of GP regression that is specifically designed for large datasets, since its computational demands scale linearly with the number of data points. This requires that the covariance matrix (the kernel) is modelled as a sum of complex exponential functions. Foreman-Mackey et al. (2017) suggest using a simple damped harmonic oscillator (SHO) as kernel in case of modelling the covariance caused by stellar activity. Mathematical details are given in Equations (9), (20), and (21) in Foreman-Mackey et al. (2017). An improved kernel for rotationally modulated stellar activity is presented by Foreman-Mackey (2018) and it consists of two SHO kernels with the period of the second being the second harmonic of the first. This dSHO kernel is parameterised by the variance σ, the oscillator (rotation) period P, the quality factor Q₀ and the difference of the quality factor of the first and second period dQ, which both represent the degree of damping, and the fractional amplitude between the two oscillators f. The resulting damping timescale is enforced to be the same for both oscillators. We used wide priors (see Table C.1) for these parameters and required that Q > 1/2 to keep the oscillators in the under-damped regime. We also allowed for additional white noise using a jitter term, that is a covariance matrix with diagonal elements only.

3.2 Nested sampling

Nested sampling (Skilling 2004) is a Bayesian inference method that can be used to estimate the logarithm of the evidence (ln 𝒵) and the posterior distributions of a model. The python package dynesty (Higson et al. 2019; Speagle 2020) is an implementation of nested sampling that can be used to perform parameter estimation and model selection. It can be used to efficiently explore complex posterior distributions, and it has a user-friendly interface that makes it easy to use in a variety of applications. One of the key features of dynesty is its dynamic nested sampling algorithm, which allows it to adapt to the structure of the posterior distribution in real-time, making it more efficient than other nested sampling methods. Additionally, the package provides a wide range of diagnostic tools and options for controlling the sampling process, enabling users to fine-tune the method to their specific needs.

We used 5000 live points, bound=multi, sample=rwalk, and stop when δ ln 𝒵 < 0.01. The model selection was made using the static sampler, while the final model was run with the dynamic sampler, optimising the posterior distributions. Multiple runs for model D2, the Keplerian three-planet model with a uniform prior for the period of the SHO kernel (see Table 3), resulted in very similar posteriors, indicating that the number of live points were sufficient.

3.3 Post-processing using SPOCK

The Stability of Planetary Orbital Configurations Klassifier (SPOCK, Tamayo et al. 2020) is a machine learning algorithm that is designed to predict the long-term stability of planetary systems over 10⁹ orbital periods of the innermost planet. The algorithm is trained on a set of numerical simulations of about 100 000 nearly co-planar planetary systems with three or more planets with masses ranging from Mars- to twice Neptune-mass. It uses a combination of physical and dynamical features to classify the stability of a given system. It has been used to predict the stability of observed exoplanetary systems (e.g. Tamayo et al. 2020; Kaye et al. 2022), and it has been shown to be accurate in identifying stable systems. Additionally, it has been used to identify new configurations of exoplanetary systems that are likely to be stable, which can help guide the search for new exoplanets.

We used SPOCK to check the resulting planetary system configurations for dynamical stability. We emphasise that we used the posterior distributions and not draws from the normal distributions given by the mean and standard deviation tabulated in Table 4.

3.4 Light curve modelling

To obtain the stellar rotation period we performed a GP fit on all the data with juliet (Espinoza et al. 2019), which provides Keplerian components by fitting RVs (Radvel, Fulton et al. 2018), transits for photometric models (batman, Kreidberg 2015), and a red noise component (GP) using the publicly available packages celerite and george (Ambikasaran et al. 2015). We used the dSHO kernel for the modelling of the stellar rotation modulation, fitting the GP hyperparameters together with a linear trend for each dataset.

4 Analysis results

4.1 Stellar rotation period

We computed the generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009) of the activity indicators and the photometric data (Fig. B.1). In Fig. B.2 we display the periodograms which have signals reaching at least a FAP level of 10% (excluding yearly and half-yearly periodicities). In addition to the signal at 1 d caused by nightly observations, the most dominant signal (FAP < 0.1%) corresponds to 130 d in one of the ASAS-SN datasets (bf cam, V filter), which is also present in the MEarth-13 data, but the latter has several significant periodicities between 100 d and 200 d. Other significant signals (FAP < 0.1%) are found in the MEarth-11 data, at 160 d, which become stronger after fitting a trend (see Section 2.2). At 200 d there are also peaks in the ASAS-SN bn and bF cam (g filter), as well as in MEarth-11, the latter also showing a dominant peak around 300 d. These long periodicities could be related to a long magnetic cycle, since stellar rotation periods have shorter timescales (Shan et al. 2024).

We used the dynamic nested sampling algorithm from dynesty (see Section 3.2 and Fig. B.2) and chose a dSHO kernel with uniform priors on the period from 1 to 200 d (the rest of the priors are listed in Table B.3). We first fitted a GP to each photometric dataset separately, obtaining only a clear detection of the GP period on MEarth-11 (~160 d, trend fitted), MEarth-13 (~46 d) and ASAS-SN bf cam (V filter, ~137 d). Multiple star spots at different longitudes can result in sub-rotational period signals, which could be the reason for the detection of shorter periods in some datsets. A possible activity cycle similar to that in our Sun could lead to non-detections of photometric variability in some datasets. In this sense, we cannot retrieve a reliable stellar rotation period, but we can make an estimate. Firstly, we do not rely entirely on the MEarth-13 data alone because of the long gap. Since MEarth-11 and ASAS-SN bf cam (V filter), which are the only remaining datasets where a periodicity has been found, give long periods (P > 120 d), we expect a rotation period of this order. In this sense, we can rely on the combination of the MEarth-11 (trend fitted) with the ASAS-SN bf cam V filter as an estimate of the stellar the rotation period. We then combined MEarth-11 (trend fitted) with ASAS-SN bf cam (V filter) and obtained a stellar rotation period of ${130}_{-14}^{+5}$ d.

The posteriors are listed in Table B.1, and we also depict the GP hyperparameters (dQ and Q) of the dSHO kernel against the periodic component in Fig. 3. The likelihood and number of posterior samples favour a rotation period of ~125 d. A better sampled photometric dataset (i.e. with no block separation, sampled for at least three probable rotation periods) is required in order to have an unambiguous detection.

The spectroscopic activity indicators, namely the Hα equivalent width probing the chromospheric activity, the chromatic index indicating RV variations do to changing spot contrast, and the dLW characterising line shape variations (2.1), do not show a significant signal near the potential rotation period or its harmonics. Only the latter shows a weak signal at about 125 d (Fig. B.2).

The inferred rotation period of ~125 d is compatible with the low rotational velocity v sin i < 2 km s⁻¹ (Dreizler et al. 2020) and the estimate age of older than 7.0 Gyr (West et al. 2008) and in agreement to the estimated rotation period from the previous analysis (Dreizler et al. 2020). As discussed below, the RV data also show a signal at 125 d, which we interpret as the rotational induced RV red noise. The corresponding modulation is shown in Fig. 7 as a thin grey line. Its amplitude is below 2 m s⁻¹.

4.2 Planetary model selection

In a first analysis step, we selected the model with the highest log 𝒵. As a selection criterion, we used the difference in the natural logarithm of the Bayesian evidence Δ ln 𝒵. Following Trotta (2008), we used the threshold of Δ ln 𝒵 > 5 as strong evidence and Δ ln 𝒵 > 2.5 as moderate evidence for a more complex model being better than a simpler model. For the model selection we used Keplerian models since the deviation between Keplerian and Newtonian models are rather small so that we do not expect significant differences in the model selection caused by this simplification.

Starting with our base model which allows for offsets and jitter terms for the two HARPS and the ESPRESSO datasets separately (Model A0, Table 3), we recursively added a Keplerian orbit at the frequency with the highest power in the GLS periodogram (Zechmeister & Kürster 2009) of RV residuals until 𝒵 was no longer increasing (Models B0 to F0, Table 3) and no peak in the periodogram was above the 1% false alarm probability. We used narrow priors for the three known planets with a width of the frequency resolution given the length of the dataset. After adding five Keplerian signals to the model, no significant power was left in the GLS periodogram of the residuals. The first three of these Keplerian models are the already known planets at 6.68, 3.20 and 13.1 d (Models B0 to D0) from Dreizler et al. (2020). The signal at 55 d has been identified as a planet candidate by Dreizler et al. (2020). In Fig. 4 we show periodograms where we have successively subtracted the planet signature and the red-noise modulation.

Instead of interpreting the 55 d signal as a planet, the long period signals were modelled with a SHO and dSHO kernel. This was motivated by the varying period and power of the Stacked Bayesian GLS (s-BGLS) periodogram (Mortier & Collier Cameron 2017) of signals in the long-period range above 50 d (Fig. 5), indicating activity induced radial velocity variations. The 55 and the 125 d periods are also close to the potential rotation period or its second harmonic (Sect. 4.1). In particular, the strength of the 55 d signal declines after N_obs ~100 as more observations were added. The comparison shown in Table 3 reveals that the more complex dSHO kernel is not better than the SHO kernel when using wide priors for the rotation period. The resulting periodicity of about 125 d is close to the rotation period (Table C.1, Fig. C.2) estimated from the analysis of the photometry (Sect. 4.1). A normal prior around the potential stellar rotation period resulted in a slightly lower 𝒵 compared to a wide uniform prior. Models with Gaussian process regression are significantly better (Δ ln 𝒵 > 5) compared to those without for three and four planets (Models D0 versus D2 and E0 versus E2). In the case of the five planet models (F0 versus F2), all have nearly identical evidence. Together, this indicates that we can confirm the three known planets of GJ 1061 (Dreizler et al. 2020). The additional signals at 55 d and 125 d are more likely due to stellar rotation, although the five-planet model (F0) is statistically moderately better than the three-planet model including the SHO kernel (D2) to model the rotational RV variability. No additional planets can be detected in the current data.

Having established the best Keplerian model, we then modeled the planet-planet interaction. We did this by using Newtonian models (see Section 3) for models with three to five planets using a stellar mass of 0.12 M_⊙, shown in the lower block in Table 3. We first checked that the parametrisation of the eccentricity and the longitude of periastron has only a minor impact on the derived model parameters. The final model wss therefore run with $h^{\mathrm{\prime }}=\sqrt{(}e)\sin \omega$ and $k^{\mathrm{\prime }}=\sqrt{(}e)\cos \omega$ using uniform priors within [−0.4, 0.4], the mean longitude λ was sampled uniformly within [0, 2 π] employing periodic boundaries (see also Table 4). We then compared the results from a co-planar (Model D4, D5, D6) to a full three-dimensional (3D) model, where a common inclination and a fixed longitude of the ascending node of Ω = 0 for all planets were fitted in the former. We also used three different GP kernels, the SHO kernel with uniform and normal priors, as well as the dSHO kernel with uniform priors. Like in the case of the Keplerian models, the SHO kernel results in better 𝒵 values. The Newtonian models are always significantly better than the corresponding Keplerian models. The GP kernels do have a slight impact on the posteriors of the planet masses (see Fig. C.1). We therefore used a weighted mean for deriving the mass uncertainties.

The results presented in Table 4 show that the mutual inclinations derived from the 3D-model are identical within the uncertainties. In combination with the fact that the co-planar model has a significantly higher 𝒵, identifying this as the statistically preferred model, we concluded that deviation from a co-planar configuration are too small to be detected, and we can thus only constrain the inclination of a co-planar configuration. Nevertheless, the planet masses are very similar between the two solutions, indicating that a co-planar model does not introduce major systematic uncertainties.

To evaluate the impact of the ESPRESSO data, we also used the co-planar model with an identical set-up of model D5 but restricted it to the HARPS data. The results are listed as a third block in Table 4. For the two inner planets in particular, the inclusion of ESPRESSO data led to a reduction in the uncertainty of the mass determination from about 20% down to 10% as well as a systematic decrease due to a shift in the inclination. When the uncertainty of the modelling of the stellar noise was accounted for, the error bars for the planetary masses were slightly larger, especially the upper limits for planets c and d are higher. Nevertheless, their masses are clearly below 2 M_⊕ and hence likely rocky.

Our final best model is the co-planar three-planet model, including an SHO kernel to model the rotation induced RV variability (D4 and D5 in Table 3). The two models differ in the priors for the rotation period (normal versus uniform, see Table 3) and show nearly equal values for ln 𝒵 and we selected model D5 with the wider uniform prior as the final best model. The posteriors distribution and the best parameters are those in the third and fourth column of Table 4. They are also displayed as corner-plots in Fig. C.2. The main result is the determination of the masses of the planets due to the planet-planet gravitational interaction. This is clearly visible in the phase folded RV data overlaid with the model shown in Fig. 6, where we can see the orbital precession and the resulting phase shift, especially in the two inner planets. The comparison of the model including the 1σ uncertainties with the RV data as function of time is shown in Fig. 7. The very good match visible between the data and model provides confidence in the mass determination from the planet-planet interaction. This is possible due to the ~20-year baseline of HARPS data and high-precision RV measurements of the ESPRESSO data.

Using the posterior parameter distribution, we used SPOCK to derive a probability for each configuration to be stable. The distribution of the stability probability of Model D5 has a median of 0.74 ± 0.1, 99.5% are above the stability threshold used by Tamayo et al. (2020). The probabilities for all parameter configurations were then taken as weights and stability corrected posteriors were then obtained. The parameter distributions were unaffected by the stability correction. We therefore conclude that the derived planetary architecture of the GJ 1061 system is dynamically stable.

Recently, Luque & Pallé (2022) reported that planetary density and not radius is the optimal parameter to characterise the internal composition of sub-Neptunian planets. Using their results, showing all planets of M stars below 2 M_⊕ having Earth density, we can infer that the three planets orbiting GJ 1061 are very likely rocky planets as opposed to water-rich or gas-rich. Assuming an Earth bulk density, the planetary radii can be inferred and are listed in Table 4. GJ 1061 c and d are therefore likely rocky planets in the optimistic and conservative habitable zone as derived by Kopparapu et al. (2014). It makes GJ 1061 c and d the third and fourth closest rocky habitable zone planets, after Proxima b and Ross 128 b.

These results make the GJ 1061 system an optimal target for further characterisation instruments combining high-resolution spectroscopy with Adaptive Optics assisted high spatial resolution like the ArmazoNes high Dispersion Echelle Spectrograph (Palle et al. 2023) or the Planetary Camera Spectrograph (Kasper et al. 2021) and with space missions such as the Habitable Worlds Observatory (HWO) or the complementary Large Interferometer For Exoplanets (LIFE) mission (Quanz et al. 2022a,b) in thermal emission. For the latter mission, the planets of GJ 1061 are three out of 212 planets detectable (signal-to-noise > 7 in less than 100 h) with the reference configuration of LIFE for host stars within 20 pc (Carrión-González et al. 2023).

Fig. 3 Hyperparameters of the GP model for the photometric data: The quality factor Q versus the rotation period (top) and the difference in quality factor between primary and secondary oscillator (bottom). The colour indicates the likelihood.

Fig. 4 Periodograms of models with an increasing complexity: from top to bottom: zero, one, two, and three planets as well as three planets + SHO kernel. Red arrows mark the subtracted signal for the following periodogram, the blue arrow indicates the potential stellar rotation period at ~125 d, the light blue arrow points at the 55 d signal. The labels indicate the model which is subtracted from the observation ‘O’ (see Table 3).

Fig. 5 Stacked Bayesian periodogram (s-BGLS) of the HARPS post-fibre upgrade data after subtraction of three planets (model D0).

5 Summary

We have validated the rocky nature of the three planets orbiting the nearby low-mass star GJ,1061 using 198 RV measurements collected with HARPS and ESPRESSO. The dense coverage, in part due to the RedDots campaigns (PI: Jeffers, ESO ID 0101.C-0516(A) and 0112.C-2112(A)), combined with the long observational baseline and the exceptional precision of the ESPRESSO data, allowed us to detect planet-planet gravitational interactions. By modelling these interactions with N-body integrated orbits, we achieved a more comprehensive and accurate characterisation of the planetary system.

We determined the masses of all three planets orbiting GJ, 1061 with a precision of about 10%. Each of the planets has a mass below 2 M_⊕, indicating that they are likely rocky planets with Earth-like bulk densities. The mutual inclinations suggest a co-planar configuration. Such well-characterised, nearby planetary systems are excellent candidates for future direct atmospheric investigations.

Additionally, using photometric time series and spectroscopic activity indicators, we determined the stellar rotation period to be approximately 125 days. Notably, two of the planets, planets c and d, are the third and fourth closest rocky planets to our Solar System that orbit within the optimistic habitable zone of their host star. Continued long-term monitoring, supported by high-precision RV measurements from instruments like ESPRESSO, would enable the determination of dynamical masses in more compact planetary systems orbiting M-type stars.

Fig. 6 Data versus model D5 (Table 3) for planets b, c, and d from top to bottom: to show the temporal evolution of the planet’s orbits, the colour of the data and the model changes from blue to red over the observational baseline to visualise the non-Keplerian orbit due to the planet-planet interaction noticeable as a phase shift. A similar figure, but separated into observing seasons is shown in Fig. D.2.

Fig. 7 Data versus Model D5 (Table 3) as a function of time. The grey shaded region comprise the 1σ uncertainty of the model, the slowly varying line is the contribution from the rotational modulation.

Acknowledgements

S.D. acknowledges support from the Deutsche Forschungsgemeinschaft under Research Unit FOR2544 “Blue Planets around Red Stars”, project no. DR 281/39-1. S.D., S.V.J., F.L. and Y.T. acknowledge the support of the German Science Foundation (DFG) priority programme SPP 1992 ‘Exploring the Diversity of Extrasolar Planets’ (DR 281/37.1, JE 701/5-1, TS 356/3-1). C.A.H. and J.R.B. acknowledge support from STFC under consolidated grant ST/X001164/1; H.R.A.J. acknowledges support from STFC consolidated grant ST/R000905/1. E.R. acknowledges financial support from the Spanish Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación through projects PID2019-107061GB-C64, PID2022-137241NB-C43, and the Centre of Excellence ‘Severo Ochoa’ Instituto de Astrofísica de Andalucía (grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033). This paper makes use of data from the MEarth Project, which is a collaboration between Harvard University and the Smithsonian Astrophysical Observatory. The MEarth Project acknowledges funding from the David and Lucile Packard Fellowship for Science and Engineering, the National Science Foundation under grants AST-0807690, AST-1109468, AST-1616624 and AST-1004488 (Alan T. Waterman Award), the National Aeronautics and Space Administration under Grant No. 80NSSC18K0476 issued through the XRP Program, and the John Templeton Foundation. ASAS-SN is funded in part by the Alfred P. Sloan Foundation under grant G202114192. Making ASAS-SN light curves public is primarily funded by grants GBMF5490 and GBMF10501 This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. The Asteroid Terrestrial-impact Last Alert System (ATLAS) project is primarily funded to search for near earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogs from the survey area. This work was partially funded by Kepler/K2 grant J1944/80NSSC19K0112 and HST GO-15889, and STFC grants ST/T000198/1 and ST/S006109/1. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen’s University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory, and The Millennium Institute of Astrophysics (MAS), Chile.

Appendix A ESPRESSO RV data

In Table A.1 we list the radial velocity measurements from ESPRESSO using serval, described in Section 2.1.

Appendix B Photometric monitoring and analysis of rotation period

The information about the instruments as well as the observing campaigns of the photometric monitoring of GJ 1061 described in Section 2.2 is listed in Table B.2, and the time series and corresponding periodograms are displayed in Fig. B.1 B.2 respectively.

Fig. B.1 Time series of nightly binned photometric monitoring of GJ 1061.

Fig. B.2 Left: GLS periodograms of photometric time series and the spectroscopic activity indicator differential line width (dLW). Other activity indicators do not show significant signals. The blue and yellow line mark the year and half year period. The red line indicates the potential rotation period at 125 d. Right: the window function of the time series.

Appendix C Priors and posterior distributions

Fig. C.1 Posteriors for the planetary masses from models D4 to D6 using three different kernels (Table 3) and the weighted combination of the three in black.

Fig. C.2 Corner plot for the best-fit model D5, see Table 3.

Appendix D Osculating orbital elements

In addition to Fig. 6, we show here the RV data separated in seven chunks for the three planet. The non-Keplerian orbit result in a variation of the order of a minute of the orbital period, as well as in an apsidal precession. The latter we also display as a function of time for planetary system parameters randomly drawn from the posteriors.

Fig. D.1 The change of the longitude of the pericentre (equal to the argument of pericentre in a co-planar planetary system) for planets b, c, and d from left to right.

Fig. D.2 Alternative visualisation of the non-Keplerian orbit (see Fig. 6) for planets b, c, and d from left to right. The upper four panels show the spares HARPS observations overt the first ~ 3000 days. The next two show the two densely sampled observing seasons with HARPS between 5000 and 5600 days, the lower panel shows the recent ESPRESSO data.

References

All Tables

All Figures

	Fig. 1 ATLAS c-band light curve of GJ 1061. Some error bars are smaller than the points.
In the text

	Fig. 2 Lomb-Scargle periodogram of the ATLAS c-band photometry of GJ 1061. The green-line is the p = 0.01 false-alarm threshold calculated using, and the blue curve is the maximum of the ensemble of 100 periograms of the randomly scrambled photometry. The red dotted lines mark the lunar synoptic period and its 1:2 and 3:2 harmonics.
In the text

	Fig. 3 Hyperparameters of the GP model for the photometric data: The quality factor Q versus the rotation period (top) and the difference in quality factor between primary and secondary oscillator (bottom). The colour indicates the likelihood.
In the text

Fig. 4 Periodograms of models with an increasing complexity: from top to bottom: zero, one, two, and three planets as well as three planets + SHO kernel. Red arrows mark the subtracted signal for the following periodogram, the blue arrow indicates the potential stellar rotation period at ~125 d, the light blue arrow points at the 55 d signal. The labels indicate the model which is subtracted from the observation ‘O’ (see Table 3).

In the text

	Fig. 5 Stacked Bayesian periodogram (s-BGLS) of the HARPS post-fibre upgrade data after subtraction of three planets (model D0).
In the text

Fig. 6 Data versus model D5 (Table 3) for planets b, c, and d from top to bottom: to show the temporal evolution of the planet’s orbits, the colour of the data and the model changes from blue to red over the observational baseline to visualise the non-Keplerian orbit due to the planet-planet interaction noticeable as a phase shift. A similar figure, but separated into observing seasons is shown in Fig. D.2.

In the text

	Fig. 7 Data versus Model D5 (Table 3) as a function of time. The grey shaded region comprise the 1σ uncertainty of the model, the slowly varying line is the contribution from the rotational modulation.
In the text

	Fig. B.1 Time series of nightly binned photometric monitoring of GJ 1061.
In the text

	Fig. B.2 Left: GLS periodograms of photometric time series and the spectroscopic activity indicator differential line width (dLW). Other activity indicators do not show significant signals. The blue and yellow line mark the year and half year period. The red line indicates the potential rotation period at 125 d. Right: the window function of the time series.
In the text

	Fig. C.1 Posteriors for the planetary masses from models D4 to D6 using three different kernels (Table 3) and the weighted combination of the three in black.
In the text

	Fig. C.2 Corner plot for the best-fit model D5, see Table 3.
In the text

	Fig. D.1 The change of the longitude of the pericentre (equal to the argument of pericentre in a co-planar planetary system) for planets b, c, and d from left to right.
In the text

	Fig. D.2 Alternative visualisation of the non-Keplerian orbit (see Fig. 6) for planets b, c, and d from left to right. The upper four panels show the spares HARPS observations overt the first ~ 3000 days. The next two show the two densely sampled observing seasons with HARPS between 5000 and 5600 days, the lower panel shows the recent ESPRESSO data.
In the text