New Light on Dark Extended Lenses with the Roman Space Telescope

Dark matter (DM) is the predominant constituent of all matter in the Universe. However, it has been detected only through its gravitational interactions and its nature remains mysterious (e.g., Gelmini 2015; Bertone & Hooper 2018). Significant efforts have been devoted to searching for both particle DM (e.g., Bertone et al. 2005) as well as macroscopic DM consisting of localized pointlike objects, the canonical example being primordial black holes (PBHs; see, e.g., Sasaki et al. 2018; Carr et al. 2021; Green & Kavanagh 2021; Escrivà et al. 2022 for a review). In the broad mass range of ∼10⁻¹⁰–10 M_⊙, microlensing surveys by Subaru Hyper-Suprime Cam (Niikura et al. 2019b), OGLE (Niikura et al. 2019a), and MACHO/EROS (Alcock et al. 1998) have dominated the sensitivity reach for probing macroscopic DM.

Beyond pointlike macroscopic DM, many motivated theories instead predict DM to be composed of extended macroscopic configurations, which have received less attention. Examples include DM halo substructures (Erickcek & Sigurdson 2011; Barenboim & Rasero 2014; Fan et al. 2014; Dror et al. 2018) and axion miniclusters with axion stars (Kolb & Tkachev 1993; Eggemeier & Niemeyer 2019). Gravitational lensing has been identified as a promising method for probing extended DM substructures (Zackrisson & Riehm 2010), with recent analyses demonstrating that existing microlensing surveys can sensitively constrain a variety of models (Fairbairn et al. 2018; Croon et al. 2020a, 2020b; Cai et al. 2023), and probes complementary parameter space to other techniques, such as halometry from astrometry (Van Tilburg et al. 2018).

The launch of the Nancy Grace Roman Space Telescope (Roman; Spergel et al. 2015) in 2027 will usher in a new era of astronomy. Among its high-priority science requirements is the Galactic Bulge Time Domain Survey (GBTDS), a microlensing survey with unprecedented sensitivity to nonluminous astrophysical bodies. Though this survey is nominally intended to provide a census of Galactic exoplanets, it also has the potential to provide a new window into DM at macroscopic scales. While its sensitivity to pointlike DM sources has been discussed previously in the literature (Pruett et al. 2022; DeRocco et al. 2024; Fardeen et al. 2023), no prior study has addressed its sensitivity to extended DM distributions.

In this work, we perform the first sensitivity analysis of the Roman Space Telescope to extended DM structures. As we demonstrate, the GBTDS will be able to probe scenarios with such substructures many orders of magnitude below existing constraints, reaching fractional contributions to the total DM energy density as low as f_DM ∼ 10⁻⁶ and exploring new, motivated regions of parameter space.

Extended DM configurations arise in many theories beyond the Standard Model and manifest in a variety of forms. We analyze Roman's sensitivity for three qualitatively distinct, well-motivated examples of extended DM distributions. Here, we consider each resulting extended lens configuration to be radially symmetric and characterized by its compactness, quantified in terms of R₉₀, the radius in which 90% of the mass of the lens is enclosed. Though we choose to focus on three fiducial models, our analysis can be readily extended to account for other DM sources.

2.1. Navarro–Frenk–White Subhalos

A hierarchy of DM halos ranging from galactic masses down to Earth-mass scales is a signature prediction of standard Λ-Cold Dark Matter (ΛCDM) cosmology. There is a large uncertainty on the low-mass end of this spectrum; hunting for these smaller subhalos provides a key test of the ΛCDM paradigm. Beyond this, various cosmological theories that extend the Standard Model predict DM subhalo formation (Erickcek & Sigurdson 2011; Barenboim & Rasero 2014; Fan et al. 2014; Dror et al. 2018) at a range of scales. Characterizing the abundance of low-mass DM subhalos could therefore provide insight into physics beyond the Standard Model.

Here, we choose to model subhalos via a characteristic Navarro–Frenk–White (NFW; Navarro et al. 1996) density profile that is known to fit simulations and galactic-scale observations:

$$\begin{eqnarray}&&{\rho }_{\mathrm{NFW}}(r,{\rho }_{s},{R}_{s})=\displaystyle \frac{{\rho }_{s}}{(r/{R}_{s}){\left(1+(r/{R}_{s})\right)}^{2}},\end{eqnarray} \tag{ 1 }$$

