© The Authors 2026

1 Introduction

With the number of confirmed detections growing steadily, the field of exoplanet research is in the process of branching out into a diverse set of subdisciplines, ranging from formation modelling to the atmospheric characterisation of individual objects. Another such ongoing pursuit is the search for low-mass companions to exoplanets, bodies conventionally called exomoons¹. With the exception of Mercury and Venus, every planet (and even some dwarf planets) in our Solar System harbours at least one moon. It is thus plausible that planets orbiting other stars are accompanied by their own satellites. A successful detection of such an object in an exoplanetary system would be the first step towards assessing the prevalence of exomoons around different types of host planets and building a population-level moon sample, which can help calibrate our models of moon formation, whether via capture (Hansen 2019), impact (Barr & Bruck Syal 2017) or in situ within the circumplanetary disc (e.g. Canup & Ward 2006; Cilibrasi et al. 2021; see Cugno et al. 2024 for the detection of a disc surrounding GQ Lup B). Such knowledge will provide vital clues that can inform formation and evolutionary models (e.g. Peale & Canup 2015; Batygin & Adams 2025), the selection of follow-up targets deemed potentially habitable (Martínez-Rodríguez et al. 2019; Dobos et al. 2022) and the interpretation of putative biosignature detections in the future (Rein et al. 2014).

The different detection methods devised and pursued since the search for exomoons began in the wake of the first exoplanet detections (see e.g. Sartoretti & Schneider 1999) include searching for anomalous features in planetary transit observations (e.g. Szabó et al. 2006; Simon et al. 2007; Kipping 2009), identifying the transit signatures of moons themselves as they eclipse free-floating planets (e.g. Limbach et al. 2021, 2024), measuring planetary radial velocity (RV) modulations (e.g. Ruffio et al. 2023; Vanderburg & Rodriguez 2021), monitoring for microlensing events (e.g. Liebig & Wambsganss 2010), probing for spectroastrometric signals (e.g. Agol et al. 2015; van Woerkom & Kleisioti 2024), attempting to image them directly (e.g. Lazzoni et al. 2020), and hunting for so-called Exo-Io signatures around hot Jupiters (Oza et al. 2019). So far, however, no bona fide detection has been made (candidates include Kepler-1625 b-i and Kepler-1708 b-i; Teachey & Kipping 2018; Heller et al. 2019; Kreidberg et al. 2019; Teachey et al. 2020; Kipping et al. 2022; Heller & Hippke 2024; see also Kenworthy et al. 2023).

More recent studies simulating the survival rate of moons around planets at different orbital separations from their hosts are suggestive of why there has not yet been a confirmed detection. The sensitivity of transit studies – inherently biased to short orbital periods and hence small separations between host and planet – appears to be limited to a parameter space where moons are prone to being ejected or tidally disrupted (see e.g. Dobos et al. 2021). To put it another way, exomoons are unlikely to survive in orbit around the planets accessible to transit observations. On the other hand, these simulations suggest that moons around far-out, long period exoplanets of the kind accessible to direct imaging studies have a high likelihood of remaining in orbit for several gigayears. Here, however, our current contrast and resolution capabilities limit the achievable sensitivity of direct detection attempts and other methods such as spectroastrometry or moon transit monitoring are still to yield a definitive detection.

A ‘middle ground’ between the two extremes is to hunt for moons around planets orbiting their host stars every few hundred days. In this window of opportunity, we expect the survival rate to be large enough to make a discovery possible, while also possessing the means of performing such a detection in the first place. As we shall show in this paper, precisely monitoring a given planet’s position over time can reveal deviations from the expected Keplerian trajectory, that is a higher-order epicyclic or ‘wobble’ motion, induced by the presence of an orbiting moon. Such a detection would amount to applying the astrometric method of planet detection by monitoring the periodic movement of a star to a planet–moon pair. This technique represents a readily applicable strategy towards making the first extrasolar satellite detection a reality. The convenient possibility of constraining the mass and orbital elements of the moon from the astrometric time-series renders this method especially compelling.

At the heart of this new drive towards dynamic moon detection lies the unprecedented astrometric precision of the near-infrared interferometer GRAVITY at ESO’s Very Large Telescope Interferometer (VLTI; GRAVITY Collaboration 2017). With an astrometric precision of 50 μas (GRAVITY Collaboration 2019), the instrument has already facilitated remarkable constraints on the dynamical masses and orbital solutions of a large sample of directly detected planets. Notably, GRAVITY enabled an independent dynamical detection and characterisation of β Pic c based solely on its perturbations to the orbit of β Pic b (Lacour et al. 2021) and provided the key data set used to resolve the first known brown dwarf, Gliese 229 B, into two near equal-mass binary components (Xuan et al. 2024; Whitebook et al. 2024).

Here, we present an analytical modelling effort that aims to determine our current astrometric exomoon detection sensitivity and explore how future interferometric instruments and facilities might aid the hunt for these elusive objects. This paper is structured as follows: in Sect. 2 we outline how to calculate the astrometric signal amplitude of a moon and to simulate its variation as a function of time. Sect. 3 deals with current and future capabilities of astrometric exomoon detection on the basis of several exemplary systems and the case of a putative moon around β Pic b in particular. The modelling of detection attempts based on multiple astrometric epochs is described in Sect. 4 before the derived sensitivity curves are presented in Sect. 5. We conclude our study in Sect. 6.

