Astrometric Weak Lensing with Gaia DR3 and Future Catalogs: Searches for Dark Matter Substructure

Small-scale dark matter structures lighter than a billion solar masses are an important probe of primordial density fluctuations and dark matter microphysics. Due to their lack of starlight emission, their only guaranteed signatures are gravitational in nature. We report on results of a search for astrometric weak lensing by compact dark matter subhalos in the Milky Way with Gaia DR3 data. Using a matched-filter analysis to look for correlated imprints of time-domain lensing on the proper motions of background stars in the Magellanic Clouds, we exclude order-unity substructure fractions in halos with masses $M_{l}$ between $10^{7} \, M_{\odot}$ and $10^{9} \, M_{\odot}$ and sizes of one parsec or smaller. We forecast that a similar approach based on proper accelerations across the entire sky with data from Gaia DR4 may be sensitive to substructure fractions of $f_{l} \gtrsim 10^{-3}$ in the much lower mass range of $10 \, M_{\odot} \lesssim M_{l} \lesssim 3 \times 10^{3} \, M_{\odot}$. We further propose an analogous technique for stacked star-star lensing events in the regime of large impact parameters. Our first implementation is not yet sufficiently sensitive but serves as a useful diagnostic and calibration tool; future data releases should enable average stellar mass measurements using this stacking method.


INTRODUCTION
The arrival of Gaia's second data release (DR2) heralded a new era in understanding the dynamics of the Milky Way (MW).The remarkable size and precision of this astrometric data set allowed scientists to pursue questions that had been previously unanswerable.As a result, we have learned of new stellar streams Mateu et al. (2018); Ibata et al. (2019Ibata et al. ( , 2021)); Li et al. (2022) and other galactic substructure Belokurov et al. (2018); Helmi et al. (2018); Koppelman et al. (2018); Myeong et al. (2018); Naidu et al. (2020).The extremely accurate astrometry of Gaia's Early Data Release 3 (EDR3) has revealed the 5 as/y global pattern in the proper motions of extragalactic sources imprinted by the acceleration of the Solar System through aberration of light Gaia Collaboration et al. (2021), and has produced the most accurate construction of the celestial reference frame, based on 1.6 million quasar sources Klioner et al. (2022).Additionally, the enhanced quality of the EDR3 astrometric data has allowed for a comprehensive study of the distance and kinematic properties of the Magellanic Clouds Luri et al. (2021) (geometrically anchoring the cosmic distance ladder), as well as for the exploration of the non-axisymmetric features of our own Galaxy Drimmel et al. (2022).
One exciting possibility with astrometric measurements of extraordinary precision is to explore the potential of time-domain astrometric lensing in the study of gravitating bodies.As these objects move across the sky, they can induce a time-dependent shift in the positions of background objects, yielding potentially observable effects.Van Tilburg et al. (2018), hereafter Paper I, introduced different classes of observables to look for nonluminous structures in the Milky Way, solely via the time-varying astrometric imprints on background luminous sources.Their detection strategies laid the groundwork for halometry, creating a flexible framework to characterize the properties of a broad spectrum of DM substructures within the Galactic halo.Capitalizing on one of these techniques to tease out subtle timedomain astrometric lensing effects, Mondino et al. (2020), hereafter Paper II, performed the first search for individual DM subhalos in the MW, by implementing a matched-filter template of local lensing distortions to the proper motions of Magellanic Cloud stars in the Gaia DR2 data.This proof-of-concept analysis constrained parsec-sized (or smaller) lenses in the 10 7 -10 8  ⊙ mass range.Although this technique exhibited sensitivity to more dilute subhalos, which could have otherwise escaped detection with photometric lensing probes, the results are currently statistics-limited.With increased observational periods and resolved systematics, the sensitivity of the search is expected to improve steeply over time, with the capability to deliver parametric leaps in reach with data from ongoing and planned astrometric surveys.
In this paper, we perform an improved search using proper motion templates for lensing by DM subhalos on Gaia DR3 data following Paper II, with the expected sensitivity enhancements borne out in this data release.We also outline extensions of this technique to parallax and acceleration templates; forecasts for the latter show exceptional promise for lighter subhalos, in the range 10  ⊙ -10 6  ⊙ .Since acceleration measurements are currently absent from the Gaia archive, we demonstrate a method to extract them by combining data from the DR2/DR3 releases; we will perform an acceleration template analysis on this catalog in follow-up work.While the proper motion template analysis presented in this paper will eventually hit a noise floor from intrinsic stellar motion dispersion -at least on the Magellanic Clouds (MCs) -the acceleration analysis is expected to be statistics-limited for the foreseeable future, with even steeper improvements expected after longer observation times.
We furthermore suggest that a stacked version of our template methods can look for collective star-star lensing in the astrometric, nontransient regime.To this end, we construct a catalog of optical doubles of foreground and background stars at substantially different line-of-sight distances, but at (accidentally) small angular separations.We execute the first star-star astrometric lensing analysis based on proper motion templates, and do not (yet) find evidence for this effect, in line with expectations.In the process, we identify several potential biases and systematics in Gaia's astrometric data for stellar pairs at very close separations.The inclusion of stellar accelerations in a future analysis of this type should yield a positive detection by Gaia DR4, potentially establishing an estimator for the average mass of the foreground stellar population, and serving as a calibrator for the DM searches.
This paper is organized as follows.In section 2, we investigate the full potential of searches for DM substructure using template observables.We give a review of their definitions in section 2.1, estimate current and future sensitivities in section 2.2, and present an updated limit from Gaia DR3 data in section 2.3.In section 3, we describe the new application of template observables to measure astrometric lensing in a collection of optical doubles.We present a test statistic that combines the three template observables in section 3.1, detail the sample selection and data processing in section 3.2, and summarize our findings in section 3.3.We conclude in section 4. Additionally, we discuss possible improvements and limitations of our test statistic for dark subhalo searches in appendix A, we compute for the first time the lensing-induced parallax and define the corresponding template observable in appendix B, and outline the extraction of accelerations from DR2/DR3 data with an estimate of their uncertainty in appendix C.
We publicly release our data analysis pipeline for the DM search and the star-star lensing analysis, allowing their use by other groups and on other astrometric data sets (not necessarily from Gaia).The code used to obtain the results of this study is available on GitHub and a link below each figure (ƭ) provides the Python code with which it was generated.Associated data are publicly available at this link.

