© The Authors 2026

1 Introduction

Co-orbital objects, also known as Trojans, are celestial bodies that share a 1:1 mean-motion resonance with a planet. They typically reside near the system’s Lagrange points, L₄ and L₅, approximately 60° ahead of and behind the planet along its orbit. Depending on their orbital behavior relative to the planet, coorbital bodies can follow two main types of trajectories: tadpole orbits and horseshoe orbits (Garfinkel 1976; Erdi 1977; Garfinkel 1978; Niederman et al. 2020; Raymond et al. 2023). In a tadpole orbit, the object librates around either the L₄ or L₅ point with a small angular amplitude, forming a path that resembles a tadpole shape in a rotating reference frame. In contrast, objects on horseshoe orbits exhibit much larger oscillations in angular separation, traveling along a broad arc that encompasses the L₃, L₄, and L₅ points, effectively tracing a horseshoe-shaped path relative to the planet. Tadpole orbits are typically more robust than horseshoe orbits, because the horns of the horseshoe approach too close to the planet (in terms of the planet’s Hill radii) for more massive planets (Laughlin & Chambers 2002; Cuk et al. 2012). Early Trojan progenitors on horseshoe orbits can be captured into stable tadpole configurations during Jupiter’s gas accretion and mass growth (Fleming 2000; Raymond et al. 2023).

The existence of Trojan planets in exoplanetary systems, known as exotrojans, has significant implications for our understanding of planet formation and migration processes (Chiang & Lithwick 2005). Similar to their counterparts in our Solar System, exotrojans are thought to be by-products of early planetary evolution. Thus, studying these celestial bodies can offer valuable insights into the initial stages of planetary system development (Lillo-Box et al. 2018a). Several theories have been proposed to explain the formation of exotrojans. Two prominent scenarios exist: the in situ formation scenario, where exotrojans are remnants of early planet formation (Laughlin & Chambers 2002; Chiang & Lithwick 2005; Lykawka & Horner 2010), and the late-capture scenario, where they are captured during planetary migration (Laughlin & Chambers 2002; Chiang & Lithwick 2005; Beaugé et al. 2007; Lyra et al. 2009).

Various techniques have been developed to detect exotrojans. One prominent method, introduced by Ford & Gaudi (2006), combines transit and radial velocity (RV) observations to search for exotrojans and constrain the mass of a possible co-orbital body. Another powerful approach involves analyzing transit timing variations (TTVs), which are deviations from a planet’s strictly periodic transit schedule. The gravitational pull from a co-orbital companion perturbs the transiting planet’s orbit, causing its transits to occur slightly earlier or later than predicted. Ford & Holman (2007) first explored this effect for exotrojans, demonstrating that even a sub-Earth-mass companion could produce a TTV signal detectable by ground-based observatories. Subsequent studies, such as Haghighipour et al. (2013), reinforced this by showing that an Earth-mass exotrojan could induce TTVs with amplitudes of up to a few hours, well within the sensitivity of current ground-based telescopes.

Previous systematic searches for Trojans in Kepler data have been conducted by Janson (2013). To extend the search for coorbital companions such as exotrojans and exomoons, Leleu et al. (2017) introduced the α-test method. More recently, the RIVERS project (Leleu et al. 2021; Leleu et al. 2022) further advanced the detection of co-orbitals using TTVs. Building on these methods, the TROY Project (Lillo-Box et al. 2018a,b; Balsalobre-Ruza et al. 2024) systematically searches for exotrojan companions using a combination of transit, RV, and TTV techniques.

The TROY Project has systematically explored the population of Trojan bodies in exoplanetary systems. In its latest study, Balsalobre-Ruza et al. (2024) analyzed 95 transiting planets across 84 systems (focusing on low-mass stars) using the α-test method, setting upper mass limits for potential co-orbital companions and ruling out exotrojans heavier than Saturn. A strong candidate was identified in the GJ 3470 system, which marks the potential first detection of an exotrojan companion in an exoplanetary system. Furthermore, recent advancements in direct imaging have enabled tentatively probing co-orbital submillimeter emission in young systems (Balsalobre-Ruza et al. 2023; Long et al. 2022). These efforts demonstrate the growing capabilities in detecting and characterizing exotrojan companions.

