© The Authors 2025

1 Introduction

Exoplanetary astronomy promises to be more fruitful than ever, with the number of confirmed planets recently having exceeded 5800 (Akeson et al. 2013; accessed Jan 2025). The rate of detection is still on the rise owing to many transit-photometry surveys, most notably Kepler (Borucki et al. 2010) and TESS (Ricker et al. 2015). Although now behind in number of detections, the method of radial velocities (RVs), which once initiated the search for planets with Mayor & Queloz (1995), continues to deliver. The RV method remains a favoured technique because: (i) it is less selective in terms of system geometry; (ii) it is convenient for follow-up confirmation of planet candidates, which is the case for many Kepler and TESS targets; and (iii) it allows for the masses of exoplanets to be constrained, which provides information on their interiors, as long as their radii are determined. Moreover, deep RV searches are very desirable for the nearest stars, since only a small fraction of planets are generally expected to transit.

The RV method measures the gravitational influence on a star from one or more other bodies. This suggests that M dwarfs – the least massive and most prevalent main-sequence stars in the Galaxy – are the theoretical favourites for this technique. Occurrence-rate surveys have indicated that M dwarfs host at least one planet on average (Dressing & Charbonneau 2015; Tuomi et al. 2019; Hsu et al. 2020; Sabotta et al. 2021; Perger et al. 2023). In practice, however, M dwarfs stubbornly mask potential planetary RV signals through physical processes inherent to late-type stars, most notably stellar activity on both short and long timescales (e.g. Borgniet et al. 2015 and Dumusque et al. 2011, respectively). In this work, we examine the RV behaviour of one such M dwarf, GJ 3998 (α = 17^h16^m00.6^s, δ = 11^◦03′27.6″; Gaia Collaboration 2023, hereafter GDR3). The distance to GJ 3998 has been constrained to 18.2 pc (GDR3), and its stellar parameters have been reasonably well established (Table 1). Most recently, GJ 3998 has been studied by the HADES RV programme that was carried out at the HARPS-N spectrograph at the 3.58 m Telescopio Nazionale Galileo (TNG; Cosentino et al. 2012). Affer et al. (2016), hereafter A16, modelled the stellar activity that GJ 3998 exhibited, and found evidence of two planetary companions around it: GJ 3998 b, with a period near 2.65 d, an eccentricity of zero, and a minimum mass of 2.47 ± 0.27 M_⊕; and GJ 3998 c, with a period near 13.74 d, an eccentricity of ${0.049}_{-0.034}^{+0.052}$, and a minimum mass of ${6.26}_{-0.76}^{+0.79}$ M_⊕. The rotational period of GJ 3998 was also determined by A16 to be $P_{\text{rot }}={31.8}_{-0.5}^{+0.6}$ d; soon after, a photometric confirmation followed (Giacobbe et al. 2020; 32.9925 d).

Since the discovery of GJ 3998 b and GJ 3998 c, the exoplanetary detection field has been enjoying the fruits of novel advancements. These include the use of multi-dimensional Gaussian-process regression (Rajpaul et al. 2015; Delisle et al. 2020, 2022) and more sophisticated RV determination techniques (e.g. Artigau et al. 2022). Such advancements can lead to the discovery of previously uncaught planetary signals. This circumstance and the proximity of this system constitute our main motivation for a more modern RV analysis of GJ 3998 that includes contemporary data. This work continues as follows.

Section 2 describes the available observational data relevant to our study. Section 3 narrates the generalised modelling approach that arose over the course of our investigation. Section 4 contains our analysis of available velocimetry. We discuss our results in Sect. 5, and we summarise them in Sect. 6.

2 Observations

2.1 HARPS-N velocimetry

HARPS-N is a fibre-fed high-resolution échelle spectrograph installed at the 3.58 m TNG at the Roque de los Muchachos Observatory, Spain (Cosentino et al. 2012). It has a resolving power of R ∼ 115 000 over a 380–690 nm passband. HARPS-N is contained in a vacuum vessel to limit spectral drifts induced by temperature- and air pressure variations. This ensures its stability with the course of time. HARPS-N comes with its own pipeline that provides wavelength-calibrated spectra, as well as higher-level data products, including: reduced two- and one-dimensional wavelength-calibrated spectra, crosscorrelation functions (CCFs), CCF bisector profiles and CCF- derived RVs. Like for other HADES targets, GJ 3998 was not observed with its own simultaneous thorium-argon (ThAr) lamp calibration. Because of this, our data is not corrected for the instrumental drift between the star fibre and the referencecalibration fibre. This drift is expected to have a mean variance of 1.00 m s⁻¹ (Perger et al. 2017). Thus, we anticipate a slight inflation in the RV model jitter.

2.1.1 Radial velocities

HARPS-N CCF RVs form a good baseline for analysis, but RV determination can be rendered even more precise through other techniques on raw spectra. The line-by-line method (LBL method; Artigau et al. 2022) offers a new mode of RV extraction that uses the full spectral information to measure velocities and that handles outliers in a statistically consistent manner. In addition, the LBL method provides with many complementary metrics that give insights of stellar activity, such as differential line widths, or even RVs measured from specific wavelength ranges.

We reduced 206 available raw HARPS-N spectra through the YABI portal hosted by the Centro Italiano Archivi Astronomici. YABI (Hunter et al. 2012; Borsa et al. 2015) implements the classical HARPS-N pipeline and computes CCFs for a given template mask. We used an M2-type mask to obtain CCF RVs. Then, we fed all raw spectra to the LBL method with automatic telluric cleaning and an effective-temperature estimate informed by Table 1¹. We moved forwards with LBL RVs and several curated activity indicators (Sect. 2.1.2). One point near 2 457 124.7 BJD was rejected by LBL because it did not meet its S/N > 8 requirement. We discarded another point near 2 456 694.7 BJD on account of being a strong outlier in three activity indicators. This returned a dataset of 204 points that spans from May 2013 to Sep 2023 and has a median temporal spacing of 1.0 d. This final dataset, which we carried forwards in subsequent modelling and analysis, is 50% greater in size than A16 (N = 136). It has an LBL RV root mean square (rms) of 4.28 m s⁻¹, which is noticeably lower in comparison to CCF RVs (4.98 m s⁻¹). Moreover, the RV precision appears to have improved significantly – we observe a median LBL RV uncertainty of 1.02 m s⁻¹, compared to 1.71 m s⁻¹ for CCF RVs.

2.1.2 Activity indicators

We took to considering four metrics that quantify stellar activity in their own way: CCF FWHM, the Mount Wilson S-index, the Hα line and the Na I D lines. We collectively refer to those metrics as ‘activity indicators’.

Active regions perturb the flux- and velocity fields of the stellar disc, thereby changing the shape of observed lines. Consequently, the CCF also changes in shape. Such perturbations are quantified by tracking the evolution of CCF width, depth, and symmetry, and their strength can be related to the coverage of the active regions, their contrast and the stellar v sin i. In this work, we use the CCF full width at half maximum (FWHM) as a primary activity indicator. The CCF FWHMs are automatically provided by the HARPS-N DRS; we took them and then colour- corrected them order by order following Suárez Mascareño et al. (2023). Our considered CCF FWHM time series are described by an rms of 5.77 m s⁻¹ and a median uncertainty of 3.42 m s⁻¹.

The emission intensity of the cores of Ca II H&K lines is a reliable proxy for the strength of the stellar magnetic field, and thereby of the stellar-rotation period, P_rot, for late-type stars (Noyes et al. 1984; Suárez Mascareño et al. 2015). We used the Mount Wilson S-index, $$S_{\text{MW}}=1.111\times \frac{{\tilde{N}}_{\text{H}}+{\tilde{N}}_{\text{K}}}{{\tilde{N}}_{\text{R}}+{\tilde{N}}_{\text{V}}}+0.0153,$$(1)

where ${\tilde{N}}_{\text{H}}$ and ${\tilde{N}}_{\text{K}}$ are triangular passbands centred at 3968.470 Å and 3933.664 Å, with a common FWHM of 1.09 Å; while ${\tilde{N}}_{\text{R}}$ and ${\tilde{N}}_{\text{V}}$ are rectangular passbands centred at 4001.07 Å and 3901.07 Å, with a common width of 20Å (Vaughan et al. 1978). Our considered measurements of S_MW are described by an rms of 0.133 and a median uncertainty of 0.077.

The Hα line and the Na I D lines are also good activity indicators for low-activity, low-mass stars. We computed Hα and Na I D fluxes by the definitions of Gomes Da Silva et al. (2011) and Díaz et al. (2007), respectively. In our considered measurements, the Hα indicator is characterised by a rms of 2.44 × 10⁻² and a median uncertainty of 0.27 × 10⁻². This is potentially indicative of strong activity relative to the instrumental precision. The Na I indicator has a rms of 4.71 × 10⁻³ and a median uncertainty of 1.53 × 10⁻³. Uncertainties of all activity indicators except CCF FWHM come from photon-noise propagation in their algebraic expressions.

2.2 HARPS velocimetry

HARPS is a high-resolution échelle spectrograph located at the 3.6m ESO telescope, La Silla Observatory, Chile (Mayor et al. 2003). The public HARPS RV database by Trifonov et al. (2020) contains 6 RV measurements of GJ 3998. Acquired measurements span from April 2008 to July 2014, with a median temporal difference of 1.0 d. Compared to the HARPS-N collection, HARPS data are too small in number and too isolated in epoch to justify an inclusion in our analysis.

2.3 ASAS-SN photometry

ASAS-SN is an automated photometric programme that looks out for supernovae and other transient events (Kochanek et al. 2017). It is comprised of 24 ground-based 14cm telescopes that are grouped in fours at six sites, all part of the Las Cumbres Observatory Network (Brown et al. 2013). Sites are located at: Haleakalā Observatory, HI, USA; two at Cerro Tololo Observatory, Coquimbo, Chile; South African Astronomical Observatory, Sutherland, South Africa; McDonald Observatory, TX, USA; Ali Observatory, Tibet, China. This heterogeneous distribution of the telescope network and the large number of unit telescopes allow ASAS-SN to survey the entire visible sky every night, with an uptime that is less prone to weather conditions. We present a short self-contained photometry analysis in Appendix A. Its key takeaway is that we find no signals near periods of discussed planetary companions.

3 Modelling and inference

We now describe our mode of operation that is generalised for N data sources (e.g. HARPS-N, HARPS etc.) and M physical quantities (e.g. RV, FWHM etc.). We denote measurements y_i,j, their associated uncertainties ∆y_{i, j}, and the times of measurement t_i,j, where i ∈ [0..N − 1] and j ∈ [0..M − 1]. We imposed two requirements: (i) ∀_{i, j} : avg y_{i, j} = 0, (ii) max t_{i, j} = 0. The first condition establishes common ground for the introduction of instrument offsets (Sect. 3.1). The second condition minimises uncertainty on published ephemerides, and avoids complex degeneracies between parameters of periodic long-term corrections (Sect. 3.2). Some additional preliminary steps need to be made in order to clean data of instrumental or physical trends. Then, the stellar activity was modelled together with potential planetary signals. The following subsections narrate the steps of pre-processing, stellar-activity modelling and the optional planetary modelling. Our models assume that stellar activity influences all considered physical quantities, while potential planets act on RV only (conventionally assigned to j = 0).

3.1 Source offsets and jitters

In cases where spectra were acquired from different instruments, we accounted for non-common zero-points through the substitution $$y_{i,j}\to y_{i,j}-O_{i,j}{}_{\forall j,}^{\forall i\ne 0,}$$(2)