DARK MATTER LENSING
Modern cosmology has cemented the role that Dark Matter (DM) plays in shaping the properties and growth of structures in the Universe.While the particle nature of DM is unknown, it is determined to be non-relativistic from  ≲ 10 7 and onward, and to have scaleinvariant adiabatic perturbations, which grow under the influence of gravity.These perturbations eventually form halos that merge hierarchically, yielding a nearly scale-invariant spectrum of structures with masses  halo ≳ 10 9  ⊙ .Clustering on smaller scales is largely unconstrained nonetheless.A variety of mechanisms -increased primordial perturbations, phase transitions, attractive interactions or dissipative dynamics -can lead to enhanced clustering at smaller scales far beyond what is predicted in ΛCDM.The main observational challenge is to detect these low-mass subhalos, which are entirely dark due to their lack of star formation or even baryonic content, indirectly through their gravitational interactions.Proposed detection methods (see Bechtol et al. (2022) for a recent overview) include strong lensing Vegetti et al. (2010); Hezaveh et al. (2016); Birrer et al. (2017), stellar streams Bonaca et al. (2019), lensing of gravitational waves Dai et al. (2018), pulsar timing arrays Dror et al. (2019); Ramani et al. (2020), and time-domain astrometric weak lensing (Paper I).In this work, we focus on the last method.
As dark subhalos move across the celestial sphere, the apparent locations of background luminous objects will undergo time-dependent shifts due to weak gravitational lensing.Specifically, the proper motion, acceleration, and parallax of the foreground subhalo will imprint correlated lensing-induced shifts to the proper motion, acceleration, and parallax of any background objects in the angular vicinity of the lens.While the lensing-induced motion of any individual background source is impossible to detect in the large-impact-parameter regime (at small impact parameters, searches for astrometric transients can be powerful Chen et al. (2023b)), the collective, correlated lensinginduced motions have a characteristic pattern that can be teased out statistically on multiple sources through a "template" analysis.These template observables are matched filters to the lensing effects on the background linear motion, linear acceleration, and parallax.
Previous works have considered microlensing effects over multiple background sources.For example, Gaudi & Bloom (2005) investigated the possibility of detecting outer solar system objects using astrometric microlensing induced by the lens parallax motion, while Di Stefano (2008) suggested to look for nearby and fast-moving lenses through a combination of photometric and astrometric lensing effects on multiple stars.Both works considered lenses that span large angles on the sky during the observation time, while here we focus on the regime of small relative changes in the lens-source impact parameters, where "template" observables can be defined.

Template observables
Template observables are designed to search for local, correlated astrometric weak lensing effects induced by extended dark lenses on multiple luminous sources.Proper motion and angular acceleration templates were first introduced in Paper I and an optimized version of the proper motion template was used in Paper II to search for DM-induced signals on Magellanic Clouds stars from Gaia's second data release (DR2).Since template observables are the main focus of this work, we summarize here the relevant existing results.
In the weak lensing regime, the lens-source angular separation is much larger than the Einstein radius of the lens and one image of the source is completely resolved.In this case, a gravitational lens  with mass   which appears to be near a background luminous source  induces an apparent angular deflection to the source's true position on the plane perpendicular to the line of sight, where  is Newton's gravitational constant,  the speed of light, and   the physical size of the lens.We also define the angular size   ≡   /  of the lens.The angular separation (impact parameter) between the source and the lens is denoted as   ≡   −   .The function   (  ) ≡   (  )/  is the normalized lens mass enclosed within the impact parameter, and we assume the source distance to be much larger than the lens distance   .
In general, the impact parameter changes over time due to the relative motion of the lens, the observer, and the source -  =   () -leading to a time-varying lensing angular deflection, which results into an apparent motion of the background source.When the relative change in impact parameter over the total observation time is small, |Δ  | ≪ |   |, the time-dependent lensing correction can be extracted through a series expansion of Eq. ( 1).The dominant contribution to the time-dependent impact parameter comes from the linear proper motion   of the lens with 3D physical velocity v and line-of-sight direction D , parametrized as: with  = 0 corresponding to the time of closest approach, at angular impact parameter   =  0  .For    ≪  0  over the observation time , the leading-order effect is then a lens-induced proper motion with a dipole-like profile: The next-to-leading order correction is a lens-induced angular acceleration with a quadrupole-like profile For a given lens mass distribution , the two-dimensional spatial profiles of the time-domain lensing distortions,  and , depend only on the position , the characteristic angular scale   , and the velocity direction μ .
The universal nature of the lensing distortion patterns from Eqs. ( 4) and ( 6) allows us to define the template observables T  and T  , which quantify the overlap between the measured proper motions {  } and accelerations {  }, respectively, of a field of stars with the expected distribution in the presence of a lensing signal.For a candidate lens at location   , with angular size   , and velocity direction μ , we define for stars at locations {  } in the surroundings of   with measured proper motion variance  , .We also introduce the convenient normalization factor N  .Analogous statistics T  and N  can be defined for lensing-induced accelerations.
If the mean motions of the background sources are zero (or subtracted) and the errors are Gaussian and uncorrelated, maximizing T is equivalent to maximizing the log-likelihood ratio of the hypothesis of a local signal (i.e. the presence of a lens in a given location with corresponding angular scale, and velocity direction) versus that of the null hypothesis.When searching for a population of dark lenses in front of a stellar target, the lens properties are usually unknown and need to be marginalized over, requiring a refined version of the T test statistic, such as the global test statistic R introduced in Paper II, utilized in the DR3 analysis in Sec.2.3, and discussed in appendix A. Nevertheless, the local template from Eq. ( 7) still captures the relevant features and its simple expression provides an easy analytic estimate of the local signal-to-noise ratio (SNR) and thus the best stellar targets in the search for a lensing signal from a dark subhalo.
By construction, the lens-induced angular acceleration in Eq. ( 5) is a subleading contribution to the lensing correction, and one could naively think that the proper motion signal always dominates.However, the sensitivity of each observable strongly depends on the properties of the measured luminous sources and is ultimately limited by the intrinsic variance of the stellar motions.We will compare the sensitivity of proper motion and acceleration templates to compact dark lenses and show their complementarity to cover different regions of the parameter space in section 2.2.
Other types of lensing corrections can arise from a lens motion which is different from the linear motion considered so far.For example, the source-lens impact parameter also changes over time due to parallax, as   () ≃   (cos , sin  ecl sin ), where   is the lens parallax,  ≃ 2/y is the orbital angular velocity, and  ecl is the ecliptic latitude of the lens -neglecting the background source's parallax and assuming for simplicity that the observer moves in a circular orbit.Expanding the angular shift in Eq. ( 1) for   ≪  0  , the lens-induced parallax of the luminous source can be obtained and the corresponding template observable can be defined.The derivation of this anomalous parallax contribution and estimates for the sensitivity of the parallax template observable, compared to the proper motion template, are presented in appendix B. When the astrometric measurements are statistically limited, T  is typically more sensitive than T  , as the total distance travelled by a dark lens with velocity   ∼ 10 −3  is usually much larger than the parallax displacement of one AU, especially over multi-year observations.However, when the proper motion dispersion is limited by the intrinsic stellar motion, T  and T  can have comparable sensitivities (see appendix B).
For other types of lenses, such as outer Solar System planets, which move slowly with respect to the observer, the anomalous parallax is the dominant effect Gaudi & Bloom (2005); Paper I. We conclude this section by highlighting the difference between the lensing regime in which the template approximation is valid and where a different treatment is needed.As pointed out earlier, template observables correctly capture the astrometric lensing corrections when the arc on the sky spanned by the lens motion during the survey is smaller than the impact parameter   .For lenses with a finite characteristic size   =   /  , the SNR for the lensing-induced proper motion is maximized at   ∼   , unless the lens has a very cuspy inner profile (see Sec. 2 of Paper I) or is effectively point-like, i.e. when   is smaller than the typical star-star separation, which for Gaia is ≳ 0.7 ′′ Fabricius et al. (2021).Therefore, the template searches in this work are designed for lenses with   ≳ 0.003 pc ( ⊥ /10 −3 )(/10 y) for the proper motion and acceleration effects, and   > AU for anomalous parallax.When this condition is not satisfied, e.g. for very compact, cuspy, or fast-moving lenses, the astrometric lensing effect manifests as a transient perturbation to the trajectory of individual (or multiple) stars that can be searched for using other observables, referred to as mono-blips (or multi-blips) in Paper I -see Chen et al. (2023b) for a projection of the mono-blip observable in (mock) Gaia DR4 data.

