© The Authors 2024

1 Introduction

The successful search for exoplanets has been a remarkable achievement in modern astronomy, resulting in the discovery of over 5000 planets outside the solar system. Furthermore, advancements in radial velocity (RV) surveys have led to an increasing number of low-mass planet detections. A small subset of these known planets holds special interest: rocky planets located within the habitable zone of their host star, where conditions could potentially support liquid water on their surfaces (Kasting et al. 1993; Kopparapu et al. 2013). The Habitable Exoplanets Catalog currently lists only 24 Earth-sized planets in the conservative sample of potentially habitable exoplanets¹.

It is noteworthy that the majority of these planets have been found orbiting M dwarfs. Among them, the TRAPPIST-1² system (Gillon et al. 2017), with four planets in the list, is the only system with precisely determined planetary masses and radii. Without the requirement of being within the habitable zone, the catalog of Transiting M dwarf Planets (TMP)³ lists 21 potentially Earth-like planets (M_planet < 2 M_⊕) including the seven rocky planets of the TRAPPIST-1 system with mass determinations from RVs or transit timing variations.

Multi-planet systems are also of great interest as they provide valuable insights into planet formation and evolution. As of April 2023, a total of 850 multiple planet systems have been discovered⁴. Lin et al. (2021) found that the number of planets in these systems is related to the size of the protoplanetary disk. Additionally, these authors show that the timescale for protoplanet appearance plays a crucial role in determining the planet configuration arising from resonant trapping. A shorter timescale results in a larger number of formed planets, which become trapped in more closely spaced resonances. Their simulations also indicate that resonant planets are generally not formed around stars with masses larger than about 0.4 M_⊙. For a comparison of these predictions with observations, we should therefore aim for the most complete knowledge of the number of planets in multi-planet systems. Another important aspect seems to be the description of the gas disk interaction with the planetary embryos, as discussed by Sánchez et al. (2022). The similarity in planet properties and the arrangement of their orbits, forming a chain of mutual resonant configurations, led Ormel et al. (2017) to conclude that the TRAPPIST-1 system and its planetary architecture can be well explained by the growth of planets driven by pebbles at the water ice line, followed by inward migration. As discussed by Drążkowska et al. (2023), such planetary systems offer valuable constraints on planet formation scenarios, thereby motivating further characterization efforts.

Earlier studies have indicated that low-mass stars often host Earth mass planets at relatively short orbital periods (Dressing & Charbonneau 2013, 2015). In such systems, planets within the habitable zone are located closer to the star, making them more accessible for investigations. Moreover, M dwarfs, which represent the most common class of stars in the solar neighborhood (e.g., Reylé et al. 2021; Golovin et al. 2023), have extremely long-lived main sequence phases, making them excellent candidates for studying potentially habitable systems. Planets around nearby M dwarfs are also the most suitable targets for future very high contrast and spacial resolution imaging, allowing to probe exoplanet atmospheres. The planetary atmosphere investigations will be possible with ground-based instruments in the era of extremely large telescopes in reflected light with instruments 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) which has been recommended to NASA as the next flagship mission by the US Astro 2020 Decadal Survey report (National Academies of Sciences, Engineering and Medicine 2023) or the complementary Large Interferometer For Exoplanets (LIFE) mission (Quanz et al. 2022a,b) in thermal emission. Carrión-González et al. (2023) investigated the detectability of exoplanet atmospheres around stars within 20 pc. From the 212 planets detectable (signal-to-noise > 7 in less than 100 h) with the reference configuration of LIFE, 49 can also be detected with the notional HWO, 163 with LIFE only. From the 38 LIFE-detectable planets in the habitable zone, 13 are below five Earth masses.

One particularly remarkable system is Teegarden’s Star, an M7.0 dwarf (Alonso-Floriano et al. 2015) discovered by Teegarden et al. (2003). This is the 25th nearest star to the Sun⁵ (Reylé et al. 2021), at a distance of only 3.831 pc. Another significant attribute of Teegarden’s Star is its relatively low magnetic activity compared to most late-M dwarfs, which enhances its potential as a target for the search for extraterrestrial life. In 2019, Zechmeister et al. (2019) presented evidence of two Earth mass planets within the potentially habitable zone of Teegarden’s Star, with orbital periods of 4.91 and 11.4 days, respectively. Both planets are among the 13 low-mass habitable-zone planets detectable with LIFE. These were the first, and at present still are, the only planets detected around an ultra-cool dwarf (spectral type later than M 7.0 V) using radial velocities. Subsequent data collection allows us now to search for weaker signals in the system so that a more complete inventory of low-mass planets can be obtained, providing more reliable constraints on planetary architectures and properties around very low-mass stars. Here we report the discovery of a third planet orbiting Teegarden’s Star and find evidence for further suggestive signals that could point to an even larger number of planets in the system, thereby resembling TRAPPIST-1.

First, we introduce the RV measurements, the spectroscopic activity indices, and transit data in Sect. 2. In Sect. 3, we describe the methods used for analysis. The stellar parameters and results are presented and discussed in Sects. 4 and 5, respectively.

2 Observations and data products

2.1 Radial velocity data

RV measurements were obtained with four different instruments described below. Based on the instrument-specific data reduction we computed the RVs from all instruments with serval⁶ (Zechmeister et al. 2018).

The CARMENES instrument (Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical Echelle Spectrographs) at the 3.5 m telescope at the Calar Alto Observatory in Almería, Spain, is a dual-channel spectrograph that operates at both optical (VIS, 0.52–0.96 μm) and near-infrared (NIR, 0.96–1.71 μm) wavelengths. The average resolving power for the two wavelength regions is ℜ = 94 600 and ℜ = 80 400, respectively. We used the available 262 measurements (July 2023) of the RV of Teegarden’s Star from the guaranteed time (GTO) and legacy project observations of the CARMENES project (Quirrenbach et al. 2014; Ribas et al. 2023). Compared to the 239 measurements with nightly zero point corrections (Trifonov et al. 2018) used by Zechmeister et al. ((2019), we now added 14 more published by Ribas et al. (2023). Additionally to the GTO data, nine measurements were taken in the last observing season (July 2022 to March 2023) as part of the CARMENES Legacy-Plus project. Due to the larger scatter of the NIR data and the small signals expected, we restricted the analysis to the VIS data of CARMENES. Bauer et al. (2020) compared the scatters of the VIS and NIR RV curves for an M dwarf of roughly the same J magnitude as Teegarden’s Star. During part of the time during which these data where acquired, the NIR channel had a peak-to-peak instrumental drift substantially larger than that of the VIS channel, and suffered from instrumental changes during maintenance operations of the active cryogenic cooling system. We like to note that recent instrumental updates have significantly reduced this effect. The impact of the telluric lines on the RV measurements has been tested using CARMENES spectra, which were cleaned from telluric contamination (Nagel et al. 2023), but we found very small differences compared to the standard results from serval. Nevertheless, we used the telluric-corrected RV data.

ESPRESSO (Pepe et al. 2021) is the Echelle Spectrograph for Rocky Exoplanets and Stable Spectroscopic Observations installed at the incoherent combined coudé facility of the Very Large Telescope (VLT). It is an ultra-stable fiber-fed échelle high-resolution spectrograph (ℜ = 140 000, 0.38–0.79 μm in high resolution single unit telescope mode). We obtained 11 ESPRESSO measurements in September 2019.

MAROON-X (Seifahrt et al. 2016) is a stabilized, fiber-fed high-resolution (ℜ = 85 000) spectrograph mounted at the 8.1 m Gemini North telescope on Mauna Kea, Hawai’i, USA. It has blue and red arms, which encompass 0.50–0.68 μm and 0.65–0.92 μm, respectively. During an observation, both arms are operated simultaneously. MAROON-X obtained nine, seven, and eight measurements with the blue and red arms in three observing runs in August and October 2021, and in August 2022, respectively.

The Habitable Zone Planet Finder (HPF) is a stabilized fiber-fed near-infrared spectrograph (0.9–1.8 μm) for the 10 m Hobby-Eberly Telescope (HET, Ramsey et al. 1998) with a resolution of ℜ = 53 000 (Mahadevan et al. 2012). Between September 2019 and October 2021, a total of143 measurements were collected, with each night’s observations organized in pairs of adjacent measurements. The spectra were extracted with the HPF pipeline HxRGproc (Ninan et al. 2018; Kaplan et al. 2019; Metcalf et al. 2019).

Overall, we have 228 new measurements extending the time baseline by about 3 yr compared to Zechmeister et al. (2019). In case of multiple observations per night, we used nightly averages with the nightly standard deviation as uncertainty for the night. We also excluded outliers using a 5σ-clipping procedure. This resulted in 355 measurements after σ-clipping and nightly binning for the full dataset (Table 1). The original RV data, that were un-binned and un-clipped, are available at CDS. A short section is shown in Table A.1.

2.2 Spectroscopic indices

From the CARMENES and MAROON-X data, we extracted several spectroscopic indices using serval. Here we make use of the chromatic index (CRX), the differential line width (dLW), and the emission strength of the hydrogen Hα line. In short, the CRX index measures the wavelength dependence of the RVs determined in each spectral order. While planets do not have a color-dependent RV, stellar activity signals from spots and plagues are likely to show them due to the surface temperature modulations. The dLW is similar to the full width half maximum (FWHM) index from cross-correlation techniques. It measures line profile variations, which can either originate from instrument changes or from stellar activity. The equivalent width of the hydrogen Hα line is notoriously difficult to determine in M dwarfs due to the absence of a continuum free of spectral lines. A slowly rotating template stars for each M spectral subclass is therefore use as reference. For details about the indices we refer to Zechmeister et al. (2018) and Schöfer et al. (2019).

Although the spectra were cleaned from telluric lines (Nagel et al. 2023) and wavelength regions with strong telluric contamination were masked for RV measurements, residual telluric lines may still have an impact. We simulated the impact of telluric lines on the RV determination using an appropriate synthetic spectrum and adding a telluric spectrum shifted by the barycentric velocity of each observing date. No planetary signal was added. The result was a time series of simulated spectra, which were then analyzed using serval to determine the RVs. Without a planetary RV signal in the simulated data a measured RV is then due to telluric lines, resulting in a telluric contamination index.

2.3 Photometric data

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014, 2015) observed Teegarden’s Star in sectors 42 and 44 in full frame image mode and in sectors 43, 70, and 71 with 120s cadence. The uncertainty of the normalized flux in this light curve is 0.002, which motivated the search for transit signals of small planets (see Sect. 4.4).