where O_{i, j} is the offset of source i, in the quantity j. HARPS-N is the zero-point and the only source in this analysis.

Beyond instrumental precision, measurements come with additional uncorrelated noise that we are in principle unable to model. We accounted for this by adding jitter terms in quadrature to measurement uncertainties; that is, through the substitution $$\text{Δ}{y}_{i,j}\to \sqrt{{\left(\text{Δ}{y}_{i,j}\right)}^2+J_{i,j}^2}{}_{\forall j,}^{\forall i,}$$(3)

where J_{i, j} is the instrumental jitter of source i in quantity j.

3.2 Long-term correction

Radial velocity extraction pipelines sometimes provide with data that follow strong trends. Said trends may be physically motivated (e.g. secular acceleration or long-term stellar cycles), but they may also be systematics uncaught by reduction pipelines (e.g. instrumental drifts in RV). We applied corrections through the substitution $$y_{i,j}\to y_{i,j}-f_j\left(t_{i,j}\right){}_{\forall j,}^{\forall i,}$$(4)

where f_j is the long-term function (LTF) for quantity j. LTFs may or may not share parameters across quantities.

3.3 Planetary signatures

Searches for planetary signatures require comparisons between models with and without planets. For models with planets, we subtracted each planetary signal in RV through the substitution $$y_{i,0}\to y_{i,0}-k_{\text{rv}}\left\{\cos \left[\omega +v\left(t_{i,0}\right)\right]+e\cos \omega \right\}\forall i,$$(5)

where k_rv is the RV semi-amplitude of the signal, ω is the argument of periastron, ν is the true anomaly function, and e is the planetary eccentricity. Planetary orbits are characterised by five parameters: k_rv, e, ω, the period, P, and the phase, φ. We defined φ such that max(t_{i, j}) has a normalised phase φ relative to the latest inferior conjunction, as was defined in Kane et al. (2009) under the name ‘midpoint of primary transit’. Inferior-conjunction ephemerides are therefore located at $$\epsilon =(n-\phi )P\forall n\in \mathbb{Z}.$$(6)

To solve Eq. (5), we utilised the EKEPL 1 procedure by Odell & Gooding (1986), which was reported to achieve precision better than 10⁻¹² rad after five iterations. Nevertheless, it is sometimes favourable to approximate Eq. (5) to $$y_{i,0}\to y_{i,0}+k_{\text{rv}}\sin \left[\frac{2\pi \left(t_{i,0}+P\phi \right)}{P}\right]\forall i,$$(7)

which is much faster to compute, and which introduces only three parameters: k_rv, P and φ. We hereafter refer to Eq. (7) as the planetary-sine approximation, or the circular-orbit approximation. Here, φ is defined as in Eq. (5). We note that this preservation of meaning required a change in sign relative to Eq. (5).

3.4 Stellar-activity modelling

Gaussian processes (GPs) lay the foundation of our stellar- activity modelling. Formally, a GP is defined as ‘a collection of random variables, any finite number of which have a joint Gaussian distribution’ (Rasmussen & Williams 2006). As GPs offer non-parametric approaches to fit data, observations themselves determine the best functional form of the model, no matter how intricate (Haywood 2015). This makes GP-based frameworks convenient tools for analysing stellar activity in RV, because they require no prior knowledge of associated active regions. GPs have already demonstrated to reap notable successes in planet detections (e.g. Haywood et al. 2014). However, the main strength of GPs is also their greatest weakness – at times, they are flexible enough to overfit data, and in turn to suppress signals of physical nature (Suárez Mascareño et al. 2023).

Stellar activity can be modelled as a stochastic process that is characterised by a periodic component (stellar rotation), as well as an exponential component (timescale of coherence of the periodic component). For this purpose, we employed the squared-exponential periodic (SEP) kernel², which is standard treatment of our scientific problem (e.g. Haywood et al. 2014; Angus et al. 2018). The SEP kernel has the form $$k(\text{Δ}t)={\kappa }^2\exp \left[-\frac{\text{Δ}{t}^2}{2{\tau }^2}-\frac{{\sin }^2(\pi \text{Δ}t/P)}{2{\eta }^2}\right],$$(8)

where κ² is the maximum amplitude of the kernel, τ is the evolutionary timescale, P is the kernel period that we associate with the stellar-rotation period P_rot, and η describes the kernel feature complexity. The hyperparameter η is parametrised differently in the literature, and its variations have different names; for example, lengthscale or harmonic complexity. We hereafter refer to η as the ‘sinescale’ by virtue of its algebraic analogy to the timescale τ. The S+LEAF library (Delisle et al. 2020, 2022) provides with many approximations of the SEP kernel, two of which we work with: the Matérn 3/2 exponential periodic (MEP) kernel³ and the exponential-sine periodic (ESP) kernel.

GP kernels are used in different regimes. One of the commonly used ones, the multi-dimensional regime, tries to fit one GP kernel on all physical quantities at the same time. This is done by solving the system of equations $$|\begin{array}{ll}{y}_{i,j}=A_{j}G\left(t_{i,j}\right)+B_j{( \text{d}G/\text{d}t)}_{t=t_{i,j}}\hfill & \forall i,\hfill \\ ⋮\hfill & \forall j,\hfill \end{array}$$(9)

where A_j and B_j are quantity-specific fit parameters, and G(t_{i, j}) is the kernel estimation of the stellar activity at a given time. The multi-dimensional regime follows the FF′ formalism that was described in Aigrain et al. (2012) and extended in Rajpaul et al. (2015). This regime is a common approach for the disentanglement of stellar activity in the velocimetry of late-type dwarfs, such as K2-233 (Barragán et al. 2023) and GJ 9827 (Passegger et al. 2024), and it is useful to break degeneracies in cases where there is a significant correlation between physical quantities (Suárez Mascareño et al. 2020). Rajpaul et al. (2015) discussed that the FF′ formalism expects only some activity indicators to be ‘gradiented’; that is, to have a non-zero B_j term in Eq. (9) for j, 0. We used the S+LEAF MultiSeriesKernel procedure to realise the multi-dimensional regime. We set B_j = 0 for each non-gradiented quantity and excluded it from inference.

3.5 Parameter priors

We use a set of default parameter priors unless stated explicitly. This establishes a common baseline for comparison, and lightens the description of the many models we present.

Our default LTF parameter priors are informed by the solutions of the Levenberg-Marquadt algorithm (LMA; Levenberg 1944; Marquardt 1963). We ran an instance of the LMA that returns a mean µ_lm and an uncertainty σ_lm for each LTF parameter. For these parameters, we then used the priors 𝒩 (µ_lm, 200σ_lm). The large scale factor of σ_lm accounts for the minute errors that the LMA tends to return. Our default offset priors are based on the common standard deviation of all y_{i, j}, which we denote σ_j. The default prior on offsets is O_i,j = 𝒩(0, 5σ_j). This gives wide enough priors in most cases. For jitters, we used J_i,j = 𝒰_log(10⁻³, 10³), which exploits that jitters are positive by definition.

The default priors on stellar-activity- and planetary parameters were chosen to be non-restrictive. Planets have default RV semi-amplitude priors of 𝒰_log(10⁻³, 10³) m s⁻¹ and period priors of 𝒰_log(1, 10³) d. The remaining parameters, φ, e, and ω, have uniform priors that are defined by their domains. For SEP kernel hyperparameters, we used timescale priors of 𝒰_log (40, 10⁴) d, period priors of 𝒰_log (2, 200) d, and sinescale priors of 𝒰_log (10⁻², 10²). Our procedure fits one kernel on all quantities, and the latter may have positive or negative correlations; therefore, A_j, B_j ∈ ℝ. On the other hand, non- gradiented activity indicators have a single parameter, A_j. In these cases, we constrained A_j to be non-negative so as to avoid a two-fold degeneracy in Eq. (9). We used A_j = 𝒰_log(10⁻³, 10³) and B_j = 𝒰 (−10³, 10³) for gradiented activity indicators; and A_j = 𝒰_log (−10³, 10³) and B_j = 0 for non-gradiented ones.

3.6 Inference

Markov chain Monte Carlo methods can efficiently sample across the high-dimensional parameter spaces of our models. Skilling (2004) noted that while directly sampling from the likelihood function, ℒ, becomes exponentially more expensive, sampling uniformly within a bound ℒ > k is much easier – and increasing k iteratively upon reaching convergence can be used to evaluate the posterior. This strategy is now known as nested sampling. Its formalism allows it to have: (i) lower computational complexity; (ii) the ability to explore multi-modal distributions of arbitrary complex shapes; (iii) the ability to numerically evaluate the Bayesian evidence, Z. The last point, through the Bayes factor, Z₁ /Z₂, allows for two models to be compared. The Bayes factor describes the relative evidence yielded by data between two models (Morey et al. 2016). We express the Bayes factor in terms of the difference in log-evidence ∆ ln Z between models. A value of ∆ ln Z = 5, for example, would imply that data favours the model at hand by a factor of e⁵ ≈ 150. In the literature, a necessary condition for planetary detection is typically ∆ ln Z ≥ 5 in favour of the model with a planet.

We used ReactiveNestedSampler, a nested-sampling integrator provided by ULTRANEST (Buchner 2021). In every inference instance of N_param model parameters, we required 40N_param live points and a region-respecting slice sampler that accepted 4N_param steps until the sample was considered independent.

4 Analysis

4.1 Features of HARPS-N velocimetry

The correlation matrix between CCF RVs, LBL RVs and activity indicators reveals two qualities of our dataset (Fig. C.1). Firstly, CCF and LBL RVs are very well correlated with one another (R = +0.93), and their head-to-head correlation with activity indicators is quite similar (|∆R| ≤ 0.06). We interpret this as a good validation of LBL RVs. Secondly, there is a positive four-way correlation between FWHM, S-index, Hα and Na I (R ≥ 0.43).

Figure 1 shows the generalised Lomb-Scargle periodograms (GLSPs; Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009; VanderPlas & Ivezić 2015) of LBL RV and all activity indicators in a 1–4000 d period interval. There are visible longterm trends in all activity indicators (Figs. 1d,g,j,m). Second- order polynomial functions handle well these trends, and they help in revealing strong 400–700 d peaks in activity-indicator GLSPs (Figs. 1e,h,k,n). These peaks appear to preserve their location and strength in the GLSP of time series without a second-order correction, as well in the GLSP of time series excluding points after 2458 200 BJD. This could be indicative of a magnetic cycle – and as such, it may potentially need to be included in the LTF of our model. We later discuss how our LTF was adapted to handle an inclusion of a cycle (Sect. 4.2), and to what extent the evidence supports such an inclusion (Sect. 4.4).

At the shorter-period end, we identify strong peaks in RV and all activity indicators near 30.7 d (Figs. 1b,g,l,q,v). We interpret these signals to correspond the stellar rotational period $P_{\text{rot }}={31.8}_{-0.5}^{+0.6}$ d measured by A16. The strongest peak in the RV GLSP lies at P_rot (Fig. 1b). There, peaks of two other known objects draw attention to themselves: GJ 3998 b (2.65 d and its 1 d alias near 1.5 d), as well as GJ 3998 c (13.7 d). However, signals of periodicities close to the planets are observed in some activity indicators: a peak in S-index near GJ 3998 b (2.80 d; Fig. 1h), and peaks in all activity indicators near GJ 3998 c (13.6 d and 14.1 d; Figs. 1e,h,k,n).