Sensitivity estimates and projections
Proper motion and acceleration template searches in current and upcoming Gaia data can be sensitive to signals from Galactic compact dark lenses, thus probing parts of previously unexplored parameter space of dark matter substructure.For the sensitivity forecasts, we briefly review here the estimates of the local SNR for the template observables and refer the reader to Sec. 4.2 of Paper I for a more extensive derivation.
The observed stellar population is assumed to have zero or subtracted background motion and uncorrelated Gaussian noise, such that ⟨T ⟩ noise = 0 and ⟨T 2 ⟩ noise = N 2 .On a field of background stars with angular number density Σ 0 , one expects that up to an O (1) numerical factor that depends on the characteristic lens density profile.For a template that perfectly matches the true lens properties, ⟨T ⟩ signal =  ,  N 2 , where Therefore, the SNR for the template observable is This expression makes evident the fact that the proper motion and acceleration observables are the two leading terms in the Taylor expansion in   /  =  ⊥ /  of the time-domain lensing signal.The Projected sensitivity with 10 y of Gaia Figure 1.Sensitivity projections for proper motion and angular acceleration templates in the parameter space of lens fractional abundance   =   / DM versus mass   .All curves correspond to local SNR T = 1 (optimistically ignoring the look-elsewhere effect), with the signal-to-noise ratio given in Eq. ( 9), taking  ⊥ = 208 km/s and the smallest lens distance from Eq. ( 10).
The proper motion template (blue contours) is evaluated for the minimum lens size given in Eq. ( 11) (solid) and a fixed value of   = 1 pc (dash-dotted), using the LMC as stellar target (Σ 0 = 5 × 10 8 rad −2 , ΔΩ = 0.02 rad 2 ,   = 0.2 mas/y,  stars = 50 kpc).The acceleration template (red contours) is evaluated for the minimum lens size using Galactic Disk stars (Σ 0 = 4.5 × 10 9 rad −2 , ΔΩ = 0.2 rad 2 ,  stars = 3 kpc), taking   7 = 4.5 as/y 2 (solid) and   Δ = 300 as/y 2 (dotted) -assuming per-epoch positional accuracy of 200 as.ƭ SNR in Eq. ( 9) is dominated by the most massive, fastest, nearby, and compact lens in the parameter combinations shown.The variation in  ⊥ =     is independent of the other lens properties and given by the DM virial velocity dispersion in the Galaxy ⟨ ⊥ ⟩ = √︁ /2 ,DM , with  ,DM ≃ 166 km/s.We assume the lens population is distributed according to the Galactic density profile   =    DM , where   is the (constant) fraction in DM substructure in halos of mass   .The closest lens in front of a stellar target with angular number density Σ 0 is expected at a distance of where  DM (and thus   ) is approximated to be constant up to  ,min .
The smallest lens size  ,min allowed in the template regime is fixed by the lens displacement over the observational time  ⊥  =      (see section 2.1).However, the sensitivity will not improve further for lenses that are effectively point-like, i.e. when their angular size is smaller than the typical angular separation ∼ 1/ √ Σ 0 of background sources.Therefore, In Fig. 1 we show SNR = 1 sensitivity projections for proper motion and angular acceleration templates for 10 years of Gaia observations, using the above values of  ,min and  ,min in Eq. ( 9).The sensitivity depends on the proper motion and acceleration dispersions,   and   , of the observed stellar populations.If these were only instrumentally limited,  obs observations with individual position uncertainty    would result in   ∼    /( √  obs ) and   ∼    /(2 √  obs ); we estimate the numerical factors of this parametric dependence in appendix C. In this case, SNR T  /SNR T  =  ⊥ /  < 1, i.e. the proper motion template would always perform better, as expected, since it is the leading-order effect when  ⊥  <   .The instrumental precision for most of the stars observed by Gaia is already below the intrinsic proper motion dispersion, or will be so by the end of the mission.On the other hand, a measurement of angular accelerations by Gaia would only be statistically limited because the intrinsic accelerations -from the Galactic potential, wide binary companions, or exoplanets -are typically far below the survey's sensitivity on an individual distant star, and are uncorrelated among nearby stars (a relatively small fraction of nearby and/or bright stars do have detectable intrinsic accelerations from binary companions Halbwachs et al. (2023); Penoyre et al. (2022); Ranalli et al. (2018)).Consequently, we see in Fig. 1 that acceleration templates offer the best prospects for measuring a compact DM subhalo lensing signal throughout a large portion of the parameter space, particularly for masses below ∼ 10 6  ⊙ .
For each observable, we consider the best stellar target, i.e. the one that maximizes the figure of merit √ Σ 0 / ,  .The best performance is obtained by a T  search on the Large Magellanic Cloud (LMC) and a T  search on Galactic Disk stars.An additional advantage of using angular accelerations is the possibility of using a larger star sample, spanning a wider portion of the sky, while minimizing intrinsic proper motion dispersion requires choosing stars that are further away (but still in a dense field).The SNR scaling obtained above breaks down at large lens masses, when Poisson fluctuations of the number of lenses in the field of view (FOV) become important.For our sensitivity projections of Fig. 1, we require more than 3 lenses in front of the stellar target.
Currently, angular accelerations for sources observed by Gaia are not directly available.The optimal way to measure the acceleration is to include it as an additional parameter to the astrometric fit.However, before the full time series of Gaia's observations become available with DR4, an indirect measurement of accelerations can be obtained by combining position and proper motion measurements at different times from the DR2 and DR3 catalogs.In appendix C, we outline a method to construct such an acceleration catalog together with an estimate of the expected statistical uncertainty.The implementation of this method and a careful investigation of the systematics involved is left to future work Chen et al. (2023a).Applying a template search for DM lensing effects on this catalog would already produce interesting results, as shown in Fig. 1, where we refer to the expected uncertainty on the derived acceleration as   Δ , while   7 refers to the statistically-limited uncertainty from a 7-parameter astrometric fit.The possible applications of an acceleration catalog for Gaia's sources of course go beyond DM searches with astrometric weak lensing and include measuring collective star-star lensing effects (see section 3), mapping out the Milky Way potential Buschmann et al. (2021), and searching for ultra-low frequency gravitational waves Pyne et al. (1996); Book & Flanagan (2011); Klioner (2018).

