The CARMENES search for exoplanets around M dwarfs A sub-Neptunian mass planet in the habitable zone of HN Lib

We report the discovery of HNLibb, a sub-Neptunian mass planet orbiting the nearby ( d ≈ = 6 . 25pc) M4.0V star HNLib detected by our CARMENES radial-velocity (RV) survey. We determined a planetary minimum mass of M b sin i = 5 . 46 ± 0 . 75 M ⊕ and an orbital period of P b = 36 . 116 ± 0 . 029d, using ∼ 5yr of CARMENES data, as well as archival RVs from HARPS and HIRES spanning more than 13yr. The ﬂux received by the planet equals half the instellation on Earth, which places it in the middle of the conservative habitable zone (HZ) of its host star. The RV data show evidence for another planet candidate with M [c] sin i = 9 . 7 ± 1 . 9 M ⊕ and P [c] = 113 . 46 ± 0 . 20d. The long-term stability of the signal and the fact that the best model for our data is a two-planet model with an independent activity component stand as strong arguments for establishing a planetary origin. However, we cannot rule out stellar activity due to its proximity to the rotation period of HNLib, which we measured using CARMENES activity indicators and photometric data from a ground-based multi-site campaign as well as archival data. The discovery adds HNLibb to the shortlist of super-Earth planets in the habitable zone of M dwarfs, but HNLib[c] probably cannot be inhabited because, if conﬁrmed, it would most likely be an icy giant.