The presence of such activity signatures near planetary periods is problematic. It may suggest stellar-activity origin for these RV signals which had been considered planetary until now. However, activity signals can be independent from the RV signal they neighbour; for example, they can be aliases of activity signatures at other periods. The window function (WF) of our dataset reveals strong signals near 1 yr and its harmonics (Fig. 1p)⁴. We expect that the common harmonics of a true signal, P_signal, and a WF peak at P_wf come forth as aliases through $$P_{\text{alias }}^{-1}=\left|P_{\text{signal }}^{-1}\pm P_{\text{wf}}^{-1}\right|,$$(10)

as was described in Dawson & Fabrycky (2010). We can therefore identify the most prominent WF peaks, and harmonically combine them with RV peaks to check if the latter are independent from their neighbouring activity peaks. We do this exercise for GJ 3998 b and GJ 3998 c in Appendix B, and show that their closest periodicities in activity can be explained by aliasing. We later continue to test possible stellar-activity affiliations to all planetary signals (Sect. 4.5.3).

One more peak catches the eye in the RV GLSP in Fig. 1b: a 41.7 d signature that we isolate in Fig. 1c. This signal has been discussed by A16, who noted its significance not only in RV, but in the S-index and Hα. The authors considered this peak to be a byproduct of stellar activity, with the signal likely arising due to stellar differential rotation. In our time series, we do not find significant signals at 41.7 d in any activity indicator. However, we do detect strong signals in S-index (Fig. 1i) and Na I (Fig. 1o), which both lie near 41.0 d instead. Figure 2 shows our test on the dependency between the strong 41.7 d RV signal and P_rot (measured 30.7 d from the GLSP). These two signals inject aliases close to one another through the same 1/3 yr alias (Fig. 2a; upward teal triangle). Our 3800 d baseline, however, allows us to resolve the signal and the alias well. We see a clear signal at 41.7 d that cannot come from an alias of the 30.7 d one (i.e. P_rot). Rather, there are two independent signals at 30.7 d and 41.7 d, with each signal having its own family of aliases. There are additional RV structures in the intervals 33–34 d and 37–38 d that may be explained as other aliases of 30.7 d and 41.7 d. Finally, we remark that one of the 1/3 yr aliases of 30.7 d comes close to 41.0 d, roughly where we observed peaks in S-index, Hα and Na I (Figs. 1i,l,o).

In order to explore the origin and consistency of the 41.7 d RV signal, we constructed stacked GLSPs similar to Mortier & Collier Cameron (2017), but in terms of the false-alarm probability (FAP; Baluev 2008). Then, we sought the strongest peak within 0.2 d from 41.7 d and computed its FAP, measurement by measurement. Figure 3 shows this FAP evolution in RV and activity, and can be considered similar to a 1D reduction of Fig. 2 in A16. Figure 3 reminds of features that are also seen in A16: sudden gradual increase in signal from measurements 70–80, with a FAP peak near measurement 130. We marked the boundary of A16 relative to our dataset in terms of time. Until the boundary, the strongest 41.5–41.9 d S-index signal had a <0.1% FAP. As soon as new measurements enter the dataset, its significance falls to >1% and even >10% FAP. This behaviour is mimicked by Na I, but to a lesser extent. CCF FWHM and Na I never attain <10% FAP signals. In contrast, LBL RVs gradually increase in significance from about measurement 70 until the very end.

The significance of 41.5–41.9 d activity signals eventually diminishes. Yet, the strength of these signals until the end of A16 begs an explanation. At that point, the baseline of measurements is almost 900 d, or less than a quarter of our current dataset (3800 d). With every additional measurement, this baseline increased little by little. Subsequently, the 41.7 d peak and its neighbouring P_rot alias would become narrower and narrower until they completely separated – and the P_rot alias would converge outside the narrow 41.5–41.9 d interval. It is the gradual departure of the P_rot alias away from the 41.5–41.9 d region that causes a decline in significance for activity signals. GLSPs of broader period intervals (e.g. within 2 d from 41.7 d) show qualitatively the same FAP evolution, but the drop in significance for activity-index signals takes place at a larger number of measurements. This is consistent with expectations – the P_rot alias needs a larger baseline to converge away from a wider period interval.

In an idealised case of negligible stellar activity, the FAP of a planetary signal will rise indefinitely in a stacked GLSP by design. In reality, the steady trend of increasing RV significance in Fig. 3 is marred by a sudden drop between measurements 135–170. GLSPs of broader period intervals preserve the shape of this FAP evolution, suggesting that the strong RV signal is not drifting away from 41.7 d. We found a transfer of powers from the 41.7 d signal to its 1 yr alias (Fig. C.2; dashed white rectangles). This transfer occurs approximately between measurements 135–170, and would explain why we see such a strong deviation from the expected growth in significance. We constructed stacked GLSPs for GJ 3998 b and GJ 3998 c as well for comparison (Figs. C.3 and C.4). The FAP evolution of GJ 3998 c has a similar drop in significance compared to GJ 3998 d.

Fig. 1 Time series of mean-subtracted: (a) LBL RV, in m s−1; (d) CCF FWHM, in m s−1; (g) S-index; (j) Hα; (m) Na I. Measurements are given in hollow blue markers. The temporal coverage of A16 is marked by a dashed black line. The time series appear to have long-term trends, which we remove with a second-order polynomial fit (solid black line) for the purposes of exploration. (b,e,h,k,n) Associated wide-period GLSPs of LBL RV and activity indicators after subtracting the polynomial fit. Three FAP levels: 10%, 1% and 0.1%, split GLSP ordinates in bands of different colour. We highlight periodicities of the two planets discovered by A16: Pb and Pc, as well as the stellar rotational period Prot (dashed grey lines). (c,f,i,l,o) Zoomed-in GLSPs that focus near 41.7 d. (p) WF of data. (q) Zoomed-in WF that focuses near 41.7 d.

Fig. 2 Alias analysis of long-term detrended LBL RV. We focus on the aliases of the peaks at: (a) 30.7 d, the Prot peak found in the raw-data GLSP; (b) 41.7 d. (c) The WF reveals five prominent peaks at: 120.6 d (upward teal triangle), 194.0 d (pink hexagon), 339.5 d (blue circle), 369.6 d (green square) and 410.9 d (downward red triangle). Relevant aliases in data periodograms are plotted with dashed lines of a corresponding colour and marker. There are other strong peaks in the WF that do not create aliases in the specified period range (hollow grey circles).

Fig. 3 FAP evolution of the strongest peak within 0.2 d of the 41.7 d signal for long-term detrended: LBL RV (blue circles), CCF FWHM (green squares), S-index (downward red triangles), Hα (upward teal triangles) and Na I (pink hexagons). The temporal coverage of A16 is marked by a dashed black line.

4.2 Model comparison

We adopted the framework in Sect. 3 to build a model grid that probes: (i) different stellar-activity kernels; (ii) planetary systems of different number and type, hereafter ‘planetary configurations’; and (iii) the possible presence of a sine component in the LTF. For stellar activity, we chose to go with three approximations of the SEP kernel: MEP, ESP with 2 harmonics (hereafter ESP2), and ESP with 3 harmonics (hereafter ESP3). The rank of said kernels increase in this order, and they can be regarded to grow in complexity likewise. In the design of different planetary configurations, we restricted ourselves to cases where inner planets may have circular orbits, while outer planets may have Keplerian orbits – but not vice versa. We started out by modelling for two to three planets. This resulted in 7 unique configurations: (i) two circular planets; (ii) circular inner and Keplerian outer; (iii) two Keplerian planets; (iv) three circular planets; and so on. In addition to those 7, we included a planet-free configuration and a one-circular configuration, so as to confirm the inner two planets discovered by A16. This gave 9 planetary configurations in total. In summary, the 3 stellar-activity kernels (MEP, ESP2, ESP3), the 9 planetary configurations and the presence/absence of a sine component in the LTF resulted in a grid of 54 models.

Models were fitted on the conjunction of LBL RV & CCF FWHM time series. We modelled stellar activity in the multidimensional regime, with FWHM being non-gradiented. We used our default SEP hyperparameter priors (Sect 3.5), even though P_rot had already been constrained in A16. For the three modelled planets, we used the period priors: 𝒰 (2, 3) d, 𝒰 (13, 14) d and 𝒰 (30, 50) d. The priors of the first two planets are informed by A16, and the third-planet prior tests for the 41.7 d signal, while allowing it to diverge to P_rot. Planetary semi-amplitudes were assigned priors of k_rv = 𝒰 (0, 5) m s⁻¹ ; and in Keplerian orbits, they took eccentricity priors of e = 𝒰 (0, 1). For models with a sine component in the LTF, we used $$f_j\left(t_{i,j}\right)=k_{\text{cyc},j}\sin \left[\frac{2\pi \left(t_{i,j}+P_{\text{cyc}}{\phi }_{\text{cyc},j}\right)}{P_{\text{cyc}}}\right]+{\alpha }_j{t}_{i,j}^2+{\beta }_j{t}_{i,j}+{\gamma }_j,$$(11)

where k_{cyc, j} and φ_{cyc, j} are the cycle amplitude and cycle phase of quantity j, and P_cyc is the common cycle period. We assigned amplitude priors of 𝒰 (0, 40) m s⁻¹, phase priors of 𝒰 (0, 1) and common-period priors of 𝒰_log (200, 2000) d. For models with a sine-free LTF, we set k_{cyc, j} = 0 and excluded k_{cyc, j}, φ_{cyc, j} and P_cyc from sampling.

Figure 4 provides with the ln Z results of our model-grid search. Entries are arranged such that planetary configurations increase in complexity downwards, while stellar-activity kernels increase in complexity rightwards. Dashed lines split the model grid in number of planets (from zero to three); as well as whether a sine term was present in the LTF. We add a chess-like coordinate system at the grid borders, which we later use to refer to a particular model for brevity. Finally, we highlight the most appropriate (‘best’) model in our analysis: d4, a model with a sine component in the LTF, three circular planets and a MEP stellar-activity kernel. We carry model d4 over to a complete analysis, and we discuss the implications of its results in Sect. 5.

Later in this work, we compare groups of models that share common features; for example, planet-free models (a9-f9) against one-planet models (a8-f8). We are generally interested in the improvement from a group of less complex models, A, to a group of more complex models, B. We quantify said improvement through two quantities that relate to the Bayes factor. The first quantity is the average improvement in evidence, $${(\text{Δ\hspace{0.17em}ln\hspace{0.17em}}Z)}_{\text{avg }}=\mathrm{avg}\left(\ln Z_{\text{B}}\right)-\mathrm{avg}\left(\ln Z_{\text{A}}\right),$$(12)

with which we report the two groups that formed these averages. The second quantity is the pessimistic improvement in evidence, $${(\text{Δ\hspace{0.17em}ln\hspace{0.17em}}Z)}_{\text{pes}}=\min \left(\ln Z_{\text{B}}\right)-\max \left(\ln Z_{\text{A}}\right),$$(13)

with which we report the models that yielded the ln Z extrema of each respective group. (∆ ln Z)_pes is defined to favour models from A, so it is naturally more conservative than (∆ ln Z)_avg. We use both quantities in tandem, but we regard (∆ ln Z)_pes to better address the true ln Z uncertainty that has been observed in the literature (e.g. Nelson et al. 2020). This comes at the cost of increased false negatives when we try to add new features in models.

