© The Authors 2024

1 Introduction

The Solar System hosts a variety of large moons, presenting environments as diverse and scientifically interesting as its planets. As the count of exoplanets found so far is in the thousands, it is natural to ask whether these extrasolar planets hide an equal number of extrasolar moons. Secondary or population-level effects of such moons have already been found (e.g. Kenworthy & Mamajek 2015; Hippke 2015; Teachey et al. 2018; Oza et al. 2019; Saillenfest et al. 2023), and so the question is not whether these moons exist, but rather what they look like, and in what systems and around which planets can they be found.

The detection and characterisation of such satellites accompanying extrasolar planets holds the potential to further our understanding of planet formation, evolution, and habitability. Different moon and planet formation scenarios and their outcomes have been linked to predictions of or requirements on, for example, the time of formation of the satellites (Cilibrasi et al. 2018), their size and mass relative to their host (Nakajima et al. 2022; Canup & Ward 2006; Hansen 2019), circumplanetary disc (CPD) composition (Batygin & Morbidelli 2020; Oberg et al. 2023), host magnetosphere (Canup & Ward 2006), host mass (Oberg et al. 2023), instellation or host migration history (Heller & Pudritz 2015a,b), and orbital properties (Li et al. 2020). While the Solar System is host to a large number of moons, the limited scenario we are presented with cannot conclusively bear witness to any of these analyses: validation of these formation studies would require the detection and characterisation of satellites of planets outside the Solar System.

Moreover, moons with the mass of Mars or greater (which are not found in the Solar System) are promising targets for habitability studies (Lammer et al. 2014; Williams et al. 1997; Dobos et al. 2022; Forgan & Dobos 2016), and the observation-ally confirmed existence of subsurface oceans on Titan (Beuthe 2015; Bills & Nimmo 2011; Baland et al. 2011), Enceladus (Beuthe 2016), among others, and the likely current or past existence of such oceans on Triton (Nimmo & Spencer 2015; Schenk et al. 2021; McKinnon & Kirk 2014; Gaeman et al. 2012) and other Solar System bodies (e.g. Hussmann et al. 2006; Burnett & Hayne 2023; Rovira-Navarro et al. 2023; Nimmo et al. 2016; Bagheri et al. 2022; Bierson & Nimmo 2022; Nimmo & Pappalardo 2016) hold promise for the potential habitability of such worlds, even at sizes observed in the Solar System. This makes Solar System-like satellites objects of interest concerning habitability in their own right, not just as extensions of the population of small icy and rocky planets orbiting stars directly.

Current searches for exomoons are limited primarily to those using transit timing variations (TTVs) and transit duration variations (TDVs), which can yield orbital and mass information on the moon (Kipping 2009a,b; Heller et al. 2014; Teachey et al. 2018, 2020; Fox & Wiegert 2021; Kipping & Yahalomi 2022; Kipping et al. 2022). For this method, unfortunately, the effects of even large moons are degenerate with those produced by other, potentially unobserved planets in the system (Fox & Wiegert 2021; Kipping & Teachey 2020). Consequently, candidate moons are met with skepticism (Teachey et al. 2020; Fox & Wiegert 2021; Kipping 2020; Tokadjian & Piro 2022), especially as they have thus far required invoking eccentric formation scenarios (Hamers & Portegies Zwart 2018; Hansen 2019; Kipping et al. 2022). Other detection and characterisation methods have been proposed, from those directly analogous to those used in planet detection such as radial velocity measurements of the host planet (Ruffio et al. 2023; Vanderburg et al. 2018) or transit spectroscopy (Kaltenegger 2010; Limbach et al. 2021) to the detection of thermal excesses in direct imaging data (Limbach & Turner 2013; Kleisioti et al. 2021, 2023) or microlensing (Han & Han 2002; Han 2008; Hwang et al. 2018). These also face major limitations, however. The first two methods, with current instrumentation, can unfortunately only detect binary-like satellites in terms of size and mass (Lazzoni et al. 2022), whereas the third requires exceptional tidal heating rates. Microlensing detections are isolated events, and therefore are difficult to confirm and cannot be followed up on.

In this analysis, we build upon previous work by Cabrera & Schneider (2007), who discussed the possibility of detecting extrasolar satellites by astrometry of the planet or photometric phenomena induced by mutually eclipsing and dimming events, and Agol et al. (2015), who put these two together and discussed the possibility of detecting extrasolar satellites by use of a method called spectroastrometry, in which one measures the differential on-sky position between the light originating from a system in different filters. A shift in the centre of light between two filters cannot be the result of a point-source or symmetric object (e.g. rings), and so explanation of such a shift in a planet signal requires the presence of a (potentially unresolved) second object. This method has previously found use in characterisation of close-separation stellar binaries and active galactic nuclei at scales well below the diffraction limit, as the centroid can be measured to precision better than the diffraction or seeing limits (Bailey 1998a,b; Porter et al. 2004; Whelan & Garcia 2008).

Bailey (1998b) already predicted that a large telescope using low-order adaptive optics to achieve image sizes of ~0.l arc second in the IR should be able to use spectroastrom-etry to make measurements to ~ 100 micro-arcsec’, which is roughly equivalent to the separations of Solar System-like moons at distances of the order of tens of parsecs. We derive improved estimates of the signal and noise for spectroastrometry of photometric points, and we show that the upcoming class of extremely large telescopes (ELTs), with the European ELT as archetype, will indeed be capable of detecting and characterising nearby mild-to-warm tidally heated exomoons (THEMs) of sizes and separations such as those found in our Solar System and compatible with moon formation theory, with relatively little regard for inclination or orientation of the satellite orbit.

In contrast, satellites with the expected mass ratios of ~10⁻⁴ compared to their host predicted by Canup & Ward (2006) at medium-to-wide planet-moon separations will remain elusive even with future instrumentation for most other methods (Ruffio et al. 2023; Lazzoni et al. 2022). Additionally, those other methods impose a restrictive orientation of the orbit of the satellite (for planet-transiting moons or radial velocity detections; Lazzoni et al. 2022) or require a fortituous transiting orientation of the host planet with respect to its star (Kipping 2009a,b). We show that spectroastrometry only imposes a weak orientation preference on observations that is complementary to that required for the radial velocity and planet-transit methods.

We structure the analysis as follows: we begin by presenting the spectroastrometric signal, noise, and some their basic and derived properties in Sects. 2.1 to 2.5, followed by a description of our benchmark scenario in Sect. 2.6. The results for this benchmark scenario are presented in Sect. 3, and the limitations of our analysis, the repercussions with respect to current-generation telescopes, and a comparison with other moon detection and characterisation methods are presented in Sect. 4. We conclude matters in Sect. 5.

2 Methods

In the following, we present the framework upon which our analysis is based. It is structured as follows: in Sects. 2.1 and 2.2 we introduce the spectroastrometric signal and some of its basic properties, followed by a discussion on the contaminating noise sources in Sect. 2.3 and the consequences on observation design in Sect. 2.4. We then derive an expression for the minimum flux a moon must have to be detectable at a given signal-to-noise ratio given any instrument specifications in Sect. 2.5, and introduce the model scenario used to evaluate the efficacy of spectroastrometry in detecting Solar System-like moons in Sect. 2.6.

2.1 The centroid and the spectroastrometric signal

In defining the spectroastrometric signal, we take an approach that differs from that taken by Agol et al. (2015); this framework allows us to compute several additional error estimates. Motivated by the idealised mathematical limit of the photon-count averaged centre of light on a detector D for ever smaller pixels, we shall define the centroid for an observation over a time interval T_F through a filter F, c_F, as the photon count-weighted integral over all positions c over the course of an observation O, during which a total of N_F photons are counted as follows: (1)

The differential photon count hitting the detector for each position c, dN_F(c), can be rewritten, accounting for the fact that the incoming photon rate is a function of the detector location it hits and of time, t. We choose to ignore the effect of the photon path through the instrument, such that we can then write dN_F(c) = I_F(c, t) dΩ dt dS with I_F the photon intensity through the filter F originating from the point c on the celestial sphere at time t, dΩ the infinitesimal area on the detector that the photons hit (expressed as solid angle on the celestial sphere), and dS the infinitesimal area of the telescope aperture that the photon field passes through. If we furthermore assume that the total photon flux originating from the part of the sky we observe through the filter F, F_F = ∫_D I_F (c, t)dΩ, is constant throughout the observation such that N_F = F_FT_FS, we obtain: (2)

where dS was taken outside the integral as we had assumed the photon flux on the detector to be independent of the photon path through the instrument. The outermost integral is ostensibly a time-average, while (3)

is the photon flux-averaged centre of light at time t. We have thus shown that the centroid as measured from an observation is the time-averaged value of the instantaneous centroid c_F (t) over the duration of that observation: c_F = 〈c_F(t)〉_T.