2 The astrometric signal

The concept of an astrometric exomoon detection essentially amounts to a re-application of the astrometric method of planet detection. As we shall see, the main differences and additional challenges with respect to the planet detection science case are the circumstance that the observable perturbation amplitudes are orders of magnitude smaller and the fact that we are required to perform the position measurements on the planets themselves. To aid future target selection and observation planning, it is therefore crucial to understand the behaviour of the observable as a function of different parameters. Here, we briefly introduce the astrometric signal amplitude and variability in time.

2.1 Astrometric signal amplitude

To assess the detectability of a given star–planet–moon configuration, we are interested in the planet’s deviation from a Keplerian orbit around the host as a consequence of being orbited by a moon. Projected onto the plane of the sky, one can picture this perturbation as a time-dependent offset from the expected planetary position. While varying periodically in right ascension (RA) and declination (Dec) according to the orbital parameters of the planet–moon system, we are initially mostly interested in the maximum radial deviation, a quantity we shall call the astrometric signal amplitude A and which we can estimate via

$$A={\left(\frac{G}{4{\pi }^2}\right)}^{1/3}{\left(\frac{P}{M_{\mathrm{p}\mathrm{l}}+M_{\mathrm{m}}}\right)}^{2/3}\frac{M_{\mathrm{m}}}{d}=\frac{a_{\mathrm{m}}}{d}\frac{M_{\mathrm{m}}}{M_{\mathrm{p}\mathrm{l}}+M_{\mathrm{m}}},$$(1)

where G is the gravitational constant, P is the orbital period of the planet–moon orbit, d is the distance of the target system, and M_pl and M_m are the masses of the planet and moon, respectively (Sahlmann et al. 2013; Quirrenbach 2010). For the second equality, we used Kepler’s third law to substitute P with an expression that includes the semi-major axis of the moon, a_m. In general, Equation (1) is an approximation that holds only for circular moon orbits.

To reiterate, the quantity A is an amplitude and is therefore only a snapshot of the deviant wobble movement exhibited by the planet. In this sense, it is comparable to the RV semi-amplitude, K. This analogy to the RV method for exoplanet detection extends even further: when hunting for planets using the RV method, we do not rely on a single measurement. To sample a large fraction of the phase-folded RV curve, many observations often spanning several orbital periods are required. This ensures two things: (1) ideally, the full radial velocity amplitude is registered and (2) the large number of epochs increases the detection significance. The same is true for the astrometric method of exomoon detection: while we do not know at which moment in time the planet shows the largest deviation and the measurement precision of our current instrumentation may not suffice to make a significant detection based on only a small number of observations, an extensive time-series of multiple astrometric epochs will solve both of these problems. We therefore require not just an estimate of the perturbation amplitude as provided by Equation (1), but also its behaviour in time. In the following, we shall briefly explain how this modulation can be computed analytically.

2.2 Simulating the signal variation in time

To investigate the behaviour of the orbital perturbation experienced by the planet as a function of different parameters and thus its detectability with current and future instruments in more detail, we need to analytically describe the time variation of the signal. To this end, we defined two models.

2.2.1 The star–planet model

First, we implemented an analytical two-body orbit model that is capable of integrating the trajectories of the host star and planet around their mutual centre of mass as a function of nine orbital parameters, namely the system’s parallax, ϖ, the stellar and planetary masses, M_s and M_pl, respectively, the semi-major axis of the planet, a_pl, its inclination, i_pl, eccentricity, e_pl, longitude of ascending node, Ω_pl, argument of periastron, ω_pl, and time of periastron passage, t_{peri, pl}. In broad terms, this is an algorithm that converts a given time of observation to a mean anomaly, solves Kepler’s equation iteratively, performs a series of three-dimensional space rotations and scales the resulting stellar position by the planet–star mass ratio to obtain the planetary position. Performed for a series of time steps, this procedure yields the orbital trajectory of the two bodies as projected onto the plane of the sky. Clearly, this model forms the baseline case of our investigation: a system devoid of a moon exhibiting perfect Keplerian two-body orbits.

2.2.2 The star–planet–moon model

On top of the star–planet model, we needed to simulate the system’s behaviour under the influence of a moon in orbit around the planet. To this end, we employed the same framework as above but extended it by seven additional orbital parameters: the moon mass, M_m, semi-major axis, a_m, inclination, i_m, eccentricity, e_m, longitude of ascending node, Ω_m, argument of periplanet, ω_m, and time of periplanet passage, t_{peri, m}. Thus, the star–planet–moon model incorporates 16 free parameters.

Whereas in the star–planet model we assumed the dynamical mass of the companion object or system to be entirely locked in the planet, we are here supposing that it is actually the sum of the planet and moon masses. Thus, the trajectory obtained from the star–planet model can be reused as the trajectory of the planet–moon barycentre. Superimposing the trajectories of the planet and moon in their own reference frame onto this barycentric trajectory yields the space motions of the planet and moon around the barycentre of the entire system. While this nested strategy avoids cumbersome three-body numerical integrations, it neglects higher-order gravitational interaction terms between star, planet, and moon. Such a simplification appears reasonable, however, given that we intend to apply this model to stable systems exhibiting wide planetary orbits and for only negligible timescales compared to long-term orbital disruption and ejection processes. As Kipping (2010) points out, this approach is justified provided the considered semi-major axis of the moon is smaller than 0.531 times the planet’s Hill radius, R_{Hill, pl}. The latter is defined as