As a nearby ultracool dwarf, Teegarden’s Star has been observed by the SPECULOOS (Search for habitable Planets EClipsing ULtra-cOOl Stars) project, which aims at finding transiting planets around such stars (Burdanov et al. 2018; Delrez et al. 2018; Sebastian et al. 2021). In particular, Teegarden’s Star was observed by SPECULOOS telescopes located in the northern hemisphere, namely Artemis (Burdanov et al. 2022) and SAINT-EX (Demory et al. 2020), which both can reach photometric precisions for a few minutes sampling of ~0.1% and mid-transit times accuracies of ~1 min (see e.g., Delrez et al. 2022; Pozuelos et al. 2023). In total, Teegarden’s Star was observed for 142 h over 35 nights between August 2021 and December 2022 (see Sect. 5.4).

3 Methods

The RVs were analyzed using a multi-planetary model employing Keplerian orbits. The Keplerian model was parameterized with the semi-amplitude K, orbital period and instead of eccentricity e and the argument of periastron ω, and the mean longitude λ, while the inclination i and the longitude of the ascending node Ω remain undetermined in the Keplerian model. Without limiting the generality, we therefore set i = 90° and Ω = 0°. This model was implemented in a series of scripts in python language.

As a cross-check, we used the versatile modeling tool juliet (Espinoza et al. 2019), which implements RV fits using RadVel (Fulton et al. 2018), a python package for modeling Keplerian orbits in RV time series. In addition to the aforementioned parameters (K, P, h, and k), the Keplerian model in juliet was parameterized with the time of periastron passage T, since juliet currently does not support the mean longitude λ as a parameter. We compared the results and the estimated values of the logarithm of the Bayesian evidence (ln 𝒵) between the models calculated with juliet and those previously described in order to ensure that the results were consistent. We also used Exo-Striker (Trifonov 2019), which allows efficient testing and visualization of modeling approaches.

We complemented the planetary models of varying complexity with a Gaussian process (GP) approach to model possible contributions from stellar activity, particularly in the form of rotational modulation. The posterior distributions of the parameters and hyper-parameters were determined with a nested sampling algorithm, which also provides the Bayesian evidence that was used to identify the best model. Finally, the parameters of the best model were checked for orbital stability following a dynamical analysis.

In an alternative approach, we used the ℓ₁-periodogram⁷ (Hara et al. 2017) to identify significant signals in the RV data. Based on an input frequency grid and a model for the covariance of the data, the ℓ₁-periodogram searches for a representation of the data with a small number of sinusoidal signals. The algorithm simultaneously tests all frequencies of the input grid. Following the tutorial from the git repository, we also tested the robustness of the noise model. A grid with various amplitudes for white noise, red noise, and correlated noise as well as periods and time scales for red noise and correlated noise provided a cross validation score for each covariance matrix.

3.1 Gaussian process

GP regression is a nonparametric Bayesian method used for modeling complex functions by assuming that the underlying function is a sample from a Gaussian distribution. The python package celerite (Foreman-Mackey et al. 2017) provides an efficient implementation of GP regression, specifically designed for large datasets, by modeling the covariance matrix as a sum of complex exponential functions. For modeling stellar activity, a simple damped harmonic oscillator (SHO) kernel is suitable (details are given in Eqs. (9), (20), and (21) of Foreman-Mackey et al. 2017) We use an improved kernel for rotation-modulated stellar activity composed of two SHO kernels (here called dSHO) from Foreman-Mackey (2018). The hyper-parameters for both SHO kernels include the scale of the activity induced noise σ₀ for each instrument, the quality factor of the oscillator Q₀, the scale ratio between the two oscillators f, the difference of the quality factors dQ, and the period of the underlying process P. The period of the second oscillator was fixed to the second harmonic of the first in order to also account for stellar activity signals at half of the rotation period. Additional white noise was incorporated using a jitter term. Following the discussion of Blunt et al. (2023) about potential overfitting with GP kernels of too much flexibility, we also tested a dSHO kernel, where the damping time scale of the two oscillators is forced to be identical. We also tested the effect of a joint covariance matrix for all RV datasets using their modified version of RadVel. The impact on the planet parameters of both is negligible. This is not unexpected since the correlated noise impact is small (see below).

3.2 Nested sampling

Nested sampling is a Bayesian inference method for estimating the logarithm of the evidence (ln 𝒵) and the posterior distributions of a model. The python package dynesty (Speagle 2020) is an implementation of nested sampling that is efficient in exploring complex posterior distributions and has a dynamic nested sampling algorithm. It provides diagnostic tools and options for controlling the sampling process. We used 5000 live points and stopped when δ ln 𝒵 < 0.01. We used the multi mode and the rwalk sampling. The priors for all of the parameters are listed in Tables 2 and G.1.

3.3 SPOCK

SPOCK (Stability of Planetary Orbital Configurations Klassifier; Tamayo et al. 2020) is a machine learning algorithm that was designed to predict the long-term stability of planetary systems. The algorithm is trained on a set of numerical simulations of planetary systems and uses a combination of physical and dynamical features to classify the stability of a given system. SPOCK has been used to predict the stability of observed exoplanetary systems (e.g., Tamayo et al. 2020, 2021; Kaye et al. 2022), and has been shown to be highly 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 and validation of new exoplanets (Wittrock et al. 2023). We used SPOCK to check the dynamical stability of the orbital solutions (see Sect. 4.5).

