The CARMENES search for exoplanets around M dwarfs. A long-period planet around GJ 1151 measured with CARMENES and HARPS-N data

Detecting a planetary companion in a short-period orbit through radio emission from the interaction with its host star is a new prospect in exoplanet science. Recently, a tantalising signal was found close to the low-mass stellar system GJ 1151 using LOFAR observations. We studied spectroscopic time-series data of GJ 1151 in order to search for planetary companions, investigate possible signatures of stellar magnetic activity, and to find possible explanations for the radio signal. We used the combined radial velocities measured from spectra acquired with the CARMENES, HARPS-N, and HPF instruments, extracted activity indices from those spectra in order to mitigate the impact of stellar magnetic activity on the data, and performed a detailed analysis of Gaia astrometry and all available photometric time series coming from the MEarth and ASAS-SN surveys. We found a M$>$10.6 M$_{\oplus}$ companion to GJ 1151 in a 390d orbit at a separation of 0.57 au. Evidence for a second modulation is also present; this could be due to long-term magnetic variability or a second (substellar) companion. The star shows episodes of elevated magnetic activity, one of which could be linked to the observed LOFAR radio emission. We show that it is highly unlikely that the detected GJ 1151 b, or any additional outer companion is the source of the detected signal. We cannot firmly rule out the suggested explanation of an undetected short-period planet that could be related to the radio emission, as we establish an upper limit of 1.2 M$_{\oplus}$ for the minimum mass.


Spectroscopic data
The CARMENES instrument consists of a visible (VIS) and a NIR spectrograph with wavelength coverages of 0.52-0.96µm and 0.96-1.71µm, and resolutions R of 94 600 and 80 400, respectively (Quirrenbach et al. 2014).The VIS data used in this study consist of the 70 epochs already published in Perger et al. (2021b) along with 31 additional observations.It is distributed over seven seasons, including 7 epochs from February 2016 to June 2016, 1 lone data point from January 2018, and 93 further epochs from February 2020 to April 2022, with a small gap from August to November.The basic reduction of the échelle spectra was done with caracal (CARMENES Reduction And CALibration software, Caballero et al. 2016).We further corrected the CARMENES VIS spectra for telluric lines using the method explained by Nagel et al. (2022).Radial velocities were then extracted with the serval code considering nightly zero points and instrumental drifts (see Sect. 2.2 for details).This dataset has the largest baseline of observations, relatively small uncertainties, and the lowest RV root mean square (rms).Two recent data points could not be considered in the calculation of the RVs because they showed signal-to-noise ratio (S /N) < 3.One further data point is an obvious RV outlier and another shows uncertainties larger than 5 σ.In the following, 97 CARMENES VIS epochs were considered.We list the measurements for each epoch in Table C.2.  (Gliese & Jahreiß 1979), G 122-49 (Giclas et al. 1971), or Karmn J11509+483 (Caballero et al. 2013).

Parameter
Value Ref.
HARPS-N is a spectrograph observing at visible wavelengths, which is installed at the Nasmyth B focus of the 3.58 m Telescopio Nazionale Galileo (TNG) at the Roque de Los Muchachos Observatory in La Palma, Spain (Cosentino et al. 2012).HARPS-N has a wavelength coverage of between 0.38 and 0.69 µm and a resolving power of R ∼115 000.We use 21 publicly available data points already used in the previous studies along with 25 newly observed epochs from the A44TAC_15 program.The dataset then consists of 46 epochs distributed over four seasons in two parts of approximately 3 months, from December 2018 to February 2019 and from February 2022 to April 2022.We extracted the RVs with both terra and serval pipelines and show the basic properties in Table 2 for completeness, but used the reduction by terra, because it shows the lower rms and uncertainty.This comparison is not done for the CARMENES spectra because serval includes a variety of corrections especially designed for the spectra of the instrument.As we can observe in Fig. 1, the two different parts show a large offset of around 15 m s −1 , resulting in the relatively large rms.We list the measurements for each epoch in Table C.1.Notes.Number of observations N, RV root mean square (rms) of the observed data, average internal RV uncertainties ⟨σ RV ⟩, time baseline T ; and average time interval between observations ∆t for the RV data for the different instruments and the combined dataset.
Fig. 1.RV time series (left panels) and periodograms for the individual RV datasets coming from CARMENES, HPF, and HARPS-N (from top to bottom).We show the behaviour at longer periods >50 days (middle panels) and down to 1 day (right panels).We mark the rotation period of the star at 140±10 d in blue, the orbital period of the proposed planet in green, and the prominent peaks in grey.
HPF is an NIR spectrograph installed at the 10 m Hobby-Eberly Telescope at the McDonald Observatory in Texas, USA (Mahadevan et al. 2012).It has a wavelength coverage from 0.8 to 1.27 µm and a resolution of R ∼55 000.Mahadevan et al. (priv. comm.)kindly provided us with the data used in their previous study of the star (Mahadevan et al. 2021).These data consist of 25 epochs acquired from March 2019 to June 2020, computed from the reduced 1D spectra using serval.We give an overview of the basic properties of the data in Table 2 and illustrate them in Fig. 1.As we can see, the HPF data have the shortest time baseline T , the shortest average sampling between observations ∆t, and the largest average uncertainty ⟨σ RV ⟩.
We further analyse 60 spectropolarimetric observations from ESPaDOnS 1 (Echelle SpectroPolarimetric Device for the Observation of Stars;Donati 2003).This is a high-resolution échelle spectrograph operating in the wavelength range from 0.37 to 1.05 µm with a resolution of R = 68 000, which is mounted at the CFHT (Canada-France-Hawaii-Telescope).GJ 1151 was observed over the course of six different nights and on ten different occasions from January to March 2017 (four Stokes I observations for one Stokes V measurement -the polarisation filter is in different angles).The dates are contemporaneous    Notes.MEarth combined data are shown under the "combined" column.⟨σ phot ⟩ refers to the photometric uncertainty.See notes in Table 2 for more details.
with the RV monitoring observations.We use the ESPaDOnS spectra to access the longitudinal stellar magnetic field and to extract log R ′ HK values (see Table 1) following the procedure in Perdelwitz et al. (2021).The RVs extracted from the spectra are not used because they show large uncertainties of approximately 15 m s −1 .