where ρ_s and R_s denote the scale density and scale radius, respectively. While the total mass of this profile is formally divergent, we cut off the distribution at 100 R_s in keeping with the existing literature (Croon et al. 2020a). The compactness of the lens can be recast in terms of R₉₀, which is ≈69 R_s for this mass distribution.

While we choose to employ this distribution for our subhalos, our analysis can be readily extended to other profiles, such as those associated with secondary infall (Bertschinger 1985).

2.2. Halo-dressed PBH-like Configurations

The unprecedented sensitivity of the GBTDS to pointlike objects such as PBHs has been recently established in Pruett et al. (2022) and DeRocco et al. (2024) and PBH microlensing in this mass range can test a variety of intriguing theoretical models (e.g. Kusenko et al. 2020; De Luca et al. 2021; Sugiyama et al. 2021; Kawana et al. 2023; Lu et al. 2023; Inomata et al. 2023). However, when PBHs constitute a subdominant fraction of DM, they can become surrounded by massive DM halos formed from the accretion of background DM via self-similar infall (Bertschinger 1985); such a PBH is said to be "dressed" by its enclosing halo. The growth of this halo occurs primarily during a matter-dominated era.

Numerical analyses (Bertschinger 1985; Berezinsky et al. 2013; Adamek et al. 2019; Serpico et al. 2020; Boudaud et al. 2021) indicate that the resulting DM halos take a universal form well described by the falling power law,

$$\begin{eqnarray}&&{\rho }_{\mathrm{dPBH}}(r,{\rho }_{s},{R}_{s})={\rho }_{s}{\left(\displaystyle \frac{r}{{R}_{s}}\right)}^{-9/4},\end{eqnarray} \tag{ 2 }$$

where ρ_s and R_s are scale parameters. As with the NFW subhalo case, we can recast the scale radius R_s into R₉₀ ≈ 86.9 R_s . The mass of the halo M_h and central black hole M_BH are related via (Mack et al. 2007; Ricotti et al. 2008)

$$\begin{eqnarray}\begin{array}{rcl}{M}_{h} & \simeq & \,3\left(\displaystyle \frac{1000}{1+{z}_{c}}\right){M}_{\mathrm{BH}}\\ {R}_{h} & \simeq & \,0.019\,\mathrm{pc}{\left(\displaystyle \frac{{M}_{h}}{1\,{M}_{\odot }}\right)}^{1/3}\left(\displaystyle \frac{1000}{1+{z}_{c}}\right),\end{array}\end{eqnarray} \tag{ 3 }$$

where we take z_c ≃ 30 as the redshift at which the growth of the halo becomes inefficient due to interactions with nonlinear cosmic structures.

We note that for masses of interest to the Roman GBTDS (M_h ≈ 10⁻⁹–10 M_⊙), Equation (3) implies that R₉₀ ≫ R_E, where R_E is the Einstein radius (see Section 3 below). Therefore, dressed PBHs with such compactness will ultimately be challenging to discriminate from bare PBHs by Roman. As stressed in Ricotti & Gould (2009), microlensing lightcurves can efficiently distinguish dressed PBHs if the halo mass within an Einstein radius is at least a significant fraction of PBH mass, a criterion that is not met for these values of R₉₀. (We note that alternative lensing methods have been suggested to test dressed PBH scenarios, especially around the solar mass range, including lensing of cosmological fast radio bursts (Oguri et al. 2023) and gravitational waves (Gil Choi et al. 2023).)

While this demonstrates that very diffuse DM objects will be difficult to probe with Roman, we adopt a similar strategy to Cai et al. (2023) and present results for PBHs enclosed by "dressing" for which the extent of the halo has been decoupled from the PBH mass for illustration. We therefore call this scenario "dressed PBH-like," to indicate its phenomenological nature.

2.3. Axion/Boson Stars

Axions and, more generally, axion-like particles (ALPs), are pseudoscalar fields that are theoretically well-motivated (see, e.g., Adams et al. 2022 for review) and appear in many extensions of the Standard Model. In the early Universe, both ALPs and new scalar degrees of freedom have been shown to condense into macroscopic bound solitonic stars (see, e.g., Kaup 1968; Lee & Pang 1989; Seidel & Suen 1994; Schunck & Mielke 2003; Liebling & Palenzuela 2023). Among other formation scenarios, this is expected to occur generically within the cores of DM substructures such as diffuse axion miniclusters (Kolb & Tkachev 1993; Eggemeier & Niemeyer 2019). Axion stars have been associated with a broad range of observational signatures, such as bosenova explosions (e.g., Eby et al. 2016, 2022; Helfer et al. 2017; Levkov et al. 2017; Arakawa et al. 2023) as well as microlensing (Fujikura et al. 2021).