4 Results

4.1 Stellar properties

The stellar atmospheric parameters (T_eff, log g, and [Fe/H]) of Teegarden’s Star in Table 3 were derived from the co-added CARMENES spectra with the STEPARSYN⁸ code (Tabernero et al. 2022) using the line list and model grid described by Marfil et al. (2021). This update led to small modifications compared to the values reported by Zechmeister et al. (2019). All planetary masses listed in this paper are derived relative to the updated stellar mass of Teegarden’s Star, and the uncertainties consider the error propagation from the stellar mass.

Based on spectroscopic indices (see Sect. 2.2), the rotation period of Teegarden’s Star was determined as 96 days by Lafarga et al. (2021). A similar period of 100 days was found in the HPF data by Terrien et al. (2022). A cluster analysis of all spectral activity indicators revealed 98 days as the most likely stellar rotation period (Kemmer et al. 2023), aligning well with the determinations made by Lafarga et al. (2021) and Terrien et al. (2022) using spectroscopic activity indicators and Zeeman signatures, respectively. This agrees with earlier estimates of the rotation period reported by Zechmeister et al. (2019) using results from West et al. (2015), Newton et al. (2017), and Jeffers et al. (2018). The long rotation period also matches well with the low stellar activity indicating an age of around 8 Gyr as discussed in more detail by Zechmeister et al. (2019).

In Fig. 1, we present the periodograms of the CARMENES RV data, along with selected spectroscopic indices. Notably, the potential rotation period at 96 days and its 1-yr alias at 79 days are evident in the RVs (first panel), the CRX (second panel), and dLW (third panel) indices. However, the spectroscopic indices and the RVs share more long-periods (320 days, 172 days, 120 days), with the most prominent one among them in the RV data occurring at 172 days.

The power observed at 320 days and 120 days can be attributed to spectral leakage from the window function (see the bottom panel in Fig. C.1). This spectral leakage also contributes weakly to the power at 96 days. Consequently, we may interpret the 172 days signal as the true rotation period, with the other long-period variabilities being considered as alias signals. However, such a long rotation period would make Teegarden’s Star an outlier in terms of slow rotation. Newton et al. (2018) list 281 M stars with rotation periods determined from photometric monitoring using MEarth South (Irwin et al. 2015), including 31 stars below 0.12 M_⊙. Only one of them has a rotation period longer than 172 days, five have rotation periods above 150days (see Fig. 4 in Newton et al. 2018). A similar result is presented by Shan et al. (2024). In their literature compilation complemented with 129 new determinations using photometric monitoring within the CARMENES project, there is only one star with comparably low mass and a rotation period of 178 days, two in the range of 150days, all other 259 M stars are reported with lower rotation periods. Therefore we conducted further investigations into the signal at 172 days.

First, we fit Keplerian orbits to the RV data and to the spectroscopic indices using a uniform prior for the orbital period in the range of 150 days to 190 days. The result is shown in Fig. D.2. In case that residual telluric contamination were the cause of the 172 days signal in the RV data, we would expect that the period and phase of the RV data and the telluric contamination index agree, which is not the case. As a result, we excluded contamination from residual telluric lines. Since the CARMENES data cover about 14 cycles of the 172 days period, the difference in period and phase seems to be robust.

The results also indicate that the signal might be related to stellar activity, since the dLW and CRX indices show signals at consistent periods (Fig. D.2). We therefore calculated the stacked Bayesian Lomb-Scargle periodogram (sBGLS; Mortier & Collier Cameron 2017) of the observed RVs, displayed in Fig. E.1 (top). The variability of the power could be an indication of non-coherence and, hence, possibly interpreted as caused by stellar activity. In order to check the impact of the irregular sampling and spectral leakage, we injected a coherent signal corresponding to a Keplerian orbit of 172 days using the parameters from Table G.2 into the residuals of our overall best-fit model (model E, see Table 4). This synthetic dataset represents a single coherent signal at 172 days with the noise and sampling characteristics of the observation. The sBGLS periodogram is shown in the bottom panel of Fig. E.1. The power of the 172 days signal is also variable, despite being coherent. Since both sBGLS periodograms show variable power at 172 days in the observed and simulated data, the power fluctuations are therefore attributed to variable spectral leakage rather than to non-coherence. The 172 days signal is therefore compatible with being coherent over 2500 days. It remains unclear whether and according to what mechanism the signal at 172 days is related to stellar activity or other subtle effects such as remaining sky emission or detector artifacts. In this study we keep 96 days as the rotation period and discuss a potential planetary origin for the 172 days signal in the following section.

Inspection of the TESS data does not provide information about the stellar rotation period. Using the Simple Aperture Photometry (SAP) data we cut out regions with obvious excess instrumental noise. The light curve of the sectors 43 and 70 each show a small trend of about 0.01 flux change over the sector length, in sector 43 it is a positive, in sector 70 a negative trend. No trend is visible in sector 71. The sectors are separated by 2 yr the gap is therefore too large to obtain a meaningful fit for a periodic variation. Within each sector, variability at the level of the rms on time scales below 10 days can be seen. Sector 43 additionally shows two small flares. The TESS light curves presents Teegarden’s Star as a photometrically very quite star. This well matches with the conclusion from Zechmeister et al. (2019), where occasional flares as well as an overall low activity level were concluded from Hα-variability.

The spectroscopic activity indices as well as the RV data show a trend at a timescale of the length of the dataset. Fitting a sinusoidal function reveals a period of at least the length of the dataset (~2500 days, Fig. F.1). Since the posterior distributions of the period and phase overlap, this may constitute an indication of an activity cycle of unknown duration.

Fig. 1 Periodogram of RV data and spectroscopic activity indicators from CARMENES VIS. From top to bottom: RV data, chromatic index (CRX), differential line width (dLW), Ha, and the telluric contamination index. The blue triangles indicate the rotation period at 96 days and its 1-yr alias at 79days, the red triangles indicate the periods of the other signals (4.9 days, 11.4days, 26 days, and 172 days). The dashed horizontal lines show the false alarm probability levels at 10%, 1% and 0.1%.

Fig. 2 Periodogram of the CARMENES VIS data. From top to bottom: signals are subsequently subtracted, with red color denoting potential planetary signals and blue indicating the rotation period. Model names correspond to those in Table 4. The 1 day alias period of the signals is marked with a triangle in lighter color. The dashed horizontal lines are the false alarm probabilities at 10%, 1% and 0.1%. From top to bottom: original data, model A, model D, model E, and model G.

4.2 Analysis of CARMENES data

The measurements from CARMENES VIS are the ones with the longest time baseline, in addition to having a quite dense sampling and relatively low uncertainties. We therefore first analyzed the CARMENES VIS data separately and later turned to the analysis of the combined dataset.

The periodogram of the RV data (Fig. 2, top panel) shows significant peaks from the two already reported planets, namely Teegarden’s Star b and c (Zechmeister et al. 2019). By subtracting these two signals (model A, see Table 4), the periodogram of the residuals shows the strongest power at 172 days (Fig. 2, second panel). The inspection of the periodogram of the spectroscopic activity indicators (Fig. 1, see Sect. 4.1) reveals that this period is present there as well. Nevertheless, we fit it with a Keplerian orbit (model B), which then leads to the next highest peak at 96days, which we relate to stellar rotation (Sect. 4.1). Subtracting this signal (model D) results in a peak at 26 days (Fig. 2, third panel). We interpret this as the signal of a third planet, Teegarden’s Stard. The Bayesian evidence in 𝒵 of model E increases by more than 5, which indicates that the more complex model is superior (Table 4, fourth column). We also tested whether adding the 26 days signal to model A (model C) results in an improved Bayesian evidence. The difference of Δ in 𝒵 = 3.8 indicates a moderate improvement.