Data extraction from spectroscopic observations
There are two main methods to extract RVs from high-resolution échelle spectra.The binary-mask technique was developed early on and implemented in the data-reduction system (DRS) of both HARPS and HARPS-N.This method computes the Doppler shift via the cross-correlation function (CCF) between the stellar spectrum and that of a predefined binary mask.The resulting shape of the CCF is measured through its moments, namely its width and height (FWHM -full width half maximum, and CON -contrast), and its asymmetry, the bisector inverse span (BIS).Those indices have been shown to be influenced by the effects of stellar activity on the spectra.The CCF technique (Simkin 1974;Rajpaul et al. 2021;Klein et al. 2022) is also implemented by the raccoon (Radial velocities and Activity indicators from Cross-COrrelatiON with masks, Lafarga et al. 2020) code, which we also use here.The raccoon pipeline for GJ 1151 uses a mask of spectral type (M4.5 V) derived from the CARMENES observations.
A second method, dubbed template matching, has proven to yield better results than the CCFs (Anglada-Escudé & Butler 2012), especially for low-mass stars, where many of the abundant spectral lines are blended.In such cases, a template spectrum of high S/N is constructed by combining all the observed spectra of a target.Each observed spectrum is then optimally matched to this template.This method constitutes the basis of the terra (Template-Enhanced Radial velocity Re-analysis Application, Anglada-Escudé & Butler 2012) and serval (SpEctrum Radial Velocity AnaLyser, Zechmeister et al. 2018) codes.
The serval code includes the calculation of two activity indicators, the differential line width (dLW), measuring the width of the average line, similar to CCF FWHM, and the chromatic index (CRX), which is the slope over the logarithmic wavelength of the RVs calculated for every order of the échelle spectrum (Baroch et al. 2020).Further, the CARMENES processing pipeline delivers measurements for a variety of line features sensitive to stellar activity by comparing their fluxes with the continuum values.Those include the Hα line at 6563 Å, which can trace the ionised hydrogen content, and the Ca II infrared triplet (CaIRT).For this index, the three ionised calcium emission lines of wavelengths 8498, 8542, and 8662 Å are measured.We then create the CaIRT index as the average over the three individual measurements.We further calculate the sodium line indices Na I D2 and Na I D1 at vacuum wavelengths of 5890 and 5896 Å, and use the average to create the Na I D index for the line doublet.
We list the CCF activity indices calculated by the raccoon code, the described indices calculated by serval, and the Hα index in Table C.2 for the CARMENES spectra and in Table C.1 for HARPS-N.Additionally, we show the calcium and sodium indices calculated from the CARMENES spectra.Five HARPS-N data points including publicly available spectra showed problems in the calculation of the CCF due to low S/N.

Photometric data
MEarth-North is an array of eight 40 cm telescopes located at the Fred Lawrence Whipple Observatory (FLWO) on Mount Hopkins, Arizona, USA (Charbonneau et al. 2008), performing photometric monitoring of the closest late M dwarfs since 2008.Our target has been observed with telescope 2, producing two different datasets from 2008 to 2010 and from 2011 to 2020, and with telescope 7 from 2011 to 2020.We provide an overview of the basic properties of the datasets in Table 3.We also show the combined data in Fig. 2  ASAS-SN (All-Sky Automated Survey for SuperNovae) is a survey instrument composed of five stations, two of which are located at Cerro Tololo International Observatory (CTIO, Chile), and one each at Haleakala Observatory (Hawaii), McDonald Observatory (Texas), and South African Astrophysical Observatory (SAAO, Sutherland, South Africa; Kochanek et al. 2017).As its automated data pipeline is not optimised for the observation of high-proper-motion stars, we calculated new coordinates prior to the download of the light curves for each year of observations.We calculated those coordinates at the middle of every season of observations (April) following the procedure by Trifonov et al. (2021).We analysed observations in the V band from July 2013 to January 2019 and in the g band from September 2017 to June 2022.The latter set includes a linear trend (43.5 ± 4.8) × 10 −6 mag d −1 , which we subtracted in order to improve the sampling of the period range we are interested in.We show the nightly binned and 3σ-clipped data in Fig. A.1 and their basic properties in Table 3.
We also inspected the TESS mission photometry of GJ 1151, which was observed from 23 February to 17 March 2020 (Sector 22) and from 30 January to 24 February 2022 (Sector 48).As shown already by Pope et al. (2021), no rotation A50, page 4 of 17 modulations (the rotation period is not covered by the 27day TESS sectors), planetary transits, or obvious flares can be found.
Furthermore, we collected photometry at 37 epochs (August 2019 to May 2022) from 320 observations from the Gaia DR3 database (Gaia Collaboration 2021) after performing a retrieval on its 'Variable dataset of Solar like stars' (stars with spectral type later than F5).We could not identify any significant signals in the noisy periodogram of the low-cadence Gaia data.

