© The Authors 2025

1. Introduction

An important piece in the puzzle of planetary system formation and evolution is the absence of short-period super-Earths and sub-Neptunes in the Solar System. In recent decades, these planets have been revealed as very common planetary types in our Galaxy (e.g. Zhu et al. 2018). The giants Jupiter and Saturn likely played a significant role in shaping the inner Solar System, and might constitute the root cause for the absence of super-Earths. This proposition can be further investigated in the exoplanet context, exploring the correlation between outer cold Jupiters (CJs) and short-period small planets (SPs) defined with 1 < M_p < 20 M_⊕ and P ≲ 100 d.

Large radial velocity (RV) surveys have investigated the CJ-SP connection. While a positive correlation was initially reported among CJs and SPs (e.g. Zhu et al. 2018), this result was later challenged (Barbato et al. 2018; Bonomo et al. 2023). Recent studies analysing larger samples of systems suggest a tentative dependency of this correlation with the stellar host metallicity (e.g. Bryan & Lee 2024; Bonomo et al. 2025). To date, these studies face two major limitations: small number statistics and limited systems completeness. The former can only be alleviated in the near future by discovering many new planets, employing, for instance, Gaia Data Release 4 (DR4) for long-period giant planets and PLATO for new systems of SPs. The completeness, on the other side, is limited by the precision and sampling of the RV surveys. A deeper completeness increases the precision in occurrence rates, which is crucial to reveal potential correlations between CJs and SPs. Large observational efforts paired with novel techniques of stellar activity mitigation would be required to improve completeness.

While this observational effort is ongoing, there is a straightforward way to tighten the RV detection limits within a given dataset, which involves the study of orbital stability. Based on the hypothesis that observed systems cannot have strongly unstable orbits (the probability of detecting a system transitioning to instability is very small), a planet can be considered detected if it fulfils two essential conditions: the signal in the data is statistically significant and it does not compromise the orbital stability of the system. Some works have included this dynamical constraint a posteriori on top of the data-driven detection limits, using simplified stability criteria (e.g. Bonomo et al. 2017; Barbato et al. 2023). Nevertheless, there is currently no systematic inclusion of the stability constraint on the derivation of detection limits for CJ-SP studies or in general exoplanet contexts. Therefore, gaining access to an easy-to-use tool that derives such constraints is timely.

In this Letter, we introduce the Algorithm for the Refinement of DEtection limits via N-body stability Threshold (ARDENT). It is an open-source Python package for the computation of dynamical detection limits; namely, detection limits that combine observational and dynamical constraints. Their computation is straightforward and can be applied to a wide range of systems. In Sect. 2, we introduce the code and describe how the stability constraint is added to the data-driven detection limits. In Sect. 3, we illustrate ARDENT based on TOI-1736. This is a two-planet system composed of a 7-day SP and 570-day CJ on an eccentric orbit. From this analysis, we derive the dynamical detection limits of TOI-1736, and sum up the advantages that ARDENT offers in Sect. 4.

2. ARDENT

2.1. General strategy

The completeness of planetary systems from RV surveys is quantified via mass detection limits. Commonly, these are computed from injection-recovery tests with a various sample of hypothetical planets (e.g. Barbato et al. 2018; Rosenthal et al. 2021; Bryan & Lee 2024). The mass above which a certain percentage of injected planets can be recovered (e.g. 95%) defines the 95% detection limits. Above these limits, we can firmly rule out the existence of additional planets within our threshold confidence level. In this work, we name these limits the data-driven detection limits, since they are directly related to the data at hand.

However, these limits are incomplete because they do not validate the orbital stability of the injected planets. To rule out a potential planet, one of the two following conditions must be met: (i) if it exists, the planet signal should have been confidently detected in the RV data already; or (ii) its orbit would be dynamically non-viable. In other words, below the data-driven detection limits, any unstable area of the period-mass space should be further excluded and the detection limits would then be adapted. While orbital stability is a well-established constraint to refine the orbital solutions of planetary systems (e.g. Tamayo et al. 2021; Stalport et al. 2022), ARDENT aims to export this strategy to the computation of detection limits.

The detectability of a planet in the data is supposed independent from its dynamical plausibility. Therefore, the probability of missing an additional planet can be expressed as

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{miss}}=(1-{\mathcal{P}}_d)\times {\mathcal{P}}_s,\end{array}$$(1)

where 𝒫_d is the probability of detecting the planet in the data and 𝒫_s is the probability that the resulting planetary system is stable. No matter the detectability of a hypothetical planet, it cannot exist if the resulting system is fully unstable (𝒫_s = 0 meaning 𝒫_miss = 0). Conversely, if 𝒫_s > 0, the injected planet cannot be missed only if its detection rate is maximal (𝒫_d = 1 meaning 𝒫_miss = 0). The aim of ARDENT is to compute the mass detection limits under the conditions of detectability and dynamical plausibility.