We note that for an I_F(c, t) that is point-symmetric about some point c₀(t) we have c_F(t) = c₀(t): motivated by this observation, we shall now assume that we are observing a planet-moon system (we note, however, that computation of the centroid of an observation according to Eq. (1) does not require this assumption), such that we can decompose the intensity originating from any given on-sky location into that originating from the planet (denoted by the subscript p) and that originating from the moon (with the corresponding subscript m). Then, by linearity of the integral: (4)

where F_Fp (respectively F_Fm) is the flux in the filter F due to the planet (respectively moon). We observe that so long as the point-spread function (PSF) for our telescope is point-symmetric, we have that the components I_Fp(c, t) and I_Fm(c, t) of our intensity are respectively point-symmetric about the on-sky position of the planet and moon, which we shall denote c_p and c_m, such that finally: (5)

analogous to Eq. (4) in Agol et al. (2015). We must note, though, that here we have shown that filter fluxes must be expressed in photon fluxes, not energy fluxes, for this relation to hold. For notational brevity, we shall henceforth write c_m(t) − c_p(t) = c_mp(t) for the on-sky projected angular separation of the two bodies. In imitation of the approach Agol et al. (2015) take to forego the necessity of a reference position, we then define the spec-troastrometric signal S_M,P as the absolute difference between the centroid location in a filter M in which we expect to observe the moon (the ‘moon filter’ ; not to be confused with the notation for absolute magnitude used in other literature) and another filter P in which we expect the planet to be dominant (the ‘planet filter’); we remark that we also differentiate between the periods over which the two observations were made, T_M and T_P respectively: (6)

where we have assumed that the time-averaged position of the planet over both the periods T_M and T_P is identical. As we shall encounter these quantities more often, let us denote the fraction of the flux in the band M (respectively P) due to the moon, (respectively ), as ƒ_Mm (respectively f_Pm).

2.2 The spectroastrometric signal for closed Keplerian orbits

For 〈c_mp(t)〉, an analytic solution exists in terms of known quantities from orbital mechanics for all closed Keplerian orbits; its derivation is described in Appendix A. One can then show that the expected value of the signal for a given observation of a moon with unconstrained inclination i and ω (i.e. a flat prior on both i and ω) is given by: (7)

with γ a parameter describing our knowledge of the inclination i and orientation (through ω) of the moon; if i and ω are unconstrained (i.e. uniformly distributed), we have p ≈ 0.842 if instead the moon is positioned in the worst possible orientation (edge-on i.e. i = π/2), we have p = 2/π; in the best possible orientation (face-on i.e. i = 0) we have p = 1. (respectively ) is the non-dimensionalised time-averaged in-orbit position over the time period T_M (respectively T_P) of the moon, and d is the system-observer distance. A derivation of Eq. (7) and the values of the quantity p as well as a more extensive derivation of (and ) are provided in Sect. C; an illustration of the geometry involved is given in Fig. 1.

It will suffice for now to know that it can be shown that the time-averaged in-orbit position depends on the orbital properties of the moon according to Eq. (8) (for a more detailed discussion of this result, the reader is referred to Appendix A for the derivation and Appendix B for a discussion on its properties): (8)

where P (not to be confused with the filter P, which only appears in subscripts) and e are the period and eccentricity of the orbit of the moon, T is the timespan over which we observe in the filter of interest, and E₀ and E₁ are the eccentric anomaly at the start and end of the observation in that same filter, respectively. For a given time since periapsis, E₀ and E₁ can be calculated by solving Kepler’s equation, for which a variety of solution methods exist; we employ the non-iterative method described by Markley (1995).

Fig. 1 Illustration of the geometry underlying the calculation of the spectroastrometric signal in the face-on case (p = 1).

2.3 Sources of spectroastrometric noise

There are, however, several sources of noise that will contaminate the signal. In the following, we shall describe relevant noise sources and their corresponding expressions, as follows: (1) photon shot noise, (2) pixel noise, (3) background and instrument flux noise and (4) pointing noise. We (will) note that each of these noise sources acts upon each of the measured centroids individually (that is, per filter) and independently. We will therefore consider in the following the noise that each of these sources imparts upon the measured centroid in a filter, noting that these noise sources are present in both the filter M and the filter P. We can then combine the two into the resulting total noise for the spectroastrometric signal in Sect. 2.3.5.

2.3.1 Photon noise

The PSF of the telescope is not just the reason by which we cannot resolve the moon and planet directly on the exposure; the resulting spread in photon arrival locations causes an inherent and unavoidable variance in the measured centroid, given that we can only ever take a finite sample of the PSF. In Appendix D.1, it is shown that under the condition of point symmetry of the PSF we have for the photon noise σ_PN: (9)

where σ_PSF is the standard deviation of the PSF for a single sample (i.e. a single photon) and N is the number of photons observed in the filter. For a diffraction-limited telescope, the Gaussian approximation of the Airy disc used by Agol et al. (2015), where σ_PSF = 0.45λ_c/D (with λ_c the central wavelength of the observation filter and D the telescope diameter), gives a useful general expression, but where available an estimate from PSF shape models would of course be preferred. As we use the ELT as archetype for the full class of extremely large telescopes (for which PSF models are not yet available), we will use the expression by Agol et al. (2015) in our simulations, but prefer to retain σ_PSF in expressions.

2.3.2 Pixel noise

Another source of noise derives from the fact that we are not probing the actual on-sky intensity distribution, but a pixelated version thereof. The finite size of the pixels on our detectors means that we lose some information about each detected photon. An upper bound for the noise that this introduces can be derived (see Appendix D.2) to be given by (10)

with α the pixel width, and N again the total number of photons observed in the filter.

2.3.3 Background and instrument flux noise

The noise in background and instrument flux will also give rise to noise in the measured centroid. In Appendix D.3, we show that an upper bound for this contribution is given by (11)

under the assumptions that the planet detection in the flux is ≥5σ, that we sample the centroid from the detector region encompassing the 3σ-region of the PSF (which contains >97% of the incoming flux) and that the noise is uncorrelated between pixels.

2.3.4 Pointing noise

The imperfect pointing accuracy of the telescope will also impact the extracted centroid; we can account for this by assuming that each observed photon has a further random offset dictated by the pointing accuracy of the telescope. In that case, the mathematical framework works out precisely as for the PSF, so long as the timescale on which the telescope pointing is affected by inaccuracies is significantly lower than the total time of the exposure, such that photons can be assumed to be independently affected. In the case of the ELT, the uncompensatable random errors vary on subsecond timescales (Rodeghiero et al. 2021), such that we deem this a reasonable assumption. METIS will achieve fine-guiding accuracies below 0.02λ/D (Brandl et al. 2021), which translates to a worst-case pointing offset below σ_PO ≈ 1 mas, such that we can expect the pointing inaccuracy σ_p to be at the very greatest (12)

This term is negligible for the ELT, but may be important or even limiting for space telescopes in particular.

This treatment then accounts for random pointing noise; systematic pointing offsets can be removed by referencing against the reference objects used. If possible with the given coron-agraphic system, the central star or another bright co-moving object would be a suitable candidate. Otherwise, a background object can be used. As such objects are in general far brighter than the planet, the corresponding spectroastrometric noise is likely to be negligible. If the telescope reacquisition pointing accuracy between different filters allows for it, referencing against other objects may even become unnecessary. This may introduce systematic bias if there is an unknown offset between filters (e.g. due to imperfect adaptive optics), however, and so this should be carefully accounted for.

2.3.5 Total noise

The fact that all noise sources go as N^−1/2 motivates us to combine them into one term, such that we can describe the total noise in a filter as σ_tot = σN^−1/2. We can consider the noise σ to be an ‘effective noise per photon’, and under the assumption that all noise sources are independent we can calculate σ to be (13)

It is instructive to remark that σ for the both filters, σ_M and σ_P, is a property inherent to a given filter on a given telescope, but it does not depend on the object to be observed. The combined noise σ affects the measured centroids of both of the filters in question. Assuming that these two measured centroids are statistically independent, the total noise in the measured spectroastrometric signal, σ_S, is then: (14)

with ε being the achromatic efficiency of the telescope and S_eff its effective surface area, where we have used that N_M = S_eff εT_MF_Mp/(1 − f_Mm) (an analogous expression of course holds for N_P).

Note that it follows directly from Eq. (14) that the best planet filter P to observe in so as to minimise the noise is the one in which is maximised; the same conclusion holds for the moon filter F, so long as the moon is expected to be sufficiently luminous in this filter, too.

2.4 Consequences for the observation time allocation

From Sect. 2.3, we can conclude an important guideline for observations: namely, one can show that given a fixed total observation time T = T_M + T_P, there exists an allocation between T_M and T_P such that σ_S as given by Eq. (14) is optimised. This time allocation can be shown to satisfy the relation (15)