Fig. 4 Bayesian-evidence comparison in our model grid. (a) Comparison between different planetary configurations and stellar-activity kernels, as well as for a potential presence of a sine component in the LTF (Sect. 4.2). Planetary configurations include planets with circular (C; e = 0) or Keplerian (K; e ≥ 0) orbits. For stellar activity, we utilised the MEP kernel and the ESP kernel with 2 or 3 harmonics. The model that we further elect for analysis assumed a sine in the LTF, three circular orbits and a MEP kernel (d4, red border; ln Z = −1114.4). We give the Bayesian factor ∆ ln Z of remaining models relative to model d4. (b) ∆ ln Z comparison of GP-free planet-free models that only model the LTF. We use this to break the tie between inclusion/exclusion of a sine in the LTF.

4.3 Evidence of two inner planets

There are striking differences in evidence between planet-free models and two-planet models (a9-f9 vs. a5-f7; Fig. 4a). If we consider models with sine-free LTFs, (∆ ln Z)_pes = 32.5 (b9 vs. c5), implying that data favours two planets by a factor of e^32.5 > 10¹⁴. The inclusion of a sine component in the LTF marginally hinders the improvement to a two-planet model – we report a minimum (∆ ln Z)_pes = 30.3 (e9 vs. f5). One-circular models serve as an intermediate point in our analysis, and their ln Z values fit nicely within our expectations for a smooth increase in evidence. On these grounds, we support the interpretation of A16 and validate the planetary nature of GJ 3998 b and GJ 3998 c. We note strong evidence in favour of circular orbits. For models with sine-free LTFs, and using the two-circular model group for reference, we report the following improvements: (∆ ln Z)_avg = −3.5 to a circular inner and a Keplerian outer planet (a7-c7 vs. a6-c6); as well as (∆ ln Z)_avg = −4.6 to a two-Keplerian model (a7-c7 vs. a5-c5). For models with a sine component in their LTFs, and using the same group for reference, we report: (∆ ln Z)_avg = −4.3 (d7-f7 vs. d6-f6); as well as (∆ ln Z)_avg = −4.6 (d7-f7 vs. d5-f5). A16 published a circular GJ 3998 b and a Keplerian GJ 3998 c, but we point out that the eccentricity posterior of the latter was already compatible with the zero (${0.049}_{-0.045}^{+0.094};2\sigma$).

4.4 Evidence for a long-term cycle

Our data supports GJ 3998 b and GJ 3998 c regardless of whether the model LTF contains a sine component (Sect. 4.3). It is however not evident from the model grid whether data favours a sine component in the LTF on the whole. For this reason, we prepared two intermediate models that simply fit a LTF instead of removing it from the data; consequently, no GPs or planets were included. Figure 4b compares them head-to-head, aside from the main model grid. The inclusion of a sine component in the LTF increases the evidence by ∆ ln Z = 17.1.

The (∆ ln Z)_avg improvements from models with a sine-free LTF to models with a sine in their LTF are: −3.7 for planet-free models (a9-c9 vs. d9-f9); −7.1 for one-circular models (a8-c8 vs. d8-f8); −4.1 for two-planet models (a5-c7 vs. d5-f7); −2.8 for three-planet models (a1-c4 vs. d1-f4). Unlike the comparison in Fig. 4b, data do not favour a sine in the LTF as soon as a GP is involved. This suggests that our stellar-activity kernels fit for this long-term trend, which is beyond their intended purpose. It is therefore best to use the results from Fig. 4b to break the tie between LTFs with and without sines. From now on, we continue to consider only models that include a sine in their LTF (i.e. d1- f9), unless stated otherwise.

We now turn to the posteriors of the sine components in LTFs. Figure C.5 compares the posterior distributions of the period P_cyc, the RV amplitude k_{cyc, 0} and the FWHM amplitude k_{cyc, 1} between all relevant models. Distributions are normalised to have equal height for the sake of readability. Posteriors agree well across two- and three-planet models. They give an overall impression of a P_cyc ≈ 300 d cycle with a well-defined RV amplitude of 0.8–1.4 m s⁻¹ and a moderately defined FWHM amplitude of 0.4–3 m s⁻¹. ASAS-SN Sloan g photometry is characterised by strong peaks near P_cyc and 2P_cyc as well (Appendix A; Fig. A.2). This may suggest a photometric manifestation of the same cycle.

Beyond that, the full picture is very different. Planet-free and one-circular models catch the 300 d periodicity, but a 500–600 d mode tends to dominate in the posterior, and sometimes a third mode near 2000 d manifests itself. We also observe differences between kernels – for example, ESP-kernel k_{cyc, 1} posteriors tend to have heavy left tails overall. The latter phenomenon is especially pronounced for ESP3 kernels, for which the mode does not seem to be well separated from the zero. This may come from the higher complexity of the ESP3 kernel, which tries to assign more long-term features to stellar activity. Finally, it may be the data that suggest a poorly constrained FWHM amplitude – we remind the reader that RV and FWHM measurements have median uncertainties of 1.02 m s⁻¹ and 3.42 m s⁻¹, respectively.

4.5 Evidence for a third planet

Three-planet models confidently stand out from their two-planet counterparts. This is true for each stellar-activity kernel taken individually. We report the following (∆ ln Z)_avg from two-planet to three-planet models with sines in their LTFs: 11.2 for MEP (d5-d7 vs. d1-d4), 14.0 for ESP2 (e5-e7 vs. e1-e4) and 14.0 for ESP3 (f5-f7 vs. f1-f4). For the same groups, we report (∆ ln Z)_pes of 1.0, 1.8, and 4.3, respectively. We can compare two- and three- planet models on a more equal footing by considering models with all-circular and all-Keplerian orbits. For all-circular orbits, the ∆ ln Z improvement is 14.9 for MEP, 14.9 for ESP2 and 14.1 for ESP3. For all-Keplerian orbits, we report ∆ln Z of 3.2, 17.6, and 8.7, respectively. All but one of these reported improvements are in principle enough to support a detection in an agnostic manner (∆ ln Z ≥ 5; Jeffreys 1946).

To test the numerical stability of third-planet solutions, we took all three planet-models, including ones with sine- free LTFs – and then individually compared the third-planet orbital-parameter posteriors between themselves. We list the comparison between P_d, φ_d, k_{rv, d}, e_d and ω_d posteriors in Tables C.1, C.2, C.3, C.4, and C.5, respectively. All models successfully recover an RV signature of period 41.77–41.79 d and of semi-amplitude 1.74–1.91 m s⁻¹. Posteriors agree in uncertainty for each orbital parameter, and for any model pair.

Among three-planet models with sines in their LTFs, the three-circular MEP model is the simplest, both in terms of planetary- and in stellar-activity complexity. Some other models fare better in evidence, but not by much (∆ ln Z ≤ 2.3). We therefore select d4 as our best-fitting model, or best model for short. Figure 5 highlights the agreement between model d4 and data in a selected dense interval with 32 clustered data points. As d4 models for three planets, the model fit time series expectedly shows much more complexity in RV. Nevertheless, data appear to be fitted well within uncertainty in both quantities. The temporal sampling of our measurements is non-trivial, as is hinted at by the WF (Fig. 1p). There are clusters of more than ten measurements within 20 d (e.g. 2 457 140–2 457 160 BJD), but also lone measurements that are isolated in time by several tens of days (e.g. near 2456 750 BJD). We inspected the model fit in the full time series, and the MEP kernel appears to have no problem with modelling activity on short and long timescales.

Figure 6 supplies with the RV and FWHM GLSPs of the residual time series, having accounted for the model jitter. Both periodograms show no significant periodicities except a 3% FAP RV signal at 1.30 d and a 5% FAP RV signal at 4.28 d. At the same time, they do not hint of overfitting. Figure 7 displays the parameter posteriors of the three planets, and Fig. 8 shows their isolated RV signals from modelled trends and stellar-activity. Data agrees well with each modelled sinusoidal signal, including the novel GJ 3998 d. We observe ordinary residual distributions and no apparent trends with phase.

Table C.6 compares our posteriors with the published results by A16. Although not all common parameters could be compared between A16 and our work, we find reasonable agreement for the ones that do. Our RV jitter $J_{0,0}={1.29}_{-0.14}^{+0.15}$ m s⁻¹ agrees with A16 (${1.19}_{-0.14}^{+0.11}$ m s⁻¹), but is marginally larger. This could be indicative of the lack of instrumental-drift correction (Sect. 2.1). Our stellar-rotation period P_rot = 30.2 ± 0.3 d somewhat agrees with A16 (${31.8}_{-0.5}^{+0.6}$ d). The posteriors of GJ 3998 b and GJ 3998 c are in excellent agreement with A16, with one exception: our inferior-conjunction ephemerides have uncertainties that are three to seven times as large. This may come from the choice of ephemeris sampling: we sampled in normalised phase, whereas A16 did so in epoch directly. Our third-planet posteriors appear well defined in model d4. The RV semi-amplitude of this new planet (${1.74}_{-0.25}^{+0.26}$ m s⁻¹) stands nearly 7σ away from the zero. Its phase is well defined, and has a similar absolute uncertainty compared to the other two planets.

In addition to inferior-conjunction ephemerides, we derived three other quantities in Table C.6: semi-major axes a, minimum planetary masses m sin i and mean insolation fluxes Φ. The first two quantities necessitate knowledge of the stellar mass M_⋆, while the third quantity requires M_⋆ and the stellar luminosity L_⋆. Through stellar parameters from Table 1, we get compatible results for the semi-major axes and the minimum masses of the inner two planets. For the third planet, we obtain $m_{\text{d}}\sin i_{\text{d}}={6.07}_{-0.96}^{+1.00}$ M_⊕. We report ${\text{Φ}}_{\text{b}}={48}_{-9}^{+11}{\text{Φ}}_{\oplus },{\text{Φ}}_{\text{c}}={5.4}_{-1.0}^{+1.3}{\text{Φ}}_{\oplus }$ and ${\text{Φ}}_{\text{d}}={1.2}_{-0.2}^{+0.3}$ Φ. We discuss the implications of these results in Sect. 5.1.

We assessed the significance of GJ 3998 b, GJ 3998 c and GJ 3998 d following Hara et al. (2022), hereafter H22. Their formalism uses the posterior distribution of an inference run to compute the probability of having no planets within a given element in the orbital-frequency space. H22 refers to this metric as the false inclusion probability (FIP). We re-ran model d4 with the following modification in priors: P_b, P_c, P_d = 𝒰(2, 50) d. We did so in order to encompass the full period range 2–50 d, which had not been fully covered by our individual priors. The formalism by H22 works in Bayesian-evidence terms, at a stage where the parameter space has been already marginalised. Therefore, having identical P_b, P_c, P_d priors should not perturb the FIP framework. We computed the FIP between 2–50 d for an angular- frequency step ∆ω = 2π/5T_span, where T_span is the observation span. Figure 9 displays our FIP calculation as a function of orbital period. Our new discovery, GJ 3998 d, is assigned a FIP of around 2 × 10⁻⁵, which passes the 10⁻² threshold suggested by H22.

Fig. 5 Raw time series (blue points) against our best model (d4; solid lines) in a selected temporal region. Data points come with two error bars: one assigned by the instrument (blue), and another that includes the model jitter (grey; Sect. 3.1). (a) LBL RV in m s−1, where the full model (black) is split into the stellar-activity GP component (green), the three-planet component (red) and the LTF (teal). (b) CCF FWHM in m s−1, where the full model (black) is split into the stellar-activity GP component (green) and the LTF (teal). Every component is offset by an arbitrary amount, labelled in the plots. The 1σ confidence intervals of the GP are given in shaded bands.