The ExoEcho (Exoplanet Ephemerides Change Observations) project aims to study exoplanet photodynamics using precise transit timing data from ground- and space-based telescopes (Wang et al. 2023; Zhang et al. 2024b; Ma et al. 2025; Wang et al. 2025). As part of the ExoEcho project, this study investigates the possible existence of Trojan companions in hot-Jupiter systems using TTV observations from the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014). TESS, launched in 2018 as a successor to the Kepler space telescope (Borucki et al. 2010), has significantly advanced exoplanet research with its high-precision, all-sky photometric survey. Designed to expand our understanding of exoplanets through the transit method, TESS surveys nearly the entire sky, including the continuous viewing zones at the ecliptic poles. It yields a photometric precision of approximately 700 ppm for a 12th-magnitude star over a 1-hour integration (Ricker et al. 2014). These capabilities make TESS an ideal instrument for detecting exoplanets using transit techniques and exploring the broader landscape of exoplanetary science through photodynamics techniques. While hot Jupiters are not dynamically favorable for retaining primordial Trojans (Leleu et al. 2019), they serve as accessible targets for observational searches due to their high transit cadence.

This paper is structured as follows. Section 2 presents the TESS photometric observations of 260 hot-Jupiter systems used in our analysis. Section 3 details the method for deriving upper mass limits for potential exotrojans and presents our results. We discuss the implications of these findings in Sect. 4, and Sect. 5 summarizes the main conclusions of this study.

Fig. 1 Radius and period properties of the sample of 260 hot-Jupiter systems studied in this work. Most have orbital periods ranging from 2 to 6 days and radii between 0.8 and 1.6 RJ. The numbers of transits are shown in different colors.

2 Data and observation

The data used in this study were obtained from observations by TESS. We selected confirmed transiting hot Jupiters with published mass determinations based on the following criteria: an orbital period shorter than 10 days and a planetary mass exceeding 0.25 M_J. To ensure a robust sample, only systems with at least three observed transits were included. This selection yields a sample of 260 confirmed hot-Jupiter systems with RV-based mass measurements, representing one of the largest datasets used in similar studies to date. In Fig. 1, we show the distribution of the planetary orbital parameters, which indicates that most systems in our sample have orbital periods between 2 and 6 days and planetary radii ranging from 0.8 to 1.6 R_J.

2.1 Data acquisition

We used Astroquery (Ginsburg et al. 2017) to access TESS light curve products for hot Jupiters produced by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) from the Mikulski Archive for Space Telescopes (MAST; Brasseur & Donaldson 2019; MAST Team 2021). The analysis involved using pre-search data conditioning simple aperture photometry (PDCSAP) light curves from SPOC (Stumpe et al. 2012; Smith et al. 2012; Stumpe et al. 2014; Jenkins et al. 2016) with a 2-minute cadence, where common instrumental systematic errors were removed. Figure 1 also shows the distribution of the number of transits in the TESS observations used for our sample of 260 hot-Jupiter systems. Most systems have fewer than 100 transits, while some exceed this count. In particular, Qatar-10 b stands out, with more than 200 transits due to its short period (1.6 days) and its specific location in the sky (Alsubai et al. 2019).

Fig. 2 TTV measurements for WASP-93 b derived using a linear ephemeris model. The blue square represents the combined TTV obtained from TESS sectors 17 and 57. The insets show individual TTV measurements from each observation within these sectors, offering a zoomed-in view. The overall RMS for the TTVs is 32 seconds, with individual sectors showing an RMS of 31-32 seconds. Epoch 0 corresponds to the reference time of inferior conjunction at BJD 2458 765.64897.

2.2 Transit modeling and TTV extraction

We began by fitting the stitched light curve from all sectors of each hot Jupiter with the PyTransit package (Parviainen 2015), which employs a Markov chain Monte Carlo (MCMC) algorithm, to obtain preliminary estimates of the transit ephemerides. We fitted a Mandel & Agol (2002) model to the phase-folded light curve of each TESS sector. The parameters used for the transit model included the inferior conjunction T₀, orbital period P, stellar density ρ_* (in g cm⁻³), impact parameter b, and planet-to-star area ratio (R_p/R_*)². We assigned normal priors to T₀ and P, and uniform priors to ρ_*, b, and (R_p /R_*)². The eccentricity and argument of periastron were fitted using the parameterization of e cos w and e sin w. Quadratic limb-darkening coefficients were fitted using q₁ and q₂ parameterization (Kipping 2013) with uniform priors.

Most MCMC priors were derived from online catalog parameters (e.g., the means and standard deviations of the normal priors for period and T₀ match the values from the respective discovery papers). The prior for the mid-transit time of each event was centered on a generic linear ephemeris derived from literature values, allowing each transit to be modeled individually. We used the PyTransit implementation of the differentialevolution global optimization algorithm, followed by an MCMC sampling. We ran 100 chains with 2000 steps each. After the MCMC analysis, we checked the corner plot to ensure that the chains had converged.