The approximation is justified by noting that for a suitable planet filter P we ought to expect that ƒ_Pm ≈ 0, and if the moon does not fully dominate in M, is also reasonable. As this approximated optimal time allocation requires no a priori knowledge of any tentative moon, we find it most representative to perform our evaluation of the method using this time allocation. We wish to note that this optimal time allocation is equally valid for the noise estimate used by Agol et al. (2015), though they do not appear to have taken note of it.

2.5 Minimum required moon flux

In general, unfortunately, we have no a priori knowledge of the expected flux of moons. Moreover, we are interested not in the signal-to-noise ratio produced by a specific moon, but rather we would like to answer the converse question: the luminosity a moon ought to have such that it is detectable with at least a given signal-to-noise ratio S /N. One can take Eqs. (7) and (14) to arrive, after some manipulation, at a criterion for detectability: (16)

where (17) (18)

This is a quadratic form in ƒ_Fm, whence one can derive the following minimum flux requirement for detectability at the specified signal-to-noise ratio: (19)

with (20) (21) (22)

This, should be stressed, holds for all closed Keplerian orbits, given the expression for in the form of Eq. (8), so long as D > 0. This latter constraint follows when writing Eq. (16) as a form that is quadratic in F_Mm/F_Mp, and it is satisfied roughly when C_P ≲ 1 such that (23)

Note that the moon fluxes in either band are, of course, both unknown quantities; as ƒ_Pm is relatively low (close to zero) if a suitable planet filter P was chosen, it is more informative (and receptive to implicit computation) to put constraints on F_Mm as a function of ƒ_Pm than the other way around. If one can relate ƒ_Mm and ƒ_Pm to a set of defining parameters (say, for example, the surface temperature and radius of a moon), it is then possible to solve the constraint posed by Eq. (19) implicitly, given a set of orbital and observation parameters.

2.6 Model scenario

To assess the limits of what types of Solar System-like moons may be detectable in nearby systems in the near-future using ELT-class telescopes, we explore a benchmark scenario. We take METIS as reference instrument, given that its capability to detect Earth-sized planets around nearby stars in thermal emission has been well-established (Brandl et al. 2021); it should therefore not be problematic to reach a 5σ detection (in flux) of a giant planet, so as to satisfy the assumptions required for Eq. (11). The telescope model parameters are summarised in Table 1. For the moon filter M, we take the N2-filter as it is the longest-wavelength continuum filter on METIS, which will therefore be best-suited to capture low-temperature (<300 K) thermal emission. For the planet filter P, by contrast, we choose the near-infrared filter M’, which captures a region in which giant planets such as Jupiter, Saturn (Roman 2023) and higher-mass exoplanets (Phillips et al. 2020) are luminous, while remaining stiff to the types of temperatures expected globally for tidally heated moons (see e.g. Dobos & Turner 2015). As an example, the used spectrum for ε Indi Ab (which we shall justify shortly) is shown with a 350 K, Earth-sized moon (representing an upper bound on the hottest, largest plausible moons) for comparison in Fig. 2: it is clear that the planet still dominates wholly in M', while the moon has an appreciable contribution in N2. METIS requirements mandate that its optical system is diffraction-limited (Brandl et al. 2021), which we will therefore take as an assumption on the PSF-width. Atmospheric transmission is set to an achromatic 80%, based on a representative scenario simulated using ESO’s SkyCalc (Jones et al. 2013; Noll et al. 2012). We set an observation time of 6h, representing a length of night that can be reasonably expected to occur regularly at the ELT site (Lombardi et al. 2009), though we note that with a scaling of the noise with T^−1/2, our results are relatively stiff to the total observation time.

We model a reference moon-planet system based on ε Indi Ab, the nearest Jupiter analogue known imageable by JWST (and thus, presumably, ELT), and its nominal parameters as determined by Feng et al. (2019). The parameters necessary for this analysis are given in Table 1; based on this planet mass and age, a spectrum is interpolated from the chemical equilibrium atmosphere-spectra in the publically available atmosphere library ATMO2020 by Phillips et al. (2020). As their simulated atmospheres only provide data for planets of this mass up to ages of ~3 Gyr, we cannot take an age in the age range of 3.7–4.3 provided by Feng et al. (2019), though we note that preliminary results using ATMO2020 show that in this case the spectroastro-metric signal-to-noise ratio is increasing with age of the planet (i.e. with decreasing temperature): this makes an age of 3 Gyr conservative.

In imitation of Solar System moons (with the notable exception of Titan), we assume the moon to be an airless icy or rocky body. In such a case, the observed brightness temperature in the infrared agree relatively well with the surface temperature (to within a couple ~10 K), regardless of surface material (see e.g. Hu et al. 2012; Whittaker et al. 2022); integrated brightness temperatures of Jupiter’s icy moons Ganymede and Callisto agree relatively well with their expected surface temperatures (Squyres 1980). We therefore choose to model the moon as a black body, such that the only free parameters that remain to explore are (1) its size (but notably not directly its mass), (2) its surface temperature and (3) its orbital properties. Each of these we vary over ranges made plausible by a combination of observed Solar System moons and moon formation theory, so as to explore the signal-to-noise ratio of the spectroastrometric signal. In particular, we vary the size between an Io-radius and an Earth-radius, which, allowing the satellite to be either icy or rocky, covers the full range of expected moon masses of 10⁻⁵–10⁻³ host masses for large satellites proposed in literature (e.g. Canup & Ward 2006; Cilibrasi et al. 2018; Moraes et al. 2018). We vary the semi-major axis over the range seen for the Galilean satellites.

While analysis of time-series observations of the system under the assumption of a moon on a Keplerian orbit may allow one to produce higher-significance detections (see Sect. 4.2.2 in Agol et al. 2015 for an example), we prefer to analyse solely the spectroastrometric signal in one observation, as it requires no assumptions on the nature of the signal. A high-significance detection of a spectroastrometric signal would therefore provide incontrovertible evidence of a signal of astrophysical origin, regardless of whether it is a moon or not; further analysis and conclusions on the nature of this signal can then follow.

Fig. 2 Spectrum for ε Indi Ab, generated from ATMO2020 assuming a 3.25 MJup, 3 Gyr planet, and equilibrium chemistry; the chosen M' and N2 bands are marked. Also shown is the extreme scenario of a 350 K, Earth-sized blackbody moon, for which it is clear that the planet still dominates by ~2 orders of magnitude in M', while the moon contributes an appreciable fraction of the flux in N2, such that a spectroastrometric signal may be observed between the two filters. We note that the region from 5.5 to 7.5 µm, though it has promising absorption regions, is inaccessible from the ground.

3 Results

We structure our results as follows; we examine the minimum temperature required for detectability in circular orbits as a function of semi-major axis in Sect. 3.1, and perform a first-order analysis of the required eccentricities for these temperatures under the assumption that the detectability is only marginally affected by low eccentricities. We validate this assumption in Sect. 3.2, and finally perform an analysis of the attainable signal-to-noise ratios as a function of moon blackbody temperature and distance for two sample moons of given size and semi-major axis in Sect. 3.3.

3.1 Minimum moon temperature for detectability

Figure 3 shows the minimum required blackbody temperature of a moon for various sizes (corresponding to Io, Mars and the Earth) to be detectable at 5σ in a circular orbit around ε Indi Ab, assuming a total of 6 h of observation time. Also shown is a first-order estimate of the eccentricities required for each of the various size bodies at each semi-major axis to reach the given minimum temperature by assuming radiative equilibrium of a blackbody moon experiencing viscoelastic dissipation as given by Segatz et al. (1988) and using the value of −Im(k₂) given for Io by Lainey (2016). For comparison, the semi-major axes and eccentricities of the seven largest Solar System moons are marked; as Triton currently has a circular orbit but is known to have migrated to its current position after capture over a timespan of ~1 Gyr (Ross & Schubert 1990; McKinnon & Kirk 2014), its constant-angular momentum post-capture evolution as described by Ross & Schubert (1990) is also drawn. The value of f_Mm corresponding to these minimum temperatures is also shown; over the explored semi-major axis values, f_Pm for the minimum detectable temperatures is negligible, and so the required values of f_Mm coincide. For this same reason, the spectroastrometric noise is nearly constant across all semi-major axes at ~ 0.0026 mas, corresponding to a minimum detectable spectroastrometric signal of ~ 0.013 mas.

3.2 Eccentricity effects

Figure 4 shows the minimum moon blackbody temperature for two sample moons (that will be explored and justified in further detail in Sect. 3.3) at an ε Indi Ab-like system-observer distance at a variety of eccentricities and orbital phases covering all possibilities for the inclination-averaged case (p ≈ 0.842). The green and blue lines, indicating the minimum required temperature for e = 0 and ±10 K deviations from that contour, show that the required temperature for moons with Solar System-like eccentricities (e ≲ 0.1) is in general well-approximated by the required temperature for the circular case throughout the orbit; additionally, the required temperature for detectability is in fact in all cases lower than for the circular case over > 50% of the orbit.