We checked whether the five signals exhibit any crosstalk due to spectral leakage in the window function of the measurement time series. The potential spectral leakage can be investigated in Fig. C.1. The three planetary signals at 4.9 days, 11.4days, and 26 days are not affected by crosstalk with any other signals. We note that the 11.4 days signal produces some power close to, but not exactly at, 26 days. Crosstalk between the two signals is therefore not expected. It should be noted that the 172 days signal has a side lobe very close to 96 days. However, the removal of the 172 days signal eliminated the aliases at 320 days and 120 days but did not affect the signal at 96 days. As a result, these two periodicities at 96 days and the 172 days seem to be independent but possibly affected by spectral leakage.

After subtracting the five signals, no remaining peak reaches a power above the 10% false alarm probability level in the periodogram, but two signals appear to be stronger than the adjacent noise. The one at 7.7 days is especially suggestive because it lies close to a 3:2 period commensurability chain with the two known planets (at 4.9 days and 11.4 days). Again, the Bayesian evidence grows (model F), at least by more than 2.5, so that the model F including an additional planet candidate at 7.7 days is moderately preferred.

There are two more peaks of similar strength in the periodogram, one at 1.104days, the other its one-day alias at 10.6 days. If we associate the signal to a planet, an orbital period of 10.6 days can be ruled out since this would make the system dynamically unstable due to the close vicinity of planet c at 11.4 days. The 1.104 days peak would therefore correspond to the true signal if it were a planet. Including a Keplerian at 1.104 days (model G) instead of the 7.7 days (model F) would also moderately increase the Bayesian evidence compared to model E. The difference between the Bayesian evidences of the six-signal models is insignificant, but it is tempting to accept the 7.7 days candidate signal as real since it would make the Teegarden’s Star planetary system dynamically compact and remarkably similar to the TRAPPIST-1 system.

In addition to fitting the potentially activity related signals above ~ 70 days with Keplerian orbits, we further tried to use the GP dSHO kernel described in Sect. 3.1 to model the variability in that range. The Bayesian evidence difference between the two-planet model with GP (model H) and the corresponding four-signal model without GP (model D) as well as the three-planet model with GP (model I) compared to the five-signal model without GP (model E) results in a higher evidence for the models without GP, which would support the interpretation of all five signals as being associated with planetary orbits. However, due to the presence of the 96 days and 172 days signals in the activity indices (Sect. 4.1), we prefer not to make a claim for the presence of further planets in the system at this point. The signals potentially caused by stellar activity are better represented by variability coherent at least over the duration of the observations indicating a very stable activity pattern on the stellar surface. We also tested a model with four planets and stellar rotation at 96 days modeled with a dSHO kernel (model J). This is superior compared to one where all long-period signals are modeled by the dSHO kernel (model I). Using CARMENES data alone, model J is moderately preferred. It is worthwhile to mention that the resulting amplitudes of the GP (dSHO kernel with and without enforced identical damping time scales as well as a single oscillator SHO kernel) in model I and J are very close to the amplitudes of the corresponding Kepler models for the 172 days (1.3 m s⁻¹) and 96 days (1.0 m s⁻¹) signals, respectively. The high Q value corresponds to a correlation time scale of about 15 000days, which is larger than the length of the observations. The SHO or dSHO kernels therefore mimic a coherent Keplerian model.

Using the ℓ₁-periodogram (Hara et al. 2017) supports our results from the signal detection using Keplerian models of increasing complexity. The ℓ₁-periodogram with the logarithm of the Bayes factor is shown in Fig. B.1. The five signals of model E are also present in the ℓ₁-periodogram and indicated as significant according to their Bayes factor. The long-period signal interpreted as activity cycle (Sect. 4.1, Fig. F.1) is also present and significant in the ℓ₁-periodogram.

4.3 Analysis of the full dataset

In a second step of the analysis, we used the full dataset and reran the models from Sect. 4.2. We would like to point out that the datasets from ESPRESSO and MAROON-X contain only a relatively small number of measurements. Since the ESPRESSO data overlap with the CARMENES data, the instrumental offset can be rather well constrained. The MAROON-X data happen to correspond to epochs that fall within a long observing gap in the CARMENES time series. However, the HPF dataset overlaps with the CARMENES observations and with two of the three MAROON-X campaigns. The disadvantage of adding multiple short duration datasets, each with an individual instrumental offset, is therefore mitigated by mutual overlaps between the datasets except for the third campaign of MAROON-X.

We investigated the impact of the additional data on the planetary parameters. The results are illustrated in Fig. 3 for the orbital period and RV amplitude. We note the slight improvement of the parameter determination through lower uncertainties (e.g., the orbital period of planet b) and mutual consistency of the parameter determinations within the uncertainties. Without the CARMENES data (gray histogram), the signals of the known planets b and c can be detected, but with significantly larger uncertainties because the time coverage is much smaller. The signal at 26 days and 172 days cannot be detected without CARMENES. The ESPRESSO and MAROON-X observations are too short each to cover at least one period, the uncertainty of the mutual offsets therefore interferes with the detection of the signals. The scatter in the HPF data connecting them is too large to compensate that.

In Table 4 (fifth and sixth columns) the Bayesian evidences and the differences to our preferred model are listed. While the CARMENES data alone would provide statistically significant support for an additional planet between planet b and c at 7.7 days, the full dataset does not. The short period signal at 1.104 days is also insignificant in the full dataset. The CARMENES data also gives preference to a model with four planets and stellar rotation at 96 days modeled as dSHO (model J) compared to a three-planet model and stellar activity variability at periods above ~70 days (model I). This, however, is not the case for the full dataset.

The model containing five Keplerian signals (model E) is our preferred final best model. Three of the signals, those at 4.9 days, 11.4 days, and 26 days, are interpreted as arising from the Keplerian motions of planet b, c, and d in the system, while the signals at 96 days and 172 days are interpreted as likely to be caused by stellar activity. Nonetheless, we emphasize that the latter ones are best represented by coherent variability over the entire time baseline of our observations, thus indicating a long-lived stellar activity pattern. The inclusion of the 26 days signal (planet d) leads to an improved Bayesian evidence, independent of the sequence in which signals are added to the models (Table 4).

The adopted fits are presented in Figs. 4 and 5 over time in intervals of half a year and in orbital phase for each planetary signal, respectively. The parameters and their uncertainties as well as the priors are listed in Tables 2 and G.1 for the planetary and data parameters, respectively. Table 2 also contains a selection of derived planetary parameters, such as the minimum planetary mass, the equilibrium temperature (assuming an albedo of 0.3), and the instellation. In the Appendix, we also list the parameters for the additional two signals that are likely due to stellar activity but better represented by Keplerian functions (Table G.2). For the 172 days signal we also list a minimum planetary mass in case the further observations reveal a planetary origin of this signal (and it would become Teegarden’s Star e). Finally, the posterior distribution and correlations for all fit parameters are displayed in Figs. H.1–H.4.

As a final check, we followed the approach of Blunt et al. (2023) and tested the robustness of the final model against overfitting. 80% of the RV measurements of each dataset were used to train model E again. The remaining 20% of the data were kept as control set. Using three random selections of the training set the combined weighted rms of the residuals of the training and control set, both calculated with the best-fit model parameters derived from the training set, are nearly identical (1.79 m s⁻¹ versus 1.78 m s⁻¹). This indicates a robust model with reliable predictive properties. This is not unexpected since adding new RV data left all planet parameters unchanged as already shown in Fig. 3.

Fig. 3 Comparison of the posterior distributions of the three planets (planet b, c, and d from top to bottom) and the signal of an unclear nature at 172 days depending on the datasets used for the analysis.

Fig. 4 Full dataset overlayed with the best-fit model (model E, see Table 4). The error bars contain the jitter contribution listed in Table G.1.

Fig. 5 Phase-folded full dataset overlayed with the best-fit model (model E, see Table 4). From top to bottom: planet b, c, d. For each panel, the contribution from the other signals has been filtered out. The error bars contain the jitter contribution listed in Table G.1.

4.4 Transit search and detection limits