Periodogram analysis
To search for periodic signals in our data, we used the Generalised Lomb Scargle (GLS) periodogram as described in Zechmeister & Kürster (2009).The significance of the signals in the periodograms is computed with the false-alarm probability (FAP) using the bootstrap method (Efron & Tibshirani 1985).We consider the signals to be tentative or significant if their FAP is below the 1% or the 0.1% level, respectively.We perform a pre-whitening process, which consists of fitting the data with a sine curve with the same period and phase as the most significant peak in the periodogram.This signal is removed and a re-analysis is performed on the residuals until no further significant signals are found in the data.

Model comparison
As an alternative approach for analysing the data, we searched for a best-fitting model through Bayesian statistics (Trotta 2008), as implemented in juliet (Espinoza et al. 2019).Starting with the simple combination of datasets (null model, e.g. for the MEarth dataset combination), we calculate different models of increasing complexity and accept a more complex model over a simpler one as significant if the difference of the logarithm of the Bayesian evidence (or log-evidence) between the more complex and simpler model is ∆ ln Z > 5, as discussed in Trotta (2008).The Bayesian evidence is the integral of the maximum likelihood estimator over the prior volume, which is the parameter space over which we evaluate our models and that we choose to be as large as possible.We calculate ∆ ln Z with juliet using the dynamic nested sampler dynesty, which iteratively evaluates the integral of the maximum likelihood over different live points.Those with poorest likelihood are replaced by randomly chosen points that improve the Bayesian evidence, until the improvement is below a certain threshold.As a dynamic nested sampler, the number of live points that are used is also updated on the fly, with a starting number of 1000.This methodology allows us to achieve precise posterior distributions of every parameter of the model within the limits of the introduced priors thanks to the efficient exploration of the prior space.

Model definitions
We fit a variety of different models, M, to our data.In the following, we describe each of the components from which the models are constructed including their respective parameters and hyper-parameters, which we generally evaluate from wide uninformative priors.
-Keplerian and sinusoid curve: with τ = t − t 0,i , which is the observed time minus an initial time t 0,i of model i, with K i being its amplitude, P i the period, ω i the argument of periastron, e i the eccentricity, and ν i the true anomaly of model i.The periods are evaluated from 1 d up to two times the baseline of our observations, 2T .We apply a uniform prior on the log-scale of the observed times so that we can equally explore the different timescales of our problem.The amplitude is uniformly explored up to 2rms RV , the eccentricity is explored with a uniform prior from 0 to 1, and ω from 0 to 360 deg.A sinusoid is modelled in the case of e i = 0. Other tests with wider priors on P i and K i were explored in order to avoid missing very eccentric signals.These solutions at periods much larger than T were highly dependant on the offsets between datasets.
-Individual offsets, linear and quadratic terms: where γ j is the offset of the j-th dataset (e.g. for different instruments CV, HA, HP for CARMENES, HARPS, HPF), and .γ and ..
γ are the linear and quadratic coefficients of a parabola.We apply uniform priors on the offsets from −2rms RV to 2rms RV .
For the linear and quadratic term, we set uniform priors from −2rms RV T −1 to 2rms RV T −1 and −2rms RV T −2 to 2rms RV T −2 , respectively, where T is the time baseline of the observations.-Individual jitter: Every dataset is assigned an individual jitter term σ j representing the relative significance between different datasets j and accounting for possible underestimations of the individual uncertainties.The jitter is added quadratically to the individual errors σ RV of each epoch in the calculation of the model likelihoods.Similarly to the offsets, we apply uniform priors of up to 2rms RV .
-Gaussian processes: Our modelling of the stellar contribution to the datasets considers the usage of a specific kernel S in the application of a Gaussian process regression.Those kernels are able to reproduce the effect of inhomogeneities on the stellar surface, and have a quasi-periodic form, where the signal of the rotation is allowed to change its parameters, namely amplitude and phase, over a certain time interval.We applied the quasi-periodic plus cosine (QPC) kernel, as now included in juliet, with the following form (Perger et al. 2021a): where λ is the exponential decay representing the average lifetime of starspots, K = h 2 1 + h 2 2 is the average amplitude of the signal, w 0 = 0.31 is a constant, and P represents the rotational period of the star.We set a log-uniform prior on λ from 60 d to the time baseline T of our observations.For P, we use wide uninformative priors.