Fig. 6 GLSPs of the residual time series of our best model (d4), accounting for the model jitter, for: (a) LBL RVs, (b) CCF FWHMs.

Fig. 7 Prior (grey, cross-hatched) and posterior (blue, dot-hatched) distributions of planetary parameters in our best model (d4). Uncertainties in our posteriors reflect the 16th and the 84th percentiles.

Fig. 8 (a,b,c) Activity-corrected, phase-folded time series (blue error bars) against the planetary fit of our best model (d4; solid black line). (d,e,f) Residual time series after the planetary fit. (g) Distribution of residuals (solid blue line).

4.5.1 Modelling with RV and other activity indicators

We ran model d4 on RV in conjunction with other activity indicators: RV&S-index, RV&Hα, and RV&Na I. If the 41.7 d signal comes from stellar activity, then model d4 would fail to describe the observed time series, or would converge to a different result. We ran these models again in the multi-dimensional regime, with activity indicators being non-gradiented in the setup, similar to FWHM. We used the same priors as in model d4, but with one alteration: we used J_0,1 = 𝒰_log 10⁻⁴, 1 for Hα and Na I to account for their small absolute uncertainties. Finally, we repeated model d4, but using CCF RVs instead of LBL RVs. This was done to check if the detection of a third planet could have been made through CCF data alone.

Table C.7 compares the posteriors of model d4 against those run variations. We could only find four disagreements between comparable parameters. The first two disagreements come from the CCF RV & CCF FWHM model: (i) the LTF sine period P_cyc tries to converge near its 2000 d prior boundary, (ii) P_d contains two secondary modes approximately at 31.1 d and 31.4 d, close to the 1/3 yr alias of the 41.7 d signal. In third place comes the very unconstrained $P_{\text{cyc}}={1384}_{-945}^{+502}$ d of the LBL RV & Na I model. Finally, the LTF of the LBL RV & S-index model has flat phase posteriors.

Among these four disagreements, only the second relates to a planetary parameter. It is not fundamental and it takes place in a dataset where RV is of known lower precision. The CCF RV & CCF FWHM model has posteriors that are unsurprisingly less constrained than the ones from model d4. Nevertheless, GJ 3998 d is detected even with CCF RVs – we measure an RV semi-amplitude of ${1.52}_{-0.45}^{+0.41}$ m s⁻¹, more than 3σ away from the zero.

Fig. 9 FIP-period plot of our best model (d4; solid black line) against the 10−2 threshold suggested by Hara et al. (2022) (dashed black line).

4.5.2 Modelling with RVs from other pipelines

We used SERVAL (Zechmeister et al. 2018) to compute a third RV data source, with which we can validate CCF and LBL RVs. We reduced all available raw HARPS-N spectra through SERVAL, and we took from there only RVs and its associated uncertainties. Then, we crossmatched our dataset with SERVAL RVs by time, thereby rejecting the same outliers. This query returned 204 datapoints with a median temporal spacing of 1.0 d. SERVAL RVs have an rms of 4.25 m s⁻¹, which approaches LBL RV rms (4.28 m s⁻¹). These new RV measurements are less precise, with a median uncertainty of 1.19 m s⁻¹ (against 1.02 m s⁻¹ in LBL RVs). We ran a duplicate of model d4, but on SERVAL RVs & CCF FWHM data. We found no differences from our d4 results (Fig. C.6, Table C.6) that warrant discussion.

4.5.3 Modelling activity signals near planetary periods

Planets manifest themselves in RV only, not in activity indicators. Therefore, strong activity signals near orbital periods of planets should prompt scrutiny. This concept has been fundamental for the validation of detections since the early days of exoplanetary discovery (e.g. Queloz et al. 2001), and is still being used to refute planetary candidates (e.g. Robertson et al. 2015). In fact, GJ 3998 has already been subject of similar analysis. Recently, Dodson-Robinson et al. (2022) analysed the historical A16 data and shared concerns that the signals attributed to GJ 3998 b and GJ 3998 c may come from stellar activity.

Not all cases, however, are easy to resolve. For example, the stellar rotational period of GJ 3998 is P_rot = 30.2 ± 0.3 d, and the orbital period of GJ 3998 c is $P_{\text{c}}={13.727}_{-0.004}^{+0.003}$ d; in other words, near P_rot/2 (Table C.6). Figure 1 was earlier discussed to reveal strong peaks in activity near P_c. This would either imply that: (a) the 13.7 d RV signal is indeed a planet and the 13.7 d activity signals are harmonics of P_rot; or that (b) RV and activity signals come from the same stellar-activity component, which somehow escapes our modelling. Because case (b) would invalidate a planet, we found it necessary to devise a general test for all planets in the system.

Suppose there is an excess signal at period P in RV and activity, which our modelling fails to address in all physical quantities. Therefore, a model that accounts for stellar activity and one planet should erroneously guide the planetary period to P. However, because the excess signal manifests in activity too, if we fit a planet-like component in activity instead of RV, the period of this component should also be erroneously guided to P. We can therefore fit variants of model d4 with three sine terms in a chosen activity indicator, and check if the posteriors of those terms compare to the three planets in model d4. If posteriors show similarity for any planet, case (b) holds and that planet must be invalidated.

We searched for three sinusoidal signals in variants of model d4 that were individually fit on each activity-indicator time series: on FWHM data only, on S-index data only and so on. We used the same priors as in model d4, with the same adjustments as in Sect. 4.5.1: J_0,0 = 𝒰_log 10⁻⁴, 1 for Hα and Na I to handle the scale of their uncertainties. Figures 10, 11, 12, and 13 compare the three-planet posteriors between model d4 and models in: FWHM, S-index, Hα, and Na I, respectively. None of these comparisons show evidence to suggest a stellar-activity origin for any of the three planets. A more thorough discussion on the results follows.

A baseline proxy for the significance of a signal is its semiamplitude posterior relative to the median uncertainty of the data. We remind the reader that all three planets in model d4 have RV semi-amplitudes that are well separated from the median RV uncertainty (1.02 m s⁻¹ ; dashed blue lines in g,h,i panels of Figs. 10–13). Some of the fitted activity sine terms have inklings of modes in semi-amplitude near the median uncertainty – but they all have other clear disagreements with d4 planet posteriors. These sine terms were designed to seek activity signals in the priors of: (i) GJ 3998 d in S-index, Fig. 11i, which converges to a period between 30–31 d, close to P_rot; (ii) GJ 3998 c in Hα, Fig. 12h, which converges to a weakly defined period and phase; (iii) GJ 3998 d in Hα, Fig. 12i, which has two period modes: 30–31 d and 48–50 d, both far away from P_d; (iv) GJ 3998 d in Na I, Fig. 13i, which converges to a period between 30–31 d, close to P_rot.

On the whole, we do observe a few faint structures in the posteriors of the fitted sine terms in activity. These structures may indeed come from residual activity that had not been picked up by the stellar-activity kernel, or from the specific WF of our data – but none of them come close to the well-constrained planetary parameters of model d4. On these grounds, we conclude that this test validates the planetary nature of GJ 3998 b, GJ 3998 c, and GJ 3998 d.

5 Discussion

5.1 The planetary system of GJ 3998

We confirm the existence of planets GJ 3998 b and GJ 3998 c that were discovered by A16. Additionally, we assign the 41.7 d RV signal to a third planet, GJ 3998 d. We constrain the orbital periods of GJ 3998 b, GJ 3998 c and GJ 3998 d to ${2.65033}_{-(19)}^{+(22)}$ d,⁵ ${13.727}_{-0.004}^{+0.003}$ d, and 41.78 ± 0.05 d, respectively (Table 2). These planets have minimum masses of ${2.50}_{-0.29}^{+0.30}{\text{M}}_{\oplus },{6.82}_{-0.75}^{+0.78}{\text{M}}_{\oplus }$, and ${6.07}_{-0.96}^{+1.00}$ M_⊕, respectively. All planets have well-defined posteriors (Table. 2, Fig. 7), and HARPS-N data matches nicely our best-model fit at their periods (Fig. 8). Our best-model residual time series reveal no significant signals in RV nor FWHM (Fig. 6). All planets likely have minute eccentricities, but our current data and modelling do not permit us to provide lower limits.

We computed the mass-period and the mass-flux relationships of known exoplanets that have had their minimum masses measured independently, and that have relative mass uncertainties smaller than one third (Akeson et al. 2013). Figure 14 presents said diagram for the period range 0.3–500 d, the flux range 0.2–900 Φ_⊕ and the mass range 0.3–50 M_⊕ against our posteriors of the three planets in the GJ 3998 system. With regards to the mass-period diagram, GJ 3998 c and GJ 3998 dlie in the bulk of discovered planets. GJ 3998 b somewhat strays from the main distribution (Fig. 14c). It is in the mass-flux diagram, however, where our three planets stand out – the minimum masses of all put them in less explored fields in the mass-flux space (Fig. 14d). This is especially the case for GJ 3998 d, which turns out to be one of the few known planets that receive an Earth-like flux. In terms of the exoplanetary statistics of similar-mass stellar hosts (±0.1 M_⊙), GJ 3998 d would be the second planet of Earth-like flux after GJ 3293 b (Astudillo-Defru et al. 2015).

The favourable flux of GJ 3998 d raises questions about its potential habitability. A traditional choice to assess the habitability of a planet is through Kopparapu et al. (2014). Their work provides the runaway- and maximum-greenhouse flux limits, which we used to define a conservative HZ (cHZ). In the same fashion, the authors also provide recent-Venus and early- Mars flux limits, which we take as boundaries of an optimistic HZ (oHZ). We used M_⋆, log₁₀ L, T_eff from Table 1 and generated 10⁵ estimations of 5 M_⊕ cHZs and oHZs. Then, we sorted those estimations by orbital period and plotted them against the best-model periods of our three planets. Figure 15 displays the results of this exercise. Our measured P_d crosses the cHZ, the oHZ, and the uninhabitable zone (uHZ). The length of the P_d line in each zone is proportional to the probability of said planet residing in said zone. Through our HZ estimations, we report the following probabilities of P_d belonging in the cHZ, the oHZ, and the uHZ: 17%, 67%, and 16%, respectively. We stress that our form of assessment assumes a rocky-planet composition for GJ 3998 d. This would imply a true mass of 5 M_⊕ at most. Having assumed a random orbit orientation, the probability of this mass inequality is 5%.

Fig. 10 Planetary-parameter posteriors of our best model (d4; Fig. 7; blue dot-hatched) against the posteriors from our FWHM model with three activity sine terms (green square-hatched). Parameters share period- and phase priors across models. We use semi-amplitude panels (g,h,i) to show the median uncertainty of RV and FWHM measurements (dashed lines in respective colour).

Fig. 11 Identical to Fig. 10 but making a comparison between our best model (d4; Fig. 7; blue dot-hatched) and our S-index model with three activity sine terms (pink star-hatched). To ease comparison, semi-amplitudes and the median uncertainty in Na I were scaled by 20.

Fig. 12 Identical to Fig. 10 but making a comparison between our best model (d4; Fig. 7; blue dot-hatched) and our Hα model with three activity sine terms (yellow circle-hatched). To ease comparison, semi-amplitudes and the median uncertainty in Hα were scaled by 300.