Search for Dark Matter Substructure with Gaia DR3
The proper motion template observable described in the previous sections has already been successfully applied to Gaia DR2 data to search for lensing signals induced by galactic DM subhalos on MCs stars in Paper II.Here we repeat the same analysis with the improved astrometry of Gaia DR3 Lindegren et al. (2021); Luri et al. (2021).We closely follow the procedure described in Paper II; we only briefly review it here and refer the reader to the original paper for more details.
Secondly, the data are cleaned in two steps: (i) The background mean proper motion field computed with a Gaussian distance kernel of radius 0.1 • is subtracted from the stars' proper motion; the smoothed angular number density map is computed with a Gaussian distance kernel of radius 0.1 • and overdense pixels of about 0.014 • are removed if their density is 2.5 times larger than the average to reduce contamination from star clusters; proper motion outliers at more than 5 are removed and the following additional cuts on the parallax and the quality of the astrometric fit are imposed Lindegren et al. (2021): parallax/parallax_error < 2, ruwe < 1.4, ipd_gof_harmonic_amplitude < 0.4, ipd_frac_multi_peak < 40, ipd_frac_odd_win < 40.
(ii) The background mean proper motion field computed with a Gaussian distance kernel of radius 0.06 • is subtracted and proper motion outliers at more than 3 are removed, repeating the procedure 3 times; the effective proper motion dispersion as a function of the stars' G magnitude and distance from the center of each cloud are computed using G magnitude bins of size 0.1 and radial bins of size 1 • to group the stars.
The optimized global test statistic R from Eq. ( 9) in Paper II, obtained by maximizing the likelihood ratio of the signal hypothesis over background, is evaluated on the end products of the procedure described in (ii).The same clean data are also used as input for the data-driven noise + signal simulations as described in Paper II.The processing step in (ii) is performed on the simulated data, before the evaluation of the test statistic.For each simulation, we fix the lens population parameters {  ,   ,   ≡   / DM }, where  DM is the Milky Way DM halo density profile, which we model as an NFW profile  DM = 4  /[(/  )(1 + /  ) 2 ] with scale radius   = 18 kpc and density at the scale radius of   = 0.003  ⊙ /pc 3 . 1The distribution of the test statistic from 150 T µ limits and projections for compact lenses (r l = 1 pc) The solid blue (red) line shows the 90% (50%) CL limit from the analysis on Gaia DR3 data, while the dashdotted lines correspond to the previously obtained results from Paper II.In light gray we show the SNR T = 1 curves, where SNR T is given in Eq. ( 9),   =  ,min from Eq. ( 10),   = 1 pc, Σ 0 = 5 × 10 8 rad −2 , ΔΩ = 0.02 rad 2 , and  ⊥ = 208 km/s.The proper motion dispersion   is fixed to 1.75 mas/y to anchor the SNR T = 1 line to the 50% CL DR2 limit, scaled by a factor of 0.5 to obtain the statistics-limited optimistic DR3 improvement (see Eq. ( 12)), and fixed to an intrinsic dispersion of 0.2 mas/y for the EOM projections.Differently from Fig. 1, here an estimate of the look-elsewhere effect is included in the analytic projections by anchoring the DR2 line to the observed limit and taking into account the LMC average effective dispersion of 1.54 mas/y in DR2 (see Fig. 4) ƭ simulations per point is then compared to the observed value in the data and the parameter space point is excluded at 90%(50%) CL if R data (  ,   ,   ) < R 90% (50%) sim (  ,   ,   ).Following the procedure described above, we can exclude the presence of dark lenses in front of the MCs.The results are shown in the parameter space of fractional lens abundance   versus lens mass   , for compact lenses (  = 1 pc) in Fig. 2 and point-like lenses (  = 10 −3 pc) in Fig. 3.We compare the updated limits obtained with the analysis on Gaia DR3 stars with the ones obtained with the analysis on Gaia DR2.The improvement in the constraints comes from the expected improvement on the stellar proper motion measurements due to the longer observation time, As shown in Fig. 4, this improvement is borne out in Gaia's reported error  ,Gaia ≡ ⟨ 2 , ⟩ 1/2 for all the stars in the LMC.However, the relevant quantity for the proper motion template evaluation is the effective proper motion dispersion  ,eff ≡ ⟨(  −  mean ) 2 ⟩ 1/2 observed in the stellar population, which improves by a factor of 2 only for stars with 19 ≲  ≲ 20.The dispersion at lower and higher G magnitude values appears not to be statistics-limited and therefore has a somewhat less favorable scaling with time.
In Figs. 2 and 3, the observed DR2 and DR3 limits on dark lenses upper limits on   that are larger than unity (but less than 2 or 3) in our conservative model can be physical.The thick lines correspond to Gaia DR3 data, while the thin lines for  ,eff and  ,Gaia correspond to Gaia DR2 data rescaled by a factor of 0.5 that corresponds to the expected statistical improvement with time (see Eq. ( 12)).ƭ are compared with the optimistic improvement expected for statisticslimited proper motion variance, using the SNR scaling derived in section 2.2.The expected statistic-limited improvement on the fractional abundance of the dark lenses is fixed   (Fig. 2) ,min (Fig. 3) The projected sensitivity for Gaia's end of the mission (EOM) results is also shown, assuming that the effective dispersion will be completely dominated by the MCs' intrinsic proper motion  ,eff =  ,intrinsic ≃ 0.2 mas/y.The global template test statistic R used in the analysis above relies on the strongest signal produced by an individual lens in front of the stellar target.In appendix A we discuss a possible generalization that includes multiple lenses and show how it is not expected to improve over current results -with our conservative treatment for noise modeling -because the signal decreases faster than the noise for subleading lenses (due to the look-elsewhere effect).Further investigation of an optimal test statistic that leverages on multiple lenses while also not being hurt by strong look-elsewhere effects go beyond the scope of this work.
Several other gravitational probes of massive compact objects which make up all or a fraction of the DM abundance have been investigated in the literature.The most well-studied candidates are primordial black holes (PBHs), which are strongly constrained in the mass range 10 6 -10 9 M ⊙ (see for example Carr et al. (2021); Bird et al. (2023)).Some of the existing bounds that would rule out all the parameter space of Fig. 3 for point-like lense only apply specifically to PBHs; these are effects that rely on accretion, such as CMB spectral distortions and anisotropies Ricotti et al. (2008); Ali-Haïmoud & Kamionkowski (2017) and direct X-ray observations Inoue & Kusenko (2017).Dynamical effects on the MW Xu & Ostriker (1994), globular clusters Moore (1993); Brandt (2016), dwarf galaxies Lu et al. (2021); Takhistov et al. (2022) 2 , and wide binaries Yoo et al. (2004); Quinn et al. (2009) would probably also be induced by the most compact dark lenses considered here; however in order to derive robust constraints from these effects, dedicated galactic simulations and direct comparison to data are required to obtain reliable constraints from dynamical friction and heating effects induced by compact (but finite size) subhalos Arvanitaki et al. (2020).PBHs constraints based on the induced formation of cosmic structure incompatible with LSS Carr & Silk (2018) and Lyman- forest observations Murgia et al. (2019) are model dependent, as they are sensitive to the cosmic history and formation time of the subhalos.Finally, traditional photometric microlensing surveys rule out ultracompact lenses with masses below about 10 3  ⊙ Blaineau et al. (2022), which are lighter than the ones considered here, and furthermore require much smaller subhalo scale radii (smaller than their Einstein radii) for the strong-lensing signal to be unsuppressed.