$$R_{\mathrm{H}\mathrm{i}\mathrm{l}\mathrm{l},\mathrm{p}\mathrm{l}}=a_{\mathrm{p}\mathrm{l}}{\left(\frac{M_{\mathrm{p}\mathrm{l}}}{3M_{\mathrm{s}}}\right)}^{1/3}.$$(2)

2.2.3 Extracting the wobble

Comparing the trajectories simulated using the star–planet and the star–planet–moon models yields the orbital deviation resulting from the presence of a moon. This perturbation can be visualised by plotting the orbital trajectories of the planet for both cases. The left hand panel of Fig. 1 shows the resulting planetary wobble for a fictitious example system. The shape of the wobble depends on the system configuration and can be deconstructed into RA and Dec components presented in the right hand panels of Fig. 1. Unlike the simplified description of the astrometric signal amplitude given in Equation (1), the star–planet and star–planet–moon models that underpin the parametrisation of the visualised wobble are capable of handling eccentric planet and moon orbits.

3 Expected signals in the context of measurement uncertainties

While, in general, an exomoon cannot be detected by a single astrometric measurement, it can be instructive to compare to our current astrometric measurement capabilities the theoretical perturbation in planet position as a consequence of an orbiting body. To illustrate this, Table 1 lists multiple known exoplanets and the astrometric signal amplitudes they would present for a series of assumptive companion moons. Here, we have varied both the mass and semi-major axis of the moon. While the former assumes three distinct values for each case, for the reasons outlined in Sect. 2.2.2, the latter is coupled to the respective Hill radius of the planet, R_{Hill, pl}. Apart from ensuring the applicability of the nested two-body framework we employed as an approximation of the three-body situation we are in reality confronted with, restricting the semi-major axis of the moon to values below a constant critical fraction, χ_crit, of the respective Hill radius also safeguards against modelling system configurations that are dynamically unstable. Amongst other things, χ_crit depends on the eccentricity and obliquity of the planetary orbit (e.g. Holman & Wiegert 1999). Different estimates of the critical distance fraction include χ_crit ≈ 0.36 and χ_crit ≈ 0.49 for prograde (Barnes & O’Brien 2002; Domingos et al. 2006) and χ_crit ≈ 0.93 for retrograde orbits (Domingos et al. 2006). To meet both the two-body approximation and the dynamical stability conditions, we therefore chose the two representative values of χ = 0.5 and 0.1 to compute the case-dependent semi-major axes according to a_m = χR_Hill,pl.

The demonstrated and predicted astrometric measurement uncertainties associated with a set of current and future instruments listed in the bottom of Table 1 place the computed signal amplitudes for the different system configurations into context. At present, GRAVITY is the only instrument capable of robustly resolving the perturbations predicted for a large subset of cases presented in Table 1. Classical imaging instruments, on the other hand, cannot be expected to be capable of detecting any but the most extreme cases, where the moon-to-planet mass ratios are unlikely large.

Clearly, the system configurations envisaged in Table 1 are entirely imaginative and serve to illustrate how the signal strength changes for a set of exemplary cases. If one were to attempt an astrometric exomoon detection in earnest, however, one would ideally want to avoid having to resort to a blind search around an arbitrary planet. This is due to the prolonged period of time throughout which astrometric epochs need to be collected for in order to achieve a positive detection. Instead, an additional line of evidence that hints towards the potential existence of an orbiting body should be the starting point of an efficient astrometric moon hunt. In this regard, the case of β Pic b warrants further attention. Recently, Poon et al. (2024) have argued that the planet’s potentially increased obliquity (to be confirmed by future observations) could be explained by the perturbing influence of an exomoon. They also give estimates as to the moon mass and semi-major axis required to induce the observed planetary obliquity. Hereafter, we use the shorthand ‘P24-moon’ to refer to objects that comply with these constraints. Prior to the suggestion of the P24-moon, Macias et al. (2024) have already ruled out moons more massive than 0.25 M_Jup.

To understand whether our current instrumentation is capable of picking up on the astrometric signal of such a P24-moon, we can compute the expected amplitude in the vicinity of the moon mass and semi-major axis constraints using the models defined in Sect. 2.2. This region of the parameter space is visualised in the right panel of Fig. 2. It demonstrates that a P 24-moon would induce an astrometric wobble amplitude of 10 to 20 μas. Despite its unprecedented astrometric precision of 50 μas, such a perturbation would likely be difficult to detect using GRAVITY. However, plans are already underway to improve the astrometric accuracy of the existing GRAVITY instrument by addressing certain limitations (e.g. optical aberrations in GRAVITY Collaboration 2021), as well as by halving the operational wavelength – observing in the Y- and J-bands instead of the K-band. These efforts are part of the VLTI/PLANETES technological development programme at ESO². A future instrument developed under this programme would be expected to achieve a precision of 10 μas (as predicted in Lacour et al. 2014), providing the capability to search for a P 24-moon around β Pic b. Comparability of the expected signal amplitude with the instrumental measurement uncertainty implies a single-epoch signal-to-noise comparable to 1. In Sect. 4, we shall see how the detection significance increases when acquiring multiple epochs.

