Teegarden’s Star revisited A nearby planetary system with at least three planets ⋆ , ⋆⋆

The two known planets in the planetary system of Teegarden’s Star are among the most Earth-like exoplanets currently known. Revis-iting this nearby planetary system with two planets in the habitable zone aims at a more complete census of planets around very low-mass stars. A significant number of new radial velocity measurements from CARMENES, ESPRESSO, MAROON-X


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, ⋆ Based on observations collected at the European Southern Observatory under ESO programme(s) 0103.C-0152(A).Table A.1 is only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgibin/qcat?J/A+A/.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 exoplan-ets1 .It is noteworthy that the majority of these planets have been found orbiting M dwarfs.Among them, the TRAPPIST-12 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)3 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 4 .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 25 th nearest star to the Sun5 (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 Sect. 4 and Sect.5, respectively.

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 serval6 (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 nearinfrared (NIR, 0.96-1.71µm) wavelengths.The average resolving power for the two wavelength regions is R = 94 600 and R = 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 Notes. (a) For the number of measurements those of the original (o), the 5 σ clipped (c), and daily binned data (b) sets are listed.The rms was calculated using the residuals after subtracting the best-fit model.The internal precision (the median of the uncertainties) was measured after the binning.
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 (R = 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, fiberfed high-resolution (R = 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 fiberfed near-infrared spectrograph (0.9-1.8 µm) for the 10 m Hobby-Eberly Telescope (HET, Ramsey et al. 1998) with a resolution of R = 53 000 (Mahadevan et al. 2012).Between September 2019 and October 2021, a total of 143 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 three years 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.

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 therefor 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.

Photometric data
The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014Ricker et al. , 2015) ) observed Teegarden's Star in sectors 42 and 44 in full frame image mode and in sectors 43, 70, and 71 with 120 s 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 midtransit 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).

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 P, h = √ e sin ω and k = √ e cos ω 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 Z) 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 ℓ 1 -periodogram7 (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 ℓ 1 -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.

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 Equations 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 σ 0 for each instrument, the quality factor of the oscillator Q 0 , 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).

Nested sampling
Nested sampling is a Bayesian inference method for estimating the logarithm of the evidence (ln Z) 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 Z < 0.01.We used the multi mode and the rwalk sampling.The priors for all of the parameters are listed in Tables 4 and G.1.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. 2020Tamayo et al. , 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).

Stellar properties
The stellar atmospheric parameters (T eff , log g, and [Fe/H]) of Teegarden's Star in Table 2 were derived from the co-added CARMENES spectra with the SteParSyn8 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 d by Lafarga et al. (2021).A similar period of 100 d was found in the HPF data by Terrien et al. (2022).A cluster analysis of all spectral activity indicators revealed 98 d 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), andJeffers 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 d and its 1-year alias at 79 d 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 d, 172 d, 120 d), with the most prominent one among them in the RV data occurring at 172 d.
The power observed at 320 d and 120 d 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 d.Consequently, we may interpret the 172 d 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 d, five have rotation periods above 150 d (see Fig. 4 Newton et al. 2018).A similar result is presented by Shan et al. (2024).In their literature compilation complemented with 129 new determinations using 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 d to 190 d.The result is shown in Fig. D.2.In case that residual telluric contamination were the cause of the 172 d 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 d 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 d using the parameters from Table G.2 into the residuals of our overall best-fit model (model E, see Table 3).This synthetic dataset represents a single coherent signal at 172 d 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 d signal is also variable, despite being coherent.Since both sBGLS periodograms show variable power at 172 d in the observed and simulated data, the power fluctuations are therefore attributed to variable spectral leakage rather than to noncoherence.The 172 d signal is therefore compatible with being coherent over 2 500 d.It remains unclear whether and according to what mechanism the signal at 172 d is related to stellar activity or other subtle effects such as remaining sky emission or detector artifacts.In this study we keep 96 d as the rotation period and discuss a potential planetary origin for the 172 d 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 two years 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 d 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 (∼ 2 500 d, 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.

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  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 d and 11.4 d).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.104 d, the other its one-day alias at 10.6 d.If we associate the signal to a planet, an orbital period of 10.6 d can be ruled out since this would make the system dynamically unstable due to the close vicinity of planet c at 11.4 d.The 1.104 days peak would therefore correspond to the true signal if it were a planet.Including a Keplerian at 1.104 d (model G) instead of the 7.7 d (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 d 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 d 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 twoplanet model with GP (model H) and the corresponding foursignal 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 d and 172 d 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 d 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 d (1.3 m s −1 ) and 96 d (1.0 m s −1 ) signals, respectively.The high Q value corresponds to a correlation time scale of about 15 000 d, which is larger that the length of the observations.The SHO or dSHO kernels therefore mimic a coherent Keplerian model.
Using the ℓ 1 -periodogram (Hara et al. 2017) supports our results from the signal detection using Keplerian models of increasing complexity.The ℓ 1 -periodogram with the logarithm of the Bayes factor is shown in

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 mul- tiple 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 d and 172 d 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 3 (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 d, the full dataset does not.The short period signal at 1.104 d is also insignificant in the full dataset.The CARMENES data also gives preference to a model with four planets and stellar rotation at 96 d modeled as dSHO (model J) compared to a three-planet model and stellar activity variability at periods above ∼ 70 d (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 d, 11.4 d, and 26 d, are interpreted as arising from the Keplerian motions of planet b, c, and d in the system, while the signals at 96 d and 172 d 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 d signal (planet d) leads to an improved Bayesian evidence, independent of the sequence in which signals are added to the models (Table 3).
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 4 and G.1 for the planetary and data parameters, respectively.Table 4 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 d 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 to 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 −1 versus 1.78 m s −1 ).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.

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).[11.37, 11.44]  [25.70, 26.30]In this context, we scrutinized the TESS data employing the SHERLOCK pipeline9 , 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 d to 13 d.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 MATRIX10 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 d and 12 d, 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 d 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 d (yellow area in Fig. I.1) and with an efficiency better than about 75% in the period range of 7 d to 12 d.Using the consecutive sectors 70 and 71 we conducted a similar test for 15 periods between 25 d and 27 d, 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 d and c at 11.4 d orbital periods if they transited.Planet d with a 26 d 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 d 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).

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 4.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 d 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.