Fig. 13 Identical to Fig. 10 but making a comparison between our best model (d4; Fig. 7; blue dot-hatched) and Na I model with three activity sine terms (red rhomboid-hatched). To ease comparison, semi-amplitudes and the median uncertainty in Na I were scaled by 800.

5.2 Detection limits for additional companions

We assessed our detection limits for additional planetary companions with a procedure similar to Suárez Mascareño et al. (2023). We injected a random planetary-sine signal as is described in Eq. (7), following the distributions P = 𝒰_log(1,1000)d, k_rv = 𝒰_log(0.1,5) m s⁻¹, and φ = 𝒰(0, 1). Then, we used the median posterior of our best model (d4; Table C.6) to fit d4 onto the modified time series. The residual time series after this fit should be void of long- and short-term stellar activity, as well as of RV signals from the three planets – but they may or may not contain the injected signal. We then constructed a GLSP on the residual time series and computed the FAP of the strongest peak in the interval 1–1000 d of 1000 log-spaced points⁶. We repeated this operation until we obtained the FAPs of 2 × 10⁶ random planetary-sine signals.

Figure 16 displays the median FAP of all injected signals in 200 × 200 bins of the aforementioned parameter space. The results of our detection-limit grid provide with several pieces of information. Firstly, with current data, we are already sensitive to the detection of additional Neptune-mass planets at periods up to 10³ d. Secondly, there tend to be nearly no <10% FAP signals for semi-amplitudes smaller than 1 m s⁻¹ (Fig. 16a). This is in line with what we expect from the uncertainty of LBL RV data. Thirdly, we observe a significantly suppressed detection near P_rot and P_rot/2. This result, although expected, is very accentuated in our case: the detection limit goes to 3 m s⁻¹ near P_rot and degrades to even higher values for P_rot/2. Finally, we observe what seems to be a forest of narrow drops in detection between 1–4 d. This may be related to the particular WF of our dataset.

All discussed planets comfortably lie in the <10% FAP zone of the parameter space – and the posterior of GJ 3998 c is only narrowly missed by the detection drop associated with P_rot/2. Generally, detection drops can be mitigated in regimes where stellar-activity contributions in RV are not as strong. In this regard, there have been some indications that velocimetry in the near infrared may be a way in the right direction (e.g. Carmona et al. 2023). In any case, we must remark that our detection-limit test is pessimistic because it adheres to the assumptions of model d4 and the posteriors we obtained from its inference.

Fig. 14 The physical parameters of the GJ 3998 planetary system against the literature. Centre: Mass-flux and minimum mass-period diagrams of exoplanets with well-established masses and periods from Akeson et al. (2013) (red squares for M-dwarf hosts, blue circles otherwise). Against those, we plot the derived posteriors of our best model (d4; black error bars) using stellar properties from Maldonado et al. (2020). Reported uncertainties reflect the 16th and the 84th percentiles of our posteriors. Top and right sides: one-dimensional histograms of well-established planets from Akeson et al. (2013) (solid red lines for M-dwarf hosts, dashed blue lines otherwise). The locations of GJ 3998 b, GJ 3998 c and GJ 3998 d are marked by dashed black lines.

Fig. 15 Optimistic (yellow, north-east hatching) and conservative (green, north-west hatching) 5 M⊕ HZs from Kopparapu et al. (2014) against orbital periods of GJ 3998 b, GJ 3998 c and GJ 3998 d (vertical dashed lines). We show 105 separate HZ estimates that account for the uncertainty on GJ 3998 stellar parameters by Maldonado et al. (2020).

5.3 Stellar activity of GJ 3998

A16 measured log R′_HK = 5.01 for GJ 3998, which implied P_rot = 34.7 ± 6.9d through the relation in Suárez Mascareño et al. (2015). More recently, Suárez Mascareño et al. (2018), hereafter SM18, reported P_rot = 33.6 ± 3.6 d. In our best model (d4; RV&FWHM), we measure P_rot = 30.2 ± 0.3 d, which agrees with SM18, but somewhat disagrees with A16 (${31.8}_{-0.5}^{+0.6}$ d). The mismatch with A16 is likely caused by our consideration of FWHM measurements together with RVs – and ultimately, we work with more data than what was available at the time. In the variations of our best model for different activity indicators (Sect. 4.5.1; Table C.7), we measure similar values of P_rot: 30.2 ± 0.2 d for RV&S-index, 30.1 ± 0.3 d for RV&Hα, and 29.9 ± 0.3 d for RV&Na I.

The timescales of active regions do not surprise either – we report τ values of: ${180}_{-55}^{+93}$ d for the main RV&FWHM model (d4), ${321}_{-121}^{+278}$ d for RV&S-index, ${130}_{\text{−34}}^{\text{+52}}$ d for RV&Hα, and ${153}_{-41}^{+65}$ d for RV&Na I. The ratio τ/P_rot appears to be roughly between 4 and 11, depending on the activity indicator. Therefore, GJ 3998 would likely have a light-curve morphology that is not Sun-like, as is defined in Giles et al. (2017). The same work analysed the Kepler photometry of about 2200 stars with rotation periods near 10 d or 20 d. The authors found the following relation between active-region timescales and photometric rms: $$\begin{array}{ll}{\log }_{10}\tau =\hfill & 10.9252+3.0123{\text{\hspace{0.17em}log}}_{10}\text{rms}\hfill \\ \hfill & +0.5062{\left({\log }_{10}\text{rms}\right)}^2-1.3606{\log }_{10}{T}_{\text{eff}},\hfill \end{array}$$(14)

where rms is measured in mag. Taking the confidence interval $\tau ={180}_{-55}^{+93}$ from model d4, using T_eff from Table 1 and solving for the rms gives an interval of 11–21 mmag. Although Eq. (14) was derived from shorter-period stars, this rms interval somewhat compares with the raw ASAS-SN photometry (15 mmag rms in V, 34 mmag rms in g; Appendix A).

The sinescales of active regions tend to be similar and span between 0.59–0.93. We observe positive A₀, B₀, A₁ for all RV-indicator pairs. This suggests that stellar-activity modulations correlate positively across RV and activity indicators. Moreover, A₀ ≪ B₀ in all cases, which hints that stellar activity manifests itself more clearly in the RV gradient. This likely corresponds to a flux effect in what we are observing. The positive value of B₀ indicates that active regions are dark spots on the stellar surface.

Our data suggest no major differences between activity indicators in the short-term stellar activity of GJ 3998. The same could be stated for the long-term activity of GJ 3998, but with slightly reduced confidence. In our main model (d4), we measure what appears to be a sine component in the LTF, with a period $P_{\text{cyc}}={316}_{-8}^{+14} \text{d}$. This result disagrees with SM18 (1.8 ± 0.4 yr). To obtain this result, SM18 inferred double-sinusoidal models at P and P/2, following Berdyugina & Järvinen (2005). We note that the ratio between our and the SM18 P_cyc is 0.48 ± 0.11. This could mean that their double-sinusoidal model had potentially picked up the strong 316 d signal, but had potentially assigned it to P_cyc /2 – which would lead to the report of a twice as large period. This sine component in the LTF could correspond to a magnetic cycle that would be very short compared to other stars in the literature (e.g. Suárez Mascareño et al. 2016; Distefano et al. 2017; SM18). Other activity indicators inform us of a similar result. In our main model (d4), the sine component in the LTF has an RV amplitude of ${1.17}_{-0.37}^{+0.36} \text{m} {\text{s}}^{-1}$ and an FWHM amplitude of ${2.17}_{-1.12}^{+1.04} \text{m} {\text{s}}^{-1}$. We observe a somewhat significant amplitude in FWHM, but not so much in other indicators – yet, our results suggest that the sine component in FWHM, Hα and Na I are in phase⁷. This is compatible with the observed strong correlation between those indicators in the raw data (Fig. C.1).

Finally, we fitted one-dimensional variants of d4 on each activity indicator for completeness. Table C.8 supplies with their posteriors. They are overall noisier than the ones in Table C.7, but we found no differences that merit further discussion.

Fig. 16 Detection-limit grid of 2.4 × 106 random planetary-sine signals of RV amplitude 0.1–5 m s−1 and of period 1–1000 d. The grid is formed by 200 × 200 log-spaced bins. Stellar rotation complicates detection considerably for orbital periods near Prot (solid blue line) and Prot/2 (dashed blue line). We display our best-model solutions (d4; black error bars) alongside solutions of Neptune- and Earth-mass planets (dotted black lines). Reported uncertainties reflect the 16th and the 84th percentiles.

6 Conclusions

We measure GJ 3998 b, GJ 3998 c, and GJ 3998 d to have orbital periods of ${2.65033}_{-(19)}^{+(22)}\text{d},{13.727}_{-0.004}^{+0.003} \text{d}$, and 41.78 ± 0.05 d, respectively; minimum masses of ${2.50}_{-0.29}^{+0.30}{\text{M}}_{\oplus },{6.82}_{-0.75}^{+0.78}{\text{M}}_{\oplus }$, and ${6.07}_{-0.96}^{+1.00}{\text{M}}_{\oplus }$, respectively; and incident fluxes of ${48}_{-9}^{+11}{\text{Φ}}_{\oplus }$, ${5.4}_{-1.0}^{+1.3}{\text{Φ}}_{\oplus }$ and ${1.2}_{-0.2}^{+0.3}{\text{Φ}}_{\oplus }$, respectively. Our results paint a pic ture of a multi-planet system of super-Earths, where the outermost companion receives about as much incident flux as the Earth. This new discovery, GJ 3998 d, is one of the few currently known planets that come close to having an Earth-like incident flux. As such, it is a welcome addition to the statistics of planetary systems in the solar neighbourhood.

We did a chain of independent tests on the planetary nature of all three planets in the system. These tests include an alias analysis of activity peaks near orbital periods, as well as model comparisons between different: (i) stellar-activity kernels, (ii) planetary configurations, (iii) activity indicators, and (iv) RV extraction pipelines. In addition, we looked for planetary-like features in activity indicators, in case our stellar-activity modelling failed to account for such potential features. All planets passed our testing suite, and their posteriors demonstrated stability in different settings. Our detection-limit test shows that if a fourth planet existed, but somehow escaped detection, it would have likely induced an RV semi-amplitude smaller than 1 m s⁻¹ (Fig. 16). Current HARPS-N data suggests that GJ 3998 may exhibit magnetic activity with a period of ${316}_{-8}^{+14}\text{d}$, which would place it among the stars with rather short magnetic-cycle periods. This longer-term stellar activity is well observed in three activity indicators: FWHM, Hα, and Na I. Finally, we provide a tighter constraint on the rotational period of GJ 3998: P_rot = 30.2 ± 0.3 d.

Acknowledgements