Next, to precisely measure the mid-transit time for each individual transit, we segmented the TESS light curve into smaller windows, each spanning three times the transit duration and centered on the expected transit time. We then fitted the light curve within each window using PdotQuest (Wang et al. 2023). The TTVs were calculated by subtracting a linear ephemeris from the observed mid-transit times of each hot Jupiter. As an example, Fig. 2 shows the TTVs for WASP-93 b, derived from a linear model. The measured mid-transit times and TTVs for WASP-93 b are listed in Table 1 as an example for the machine-readable table available at the Centre de Données astronomiques de Strasbourg (CDS).

3 Analysis and results

3.1 TTV methodology

As demonstrated by Ford & Holman (2007), dynamical perturbations from an exotrojan librating around the L₄ and L₅ Lagrange points (as illustrated in Fig. 3) cause deviations in the strictly periodic transit times of the primary planet. These deviations manifest as a TTV signal with an approximate amplitude of $$A_{\mathrm{T}\mathrm{T}\mathrm{V}}\simeq 60\mathrm{s}\left(\frac{P}{4\mathrm{d}}\right)\left(\frac{M_{\mathrm{T}\mathrm{r}\mathrm{o}\mathrm{j}}}{M_{\oplus }}\right)\left(\frac{0.5M_{\mathrm{J}}}{M_{\mathrm{p}}+M_{\mathrm{T}\mathrm{r}\mathrm{o}\mathrm{j}}}\right)\left(\frac{K_{\mathrm{\Delta }M}}{10°}\right),$$(1)

where M_p and M_Troj denote the masses of the primary planet and the exotrojan companion, respectively, P is the orbital period of the primary, and K_ΔM represents the libration amplitude of the exotrojan. This method has inherent limitations, most notably that the TTV signal is directly coupled to the companion’s libration amplitude. Consequently, exotrojans with vanishingly small libration amplitudes - those residing at or very close to the triangular Lagrange points - produce negligible TTV signals and remain undetectable regardless of their mass. Furthermore, Eq. (1) reveals a physical degeneracy between the Trojan mass and libration amplitude, implying that our derived mass limits depend on the assumed libration.

Due to the symmetry of the restricted three-body problem in the circular approximation, Eq. (1) applies to both the L₄ and L₅ configurations. Although the phase of the TTV signal differs between a leading (L₄) and a trailing (L₅) companion, our methodology relies on the RMS of TTV residuals to establish upper mass limits. This approach is therefore insensitive to the specific location (L₄ vs. L₅) and constrains the presence of a co-orbital companion at either point capable of inducing TTVs above the noise floor. We also note that a symmetric population of Trojans at both L₄ and L₅ would dampen the net TTV signal, constituting another limitation of this technique.

We also consider the impact of calculating planet masses from RV observations. In systems containing an unaccounted massive Trojan, the observed RV amplitude reflects the combined mass of the planet and the Trojan (M_cat ≈ M_p + M_Troj)-As shown in Eq. (1), the TTV amplitude scales with the mass ratio M_Troj/(M_p + M_Troj). Since the catalog mass M_cat overestimates the true planet mass M_p, using it in the denominator leads to a reduction in the predicted TTV amplitude for a given Trojan mass. Consequently, a larger Trojan mass would be necessary to produce a detectable signal. Thus, our derived upper mass limits, which rely on these catalog masses, are actually conservative values for this upper limit.

Fig. 3 Schematic illustration of exotrojan orbits near a hot Jupiter. The figure shows both tadpole orbits, which librate around the Lagrangian points L4 and L5, and horseshoe orbits, which encompass a horseshoeshaped region from L4 to L5. The host star and hot Jupiter are indicated, and the rotating reference frame is centered on the star-hot Jupiter barycenter.

3.2 N-body simulations

We then performed N-body simulations using REBOUND (Rein & Liu 2012) to model a system that features an exotrojan companion and assess its ability to reproduce the observed data. The system’s initial conditions were derived from the best-fit parameters of the primary planet, assuming circular orbits for both the hot Jupiter and the Trojan. The exotrojan companion was placed near the 1:1 mean motion resonance (MMR) with the primary planet at the L₄ point (60° leading the planet). A further discussion of the impact of the circular orbit assumption is detailed in Sect. 4. To explore a range of libration amplitudes, we varied the initial osculating period ratio between the companion and the planet P_T/P_p from 0.9 to 1.1 in steps of 0.005. This initial period-ratio scan is an indirect way to sample the resulting libration amplitudes of near-co-orbital configurations, rather than directly specifying a physical period ratio away from 1:1 resonance. Some initial conditions in this scan correspond to chaotic, horseshoe-like, or unstable trajectories and were therefore excluded by the stability filtering described below. For each configuration, the libration amplitude K_ΔM was determined as the root mean square (RMS) of the relative mean longitude deviation from the Lagrangian point (σ_λtroj-_λplanet). We then calculated the corresponding TTV signals and derived mass constraints for each configuration.