Photometry
We show in Figs. 2 and A (Newton et al. 2016).It can be seen in the additional MEarth sets that the signal actually appears at a significantly longer period (∼140 d), and that yearly aliases of this period are visible at around 220 and 100 d.Therefore, we update the rotational period to 140±10 d, setting a conservative error bar that includes the prominent signals observed in the periodograms of the ASAS-SN and MEarth data.We performed a combined analysis of the three different MEarth datasets where we searched for the offsets and jitters for each set that maximise the log-evidence (see Sect. 3.3).Besides the main signal at around 140 d, we also identify a modulation of 500 d.We attempted to fit a GP model with a QPC kernel to the photometry, but the process is not able to converge on a specific period or lifetime of the spots.Most importantly, we do not see any periodicity at ∼390 d in the photometric data, which is the period for the RV signal of the potential planetary companion.

Radial velocities
The RVs of the different datasets are shown in Fig. 1.We can see that the periodogram of the CARMENES data contains multiple peaks with large power at >200 d, including a long-term modulation of the order of the time baseline of the observations (>2000 d).The periodogram of the HPF data does not show any significant signal except for the tentative variability at around 2 d reported by Mahadevan et al. (2021), which was later found to be caused by the poorly sampled long-period signal (Perger et al. 2021b).The periodograms of the sparsely sampled HARPS-N datasets are dominated by the window function and do not show clear periodogram peaks.No RV dataset shows significant power at the rotation period of 140 d.
We fit different models with increasing level of complexity to the three RV datasets in order to optimally combine the data and to find significant signals in the full data set.Each of the parameters of the models is evaluated over wide uninformative priors.All the models should represent either a stellar or a planetary signal imprinted on the measurements.We show the log-evidence values of all models as evaluated by the juliet code in Table 4.
We obtain a value of ln Z = −437.5 as the largest logevidence.Models with ∆ ln Z < 5 with respect to this model, ln Z < −442.5, are statistically equivalent.All such models have both a common, strong signal at ∼390 days and a second longterm signal.Given our statistical threshold, the 390d signal is best fitted by a sinusoid.The second signal cannot be as clearly identified because the periods found from both sinusoid and Keplerian models are below the Nyquist frequency of f = 0.0087 d −1 (P = 0.5T ≈ 1145 d) and do not converge on a specific period.This is also true for the hyper-parameters of the fitted QPC kernel in the GP approach, which also do not converge and which we therefore do not show in the table.As a Notes.∆ ln Z refers to the difference in log-evidence of the best model compared to any other model.The first word in each line of the first column refers to the variation with the shortest period (around 390 d), the second word to the variation with the largest period.The model we chose to be most relevant for this study is shown in bold face.result, we prefer the less complex model and fit the second signal using a parabola.We show the combined time series and the corresponding periodograms in Fig. 3.

Detection of a sub-Neptune planet
Following the criteria of an improvement in ln Z > 5, the bestperforming model is the combination of a sinusoid with a parabola (Table 4).Due to the lack of similar periodicities in any of our stellar activity tracers (Sect.4.2), we conclude that the sinusoidal modulation at 390 days is most likely produced by a planetary companion in a nearly circular orbit.The strength of the signal, its clear visibility in the RV curve, and our thorough RV extraction make it furthermore very unlikely for the signal to be introduced by any yearly effects (as reported for HARPS data in Dumusque et al. 2015).We show all fitted parameters of this model in Table 5, including the detailed priors and additionally derived parameters.A corner plot of the posterior distribution for those parameters is shown in Fig. B.1.Following those parameters, GJ 1151 c appears to be a sub-Neptune mass planet with a minimum mass of 10.6 M ⊕ , and in an orbit of non-measurable eccentricity.The orbital distance is 0.57 au, which places the planet well outside the habitable zone around GJ 1151.We further show the RV residuals after the application of the preferred model and the corresponding periodograms in Fig. 3.The RV measurements are shown in Fig. 4 after subtracting the parabola and being phase-folded to the planetary orbital period.