The density profile of an axion or boson star does not afford a closed-form analytic solution and must instead be calculated numerically by solving the coupled Schrodinger–Poisson equations. Rather than solving this for a particular microphysical model, we instead adopt a phenomenological model of the ground state solution, which has been shown to be well-described by (Schiappacasse & Hertzberg 2018)

$$\begin{eqnarray}&&{\rho }_{\mathrm{ax}}(r,{\rho }_{s},{R}_{s})=\displaystyle \frac{3{\rho }_{s}}{{\pi }^{3}{R}_{s}^{3}}{{\rm{sech}} }^{2}\left(\displaystyle \frac{r}{{R}_{s}}\right).\end{eqnarray} \tag{ 4 }$$

For this profile, we find R₉₀ ≃ 2.8 R_s and truncate the profile at 20 R_s .

In contrast to the dressed PBH scenario, simple axion models provide a mass–radius relation for stable axion stars (Schiappacasse & Hertzberg 2018; Visinelli et al. 2018; Sugiyama et al. 2023) that lead to R₉₀ ≪ R_E. Hence, this particular subclass of axion/boson star appears effectively pointlike to Roman and the GBTDS will be particularly sensitive to these scenarios.

Gravitational microlensing (Paczynski 1986), the apparent magnification of a luminous source by the gravitational field of a lensing mass, is one of the strongest observational probes of macroscopic DM in the mass range ≈10⁻¹¹ M_⊙–10 M_⊙ (Alcock et al. 1998; Niikura et al. 2019a, 2019b). This magnification depends on the impact parameter of the event u, the transverse distance in the lensing plane between the center of the source and the center of the lens. Throughout this work, we normalize the impact parameter by the Einstein radius of the lens,

$$\begin{eqnarray}&&{R}_{{\rm{E}}}=\sqrt{\displaystyle \frac{4{{GMD}}_{L}(1-{D}_{L}/{D}_{S})}{{c}^{2}}},\end{eqnarray} \tag{ 5 }$$

where M is the mass of the lens and D_S and D_L are the distance from the Earth to the source and lens, respectively. This quantity sets the typical transverse scale in the plane of the lens over which the lensing effect is appreciable, θ_E = R_E/D_L . The angular size of the source is likewise given by θ_⋆ = R_⋆/D_S , where R_⋆ is the radius of the source.

When both the source and the lens can be approximated as points (the "point-source point-lens" (PSPL) regime), the magnification curve follows a well-known analytic form given by (Nakamura & Deguchi 1999)

$$\begin{eqnarray}&&{A}_{\mathrm{PSPL}}(u)=\displaystyle \frac{{u}^{2}+2}{u\sqrt{{u}^{2}+4}},\end{eqnarray} \tag{ 6 }$$

where A is the observed magnification. In general, however, the angular size of the source and lens may not be negligible. The magnification can no longer be computed analytically and one must instead compute the point-source, finite-lens magnification A_PSFL numerically, then integrate this over the angular extent of the source. One can derive A_PSFL(u) via the lensing equation:

$$\begin{eqnarray}&&u=v-\displaystyle \frac{1}{v}\displaystyle \frac{M(v)}{{M}_{\mathrm{tot}}},\end{eqnarray} \tag{ 7 }$$

where v is the transverse distance from the center of the lens to the image, once again normalized to R_E, and M(v) is the projected mass contained within a radius v in the plane of the lens; hence,

$$\begin{eqnarray}&&\displaystyle \frac{M(v)}{{M}_{\mathrm{tot}}}=\displaystyle \frac{{\int }_{0}^{v}{dv}^{\prime} v^{\prime} {\int }_{-\infty }^{\infty }{dz}\,\rho (\sqrt{{z}^{2}+v{{\prime} }^{2}})}{{\int }_{0}^{\infty }{dv}^{\prime} v^{\prime} {\int }_{-\infty }^{\infty }{dz}\,\rho (\sqrt{{z}^{2}+v{{\prime} }^{2}})}.\end{eqnarray} \tag{ 8 }$$