Fig. 3 Minimum blackbody temperature and corresponding moon flux fraction in M for a total observation time of 6 h on ELT/METIS as a function of semi-major axis for 5σ detectability of a moon around ε Indi Ab for an Io-sized, Mars-sized, and Earth-sized moon (top), with a corresponding first-order estimate of the required eccentricity (bottom). The grey region delineates the area bounded by the edge-on and face-on cases (p = 2/π and p = 1); the red and blue lines correspond to the inclination-averaged case (p ≲ 0.842). Semi-major axes and eccentricities of the seven largest Solar System moons Ganymede, Titan, Callisto, Io, the Moon, Europa and Triton are marked by their first letter(s). The post-capture migration path of Triton assuming conservation of angular momentum (such as in Ross & Schubert 1990) is also marked to illustrate the observational potential of captured moons. The moons used for further analysis in Sects. 3.2 and 3.3 are marked in green in the top plot. For all three moons over all explored semi-major axes, the minimum 5σ-detectable temperature corresponds to a spectroastrometric signal of roughly 0.013 mas.

Fig. 4 Minimum temperature for 5σ detectability as a function of the observed orbital phase and eccentricity for an Io-like moon with an Io-like semi-major axis (top) and for an Earth-like moon with a Callisto-like semi-major axis (bottom). The green and blue contours indicate the contour corresponding to the nominal (zero-eccentricity) required temperature (as marked on the top plot in Fig. 3), and the 10 K deviations from this temperature, respectively; in both cases, these are also marked on the colourbar. We note that in all non-zero eccentricity cases the temperature required for observability is in fact lower than for the circular case over >50% of the orbit.

3.3 Applicability to distant systems

The results presented in Sects. 3.1 and 3.2 suggest that there is a set of nearby systems in which planet-like moons of directly imageable giant planets can be detected. This motivates us to explore the extent of this viable detection space for two scenarios, representing the end-points of what Lazzoni et al. (2022) deem the class of planet-like satellites (namely those formed by core accretion in the CPD, analogous to most Solar System satellites): (1) an Io-sized moon at Io-like separation from its host and (2) an Earth-sized moon at a Callisto-like separation from its host. The first requires the least number of assumptions, as it is a type of moon observed in the Solar System, and tidal heating-mechanisms through which it might reach temperatures to be luminous in the IR over observable timescales are well-established (e.g. Dobos & Turner 2015; Rovira-Navarro et al. 2021). The second would be on the larger end of what one expects to see around a host like ε Indi Ab according to moon formation studies such as Canup & Ward (2006), though more recently Cilibrasi et al. (2018) produced results that seem to suggest such masses may be attainable. Additionally, where fixed-Q tidal theory would predict that such a far-out moon should be unlikely to experience significant tidal interactions, the recently proposed paradigm of resonance locking would allow such far-out moons to experience tidal interactions still (Fuller et al. 2016; Lainey et al. 2020). These two cases should therefore in principle bracket the full range of plausible planet-like moons that might be expected to be luminous in the IR; the signal-to-noise ratio as a function of blackbody temperature and system-observer distance for these two cases are illustrated in Fig. 5.

Fig. 5 Signal-to-noise ratio for an Io-sized moon on an Io-like orbit (left) and an Earth-sized moon on a Callisto-like orbit (right) around an ε Indi Ab-like planet as a function of blackbody temperature and system-observer distance (we note that the axis scale and datum differ between the two moons). The grey, blue, and red lines mark the 1σ, 3σ, and 5σ boundaries, respectively.

4 Discussion

Figures 3 and 5 indicate that spectroastrometry may be able to provide detections of large icy moons (with surface temperatures as low as T_s ≲ 150 K) or even smaller, hot rocky moons (with surface temperatures as low as T_s ≲ 300 K) in nearby systems; such temperatures are plausible for tidally heated moons on orbits comparable to modern-day Solar System moons (e.g. Dobos & Turner 2015; Forgan & Dobos 2016; Rovira-Navarro et al. 2021) and easily exceeded by conditions as have been presumed to exist in the early Solar System (e.g. Ross & Schubert 1990), potentially over Gyr timescales (McKinnon & Benner 1990; Lunine & Nolan 1992). The possibility of observing large icy moons at low temperatures is interesting, as the large icy moons Ganymede and Callisto in our Solar System have (and are stable at) surface (and brightness) temperatures of order ~150 K (Squyres 1980; though a large part of this is contributed by solar flux), and the largest and most massive moons are thought to form at or beyond the ice line of their CPD (Heller & Pudritz 2015b), therefore being icy. The results for smaller hot moons are also of interest, as Rovira-Navarro et al. (2021) have shown that plausible scenarios may arise in which such temperatures are maintained for timescales ≳ 1 Gyr.

We will discuss these results in several contexts: in Sect. 4.1, we compare our framework against that in previous literature, and we discuss several caveats to our new formulation in Sect. 4.2. One may also wonder if perhaps current-generation infrared telescopes may be able to detect objects in this manner, too: we discuss what this would take in Sect. 4.3, and follow this up with a discussion on the applicability of spectroastrometry to moon-planet systems other than the ones we have studied in this paper in Sect. 4.4. Finally, we make a comparison of spectroastrometry against other methods in Sect. 4.5 and discuss what combinations with other methods may possibly allow unambiguous detection, characterisation and confirmation of moon candidates in Sect. 4.6.

4.1 Comparison to previous work

While we have expanded the framework provided by Agol et al. (2015) to work for closed Keplerian orbits in general (Sect. 2.2) and to include additional noise bounds (Sect. 2.3), a great deal of conclusions in their work still hold. In particular, the scaling of the signal-to-noise ratio as S /N ∝ (ϵT)^1/2d⁻² remains generally applicable, and under the assumptions that the detector pixel size is designed to sample the PSF well (i.e. α ∝ σ_PSF), that the pointing stability noise is negligible and that the telescope is diffraction-limited we also recover the telescope-diameter proportionality S /N ∝ D². Additionally, while Agol et al. (2015) do not appear to explicitly have taken note of this, (1) in their formulation the noise is a function of the spectral properties of the moon, but not of its orbit (excepting transits), (2) the spec-troastrometric signal-to-noise ratio of a single observation can be straightforwardly calculated for a given observation without requiring any a priori knowledge on the nature of the observed system and (3) there exists an optimal time allocation that minimises the noise of an observation that can be well-approximated by quantities that can be known or estimated a priori (Sect. 2.4). Each of these conclusions still hold in our expanded formulation.

There are several additions we have gained over this previous work, however. We have shown that the viability of spectroas-trometry is relatively insensitive to the orbital parameters of an exomoon with the exception of its semi-major axis and particularly high eccentricities (well exceeding those found in the Solar System). We found an expression for the best-case and worst-case inclinations as well as the expected spectroastromet-ric signal for a flat prior on moon inclinations. Additionally, one can show from Eq. (8) that in the limit of large T/P, , such that circular or small-eccentricity orbits become invisible. Further analysis of Eq. (8) (e.g. by analysing the pre-factor P/(πT) sin ((E₁ − E₀)/2) in Eq. (A.6)) shows that this already occurs for T/P ≳ 1. Hence, we should not expect to be able to observe moons with periods on the order of or shorter than our observation time in the relevant filter, both because their semi-major axis and the resulting orbital motion is small, which was already apparent from the work by Agol et al. (2015), but also because their signal is washed out by this time-averaging effect for T/ P ≳ 1. Conversely, highly eccentric moons on short-period orbits may instead become visible precisely because of this time-averaging effect. While for ground-based telescopes observation times are unlikely to exceed orbital periods, this is a plausible scenario for space-based observations, which would require far longer observation time to reach a performance similar to large-diameter ground telescopes.

Additionally, we have shown that while background and instrument noise in centroid measurements for a significant (>5σ) detection of the planet is likely to be of similar magnitude to the photon noise, the noise due to a pixelated detector (for a well-sampled PSF) and pointing stability of the telescope in question are unlikely to provide major contributions to the centroid noise, and therefore to the spectroastrometric noise.

4.2 Limitations on this formulation

A major limitation of the method as developed in this manner is the lack of an expression for speckle noise, which ostensibly does not satisfy the assumption of independence between pixels required for Eq. (11). Until this effect can be quantified, these results should be taken only as representative for moons of far-out or free-floating planets, where speckle noise is negligible, though fortunately the majority of large exomoons are expected to occur around wide-orbit planets (Dobos et al. 2021; Heller & Pudritz 2015a), beyond the snowline (Inderbitzi et al. 2020). If an expression or estimate for speckle noise in spectroastrometry can be found, we might be able to push spectroastrometry to find moons of close-in planets in reflected light, too, such as the case of an Earth-Moon analogue studied by Agol et al. (2015), though this would require a different (augmented) model of the moons’ spectral energy distribution incorporating reflected light.