AKS acknowledges the support of a fellowship from the “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/DI23/11990071. AKS, JIGH, ASM, NN and RR acknowledge financial support from the Spanish Ministry of Science and Innovation (MICINN) project PID2020-117493GB-I00 and from the Government of the Canary Islands project ProID2020010129. NN acknowledges funding from Light Bridges for the Doctoral Thesis “Habitable Earth-like planets with ESPRESSO and NIRPS”, in cooperation with the Instituto de Astrofísica de Canarias, and the use of Indefeasible Computer Rights (ICR) being commissioned at the ASTRO POC project in Tenerife, Canary Islands, Spain. The ICR-ASTRONOMY used for his research was provided by Light Bridges in cooperation with Hewlett Packard Enterprise (HPE). LA and GM acknowledge support from the agreement ASI-INAF 2021-5-HH.2-2024. IR and MP acknowledge fincancial support from Spanish grants PID2021-125627OB-C31 funded by MCIU/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, PID2020-120375GB-I00 funded by MCIU/AEI, by the programme Unidad de Excelencia María de Maeztu CEX2020-001058-M, and by the Generalitat de Catalunya/CERCA programme. IR acknowledges further financial support from the European Research Council (ERC) under the European Union’s Horizon Europe programme (ERC Advanced Grant SPOTLESS; no. 101140786). MP acknowledge support from the European Union – NextGenerationEU (PRIN MUR 2022 20229R43BH) and the “Pro- gramma di Ricerca Fondamentale INAF 2023”. EGÁ acknowledges financial support from the Universidad Complutense de Madrid and the Spanish Ministerio de Ciencia e Innovación through the project PID2022-137241NB-C4[1,4]. JM and LA acknowledge support from the Italian Ministero dell’Università e della Ricerca and the European Union – Next Generation EU through project PRIN 2022 PM4JLH “Know your little neighbours: characterising low-mass stars and planets in the Solar neighbourhood”. AKS gratefully dedicates this work to V. D. Mihov for teaching him the principles of programming beyond mere coding. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. We made use of NASA’s Astrophysics Data System, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This research has made use of the SIMBAD (Wenger et al. 2000) and VIZIER (Ochsenbein et al. 2000) databases, both operated at CDS, Strasbourg, France. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. We used the following PYTHON packages for data analysis and visualisation: ASTROPY (Astropy Collaboration 2022), MATPLOTLIB (Hunter 2007), NIEVA (Stefanov et al., in prep.), NUMPY (Harris et al. 2020), PANDAS (McKinney 2010; The Pandas Development Team 2022), SCIPY (Virtanen et al. 2020), SEABORN (Waskom 2021), S+LEAF (Delisle et al. 2020, 2022) and ULTRANEST (Buchner 2021). This manuscript was written and compiled in OVERLEAF. All presented analysis was conducted on UBUNTU machines. The bulk of modelling and inference was done on the Diva cluster (192 Xeon E7-4850 2.1 GHz CPUs; 4.4 TB RAM) at Instituto de Astrofísica de Canarias, Tenerife, Spain.

Appendix A Analysis of ASAS-SN photometry

We used the ASAS-SN Sky Patrol online service⁸ to inspect the photometry of GJ 3998 at hand. At the time of access, the field of our target had been covered by ASAS-SN from October 2012 to June 2024 (2 456 200–2 460 500 BJD). Our preliminary check revealed that there are no significant contam- inators in the field that could bias photometry. However, the total proper motion of GJ 3998 approaches 0.4″/yr (GDR3), which may become comparable to the pixel size of the detector (8″). Past works addressed this issue through a computation of many separate light curves associated with smaller temporal intervals (e.g. Trifonov et al. 2021; Damasso et al. 2023). We went forwards with a similar procedure. We did split the region 2 456 000–2 460 500 BJD, including endpoints, in timestamps at every 300 BJD. Then, for each time stamp, we computed the expected position of GJ 3998 and used the Sky Patrol service to compute an individual light curve for the temporal region defined by ±150 BJD around the time stamp. Finally, we concatenated the light curves to a single data product, and performed iterative 3σ clipping by flux for each bandpass. Our final dataset contains 1110 measurements in Johnson V and 3335 photometric measurements in Sloan g.

Figure A.1 displays the final combined ASAS-SN time series, as well as a binned scatter plot between flux and lunar separation to assess the level of lunar contamination. The rms of measurements against the median photometric uncertainty is 2.3 mJy against 0.9 mJy for Johnson V; and 2.5 mJy against 0.5 mJy for Sloan g. There are visible long-term trends in both bandpasses – we observe a gradual increase in flux in Jonhson V and a likewise decrease in Sloan g. We calculated the minimum lunar separation throughout observations to be about 29.6 deg. Nevertheless, we do not observe appreciable differences in flux in low-separation regimes (Fig. A.1b,d). We thus attribute the long-term changes to an incomplete proper-motion correction.

Figure A.2 gives the 10–1000 d GLSPs of our time series in different bands. Before us stands a bouquet of P_rot, P_cyc and 1 yr harmonics, together with some higher-order harmonics of known and unknown periodicities (see labelled signals in Fig. A.2). Compared to the stellar activity we measure in RV (Sect. 5.3), Johnson V photometry has insignificant Prot signals, a somewhat significant 2P_cyc, and two strong signatures that we could not identify: at 12.6 d and at 57 d. On the other hand, Sloan g photometry contains strong P_rot and P_cyc signatures, as well as many unidentified peaks. Those include the 12.6 d and 57 d peaks from Johnson V, another 23.8 d peak, but most notably, a mysterious 141 d signal (hereafter P₇ ; see Fig. A.2). This last signal appears to comes to together with the signatures of 2P7 and 3P7 ; as well as its two yearly harmonics at 33 d and 233 d.

What remains to be verified – and what should not be affected by a long-term proper-motion drift – is whether there are transit-like structures in the time series. We constructed box least-squares periodograms (BLSPs; Kovács et al. 2002) in the period interval between 2–1000 d and a transit-duration interval between 0.01–1 d. Figure A.3 displays the result of this exercise. We find no significant peaks at the periodicities of any planet, nor at our measured rotational period P_rot or long-term cycle P_cyc. We folded the time series for our best-model P_b, P_c, P_d (Table C.6), but we could not find no apparent structures. Repeating the exercise without the final step of sigma-clipping also gives a negative result. We conclude that current ASAS-SN data supports the argument of non-transiting GJ 3998 b, GJ 3998 c and GJ 3998 d.

Fig. A.1 (a,c) GJ 3998 ASAS-SN photometry in two bands: Johnson V and Sloan g. (b,d) The same time series against the lunar separation at the times of measurement, binned every 2°.

Fig. A.2 GLSPs of GJ 3998 ASAS-SN photometry in: (a) Johnson V, (b) Sloan g. Both periodograms are saturated with many peaks. Some of them relate to the stellar rotation period Prot (solid blue line; circle), the sinusoidal cycle Pcyc in the LTF (solid red line; downward triangle) to 1 yr (solid green line; square). We plot the first few harmonics of each of these signals in dashed lines of a corresponding colour and marker. After this initial filtering, ten significant signals remain: (1) 12.6 d, unidentified; (2) 17.1 d, unidentified; (3) 23.8 d, unidentified; (4) 27.7 d, near the common harmonic of Prot and 1 yr; (5) 33 d, near the common harmonic of P7 and 1 yr; (6) 57 d, unidentified; (7) 141 d, unidentified; (8) 233 d, near the common harmonic of P7 and 1 yr; (9) 276 d, near 2P7; (10) 393 d, near 3P7.

Fig. A.3 BLSPs of GJ 3998 ASAS-SN photometry in: (a) Johnson V, (b) Sloan g.

Appendix B Alias analysis of planets b and c

We begin our alias analysis on the neighbourhood around GJ 3998 b. There are two peaks which we test for dependency: the planetary signal in RV (2.65 d; Fig. 1c) and the 2.80 d signal in S-index (Fig. 1m). We select the most significant WF peaks in the interval 10–1000 d. Three strong WF peaks come close to ⅓ yr, ½ yr and 1 yr; these are 120.6 d, 194.0 d, and 369.6 d, respectively. For the two peaks of interest, they generate aliases within 2.61–2.83 d. We select these intervals from the long-term detrended GLSPs and show them in Fig. B.1. Each signal comes with its 1 yr aliases in its own dimension. The aliases of the 2.80 d signal are not even remotely close to the orbital period of GJ 3998 b. Equivalently, the 2.65 d planetary signal does not seem to form aliases near 2.80 d.

We continue with GJ 3998 c (13.7 d) which is surrounded by strong signals in all activity indicators. We found out P_rot/2 has a ⅓yr alias near 13.6 d, close to the orbital period of GJ 3998 c. This is indeed what we observe in Fig. B.2: the aforementioned alias fits well peaks near the orbital period of the planet (Figs. B.2d,f,j; upward teal triangle) which do not agree exactly with the RV signal (compare Fig. B.2a with Figs. B.2c,e,i). There is another strong 14.1 d signal in S-index, which cannot be well explained by aliases of P_rot/2 (Fig. B.2f). This periodicity, however, is quite far away from GJ 3998 c (13.7 d; Fig. B.2e).

Fig. B.1 Alias analysis of long-term detrended signals near 2.65 d in: (a,b) LBL RV, (c,d) S-index. The left column shows aliases of the 2.65 d signal, found in RV and attributed to GJ 3998 b. The right column shows aliases of a 2.80 d signal found in S-index. (e) The WF reveals five prominent peaks at: 120.6 d (upward teal triangle), 194.0 d (pink hexagon), 339.5 d (blue circle), 369.6 d (green square) and 410.9 d (downward red triangle). Relevant aliases in data periodograms are plotted with dashed lines of a corresponding colour and marker. There are other strong peaks in the WF that do not create aliases in the specified period range (hollow grey circles).

Fig. B.2 Alias analysis of long-term detrended signals near 13.7 d in: (a,b) LBL RV, (c,d) CCF FWHM, (e,f) S-index, (g,h) Hα, (i,j) Na I. The left column shows aliases of the 13.7 d signal, found in RV and attributed to GJ 3998 c. The right column shows aliases of half the Prot peak found in the raw-data GLSP, at 15.4 d. (k) The WF reveals five prominent peaks at: 120.6 d (upward teal triangle), 194.0 d (pink hexagon), 339.5 d (blue circle), 369.6 d (green square) and 410.9 d (downward red triangle). Relevant aliases in data periodograms are plotted with dashed lines of a corresponding colour and marker. There are other strong peaks in the WF that do not create aliases in the specified period range (hollow grey circles).

Appendix C Supplementary material

Fig. C.1 Pearson correlation matrix between LBL RV, CCF RV, and discussed activity indicators.

Fig. C.2 FAP evolution of long-term detrended RV signals between 36–44 d. We plot the location of the 41.7 d signal and its aliases within the interval. Alias symbols match Fig. B.1e, Fig. B.2k and Fig. 2c. Dashed white rectangles display the power transfer from the 41.7 d signal to its 37–38 d alias family.

Fig. C.3 FAP evolution of the strongest peak within 0.2 d of the 2.65 d signal for long-term detrended: LBL RV (blue circles), CCF FWHM (green squares), S-index (downward red triangles), Hα (upward teal triangles) and Na I (pink hexagons). The temporal coverage of A16 is marked by a dashed black line.

Fig. C.4 FAP evolution of the strongest peak within 0.2 d of the 13.7 d signal for long-term detrended: LBL RV (blue circles), CCF FWHM (green squares), S-index (downward red triangles), Hα (upward teal triangles) and Na I (pink hexagons). The temporal coverage of A16 is marked by a dashed black line.

Fig. C.5 Violin plots of all models with a sine component in their TFs that show the posteriors of: (a) the cycle period Pcyc, (b) the RV amplitude kcyc, 0, (c) the RV amplitude kcyc, 0. Violin plots are arranged vertically by planetary configuration in the same order as in Fig. 4. Plots come in three colours: blue for MEP, green for ESP2 and red for ESP3. Every distribution is accompanied by its corresponding model label.

Fig. C.6 Visual extension of Table C.6: prior (grey, cross-hatched) and posterior (blue, dot-hatched) distributions of all parameters in our best model (d4). Uncertainties in our posteriors reflect the 16th and the 84th percentiles. For the sake of brevity, we show three significant digits for posterior medians, and truncate uncertainties to match the same precision.