Activity level of GJ 1151
GJ 1151 is a slow rotator with an inferred rotation period of P = 140 ± 10 d (Sect.4.1) and projected equatorial velocity v sin i < 2 km s −1 .Its X-ray luminosity is L X ∼ 5.5 × 10 26 erg s −1 (see Table 1 for properties and related references).We calculate a low value of log R ′ HK = −5.22 ± 0.18 dex using the ESPaDOnS observations from early 2017.We find published pEW(Hα) values of 0.0 ± 0.2 Å from the CAFE spectrograph (Calar Alto Fiber-fed Echelle; Jeffers et al. 2018)     Notes.Detection status refers to the Stokes V detection diagnosis (Donati et al. 1997), and the longitudinal magnetic field was calculated following Brown et al. (2022).
With the additional epochs of observations acquired with CARMENES, we update this value to an average of −0.08 ± 0.15 Å. GJ 1151 can be classified as inactive when adopting a threshold of −0.(Donati et al. 1997), the spectropolarimetric ESPaDOnS observations were able to marginally detect its magnetic field, with one epoch being on the verge of a definite detection (see Table 6).We compute the values for the longitudinal magnetic field, B l , using the method described in Brown et al. (2022).
On the other hand, the residual RVs (after the subtraction of our best model) still present a scatter of ∼3 m s −1 .This is rather high given the median uncertainty of the CARMENES data, 1.8 m s −1 , and the reported experiences when using the CARMENES survey to investigate these kinds of targets (Ribas et al. 2023).From mid-2021 (BJD ∼ 2 459 400) onwards, the star also seemed to enter into a more active phase as can be seen in the Hα indices as depicted in Fig. 2 for CARMENES and in Fig. A.3 for HARPS-N measurements, and by the increase of the RV residual scatter in the latest observations (see Fig. 3).We also find signals of long-term variability in some activity indices and photometry, which could be induced by such an increase in magnetic activity.

Constraints on the long-term RV signal
Besides the ∼390 d signal that we classify as a Keplerian motion induced by a planetary companion, our RV data additionally show a long-term modulation with a timescale comparable to the baseline of the observations.We choose to fit a parabola to this modulation as the simplest model providing a good match to the data (Table 4).However, this does not provide any information as to the nature of the observed RV variability.
We can assume that the curvature in the RVs is due to a longperiod companion.Following Kipping et al. (2011), the approximate minimum-velocity semi-amplitude K > 7.5 m s −1 and measured RV acceleration γ = −7.54+0.59  −0.61 × 10 −6 m s −1 d −2 can be combined to provide a lower limit on the orbital period P 2 > 17.2 ± 0.7 yr and minimum mass M 2 sin i > 0.16 ± 0.10 M Jup .We expand on this possibility using astrometric observations in the following section.
The possibility of a long-period companion was tested in our modelling.The posteriors of the fit with a Keplerian could be compatible with a sub-Saturn planet on a highly eccentric orbit, e ∼ 0.5 at 1.7 au.However, the resulting period, 2268 d, is comparable to the time baseline and well above the period calculated from the Nyquist frequency.Moreover, the eccentricity value seems to depend strongly on one single measurement at BJD 2 458 140 (Fig. 3), which casts doubt on the robustness of the solution.
Giving the large scatter of 3 m s −1 of the uncorrelated residual RVs, the possibility for this modulation being of stellar origin seems high.As already explained, we considered a model with a GP and the QPC kernel (Sect.4.3) to account for the possibility that stellar surface phenomena imprint the rotational period on the RVs, but no conclusive results were obtained.However, the variations of scatter of the RV residuals and of the Hα-index measurements, and the long-term periodicities found in some activity indicators and photometric time series, indicate a stellar influence on the data.However, this influence could only depend on episodic increases in the magnetic activity level of GJ 1151.As it oscillates between more or less active phases, the star would imprint a larger or smaller scatter on the RV data, respectively.

Constraints on companions from absolute astrometry
Additional constraints on an outer companion can be provided based on information from the 320 observations from the Gaia DR3 archive.The renormalised unit weight error (RUWE) statistic for GJ 1151 is 0.99, a value well below the threshold of 1.4 (e.g.Lindegren et al. 2018Lindegren et al. , 2021)).Therefore, a single-star model describes the data satisfactorily.We further followed the approach by e.g.Belokurov et al. (2020) and Penoyre et al. (2020) to investigate the range of orbital separations and companion masses that would induce excess astrometric residuals with respect to a single-star model, therefore producing larger RUWE values.
We generated synthetic Gaia observations of GJ 1151 using the nominal astrometric parameters, and the values of observation times, scan angle, and along-scan parallax factors encompassing the mission time-span using the Gaia observation forecast tool 2 .We linearly added orbital motion effects due to companions with orbital periods in the range from 15 to 40 yr and planetary mass in the range from 0.2 to 40 M Jup , taking 2 https://gaia.esac.esa.int/gost/A50, page 9 of 17 A&A 671, A50 (2023) Fig. 5. Gaia DR3 sensitivity to companions of given mass (in M Jup ) as a function of orbital period (in year) orbiting GJ 1151.Solid, dashed, and dashed-dotted lines correspond to iso-probability curves for 90%, 95%, and 99% probability of a companion of the given properties producing RUWE > 0.99.
into account the lower limits on both parameters from the RV datasets, and assuming a primary mass of 0.16 M ⊙ from Table 1.One hundred random realisations of the remaining orbital elements were produced for each mass-period pair, all drawn from uniform distributions within their nominal intervals.We then perturbed the observations of each system with Gaussian measurement uncertainties appropriate for the case of a G = 11.7 mag star such as GJ 1151 (see e.g.Fig. 3 in Holl et al. 2022).Finally, we performed linear fits to the synthetic data using a single-star model, and recorded the fraction of systems with RUWE exceeding 0.99 for each period-mass pair.Figure 5 shows iso-probability contours in companion mass-period space corresponding to the fractions of systems with RUWE > 0.99.Already for the minimum of a 17 yr orbit, the planetary mass of a companion would be of 3(8) M Jup in order to produce a RUWE in excess of 0.99 with 90% (99%) probability.
In summary, Gaia DR3 astrometry does not have sufficient sensitivity to probe for the presence of a giant outer companion of GJ 1151 with true masses close to the lower limit set by the constraints from RVs.