STAR-STAR LENSING
The previous section expounded on a blind search for astrometric lensing effects from dark subhalos, lenses whose properties are, at best, only known at the population level.We now consider a different kind of lens target: isolated stars in the Gaia catalog that, by chance, happen to be angularly close to other, more distant stars.We will refer to such resolvable optical doubles as pairs of foreground and background stars.In this case, most properties of the lens -its location, velocity, and distance -are known, except its mass.The measurement of lensing corrections induced in the motion of the background stars can therefore be used for directly inferring the mass of the foreground (lens) star.The possibility of a mass measurement for individual stellar lenses through astrometric microlensing has been investigated in Boden et al. (1998) 2022) and a precision better than 15% in the mass measurement is expected for a dozen events Klüter et al. (2020).
We suggest a complementary approach, which aims at detecting collective weak lensing distortions for a large number of stellar pairs at wider angular separations, where the lensing event is not a transient on typical survey time scales.The cumulative effect from multiple star-star lensing corrections can be captured using our template observables and searched for in existing and upcoming astrometric data sets, including those of Gaia.
To gain intuition on the sensitivity of current and future data, we can estimate the SNR within some approximations.The signal is dominated by optical doubles at the smallest observable angular separation  min , and with the lens at the smallest possible distance  ,min .However, the template regime only applies for  ,min  min >  ⊥ .The best-case scenario is thus at the saturation of this lower bound.We can stack all of the star-star lensing events and compute collective test statistics T ★  and T ★  to arrive at SNRs analogous to those in Eq. ( 9): for the proper motion and the acceleration templates, respectively.In the above, we have assumed solar-mass lenses with a constant number density   up to  ,min -about a few hundred pc -and a uniform angular number density Σ 0 of distant background stars.The smallest star-star separation is limited by Gaia's completeness to be  min ∼ 0.7 ′′ Fabricius et al. (2021).
From Eq. ( 14), we see that a star-star lensing detection using proper motion templates on Gaia DR3 data is challenging.If the proper motion dispersion of background stars were statistically limited, it would scale as   ∝  −3/2 and the sensitivity would improve linearly with time.However, from the analysis presented in the next sections, we find that   is most likely already limited by the intrinsic stellar motion and will not decrease with time.On the other hand, a near-future detection using angular acceleration templates looks promising.In Eq. ( 15), we used the dispersion expected for accelerations inferred from the combination of Gaia DR2 and DR3 observations (see appendix C).However, a dedicated 7-parameter astrometric fit including angular accelerations could lead to: where we used a per-epoch positional accuracy of 250 as, reflecting Gaia's astrometric performances on stars with magnitude  = 17; for fainter stars with  = 19, the sensitivity worsens by about a factor of 3.6.Effectively, we predict that Gaia DR4 should produce a high-SNR measurement of collective star-star lensing.Given the scaling with observation time, our approach may lead to meaningful mass measurements at the population level (for isolated stars) by the end of the Gaia mission ( ∼ 10 y) if instrumental systematics are kept at bay.
In the following, we first describe how template observables can be easily adapted to search for star-star time-domain gravitational lensing signals.We then present an analysis applying proper motion templates on a sample of foreground-background stars selected from the Gaia DR3 catalog.As expected, a detection is not yet possible, but our implementation of the analysis pipeline provides insights on the sensitivity to the signal beyond the estimates of Eq. ( 14), identifies parts of phase space with systematic effects, and sets the stage for a future detection using stellar accelerations.The technique would enable a measurement of a weighted population average mass for the foreground stars, instead of individual masses.With sufficient sensitivity, the star-lenses could be grouped in classes of stars with similar properties, as inferred from photometric and spectroscopic information, to obtain average mass measurements of more homogeneous populations.Finally, template analyses of star-star lensing serve as a valuable calibration for searches of dark lenses that rely on similar techniques (section 2).

Combined template observables
Consider a pair of accidentally-close, but resolvable, stars on the sky at different line-of-sight distances: the foreground star, , will act as a gravitational lens on the background star, , inducing the same timedomain astrometric lensing effects as described in section 2.1.In this case, the lens is point-like, with   () = 1 in Eqs. ( 1)-( 6), and known distance, position and velocity.For each foreground-background pair, the signal is too small to be detected individually, but we can add up the contributions from many such pairs  with angular separation   <  max , where the maximum separation taken to be  max = 3 ′′3 .We define the star-star proper motion template observable and its normalization as: where    and Σ  , are the proper motion and proper motion covariance matrix of the background star in the pair .Notice that a single foreground star could be lensing multiple background stars and so appear in multiple pairs in Eq. ( 18).We define the lensing distortion as the one produced by a lens of one solar mass, where   =   −   ,   =   −   , and   () denotes the distance to the source (lens).Similar template observables can be defined for the acceleration and parallax, T ★  and T ★  .The three templates can be combined in a single test statistic The expected signal-to-noise ratio for the star-star lensing signal is SNR T ★ = ⟨  ⟩N ★ all , where ⟨  ⟩ is an average lens mass in solar mass units.An estimator for the average stellar mass is thus: If an estimator of the foreground star mass is available, i.e. as estimated from photometry or spectroscopy, it can be included as an additional weight in the lens-induced distortion of Eqs. ( 19) 4 .The Gaia catalog contains several millions of accidentally-close foreground-background pairs as defined in this work.Below, we describe our selection of the data sample and present an analysis with the proper motion template observable.We do not include the angular acceleration, due to the lack of acceleration measurements for now, nor the parallax template, due to the challenges in subtracting the background (see section 3.3 for further details).Once a catalog of angular accelerations is built, it will be straightforward to extend our analysis to include the acceleration template.

Data sample and background subtraction
To compute the above star-star lensing template, we create a catalog of optical doubles within Gaia DR35 .Any of these accidental doubles consists of a foreground star -the lens  -and a background starthe source .Specifically, we searched for all stellar pairs satisfying the two conditions |   | <  max and In other words, we selected stars within  max = 3 ′′ from each other on the celestial sphere, but at significantly different line-of-sight distance (at significance of   = 2), thereby reducing the contamination from binary stars and incorrectly-tagged pairs.
With those criteria, we identified about 61 million optical doubles.In the following sections, we describe a procedure whereby we further select stellar pairs within the original sample to ensure quality of the astrometric fit, reduce contamination from incorrectlyclassified pairs, subtract the mean proper motion from the population of background stars, and remove proper motion outliers.After applying these additional selection cuts, we retain a clean sample of about 11.4 million optical doubles that will be used to evaluate the lensing template.A sky density map of the selected stellar pairs is shown in Fig. 5.

Quality and distance cuts
To ensure a high quality of the astrometric data used in our analysis, we require the goodness-of-fit statistic RUWE to be < 1.4 Lindegren et al. (2021).Furthermore, since for most of Gaia's stars the inverse parallax is a poor estimate of the distance due to large fractional uncertainties, we make use of the probabilistic stellar distances estimated in Bailer-Jones et al. (2021) to further ensure that foreground sources are indeed in front of their background counterparts -in addition to the parallax selection described above.Bailer-Jones et al. ( 2021) provides two types of distances: geometric, combining parallax together with a direction-dependent prior on distance, and photogeometric, which additionally uses color and apparent magnitude measurements.We use their geometric distances, since accidental doubles occur in dense regions of the sky so their color might be incorrectly estimated.We impose the condition , where   () is the posterior median of the geometric distance of the lens (source), and  ,86th and  ,14th Stellar pair density denote the difference between the posterior median and the 86th and 14th percentiles for the lens and the source, respectively.The quality and distance cuts heavily reduce the size of our original sample, leaving about 25% of the original stellar pairs6 .