The computation of the data-driven detection limits in ARDENT is described in Sect. 2.2. By default, ARDENT computes the 95% detection limits. Above this limit curve, 𝒫_miss ∼ 0 no matter the dynamical plausibility, 𝒫_s, of the injected planets. Below this curve, 𝒫_miss ∼ 0 exclusively when 𝒫_s ∼ 0. Consequently, ARDENT only investigates the stability constraint below the data-driven detection limits. In the period-mass space where the detection limits are computed, the orbital stability depends on both parameters. For a given orbital period, the more massive the injected planet, the stronger the planet-planet gravitational interactions that can ultimately lead to orbital instability. Therefore, ARDENT searches for the minimal mass (below the data-driven limit curve) above which the resulting system is unstable. This limit mass defines the dynamical detection limits and pushes down the system completeness to lower masses. Figure 1 illustrates this general strategy. We cannot further rule out hidden planet candidates below the data-driven limit curve when they are not fully unstable, namely, when 𝒫_s > 0.

Fig. 1. Schematics illustrating dynamical detection limits computation in ARDENT. Empty circles represent the data-driven detection limits, above which 𝒫d ∼ 1. ARDENT models the dynamical influence of known planets in the system (diamond symbol), inducing an instability zone where additional planets are excluded. The resulting dynamical detection limits are shown in red. The inner panel illustrates the iterative mass limits update accounting for stability, with dM the difference in injected mass between two consecutive steps and δ the precision threshold (cf. Sect. 2.3.1).

2.2. Data-driven detection limits

The data-driven detection limits are the most common investigation in completeness studies. To compute them, ARDENT performs injection-recovery tests in the residual RV time series, namely, the time series that is left after subtracting the best Keplerian+activity model and that serves as input to the code. Planets are injected randomly in a range of periods, P, and RV semi-amplitudes, K. The range of K is expressed in units of root-mean-square (rms) to decrease the instrumental dependence and provide a default range valid for most applications. Furthermore, a uniform sampling of K results in a more efficient exploration of the period-mass space compared to a direct mass sampling (cf. Appendix A).

For each sampled P and K, ARDENT injects the planet with a set of N_p different orbital phase values, thereby performing N_p injection-recovery tests for every point in the P − K space. ARDENT injects each signal in the time series and attempts at recovering it via the computation of the generalised Lomb-Scargle (GLS) periodogram. The injected planet is considered as recovered if the most significant signal in the GLS periodogram has a false alarm probability (FAP) below 1% with a period estimation less than 5% away from the injected period. The detection rate of the sampled planet is thus given by its number of detections over the N_p tests. Repeating this process on the entire sample of P and K fills up the parameter space within which the detection limits are computed.

At the end of the injection-recovery tests, the period, mass (or K), and detection rate are used to compute the detection limits. ARDENT splits the sample into period bins. Given the period was sampled uniformly in log scale, each bin should contain approximately the same number of points if the number of injection-recovery tests is large enough. For each period bin, ARDENT fits for the cumulative distribution function (cdf) of the detection rate as a function of the injected mass. The fitted cdf is used to estimate the limit mass for a given detection threshold. By default, ARDENT returns the 95% detection limit. Repeating this process on each period bin leads to a detection limit curve in the period-mass space. These are the data-driven detection limits returned by ARDENT, which are used to investigate the stability constraint.

2.3. Dynamical detection limits

The dynamical detection limits add the stability constraint to the condition of detectability of a hypothetical planet. In this framework, the probability of missing a planet is null if either its detectability is granted or if it is dynamically not viable (cf. Eq. (1)). This second condition allows to exclude zones of the parameter space where additional undetected planets would be unstable. As such, exclusively below the data-driven detection limits can the dynamical constraints make a difference, where 𝒫_miss = 0 only if 𝒫_s = 0. Investigating the orbital stability of a planetary system can be computationally expensive. The motivation for a fast inclusion of the stability constraint is what shapes the process of dynamical detection limits in ARDENT.

2.3.1. Selective mass and period samplings

For a given set of orbital parameters, increasing the mass of an injected planet increases also the planet-planet gravitational interactions and the risk of orbital instabilities. Therefore, ARDENT searches for new mass limits below the data-driven curve that transition from stability to instability.