Detection limits for a short-period planet
The residuals from the best-fitting model produce a periodogram with no significant remaining signals (Fig. 3), although with a large scatter of 3 m s −1 .For this reason, a yet undetected terrestrial planet compatible with the constraints from Vedantham et al. (2020) could be responsible for the LOFAR radio emission.Assuming that the residuals are uncorrelated, we calculate the detection limits of our data following Bonfils et al. (2013) and Zechmeister et al. (2009).To test this possibility, we inject circular planetary orbits with periods from 1 to 5 d, with a linear sampling of (20 000 d) −1 for 12 different phases.We start our trials with the semi-amplitude of the best-fit sine wave at every frequency.We iteratively increase the semi-amplitude of the injected planet in steps of 0.01 m s −1 until it generates a significant signal in the periodogram of the residuals (FAP < 0.1%).The FAP is calculated using the approximate method by Baluev (2008), which is computationally less expensive than the bootstrap method used before (see Sect. 3.1).Therefore, we obtain the lower limit semi-amplitude for a significant signal to be detectable in our data given the noise level.Finally, we convert this value into an upper limit on the minimum mass of a planet to remain undetected.As we can see in Fig. 6, we determine a limiting semi-amplitude of K = 1.52 ± 0.25 m s −1 for the studied period range, which translates into upper limits on the minimum mass of a planet from 0.73 to 1.25 M ⊕ for 1 to 5d orbits, respectively.We thereby only confirm our previously published value (Perger et al. 2021b) because of the increase in the RV residual scatter in the last observational season.

LOFAR radio signal
The LOFAR radio emission detection by Vedantham et al. (2020) placed GJ 1151 in the spotlight.The signal was found by crossmatching nearby stars (d < 20 pc) from the Gaia database with the LOFAR Two-Metre Sky Survey (LoTSS, Shimwell et al. 2019) data release (a similar following study, Callingham et al. 2021, expanded such detections to a further 19 M-dwarfs).In particular, Vedantham et al. (2020) found a highly significant signal in one of the eight-hour pointings, during which the emission was detected for the full length of the observation.The emission was broadband, covering the entirety of the observed frequency range 120−167 MHz, with a nearly flat spectral shape.The emission was coherent and highly circularly polarised.All these features point to the electron cyclotron maser (ECM) mechanism, where a population of electrons group together in gyration phase and emit strongly.The resulting cut-off frequency, ν = 2.8 B[G] MHz, is a direct indication of the local magnetic field.The LoTSS upper frequency (167 MHz) then points to a local magnetic field of at least B ∼ 60 G.However, we note that despite the multiple observations of GJ 1151 with the LOFAR instrument, the source has not been detected ever since (Vedantham, priv. comm.).
With our current data and analysis, we can perform an assessment of the possible source of the low-frequency radio emission.We discuss four possible scenarios below.