In this subsection, we aim to confirm or refute the transiting nature of the planets orbiting Teegarden’s Star. To this end, we explored the public data provided by the TESS mission, which observed this star in sectors 43, 70, and 71 during the first mission extension in September 2021 and September and October 2023 using the 2 min and 10 min cadences. For our study, we used only the data corresponding to the 2 min cadence.

Neither the Science Processing Operations Center (SPOC; Jenkins et al. 2016) nor the Quick-Look Pipeline (QLP; Huang et al. 2020) have issued an alert for Teegarden’s Star, which is a hint of the non-transiting nature of this system. However, these pipelines may not detect some periodic transits if they are shallow and their S/N are below their detection thresholds. Then, the community often conducts complementary planetary searches either using space- and ground-based telescopes (see e.g., Delrez et al. 2022; Peterson et al. 2023; Trifonov et al. 2021) or using alternative custom detection pipelines (see e.g., Montalto 2023).

In this context, we scrutinized the TESS data employing the SHERLOCK pipeline⁹, presented initially by Pozuelos et al. (2020) and Demory et al. (2020), and used in several studies (see e.g., Wells et al. 2021; Van Grootel et al. 2021; Schanche et al. 2022). This pipeline allows the user to recover known planets, identify planetary candidates, and to search for new periodic signals, which may hint at the existence of additional transiting planets. In short, SHERLOCK combines several modules to (i) download and prepare the light curves from the MAST repository, (ii) search for planetary candidates, (iii) perform a semi-automatic vetting, (iv) compute a statistical validation, (v) model the signals to refine their ephemerides, and (vi) compute observational windows for ground-based observatories to trigger a follow-up campaign. We refer the reader to Pozuelos et al. (2023) for a recent application and further details.

Our first investigation covered potential orbital periods from 0.5 day to 13 days. This range guarantees the occurrence of at least two transits for any transiting planet. Simultaneously, it provides ample opportunity to identify planets b and c, in the hypothetical case that they transited. In a second trial, we searched for single transit-like events, which covered the potential scenario where the outermost planet d transits. In neither of these two searches, we found any hint of a signal that might correspond to a transiting planet in the system. All the signals detected were attributable to systematics or noise.