Here, ρ(r) is the radial mass distribution of the lens. As discussed in the previous section, this function may take many forms but is typically controlled by a scale parameter that determines the lens compactness. Here, we choose r₉₀, the radius in which 90% of the mass is enclosed, in units of the Einstein radius. Note that we adopt the convention of lowercase names for length scales in units of R_E or θ_E, e.g., R₉₀ = r₉₀ R_E.

For a given mass distribution, Equation (7) can be solved to yield v(u), the image location as a function of the impact parameter. This function is often multivalued, with different branches corresponding to different images. Therefore, one must compute the magnification associated with each image separately and sum them to find the total magnification. This yields

$$\begin{eqnarray}&&{A}_{\mathrm{PSFL}}(u,{r}_{90})=\displaystyle \sum _{i}\left|\displaystyle \frac{{v}_{i}(u)}{u}\displaystyle \frac{{{dv}}_{i}(u)}{{du}}\right|,\end{eqnarray} \tag{ 9 }$$

where i denotes a particular image. At points where the slope of v(u) becomes infinite, the magnification diverges. This feature is known as a caustic (see, e.g., Hurtado et al. 2014; Karamazov et al. 2021 for a discussion of their presence in the lightcurves of spherical extended lenses). If these caustic features can be resolved, their location and shape contain information about the underlying mass profile, providing a unique opportunity to uncover the microphysical nature of the lens. We illustrate this possibility in Figure 1. However, we leave a full characterization of Roman's sensitivity to extracting these features to future work.

Zoom In Zoom Out Reset image size

Having numerically calculated A_PSFL, the ultimate magnification curve, taking into account both finite-source effects and the extended lens, is given by an integral over the angular extent of the source. Working in polar coordinates (r, ϕ) with the origin aligned to the center of the source, this is given by (Witt & Mao 1994; Matsunaga & Yamamoto 2006; Sugiyama et al. 2020)

$$\begin{eqnarray}\begin{array}{rcl}{A}_{\mathrm{FSFL}}(u,\rho ,{r}_{90}) & \equiv & \displaystyle \frac{1}{\pi {\rho }^{2}}{\displaystyle \int }_{0}^{\rho }{dr}{\displaystyle \int }_{0}^{2\pi }d\phi \,r\,{A}_{\mathrm{PSFL}}\\ & & \times \,\left(\sqrt{{r}^{2}+{u}^{2}-2{ur}\cos (\phi )},{r}_{90}\right),\end{array}\end{eqnarray} \tag{ 10 }$$

where ρ ≡ θ_⋆/θ_E.

Given this curve, one can implicitly solve for the maximum impact parameter u_T that yields a "detectable" event, where we define detectability as exceeding some threshold magnification A_T ≡ A_FSFL(u_T ). This sets the integration range over which we compute our event rate, as discussed in the following section.

Here, we describe our methodology for estimating the expected number of extended-lens events that Roman will detect during the GBTDS. We begin by computing the differential expected rate of detectable events, which is given by (Batista et al. 2011)

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Gamma }}}{{{dD}}_{L}\,{{dt}}_{\mathrm{dur}}\,{{du}}_{0}}={f}_{\mathrm{DM}}\displaystyle \frac{2}{\sqrt{{u}_{T}^{2}-{u}_{0}^{2}}}\displaystyle \frac{{v}_{T}^{4}}{{v}_{c}^{2}}\exp \left[-\displaystyle \frac{{v}_{T}^{2}}{{v}_{c}^{2}}\right]\displaystyle \frac{{\rho }_{M}}{M}\varepsilon ({t}_{\mathrm{dur}}),\end{eqnarray} \tag{ 11 }$$

where f_DM is the fraction of DM mass density composed of the extended-lens population of interest, ρ_M is the mass density of lenses, u₀ is the impact parameter at the point of closest approach, and ε(t_dur) is the detection efficiency (Johnson et al. 2020). Here, v_c is the typical velocity dispersion of the dark lenses, calculated assuming the distribution of DM in the Galactic halo follows an NFW profile (see Equation (1)),

$$\begin{eqnarray}&&{v}_{c}(r)=\sqrt{\displaystyle \frac{{{GM}}_{\mathrm{NFW}}(\lt r)}{r}}.\end{eqnarray} \tag{ 12 }$$