Background motion subtraction and outliers removal
In order to avoid spurious signals contributing to the the proper motion template of Eq. ( 18), we should make sure that T ★  vanishes in the absence of a real lensing signal.In the limit of background stars uniformly distributed around the lenses, we expect that ⟨T ★  ⟩ noise = 0, owing to the parity symmetry of the dipole pattern of the lens-induced distortion Δ  of Eq. ( 19).This remains true even with a preferred lens velocity direction, ⟨ μ ⟩ ≠ 0, and a nonzero background star motion, ⟨   ⟩ ≠ 0, with averages taken over all pairs .However, we find that the background stars are not distributed isotropically, but more background stars are observed in specific directions, particularly at small lens-source angular separation, as shown in the upper row of Fig. 6 in equatorial, ecliptic, and galactic coordinates.The anisotropy does not seem to be (strongly) related to the foreground star's proper motion direction, as shown in the bottom left panel of Fig. 6, while the bottom right panel shows the distribution of the lenses' proper motion.We suspect that the most likely (and plausible) candidate for this anisotropy in number density can be attributed to Gaia's scanning law pattern.
In order to mitigate the effect of the observed anisotropy, we subtract the mean motion of the background stars, which ensures that ⟨T ★  ⟩ noise = 0 irrespective of the background stars' distribution.To do so, we bin the background stars in a 4-dimensional histogram based on their on-sky position, G magnitude, angular separation from their foreground counterpart, and inferred geometric distance.For the spatial pixelation, we use the nested HEALPix7 scheme at level 4 Zonca et al. (2019); Górski et al. (2005), corresponding to a resolution of approximately 3.7 deg.In G magnitude and angular separation, the bin sizes are chosen to be 1 and 0.3 ′′ , respectively.Finally, due to the large uncertainties on the stars' distances, we use a probabilistic bin assignment in 4 logarithmically-spaced bins between 1 and 10 kpc (plus 2 additional bins for stars falling below and above this range); since Bailer-Jones et al. ( 2021) do not provide the full posterior distribution, we crudely model it as a Gaussian centered on the median, and with standard deviation ( ,86th +  ,14th )/2.To ensure that we can reliably estimate the statistics in each bin, we only retain stars with at least 80% of their probability support in bins with a threshold count of 30 stars.The removal of stars falling in "sparse" bins results in the reduction of the sample with spatial features due to the HEALPix pixelation, as shown in the right panel of Fig. 5.In each bin, we compute the proper motion mean μ and the effective proper motion dispersion covariance matrix Σ  .Outlier stars with proper motion  too far from the mean μ -those with Since the lensing signal is not expected to induce any such large deviation on individual stars, this cleaning ensures the removal of spurious fast-moving stars while retaining most of the signal.This procedure is iteratively repeated 10 times, with only about 0.01% of outliers identified in the last iteration.
The final sample with the subtracted proper motion mean and the effective proper motion covariance can be used to compute the template of Eq. ( 18).Altogether, about 77% of the remaining stellar pairs survive the aforementioned cleaning procedure.We compare the estimated effective proper motion dispersion to the instrumental errors reported by Gaia as a function of the background stars' G magnitude and angular separation from their foreground counterpart in Fig. 7.The observed dispersion of the stellar population is much larger than the error reported by Gaia and nearly independent of G magnitude, showing that for the majority of the background stars in our sample the instrumental precision is well below the intrinsic proper motion dispersion, as expected.The right panel of Fig. 7 shows that background stars which are closer to their foreground counterpart are more poorly measured and have larger dispersion, since they are observed in crowded regions; on the other hand, for angular separations ≳ 1.5 ′′ the distributions are nearly flat.In the next section, we report results of the proper motion template analysis on the clean sample of optical double stars.

Unweighted proper motion template
The computation of the proper motion template from Eq. ( 18) on the cleaned sample of about 11.7 million stellar pairs leads to the results For all stellar pairs, the lens is at the origin, and the black arrow in the left panel points in the direction of the average lens velocity over all pairs.The number of background stars decreases significantly below ∼ 1.5 ′′ and it is not uniformly distributed.(Bottom left) Number density of background stars around their lens counterpart as a function of the parallel and perpendicular components of the angular impact parameter with respect to the lens velocity direction where, for all stellar pairs, the lens is at the origin; in this plane, all the lens velocities are directed along the positive horizontal axis and the anisotropy is not visible, suggesting that it is not related to the lens velocity direction.(Bottom right) Distribution of the proper motions of the star-lenses.ƭ   1. Proper motion template evaluated on the sample of 11.7 million cleaned Gaia DR3 optical doubles.The first row corresponds to the signal channel, where the observed background stars' proper motion are matched filtered with the lensing-induced dipole profile, while the other rows correspond to control channels, where no signal is expected.In each case, the template, its normalization (expected SNR) and their ratio (observed SNR) are reported.For the signal channel, the measured average lens mass is also indicated.ƭ shown in Tab. 1.All stellar lenses are equally weighted with unit Solar mass.In addition to the signal channel, which computes the overlap of the background stars' proper motion with the lens-induced dipole profile, we consider three control channels: a dipole profile rotated by 90 • , a monopole profile Δ  ∝ β , and a quadrupole profile Δ  ∝ 2 μ β • μ + β 1 − 4( β • μ ) 2 , all with the same radial scaling as the signal dipole channel.The expected SNR is similar in all channels and below unity, as expected, about a factor of 4 smaller than our naive estimate of Eq. ( 14).The dipole signal channel shows an excess of 3.75.However, at this point, we caution against interpreting such measurement as a positive detection -and attribute it to a potentially-resolvable systematic -for the following reasons.Taken at face value, the measured value of T ★  corresponds to an average lens mass of 11  ⊙ , which is incredulously high given Gaia's estimated astrophysical parameters, as discussed below.Moreover, we further test the signal hypothesis by evaluating the template on subsamples of the entire catalog of optical doubles (see Fig. 8), in which the pairs are binned according to the foreground and background distances (left and middle panel) and their angular impact parameter (right panel).The expected scaling of N ★  with   and   is borne out in the data.The peaks in T ★  /N ★  occur in bins with the smallest SNR, contrary to the expectation from a real signal.We conclude that the large T ★  value reported in Tab. 1 is most likely due to a systematic effect that has not been properly removed by our cleaning procedure and that is more pronounced in bins with poor statistics.The middle panel of Fig. 8 reveals how the anomalously-large T ★  value is mostly due to optical doubles where the lens is further than about 300 pc, where the SNR is below 0.07.This justifies removing pairs with   > 300 pc from our sample, reducing the size to about 0.814 million stellar pairs that retain most of the sensitivity.As shown in Tab. 2, the value of N ★  is almost unchanged and the measured T ★  does not show a significant fluctuation, as expected for the (statistical-)background-only hypothesis.The removal of distant lenses is also justified by the sensitivity estimates discussed above (see Eq. ( 14)), which showed how the SNR is maximized by lenses at   ≃   / min ≃ 270 pc, for   ≃ 10 −3 ,  = 3 y, and  min = 0.7 ′′ .
Recently McGill et al. (2020) pointed out that a significant fraction of microlensing events predicted in Gaia DR2 were spurious events, where the background source candidate turned out to be either a duplicate detection or a binary companion of the lens.The issue was found to mostly affect bright sources with no 5D astrometric solution.We believe it is unlikely that a similar spurious contamination affects our sample of lens-source pairs for the following reasons.First, by construction, we only select background sources with proper motion measurement and we additionally require high quality of the astro- metric solution, together with significant difference in the measured parallaxes and distances in each pair (see Sec. 3.2.1).These requirements ensure that the majority of the pairs are genuine sources that are widely separated along the line of sight.Nevertheless, we performed sanity checks on the distribution of source-lens difference in magnitude and color, following the procedure adopted by McGill et al. (2020), and found no anomalous clustering of photometrically identical pairs, particularly for the not-super-faint sources that were identified as problematic in the microlensing events.Based on these results, we are confident that the sample of about 814 thousand pairs that drives the sensitivity to the astrometric lensing signal is reliable.