Planet, stellar activity, or instrumental artifact
As already mentioned in Sect.4.2, the signal at 96 d 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 d signal is unclear.On the one hand, if the rotation period is 96 d, 172 d is sufficiently separated from twice the rotation period, and likewise if 96 d was assumed to be the 1-year alias of a 79 d 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 d signal and investigate its potential wavelength dependence.The conservative interpretation for the 172 d 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 d falls comfortably outside the lunar sidereal and synodic periods of 27.32 d and 29.53 d, 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 d 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.

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 habit-ability 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 catalog11 .While planet b has experienced a marginal ESI downgrade, planet c now has an ESI of 0.88, closely resembling Proxima b (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.

Planetary system architecture
We now compare the planetary system architecture of Teegarden's Star to that of similar planetary systems.From the NASA Exoplanet archive12 we extracted all systems fulfilling the following criteria: (i) more than one known planet, (ii) host star with M star < 0.2 M ⊙ , (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 M ⊕ .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 M ⊕ .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 9 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 YZ Ceti 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 d 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 d would correspond to a planet below 0.5 M ⊕ and an RV amplitude of about 0.5 m s −1 .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 d 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 lowmass 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.

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 Section 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 Uband 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 −α dE (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 4), to enter the abiogenesis zone defined by Rimmer et al. (2018).On average, flare energies ≳ 10 35 ergs 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 years.

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 d.The RV amplitude of 0.86 ± 0.17 m s −1 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 d, which is in good agreement with the stellar rotation period previously reported by Lafarga et al. (2021).The nature of another signal at 172 d, with an amplitude of 1.3 m s −1 , remains unclear.It can be attributed to stellar activity on a timescale longer than the rotation period or potentially indicate the presence of a longperiod 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 2 500 d 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 −1 were also identified.However, they are not detectable in the full dataset.One of these signals, with a period of 7.7 d, 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.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

Fig. 1 .
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), Hα, and the telluric contamination index.The blue triangles indicate the rotation period at 96 days and its 1-year alias at 79 days, the red triangles indicate the periods of the other signals (4.9 days, 11.4 days, 26 days, and 172 days).The dashed horizontal lines show the false alarm probability levels at 10%, 1% and 0.1%.

Fig. 2 .
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 3.The 1 d 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.
Teegarden's Star b and c(Zechmeister et al. 2019).By subtracting these two signals (model A, see Table3), the periodogram of the residuals shows the strongest power at 172 d (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 96 d, which we relate to stellar rotation (Sect.4.1).Subtracting this signal (model D) results in a peak at 26 d (Fig.2, third panel).We interpret this as the signal of a third planet, Teegarden's Star d.The Bayesian evidence ln Z of model E increases by more than 5, which indicates that the more complex model is superior (Table3, fourth column).We also tested whether adding the 26 d signal to model A (model C) results in an improved Bayesian evidence.The difference of ∆ ln Z = 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 d, 11.4 d, and 26 d are not affected by crosstalk with any other signals.We note that the 11.4 d signal produces some power close to, but not exactly at, 26 d.Crosstalk between the two signals is therefore not expected.It should be noted that the 172 d signal has a side lobe very close to 96 d.However, the removal of the 172 d signal eliminated the aliases at 320 d and 120 d but did not affect the signal at 96 d.As a result, these two periodicities at 96 d and the 172 d seem to be independent but possibly affected by spectral leakage.
Fig. B.1.The five signals of model E are also present in the ℓ 1 -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 ℓ 1 -periodogram.

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 d depending on the datasets used for the analysis.

Fig. 6 .
Fig. 5. Phase-folded full dataset overlayed with the best-fit model (model E, see Table 3).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.

Table 1 .
List of RV datasets(a).

Table 3 .
The 1 d 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.

Table 3 .
Model selection based on Bayesian evidence(a).Notes.(a)Columnsdenote the model name, a short description, the Bayesian evidence for the model using CARMENES VIS data only as well as for the full dataset, and the difference with respect to our best model.

Table 4 .
Fit and derived planetary parameters (model E).