References

All Tables

All Figures

Fig. 1 Time series of mean-subtracted: (a) LBL RV, in m s−1; (d) CCF FWHM, in m s−1; (g) S-index; (j) Hα; (m) Na I. Measurements are given in hollow blue markers. The temporal coverage of A16 is marked by a dashed black line. The time series appear to have long-term trends, which we remove with a second-order polynomial fit (solid black line) for the purposes of exploration. (b,e,h,k,n) Associated wide-period GLSPs of LBL RV and activity indicators after subtracting the polynomial fit. Three FAP levels: 10%, 1% and 0.1%, split GLSP ordinates in bands of different colour. We highlight periodicities of the two planets discovered by A16: Pb and Pc, as well as the stellar rotational period Prot (dashed grey lines). (c,f,i,l,o) Zoomed-in GLSPs that focus near 41.7 d. (p) WF of data. (q) Zoomed-in WF that focuses near 41.7 d.

In the text

Fig. 2 Alias analysis of long-term detrended LBL RV. We focus on the aliases of the peaks at: (a) 30.7 d, the Prot peak found in the raw-data GLSP; (b) 41.7 d. (c) The WF reveals five prominent peaks at: 120.6 d (upward teal triangle), 194.0 d (pink hexagon), 339.5 d (blue circle), 369.6 d (green square) and 410.9 d (downward red triangle). Relevant aliases in data periodograms are plotted with dashed lines of a corresponding colour and marker. There are other strong peaks in the WF that do not create aliases in the specified period range (hollow grey circles).

In the text

	Fig. 3 FAP evolution of the strongest peak within 0.2 d of the 41.7 d signal for long-term detrended: LBL RV (blue circles), CCF FWHM (green squares), S-index (downward red triangles), Hα (upward teal triangles) and Na I (pink hexagons). The temporal coverage of A16 is marked by a dashed black line.
In the text

Fig. 4 Bayesian-evidence comparison in our model grid. (a) Comparison between different planetary configurations and stellar-activity kernels, as well as for a potential presence of a sine component in the LTF (Sect. 4.2). Planetary configurations include planets with circular (C; e = 0) or Keplerian (K; e ≥ 0) orbits. For stellar activity, we utilised the MEP kernel and the ESP kernel with 2 or 3 harmonics. The model that we further elect for analysis assumed a sine in the LTF, three circular orbits and a MEP kernel (d4, red border; ln Z = −1114.4). We give the Bayesian factor ∆ ln Z of remaining models relative to model d4. (b) ∆ ln Z comparison of GP-free planet-free models that only model the LTF. We use this to break the tie between inclusion/exclusion of a sine in the LTF.

In the text

Fig. 5 Raw time series (blue points) against our best model (d4; solid lines) in a selected temporal region. Data points come with two error bars: one assigned by the instrument (blue), and another that includes the model jitter (grey; Sect. 3.1). (a) LBL RV in m s−1, where the full model (black) is split into the stellar-activity GP component (green), the three-planet component (red) and the LTF (teal). (b) CCF FWHM in m s−1, where the full model (black) is split into the stellar-activity GP component (green) and the LTF (teal). Every component is offset by an arbitrary amount, labelled in the plots. The 1σ confidence intervals of the GP are given in shaded bands.

In the text

	Fig. 6 GLSPs of the residual time series of our best model (d4), accounting for the model jitter, for: (a) LBL RVs, (b) CCF FWHMs.
In the text

	Fig. 7 Prior (grey, cross-hatched) and posterior (blue, dot-hatched) distributions of planetary parameters in our best model (d4). Uncertainties in our posteriors reflect the 16th and the 84th percentiles.
In the text

	Fig. 8 (a,b,c) Activity-corrected, phase-folded time series (blue error bars) against the planetary fit of our best model (d4; solid black line). (d,e,f) Residual time series after the planetary fit. (g) Distribution of residuals (solid blue line).
In the text

	Fig. 9 FIP-period plot of our best model (d4; solid black line) against the 10−2 threshold suggested by Hara et al. (2022) (dashed black line).
In the text

	Fig. 10 Planetary-parameter posteriors of our best model (d4; Fig. 7; blue dot-hatched) against the posteriors from our FWHM model with three activity sine terms (green square-hatched). Parameters share period- and phase priors across models. We use semi-amplitude panels (g,h,i) to show the median uncertainty of RV and FWHM measurements (dashed lines in respective colour).
In the text

	Fig. 11 Identical to Fig. 10 but making a comparison between our best model (d4; Fig. 7; blue dot-hatched) and our S-index model with three activity sine terms (pink star-hatched). To ease comparison, semi-amplitudes and the median uncertainty in Na I were scaled by 20.
In the text

	Fig. 12 Identical to Fig. 10 but making a comparison between our best model (d4; Fig. 7; blue dot-hatched) and our Hα model with three activity sine terms (yellow circle-hatched). To ease comparison, semi-amplitudes and the median uncertainty in Hα were scaled by 300.
In the text

	Fig. 13 Identical to Fig. 10 but making a comparison between our best model (d4; Fig. 7; blue dot-hatched) and Na I model with three activity sine terms (red rhomboid-hatched). To ease comparison, semi-amplitudes and the median uncertainty in Na I were scaled by 800.
In the text

Fig. 14 The physical parameters of the GJ 3998 planetary system against the literature. Centre: Mass-flux and minimum mass-period diagrams of exoplanets with well-established masses and periods from Akeson et al. (2013) (red squares for M-dwarf hosts, blue circles otherwise). Against those, we plot the derived posteriors of our best model (d4; black error bars) using stellar properties from Maldonado et al. (2020). Reported uncertainties reflect the 16th and the 84th percentiles of our posteriors. Top and right sides: one-dimensional histograms of well-established planets from Akeson et al. (2013) (solid red lines for M-dwarf hosts, dashed blue lines otherwise). The locations of GJ 3998 b, GJ 3998 c and GJ 3998 d are marked by dashed black lines.

In the text

	Fig. 15 Optimistic (yellow, north-east hatching) and conservative (green, north-west hatching) 5 M⊕ HZs from Kopparapu et al. (2014) against orbital periods of GJ 3998 b, GJ 3998 c and GJ 3998 d (vertical dashed lines). We show 105 separate HZ estimates that account for the uncertainty on GJ 3998 stellar parameters by Maldonado et al. (2020).
In the text

Fig. 16 Detection-limit grid of 2.4 × 106 random planetary-sine signals of RV amplitude 0.1–5 m s−1 and of period 1–1000 d. The grid is formed by 200 × 200 log-spaced bins. Stellar rotation complicates detection considerably for orbital periods near Prot (solid blue line) and Prot/2 (dashed blue line). We display our best-model solutions (d4; black error bars) alongside solutions of Neptune- and Earth-mass planets (dotted black lines). Reported uncertainties reflect the 16th and the 84th percentiles.

In the text

	Fig. A.1 (a,c) GJ 3998 ASAS-SN photometry in two bands: Johnson V and Sloan g. (b,d) The same time series against the lunar separation at the times of measurement, binned every 2°.
In the text

Fig. A.2 GLSPs of GJ 3998 ASAS-SN photometry in: (a) Johnson V, (b) Sloan g. Both periodograms are saturated with many peaks. Some of them relate to the stellar rotation period Prot (solid blue line; circle), the sinusoidal cycle Pcyc in the LTF (solid red line; downward triangle) to 1 yr (solid green line; square). We plot the first few harmonics of each of these signals in dashed lines of a corresponding colour and marker. After this initial filtering, ten significant signals remain: (1) 12.6 d, unidentified; (2) 17.1 d, unidentified; (3) 23.8 d, unidentified; (4) 27.7 d, near the common harmonic of Prot and 1 yr; (5) 33 d, near the common harmonic of P7 and 1 yr; (6) 57 d, unidentified; (7) 141 d, unidentified; (8) 233 d, near the common harmonic of P7 and 1 yr; (9) 276 d, near 2P7; (10) 393 d, near 3P7.

In the text

	Fig. A.3 BLSPs of GJ 3998 ASAS-SN photometry in: (a) Johnson V, (b) Sloan g.
In the text

Fig. B.1 Alias analysis of long-term detrended signals near 2.65 d in: (a,b) LBL RV, (c,d) S-index. The left column shows aliases of the 2.65 d signal, found in RV and attributed to GJ 3998 b. The right column shows aliases of a 2.80 d signal found in S-index. (e) The WF reveals five prominent peaks at: 120.6 d (upward teal triangle), 194.0 d (pink hexagon), 339.5 d (blue circle), 369.6 d (green square) and 410.9 d (downward red triangle). Relevant aliases in data periodograms are plotted with dashed lines of a corresponding colour and marker. There are other strong peaks in the WF that do not create aliases in the specified period range (hollow grey circles).

In the text

Fig. B.2 Alias analysis of long-term detrended signals near 13.7 d in: (a,b) LBL RV, (c,d) CCF FWHM, (e,f) S-index, (g,h) Hα, (i,j) Na I. The left column shows aliases of the 13.7 d signal, found in RV and attributed to GJ 3998 c. The right column shows aliases of half the Prot peak found in the raw-data GLSP, at 15.4 d. (k) The WF reveals five prominent peaks at: 120.6 d (upward teal triangle), 194.0 d (pink hexagon), 339.5 d (blue circle), 369.6 d (green square) and 410.9 d (downward red triangle). Relevant aliases in data periodograms are plotted with dashed lines of a corresponding colour and marker. There are other strong peaks in the WF that do not create aliases in the specified period range (hollow grey circles).

In the text

	Fig. C.1 Pearson correlation matrix between LBL RV, CCF RV, and discussed activity indicators.
In the text

	Fig. C.2 FAP evolution of long-term detrended RV signals between 36–44 d. We plot the location of the 41.7 d signal and its aliases within the interval. Alias symbols match Fig. B.1e, Fig. B.2k and Fig. 2c. Dashed white rectangles display the power transfer from the 41.7 d signal to its 37–38 d alias family.
In the text

	Fig. C.3 FAP evolution of the strongest peak within 0.2 d of the 2.65 d signal for long-term detrended: LBL RV (blue circles), CCF FWHM (green squares), S-index (downward red triangles), Hα (upward teal triangles) and Na I (pink hexagons). The temporal coverage of A16 is marked by a dashed black line.
In the text

	Fig. C.4 FAP evolution of the strongest peak within 0.2 d of the 13.7 d signal for long-term detrended: LBL RV (blue circles), CCF FWHM (green squares), S-index (downward red triangles), Hα (upward teal triangles) and Na I (pink hexagons). The temporal coverage of A16 is marked by a dashed black line.
In the text

Fig. C.5 Violin plots of all models with a sine component in their TFs that show the posteriors of: (a) the cycle period Pcyc, (b) the RV amplitude kcyc, 0, (c) the RV amplitude kcyc, 0. Violin plots are arranged vertically by planetary configuration in the same order as in Fig. 4. Plots come in three colours: blue for MEP, green for ESP2 and red for ESP3. Every distribution is accompanied by its corresponding model label.

In the text

	Fig. C.6 Visual extension of Table C.6: prior (grey, cross-hatched) and posterior (blue, dot-hatched) distributions of all parameters in our best model (d4). Uncertainties in our posteriors reflect the 16th and the 84th percentiles. For the sake of brevity, we show three significant digits for posterior medians, and truncate uncertainties to match the same precision.
In the text