Stellar emission.
Radio emission from stars is associated with strong activity and magnetic fields stronger than the upper limit of 680 G inferred by Reiners et al. (2022).It would be quite uncommon for the observed radio emission to be produced A50, page 10 of 17 by GJ 1151 itself.Also, the characteristics (duration, polarisation, bandwidth) of the observation do not fit into the standard emission processes seen on stars, whether of incoherent gyrosynchrotron origin, or from coherent plasma emission (Vedantham et al. 2020).Nevertheless, we cannot discard that the LOFAR detection is probing a yet unexplored domain of stellar radio emission at low frequencies, especially considering the apparent episodes of increased activity of the star.
Planetary emission.The observed radio signal could be produced by auroral activity on the newly discovered planet GJ 1151 c.However, planets with sub-Jovian masses are not expected to possess strong magnetic fields (e.g.Stevenson 2010, Table 1).The dependence of the emitted frequencies on the magnetic field, which further depends on the planetary mass, suggests that a magnetic field B > 60 G would need a giant planet with a mass of at least ∼700 M ⊕ (according to the scaling law by Grießmeier et al. 2007; other scaling laws provide similarly large values).This mass is more than twice the mass of Jupiter.Agreement with our orbital solution and RV semi-amplitude would require a very small, very unlikely orbital inclination angle, with a nearly face-on configuration (sin i < 0.016, i < 0.9 deg, probability ∼0.01%).Using the same approach as that described in Sect.5.4, we find that, at a separation of 0.57 au, GJ 1151 c would produce a larger RUWE than the one reported in the Gaia DR3 archive with 99% probability if its true mass is ≳1 M Jup , corresponding to i ≃ 2 deg.The upper limit on the true mass of GJ 1151 c from Gaia DR3 astrometry is not very informative, but it helps to narrowly rule out a true mass of 2 Jup for the companion, casting further doubt on the hypothesis of radio emission from the planet itself.
Planet-star interaction.GJ 1151 c could interact with its host star, giving rise to stellar auroral radio emission.We computed the radio emission arising from this star-planet interaction for both an open and a closed dipolar geometry, and for two different models of star-planet interaction: the Zarka-Lanza model (Zarka 2007;Lanza 2009) and the Saur-Turnpenney model (Saur et al. 2013;Turnpenney et al. 2018).For illustration purposes, we show the results for the open magnetic field topologue in Fig. 7.The bottom panel shows the predicted flux density as a function of orbital distance arising from a star-planet interaction involving GJ 1151 c for the case of an open magnetic field geometry of the star, following the prescriptions in Pérez-Torres et al. ( 2021), and assuming a 50 G magnetic field at the local site where the cyclotron emission is produced, and an unmagnetised planet.This magnetic field corresponds to an electron-cyclotron frequency of 140 MHz, which is the central observing frequency of the LOFAR observations.We note that the level of expected radio emission at the orbital position of GJ 1151 c is orders of magnitude below the typical detection sensitivity of LOFAR (∼0.1 mJy after a typical 4-8 h run).We further note that for any realistic magnetic field of the planet, the expected flux density level remains well below any detection limit.The top panel shows the Alfvén Mach number as a function of orbital distance and orbital period.The interaction between planet GJ 1151 c and its host star happens in the supra-Alfvénic regime, meaning that energy and momentum cannot be transferred to the star.That is, for orbital periods above approximately 100 days, all motion is in the supra-Alfvénic regime, and the flux densities drawn in the bottom panel do not apply, as the Poynting flux needed to feed the radio emission is never transferred to the star.Thus, starplanet interaction between GJ 1151 and GJ 1151 c cannot account for the observed LOFAR radio emission.For further details, we refer to Pérez-Torres et al. ( 2021).Alfv > 1; top panel) at the location of the planet (vertical dashed line).Therefore, star-planet interaction cannot account for the observed LOFAR radio emission.The planet in the bottom panel is drawn at around 0.1 mJy in the y-axis, which corresponds to the 1 σ level of the LOFAR radio observations.See main text for details.
Inner planet.A low-mass planet in a close orbit around its (weakly) magnetised host star (B = 60 G, or 30 G in the case of bright second harmonic emission) could induce the detected radio emission, as proposed initially by Vedantham et al. (2020).As calculated in the previous section, the planet would have to have a minimum mass of below 1.2 M ⊕ and would therefore remain undetected given the properties of our dataset.

Conclusions
The M4.5 dwarf GJ 1151 was monitored from 2016 to 2022 with the CARMENES, HARPS-N, and HPF instruments.The resulting datasets can be used to characterise both the planetary system of the star as well as the effects arising from stellar magnetic activity.
We report the discovery of a planet compatible with sub-Neptune mass, M sin i = 10.62 +1.31  −1.47 M ⊕ , in a circular orbit with a period of P = 389.7 +5.4  −6.5 d and a semi-major axis of a = 0.5714 +0.0053 −0.0064 au.We confirm the planetary nature of the radialvelocity signal by analysing spectroscopic activity diagnostics and also photometric variability of the star using archival data from MEarth and ASAS-SN, which have been collecting data on A50, page 11 of 17 A&A 671, A50 (2023) this target over the last 10 yr.We update the rotation period of the star to 140 ± 10 d, which is compatible with the upper limit of 2 m s −1 for v sin i.We also analyse a large number of stellar indicators derived from our spectroscopic observations with CARMENES and HARPS-N.We confirm that GJ 1151 is a star with a generally low magnetic activity level.However, the star shows episodes of larger magnetic activity where it induces variability in RVs, activity indices (e.g.Hα emission), and photometry, and shows a detectable magnetic field with ESPaDOnS observations.
We find that it is highly unlikely that GJ 1151 c, by itself or by interacting with its host star, is responsible for the detected LOFAR radio emission.Also, it seems unlikely that an outer companion could be responsible, namely either a sub-Saturn on an eccentric orbit (as suggested by the RVs) or a Jupiter-like planet with P > 17 yr (from the limits from Gaia observations).The Earth-mass planet in a close-in orbit proposed by Vedantham et al. (2020) remains a possible explanation for the radio signal, but this possibility has a restricted parameter space as it remains undetected in our data.Using the residuals from our best-fit model, we confirm the previous detection limits.Those translate into an upper limit for the minimum mass of a planet candidate at 0.73-1.25 M ⊕ for 1-5 d orbits, respectively.However, we note that the fact that LOFAR has not detected GJ 1151 again does not exclude the possibility that the observed radio emission could be due to a low-mass planetary companion.To shed light on these different scenarios, a simultaneous, multiwavelength radio campaign with different facilities and covering a time baseline of a few days (in order to follow a possible planet along its orbit) could confirm the previous LOFAR detection, and determine the frequency cut-off, which would give a direct measurement of the surface magnetic field.
Particularly surprising is the fact that the remaining RV residuals after the subtraction of both the sinusoidal and quadratic models are still of the order of 3 m s −1 .This appears to be quite high for a star with such a low activity level given the reported experience from the CARMENES survey (Ribas et al. 2023).The target will certainly continue to be observed in the future with the CARMENES instrument with the goal of better constraining the long-term modulations and the probably strong influence of the varying levels of stellar activity of GJ 1151.The apparently episodic detection of a coherent, polarised, and broad-band radio emission, although not predicted for a generally inactive star such as GJ 1151, could be the result of an unknown mechanism able to produce low-frequency emission in weakly magnetised stars, possibly compatible with at least some of the 19 signals detected in other Gaia-LOFAR matching searches (Callingham et al. 2021).
and the individual data in Fig. A.1, after the application of a daily binning and a 3σ clipping of the data.