We note that generating initial conditions based on the exact libration amplitude would be a more direct approach. For instance, varying the initial angular separation ζ from the Lagrange point (as adopted in e.g., Leleu et al. 2019) is equally straightforward to set up in REBOUND. Other alternatives, such as perturbing the initial velocity, also exist. We verified that all three approaches - varying the period ratio, the angular separation, or the velocity - produce equivalent TTV signals for a given libration amplitude. None of these methods guarantees an exact match to a desired libration amplitude due to the chaotic nature of N-body interactions, and our systematic scan of initial period ratios effectively covers the relevant range of libration amplitudes. A drawback of our simplified method is that some initial conditions with relatively large period deviations may correspond to chaotic or unstable orbits. To address this, we integrated a subsequent stability selection step (detailed in Sect. 3.3) to exclude any inherently unstable configurations from our final mass limit statistics. The duration of the simulated data was set to twice that of the TESS observations to avoid underestimating the RMS amplitudes. Figure 4 shows the simulated TTV amplitudes as a function of the exotrojan companion mass, compared to the observed data.

The RMS method has been widely used to analyze the TTV observation data of hot Jupiters in the past. For example, HAT-P-25 b is a hot Jupiter discovered by Quinn et al. (2011). Wang et al. (2018) utilized the RMS approach to establish the maximum mass limit of HAT-P-25 b’s companion and ruled out any hidden long-period perturber with M_pert > 0.5,0.3, and 0.5 M_⊕ near 1:2, 2:1, and 3:1 resonances. Figure 5 illustrates the RMS variations of simulated TTVs as a function of exotrojan mass for different libration amplitudes (K_ΔM). The intersections of these curves with the observed TTV RMS level (horizontal dashed line) define the upper mass limits corresponding to each libration amplitude. This demonstrates the degeneracy between the companion’s mass and its libration amplitude: larger libration amplitudes require lower masses to produce the same TTV signal. This mass sampling strategy aligns with RMS-based sensitivity analyses for companion mass constraints in previous TTV studies (Wang et al. 2018, 2021). The above procedure was performed using our newly developed tool, Companion Mass from TTV modeling (CMAT; Zhang et al. 2024a)¹, a Python package designed to automate the comparison between N-body simulations and observational data to constrain the mass of planetary companions (exotrojans in this case). The light curve modeling and dynamic simulation in this work closely follow the methods used in our previous paper (Zhang et al. 2024b).

To explicitly account for the uncertainties in the TTV measurements, we propagated the individual TTV errors to the final RMS calculation using standard uncertainty propagation. The uncertainty in the observed RMS (σ_RMS) was calculated as $${\sigma }_{\mathrm{R}\mathrm{M}\mathrm{S}}\approx \frac{1}{N\cdot \mathrm{R}\mathrm{M}{\mathrm{S}}_{\mathrm{o}\mathrm{b}\mathrm{s}}}{[\sum _{i=1}^N({\mathrm{T}\mathrm{T}\mathrm{V}}_i\cdot {\sigma }_i{)}^2]}^{1/2},$$(2)

where TTV_i and σ_i are the individual TTV measurements and their uncertainties and N is the number of data points. This uncertainty was then projected onto the mass-libration plane to determine the lower and upper bounds of the derived upper mass limit at the 1σ level. Additionally, we performed a χ² test to complement the RMS method. For each point in the parameter space (exotrojan mass vs. libration amplitude), we calculated the reduced χ² statistic (${\chi }_{\nu }^2$) comparing the observed TTVs with the simulated ones. This resulted in a χ² map with libration amplitude on the x-axis and exotrojan mass on the y-axis, where the color contours indicate the goodness of fit. This approach allowed us to determine the compatibility limits with statistical confidence, incorporating the observational errors.

Fig. 4 Observed TTV of the hot-Jupiter system WASP-93 b, compared with simulated TTV signals induced by co-orbital exotrojans of varying masses. Different curves represent exotrojan masses ranging from Earth-mass to sub-Neptune mass, illustrating the mass-dependent amplitude and periodicity of the induced TTV signals.

Fig. 5 RMS variations of simulated TTVs as a function of exotrojan mass for WASP-93 b, assuming different libration amplitudes. The colored curves represent different libration amplitudes (KΔM). The horizontal dashed red line indicates the observed TTV RMS. The intersections determine the upper mass limits for each specific libration amplitude. The 1σ uncertainty of the observed RMS, calculated via error propagation, determines the lower and upper bounds of the derived mass upper limit.

3.3 Stability analysis