A second limitation for close-in planets is the lack of inclusion of the planet-moon barycentre movement over the observations in our formulation. While the movement thereof between observations can be accounted for if the orbital properties of the planet about its host are known, the movement throughout the observation will likely have an effect that is worth quantifying in future studies: using the orbital parameters given by Feng et al. (2023), the worst-case angular velocity (i.e. for a face-on orbit at pericentre) for ε Indi Ab is on the order of 0.1 mas h⁻¹, which is of the same order as spectroastrometric effects. The authors deem it likely that this can be accounted for in a similar fashion to the time-averaging method employed in Appendix A, either numerically or possibly analytically, by inserting the barycentre motion into the integral as an offset. For wide-separation planets or edge-on planets observed at their most distant apparent separation, this effect is negligible, however. As our neglect of speckle noise means that wide-separation planets should be preferred to start with, for this effect to become strong will most likely require that speckle noise be dealt with, first.

Another matter that needs mentioning is the fact that there are objects that cannot immediately be distinguished from moons; background objects may, per chance, induce a spectroas-trometric signal. Follow-up observations will, however, reject such candidates without issue, given the rapid motion expected of moons about their host that such background objects will not display. Additionally, thermal emission from asymmetric ring systems around planets will cause a signal that may mimic thermal emission from a hot moon; one could foresee such asymmetric rings arising in a miniature version of the scenario invoked by András & Rieke (2020) for the debris cloud formerly identified with Formalhaut b, for example, though notably Formalhaut b was not detectable in the infrared. As such asymmetries on a Keplerian orbit should be expected to dissipate over the order of several orbits, a signal due to such occurrences should not be expected to remain stable over observational timescales; additionally, the presence of rings (asymmetric or not) may be excluded altogether by the method proposed by Lazzoni et al. (2020), for example.

Spots or asymmetric patterns on the planet may produce an astrometric signal, too: we note that the Great Red Spot on Jupiter is dim in the infrared (Ge et al. 2019), and so we should by analogy expect any such effects to stem from larger-scale variability. For Jupiter, such rotational variability is on the order of 1% in flux around the 10.5 µm region, whereas in the 5 µm region this is roughly 20% (Ge et al. 2019). For the moon filter we can thus safely disregard this variability: the planet filter does warrant a small calculation. In the worst-case scenario of a pole-on planet with all this variability concentrated on a single spot at the equator, and with a Jupiter-like radius and rotation period of 10 h, the resultant astrometric signal will be on the order of 10⁻² mas. As the 20% variability on Jupiter is the result of large-scale variability, not spot-like localised variability (Ge et al. 2019), we expect that this effect will be significantly smaller in practice. Additionally, higher-mass planets are expected to have a lower rotational period (Snellen et al. 2014), which will dampen the signal strongly if the rotational period becomes of the same order as the observation time in a filter (see Eq. (B.3)). Nonetheless, these variability effects may thus, in the worst case, be present at an order of magnitude that will contaminate the most sensitive of our results: follow-up observations in different planet-bands (perhaps with lower variability) will then provide a conclusive result on the presence or absence of a moon.

It is thus possible to exclude background objects, (asymmetric) ring systems, or planet variability either immediately or after follow-up observations. Consequently, we expect that moons are the only objects or phenomena that can be responsible for long-term spectroastrometric signals, though any putative signal will need follow-up observations for confirmation.

4.3 Spectroastrometric capabilities of current-generation telescopes

As these results show that observing tidally heated exomoons on future IR telescopes is plausible, one might wonder whether current-generation IR telescopes have the capabilities required to observe tidally heated moons in nearby systems. The prime candidate here, of course, is JWST; unfortunately, the scaling of the signal-to-noise ratio with D² suggests that the observation time required would be on the order of days. JWST would thus only be sensitive to hot moons around nearby directly imageable planets that are less luminous than ε Indi Ab, none of which are currently known.

Additionally, while the fine-guidance pointing stability of JWST is of an order similar to that of ELT, target reacquisition (which would be required when switching between filters) only has accuracies on the order of mas (Hartig & Lallo 2022), which introduces systematic error. This may be remedied if one uses astrometric observations (i.e. using a single filter) of the moon throughout its orbit instead, but in that case it still remains to be shown that the measured effect is due to motion of a moon rather than due to a third object disturbing the Keplerian motion of the planet. In fortituous circumstances, one might also be able to use sufficiently bright background objects (that can be astrometrically positioned to the required accuracy) to reference the observations in the two filters against each other, but this requires that such objects be available and fixed. We leave it to future work to determine the viability of each of these methods to perform spectroastrometric measurements using JWST.

4.4 Detectability of other moon-planet systems

In this analysis, we have considered Solar System-like moons around a planet analogous to ε Indi Ab, as it is the nearest known Jupiter analogue that is directly imageable. However, such a system need not be optimal in terms of detectability of any moons; moreover, it is likely that with JWST operational many more directly imageable planets will be discovered in nearby system that had previously eluded detection. It is thus productive to consider what other planets might be found, or perhaps what types of moons might be detectable around them.

One such category is lower-mass planets; exploratory results using the code in this paper suggest that even lower-temperature moons will be observable around lower-mass planets, as the relative flux of the moon then increases appreciable in the moon filter. If located in a suitable nearby system (such that stellar contamination is not a problem), this may thus allow the detection of true Solar System-analogues: that is to say, Solar System-analogue planets accompanied by Solar System-size and Solar System-temperature moons. Down to Saturn-masses this is hardly problematic, but for lower-mass planets (transitioning into the ice giant-regime) the spectral feature at 5 µm starts to become muted (Roman 2023; Linder et al. 2019), such that the planet will no longer outshine its moon in that band; another planet filter than M’ (which was selected in a fairly ad hoc fashion) is required for ice giants. Detection of moons around such planets would thus require a more robust filter selection procedure than we have followed in our analysis of ε Indi Ab. Hence, we propose that in anticipation of the detection of directly image-able low-mass giant planets in nearby systems by JWST further research is done on filter performance over a wide variety of planet masses and ages.

A second category of systems does not concern the type of planet, but rather the type of moon; while we have for the purposes of our analysis explicitly remained in the familiar realm of Solar System-like moons, current exomoon candidates are ostensibly unlike the moons seen in our Solar System. The candidates put forward by Teachey et al. (2018), Teachey & Kipping (2018), Lazzoni et al. (2020) and Kipping et al. (2022) are perhaps rather thought of as binary objects; this poses unique challenges for spectroastrometry in that the spectra of two gas giant planets in a mutual orbit are likely to possess similar features qualitatively over a broad range of wavelengths. As both objects are, without requiring tidal heating, likely to be relatively bright, however, their spectroastrometric signal may still be detectable, and if so, perhaps at system distances further away than Solar System-like moons. Given the importance of such objects in informing moon formation theory (e.g. Hamers & Portegies Zwart 2018; Hansen 2019; Moraes & Neto 2020), it is worth exploring whether any filter combination might be able to detect or rule out the presence of any such moons; perhaps, even, whether a filter combination may be able to detect both Solar System-like moons and binary-like objects.

4.5 Comparing against other moon detection methods

Over other methods, spectroastrometry provides an additional boon in the fact that, if the barycentre motion discussed in Sect. 4.2 is accounted for, the astrometry produced is of a sufficient level to potentially detect the orbital motion of nearby exoplanets over a single observation; even in the absence of a moon detection, the data is therefore still useful in novel ways. Another major advantage that spectroastrometry has over the current suite of standard exomoon detection efforts is the lack of degeneracies (if observed over multiple epochs), such as those induced by unseen second planets in the case of the methods proposed for transiting planets (Fox & Wiegert 2021; Kipping & Teachey 2020); as one measures explicitly the movement of the moon about its planet, these third-object effects are unambiguously removed. Finally, spectroastrometry provides the possibility of repeatable observations of nearby moons as small and nearly as cold as those observed in the Solar System, far below temperatures required for other direct imaging methods for tidally heated moons (e.g. Kleisioti et al. 2023), while being relatively insensitive to the inclination or orientation of the moon or its planet.

4.6 Synergies with other methods

Further modelling of the moon or its orbit can further enhance the significance of results and produce predictions or constraints on quantities such as the semi-major axis, period, size, host planet mass and flux distribution between satellite and host (Agol et al. 2015) and consequently even interior properties of the satellite (e.g. Kleisioti et al. 2023) that can be (independently) verified with other methods, such as radial velocity or astrometric measurements of the planet or such as TTVs/TDVs, auxiliary stellar transits or moon-transits of the planet, if the system is fortituously oriented. Especially noteworthy is the fact that at satellite inclinations where spectroastrometry fares poorest (near edge-on), transits of their host and radial velocity measurements become viable, which makes spectroastrometry a useful complement to those methods. Finally, mutual events (Cabrera & Schneider 2007) observed in the infrared might provide information on the composition of the planet and moon as well as provide independent constraints on the moon orbit (Schneider et al. 2015).