To proceed, a selective mass sampling is applied to each period bin. The code first tests the stability of an injected planet with a mass of the data-driven detection limit (cf. Sect. 2.3.2) and for an even sample of orbital phases. If it is stable for at least one phase angle, that point of the period-mass space cannot be excluded and the dynamical constraint does not provide further information than the RV data. Conversely, if the injected planet is unstable for all phase angles, it can be rejected and an updated mass limit has to be found. The code then tests the stability of a negligible mass (cf. Fig. 1, step 2). Either the result remains unstable (𝒫_s = 0), meaning no planet can exist at the period tested and the updated mass limit drops to zero; or the result is stable. In this case, ARDENT initiates a convergence process to find the mass transitioning from instability (step 1) to stability (step 2). At each new step, it injects planets with a mass in-between the nearest stable-unstable pair until a precision threshold on the mass, δ, is reached. By default, we set δ = 0.5 M_⊕, which is lower than the precision of most RV instruments. This process is repeated independently for each period bin. Hence, it can be easily multi-processed. A built-in ARDENT parameter defines the number of local cores desired to compute the dynamical constraints.

Regarding the scan in orbital periods, ARDENT employs an adaptive sampling. An initial loose-period binning (15 bins by default, equally spaced in log-scale) was used to compute the data-driven mass limits from the injection-recovery tests performed prior. Additionally, the sampling is significantly increased around each known planet. Specifically, ARDENT generates 100 more bins between 0.5P and 2P, then interpolates the data-driven mass limits at each new bin. Such dense sampling around the known planets aims to resolve the main first-order mean-motion resonances (MMRs), which are particularly relevant in compact systems of small planets. ARDENT computes the updated mass limit for each period bin.

2.3.2. Analytical and numerical stability criteria

For each stability test, analytical criteria were first employed to identify unambiguously stable or unstable systems. To proceed, the angular momentum deficit (AMD) was computed (Laskar & Petit 2017). It measures the difference between the angular momentum of the system and its equivalent with circular and co-planar orbits only. Within this framework, it is possible to define a threshold AMD below which even the worst redistribution of angular momentum among the planets would not lead to close encounters, ensuring the Hill stability (Petit et al. 2018). Alternatively, the AMD framework can be used to identify the region of the parameter space where MMRs overlap, leading to strong instability (Petit et al. 2017). Furthermore, a system is also classified as unstable if the initial planetary orbits cross but are not close to a MMR. We provide a detailed description of these analytical tests in Appendix B.

When a system’s stability cannot be clearly classified, further numerical investigations are carried out. ARDENT employs the REBOUND Python package (Rein & Liu 2012) with the WHFast symplectic integrator (Rein & Tamayo 2015) to perform numerical integrations over 10⁴ orbits of the outermost planet. If needed, a general relativity (GR) correction can be activated via REBOUNDx (Tamayo et al. 2020) and the gr force structure. If the system survived the simulation (i.e. no close encounter, nor escape), the semi-major axes drift is computed based on 2000 evenly spaced outputs, from which the stability of the system is estimated. We provide further details and a calibration of stability on this drift in Appendix C.

3. Application to TOI-1736

TOI-1736 is a solar analogue around which a transiting sub-Neptune and an outer giant planet have been found (Martioli et al. 2023). Investigating the multiplicity of inner SP systems in the presence of CJs is very informative. If the latter prevented the migration of SPs or their building blocks, a CJ-multiplicity correlation would be expected (e.g. Izidoro et al. 2015; Bonomo et al. 2025). This trend might be even stronger in eccentric CJ systems, where dynamically active past evolution with planet-planet scattering likely occurred (Huang et al. 2017). Refining the completeness of systems such as TOI-1736 is of prime importance in tackling these open questions.

Martioli et al. (2023) used 152 RV observations taken with the SOPHIE spectrograph to measure the planet masses. TOI-1736 b orbits the host star in 7.07 days with a mass and radius of 12.8 ± 1.8 M_⊕ and 2.44 ± 0.18 R_⊕. TOI-1736 c has a period and mass of 570.1 ± 0.7 days and 8.7${}_{-0.6}^{+1.5}{M}_{\mathrm{Jup}}$, and an orbital inclination $i_c={69}_{-17}^{+15}$ deg measured with Gaia DR3 excess noise (Lindegren et al. 2021). The orbital eccentricity is well constrained on TOI-1736 c (e_c = 0.362 ± 0.003), but it is not significant on planet b (e_b < 0.21, Martioli et al. 2023). The over-estimation of eccentricities is a common issue in RV observations, particularly for small-amplitude signals (Hara et al. 2019). In the context of ARDENT, this erroneously biases the detection limits by enlarging the instability zones around the planets. To prevent this issue, the parameter elower indicates the lower limit of the eccentricity compatible with the observations. When specified (recommended when uncertainties are significant), ARDENT automatically switches to this value to model the system dynamics. As such, the eccentricity of TOI-1736 b is considered to be null.