The unique potential inherent to the astrometric method is reinforced by the comparison to the detection capabilities of other methods. For one thing, strategies that hinge on observations of the planet transiting the stellar disc are not applicable to this non-transiting planet. What is more, Poon et al. (2024) point out that transits of the putative moon in front of the planet are unlikely. An alternative avenue towards eventual detection might lie in a comprehensive RV monitoring programme of the exoplanet. To gauge whether such a strategy is feasible for detecting a P24-moon, we visualised the expected RV semi-amplitude, K, as a function of moon mass and semi-major axis in the left hand panel of Fig. 2. As expected, the sensitivity gradient of the RV method is opposite to that of the astrometric technique: while both signals increase with the moon mass, the astrometric signal is stronger for larger semi-major axes, whereas the RV semi-amplitude decreases with larger separations. To probe the entire parameter space region predicted for the P24-moon requires an RV precision of approximately 75 m/s. Recent studies have already demonstrated measurements of planetary RVs to be possible when coupling an adaptive optics system with a high-resolution spectrograph. Denis et al. (2025) could constrain the RV of AF Lep b to within approximately 1 km/s using VLT/HiRISE (Vigan et al. 2024), while Parker et al. (2024) and Landman et al. (2024) achieved measurements of β Pic b with uncertainties of approximately 2 km/s and 300 m/s, respectively. Using KPIC (Delorme et al. 2021), Horstman et al. (2024) was able to obtain RV epochs of GQ Lup B with precisions ranging from 400 to 1000 m/s. Notwithstanding these promising results that will stimulate further development of the involved techniques and instruments, the attainable planetary RV precision does not yet suffice to detect the putative P24-moon.

Regardless of whether one pursues an RV or astrometry-based strategy, a single epoch will never suffice for detection. Instead, extensive monitoring of the planet is necessary to reveal its periodic reflex motion. Accordingly, the study of a simulated time series of mock epochs is the logical next step towards assessing what is required to astrometrically detect an exomoon.

Fig. 1 Orbital wobble caused by the perturbing gravitational influence of an exomoon. The left panel shows the orbital trajectory of the moon-hosting planet relative to the star at the origin for different configurations, the right panels display the corresponding astrometric signals in right ascension and declination as a function of time. For visualisation purposes, the mass ratio between moon and planet, Mm/Mpl, used in these examples is unrealistically high to the point that the system would more fittingly be described as a ‘binary planet’. The plane parallel moon case corresponds exactly to the configuration given in the top right of the left hand panel. The eccentric and inclined examples were computed using the adjusted orbital elements specified in the brackets in the bottom left of the panel. The ‘0’ in the indices of the semi-major axes indicate that these values correspond to the separation between the two bodies, that is between the planet–moon barycentre and the host star as well as between the moon and the planet, respectively. In other words, they are the sums of the semi-major axes of the respective hosts and orbiting bodies around their mutual centre of mass. Since we define the inclination of the moon to be 0° and 90° for face-on and edge-on orbits as projected onto the plane of the sky, respectively, simply increasing the inclination results in a dampening and eventual complete suppression of the astrometric signal in the declination component while the right ascension component remains unaffected. For the inclined example, we adjusted the argument of periastron to add a phase shift in the signal so as to better distinguish the different cases.

Fig. 2 Signal amplitudes in the vicinity of the putative exomoon around β Pic b suggested by Poon et al. (2024) as a function of moon mass and semi-major axis for RV studies on the left and astrometric searches on the right. The grey shaded region is the parameter space where a P24-moon is likely to reside if it is indeed real. The dashed lines indicate contours delineating the provided signal strengths.

4 Modelling multi-epoch detection attempts

To probe which kind of systems a given instrument is sensitive to, we need to establish different observation sampling strategies, implement a procedure of generating mock astrometric epochs, provide safeguards against over- or underestimating the signal of a given moon, and define a robust way of computing the significance for any combination of instrument, observing strategy and system configuration. In the following, these points will be addressed individually.

4.1 Defining observing strategies

The chosen observational cadence is of consequence for the system configurations one is most sensitive to. First, it is important to note that there is no one optimal strategy. Instead, the approach one adopts would ideally be informed by the target and the period one suspects the moon to exhibit. Assessing which sampling strategy is ideal for which system is a far-reaching problem in itself. Avoiding aliasing and achieving an efficient sampling of the astrometric curve are key concerns that need to be addressed (see Madore & Freedman 2005). Doing so is beyond the scope of this work, however. We shall instead confine our study to investigating how the attainable sensitivity of different instruments behaves as a function of the number of astrometric epochs. Secondly, most targets are not observable throughout the entire year from most locations on Earth. The chosen observatory and target system position thus impose a second constraint on the observing strategy. Finally, there might be special boundary conditions set by the operational processes and time allocation procedure of the chosen observatory. For instance, on Cerro Paranal, the VLTI observation windows are scheduled to coincide with the full moon and typically last for five days each. Since at the moment there is little to no possibility of deviating from this fixed schedule we shall consider it a hard condition in the analysis to follow. Bearing the above in mind, we define two basic strategies:

1. 12 epochs: A sequence of evenly spaced observations with a cadence of one synodic month; that is, approximately 29.5 days. Once six epochs have been obtained, we pause for a period of six months during which the target is assumed to be unobservable, before observing for another six months in a row.

2. 18 epochs: Same as above but including a second pause and third sequence of six epochs.

Naturally, these approaches only represent a small subset of all possible strategies and are mainly geared towards investigating the sensitivity gain obtained from an additional set of six astrometric epochs. Beyond the above approaches, one could adopt a higher-frequency tactic where an astrometric epoch is obtained during every night of several consecutive VLTI windows. Such a ‘rapid-fire’ strategy would be especially suited to hunt for short-period moons. As was mentioned above, however, we defer a thorough signal-processing-based identification of optimal sampling strategies for different system configurations to future works.

Fig. 3 Fiducial astrometric signal in RA and Dec as a function of time for an exemplary system configuration (Ms = 1.5 M⊙, Mpl = 10 MJup, apl, 0 = 10 AU, Mm = 1 MJup, ϖ = 50 mas, am,0 = 0.1 AU, im = 0°, em = 0; see caption of Fig. 1 for index definitions). The dark grey curve indicates the fiducial signal computed using the star–planet–moon model described in Sect. 2.2.2, while the light grey line shows the zero-signal expected in the absence of a moon. The light and dark grey circles show the fiducial momentary signals in RA and Dec when employing the 12 epoch strategy defined in Sect. 4.1. In orange we indicate the associated GRAVITY mock epochs generated according to the procedure outlined in Sect. 4.2. The right hand panels show the same data but phase-folded by the fiducial orbital period.

4.2 Creating astrometric mock epochs

Based on the chosen number of epochs, we next need to generate a mock astrometric data set. To this end, we can use the models introduced in Sect. 2 as our underlying fiducial models.

They provide the exact on-sky positions of the star, planet, and moon relative to the system barycentre for any given time of observation. By extracting the planetary position according to the two models at a series of pre-defined observing times, we obtain our fiducial planet positions. These are not the positions that our instrument will measure, however. Shifting them in both RA and Dec by a noise component drawn from the Gaussian distribution 𝒩(0, σ_A), where σ_A is the measurement uncertainty of the chosen instrument, turns the fiducial planet positions into our mock astrometric epochs. Figure 3 shows the fiducial signal and generated GRAVITY mock epochs as a function of time in RA and Dec for an exemplary system configuration.

4.3 Randomising the time of periplanet passage

The fact that the astrometric signal varies periodically over time holds the danger of significantly misjudging the achievable detection significance. In case the epochs are taken at regular intervals and the chosen cadence coincides approximately with a multiple of P_m/2, where P_m is the orbital period of the moon, one runs the risk of consistently sampling the wobble at the same angular deviation. For instance, one could conceivably ‘be lucky’ and systematically observe the planet at its largest deviation from the Keplerian trajectory or ‘be unlucky’ and catch it at just the moments when the astrometric signal is minimal. To avoid this pitfall when hunting for moons in earnest, it is therefore advisable to refrain from employing an evenly spaced observing strategy in cases where one expects the moon orbit to be of a similar period. In our simulation, however, we can additionally counteract the issue by randomising the time of periplanet passage, t_{peri, m}. This effectively entails a randomisation of the planet–moon orbital phase such that an evenly spaced observing sequence can still coincidentally align with the planetary wobble for a single run but not for the entire sample of N_rand runs. Another consequence of this manoeuvre is that for a given star–planet–moon configuration we will be left with a distribution of detection significances.

4.4 Estimating the detection significance

The last ingredient required to assess our current and future detection capabilities is a framework that allows us to quantify in which parameter spaces the star–planet–moon model performs significantly better at explaining the mock astrometric epochs than the mere star–planet model. The two models are nested in the sense that the star–planet–moon model transitions into the star–planet model when the moon mass becomes negligible. That being said, the star–planet–moon model requires more free parameters, which implies that any approach to be used to compare the performance of the two models must be capable of penalising such increased complexity. Since it computes the Bayesian evidence, a nested sampling approach (see e.g. Feroz & Hobson 2008; Buchner 2021) would enable us to account for the number of free parameters in a given model. As we are planning to compute the significance achievable by different instruments using various observing strategies over a grid of target system parameters, while also randomising the time of periplanet passage for each grid point (see Sect. 4.3), such a strategy would require inordinate computational effort. Given the general character of this work as an initial feasibility assessment, we decided to defer the implementation of a nested sampling routine to a future study and focus the sensitivity estimates on comparing the χ-squared values presented by the two models at the fiducial parameter settings. This approach presumes that the fiducial solution can be converged upon using a suitable fitting routine, an assumption that is justified for the cases we are most interested in: the ones where $\mathrm{\Delta }{\chi }_{\text{red}}^2$, the difference between the reduced χ-squared values of the two models, is significant. Starting from the fiducial parameter set we ran a gradient descent routine to ensure the minimum was stable. The χ-squared values used throughout the remainder of the analysis presented here are based on these converged-upon parameter sets.