Another detail of particular interest is the fact that spectroas-trometry provides an independent constraint on moon and planet fluxes as well as the moon orbit that may be combined with those derivable from moon-planet models such as discussed in Kleisioti et al. (2023), which can be produced from the same direct imaging data that spectroastrometry can be performed on. In combination with the results derivable from other methods, it is therefore possible to come to a near-full characterisation of the moon in question in terms of its orbit, composition, size and surface conditions.

We would therefore like to emphasise the possible capabilities of simultaneous use of spectroastrometry over multiple epochs, photometric modelling, astrometry and transit observations of moon-transits of the planet, as all of these should in principle be possible with the same series of direct-imaging observations. Given the set of mutually independent estimates for similar parameters available between the set of these, self-consistency of any candidate can be straightforwardly checked. We therefore strongly recommend that the viability of this combination is evaluated in future studies.

5 Conclusion

Spectroastrometry has previously been shown to be a promising method for detection of Earth-Moon-like systems or Earth-like moons orbiting Jovians by Agol et al. (2015). We have shown that there is a further class of satellites, tidally heated moons analogous to those in the Solar System, that are observable in nearby systems using this method with next-generation ground telescopes. Illustrated by two example systems motivated by moons observed in our current-day Solar System, we see that for nearby systems even large icy bodies or hot bodies comparable in size to those seen in our Solar System may be observable.

In showing this, we have derived an expression for the spec-troastrometric signal that covers all closed Keplerian orbits, along with the orbital motion of the moon throughout the observation, which shows that the efficacy of spectroastrometry is only weakly dependent on the orbital properties of an exomoon, except for moons with periods equal to or lower than the observation time, which are unlikely to be detected. Moreover, we have derived additional conservative noise estimates for noise due to (1) background and instrument noise, (2) a pixelated detector and (3) pointing inaccuracies. This now allows evaluation of spectroastrometry as a method on telescopes other than the ideal photon-noise limited space telescopes assumed in previous research and without requiring a priori any assumptions on the orbit or nature of any potential moon.

Acknowledgements

We would like to thank Amy Louca and Christiaan van Buchem for the productive discussions on the prospects of detecting tidally heated moons and Neptune-like moons that led to the research presented in this paper, and Matthew Kenworthy for motivating us to work towards publishing the result. We must also acknowledge the insightful comments provided by the anonymous referee, which have greatly strengthened the presented discussion of spectroastrometry.

Appendix A Derivation of the general time-averaged angular separation between moon and planet for closed Keplerian orbits

Let us suppose that we take an observation over a time period [t₀, t₀ + T] of a moon-planet system in mutual (closed) orbit around one another with Keplerian elements (α, i, e, Ω, ω, M₀) at a system-observer distance d. Here we use the astronomical conventions: the inclination is measured between the orbital angular momentum vector and the line-of-sight (such that the reference plane is oriented perpendicular to the line-of-sight) and the longitude of the ascending node is measured from celestial north eastward (i.e. anti-clockwise). To predict the measured spectroastrometric signal (Eq. 6) we then need the time-averaged on-sky position of the moon with respect to its host, 〈c_mp (t)〉: as the on-sky projection of this time-averaged position are precisely the x- and y-components in the reference coordinate system (i.e. the coordinates in the plane perpendicular to the line-of-sight along celestial north and east, respectively), we have: (A.1)

where ξ and η denote the coordinates in the orbital plane with ξ measured from barycentre to pericentre and η oriented along the latus rectum so as to have a right-handed frame; moreover, R_i(α) denotes the rotation matrix for the rotation by an angle a about the i-th coordinate axis; these rotations can be pulled out of the integral by linearity, meaning that in fact we need only calculate the average position of the moon within its orbital plane and transform this by its orbital elements to the on-sky plane in the usual way. For brevity, let us denote the integral at the end of the expression (which is in fact the position of the moon in the orbital plane averaged over time) (including a will make dependent on a only through the orbital period, as we will later see), such that: (A.2)

From basic orbital mechanics, we observe that: (A.3)

where we have used that , introducing the gravitational parameter µ; moreover, we have used that r = α(1 − e²)(1 + ecos θ)⁻¹ and that for the period P of an orbit we have (for each of these, one can consult any orbital mechanics textbook e.g. Curtis (2013)). We obtain: (A.4)

It is not immediately clear how to go about solving the integral in Eq. A.4, but in general one can integrate rational integrands of trigonometric functions through the Weierstrass substitution x = tan (θ/2) (see e.g. Stewart (2015), p. 502). We follow a similar route, but also note that for elliptical and circular orbits we might as well immediately make use of the relation between the eccentric anomaly E and the true anomaly θ, (A.5)