Introduction
In the last decade, spectrographs delivering high-precision radial velocity (RV) measurements have reached the necessary precision to detect small planets close to (and even within) the habitable zone (HZ; e.g., Kasting et al. 1993;Kopparapu et al. 2014) of late-type main-sequence stars. The predicted high number of super-Earth and Earth-like planets in the HZ that could be detected around M-dwarf stars makes these stars extremely inter-Send offprint requests to: Esther González-Álvarez and Jonas Kemmer e-mail: estgon11@ucm.es, jkemmer@lsw.uni-heidelberg.de ⋆ Both authors contributed equally to the manuscript. esting objects for planetary RV searches (Udry et al. 2007;Dressing & Charbonneau 2013;Bonfils et al. 2013;Tuomi et al. 2014;Sabotta et al. 2021). Despite the significant number of detected planets around M-dwarf primaries, we are still far from fully understanding fundamental questions, such as how these planetary systems form and evolve (Ida & Lin 2005;Raymond et al. 2007;Burn et al. 2021;Schlecker et al. 2022). Especially for the planets in the HZ, the formation history is important. This is because not only is the incoming stellar flux a crucial component in determining the actual habitability, but so are the properties of the planet itself, such as the composition of a potential atmosphere (Tarter et al. 2007;Kaltenegger 2017 In this context, the increasing number of M-dwarf stars that harbor planets in the HZ (e.g., Suárez Mascareño et al. 2023 andKossakowski et al. 2023, to cite the two most recent examples based on the radial velocity method) plays a significant role in exploring theories of planet formation and evolution.
Here, we present the discovery of a sub-Neptunian mass planet orbiting in the HZ of the nearby mid-M dwarf HN Lib. The planet's orbital parameters are derived from an RV model that includes Keplerian orbits for the planet and another planet candidate, as well as a red-noise model driven by a Gaussian process (GP) to simultaneously account for the stellar activity.
In Sect. 2, we present the analyzed ground-based observations. The properties of the host star, HN Lib, are introduced in Sect. 3. The stellar activity analysis and the determination of the rotation period of the parent star are carried out in Sect. 4. Also in this section, with the stellar rotation period constrained, we proceed to analyze the spectroscopic observations with the aim of identifying and characterizing planet candidates. The properties of the newly confirmed and discovered planet and the planet candidate orbiting HN Lib are given in Sect. 5, together with the discussion on the implications of these findings. Finally, a summary is given in Sect. 6.

Observations
2.1. Spectroscopic data 2.1.1. CARMENES HN Lib was spectroscopically observed with CARMENES between 27 January 2016 and 31 December 2020, resulting in a total of 94 RV measurements. CARMENES is installed at the 3.5 m telescope of the Calar Alto Observatory in Almería (Spain). It was specifically designed to deliver high-resolution spectra at optical (resolving power R ≈ 94,600) and nearinfrared (R ≈ 80,400) wavelengths covering the range from 520 nm to 1710 nm. CARMENES has two different channels, one for the optical (the VIS channel) and one for the nearinfrared (the NIR channel) with a dichroic at 960 nm (Quirrenbach et al. 2014). All data were acquired with integration times of 1800 s, which is the maximum exposure time employed for precise RV measurements with CARMENES.
The spectra were processed following the data flow of the CARMENES guaranteed time observations (GTO) program (Ribas et al. 2023). Raw data are reduced with the caracal pipeline (Caballero et al. 2016a). Relative RVs are extracted separately for the VIS and NIR channels using the serval software . The final RV per epoch of each channel is computed as the weighted RV mean over all échelle orders. The data from the CARMENES VIS channel are shown in the bottom panel of Fig. 1 (individual CARMENES relative RVs and their uncertainties are listed in Table C.1). The root mean square (rms) from the final fit residuals and the mean errorbars of the CARMENES RV data are 1.75 m s −1 and 1.31 m s −1 , respectively.
At high spectral resolution, the profile of the stellar lines may change due to photospheric and chromospheric activity, which has an impact on accurate RV measurements. Therefore, for CARMENES, serval provides further measurements for a number of spectral features in the instruments' wavelength range that are considered indicators of stellar activity and may have a chromospheric component in active M dwarfs: the differential line width (dLW), the Ca ii infrared triplet, Hα, and the chromatic index (CRX). The latter determines the RV-log λ correla-tion, and it is used as an indicator for the presence of stellar active regions. Additionally, as part of the data processing, we also calculate measurements of molecular absorption bands of two species: TiO and VO, together with pseudo-equivalent widths (pEW) for different indices after the subtraction of a reference star spectrum as described by Schöfer et al. (2019). An overview of the CARMENES activity indicators of HN Lib that exhibit significant peaks in the generalized Lomb-Scargle periodograms (GLS; ) is shown in Fig. 2.

Archival data
Additional relative RV data for HN Lib are available in the archives. We employed 14 RVs (Trifonov et al. 2020) from the High Accuracy Radial velocity Planet Searcher (HARPS, Mayor et al. 2003) at the ESO La Silla 3.6 m telescope and 34 RV data (Tal-Or et al. 2019) from the High Resolution Echelle Spectrometer (HIRES, Vogt et al. 1994) on the Keck 10 m telescope. These velocities span a period of ∼9 yr, almost doubling the total time coverage of our RV analysis. The mean error bars are 1.47 m s −1 for the HARPS data and 2.83 m s −1 for HIRES.

Ground-based photometry
There are data from many photometric monitoring campaigns of HN Lib available. An overview of the data used in this work and their properties can be found in Table 1.

OSN and TJO
HN Lib was observed by our team using the 90 cm Ritchey-Chrétien telescope of the Observatorio de Sierra Nevada (OSN) in Almería, Spain, and the VersArray 2k × 2k CCD camera (Rodríguez et al. 2010). The total number of frames amounts to 915 in the V filter and 903 in the R filter, where the typical exposure times were 100 s and 60 s, respectively. Each filter dataset was corrected for bias and flat fielding before the target flux was determined using differential aperture photometry. Outliers differing by more than 3σ from the mean were removed from further analysis. The nightly binned light curves are presented in the first two panels of Fig. 3.
Further photometric CCD observations were collected by us with the 0.8 m Telescopi Joan Oró (TJO, Colomé et al. 2010) at the Observatori Astronòmic del Montsec in Lleida, Spain. A total of 1910 frames with a exposure time of 25 s were obtained using the Johnson R filter of the LAIA imager, a 4k × 4k CCD with a field of view of 30 arcmin, and a scale of 0.4 arcsec pixel −1 . The images were calibrated with bias, dark, and flat-field frames with the ICAT pipeline (Colome & Ribas 2006) of the TJO. The differential photometry was extracted with AstroImageJ (Collins et al. 2017) using the aperture size that minimized the rms of the resulting relative fluxes, and a selection of the 30 brightest comparison stars in the field that did not show variability. Lastly, we removed outliers and measurements affected by poor observing conditions with a low signal-to-noise ratio. The nightly binned TJO R-band light curve is presented in the middle panel of Fig. 3.

Archival data
Photometric public data of HN Lib were also used for the analysis. We retrieved archival data from the MEarth project (Irwin et al. 2015, DR 11 1 -last two panels of Fig. 3 Pollacco et al. 2006), and the All-Sky Automated Survey for Supernovae project (ASAS-SN, Shappee et al. 2014;Kochanek et al. 2017). Díez Alonso et al. (2019) searched for rotation periods using the ASAS-3 photometric data but did not report a detection for HN Lib. To reach a higher precision for our analysis, we created nightly bins of the data from each instrument if multiple data points were taken per night.
Since the quality of the light curves of ASAS-3 and ASAS-SN is considerably worse than that of the others (cf . Table 1), we did not use them for our final analysis. For the MEarth data, we considered a photometric offset for each version of the respective instruments since they go along with the creation of a new perstar magnitude and meridian offset for the data (Newton et al. 2018). Some resulting sub-data sets have only a few data points and are not very useful due to the ambiguous offset in the fit, which is why we did not consider them in the final analysis (see also Table 1).

HN Lib
HN Lib is a nearby, single M4.0 V star located at a distance of only about 6.25 pc (Gaia Collaboration et al. 2022). Because of its relative brightness, this star has often appeared in the literature and in all-sky surveys (e.g. Schönfeld 1886; Jackson & Turner 1952;Leggett 1992;Mann et al. 2015). Table 2 summarises its stellar parameters.
HN Lib is cataloged as a variable star that exhibits variations in its luminosity due to rotation of the star coupled with star spots and other chromospheric activity (Kukarkin & Kholopov 1982;Weis 1994;Hosey et al. 2015). Furthermore, Astudillo-A&A proofs: manuscript no. HNLib  Defru et al. (2017) determined a log R ′ HK value for HN Lib of −5.58 ± 0.05 dex and estimated from activity-rotation relationships a stellar rotation period of 102 ± 11 d, which classifies our target as a slow rotator. This conclusion is in line with the small upper limits of rotational velocity, magnetic field, and X-ray luminosity (Stelzer et al. 2013;Reiners et al. 2018Reiners et al. , 2022, a new, more robust, mean log R ′ HK value of −5.23 ± 0.07 dex homogeneously computed by averaging individual values collected with HARPS, FEROS, HIRES, and ESPaDONS (Perdelwitz et al. 2021), the measurement of chromospheric lines in absorption or moderate emission (Fuhrmeister et al. 2020), and our own rotation period determination (see below).
The fundamental atmospheric parameters (i.e., T eff , log g, and [Fe/H]) of HN Lib using the CARMENES VIS and NIR template spectra were recently determined by Marfil et al. (2021) by means of the SteParSyn code 2 (Tabernero et al. 2022). The resulting values, listed in Table 2, agree at the 2σ level with those derived by other authors via spectral synthesis (e.g. Rojas-Ayala et al. 2012;Gaidos et al. 2015;Passegger et al. 2018;Rajpurohit et al. 2018;Maldonado et al. 2020) and from fits to the photometric spectral energy distribution (Cifuentes et al. 2020).
We followed the recipes of the latter authors to determine a new bolometric luminosity, L ⋆ , with updated Gaia DR3 data. We used L ⋆ and T eff to determine the mass and radius of HN Lib to be M ⋆ = 0.291 ± 0.014 M ⊙ and R ⋆ = 0.299 ± 0.009 R ⊙ with Stefan-Boltzmann law and the mass-radius relation of Schweitzer et al. (2019), which was based on detached, doublelined, double-eclipsing, main-sequence M-dwarf binaries from the literature.
The galactic space velocities UVW of HN Lib were derived using the Gaia DR3 coordinates, proper motions, and systemic radial velocity with the formulation developed by Johnson & Soderblom (1987).
HN Lib does not appear to belong to any known young stellar moving group and has kinematics typical of a stars in the galactic thin disk, indicating a likely age of 0.8-8.0 Gyr, which is in agreement with the very weak stellar activity. The proximity of HN Lib has made our target to be the subject of numerous searches for close companions (e.g., Jameson et al. 1983;Skrutskie et al. 1989;Nakajima et al. 1994;Simons et al. 1996;Oppenheimer et al. 2001;Hinz et al. 2002;Tanner et al. 2010;Jódar et al. 2013;Ward-Duong et al. 2015;Davison et al. 2015;Cortés-Contreras et al. 2017). The high-resolution images of Dieterich et al. (2012) were the most sensitive to very-low-mass stars, brown dwarfs, or planets using NICMOS on the Hubble with the F180M near-infrared filter. Dieterich et al. (2012) established companion magnitude and angular separation limits ranging from 11.5 mag at 0.4 arsec (2.5 au) to 18 mag at 3.0-4.0 arcsec (19-25 au) for companions that could be ruled out for HN Lib. These limits, together with the age of the system and the evolutionary models for very-low-mass stars and brown dwarfs with dusty atmospheres from Baraffe et al. (2015), discard any hypothetical companion with mass from 0.07 to 0.08 M ⊙ at projected orbital separations larger than 0.4 arsec, and from 0.025 to 0.06 M ⊙ at separations larger than 3.0-4.0 arcsec. Finally, as in (Caballero et al. 2022), we searched for objects with common Gaia DR3 parallax and proper motions up to a projected physical separation of 100 000 au. We found no hints of any wide potential companion.

Modeling techniques
To get a quick overview of all the different data available to us, we used the graphical user interface capabilities of  For the fitting process, juliet offers different options. We chose to use dynesty (Speagle 2020), which implements a nested sampling algorithm and supports computing Bayesian log evidence, ln Z, for the models.

Activity indicators
To identify signals that originate from stellar activity, we investigated the GLS periodograms for the range of activity indicators determined from the CARMENES VIS and NIR spectra as described in Sect. 2.1. Since our observations span several years, many of the periodograms suffer from strong signals with periods of 360 d to 380 d, as well as approximately half a year, caused by uncorrected contamination by telluric lines. For those spectroscopic activity indicators that show peaks at these periodicities, we proceeded with periodograms that were prewhitened by subtracting the best matching solution to these signals from the data. For clarity, we show only the periodograms for the activity indicators with peaks of less than 10 % false alarm probability (FAP) after the pre-whitening in Fig. 2. All shown indicators share a pair of signals with periods of ∼125 d and ∼94 d that are caused by aliasing due to the seasonal observability of HN Lib 5 . Both periods are close to the estimate of P rot,est = 102 ± 11 d based on the log R ′ HK measurement by Astudillo-Defru et al. (2017), which suggests that the underlying signal is indeed related to the stellar rotation. However, since both signals are within the tolerance of the predicted rotation period, we could not resolve the aliasing with the spectroscopic data.

Photometry
The GLS periodogram of the TJO photometric time series in Fig. 2 shows a peak at 95 d, with aliases at 74 d and 135 d due 3 https://github.com/3fon3fonov/exostriker 4 https://juliet.readthedocs.io/en/latest/index.html 5 The seasonal observability imprints a signal of about 1/365 d −1 to the window function of the data. Aliasing occurs at f a = f t ±m× f s , where f a are the alias frequencies, f t is the true underlying signal, m is an integer value and f s is the sampling frequency. Given a signal of ∼125 d, we therefore expect first-order aliasing at f a = 1/125 d + 1/365 d ≈ 1/93 d.  Signal search with GLS periodograms in the spectral activity indicators (left) and photometric data (right) of HN Lib. The data of the activity indicators were pre-whitened by subtracting the long-period ∼360 d or ∼180 d signals accordingly, whereas the photometric data were detrended using CAFM to remove the appearing long-term trends. Original periodograms are shown in grey and the corrected ones in black. We mark the 36-day signal from the RVs with a red solid line and the 113-day period with a blue solid line. The common alias pair of the activity indicators with periods of ∼125 d and ∼94 d are depicted by green dashed and solid lines, respectively. FAPs of 10, 1, and 0.1 % were calculated with the analytical expression of  and are shown with the horizontal grey dashed lines from bottom to top.
to the seasonal observations. Furthermore, a long periodic signal with a period of 263 d is visible, that could be an alias of a longer activity cycle of the star. Because of the data gap of about ∼200 d (and the consequently uncertain offset) between the two OSN observing runs and the short baseline of each data chunk, the GLS periodogram of the OSN data shows a large plateau for long periods. However, when inspected visually, the OSN and TJO photometric times series pretty much overlap (see Fig. A.4 for comparison). A periodicity of ≳ 100 d is also apparent in the periodogram of the daily binned data from the SuperWASP survey; however, despite being the most prominent peak, it is not significant in the GLS. A signal of an unknown origin at a period of 247 d is dominant in the MEarth data. In the regime of around 100 d, the strongest peak appears at 96 d, which is consistent with the TJO data.
Based on the log R ′ HK measurements and the GLS analysis of the activity indicators, we expected the rotation period to be in the range between 90 and 130 days. Consequently, we detrended the data before further analysis due to the appearing long-term trends. For this, we used wotan (Hippke et al. 2019) and Cosine Filtering with Autocorrelation Minimization (CFAM - Kip-ping et al. 2013;Rodenbeck et al. 2018), where we set the window length to 200 d. The CFAM had a particularly strong effect on the OSN data, as it removed the offset between the two observing seasons. Overall, the detrended photometric data share a common periodicity, which appears as a very significant signal in the combined GLS periodogram at a period of 97 d, with yearly aliases of first and second order on both sides.
To measure the photometric rotation period of HN Lib considering all photometric data combined, we used GP modeling. In contrast to the GLS periodograms, which are based on static sinusoidal models, GPs are able to represent also the quasiperiodic nature of stellar activity. For our purposes, we used two different kernels. The first one was a kernel that is the sum of two simple harmonic oscillators, in the following abbreviated as "dSHO" kernel (e.g. David et al. 2019;Gillen et al. 2020). A description of the implementation and parametrization of this kernel in juliet can be found in Kossakowski et al. (2021). The second kernel was the multiplication of an exponentialsine-squared kernel with a squared-exponential kernel, which is usually called the "quasi-periodic" (QP) kernel (Haywood et al.   2014; Rajpaul et al. 2015). The parametrization in juliet for this kernel was described by Espinoza et al. (2019).
Because they are far away in time from the MEarth, OSN, and TJO observations and, thus, also from our RVs, we did not include the SuperWASP data in the fits, considering that the star could have been in a different stage of the activity cycle (Baliunas et al. 1995;Díez Alonso et al. 2019;Fuhrmeister et al. 2023). We used shared hyperparameters between the instruments for the rotation period and the parameters that describe the coherence and shape of the periodic components of the kernels (i.e., the quality factor and the difference in the quality factor of the dSHO kernel and the α and Γ parameters for the quasi-periodic kernel). The amplitude of the GP was considered individually for each instrument to account for potential wavelength dependencies and different states of activity. A summary of the GP hyperparameters and the used priors can be found in Table B.1. The MEarth data used comparison stars that are not M dwarfs (Irwin et al. 2011;Newton et al. 2016), introducing systematics. Their light curves were corrected by applying a simultaneous linear detrending of these data using the common mode parameter, which is based on the observations of all M-dwarf targets in the MEarth data combined to make a lower cadence (binned) comparison light curve (Newton et al. 2016(Newton et al. , 2018. Our fit with uninformative priors on the period between 10 d and 200 d resulted in a period of 94 ± 4 d using the QP-GP and 96±2 d using the dSHO-GP. The photometry and the model from the dSHO-GP fit are shown in Fig. 3. The period from the photometric fit hence matches the ∼94 d period of the spectroscopic activity indicators (see the left panels of Fig. 2), although it is the least significant of the two occurring aliases in most of the nonpre-whitened periodograms. One explanation for it could be that both ∼125 d and ∼94 d are close to harmonics of the long-term signals with periods of approximately one year, which we removed by pre-whitening the data. Since ∼125 d matches a lower order harmonic, it is not unexpected that it is also more affected by residual signal in the data and thus increased in power in the GLS. In agreement with the log R ′ HK measurements and the spectral activity indicators, we therefore established a rotation period of 96 ± 2 d for HN Lib, the more precise value obtained from the dSHO-GP fit. In Table 3 we summarize the different values derived for the rotation period of the star using various methods (i.e., activity-rotation relationships, CARMENES spectroscopic activity indices, and photometric light curves) to ensure that all of them are consistent within their error bars. The error bars for the photometric cases come from our GP fit.

RV signal detection
We began the analysis of the RV data with a signal search using the GLS periodograms, as shown in Fig. 4. We considered a signal to be significant in the GLS periodogram if it has an FAP < 1 %. The CARMENES data (second row of Fig. 4) show two significant peaks at periods of ∼36.1 d and ∼113.1 d. When combining the CARMENES data with the HIRES and HARPS measurements (third row and following), both signals remain significant. However, multiple alias peaks appear on both sides of the 36-day signal due to the long time baseline and seasonal sampling of the combined data. Considering that this aliasing does not occur in the CARMENES data alone, we were consequently confident in the ∼36.1 d period of the 36-day signal. Further, a signal with a FAP of almost 0.1% is visible near the period of the stellar rotation period (third row of Fig. 4), as determined in the previous section. The peak is distinct from the 113-day signal and not related by aliasing, as shown by the dashed lines in Fig. 4. Subtracting the 36-day signal with a sinusoidal model considerably increases the power of the 113-day signal in the residuals (fourth row of Fig. 4), but decreases the power of the signal at the rotation period. The residuals of a simultaneous sinusoidal fit of the 36-day signal and the 113-day period do not show any further significant signals (fifth row of Fig. 4). Since the period of ∼113.1 d is close to the rotation period determined in Sect. 4.2, we assumed that it could be an imprint of the stellar activity onto the RVs and, consequently, performed a more in-depth analysis of this signal, which is presented in Sect. 4.5.
We verified that both significant RV signals are stable and coherent over the entire observational time baseline by producing the stacked Bayesian generalized Lomb-Scargle periodograms (s-BGLS, Mortier et al. 2015) shown in Fig. 5. The probability of the 36-day signal increases with time until a constant high degree is reached above one hundred observations. The aliases of 1 yr are also seen on both sides of the 36.1 d period, but with varying probability. With the increasing number of observations, the uncertainty of the period of the signal also becomes narrower, which is a behavior expected for signals with a Keplerian origin (Stock et al. 2020;Bluhm et al. 2021;Chaturvedi et al. 2022). Remarkably, the 113-day signal, for which we expected a more erratic behavior in case of a stellar activity origin, seems to be also a stable signal after reaching ∼ 80 observations. The region between 94 d and 96 d has an unstable behaviour revealed by changes in its width and probability when increasing the number of observations. As a result, it is likely to the site of stellar activity.

RV modeling and planetary parameters
We found two potentially planetary signals in the periodogram analysis: an isolated signal with a period of 36.1 d and a second signal close to the stellar rotation with a period of 113.1 d. Following this, we conducted a model comparison based on the Bayesian evidence to find the model that explains the RV data the best. Our base assumption was that the RV data can be fully explained by stellar activity (= "0P-model"). For this model, we used GP modeling as in Sect. 4.2 to account for the quasi-periodic nature of such signals. To ensure consistency with the photometry, we applied the dSHO kernel in our analysis. Therefore, we used individual GP amplitudes for each instrument to consider the different scatter in the HIRES, HARPS, and CARMENES data. Furthermore, we made use of the determined stellar rotation period and adopted a normal prior centered on our most precise measurement from the dSHO kernel of P rot ≈ 96 d, with three times the uncertainty of the photometric fit (±6 d). We compared this activity-only model with two more complex models that assume the presence of further signals in the data. The first one additionally includes the 36-day signal using a Keplerian component with a uniform period prior between 30 d and 40 d (hereafter, the "1P-model"). We set t 0 , the time of inferior conjunction, between BJD = 2457433 and 2457450 to avoid a split posterior due to the aliasing. Further, we re-parametrized the eccentricity and the argument of periastron by S 1 = √ e sin ω and S 2 = √ e cos ω to allow for a uniform sampling between −1 and 1 (Eastman et al. 2013). In another fit, we additionally considered the 113-day signal by adding a second Keplerian component to the model. For this, we used a period prior between 100 d and 120 d and t 0 between BJD = 2457433 and 2457550 (hereafter, the "2P-model"). For both 1P and 2P models, we tested circular as well as eccentric models for the Keplerian signals. For all instruments, we set individual uninformed uniform priors for the RV offset and the jitter. A summary of the used priors can be found in Table B.2.
We found that the models including the 36-day signal are generally significantly (i.e., ∆ ln Z > 5; Trotta 2008) better than the model considering only the stellar activity (see Table 4). Likewise, the models including two Keplerian signals are mostly favored against the 1P-models, except for the model in which both planetary orbits are assumed to be eccentric. The highest evidence corresponds to the two planetary model, one eccentric orbit at 36 d and another circular orbit at 113 d, plus the dSHO kernel of P rot at 96 d (2P (36 d-ecc, 113 d-circ) + dSHO-GP (96 d) ), which we chose as our final model for the further analysis (shown in Figs. 1 and 6). The inner planet orbits in a nearly circular orbit of eccentricity 0.079 0.090 0.055 . The posteriors from the fit for the planetary parameters and the GP can be found in Table 5 and Fig. A.1, and the instrumental posteriors in Table B.3. However, regardless of which model is chosen, all determined period and semiamplitude combinations of the 36-day signal, which is the focus of our analysis, agree within their uncertainties (see Fig. A.2).

Investigation of the 113-day signal
Our findings about the 113-day signal up to this point can be summarized as follows: We found a signal in the RVs that is close to, but not exactly at, the period of the stellar rotation period, which we determined consistently from photometry and spectroscopic activity indicators. Furthermore, the stellar rotation period of 96 ± 2 d and the period in question at 113 d are not related by aliasing. In fact, we see evidence for a periodic signal in the RVs that coincides with the stellar rotation period that is distinct in period from the yearly aliases of the 113-day signal. Our analysis of the signal's sBGLS periodogram and the model comparison show that it is very well described by a static Keplerian model.

Wavelength dependence
In the presence of strong stellar activity, we would expect a wavelength dependence of the RVs. A good first indicator for activity is the CRX, which expresses the RV-log λ-correlation Tal-Or et al. 2018;Jeffers et al. 2022). We used the CRX calculated from the CARMENES visual channel data to investigate the presence of stellar activity in HN Lib. However, the GLS periodogram of the CRX data does not show any peaks near the period in question (see Fig. A.3). Further,  both in the original data and after removing the 36-day signal, there is no significant correlation between the RVs and the CRX. Another method to investigate a potential wavelength dependence of a RV signal is to make use of the RVs computed for each échelle order (Bauer et al. 2020). To do this, we divided the CARMENES RVs into three roughly equal wavelength ranges and performed a fit to each of them to determine the semiamplitude of the 113-day signal. For the fit, we took the results from the best fit in Sect. 4.4 as normal priors for the 36-day signal and the dSHO-GP kernel, as well as the period of the 113-day signal. The prior of the semi-amplitude of the 113-day signal was set uniform between 0 m s −1 and 10 m s −1 . The results for the posterior distribution from those fits are depicted in Fig. 7. There are no significant discrepancies. The increase in amplitude and, at the same time, improvement in precision with increasing wavelength can be explained by the lower RV information con-

HN Lib b
The monitoring of the M dwarf HN Lib during our observing campaign within the CARMENES GTO program yielded a sub-Neptunian-mass planet discovery, namely HN Lib b. The minimum mass of the planet is M b sin i = 5.46 ± 0.75 M ⊕ , and it orbits its host star in a nearly circular orbit at a separation of a b = 0.1417 ± 0.0023 au with an orbital period of P b = 36.116 ± 0.029 d. HN Lib b's signal is unrelated to the stellar rotation period that we determined to be P rot = 96 ± 2 d and has no significant counterparts in the activity indicators. Furthermore, the stability and the coherence of the 36-d signal indicate a long-lived behavior, stable over the time span of our observations and thus providing strong arguments in favor of the planetary nature of HN Lib b. As it can be seen in Fig. 8, the planet adds to the sam- ple of small, sub-Neptunian mass planets, which make up the bulk of planets detected to orbit M dwarfs on short orbits and that are predicted to be the most abundant planet type around M dwarfs (Dressing & Charbonneau 2015;Sabotta et al. 2021;Pinamonti et al. 2022). Kopparapu et al. (2013Kopparapu et al. ( , 2014 calculated conservative estimations of the inner HZ around stars with stellar effective temperatures in the range 2600-7200 K for planetary masses between 0.1 and 5 M ⊕ . Taking into account the runaway greenhouse effect for 5 M ⊕ planets such as HN Lib b, the inner edge of HN Lib's HZ is accordingly placed at 0.101 au. With an instellation of S ≈ 0.50 S ⊕ and a semi-major axis of 0.142 au, HN Lib b thus resides in the conservative HZ of HN Lib. At a distance of only 6.25 pc, HN Lib is one of the closest systems to Earth with a planet in the conservative HZ. In the two extreme cases of the Bond albedo, A = 0.0 and A = 0.65, the T eq for HN Lib b falls in the range of 243-174 K for a non-reflecting planet and for a highly reflecting planet, respectively. In Fig. 9 we compare the position of HN Lib b with other known habitable zone sub-Neptune-type and super-Earth planets with (minimum) mass measurements around M dwarfs. However, as for all planets in this regime, the actual habitability depends on many factors. The most decisive is probably the nature of its atmosphere (Kaltenegger 2017 observations planned for future TESS sectors. Given the ambiguity of the mass-radius relation for small planets (e.g., Fulton et al. 2017;Van Eylen et al. 2018; Petigura 2020; Cloutier & Menou 2020), we thus only can hypothesize about its size and composition. Following the recent empirical analysis of small, transiting planets with mass measurements by Luque & Pallé (2022), there are three different plausible compositions given the mass of HN Lib b: a rocky planet with an Earth-like density, a water-rich planet, and a gaseous planet with an extended hydrogen atmosphere. Applying the theoretical mass-radius relationships by Zeng et al. (2019) to the different possibilities results in radii of ∼ 1.6 R ⊕ for an Earth-like planet, ∼ 2.0 R ⊕ for a water-rich planet, and ∼ 2.4 R ⊕ for 1% of hydrogen added to an Earth-like core (assuming an equilibrium temperature of 300 K in all three cases). The resulting Earth similarity index (Schulze-Makuch et al. 2011) for a rocky planet would be 0.60, comparable to TRAPPIST-1 g (Gillon et al. 2017;Grimm et al. 2018;Ducrot et al. 2020;Agol et al. 2021) or the recently discovered GJ 1002 b (Suárez Mascareño et al. 2023). For a full transit, the depth would be in the range of ∼1600-5400 ppm taking into account the different possible compositions. However, the probability that our planet transits is low (0.98 +1.02 −0.94 %). We explored the high energy environment of HN Lib b to test for the habitability conditions on the planet. There are no positive X-ray detections of HN Lib in the literature (Wood et al. 1994;Schmitt et al. 1995), but Stelzer et al. (2013) calculated an upper limit of log L X < 26.80. As an alternative approach, we used the rotation period of 96 ± 2 d to calculate the expected A&A proofs: manuscript no. HNLib X-ray emission, based on the relations of Wright et al. (2011), resulting in log L X ∼ 26.9, in approximate agreement with the aforementioned upper limit. We then used the relationships in Sanz-Forcada et al. (2011) to calculate the extreme ultraviolet (EUV) luminosity (L EUV ) corresponding to this X-ray flux, log L EUV = 27.9, and the expected (upper limit) mass loss rate in an energy-limited approach, 1.8 × 10 9 g s −1 or 0.0097 M ⊕ Gyr −1 . This mass loss rate is too low to pose any problem for the present stability of the planet's atmosphere, which is consistent with the planet's location in the potentially HZ of HN Lib. The stellar age based on the X-ray luminosity is quite uncertain (∼7 Gyr, Sanz-Forcada et al. 2011), but it is not consistent with values as low as 0.8 Gyr, which would otherwise imply X-ray activity approaching saturation, log L X ∼ 28.4, and a shorter rotation period.

HN Lib [c]
The RV observations show another significant signal at a period of ∼113 d, close to the stellar rotation period of 96 ± 2 d. Our model comparison shows that the signal is reasonably well described by a circular Keplerian orbit and, following the sBGLS analysis, seems to be reasonably stable over the whole period of our observations, spanning 21 years. When analyzing its wavelength dependency, we found no significant changes in the amplitude of the signal that would hint at its origins coming from stellar activity. However, given the proximity to the stellar rotation period, we could not fully rule out this scenario (especially due to the presence of differential rotation) and thus report it as a planet candidate. HN Lib [c] has a minimum mass of 9.7±1.9 M ⊕ and orbits HN Lib at a separation of 0.3040 ± 0.0051 au. In contrast to HN Lib b, the orbit lies outside the HZ of HN Lib. As for planet b, predictions about its composition are ambiguous because it resides in the same mass regime. However, due to the orbit beyond the ice line of HN Lib, an ice-rich planet is very likely (Burn et al. 2021). In combination with the also comparatively massive inner companion, the planetary system of HN Lib is therefore very interesting with regard to the formation and evolution of systems with planets in the HZ.

Summary
The newly discovered planet HN Lib b orbits a nearby M4.0 V star and has a minimum mass of M b sin i = 5.46±0.75 M ⊕ . Given the ambiguity of the radius in this mass range, this planet could be either rocky, a water-rich, or gaseous. The orbital period of the planet, P b = 36.116 ± 0.029 d, corresponds to a separation to its host of a b = 0.1417 ± 0.0023 au centered in the conservative HZ of HN Lib (S ≈ 0.50 S ⊕ ). HN Lib b is one of the closest planets orbiting its host in the HZ as seen from Earth and, assuming that it has a rocky composition, its Earth similarity index would be 0.60.
Additionally, our RV measurements show another significant signal with a period of P [c] = 113.46 ± 0.20 d, which is close to the stellar rotation period of P rot = 96 ± 2 d that we determined from photometric measurements -and which was additionally validated with spectroscopic activity indicators. Since we can neither confirm nor exclude that the signal is an imprint of stellar activity, although we do know is stable over the time span of our observations, we report it to be a planet candidate with a minimum mass of 9.7 ± 1.9 M ⊕ .    Tables   Table B.