To quantify the relative performance of the two models, we applied a variant of the F-test based on the Fisher-Snedecor distribution (see e.g. Eadie et al. 1971). Following Band et al. (1997), a given $\mathrm{\Delta }{\chi }_{\text{red}}^2$ value was converted to an F-value according to

$$F=\frac{\mathrm{\Delta }{\chi }^2/\mathrm{\Delta }\nu }{{\chi }_{\text{moon}}^2/v_{\text{moon}}},$$(3)

where Δν is the difference in the degrees of freedom ν = N_data − N_fp + 1 between two models. Here, N_data and N_fp are the number of data points and free parameters, respectively. To account for the fact that each astrometric epoch corresponds to a measurement in RA and Dec, we used N_data = 2 N_epoch. The quantities in the denominator of Equation (3) are the χ-squared and degrees of freedom of the more complex model; that is, the one including the moon. From the resulting F-values, we computed the corresponding p-values via its survival function. Finally, we converted the p-values into a significance via the quantile function. Both operations were handled using the respective scipy methods (Virtanen et al. 2020).

We thus possess all the tools required to collapse the χ-squared values resulting from the comparison of the two models with a set of noisy mock data obtained by a given instrument employing a certain observing sequence of a specific star–planet–moon configuration into a single detection significance.

Fig. 4 Sensitivity curves of different current and future interferometric instruments capable of astrometric measurements. The solid and dashed lines indicate the 5 σ detection limits based on 12 and 18 epochs, respectively. The grey dots in the background mark the probed grid points. The rescaled axes in the top and on the right provide a parametrisation where the fiducial target system distance, d = 20 pc, and planet mass, Mpl = 10 MJup, have been folded in. They can thus be used to estimate the moon detection significance achievable for a genuine exoplanet of interest. Three exemplary moons detectable around AF Lep b (Mpl = (3.27 ± 0.25) MJup; Balmer et al. 2025) and HR 8799 d (Mpl = (9.3 ± 0.5) MJup; Zurlo et al. 2022) are indicated by the dotted lines. The corresponding moon masses and semi-major axes are provided.

5 Multi-epoch sensitivity curves

We applied the procedure outlined above to estimate the detection significance achievable using VLTI/GRAVITY, VLTI/PLANETES, and a hypothetical future interferometric facility with a 3 km baseline (see Table 1 for the respective planetometric measurement uncertainties), each employing the two observing strategies described in Sect. 4.1. To facilitate comparability with a large fraction of the directly imaged planet population, the fictional star–planet–moon system this investigation is based on was deliberately configured as general as possible. Accordingly, we chose a 10 M_Jup planet revolving around a 1.5 M_⊙ star on a circular, face-on orbit (i = 0°) with a semi-major axis of 10 AU and at a distance of 20 pc. The moon orbit was assumed circular and co-planar with the planetary orbit. The planetary Hill radius resulting from these fiducial parameter settings amounts to approximately 1.2 AU. These parameters remained fixed throughout the grid exploration of a two-dimensional plane spanned by the moon semi-major axis and moon mass axes, which were logarithmically sampled by 15 values from 10⁻² to 1 AU and 10⁻² to 1 M_Jup, respectively. Every grid point was sampled 100 times, generating a new set of mock epochs with a randomly drawn noise term for each iteration. This procedure ensures adequate sampling of the underlying χ-squared statistic and sufficient randomisation of the time of periplanet passage which was drawn anew for each sample (see Sect. 4.3). Accordingly, we obtained a set of 100 χ-squared values per model fit performed.

At this point, one could choose between retaining the minimum, median or maximum χ-squared value (and associated significance computed via the F-test routine) depending on whether one is interested in extracting the least favourable, median or most favourable detection scenario. To arrive at a realistic assessment of our detection capabilities we opted to retain the median χ-squared value for each model and grid point. The 5 σ detection contours resulting from the different combinations of instruments and observing strategies are visualised in Fig. 4.

Besides the number of obtained epochs, the detection significance at a given grid point in Fig. 4 is mainly governed by the expected astrometric signal amplitude. In Equation (1), we found that A ∝ a_m/d and A ∝ M_m/(M_pl + M_m). By appropriately folding in the fiducial distance, d, and the total mass, M_pl + M_m, of the fictional underlying planet–moon system, we rescaled the axes in Fig. 4 providing the resulting alternate axes along the right and top sides of the panel. These rescaled axes enable comparison of the obtained sensitivity contours with genuine exoplanets that might be considered for a follow-up aimed at unveiling an exomoon.