The stability of a system with equilateral configuration - comprising a star and two planets of equal mass - was investigated by Gascheau (1843). The system is linearly stable for mass ratios satisfying $$\frac{M_{\mathrm{p}}+M_{\mathrm{T}\mathrm{r}\mathrm{o}\mathrm{j}}}{M_{\star }+M_{\mathrm{p}}+M_{\mathrm{T}\mathrm{r}\mathrm{o}\mathrm{j}}}<0.038,$$(3)

where M_p and M_Troj are the masses of the primary planet and exotrojan companion, respectively, and M_* is the mass of the host star. We used this criterion to ensure that all the orbital configurations studied are stable.

In addition to the analytical criterion, we performed numerical stability analyses using the mean exponential growth factor of nearby orbits (MEGNO) indicator (Cincotta et al. 2003). For each configuration, we integrated the system for 10000 orbits and calculated the MEGNO value. We excluded any configuration with a MEGNO value greater than 3, which indicates a chaotic or unstable orbit. While integrations spanning 10⁴ orbits (typically hundreds of years for short-period planets) are short compared to the system’s lifetime, for these close-in hot Jupiters, the libration periods around L₄ and L₅ are typically only a few dozen planetary orbits. This means 10⁴ planetary orbits already traverse hundreds of complete libration cycles, making them well suited for identifying the stability structure relevant to the mass-limit contours studied here. To assess how this baseline affects our derived mass limits, we repeated the MEGNO analysis for four representative systems (WASP-93 b, Wendelstein-1 b, WASP-20 b, and HAT-P-4 b) up to 10⁶ orbits. These representative validation cases for the present mass-limit workflow show that, although the longer integrations reclassify an additional 1522% of the marginal edge grid points as unstable, the resulting Trojan mass upper limits change by at most 1%, because the affected points lie near the outer boundary of the tadpole region. Extending the integrations to 10⁶ orbits would therefore increase the computational cost by a factor of ∼100 while providing only negligible improvement for the mass-limit contours relevant to this study. We therefore adopted 10⁴ orbits as a computationally efficient baseline for the present compatibility analysis, while noting that this comparison is not intended as a general statement about the long-term stability of all 260 systems.

We note that for some systems (particularly those with massive planets), stable configurations could not be established for libration amplitudes as large as 30°, meaning these specific configurations are intrinsically unstable and are therefore automatically excluded from the derived mass limit statistics at these high libration amplitudes. This is consistent with the restricted three-body problem, where the size of the stable region around the Lagrange points decreases as the mass ratio increases (Murray & Dermott 2000). Our final resulting upper mass limits thus specifically denote the boundaries beyond which TTVs are observationally incompatible, considering only those companion configurations that remain strictly dynamically stable. Figure 6 illustrates this approach for the representative system WASP-93 b. The MEGNO mask defines the dynamically stable region, while the 30° contour serves as a practical reference line within the tadpole zone used in our discussion of conservative mass limits.

Fig. 6 Stability and dynamics heatmap for WASP-93 b. The color scale and contours represent the libration amplitude corresponding to the initial period ratio and Trojan mass. The dashed red line marks the 30° libration contour, shown here as a practical reference line within the tadpole zone. Our conservative upper mass limits are strictly derived from the subset of configurations that remain dynamically stable according to the MEGNO mask.

3.4 Population results

Our analysis of 260 hot-Jupiter systems using CMAT indicates that the upper mass limits for exotrojan companions strongly depend on the libration amplitude. For a typical libration amplitude of 15° (close to the mean libration amplitude of Jupiter’s Trojans, ∼14°; Milani 1993; Murray & Dermott 2000), we find that, with the exception of the two massive super-Jupiters (HATS-70 b and HATS-41 b, both with M_p ≈ 10-13 M_j), all other studied systems have upper mass limits for Trojans strictly below 10 M_⊕ (with a median limit of 1 M_⊕). The higher mass limits for these super-Jupiters are primarily due to the extreme masses of the primary planets, which dampen the TTV effect of any Trojan. This represents a significant advancement in ruling out massive exotrojans in hot-Jupiter systems. We plot the cumulative distribution function (CDF) of the upper mass limit of potential exotrojans in Fig. 7. The left panel shows the results based on the RMS method, while the right panel displays the results from our χ² analysis. As demonstrated in the right panel of Fig. 7, while the χ²-based limits are more conservative than the RMS-based limits (left panel), they provide a statistically robust exclusion zone. Specifically, based on the RMS method, 50% of the hot-Jupiter systems in our study cannot support an exotrojan with a mass exceeding 1 M_⊕, and 80% cannot support one with 2 M_⊕. For the more conservative χ² analysis, limits are slightly higher, with 50 and 80% of systems excluding companions more massive than 3 M_⊕ and 6 M_⊕, respectively. Additionally, as discussed in Sect. 3, the use of catalog masses (which overestimate M_p) makes our derived limits conservative; using the true planet masses would shift the CDF curves in Fig. 7 to the left.