M_NFW(<r) is the mass enclosed within radius r of the Galactic center, where we take ρ₀ = 4.88 × 10⁶ M_⊙ kpc⁻³ and R_s = 21.5 kpc following Klypin et al. (2001). The relative source-lens velocity is given by

$$\begin{eqnarray}&&{v}_{T}=2{\theta }_{{\rm{E}}}{D}_{L}\sqrt{{u}_{T}^{2}-{u}_{0}^{2}}/{t}_{\mathrm{dur}},\end{eqnarray} \tag{ 13 }$$

where the duration of the lensing event, t_dur, is defined as the longest continuous interval for which the magnification exceeds the threshold for detection (A > A_T ). Often A_T = 1.34 is chosen as a fiducial threshold as it corresponds to u_T = r_E in the PSPL limit. However, this threshold is likely overly conservative for Roman, which is expected to achieve photometric sensitivity sufficient to resolve percent-level changes in flux (Johnson et al. 2020).

As such, we determine A_T by using more realistic detection criteria that are consistent with existing analyses of Roman sensitivity (Johnson et al. 2020). We require at least six measurements with magnifications at least 3σ above the baseline, where the baseline depends on both the source flux and blending from neighboring stars. This corresponds to

$$\begin{eqnarray}&&{A}_{T}=\displaystyle \frac{1+3\sigma -{{\rm{\Gamma }}}_{c}}{1-{{\rm{\Gamma }}}_{c}},\end{eqnarray} \tag{ 14 }$$

where ${{\rm{\Gamma }}}_{c}\equiv \tfrac{{f}_{b}}{{f}_{s}+{f}_{b}}$ is the flux contamination, indicating the fraction of the total flux contributed by the source's neighbors, and σ is the magnitude-dependent photometric precision, which we take from Figures 5 and 9 of Wilson et al. (2023), respectively. Finally, we weight the magnification-dependent threshold by the stellar population of the GBTDS using Figure 15 of Wilson et al. (2023), resulting in a population-weighted detection threshold of A_T ≈ 1.05.

The total detected event rate, Γ, is then calculated by integrating

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Gamma }} & = & 2{\displaystyle \int }_{0}^{{D}_{S}}{{dD}}_{L}{\displaystyle \int }_{0}^{{u}_{T}}{{du}}_{0}{\displaystyle \int }_{{t}_{\min }}^{{t}_{\max }}{{dt}}_{\mathrm{dur}}{f}_{\mathrm{DM}}\\ & & \times \,\displaystyle \frac{1}{{\sqrt{{u}_{T}^{2}-{u}_{0}}}^{2}}\displaystyle \frac{{v}_{T}^{4}}{{v}_{c}^{2}}\exp \left[-\displaystyle \frac{{v}_{T}^{2}}{{v}_{c}^{2}}\right]\displaystyle \frac{{\rho }_{M}}{M}\varepsilon ({t}_{\mathrm{dur}}),\end{array}\end{eqnarray} \tag{ 15 }$$

where ${t}_{\min }=90\,\mathrm{minutes}$ and ${t}_{\max }=72\,\mathrm{days}$. The lower cutoff, ${t}_{\min }$, corresponds to 6 × 15 minutes, which is necessary to satisfy our detection criterion of six measurements with A > 3σ, where 15 minutes is the nominal cadence of the GBTDS. The higher cutoff, ${t}_{\max }$, is the proposed observation season duration for the GBTDS. This has still yet to be determined, though it is constrained to be between 60 and 72 days (Wilson et al. 2023). As such, we adopt the current nominal season duration of 72 days and display how shortening the season would influence our results with the dashed curve in Figure 2. Regardless of the choice for ${t}_{\max }$, the total expected number of events is then simply given by $N=6\times {\rm{\Gamma }}{t}_{\max }$, due to the six observation seasons.

Zoom In Zoom Out Reset image size

If the astrophysical backgrounds are fully known, i.e., no detected events can be attributed to an extended lens, then we determine the maximum sensitivity of Roman to a given population under the condition of null detection. In practice, it is often not possible to determine the identity of a lens on an event-by-event basis. Therefore, our results are not a detailed projection of expected limits, but rather a demonstration of Roman's maximum sensitivity to extended lenses. The Poisson 95% confidence limit under the assumption of null detection corresponds to N ≤ 3, which yields a maximum allowed value for f_DM. We perform the calculation of the event rate using the open-source Python package LensCalcPy ⁶ (Smyth & DeRocco 2024), a flexible tool for estimating microlensing rates that has been shown to yield results consistent with full population synthesis models (DeRocco et al. 2024).