To give two illustrative examples, we included the cases of AF Lep b (M_pl = (3.27 ± 0.25) M_Jup; Balmer et al. 2025) and HR 8799 d (M_pl = (9.3 ± 0.5) M_Jup; Zurlo et al. 2022). Multi-plying the values along the rescaled semi-major axis axis with the distance of a target system yields the physical semi-major axis values corresponding to the significance contours. Likewise, using the mass of an exemplary target planet, one can convert the values along the mass ratio axis into the corresponding moon masses. Doing so, we find that 12 astrometric GRAVITY epochs suffice to detect a moon of 0.34 M_Jup at a separation of 0.59 AU around HR 8799 d with a confidence of 5 σ. Targeting the lower mass AF Lep b, 12 epochs could reveal a 0.14 M_Jup moon at a separation of 0.39 AU. Both of these hypothetical moons would reside within 0.5R_{Hill, pl} of their respective planetary hosts. Since more massive or more distant moons would induce a larger astrometric signal, they would manifest themselves in the data at even higher significances. Conversely, if none is present, 12 epochs would be adequate to rule out the existence of such moons around AF Lep b and HR 8799 d. The arrival of PLANETES at the VLTI will amplify our sensitivity such that moons of 0.06 M_Jup at 0.17 AU and 0.16 M_Jup at 0.25 AU were detectable at 5 σ around AF Lep b and HR 8799 d, respectively. Finally, a future 3 km baseline interferometric facility will be capable of detecting moons of 0.02 M_Jup at 0.06 AU and 0.05 M_Jup at 0.08 AU around the two planets at the same confidence. These scenarios, indicated by the dotted grey lines in Fig. 4, are only three example cases for two arbitrary exoplanets. Note that 0.02 and 0.05 M_Jup translate to approximately 6 and 16 M_Earth, respectively. Thus, the moons detectable with a 3 km baseline interferometric facility are comparable with Super-Earth exoplanets (Valencia et al. 2007). At higher moon masses the instruments would be capable of probing smaller semi-major axes. Similarly, at larger moon semi-major axes the detectable masses can be substantially lower. Targeting potential host planets of lesser mass than either AF Lep b or HR 8799 d in their respective circumstellar habitable zones may eventually teach us how prevalent potentially habitable exomoons are. The recently imaged planet candidate around α Cen A (Beichman et al. 2025; Sanghi et al. 2025) might provide a suitable first target as it combines the favourable characteristics of being nearby and low in mass while at the same time orbiting within the circumstellar habitable zone.

As long as they are comparable to the underlying fiducial system configuration on which the significance contours are based (planets on far-out, circular, face-on orbits), Fig. 4 can be used to estimate the attainable instrument-dependent detection significance for any planet of interest. Star–planet–moon configurations and viewing angles that differ from the underlying fiducial system such as eccentric or edge-on orbits warrant dedicated simulations. Intuitively, for edge-on or near-edge-on cases the signal amplitude would remain the same, even if the reflex motion of the planet were confined to tracing out a line as projected onto the plane of the sky. With the publication of the astrometric time series data in Gaia DR4, we can expect to confirm dozens of planet candidates using GRAVITY. The combination of Gaia astrometry and the direct detection will yield a large sample of narrowly constrained orbital solutions and precise dynamical masses (see Winterhalder et al. 2024) that will pave the way for an astrometric exomoon hunt.

Broadly speaking, based on 12 or 18 epochs, GRAVITY and PLANETES are sensitive at 5 σ to moons between 0.1 and a few per cent of the mass of their planetary host. According to Canup & Ward (2006), moons that have formed in situ, that is in orbit around the planetary host, should not exceed a moon-to-planet mass ratio of approximately 10⁻⁴. There thus does not exist an overlap between our current and near-future astrometric exomoon detection capabilities and the regime where these regular moons are expected to reside. Judging by Fig. 4, a kilometre-baseline facility might begin to tackle this limit – at least for comparatively large moon semi-major axes. So-called irregular moons, that is ones that were either captured by or formed as a result of an impact with their host, do not abide by this upper mass-ratio limit.

To assess the potential of an astrometric search for exomoons, it is of interest what expenditure in observation time such a study would entail. Based on experiences drawn from observations conducted within the framework of the ExoGRAVITY Large Programme (ESO ID 1104.C-0651 Lacour et al. 2020), obtaining one GRAVITY epoch can be assumed to require a total of 2 hours of telescope time including overheads. Data sets comprising 12 and 18 epochs, such as the ones the sensitivity curves in Fig. 4 are based on, would therefore carry a cost between 20 and 40 hours of VLTI UT-mode time. While this is a non-negligible investment, a potential moon discovery would make such a project highly worthwhile. Risk mitigation can be achieved by an informed target selection that hinges on clues as to the potential existence of an orbiting moon from auxiliary studies. The spectroscopic data products that result from a GRAVITY observation can act as contingency science cases in the event that no moon is present or detectable: time-series observations of exoplanet and brown dwarfs can be probed for variability or combined into a single high signal-to-noise spectrum (see the ExoGRAVITY spectral library presented in Kammerer et al. 2025).

Finally, a note on distinguishing between perturbations caused by a planet-orbiting moon and a star-orbiting interior planet, since, in general, both object classes would manifest themselves as an astrometric wobble. The expected periods are significantly different, however. Whereas the sought-after moons are expected to impart signals with periods of days to months, inner planets can be assumed to complete their orbits on typical timescales of several years (Lacour et al. 2021). There does exist an overlap between the two cases, however: undetected hot Jupiters would induce a stellar reflex motion that could mimic the signature of a moon since our astrometric epochs correspond to relative separation measurements between the host star and the targeted planet. Such niche cases could easily be handled by obtaining an RV data set of the star. If present, a signal of a similar period would strongly suggest an inner planet causing the observed astrometric wobble. An exomoon, on the other hand, would not impart any measurable stellar RV modulation.