In view of this negative result, we stressed the TESS data to test if the photometric precision is sufficient to detect the presence of planets b and c, or any other planets in close-in orbits with masses below the detectability limit of our RVs. For this purpose, we used the MATRIX¹⁰ code (Dévora-Pajares & Pozuelos 2022; Pozuelos et al. 2020) to search for small, nearby planets using an injection and retrieval test. The simulation was set up with 24 period bins between 0.5 day and 12 days, as well as 15 planetary radius bins between 0.5 R_⊕ and 1.3 R_⊕. Each synthetic planet in this period-radius parameter space was injected at four different phases. The data were detrended with Tukey’s biweight (or bisquare) method (Mosteller & Tukey 1977) with a 0.5 day window size. As a result, planets larger than 0.5 R_⊕ would be detected with a detection efficiency close to 100% in orbits shorter than 7 days (yellow area in Fig. I.1 and with an efficiency better than about 75% in the period range of 7 days to 12 days. Using the consecutive sectors 70 and 71 we conducted a similar test for 15 periods between 25 days and 27 days, 15 radius bins between 0.5 R_⊕ and 1.3 R_⊕ and 4 phase bins (Fig. I.2). Planet d would have been detected with high probablity if it would be a transiting planet.

A non-detection in the TESS data, therefore, excludes transiting planets with most of the parameters that we tested, especially small planets in close-in orbits. We should have found planet b at 4.9 days and c at 11.4 days orbital periods if they transited. Planet d with a 26 days orbital period would at most show one transit in a single TESS sector and up to two in the two consecutive ones. However, due to the lower geometric transit probability, planet d is unlikely to be transiting, while b and c are not. Combined together, these results confirm the non-transiting nature of any planet with an orbital period shorter than 12 days in the Teegarden’s Star system, in particular the planets b and c, supporting the preliminary results from ground-based data by Zechmeister et al. (2019).

4.5 Dynamical stability

We used SPOCK to check the resulting planetary system configurations for dynamical stability. We used the posterior distributions and not draws from the normal distributions given by the mean and standard deviation shown in Table 2. This analysis reveals that all parameter configurations drawn from the posteriors are dynamically stable. This is not surprising, since the planetary system is not dynamically compact as further discussed below (Sect. 5). In case of the confirmation of the 7.7 days signal (Sect. 4.2) as the orbital period of an additional planet, the dynamical stability of the planetary system would then require small eccentricities in order to prevent close encounters.

5 Discussion

5.1 Planet, stellar activity, or instrumental artifact

As already mentioned in Sect. 4.2, the signal at 96days is likely associated with stellar activity and we interpret it as the potential stellar rotation period (Lafarga et al. 2021). Also as described in Sect. 4, the nature of the 172 days signal is unclear. On the one hand, if the rotation period is 96 days, 172 days is sufficiently separated from twice the rotation period, and likewise if 96 days was assumed to be the 1-yr alias of a 79 days rotation period. Since the signal is coherent over the duration of the observations, the activity-induced variability would need to be very stable in time. Normally, a coherent signal without a harmonic close to the potential rotation period or to its alias would be an indication of a planetary signal. On the other hand, the existence of variability with very similar period in spectroscopic activity indices is a strong indicator of an intrinsic stellar effect rather than a Keplerian orbital motion. The wavelength dependence of the amplitude does not provide any further information to discern the signal’s nature. The long-duration datasets that expand the wavelength coverage, CARMENES NIR and HPF, have insufficient precision to determine the amplitude of the 172 days signal and investigate its potential wavelength dependence. The conservative interpretation for the 172 days signal therefore is an activity related signal, probably not connected to the stellar rotation period.

The period of planet d at 26.13 ± 0.04 days falls comfortably outside the lunar sidereal and synodic periods of 27.32 days and 29.53 days, respectively, which generate small peaks in the periodograms of the spectroscopic activity indicators (Kemmer et al. 2023). To investigate the occurrence of common periodicities detected in all activity indicators of CARMENES stars, Kemmer et al. use a clustering algorithm. A cluster close to the sidereal and synodic month was detected in 5 out of 136 stars. This could indicate that the preferential scheduling of CARMENES observations in bright or gray time imprints the lunar cycle into the window function of the sampling in some cases. The distinction between the period of planet d and the sidereal month assures that there are no concerns regarding the interpretation of the 26.13 days signal as planetary in nature. It also does not coincide with a harmonic of the rotation period nor is it affected by spectral leakage of one of the other signals (Fig. C.1). We therefore conclude that this signal is a bona fide planetary RV signal.

Fig. 6 Planetary system architectures in comparison. Left panel: semi-major axes in planetary systems in comparison. For Teegarden’s Star we show two possible configurations: the three planets as fit by our preferred model and a hypothetical configuration (labeled as Teegarden (h)) with the three planets, the candidate at 7.7 days (light blue) as well as a purely hypothetical planet between planets c and d (light gray). The size of the plot symbols in the left panel are scaled assuming a constant planet density using the masses from Tables 2 and Table G.2. Right panel: mutual Hill separations in planetary systems in comparison. The light blue line indicates the mutual Hill radius above which the orbital crossing time is larger than 109 orbits.

5.2 Habitability

With the minor adjustments in stellar and planetary parameters when compared to Zechmeister et al. (2019) and with a new planet in the system, it is worthwhile to reevaluate the habitability of these planets. Teegarden’s Star planet b continues to exhibit Earth-like characteristics, with an equilibrium temperature of 277 K, assuming an albedo of 0.3, and an instellation (1.08 S_⊙) very close to that of our Earth by the Sun. The Earth Similarity Index (ESI; Schulze-Makuch et al. 2011) has only slightly decreased to 0.90, no longer holding the highest ESI ranking according to the habitable exoplanet catalog¹¹. While planet b has experienced a marginal ESI downgrade, planet c now has an ESI of 0.88, closely resembling Proximab (Anglada-Escudé et al. 2016). In contrast, the newly discovered planet d is cold, residing on an orbit of about a month, resulting in temperatures akin to Jupiter or its icy moon Ganymede. The prospect of directly detecting the planet’s atmosphere, which may be feasible with future instruments given the relatively short distance to Teegarden’s Star, becomes profoundly intriguing. Such a discovery would offer the opportunity to investigate the atmospheres of small rocky planets across a wide range of surface temperatures.

5.3 Planetary system architecture

We now compare the planetary system architecture of Teegarden’s Star to that of similar planetary systems. From the NASA Exoplanet archive¹² we extracted all systems fulfilling the following criteria: (i) more than one known planet, (ii) host star with M_star < 0.2 Μ_⊙, (iii) planet mass measured from RVs or transit timing variations (i.e., dynamical mass), and (iv) at least one potential rocky planet M_planet < 2 Μ_⊕. This yielded a total of seven systems including Teegarden’s Star, namely the following sorted by decreasing host star mass (taken from the references of the default parameter set as listed in the exoplanet archive): LHS 1140 (c[+b]; Lillo-Box et al. 2020), GJ 1132 (b[+c]; Bonfils et al. 2018), YZ Cet (bcd; Stock et al. 2020), GJ 1002 (bc; Suárez Mascareño et al. 2023), GJ 1061 (bcd; Dreizler et al. 2020), and TRAPPIST-1 (bcdefgh; Agol et al. 2021).

The planetary system architectures are displayed in Fig. 6, where the size of the plot symbols scales with the planet radius derived assuming Earth-like bulk density. We note that the minimum planet masses are used except for TRAPPIST-1. These systems are all compact and their planets orbit within 0.1 au. Only the two with the highest host star mass, LHS 1140 and GJ 1132, have planets above M_planet > 2 Μ_⊕. A good indicator of the dynamic compactness of a planetary system is the mutual Hill radius also termed fractional orbital separation (Gladman 1993). Following Obertas et al. (2017), it is calculated from the stellar mass M_star = M, the planetary masses M_planet = m, and semi-major axes a of two neighboring planets j and j + 1 as

Using the mutual Hill radii (Fig. 6) also shows that the LHS 1140 and GJ 1132 systems seem different compared to the ones with lower host star mass. The planets of systems with more massive planets have mutual Hill radius separations above Δ ≳ 30 while the ones with low-mass, likely terrestrial, planets have Δ ≲ 20.

The mutual Hill separations can be compared to those derived from multi-planet systems detected by Kepler (e.g., Weiss et al. 2018, their Fig. 14). The median mutual Hill separation of two-planet Kepler systems, whose hosts are predominantly FGK stars, is about 25. The planets of LHS 1140 and GJ 1132 fall within this distribution. The results by Weiss et al. (2018) also indicate that the mutual Hill radii tend to be smaller if more planets are present, but only about 10 % of the systems have a Δ < 10. While the stability of two-planet systems can be described analytically (Hill 1878; Gladman 1993), the investigation of the orbital stability of systems with more planets requires extensive numerical simulations (e.g., Chambers et al. 1996; Smith & Lissauer 2009; Ormel et al. 2017; Gratia & Lissauer 2021; Rice & Steffen 2023). A common result is that the logarithm of the time of the first orbital crossing depends linearly on the mutual Hill radius. At about 13 mutual Hill radii, the orbital crossing time scale seems to be longer than at least 10⁹ orbits for all systems. Below this limit the orbital crossing time scale changes by more than one order of magnitude for small changes in mutual Hill separation. These fluctuations are due to first and second order period commensurabilities, with additional dependencies on eccentricities and planet mass ratios.

The mutual Hill radii (Fig. 6), however, reveal that there seem to be two types of configurations. All systems except TRAPPIST-1 and YZCeti have mutual Hill radius separation above the threshold Δ ≳ 13 derived from numerical simulations where the orbit crossing time scale is probably at least as long as the system age. The planetary system of Teegarden’s Star, as obtained from our analysis with three planets, is well above this threshold for the two adjacent planet pairs (Teegarden label in Fig. 6). A putative planet at about 7.7 days would, however, change this picture drastically. Planets b and c would then form pairs with this potential planet e with mutual Hill separations only slightly above 10, very similar to planet pairs in the TRAPPIST-1 system (Teegarden (h) in Fig. 6). The signal that we investigated in Sect. 4.2 at 7.7 days would correspond to a planet below 0.5 M_⊕ and an RV amplitude of about 0.5 m s⁻¹. Given the measurement uncertainties, such a low-amplitude signal is at or below our detection limit but not at all unrealistic given the masses of the TRAPPIST-1 planets, which range from 0.33 M_⊕ to 1.37 M_⊕ (Agol et al. 2021). A hypothetical low-mass planet between planets c and d at about 17 days would be even more difficult to detect, but would turn the planetary system of Teegarden’s Star into a closer, and therefore brighter, twin of the TRAPPIST-1 system, with similar very low stellar mass, and multiple low-mass planets in a tightly packed system. More stringent upper limits on missing planets could, on the other hand, help to better establish the alternative, namely a rather loosely packed architecture of the planetary system of Teegarden’s Star.

5.4 Stellar flares observed using SPECULOOS telescopes

A detailed inspection of the SPECULOOS data confirmed the non-transiting nature of the system, supporting the results obtained in Sect. 4.4 using TESS data. However, these observations revealed several stellar flares. We identified 13 stellar flares by eye: 10 clear, single-peak flares and one multi-flare region with 3 peaks. There was an additional large flare with a slow, extended tail that could not be satisfactorily fit with the Davenport (2016) flare model or combination of flares and was not considered in this flare study. To obtain accurate flare energies we used PYMC3’s MCMC sampler (Salvatier et al. 2016) and xoflares (Gilbert et al. 2022; Davenport 2016) to model the flares simultaneously with linear systematic models for the full width half-maximum (FWHM), change in x and y CCD positions, sky background and airmass. We followed the method outlined in Murray et al. (2022) to calculate the flaring rates and U-band energies and generate the flare frequency diagram (FFD) for Teegarden’s Star. We do not account for the sensitivity of our flare detection as we are considering a single object. The reduced detection sensitivity at low flare energies will result in a tail-off (or flattening of the FFD) at low energies, however, it is not visible in the FFD in this case. The power law dN(E) = kE^−adE (Gershberg 1972; Lacy et al. 1976; Hawley et al. 2014) was fit to the FFD using PYMC3, where N is the flare occurrence rate, E is the flare energy, and k and α are constants. We obtained a best-fit α of 1.84±0.05, consistent with several sample studies of mid-to-late M dwarfs (Gizis et al. 2017; Paudel et al. 2018; Raetz et al. 2020; Murray et al. 2022). The FFD indicates that Teegarden’s Star may produce flares energetic enough for a planet with 1.08 Earth irradiation, such as planet b (Table 2), to enter the abiogenesis zone defined by Rimmer et al. (2018). On average, flare energies ≳10³⁵ erg will provide enough UV energy to build up a prebiotic inventory for RNA synthesis on planet b and, by extrapolating from SPECULOOS’s low energy flare sample, these type of flares should occur at most once every 2.4 yr.

6 Summary

Through a reanalysis of 346 nightly binned RV data collected from CARMENES, ESPRESSO, MAROON-X, and HPF, we report the discovery of a third planet in the Teegarden’s Star system. This newly identified planet, Teegarden’s Star d, has an orbital period of 26.13 + 0.04 days. The RV amplitude of 0.86 + 0.17 m s⁻¹ corresponds to a minimum mass of 0.82 M_⊕. However, due to the very low mass of the primary star, this planet orbits outside the habitable zone of Teegarden’s Star.

We have also observed an additional signal at 96 days, which is in good agreement with the stellar rotation period previously reported by Lafarga et al. (2021). The nature of another signal at 172 days, with an amplitude of 1.3 m s⁻¹, remains unclear. It can be attributed to stellar activity on a timescale longer than the rotation period or potentially indicate the presence of a long-period planet with a minimum mass of 2.3 M_⊕. However, the coincidence of the period with spectroscopic indicators makes the planet hypothesis less likely. Moreover, we found evidence of very long-period variability, spanning around 2500 days or longer, in the RV data and spectroscopic activity indicators, which may indicate the presence of an activity cycle.

In the CARMENES data, two additional signals of modest significance at a level of 0.5 m s⁻¹ were also identified. However, they are not detectable in the full dataset. One of these signals, with a period of 7.7 days, would correspond to a planet with approximately half the mass of Earth, situated close to a 3:2 commensurability chain with planets b and c. If confirmed, this planet would make the Teegarden’s Star planetary system tightly packed, resembling the TRAPPIST-1 system. A thorough search for yet undetected sub-Earth mass planets is crucial to elucidate the architectural characteristics of nearby planetary systems, particularly those suitable for atmospheric characterization.

Acknowledgements

We thank the referee for constructive and timely comments. We acknowledge the support from the Deutsche Forschungsgemeinschaft (DFG) under Research Unit FOR 2544 “Blue Planets around Red Stars” through project DR 281/32-1. This publication was based on observations collected under the CARMENES Legacy+ project. CARMENES is an instrument at the Centro Astronómico Hispano en Andalucía (CAHA) at Calar Alto (Almería, Spain), operated jointly by the Junta de Andalucía and the Instituto de Astrofísica de Andalucía (CSIC). CARMENES was funded by the Max-Planck-Gesellschaft (MPG), the Consejo Superior de Investigaciones Científicas (CSIC), the Ministerio de Economía y Competitividad (MINECO) and the European Regional Development Fund (ERDF) through projects FICTS-2011-02, ICTS-2017-07-CAHA-4, and CAHA16-CE-3978, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l’Espai, Institut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the MINECO, the DFG through the Major Research Instrumentation Programme and Research Unit FOR2544 “Blue Planets around Red Stars”, the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía. Based on observations obtained at the international Gemini Observatory, a program of NSF’s NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministèrio da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This work was enabled by observations made from the Gemini North telescope, located within the Maunakea Science Reserve and adjacent to the summit of Mauna Kea. We are grateful for the privilege of observing the Universe from a place that is unique in both its astronomical quality and its cultural significance. Based on observations obtained with the Hobby-Eberly Telescope (HET), which is a joint project of the University of Texas at Austin, the Pennsylvania State University, Ludwig-Maximillians-Universitaet München, and Georg-August Universitaet Göttingen. The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. We acknowledge the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (DR 281/37-1, HE 1935/28-1, QU 113/8-1, RE 1664/20-1). We acknowledge financial support from the Agencia Estatal de Investigación (AEI/10.13039/501100011033) of the Ministerio de Ciencia e Innovación and the ERDF “A way of making Europe” through projects PID2021-125627OB-C31, PID2019-109522GB-C5[1:4], PID2021-125627OB-C32, PID2022-137241NB-C43 and the Centre of Excellence “Severo Ochoa” and “María de Maeztu” awards to the Insti-tuto de Astrofísica de Canarias (CEX2019-000920-S), Instituto de Astrofísica de Andalucía (CEX2021-001131-S) and Institut de Ciències de l’Espai (CEX2020-001058-M). This work was also funded by the Generalitat de Catalunya/CERCA programme; the Universidad La Laguna and the EU Next Generation through the Margarita Salas Fellowship from the Spanish Ministerio de Universidades under grant UNI/551/2021-May-26; ERASMUS+, the European Union programme for education, training, youth and sport; the David and Lucile Packard Foundation, the Heising-Simons Foundation, and the Gordon and Betty Moore Foundation; the Gemini Observatory; the National Science Foundation under grant 2108465 and graduate research fellowship DGE 1746045; and NASA under grant 80NSSC22K0117 and the NASA Hubble Fellowship grant HST-HF2-51519.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. The ULiege’s contribution to SPECULOOS has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007–2013) (grant Agreement n° 336480/SPECULOOS), from the Balzan Prize and Francqui Foundations, from the Belgian Scientific Research Foundation (F.R.S.-FNRS; grant n° T.0109.20), from the University of Liege, and from the ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. MG is F.R.S-FNRS Research Director. The SPECULOOS-North team gratefully acknowledges financial support from the Heising-Simons Foundation, Dr. and Mrs. Colin Masson and Dr. Peter A. Gilman for Artemis, the first telescope of the SPECULOOS network situated in Tenerife, Spain. This work is supported by the Swiss National Science Foundation (PP00P2-163967, PP00P2-190080 and the National Centre for Competence in Research PlanetS). This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Appendix A Radial velocity data