3.5 Correlations with planetary parameters

We also sought to identify correlations between the upper mass limits of Trojan companions and the hot Jupiters’ parameters shown in Fig. 8. The orbital period (P) and radius (R_p) of the planets are categorized into 15 × 6 bins, using a logarithmic scale for the orbital period. The color scheme represents the median upper mass limit in Earth mass for the exotrojans. We computed Pearson correlation coefficients between the upper mass limits of exotrojan companions and the radii and orbital periods of hot Jupiters. Our analysis does not reveal any significant correlations within the current dataset. This may suggest that the current precision of TTV measurements is primarily limited by observational noise rather than being sensitive to subtle physical correlations with planetary parameters. However, the absence of hot Jupiters in the top right corner bins of Fig. 8 suggests that the inclusion of future data in these bins could potentially impact the correlation outcomes.

4 Discussion

4.1 Comparison with previous work

In this section, we compare our findings here with those of the TROY Project (Lillo-Box et al. 2018a,b; Balsalobre-Ruza et al. 2024). Our study focuses on co-orbital objects associated with hot Jupiters, whereas TROY Project I and II (Lillo-Box et al. 2018b) concentrate on exoplanets with short orbital periods (P < 5 days), and TROY Project III (Balsalobre-Ruza et al. 2024) examines systems orbiting low-mass stars. As illustrated in Fig. 7, our results follow a trend similar to that in previous studies. However, it is important to note the methodological differences: while the TROY project models data directly with exotrojan models to perform a rigorous dynamical characterization, our study presents a compatibility analysis assessing whether TTVs are consistent with a hypothetical Trojan. Consequently, although our sample size is considerably larger (five times that of TROY project I, ten times that of TROY project II, and twice that of TROY project III), our analysis provides a constrained compatibility check rather than the full dynamical modeling of the TROY project. This larger sample size allows for a broad statistical overview, complementing the in-depth targeted studies. Specifically, 29 hot-Jupiter systems in our sample have also been studied in the TROY Project. As illustrated in Fig. 9, our results generally provide improved upper mass limits for exotrojan companions in these systems. This improvement is likely attributable to the high photometric precision and denser observational cadence of TESS, which lead to more precise midtransit time measurements compared to the ground-based data used in parts of the TROY project.

Furthermore, we note that for the comparison in Fig. 9, we adopted a typical libration amplitude of 15° (consistent with the mean libration of Jupiter Trojans; Milani 1993; Murray & Dermott 2000). These sensitivity limits are therefore biased toward this specific libration amplitude; as demonstrated in Sect. 3, a smaller libration amplitude would result in weaker mass constraints. This dependency on the dynamical configuration is a key limitation of our method compared to full dynamical fits.

In Fig. 9, we present the upper mass limits derived from our TTV analysis. The 1σ uncertainties correspond to the intersection of observed TTV RMS bounds with our simulation grids. For systems where the lower bound of the upper mass limit is unconstrained (consistent with zero mass at the 1σ level), we indicate this using vertical lines extending to the bottom of the axis. While most analyses assume a libration amplitude of 15°, we adopted 14.5° for HAT-P-20 b to ensure orbital stability. For improved visual clarity, error bars smaller than 15% of the mass limit were omitted.

Fig. 7 Cumulative probability distribution of the upper mass limit for potential exotrojans in 260 hot-Jupiter systems. Left: results from the RMS method, compared with TROY project results in red. Right: results from our χ2 analysis, offering a more conservative but statistically robust upper limit. Reference masses for Neptune, Saturn, and Jupiter are shown. Note that systems with unstable orbital configurations at the considered libration amplitudes (e.g., typically around 30° for massive planets) are excluded from the statistics. The shaded region represents the 1σ uncertainty for a libration amplitude of 15°.

Fig. 8 Median upper mass limits of potential exotrojans in 260 hot-Jupiter systems, binned by orbital period (1-30 days) and planetary radius (0.5-2.2 RJup). The colors represent the median mass limit in each bin.

4.2 Degeneracy between mass and libration amplitude

A key constraint of the TTV method is the degeneracy between the exotrojan mass and the libration amplitude (Eq. (1)). We addressed this by scanning a range of period ratios to simulate various libration amplitudes. Consequently, the mass limits presented in this work are explicitly defined as a function of the libration amplitude. While a smaller true libration amplitude would indeed enable a more massive companion to remain undetected, our results provide a comprehensive map of these compatibility limits across the dynamical parameter space.

4.3 Impact of eccentricity