Alternative templates
Gaia DR3 provides astrophysical parameters from BP/RP and RVS spectra for nearly 0.5 million sources Creevey et al. (2023), Fouesneau et al. (2023), Delchambre et al. (2023), including an inferred mass for 128 million stars.In principle, the mass estimate for the foreground stars in our sample -or a proxy for the mass, such as the luminosity or the effective temperature -can be used as an additional weight for the lensing correction in the template of Eq. ( 18).From our sample of 11.7 million clean pairs, around 2 million lenses have an estimated mass.The SNR of this subsample (weighting all the lenses equally with a mass of  ⊙ ) is reduced by a factor of 0.13 with respect to the full sample, which is 3 times worse than the naive scaling with √︁  pairs .This is likely due to the fact that inferring astrophysical parameters is harder in crowded regions, where the lensing signal is maximized.If we include the inferred mass as weights, we get N ★  = 0.03 instead of N ★  = 0.046 with unit weights, a reduction that is expected given the lenses' median mass of 0.89  ⊙ .
Given the reduced sensitivity, we refrain from using astrophysical parameters at this point.This is further justified by the well known issue of Gaia's photometry in crowded regions: the mean source photometry is more likely to be blended with other sources in all or some of the observation epochs, possibly having a significant effect on the estimated photometry Riello et al. (2021).In future data releases, the effect of crowding on the spectra will be mitigated, opening up the possibility of reliably using astrophysical parameters for stars with an accidentally-close pair that can induce a lensing signal.
The time-varying lensing effects induced by the foreground stars could also manifest as an induced parallax in the background stars that can be captured by an appropriate template observable, as discussed in Sec.2.1 and worked out in detail in appendix B. We did not include it in our current analysis due to the challenges encountered when modeling and subtracting the background, which is crucial to remove large systematic errors.Differently from the stellar proper motions, the distribution of observed parallaxes is very different from a Gaussian distribution and highly skewed.This makes it difficult to subtract the mean background parallax and model the effective dispersion, which are needed in order to evaluate the significance of the observed value of T ★  /N ★  .Since the parallax template is expected to be less sensitive than the proper motion template (or at most comparable) and the signal is currently below the detection threshold, we leave further studies of the parallax background to future work.

CONCLUSIONS
The precision astrometry provided by Gaia has opened several windows for understanding the dark matter that makes up the Milky Way halo.Aside from singular astrometric microlensing events, the precision of Gaia allows one to consider aggregate statistical observables as well.
One class of such observables is the matched-filter or template observable to look for correlated motions or accelerations of background stars induced by foreground lenses (Paper I).Here, while individual stars may not provide a signal above noise levels, fitting a combination of luminous sources can yield a signal, or provide constraints.The first search using these techniques was performed by Paper II, finding no signal as expected, even assuming all of the dark matter was in the form of such objects.
In this work, we have extended those analyses using DR3, placing limits on dense subhalos with scale radii of 1 pc and 10 −3 pc, with constraints on order-unity fractions of the dark matter in such objects down to halo masses of 2 × 10 7  ⊙ and 10 6  ⊙ , respectively.Such limits can yet improve by an order of magnitude or more in future Gaia data releases before intrinsic stellar proper motion dispersion limits further improvement.
We have also applied the velocity template to star-star pairs.While astrometric lensing by individual stars is too weak to generate a detectable signal on a companion star, it is possible that by aggregating the signals from many stars, we may see this effect.We have created a catalog of 61 million optical doubles, from which we have extracted a cleaned subsample of 11.7 million, with 814 thousand pairs in our expected signal region.As expected, we find no excess in our main signal region, while we find a small excess correlation (inconsistent with a real signal) in the full cleaned sample.The origin of the latter spurious fluctuation remains unknown, but intrinsic correlation of stellar motions may play some role.While we do not expect velocity templates to provide a measurement of masses given the intrinsic dispersion of stellar motions, acceleration templates should be able to detect a signal and perhaps serve as a tool to measure properties of stars.
Going beyond velocity templates, we have further extended the matched-filter approach to accelerations and parallax.While accelerations were previously noted as a possible template in Paper I, no sensitivity analyses were performed.We have forecast the sensitivity of acceleration templates and found they can be a powerful probe in the range of 1  ⊙ − 10 7  ⊙ .While Gaia does not provide acceleration data at this time, one can extract accelerations analytically using a combination of DR2 and DR3 data.We have forecast the sensitivity for this approach and find that even this could provide interesting limits on dark halos in the mass range 10 3  ⊙ − 10 7  ⊙ .
We have also studied for the first time the time-dependent lensing effect arising from the parallax motion of the intervening lens.The parallax motion of the lens can induce an anomalous parallax motion on a lensed object, yielding motions that have annual periodicity, but can, in general, be quite different from true parallax.As these motions are not currently included in Gaia's astrometric solution, we develop a template that extracts the component of parallax that would be found when fitting to conventional parallax.The test statistic from this template does not improve as rapidly with time as velocity or acceleration (both of which naturally improve from the extended time baseline, alone) and thus is not expected to be a significant contributor to limits or discovery in the future.However, such an effect could be useful in looking for "outlier" velocities and accelerations as suggested in Paper I, distinguishing them from non-lens sources of motion.
In summary, it is clear that template observables of time-domain lensing are approaching an exciting era, where measured velocities and accelerations will probe regions of parameter space for dark halos that are otherwise unexplored.In addition, there is the possibility to measure this effect in large sets of star-star pairs in the near future.Anomalous parallax may prove a useful tool in future analyses as Gaia continues to provide subtle insights into the nature of the dark matter in the Milky Way.
supported by the Simons Foundation.Some of the results in this paper have been derived using the healpy and HEALPix packages.This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.This work was supported in part through the NYU IT High Performance Computing resources, services, and staff expertise.where β,ecl denotes the component along the ecliptic latitude of the angular impact parameter unit vector.From the equation above, it is clear that the anomalous parallax effect is maximized at the ecliptic equator and suppressed at the poles.In analogy to the proper motion and angular acceleration templates, the parallax test statistic T  and its normalization N  can be defined as When both the proper motion and the parallax uncertainties are statistically limited,   ∝  −3/2 and   ∝  −1/2 , and the above ratio becomes AU/ ⊥  ∼ 10 −3 , for a lens velocity of  ⊥ ∼ 10 −3  and a survey duration of  = 10 yr.When instead the astrometric measurement can resolve the sources' intrinsic proper motion and parallax,   ≃  ⊥ /  and   ≃ AU/  , and the above ratio becomes  ⊥ /  , which is typically of O (1), unless we choose a stellar target with intrinsic velocity dispersion much smaller than the Milky Way's virial velocity.We conclude by noting that the lens parallax motion induces additional corrections to the 2D stellar motion, beyond the one captured by the astrometric fit, given in equation (B3).These corrections -as well as the angular acceleration of equation (5) -are not modeled by a standard 5-parameter astrometric fit (like the one described in Lindegren et al. (2012), and used by the Gaia collaboration).Such an astrometric fit would then perform poorly on lensed luminous sources, leading to large values of the  2 .Unfortunately, there are several competing effects that can degrade the quality of the astrometric solution, and sources with a bad fit are not uncommon in the Gaia data Fabricius et al. (2021); Lindegren et al. (2021).It is therefore hard to rely uniquely on the  2 observable to tease out lensing effects.A more effective strategy would be to perform a dedicated astrometric fit that includes additional parameters to model the lensing-induced acceleration and anomalous parallax corrections.This fit would presumably be computationally costly, but could be applied selectively on sources that are good candidates for a lensing event, i.e. in locations where there is a large proper motion template and/or sources with a large  2 from the 5-parameter fit.