Table A.1 lists the original RV data from CARMENES, ESPRESSO, MAROON-X, and HPF. We note that the CARMENES data slightly differ from Ribas et al. (2023), since we here use spectra cleaned from telluric contamination.

Appendix B Results of the ℓ₁-periodogram analysis

Fig. B.1 shows the ℓ₁-periodogram of the CARMENES VIS RV data using the best noise representation after cross validation of various covariance matrices. The cross validation results in no additional noise terms.

Fig. B.1 ℓ1-periodogram of the CARMENES VIS RV data. The detected signals are identified with the corresponding period in days (top) and with the logarithm of the Bayes factor (bottom).

Appendix C Crosstalk from signals and spectral leakage

Due to the complex window function of the time series of the CARMENES observations, we calculated the periodograms of the five signals with the parameters from Table 2 and Table G.2 with the sampling of the CARMENES observations. The result is depicted in Fig. C.1.

Fig. C.1 Periodogram of the Keplerian models with the parameters from Table 2. From top to bottom upper panel: 4.9 d, 11.4 d, 26 d, 96 d, and 172 d. The blue triangle indicates the rotation period at 96 d and its one-year alias at 79 d, the red triangles indicate the periods of the other signals (4.9 d, 11.4 d, 26 d, and 172 d). Lower panel: zoom into the long period range for the 96 d, and 172 d signals.

Appendix D Signal at 172 d in RV data and spectrospcopic indices