Our results are displayed in Figure 2, where the horizontal axis shows the total mass of the lens and the vertical axis shows the mass fraction of DM contributed by a population of lenses at that mass. The colors of the curves correspond to various values of R₉₀ ranging from R₉₀ = 0.1 R_⊙ (effectively a pointlike lens) up to R₉₀ = 100 R_⊙ (a highly extended lens). The filled-in regions correspond to existing constraints from other microlensing surveys (Croon et al. 2020b; Cai et al. 2023).

We find that Roman will be most sensitive to lens masses of ∼10⁻⁷ M_⊙, corresponding to event durations on the timescale of hours. In general, for more diffuse lenses (larger R₉₀), the sensitivity weakens. This is particularly noticeable at low masses in which the magnification peak departs from the quasi-singular PSPL curve (Equation (6)). At high masses, the curves converge since the lenses are effectively pointlike, even for lenses with R₉₀ = 100 R_⊙. At M ≳ 10⁻³ M_⊙, the various curves all begin to scale as M^−1/2. In this regime, the event rate is linearly proportional to the number density of lenses and the cross section for a detectable event; since the number density for a monochromatic population scales as 1/M, while the Einstein radius, which effectively defines the cross section of microlensing events in the point-source limit, scales as ∼M^1/2, the resulting rate dependence scales as M^−1/2.

The sharp cutoff at low masses corresponds to the point at which finite-source effects suppress the peak magnification sufficiently to render events undetectable. The maximum magnification in this regime scales as ${A}_{\max }=\sqrt{1+4/{\rho }^{2}},\,(\rho \gtrsim 1)$. The cutoff in mass corresponding to this condition depends on both the compactness of the lens, which changes the effective R_E for a given total mass, and the sensitivity of the instrument, which determines A_T . For a pointlike lens, this cutoff translates to

$$\begin{eqnarray}&&{M}_{\mathrm{cut}}\approx \displaystyle \frac{({A}_{T}^{2}-1){D}_{L}{R}_{\star }^{2}{c}^{2}}{16{D}_{S}^{2}G},\end{eqnarray} \tag{ 16 }$$

where it is assumed D_L ≪ D_S . This corresponds to M_cut ≈ 8.7 × 10⁻¹⁰ M_⊙ for A_T = 1.05, D_L = 1 kpc, D_S = 8.5 kpc, and R_⋆ = R_⊙, which coincides well with the corresponding cutoff in the figure.

The dashed curve shows the sensitivity for R₉₀ = 0.1R_⊙ if the Roman observational season duration is shortened from 72 to 60 days (see Section 4). Since the peak sensitivity to extended lenses lies at ∼hour-long event durations, this change does not have an appreciable effect on our results, with sensitivity weakening by a factor of (72/60) ≈ 1.2 at low masses and (72/60)⁴ ≈ 2 at high masses due to the reduced phase space.

A wide variety of motivated theoretical models predict the formation of extended substructures in the dark sector. For the first time, we have comprehensively analyzed the sensitivity of the upcoming Roman Space Telescope to detecting such substructures in the Galactic Bulge Time Domain Survey. Our results, illustrated for three well-motivated but qualitatively distinct extended DM targets, show that Roman will be able to explore a broad parameter space at up to 4 orders of magnitude lower fractional DM mass contributions than existing constraints. Furthermore, the methodology outlined in this Letter can be readily extended to other sources. We note that the presence of caustic features in lightcurves will provide future opportunities to discriminate extended substructures from pointlike DM candidates, providing Roman with the potential to not only detect DM substructures but to probe their microphysical nature as well. We leave further exploration of this exciting possibility to future work. Furthermore, our work strongly motivates the GBTDS to adopt the maximal season duration subject to mission constraints (72 days); this result has immediate relevance for other targets in this mass range, such as free-floating planets and isolated black holes as well.

W.D. and N.S. acknowledge the support of DOE grant No. DE-SC0010107. V.T. acknowledges the support by the World Premier International Research Center Initiative (WPI), MEXT, Japan and the JSPS KAKENHI grant 23K13109.