APPENDIX C: ACCELERATIONS FROM THE DR2-DR3 POSITION OFFSET
The astrometric solutions provided by the Gaia collaboration refer to a certain astrometric epoch, chosen so as to minimize the correlations between the position and proper motion fitting parameters Brown et al. (2021).The reference epochs for DR2 and DR3 catalogs are J2015.5 and J2016.0,respectively, and this offset of 0.5 y has to be taken into account when comparing source positions between the two releases.
Maps of the mean positional difference between the DR3 and DR2 Gaia catalogs -having corrected for the difference in epoch -reveal discrepancies that begin at the sub-mas level fluctuating up to 2 mas, and could be attributable to DR2 positional uncertainties Fabricius et al. (2021).Since the astrometric solutions currently presented by the Gaia collaboration do not capture the full range of dynamics that may govern apparent stellar motions, these mismatches could also encapsulate accelerations, including of lens-induced, apparent accelerations in stars that have been fit only for parallax, position, and proper motion.
Reconstructing angular accelerations for stars by combining their astrometric measurements at different epochs was made possible through the cross-calibration of the Hipparcos and Gaia catalogs in Brandt (2018Brandt ( , 2021)).The clear advantage in combining observations from these two missions is the long baseline of ∼ 24 y, while the drawback is the limited number of available sources, just over 115, 000.
Lacking access to Gaia's transit timing information for each source, which would allow us to fit for nonlinear trajectories and look for lens-induced accelerations, we present here the derivation of an alternative estimator for calculating stellar accelerations {  }, by leveraging the difference in epoch between different data releases.By cross-matching the DR2 and DR3 Gaia datasets, we can leverage the larger number of matched sources and create a catalog that can serve as a look-up table for stellar accelerations, albeit potentially being limited by the shorter observational baseline of ∼ 0.5 y between the two releases.A machine-learning algorithm that takes advantage of this difference in Gaia's observational timelines has generated a catalog of nearby accelerating star candidates at high statistical significance Whiting et al. (2023).We propose a complementary analytic approach that can be used on all sources with a 5-parameter astrometric solution.
We simplify Gaia's astrometric fit by only considering the stars' 2D motion in the plane perpendicular to the line of sight.We further assume uncorrelated uncertainties    per individual observation in both angular directions, and that the position of each source is recorded at regular time intervals with frequency  .Furthermore, each source is observed for a total time  and Gaia's best-fit parameters { θG , μG , πG } are given at the midpoint of the observations.The astrometric fit will be the result of minimizing the following test statistic

Figure 5 .
Figure 5. Angular number density of foreground stars in the sample of optical doubles selected from Gaia's DR3 on a pixel size of about 0.45 deg in Galactic coordinates.The visible large-scale structure is caused by extinction dropping the background source density in the optical wavelengths.(left) Sample of 15.1 million pairs after the cuts on RUWE and stellar distances described in section 3.2.1.(right) Sample of 11.7 million pairs after the background subtraction and proper motion outlier removal described in section 3.2.2; the visible spatial features are due to the spatial HEALPix pixelation on a scale of about 3.7 deg during the background subtraction and subsequent removal of stars falling in "sparse" bins (see the main text for further details).ƭ

Figure 6 .
Figure 6.(Top row) Number density of background stars around their lens counterpart in equatorial (left), ecliptic (center), and galactic (right) coordinates.For all stellar pairs, the lens is at the origin, and the black arrow in the left panel points in the direction of the average lens velocity over all pairs.The number of background stars decreases significantly below ∼ 1.5 ′′ and it is not uniformly distributed.(Bottom left) Number density of background stars around their lens counterpart as a function of the parallel and perpendicular components of the angular impact parameter with respect to the lens velocity direction where, for all stellar pairs, the lens is at the origin; in this plane, all the lens velocities are directed along the positive horizontal axis and the anisotropy is not visible, suggesting that it is not related to the lens velocity direction.(Bottom right) Distribution of the proper motions of the star-lenses.ƭ

Figure 7 .
Figure7.Proper motion dispersion of the background stars in the final sample of optical doubles stars in Gaia's EDR3, after quality and distance cuts, background motion subtraction and outlier removal (see sections 3.2.1 and 3.2.2). (left) Histograms of the number of stars per G magnitude bin (gray), average proper motion error reported by Gaia (blue), and average effective dispersion of the stars' proper motion in each bin (red).(right) Same as the left panel, but as a function of the background-foreground stellar angular separation.For both the effective and the Gaia error, the total proper motion variance is computed as  2  ≡  2  +  2   + 2         , where     is the correlation between the proper motion in the  and  directions.ƭ

Figure 8 .
Figure 8. (left) Expected SNR, N ★  , for the star-star proper motion template on the clean sample of 11.7 million pairs as a function of lens and source distance.(center) Same as the left panel but for the observed T ★  /N ★  .The region with   < 300 pc contains most of the sensitivity and shows the smallest fluctuations in T ★  /N ★  ; for   > 300 pc we observe high positive fluctuations around   ≃ 7 kpc which are most likely due to a systematic error that has not been removed by the cleaning procedure.(right) Expected SNR, N ★  (thick blue, left axis), and observed SNR, T ★  /N ★  (red, right axis), as function of the lens-source angular separation  ls .The thin blue line denotes the expected scaling of the SNR as √︁  pairs /(   ,eff ) and is arbitrarily normalized to approximately match the bins with the largest SNR.Here  pairs and  ,eff denote the total number of pairs and the effective dispersion averaged over all the background stars in each bin.ƭ ) −   ( θG , μG , πG |  )] 2 Limits and projected sensitivity on the fraction of dark compact lenses   with mass   and characteristic size   = 1 pc from the proper motion template analysis on MC stars.

Table
is computed as  2  ≡  2  +  2   + 2         , where     is the correlation between the proper motion in the  and  directions.ƭ

Table 2 .
Same as Tab. 1 but for a sample of 814 thousand cleaned Gaia DR3 optical doubles with lenses closer than 300 pc to the observer.ƭ The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation.Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade.canbefound by minimizing |Δ  () − Δ  π | 2 and averaging over one period.The resulting lens-induced anomalous parallax isΔ  = − 4           (  ,   ),(B3)where here   ≃  0  denotes to the zeroth-order impact parameter and the spatial profile is   (  , ) =

T
(  ,   ) = ∑︁      (  ,   )  (  ,   ) =  is the measured parallax of the i-th star and  , the parallax dispersion.Following the signal-to-noise ratio derivation of section 2.2, we can forecast the sensitivity of the parallax template usingSNR T  ≃ 4   AU      = AU/  and assumed the maximum signal by taking sin  ecl = 0, for simplicity.We can compare this with the SNR for the proper motion template