We acknowledge that assuming circular orbits is a simplification, as hot Jupiters often possess nonzero eccentricities. While our bulk analysis maintains this assumption due to the typically low eccentricities of these planets, we recognize that eccentricity can introduce systematic effects. To quantify this, we performed additional simulations using measured eccentricities for a subset of 82 systems with e > 0.05 (Fig. 10). We find that including eccentricity generally tightens the derived mass limits while preserving the overall trend. Interestingly, the impact of eccentricity on TTV amplitude is nonlinear with respect to libration amplitude, with the effect being minimized at ∼15°. This minimum coincides remarkably well with the mean libration amplitude of Jupiter’s Trojans (∼14° ; Milani 1993), although the physical significance of this agreement remains unclear. Finally, we note that for systems with moderate eccentricities, the stable point for a “quiet” configuration (small libration) shifts away from the exact 60° separation, a dynamical feature that must be accounted for in accurate definitions of libration amplitude (Namouni 1999).

4.4 Chi-square analysis

Additionally, we incorporated a χ² analysis to complement the RMS method. While the RMS approach provides stringent constraints by globally aggregating TTV scatter, it treats all data points equally, irrespective of their individual observational uncertainties. In contrast, the χ² analysis explicitly weights deviations by their measurement errors, providing a more statistically robust assessment of compatibility. This method prevents data points with high scatter and large uncertainties - which are common in ground-based follow-up observations or fainter TESS targets - from artificially inflating the apparent TTV signal. Our results show that the χ²-based mass limits are generally more conservative (i.e., higher) than those derived from the RMS method. For instance, while the RMS method excludes companions more massive than 1 M_⊕ for 50% of the systems, the χ² analysis places this threshold at 3 M_⊕. This difference underscores that while the RMS metric serves as a useful proxy for overall variability, the χ² statistic provides a more robust exclusion criterion by incorporating data quality. Despite these quantitative differences, the consistency between the overall trends of the RMS and χ² results strengthens the validity of our general conclusions concerning the scarcity of massive exotrojans in these systems.

Fig. 9 Upper mass limits of planets from TROY project I, II, III, and this work (ExoEcho). The x-axis displays planet names, and the y-axis presents mass limits on a logarithmic scale. The colors and markers distinguish each study. The red triangles represent our results assuming a libration amplitude of ∼15°. The error bars indicate the 1σ uncertainty, which is omitted when smaller than 15% of the mass limit for visual clarity. For systems where the lower 1σ boundary is unconstrained (i.e., consistent with zero mass), we represent the limit using a vertical line extending to the bottom of the axis.

Fig. 10 Cumulative probability distribution of the upper mass limit for a subset of systems with eccentricity e > 0.05. The dotted gray line represents the results assuming circular orbits, while the solid colored line shows the results using the measured eccentricities. The inclusion of eccentricity generally tightens the mass constraints. Similar to Fig. 7, systems that are dynamically unstable at the tested libration amplitudes are excluded.

4.5 Future prospects

Future work will expand this study by including less massive exoplanets in our sample. Since the TTV amplitude is related to the mass ratio between the exotrojan and the target planet (Eq. (1)), we expect to have much stronger constraints on their exotrojans for less massive target planets.

5 Conclusion

As part of the ExoEcho project, this study presents a comprehensive analysis of the compatibility limits of exotrojans in 260 hot-Jupiter systems using the TTV technique with the latest data from the TESS mission. Combining TTV data from the target planet with N-body simulations, we put constraints on the mass of exotrojans in 260 hot-Jupiter systems, mostly below 10 M_⊕, while assuming a libration amplitude of 15°. Thanks to the all-sky survey capability of TESS, our study has one of the largest sample sizes thus far among studies focusing on exotrojan detection. Our constraints on exotrojan masses in these systems, reinforced by χ² analysis, provide new statistical limits on co-orbital dynamics in the hot-Jupiter regime.

The libration of exotrojans around the Lagrangian points induces the TTVs of the hot Jupiters, which are typically nonperiodic and thus more challenging to detect in photometric observations (Janson 2013). Additionally, an inclined orbit can cause these exotrojans to avoid detection during transits, and their low masses make them difficult to identify through RV methods. Continued transit observations from TESS and future space missions such as PLATO (Rauer et al. 2024) and ET (Ge et al. 2024) will enable more sensitive mass limits, further constraining theories of exotrojan formation. Finally, it is crucial to reiterate the limitations of the TTV technique in this context. Our derived upper mass limits are sensitive to the dynamical state of the potential companion; specifically, our method is blind to exotrojans with negligible libration amplitudes, as they do not generate detectable TTV signatures. Consequently, the absence of a TTV signal excludes only those massive companions that are undergoing significant libration.

Data availability

The full machine-readable version of Table 1, containing the measured mid-transit times for all 260 systems, is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/709/A35.