6 Conclusions

In this feasibility study, we have demonstrated the potential of hunting for exomoons and binary planets by means of an astrometric time series. The unprecedented precision of interferometric instruments such as VLTI/GRAVITY renders the technique viable and sensitive to low-mass moons around far-out directly imaged planets, a regime where moons are expected to survive over long periods. While, contrary to the RV method, the astrometric signal drops for target systems at greater distances, if successfully applied, the method can provide constraints on all orbital parameters including the absolute dynamical mass of the moon, the key parameter upon which any characterisation of the satellite object depends. The first dedicated GRAVITY search for an exomoon targeting HD 206893 B is presented in Kral et al. (2026). We showed that, a next-generation instrument like VLTI/PLANETES will likely provide an astrometric measurement precision comparable to the signature induced by the putative exomoon around β Pic b as suggested by Poon et al. (2024).

The sensitivity curves obtained from modelling the astrometric signatures expected for different target system configurations quantify our current and future capabilities and provide a compelling science case for pursuing kilometre-baseline interferometric facilities. Future studies should concentrate on implementing a fitting framework based on a nested sampling routine in order to robustly assess the degree to which a genuine astrometric time-series suggests the existence of an exomoon around a given planet. Moreover, the signal-processing conundrum of how to optimally sample an astrometric signature of a given period should be examined. In particular, such an investigation would help answer the question of whether the strict monthly VLTI observation window allocation inhibits the detection of short-period signals and whether the current scheduling system should be reassessed.

On the whole, astrometric time-series observations provide an excellent means of detecting moons around gas giant exoplanets. Intriguingly, a future kilometre-baseline interferometric facility, used for monitoring gas giants that orbit within their respective habitable zones, will be capable of detecting and characterising exomoons with masses comparable to that of Earth.

Acknowledgements

SL acknowledges funding from the European Union (ERC AdG 101142746 PLANETES). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

References

All Tables

All Figures

Fig. 1 Orbital wobble caused by the perturbing gravitational influence of an exomoon. The left panel shows the orbital trajectory of the moon-hosting planet relative to the star at the origin for different configurations, the right panels display the corresponding astrometric signals in right ascension and declination as a function of time. For visualisation purposes, the mass ratio between moon and planet, Mm/Mpl, used in these examples is unrealistically high to the point that the system would more fittingly be described as a ‘binary planet’. The plane parallel moon case corresponds exactly to the configuration given in the top right of the left hand panel. The eccentric and inclined examples were computed using the adjusted orbital elements specified in the brackets in the bottom left of the panel. The ‘0’ in the indices of the semi-major axes indicate that these values correspond to the separation between the two bodies, that is between the planet–moon barycentre and the host star as well as between the moon and the planet, respectively. In other words, they are the sums of the semi-major axes of the respective hosts and orbiting bodies around their mutual centre of mass. Since we define the inclination of the moon to be 0° and 90° for face-on and edge-on orbits as projected onto the plane of the sky, respectively, simply increasing the inclination results in a dampening and eventual complete suppression of the astrometric signal in the declination component while the right ascension component remains unaffected. For the inclined example, we adjusted the argument of periastron to add a phase shift in the signal so as to better distinguish the different cases.

In the text

Fig. 2 Signal amplitudes in the vicinity of the putative exomoon around β Pic b suggested by Poon et al. (2024) as a function of moon mass and semi-major axis for RV studies on the left and astrometric searches on the right. The grey shaded region is the parameter space where a P24-moon is likely to reside if it is indeed real. The dashed lines indicate contours delineating the provided signal strengths.

In the text

Fig. 3 Fiducial astrometric signal in RA and Dec as a function of time for an exemplary system configuration (Ms = 1.5 M⊙, Mpl = 10 MJup, apl, 0 = 10 AU, Mm = 1 MJup, ϖ = 50 mas, am,0 = 0.1 AU, im = 0°, em = 0; see caption of Fig. 1 for index definitions). The dark grey curve indicates the fiducial signal computed using the star–planet–moon model described in Sect. 2.2.2, while the light grey line shows the zero-signal expected in the absence of a moon. The light and dark grey circles show the fiducial momentary signals in RA and Dec when employing the 12 epoch strategy defined in Sect. 4.1. In orange we indicate the associated GRAVITY mock epochs generated according to the procedure outlined in Sect. 4.2. The right hand panels show the same data but phase-folded by the fiducial orbital period.

In the text

Fig. 4 Sensitivity curves of different current and future interferometric instruments capable of astrometric measurements. The solid and dashed lines indicate the 5 σ detection limits based on 12 and 18 epochs, respectively. The grey dots in the background mark the probed grid points. The rescaled axes in the top and on the right provide a parametrisation where the fiducial target system distance, d = 20 pc, and planet mass, Mpl = 10 MJup, have been folded in. They can thus be used to estimate the moon detection significance achievable for a genuine exoplanet of interest. Three exemplary moons detectable around AF Lep b (Mpl = (3.27 ± 0.25) MJup; Balmer et al. 2025) and HR 8799 d (Mpl = (9.3 ± 0.5) MJup; Zurlo et al. 2022) are indicated by the dotted lines. The corresponding moon masses and semi-major axes are provided.

In the text