Fig. 2 .
Fig. 2. Photometric and spectral index time series.Significant signals are marked.We show the full dataset (left panels) and their periodograms for longer periods (>50 d, middle panels) and for all periods (right panels).The grey vertical dashed lines mark the significant peaks.We show the period of the planet in green and the rotational period of the star at 140±10 d in blue.For the combined MEarth data, we show the different individual sets in colours.The red dashed horizontal lines show the FAP levels of 10, 1, and 0.1%.A50, page 5 of 17 .1 the data and periodograms for different individual and combined photometric datasets.In the ASAS-SN data, we identify long-term modulations at 490 and 1055 d, and prominent signals around 140 d and 210 d, which may be connected by yearly aliasing, that is (140 −1 −365 −1 ) −1 .The MEarth dataset from 2008 to 2010 is the basis for the A50, page 6 of 17 previously published rotational periods of GJ 1151 of 125 d (Díez Alonso et al. 2019) and 117.6 d

Figure 2
Figure 2 shows the spectroscopic indices that have significant signals calculated from CARMENES and HARPS-N spectra using serval, raccoon, and terra as discussed.All other indices are shown in Figs.A.2 and A.3.The distribution of peaks in the periodograms of the FWHM, Hα, and Na I D indices from CARMENES observations is similar to that found in the photometry.All of these indices show long-term variability at >1000 d, and possibly the two yearly aliases of the 140d rotation period at around 230 and 100 d.On the other hand, no index shows a clear signal at 390 d or at 140 d.

Fig. 3 .
Fig. 3. RV time-series data (left panels) and periodograms (red dashed lines as the 0.1, 1, and 10% FAP levels) of the best-fitting model consisting of a sinusoid and a parabola.The middle panels show the low-frequency portion of the periodograms (P > 50 d) and we mark the rotational period at 140, including its uncertainty of ±10 d, with the blue vertical area, and the planetary period at 390 d with the green vertical line.In the right panel, the full periodogram is shown.In the top panels, we show the data as observed (green: CARMENES, red: HPF, blue: HARPS-N) and indicate the best-fitting model by the black curve.The middle row shows the RV residuals after the subtraction of the parabola.The bottom row shows the residuals of the best fit and we provide the RV rms.Table 5. Parameters for the best-fit model of a circular orbit and a second signal represented by a second-order polynomial.Parameter Units Sinusoid + parabola Priors model Fitted parameters P 1 d 389.7 +5.4

Fig. 4 .
Fig. 4. Phase-folded RV data of the sub-Neptune planet orbit for the sinusoid + parabola model with the parabola term subtracted.Symbols are coloured as in Fig. 3. Black dots show the RV average of each 0.1 phase bin.

Fig. 6 .
Fig. 6.Detection limits of the RV residuals for the period of the assumed low-mass planet responsible for the LOFAR signal, after subtracting the best-fitting model consisting of the effect of GJ 1151 c and an additional quadratic polynomial.The black lines correspond to the minimum masses of the planetary signals as recovered by the simulations.The red line corresponds to the minimum mass of the average limiting semiamplitude along with its standard deviation limits in blue.

Fig. 7 .
Fig. 7. Expected flux density for auroral radio emission arising from a star-planet interaction in the system GJ 1151 -GJ 1151 c as a function of orbital distance.The interaction is clearly in the supra-Alfvénic regime (i.e.Alfvén Mach number M A = v rel • v −1Alfv > 1; top panel) at the location of the planet (vertical dashed line).Therefore, star-planet interaction cannot account for the observed LOFAR radio emission.The planet in the bottom panel is drawn at around 0.1 mJy in the y-axis, which corresponds to the 1 σ level of the LOFAR radio observations.See main text for details.

Table 1 .
Stellar parameters for GJ 1151

Table 2 .
Basic characteristics of the RV data for the different instruments and the combined dataset.

Table 3 .
Basic characteristics of the nightly binned and 3σ-clipped photometric data.

Table 4 .
Bayesian evidences of the different models calculated for the three combined RV datasets from CARMENES, HARPS-N, and HPF instruments.
Hα /L bol = −4.75.Reiners et al. (2022) report an upper limit to the average surface magnetic field of GJ 1151 of B < 680 G, comparing (for hundreds of M dwarfs) the observed broadening of different magnetically sensitive lines to the polarised radiative transfer calculations.With the Stokes V (circularly polarised) detection diagnosis Reiners et al. (2018b)21)to separate magnetically active from inactive stars; but Hα emission was detected byReiners et al. (2018b), with log L