Acknowledgements

We express our gratitude to the TESS mission for its outstanding contribution to exoplanet science. The precise transit-time measurements from TESS were instrumental in our study of TTVs in hot Jupiters. We thank NSFC grants 12073092 and 12103098, the Natural Science Foundation of Hainan Province (Grant Nos. 125RC632 and 424QN217), the National Key R&D Program of China (2020YFC2201400), the science research grants from the China Manned Space Project (No. CMS-CSST2021-B09), and the Earth 2.0 Space Mission program from SHAO for their financial support. Additionally, we appreciate the valuable assistance of the NASA Exoplanet Archive, which granted us access to a vast amount of observational data and resources that greatly aided our research efforts. We also thank open-source software developers for their contributions to the scientific community.

References

All Tables

All Figures

	Fig. 1 Radius and period properties of the sample of 260 hot-Jupiter systems studied in this work. Most have orbital periods ranging from 2 to 6 days and radii between 0.8 and 1.6 RJ. The numbers of transits are shown in different colors.
In the text

Fig. 2 TTV measurements for WASP-93 b derived using a linear ephemeris model. The blue square represents the combined TTV obtained from TESS sectors 17 and 57. The insets show individual TTV measurements from each observation within these sectors, offering a zoomed-in view. The overall RMS for the TTVs is 32 seconds, with individual sectors showing an RMS of 31-32 seconds. Epoch 0 corresponds to the reference time of inferior conjunction at BJD 2458 765.64897.

In the text

	Fig. 3 Schematic illustration of exotrojan orbits near a hot Jupiter. The figure shows both tadpole orbits, which librate around the Lagrangian points L4 and L5, and horseshoe orbits, which encompass a horseshoeshaped region from L4 to L5. The host star and hot Jupiter are indicated, and the rotating reference frame is centered on the star-hot Jupiter barycenter.
In the text

	Fig. 4 Observed TTV of the hot-Jupiter system WASP-93 b, compared with simulated TTV signals induced by co-orbital exotrojans of varying masses. Different curves represent exotrojan masses ranging from Earth-mass to sub-Neptune mass, illustrating the mass-dependent amplitude and periodicity of the induced TTV signals.
In the text

Fig. 5 RMS variations of simulated TTVs as a function of exotrojan mass for WASP-93 b, assuming different libration amplitudes. The colored curves represent different libration amplitudes (KΔM). The horizontal dashed red line indicates the observed TTV RMS. The intersections determine the upper mass limits for each specific libration amplitude. The 1σ uncertainty of the observed RMS, calculated via error propagation, determines the lower and upper bounds of the derived mass upper limit.

In the text

Fig. 6 Stability and dynamics heatmap for WASP-93 b. The color scale and contours represent the libration amplitude corresponding to the initial period ratio and Trojan mass. The dashed red line marks the 30° libration contour, shown here as a practical reference line within the tadpole zone. Our conservative upper mass limits are strictly derived from the subset of configurations that remain dynamically stable according to the MEGNO mask.

In the text

Fig. 7 Cumulative probability distribution of the upper mass limit for potential exotrojans in 260 hot-Jupiter systems. Left: results from the RMS method, compared with TROY project results in red. Right: results from our χ2 analysis, offering a more conservative but statistically robust upper limit. Reference masses for Neptune, Saturn, and Jupiter are shown. Note that systems with unstable orbital configurations at the considered libration amplitudes (e.g., typically around 30° for massive planets) are excluded from the statistics. The shaded region represents the 1σ uncertainty for a libration amplitude of 15°.

In the text

	Fig. 8 Median upper mass limits of potential exotrojans in 260 hot-Jupiter systems, binned by orbital period (1-30 days) and planetary radius (0.5-2.2 RJup). The colors represent the median mass limit in each bin.
In the text

Fig. 9 Upper mass limits of planets from TROY project I, II, III, and this work (ExoEcho). The x-axis displays planet names, and the y-axis presents mass limits on a logarithmic scale. The colors and markers distinguish each study. The red triangles represent our results assuming a libration amplitude of ∼15°. The error bars indicate the 1σ uncertainty, which is omitted when smaller than 15% of the mass limit for visual clarity. For systems where the lower 1σ boundary is unconstrained (i.e., consistent with zero mass), we represent the limit using a vertical line extending to the bottom of the axis.

In the text

Fig. 10 Cumulative probability distribution of the upper mass limit for a subset of systems with eccentricity e > 0.05. The dotted gray line represents the results assuming circular orbits, while the solid colored line shows the results using the measured eccentricities. The inclusion of eccentricity generally tightens the mass constraints. Similar to Fig. 7, systems that are dynamically unstable at the tested libration amplitudes are excluded.

In the text