We show the comparison between the correlated noise model using draws from the posteriors of the hyperparameters for the dSHO kernel from model I in blue (Fig. D.1 abd the best-fit Keplerian model for the 172 d signal from model Ε (red). The data points are the offset-corrected data for all instruments, the error bars are adjusted to the fit jitter value. The contribution of planets b, c, and d have been subtracted. The correlated noise model shows a very long coherence and closely resembles the Keplerian model in amplitude and phase. This similarity of the signals can also be see comparing the hyperparameter in Table G.1 with the RV amplitude of the 172 d signal in Table G.2. The quality factor Q given in Table G.1 indicate a damping time scale on the order of twice the observing base line. If interpreted as spot life time, Teegarden’s Star would have a spot pattern stable for at least about 4000 d or 40 rotation periods.

In order to investigate whether the signal at 172 d in the RV data is due to stellar activity, instrumental effects, or a planet, we fit various datasets with identical priors and compared the output. In Fig. D.2 we show the posterior distributions for the periods and mean longitudes, which can be interpreted as phases in this context.

Fig. D.1 The noise model in frequency and time domain. Top: Periodogram of the correlated noise model of the best-fit of model I. The peri-odogram shows peaks both at 172 d (red) and the 96 d (blue) signal. Bottom: Comparison between the correlated noise model using draws from the posteriors of the hyperparameters for the dSHO kernel from model I (blue) and the best-fit Keplerian model for the 172 d signal from model E (red). The data points are the offset- and jitter-corrected measurements for all instruments where the contribution from planets b, c, and d have been subtracted.

Fig. D.2 Posterior distribution for a fit of a Keplerian model of the CARMENES RV data (RV) as well as spectroscopic indices: Telluric line contamination (telluric), differential line width (dLW), chromatic index (CRX), and Hα index (Hα). Posteriors for the orbital periods (left) and mean longitudes (right) are shown.

Appendix E Stacked Bayesian generalized Lomb-Scargle periodogram

In Fig. E.1 we compare the sBGLS periodogram of the RV data where the dominating signals of planets b and c were removed with a simulation of a coherent signal of 172 d period. A coherent 172 d signal with properties taken from Table G.2 was added to the residuals of the cArMENES RV data, where the five signals had been removed.

Fig. E.1 Stacked Bayesian GLS periodogram in the long-period range after subtraction of the signals from Teegarden’s Star b and c (model A, top) and of a simulated 172 d Keplerian orbit (bottom).

Appendix F Very long-period variations

The RV data and the spectroscopic activity indices show very long-period variations on the time scale of the total length of the dataset. We used identical priors (K: 𝒰[0, 3]m s⁻¹, P: 𝒰[500, 5000] d) to fit long-period signals and compare with the results illustrated by Fig. F.1. Given the time base of our datasets of 2596 d, the posterior distributions of the periods are consistent with each other. The phase shift between activity induced radial velocities and activity indicators is known for example from solar observations (Collier Cameron et al. 2019).

Fig. F.1 Posterior distributions for the orbital periods (left) and mean longitudes (right) for a fit of a Keplerian model of the CARMENES spectroscopic indices.

Appendix G Additional tables

In Table G.1 we list the priors and posteriors of the nine RV datasets, the jitter and offset parameters for each instrument as well as a linear and quadratic term. The two signals in the RV data, which are likely not due to orbiting planets, are listed in Table G.2.

Appendix H Corner plot

The corner plot of the best model (model E from Table 4) is shown in Fig. H.1. It depicts the posteriors and the covariances of all parameters used. The results are also listed in Tables 2, G.1, and G.2.

Fig. H.1 Corner plot of all fit parameters. From left to right: Jitter terms for all instruments, Keplerian parameters for 5 signals, RV offsets for all instruments as well as the linear and quadratic trend terms.

Fig. H.2 Corner plot of parameters for planet b.

Fig. H.3 Corner plot of parameters for planet c.

Fig. H.4 Corner plot of parameters for planet d.

Appendix I TESS transit detection

Here we present the result of the planet injection test using MATRIX¹³ for TESS sectors 43, 70, and 71. Planets with orbital periods between 1 d and 12 d and radii between 0.5 R_⊕ and 1.3 R_⊕ were injected in TESS sector 43 with four different phase values. The recovery fraction is shown in Fig. I.1. The potential planet labeled e and f corresponds to the signal at 7.2 d and 1.1 d (models F and G) discussed in Section 4.4.1. Using the two consecutive sectors 70 and 71, we tested the detectability of planet d injecting planets with periods between 25 d and 27 d and radii between 0.5 R_⊕ and 1.3 R_⊕ (Fig. I.2). The planet radii shown in the two figures assume Earth density, the assumed masses are the minimum masses and the corresponding uncertainty.

Fig. I.1 Planet injection and recovery test for the planets and potential planet candidates with periods below 12 d using MATRIX and TESS sector 43.

Fig. I.2 Planet injection and recovery test for the planet d using MATRIX and TESS sectors 70 and 71.

References

All Tables

All Figures

Fig. 1 Periodogram of RV data and spectroscopic activity indicators from CARMENES VIS. From top to bottom: RV data, chromatic index (CRX), differential line width (dLW), Ha, and the telluric contamination index. The blue triangles indicate the rotation period at 96 days and its 1-yr alias at 79days, the red triangles indicate the periods of the other signals (4.9 days, 11.4days, 26 days, and 172 days). The dashed horizontal lines show the false alarm probability levels at 10%, 1% and 0.1%.

In the text

Fig. 2 Periodogram of the CARMENES VIS data. From top to bottom: signals are subsequently subtracted, with red color denoting potential planetary signals and blue indicating the rotation period. Model names correspond to those in Table 4. The 1 day alias period of the signals is marked with a triangle in lighter color. The dashed horizontal lines are the false alarm probabilities at 10%, 1% and 0.1%. From top to bottom: original data, model A, model D, model E, and model G.

In the text

	Fig. 3 Comparison of the posterior distributions of the three planets (planet b, c, and d from top to bottom) and the signal of an unclear nature at 172 days depending on the datasets used for the analysis.
In the text

	Fig. 4 Full dataset overlayed with the best-fit model (model E, see Table 4). The error bars contain the jitter contribution listed in Table G.1.
In the text

	Fig. 5 Phase-folded full dataset overlayed with the best-fit model (model E, see Table 4). From top to bottom: planet b, c, d. For each panel, the contribution from the other signals has been filtered out. The error bars contain the jitter contribution listed in Table G.1.
In the text

Fig. 6 Planetary system architectures in comparison. Left panel: semi-major axes in planetary systems in comparison. For Teegarden’s Star we show two possible configurations: the three planets as fit by our preferred model and a hypothetical configuration (labeled as Teegarden (h)) with the three planets, the candidate at 7.7 days (light blue) as well as a purely hypothetical planet between planets c and d (light gray). The size of the plot symbols in the left panel are scaled assuming a constant planet density using the masses from Tables 2 and Table G.2. Right panel: mutual Hill separations in planetary systems in comparison. The light blue line indicates the mutual Hill radius above which the orbital crossing time is larger than 109 orbits.

In the text

	Fig. B.1 ℓ1-periodogram of the CARMENES VIS RV data. The detected signals are identified with the corresponding period in days (top) and with the logarithm of the Bayes factor (bottom).
In the text

Fig. C.1 Periodogram of the Keplerian models with the parameters from Table 2. From top to bottom upper panel: 4.9 d, 11.4 d, 26 d, 96 d, and 172 d. The blue triangle indicates the rotation period at 96 d and its one-year alias at 79 d, the red triangles indicate the periods of the other signals (4.9 d, 11.4 d, 26 d, and 172 d). Lower panel: zoom into the long period range for the 96 d, and 172 d signals.

In the text

Fig. D.1 The noise model in frequency and time domain. Top: Periodogram of the correlated noise model of the best-fit of model I. The peri-odogram shows peaks both at 172 d (red) and the 96 d (blue) signal. Bottom: Comparison between the correlated noise model using draws from the posteriors of the hyperparameters for the dSHO kernel from model I (blue) and the best-fit Keplerian model for the 172 d signal from model E (red). The data points are the offset- and jitter-corrected measurements for all instruments where the contribution from planets b, c, and d have been subtracted.

In the text

	Fig. D.2 Posterior distribution for a fit of a Keplerian model of the CARMENES RV data (RV) as well as spectroscopic indices: Telluric line contamination (telluric), differential line width (dLW), chromatic index (CRX), and Hα index (Hα). Posteriors for the orbital periods (left) and mean longitudes (right) are shown.
In the text

	Fig. E.1 Stacked Bayesian GLS periodogram in the long-period range after subtraction of the signals from Teegarden’s Star b and c (model A, top) and of a simulated 172 d Keplerian orbit (bottom).
In the text

	Fig. F.1 Posterior distributions for the orbital periods (left) and mean longitudes (right) for a fit of a Keplerian model of the CARMENES spectroscopic indices.
In the text

	Fig. H.1 Corner plot of all fit parameters. From left to right: Jitter terms for all instruments, Keplerian parameters for 5 signals, RV offsets for all instruments as well as the linear and quadratic trend terms.
In the text

	Fig. H.2 Corner plot of parameters for planet b.
In the text

	Fig. H.3 Corner plot of parameters for planet c.
In the text

	Fig. H.4 Corner plot of parameters for planet d.
In the text

	Fig. I.1 Planet injection and recovery test for the planets and potential planet candidates with periods below 12 d using MATRIX and TESS sector 43.
In the text

	Fig. I.2 Planet injection and recovery test for the planet d using MATRIX and TESS sectors 70 and 71.
In the text