We retrieved the public SOPHIE RV data, and fit the time series with a model including two Keplerians and a linear drift¹, similarly to the approach presented in Martioli et al. (2023). Injection-recovery tests were performed on the residual RVs to compute the data-driven mass limits in the period range [2, 600] days. In total, 2000 samples in the P − K space were generated, each composed of ten tests with different orbital phases with an even spread in [ − π, π[. We retrieved the 95% data-driven detection limits from these tests, which were then used together with the known planet parameters to compute the dynamical detection limits. Using the default ARDENT settings and four cores on a MacBook Pro equipped with the M2 chip, the dynamical limits computation took about 20 minutes. The result is presented in Fig. 2. Appendix D details the command lines.

Fig. 2. Data-driven (in brown) and dynamical (in yellow-red) detection limits of TOI-1736 in the period-mass space, derived from SOPHIE RVs. Known planets are shown in black (triangle symbol if the mass is beyond the plotted range). Additional planets may only exist in the red-shaded area. A zoomed-in window around the inner planet illustrates nearby low-order MMRs (green: 6:5; gray: 7:6; blue: 8:7).

Strikingly, the outer CJ prevents any additional inner planet to exist beyond an orbital period of ∼150 days. This provides an upper limit to the inner system. Additionally, a small instability zone is observed in the neighbourhood of TOI-1736 b, also impacting the detection limits. The effect of nearby MMR overlap is visible, as highlighted by the vertical dashed lines. Overall, the inclusion of the stability constraint significantly strengthens the completeness of the system. Furthermore, this result is of high importance for RV follow-up programs. Indeed, by reducing the period window of potential additional planets, the observational strategy can be adapted consequently. An alternative view of these detection limits, in the orbital plane, is presented in Appendix D. Finally, we also applied ARDENT in the period range between 100 and 3000 days. While data-driven mass limits are saturated at periods larger than the RV baseline (∼900 days), the stability constraint informs the range of forbidden periods around the giant TOI-1736 c. We find that no additional planet below 1500 days can exist in the system.

4. Conclusions

The impact of CJs on SP formation and evolution remains an open question, placing the Solar System in the exoplanet context. Soon, data from Gaia DR4, PLATO, and subsequent RV follow-ups are likely to dramatically increase the number of SP-CJ systems. In this context, benefiting from the strongest systems completeness with the data at hand will be key to reveal CJ-SP (anti-)correlations. It is therefore timely to have a tool that computes the detection limits with the stability constraints. In this Letter, we introduce ARDENT, a publicly available Python package that derives both data-driven and dynamical detection limits. ARDENT is designed to be fast, thanks to a selective exploration of the period-mass space and the combination of analytical and numerical stability criteria. To test the method, we analysed TOI-1736, a system containing a 7-day sub-Neptune and a 570-day eccentric cold giant, from the publicly available SOPHIE RVs. This study highlighted the advantages of ARDENT. Within a limited amount of computation time, we were able to strengthen the system completeness and provide RV follow-up programs with a reduced period window for potential additional planets. The perspective of employing ARDENT on a sample of CJ systems is expected to improve the precision in SP occurrence rates, a crucial step to revealing SP-CJ correlations. While CJ systems are well-suited for ARDENT, the code can be used on a diverse range of system architectures, including compact systems or eccentric planets. The package is user-friendly, and can be automatised to derive detection limits on a large sample of systems².

Acknowledgments

We thank the referee for insightful feedback on the code and manuscript. M.S. acknowledges financial support from the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number C4000140754. L.N. acknowledges financial contribution from the INAF Large Grant 2023 “EXODEMO”. V.V.G. is an F.R.S.-FNRS Research Associate.

References

Appendix A: Sampling planets for the injection-recovery tests

To compute the data-driven detection limits, ARDENT proceeds to injection-recovery tests in the period versus RV semi-amplitude space. It generates a sample of injected planets by drawing P and K from log-uniform and uniform distributions, respectively. The range of K is expressed in units of rms (by default, [0.1; 1.2] rms inside which K is drawn uniformly). As such, it does not need adjustments from an instrument to another, but the user can modify this input range if desired.

Interestingly, the data-driven detection limits in K do not depend on the orbital period, but only on the timespan of the observations. In other words, the detection limit along K is a constant. This is illustrated in Fig. A.1, which presents the ARDENT output of the data-driven K detection limits for TOI-1736 analysed in this work. Each point in this period-K space corresponds to a set of ten different injected planets, each of them with a different orbital phase. The colour code depicts the detection rate as the number of detected planets over the ten different tests per (P, K). The detection limits stand at a constant value of K, until the RV timespan becomes an intrinsic limitation for the detectability.

Fig. A.1. TOI-1736 injection-recovery tests in the period versus K space. The detection limits (here the 50 and 95% limits are represented) are constant in K, until the RV timespan becomes too short compared to the injected periods. Consequently, a uniform sampling in K is efficient to sample the zone close to the detection limits.

Since detection limits are constant in the (P,K) space (considering long enough RV time spans), a uniform sampling in K is particularly efficient to sample the injected planet mass. Indeed, it insures that the sampled mass is always close to the desired detection limit. After converting the K distribution into planet mass, we obtain the plot shown in Fig. A.2 (which is also a direct output of ARDENT). This distribution of planets in the period-mass space samples densely the useful region surrounding the detection limits, as a direct consequence of the sampling in K. From this period-mass sample, the data-driven mass limits are computed by fitting the cumulative distribution function (cdf) of the detection rate as a function of injected mass, for each period bin. Because this sampling is dense around the transition between 0 and 100% detection rates, the cdf is also well sampled.

Fig. A.2. TOI-1736 injection-recovery tests in the period-mass space, obtained after converting the K distribution (cf. Fig. A.1) into mass. Because of the efficient sampling in K, this mass distribution is also efficient to sample the detection limits.

Appendix B: Analytical stability criteria

ARDENT employs analytical criteria to estimate the orbital stability of planetary systems. These enable us to identify long-term stable and short-term unstable systems. If a clear answer about stability cannot be provided, the code initiates numerical investigations. In this section we expand on the description of the analytical tests.

B.1. Orbit crossings

ARDENT checks for crossing orbits on all the planet pairs in the system via the simple expression

$$\begin{array}{c}\hfill a_i(1+e_i)<a_o(1-e_o),\end{array}$$(B.1)

where the subscripts i and o refer to the inner and outer planets of the pair, and a and e are the initial semi-major axis and orbital eccentricity, respectively. When this condition is violated, the orbits cross. Unless the planet pair is in mean-motion resonance (MMR), the orbit crossing will ultimately lead to planet-planet close encounters and rapid instability. However, inside a MMR, the orbital phases can be protected from close encounters in crossing orbits (the example of the Neptune - Pluto pair in the Solar System). Therefore, crossing orbits are indicative of rapid instability only if the pair is not in resonance.

To systematically check for MMR, we adopt the simple resonance criterion of Dai et al. (2024) based on the distance from the exact period commensurability. The latter is defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }\equiv \frac{P_{\mathit{out}}/P_{\mathit{in}}}{p/q}-1,\end{array}$$(B.2)

where P_in and P_out are the orbital period of the inner and outer planet of the pair, and p and q are small integers such that p − q gives the order of the resonance. Based on the Kepler sample of multi-planet systems, Dai et al. (2024) define a distance threshold below which the planet pair is classified as resonant. In particular, for first-order MMRs (q=p-1), -0.015 < Δ < 0.03. For second and third-order MMRs (q=p-2 and q=p-3), -0.015 < Δ < 0.015. Within ARDENT, this criterion is used to check for the main low-order resonances:

• Order 1 MMRs: 2:1, 3:2, 4:3, 5:4, 6:5

• Order 2 MMRs: 3:1, 5:3, 7:5

• Order 3 MMRs: 4:1, 5:2

As noted by Dai et al. (2024), a significant amount of planet pairs satisfying this criterion do not have librating resonant arguments (yet, they likely underwent resonant capture in the past). This criterion is therefore optimistic in considering planet pairs in resonance. In the case of ARDENT, it implies that these close-to-MMR planet pairs will not be automatically rejected from the orbit crossing criterion, and will go through a numerical investigation. Very likely, orbital instability will then rapidly set in these non-resonant crossing systems.

B.2. AMD framework

In case the system passed the orbits crossing test, another stability criterion is applied based on the AMD of the system (Laskar & Petit 2017). The AMD measures the difference between the angular momentum of the system, and the one it would have if all the orbits were circular and co-planar. Hence it informs about the dynamical excitation of the system. The AMD of a system is given by summing the AMDs of its planets in the heliocentric coordinates (Laskar & Petit 2017):

$$\begin{array}{c}\hfill AMD={\displaystyle \sum _{k=1}^n}{\mathrm{\Lambda }}_k(1-\sqrt{1-e_k^2}cos{i}_k),\end{array}$$(B.3)

where ${\mathrm{\Lambda }}_k=M_k\sqrt{G{M}_{\star }{a}_k}$, M_k and M_⋆ are the planet and stellar masses, G is the gravitational constant, a_k is the semi-major axis of planet k, e_k, and i_k are the orbital eccentricity and inclination, respectively.

Within this framework, a system is long-term stable if any redistribution of the angular momentum among the planets cannot lead to close encounters, compatible with the Hill stability (Petit et al. 2018). Therefore, it is possible to define a critical AMD C_c^H below which this condition is met and the system is Hill-AMD stable. This Hill-AMD criterion was investigated for each consecutive planet pair, and the system is AMD-stable if all the pairs are. Considering two planets with semi-major axes a₁ and a₂ such that $\alpha =\frac{a_1}{a_2}<1$, the critical Hill-AMD is given by (Petit et al. 2018)

$$\begin{array}{c}\hfill C_c^H=\gamma \sqrt{\alpha }+1-{(1+\gamma )}^{3/2}\sqrt{\frac{\alpha }{\gamma +\alpha }(1+\frac{3^{4/3}{ϵ}^{2/3}\gamma }{{(1+\gamma )}^2})},\end{array}$$(B.4)

where $\gamma =\frac{M_1}{M_2}$ and $ϵ=\frac{M_1+M_2}{M_{\star }}$. The condition of Hill stability is verified if for every consecutive planet pair

$$\begin{array}{c}\hfill {\beta }_H\equiv \frac{AMD/{\mathrm{\Lambda }}_2}{C_c^H}<1.\end{array}$$(B.5)

Alternatively, for close planet pairs, the resonance web becomes dense and a critical AMD C_c^MMR based on the overlap of first-order MMRs can be defined (Petit et al. 2017). Above this limit, MMR overlaps induce rapid instability. The critical AMD in the context of MMR overlaps is given by (Petit et al. 2017)

$$\begin{array}{c}\hfill C_c^{\mathit{MMR}}=\frac{\gamma }{1+\gamma }\frac{g^2}{2},\end{array}$$(B.6)

where $g=\frac{3^4{(1-\alpha )}^5}{2^{9}rϵ}-\frac{32rϵ}{9{(1-\alpha )}^2}$ and r = 0.80199. The planet pair is AMD-unstable if

$$\begin{array}{c}\hfill {\beta }_{\mathit{MMR}}\equiv \frac{AMD/{\mathrm{\Lambda }}_2}{C_c^{\mathit{MMR}}}>1.\end{array}$$(B.7)

Petit et al. (2017) derived the threshold semi-major axes ratio α_R above which the MMR overlap criterion is to be considered. Below this threshold, the Hill criterion prevails.

To synthesise, the AMD of consecutive planet pairs is computed and compared to C_c^H or C_c^MMR depending on the semi-major axes ratio α = a₁/a₂ of the pair:

• α < α_R: stable system if β_H< 1

• α > α_R: unstable system if β_MMR> 1

In conclusion, several analytical criteria are used in ARDENT to classify the stability of planetary systems. The systems are classified as

1. unstable if (i) orbits cross far from MMRs, or (ii) if the systems are located in the MMR overlap region (β_MMR> 1);

2. stable if they pass the Hill-AMD stability criterion (β_H< 1);

3. undefined otherwise. In this case, a numerical investigation is initiated.

Finally, it is worth noting that in systems containing an eccentric outer giant planet such as in TOI-1736, the majority of the injected planets are either unstable due to the orbits crossing criterion, or need further numerical investigation. Fewer injected planets are Hill-AMD stable, the reason being that any planet pair including the outer eccentric giant has a large AMD which penalises the Hill stability.

Appendix C: Calibration of the threshold in semi-major axis drift

ARDENT classifies the stability of planetary systems based on analytic criteria at first, and via numerical simulations if needed. The latter consist of short integrations over 10⁴ orbits of the outermost planet, after which the level of chaos is estimated. Chaos is indicative of orbital instability (e.g. Stalport et al. 2022) before that a close encounter or planet escape occurs. Therefore, it is an efficient way to speed up computation time. To estimate chaos in planetary systems, ARDENT measures the drift in semi-major axes. Under secular, non chaotic dynamics, the averaged semi-major axes are conserved quantities and do not vary on long time scales. In contrast, they drift significantly in chaotic dynamics. Therefore, the measure of the semi-major axes drift offers an estimation of chaos and, indirectly, orbital stability. As was shown by Tamayo et al. (2016), the study of semi-major axes provides a reliable stability indicator. They trained a machine on several orbital features based on short (10⁴ orbits) integrations to recognise patterns of instability, and compared those results with a control set derived from long-term (10⁷ orbits) integrations. The variation in semi-major axes was ranked among the top features predicting (in)stability.

In ARDENT, we define the relative drift in semi-major axis as

$$\begin{array}{c}\hfill \mathrm{\Delta }a=ma{x}_i\left(\frac{|a_{2,i}-a_{1,i}|}{a_{0,i}}\right),\end{array}$$(C.1)

with a_1, i and a_2, i as the median semi-major axis of planet i over the first and second halves of the integration, respectively, and a_0, i as the initial semi-major axis. The semi-major axis drift of the whole system is defined as the maximum relative drift amongst all planets. At the end of the numerical simulation, the planetary system is classified as unstable if either the integration did not reach the end (escape or close-encounter), or if one of the planetary orbits drifted by more than a certain threshold Δa > Δa_max.

This threshold drift Δa_max is an argument of the ARDENT’s function computing the dynamical detection limits, and hence it can be easily adjusted. However, we carried out experiments to calibrate this threshold and motivate its default value in ARDENT. We generated three families of 1000 artificial 3-planet systems each, composed of two inner SP and an outer CJ. In all systems, the period and mass of the innermost and outermost planets are fixed at [10, 500] days and [5 M_⊕, 5 M_Jup], respectively. Regarding the middle planet, its period and mass are randomly drawn from log-uniform distributions in [11, 200] days and [0.1, 10] ×M_b, respectively. The inner SPs have null eccentricities. In the first set of systems, the CJ has a null eccentricity e_CJ = 0 (the null eccentricity set). In the two other samples, e_CJ = 0.3 and 0.7 (the intermediate and high eccentricity sets, respectively).

For every planetary system, we carried out two numerical simulations. We ran a short simulation over 10⁴ × P_outer, similar to the ARDENT simulations time scale, after which Δa was computed. Secondly, we ran a long simulation over 10⁷ years for reference³. With an arbitrary drift threshold Δa_max, we could estimate each system stability from the short simulations and compare this classification to the real answer coming from the control sample of long integrations. Systems that were classified as unstable but survived the long integration were defined as false negatives (FN). Alternatively, systems wrongly classified as stable were defined as false positives (FP). We investigated the amounts of FN and FP for a various range of drift thresholds Δa_max, and present the result in Fig. C.1. As Δa_max increases, the quantity of FN decreases, while FP increases.

In ARDENT, we opted to prioritise minimising the FN over the FP. Indeed, FN are systems wrongly classified as unstable, leading to their rejection. The resulting risk is to over-estimate the dynamical constraints, excluding potential valid zones of the period-mass space. On the other side, a large amount of FP leads to poorer constraints on the dynamical detection limits, because many systems are wrongly kept after applying the stability constraint. While this needs to be avoided, it does not lead to wrong results as opposed to large FN numbers. Imposing a rate of FN smaller than 5%, and minimising the rate of FP accordingly, we converge to a drift threshold Δa_max = 0.025%. It is represented by the diamond symbol in Fig. C.1. At this stability threshold, the rate of FP is ∼50%. This is a conservative approach that favours the reliability of the results over their ability to refine the updated mass limits. We recommend this strategy, and define the default drift threshold in ARDENT as Δa_max = 0.025%. However, we emphasise that the Δa_max parameter in ARDENT can be easily modified by the user (via the argument max_drift_a in the function ARDENT_DetectionLimitStab). Importantly, the calibration experiments presented here assumed a certain system architecture, restricting their validity. The optimal stability threshold may differ for significantly different system architectures. The best approach consists in testing various drift thresholds or calibrating this drift again, when ARDENT is to be employed on systems with different architectures than the CJ-SP case presented here.

Fig. C.1. Calibration of the numerical stability threshold employed in ARDENT. This plot results from short and long numerical integrations on three sets of systems, differing from each other with the eccentricity of the outer CJ. The relative amounts of FP and FN are shown for all sets combined, and for a various range of threshold drifts in semi-major axis Δamax indicated by the colorbar. The diamond symbol corresponds to the default ARDENT value of Δamax = 0.025%. This conservative choice favours minimising the amount of FN.

Appendix D: ARDENT command lines for TOI-1736

This section illustrates ARDENT practically, and details the command lines used to produce the results presented in Sect. 3 on TOI-1736. The residual SOPHIE RVs, obtained after removing the best fit model, were derived prior to run this code. These residuals, together with an exhaustive overview of ARDENT functions, are available at https://github.com/manustalport/ardent. Noteworthy, the functions computing the data-driven and dynamical detection limits are distinct. As a result, it is possible to skip the computation of data-driven detection limits if they are already available, and to derive the dynamical constraints straight away. We refer to the link above for further details.

import ardent
import numpy as np

data = np.genfromtxt( \
    ’TOI1736_sophie_residualRV.dat’, \
    skip_header=int(2))
jdb, rv, rv_err = data[:,0:3].T

# Initialise ardent
vec = ardent.ARDENT_tableXY(jdb, rv, rv_err)

# Define an output folder that will contain the
# results of this ardent run
vec.ARDENT_Set_output_dir(’TOI-1736’)

# Initialise the star and planets
vec.ARDENT_AddStar(mass=1.08,starname=’TOI-1736’)
vec.ARDENT_AddPlanets(p=570.3, k=200.4, \
    mean_long=244.0, e=0.362, omega=164.3, \
    asc_node=0.0, inc=70.0)
vec.ARDENT_AddPlanets(p=7.065, k=4.34, \
    mean_long=247.2, e=0.21, elower=0.0, \	
    omega=0.0, asc_node=0.0, inc=90.0)

# Data-driven detection limits
vec.ARDENT_DetectionLimitRV(rangeP=[2., 600.])

# Dynamical detection limits
vec.ARDENT_DetectionLimitStab(NlocalCPU=4)

# Plot the detection limits
vec.ARDENT_Plot_StabDL(inset_plot=True, \
    zoom_xlim=[5.6,9.0], zoom_ylim=[0,16])

# Plot these limits in the orbital plane
vec.ARDENT_FinalPlot()

The last line above produces an alternative representation of the detection limits in the orbital plane. It is presented in Fig. D.1, where the planetary orbits are drawn on top of the detection limits. This is visually informative to measure the impact of the stability constraints on the detection limits, and sense the extent of the zone where potential unseen planets could exist (e.g. Brewer et al. 2020). The left plot shows the data-driven detection limits, to be compared with the plot on the right presenting the dynamical detection limits. These limits were computed in a period range between 2 and 600 days; hence, the coloured areas do not extend until the borders of the plots. The impact of the eccentric giant planet is very effective in pushing down the detection limits. Any follow-up RV survey must focus their search for additional planets in the inner part of the system only, within the red and blue zones of the plot on the right.

Fig. D.1. TOI-1736 detection limits presented in the orbital plane. The star and the two known planetary orbits are represented. The green shaded region is the predictive habitable zone, whose extent is computed from Kopparapu et al. (2013). The concentric discs illustrate the detection limits, with a colour code depicting the 95% limit mass. Left: Data-driven detection limits, based on the SOPHIE RVs only. Right: Dynamical detection limits, which add the constraint of orbital stability.

All Figures

Fig. 1. Schematics illustrating dynamical detection limits computation in ARDENT. Empty circles represent the data-driven detection limits, above which 𝒫d ∼ 1. ARDENT models the dynamical influence of known planets in the system (diamond symbol), inducing an instability zone where additional planets are excluded. The resulting dynamical detection limits are shown in red. The inner panel illustrates the iterative mass limits update accounting for stability, with dM the difference in injected mass between two consecutive steps and δ the precision threshold (cf. Sect. 2.3.1).

In the text

Fig. 2. Data-driven (in brown) and dynamical (in yellow-red) detection limits of TOI-1736 in the period-mass space, derived from SOPHIE RVs. Known planets are shown in black (triangle symbol if the mass is beyond the plotted range). Additional planets may only exist in the red-shaded area. A zoomed-in window around the inner planet illustrates nearby low-order MMRs (green: 6:5; gray: 7:6; blue: 8:7).

In the text

	Fig. A.1. TOI-1736 injection-recovery tests in the period versus K space. The detection limits (here the 50 and 95% limits are represented) are constant in K, until the RV timespan becomes too short compared to the injected periods. Consequently, a uniform sampling in K is efficient to sample the zone close to the detection limits.
In the text

	Fig. A.2. TOI-1736 injection-recovery tests in the period-mass space, obtained after converting the K distribution (cf. Fig. A.1) into mass. Because of the efficient sampling in K, this mass distribution is also efficient to sample the detection limits.
In the text

Fig. C.1. Calibration of the numerical stability threshold employed in ARDENT. This plot results from short and long numerical integrations on three sets of systems, differing from each other with the eccentricity of the outer CJ. The relative amounts of FP and FN are shown for all sets combined, and for a various range of threshold drifts in semi-major axis Δamax indicated by the colorbar. The diamond symbol corresponds to the default ARDENT value of Δamax = 0.025%. This conservative choice favours minimising the amount of FN.

In the text

Fig. D.1. TOI-1736 detection limits presented in the orbital plane. The star and the two known planetary orbits are represented. The green shaded region is the predictive habitable zone, whose extent is computed from Kopparapu et al. (2013). The concentric discs illustrate the detection limits, with a colour code depicting the 95% limit mass. Left: Data-driven detection limits, based on the SOPHIE RVs only. Right: Dynamical detection limits, which add the constraint of orbital stability.

In the text