(again, see e.g. Curtis (2013). In this form we effectively perform the Weierstrass substitution x = tan (θ/2) followed by the substitution which after some algebra reduces expression A.4 to a set of standard integrals which can be evaluated to yield: (A.6)

where we define E₁ = E(t₀ + T) for notational convenience, and note that in fact and , such that all semi-major axis-dependency that was still contained in P can be parametrised in terms of period phase. By using the eccentric anomaly-true anomaly relation rather than solely the Weierstrass substitution we have our result solely in terms of quantities from orbital mechanics that are either known or that can be obtained from Kepler’s equation for which efficient, well-known solution computation methods exist (e.g. Mikkola 1987; Markley 1995; Fukushima 1997; Zechmeister 2018).

Appendix B Properties of the time-averaged angular separation

We remark that up until e ≈ 0.43 to within 10%, and that up until e ≈ 0.77 we have to within 10%. For the vast majority of plausible eccentricities, therefore, is well-approximated (i.e. to within 10%) as: (B.1)

As E₁ − E₀ ~ 2πT/P (from Kepler’s equation), we see that for eccentricities e ≲ 0.4 the dynamic part (comprising the first two terms) of is dominated by a term of magnitude pointing roughly in a direction , with a secondary (correction) term of order ~ e/4sinc(2πT/P) in the direction , where we adopt the convention that sincx = x⁻¹ sin x. We thus see that this second term is only important when T ≳ P and e is close to 1 , but also that its magnitude is never greater than the last, static component, which starts to play an important role for eccentricities of order e ≳ 0.1. A useful approximation for first-order estimate purposes (though too inaccurate for computational purposes) for is therefore given by: (B.2)

In the circular limit (e = 0) this becomes exact, giving the particularly simple form: (B.3)

where we have defined θ₀ = θ(t₀).

An interesting conclusion can be drawn from these expressions. As we observe the moon at an arbitrary time (i.e. we can consider the mean anomaly at which we observe the moon to be a uniformly distributed random variable) we can compute that the expected value of the component of 〈c_mp(t)〉 along the major axis of the orbit has magnitude (that along the latus rectum, as one should expect given symmetry, is zero), regardless of the observation time. In fact, in the limit as T/P grows large, the dynamic terms vanish and extremely eccentric orbits when averaged over time will produce a consistent offset with respect to the planet. Elliptic orbits will thus produce a bias along their major axis. Circular (or near-circular) orbits, as expected, will have a mean position that is centred on the moon-planet barycentre, and so to yield a signal we must observe them over a sufficiently short period that their signal is not averaged out. In practice, this means that it will be difficult to observe moons that require an observation time greater than their period, except for particularly eccentric orbits.

Appendix C Computing the spectroastrometric signal for closed Keplerian orbits

To obtain a less unwieldy expression for the spectroastrometric signal, we assume that the observations are taken back-to-back, first the filter M and then the filter P. In that case, we obtain for the definition of the spectroastrometric signal (Eq. 6) that: (C.1)

While this expression looks as though we have not made any progress with respect to Eq. 6 aside from a rearranging of terms, recall that in the form of Eq. A.6 we have an analytical expression for and . If we then define θ_M (respectively θ_P) to be the direction of the vector (respectively ) in the ξ − η plane, we find upon a set of applications of the sine and cosine laws in the triangle formed by and their difference that: (C.2)

where we will call the ‘dimension-less signal’; the signal if we had observed the orbit face-on in units of : (C.3)

and α is the unique angle satisfying simultaneously (C.4) (C.5)

or equivalently, satisfying the relation: (C.6)

where atan2(y, x) is the 2-argument arctangent. We note that , but that it is not a function of i and ω. Hence, all dependency on i and ω is now expressed explicitly. Moreover, we observe that for ƒ_Mm ⪢ ƒ_Pm we have that α ≈ 0, and that for moons with e = 0 we have that α is fixed throughout its orbit.

This is nearly in a computable form, except for the fact that we do not know a priori the inclination nor the orientation of the orbit of the moon (as we are trying to establish whether there is a moon in the first place). To ameliorate this, let us take a flat prior on the argument of periapsis, ω, which takes values in [0,2π], and marginalise over this prior. To include circular orbits, we will take ω to be the reference point on the orbit for e = 0; this, of course, is equally well-suited to a flat prior and so follows the same mathematical argument. For a face-on orbit (where the line of nodes is degenerate, and so ω is undefined), the same argument holds if we measure ω from celestial north instead. This reasoning thus includes all closed orbits.

We justify the flat prior by noting that the argument of peri-apsis is subject to precession, which for the major satellites of Saturn, for example, is known to occur on timescales of only days to thousands of years (Jacobson 2022); moreover, theory predicts for the apsidal precession rate a value independent of orientation (Greenberg 1981), such that the argument of periapsis effectively becomes a random variable that is uniformly distributed in time (of course, these timescales are in general not short enough to warrant treating it as such for consequent observations of the same moon). Therefore we can marginalise over ω, and obtain for S_M,P (in a small abuse of notation we shall prefer to denote the signal marginalised over ω and later i as S_M,P, too); (C.7)

where is the complete elliptic integral of the second kind, which has no solution in terms of elementary functions but has been well-studied and is included in most numerical programming toolkits. If one has knowledge on the inclination of any putative moon, this expression suffices, and there are good arguments to see why one would; formation models predict, for example, that satellites that formed in situ will have negligible inclinations with respect to their host due to inclination damping in the circumplanetary disc (see e.g. Moraes et al. 2018), and most major regular Solar System moons do indeed have near-equatorial inclinations (Musotto et al. 2002; Jacobson 2022) (though Triton forms a notable exception; see e.g. McKinnon & Kirk 2014); however, out of the Solar System giant planets only Jupiter has an equatorial plane relatively close to the ecliptic (Archinal et al. 2018), and it is predicted to leave this state over time (Saillenfest et al. 2020). Currently no population-level information about the obliquity of extrasolar planets is known, and recent work has shown that mechanisms exist by which planets may be tilted to obliquities near perpendicular to their orbital angular momentum by migrating moons (Saillenfest et al. 2022, 2023), tilting the plane of their satellite system simultaneously.

Until such population-level obliquity information becomes available (which could consequently inform moon searches on a per-planet level), then, one will not know the (likely) inclination of the moon and so we must go one step further. Taking a flat prior on the inclination over i ∈ [0,π], we can marginalise over the inclination of the moon, finally arriving at the expected value of the signal for a given semi-major axis and period: (C.8) (C.9) (C.10)

where is the complete integral of the first kind; the expression for the integral of E in terms of (that is, Eq. C.9) is due to Byrd & Friedman (1971), Eq. 531.02. This shows that the overall effect of the inclination and orientation of the orbit, while not negligible, is also not too severe. In the best case, with sin i = 0, we have (independent of ω), but even in the very worst case for an edge-on orbit with sin i = 1, we find that such that the signal is reduced to ~ 63.7% its best-case value. We do remark that in this case we have maintained the marginalisation over ω; for a particularly unfortunate realisation of both i and ω it is indeed possible that the signal at a given time is lower, yet even for sin i = 1 the moon will spend the greater part of its time in the region where this multiplicative constant is between and 1. In the circular case, we note that the dimensionless signal Ŝ_M,P is independent of the time since periapsis, and so the marginal-isation over the argument of periapsis effectively amounted to a marginalisation over mean anomaly (i.e. time), too. In that case the expected value in time for the edge-on case is precisely a fraction 2/π of the face-on value, as a result of the orbital sampling effect previously noted in another context by Heller (2014); the moon spends more time at a greater projected distance than it does close to conjunction.

Overall, when observing a random system with no a priori knowledge on the inclination nor orientation of any potential moon, we expect to observe any signal at ~ 84.2% of its best-case value. For now we shall parametrise this value as γ, such that with (i.e. of order unity). We then remark that for a circular orbit, the time-dependency of the signal has dropped out during our marginalisation over the argument of periapsis already; for an elliptical orbit, we could marginalise over t₀ ~ U([0, P]), still, but choose not to, as the resulting integral must be evaluated numerically and does not give any particular insight. The expression that we will use for now thus becomes: (C.11)

where we will explicitly note that such that now only the shape of the orbit (in the form of a and e), the time at which we observe with respect to last pericentre passage t₀ and the spectral characteristics ƒ_Mm and ƒ_Pm are of interest, and all other quantities are known or have been marginalised out.

Appendix D Derivation of the noise sources

Most of the noise sources as summarised in Sec. 2.3, to the authors’ knowledge, have not yet appeared in previous literature on spectroastrometry, as that has mostly been aimed at spec-troastrometry using spectra, not photometry (hence, of course, the name), of high-S/N binary sources (see e.g. Bailey 1998b,a; Porter et al. 2004; Whelan & Garcia 2008), which do not require such precise noise estimates. Therefore, we feel this warrants a more detailed explanation of their mathematical origin.

Appendix D.1 Derivation of the photon noise

Though an expression for the photon noise appears in the work by Agol et al. (2015), we feel that a more rigorous explanation is warranted in the context of the more rigorous derivations we give of the other noise sources. Let us assume we observe the source in a filter over a given time period, in which we observe N photons, of which ƒ_mN are due to the moon and (1 − ƒ_m)N due to its host planet (as we will only be discussing a single filter, we will forego the corresponding subscript in the moon flux fraction). For the sake of convenience, let us label the photons i ∈ [1, N] with the photons i ∈ [1, ƒ_mN] originating from the moon, and the photons i ∈ [ƒ_mN + 1, N] originating from the planet. If we then denote the measured location of origin of the photon i as c_i, we can denote total the variance due to the PSF photon noise in the centroid σ_PN as: (D.1)

where c_F is the centroid position over the measurement, which is given by Eq. 5; as we have already solved the effect of the bodies’ movement with time in Secs. A and C, we will neglect this effect for now such that c_F = ƒ_mc_m + (1 − ƒ_m)c_p. We then observe that we can write (c_i − c_F) differently depending on whether the photon i originates from the moon or planet; (D.2)

where we recall that c_mp = c_m − c_p. We can therefore divide the summands of Eq. D.1 into several cases: (1) i = j ∈ [1, ƒ_mN], (2) i = j ∈ [ƒ_mN + 1, N], (3) i ≠ j ∈ [1, ƒ_mN], (4) i ≠ j ∈ [ƒ_mN + 1, N], (5) i ∈ [1, ƒ_mN], j ∈ [ƒ_mN + 1, N] and finally (6) j ∈ [1, ƒ_mN], i ∈ [ƒ_mN + 1,N]. Of course, it must be noted that these cases appear in pairs: (1) and (2), (3) and (4), (5) and (6), the calculation of which is nearly identical. Let us then set for brevity S (i, j) = 𝔼 [(c_i − c_F) · (c_j − c_F)]. If we then (crucially) assume that the PSF of our telescope is point-symmetric such that 𝔼 [(c_i − c_b) · c_mp] = 0 for both b = m and b = p, we obtain for each of the six cases the following; (D.3)

where we have used the fact that for both b = m and b = p, where σ_PSF is the photon shot noise for a single photon (encapsulated in the PSF distribution), and moreover we have used the fact that any two photons are independent. In Eq. D.1 we note then that there are f_mN copies of case (1), (1 - f_m)N copies of case (2), f_mN(f_mN − 1) copies of case (3), (1 − f_m)N((1 − f_m)N − 1) copies of case (4) and a total of 2 f_m(1 − f_m)N² copies of case (5) and (6) combined. Hence, upon summing as many copies of each term we find that indeed such that , in agreement with Agol et al. (2015). What we have gained, however, is the notion that all we need is a point-symmetric PSF, which we already required to show that the measured centroid converges to the geometric centroid in Sec. 2.1.

Appendix D.2 Derivation of the pixel noise

With a robust estimate for the photon noise which accounts for the finite number of photons we count, we are now licensed to discretise the integral in Eq. 1 over the photons. Hence, the ‘actual centroid’ for the photons observed over the course of the observation (that is, the average location of origin), which we will call c_ph, becomes: (D.4)

where c_ph,i· is the ‘centroid’ (i.e. location of arrival on the detector) of the i-th photon that was detected overall while c_ph,i(p) is the location of arrival of the i-th photon that was detected on a pixel p; the reason for this re-indexing will become clear shortly. N is the total number of photons that have struck the detector in this filter, and N(p) is the number of photons that struck pixel p. However, with conventional detectors we do not record the location of arrival of the photons; instead, we record that it arrived in a given pixel on our detector. Hence, the total ‘pixelated’ centroid that we compute from the observation data, c_px,tot, is given by: (D.5)

where c_px(p) is the location of arrival we assign to photons that arrive in the pixel p (we note that this need not be the centre a priori, but we shall see that this follows immediately if we require convergence regardless of the on-sky intensity distribution). We are interested in the ‘pixel noise’: that is, the noise caused by the fact that we assign all photons in the pixel p the location of arrival c_px(p) rather than their actual location of arrival c_ph,i·(p), however: we therefore compute the difference of c_ph and c_px,tot: (D.6)

It follows upon taking expectations by linearity that 𝔼 [c_ph − c_px,tot] = 0 is satisfied regardless of the intensity distribution when 𝔼 [c_ph.i(p) − c_px(p)] = 0 for all i and p individually. As the on-sky intensity distribution of our source ought not care about the orientation or pixel size of our detector, we should expect that c_ph.i(p) − c_px(p) is uniformly distributed throughout the pixel for each p, such that 𝔼 [c_ph.i(p) − c_px(p)] = 0 if and only if c_px(p) is precisely the centre of the pixel, which is therefore a requirement for convergence we must take into account. The pixel noise σ_px then follows precisely from taking the expectation value of the square of Eq. D.6: (D.7)

where the last equality follows from expansion of the square of the sum over p into a double-indexed sum over the cross-products of (p, q) (introducing the second index q) and we have defined an auxiliary function S (p): (D.8)

As we assume any two photons to be statistically independent, any series of photons landing in one pixel must equally well be independent from those landing in another pixel. Hence, 𝔼[S (p)S (q)] = 𝔼[S (p)]𝔼[S (q)]; moreover, . As a result, all (p, q) cross-terms in Eq. D.7 vanish, and we find: (D.9)

where we can again rewrite the square of a sum as an expansion into a double sum over cross-terms. As all photons are independent and 𝔼 [c_ph.i(p) − c_px(p)] = 0 all cross-terms vanish again and so we are left with (D.10)

for which we can produce a worst-case upper bound (that is, a conservative estimate) by noting that regardless of the intensity distribution of the source, (where α is the pixel width), as no point in the pixel can be further away from the pixel centre than its vertex, at a distance . A conservative estimate for the pixel noise immediately follows, as then: (D.11)

such that . A less conservative but reasonable estimate can be obtained by assuming that the c_ph,i are uniformly distributed through the pixels (as the intensity distribution should be independent of the orientation and pixel size of our detector); in that case, one obtains instead. As the former estimate is guaranteed to be conservative, though, we shall maintain . Given that, as mentioned, the intensity distribution and our detector properties ought to be independent, it is fair to assume that this effect is independent from the photon shot noise.

Appendix D.3 Derivation of the background and instrument flux noise

A third effect that will affect the centroid we observe is the inherent random fluctuations in the noise within the area on the detector that we sample. Effectively, this is nothing other than the photon noise term for the noise, but as the noise is evenly distributed over the background sky (or rather, we assume it to be) this term is not contained in our expression for the photon noise from Sec. Appendix D.1, nor is the derivation fully analogous. We shall therefore derive it separately.

Let us consider the centroid as measured over a noisy observation, c_m, where the noise-corrected observed photon count over each pixel p can be expressed as N(p) = N_nƒ (p) + N_n(p) − µ_n where N(p) is the photon count due to the source (i.e. planet or moon) in the pixel p, N_nƒ is the noise-free photon count in that pixel and N_n(p) is the photon count due to the noise in that same pixel; we will assume that N_n(p) ~ Pois(µ_n) for all p, where µ_n is the mean noise per pixel. Correspondingly, we shall take N, N_nƒ and N_n (without the argument p) to be the total number of noise-corrected (i.e. measured), noise-free and noise photons over all pixels that we sample. We will moreover assume that µ_n has been obtained from an analysis of the source-free parts of the exposure and is therefore known. As we have taken into account the discretisation of our centroid calculation in the photons (through the photon noise) and over the pixels (through the pixel noise) we can now write for the noisy centroid c_m that (D.12)

for which we notice that for large N (when N ⪢ |N_n − N_pxµ_n|), we have (D.13)

with c_nƒ the noise-free centroid. Hence, we have in that case that (D.14)

such that we obtain (1) by taking expectations and using the fact that 𝔼 [N_n(p) − µ_n] = 0 that the measured centroid is unbiased (i.e. still converges to the true centroid) and (2) upon squaring and then taking expectations an expression for the variance due to the background and instrument flux noise, : (D.15)

whence we have, assuming that the noise fluxes in any two pixels are independent, that (D.16)

but as N_n(p) ~ Pois(µ_n), we have that the variance of N_n(p) (which is precisely what the expectation term in Eq. D.16 is) is equal to µ_n, such that we have (D.17)

where we have used that N_n ≈ µ_nN_px. If we then set the area per pixel α² = ΔΩ and the total area covered by all pixels in our sample region Ω = N_pxΔΩ, we have (D.18)

We then recognise (1) the ratio N/N_n as the signal-to-noise ratio of the planet detection in the exposure, which we shall set to 5 as a lower bound, and (2) the sum as the (polar) second moment of area of the pixel region over which we sample with respect to the origin of our coordinate system J = ∫ ∫ c² dΩ; we should therefore set our coordinate system at the centre of the region over which we sample. Moreover, this region ought to minimise the quantity J/Ω; hence, there is clearly a trade-off between sampling a sufficiently large region so as to achieve the greatest photon count from the planet and moon on the one hand, and minimising the ratio J/Ω on the other hand. A full numerical analysis could potentially yield an optimal area, but we shall simply take at least the 3σ-enclosing region for the PSF (which contains > 97% of the flux). As we have a pixelated region, we cannot achieve the optimal shape for the lowest J/Ω (which would be a circle), but we can give an ‘upper bound’ or worst-case shape that will always include the 3σ-enclosing region, which is then a square of sidelengths [6σ_PSF/α]α (the ceiling function is necessary, as we can only sample in full pixels), such that we have an upper bound for J/Ω = [6σ_PSF/α]²α²/6. Hence, we have (D.19)

whence (D.20)

which we shall therefore adopt as a conservative upper bound. We note that while this is an upper bound, for a complete measurement one can in practice calculate an exact estimate for this term a posteriori (as then the signal-to-noise ratio of the detection in the flux and the sampling region are known exactly).

References

All Tables

All Figures

	Fig. 1 Illustration of the geometry underlying the calculation of the spectroastrometric signal in the face-on case (p = 1).
In the text

Fig. 2 Spectrum for ε Indi Ab, generated from ATMO2020 assuming a 3.25 MJup, 3 Gyr planet, and equilibrium chemistry; the chosen M' and N2 bands are marked. Also shown is the extreme scenario of a 350 K, Earth-sized blackbody moon, for which it is clear that the planet still dominates by ~2 orders of magnitude in M', while the moon contributes an appreciable fraction of the flux in N2, such that a spectroastrometric signal may be observed between the two filters. We note that the region from 5.5 to 7.5 µm, though it has promising absorption regions, is inaccessible from the ground.

In the text

Fig. 3 Minimum blackbody temperature and corresponding moon flux fraction in M for a total observation time of 6 h on ELT/METIS as a function of semi-major axis for 5σ detectability of a moon around ε Indi Ab for an Io-sized, Mars-sized, and Earth-sized moon (top), with a corresponding first-order estimate of the required eccentricity (bottom). The grey region delineates the area bounded by the edge-on and face-on cases (p = 2/π and p = 1); the red and blue lines correspond to the inclination-averaged case (p ≲ 0.842). Semi-major axes and eccentricities of the seven largest Solar System moons Ganymede, Titan, Callisto, Io, the Moon, Europa and Triton are marked by their first letter(s). The post-capture migration path of Triton assuming conservation of angular momentum (such as in Ross & Schubert 1990) is also marked to illustrate the observational potential of captured moons. The moons used for further analysis in Sects. 3.2 and 3.3 are marked in green in the top plot. For all three moons over all explored semi-major axes, the minimum 5σ-detectable temperature corresponds to a spectroastrometric signal of roughly 0.013 mas.

In the text

Fig. 4 Minimum temperature for 5σ detectability as a function of the observed orbital phase and eccentricity for an Io-like moon with an Io-like semi-major axis (top) and for an Earth-like moon with a Callisto-like semi-major axis (bottom). The green and blue contours indicate the contour corresponding to the nominal (zero-eccentricity) required temperature (as marked on the top plot in Fig. 3), and the 10 K deviations from this temperature, respectively; in both cases, these are also marked on the colourbar. We note that in all non-zero eccentricity cases the temperature required for observability is in fact lower than for the circular case over >50% of the orbit.

In the text

	Fig. 5 Signal-to-noise ratio for an Io-sized moon on an Io-like orbit (left) and an Earth-sized moon on a Callisto-like orbit (right) around an ε Indi Ab-like planet as a function of blackbody temperature and system-observer distance (we note that the axis scale and datum differ between the two moons). The grey, blue, and red lines mark the 1σ, 3σ, and 5σ boundaries, respectively.
In the text