A wealth of information about our universe and galaxy is contained in the spectrum of its density fluctuations and the gravitational influence they exert on other objects. All evidence for dark matter (DM) is, so far, of this kind: gravitational back-reaction on the cosmic microwave background, large-scale structure formation, cluster-and galaxy-scale velocities, and weak gravitational lensing on extra-galactic scales. From these and other indirect gravitational probes, we have learned about our cosmological history and the properties of DM and astrophysical systems on large scales.
There also exists a "dark world" on small scales. Most types of compact objects, such as astrophysical black holes (BHs), neutron stars, white dwarfs, brown dwarfs, and planets generically emit or reflect too little electromagnetic radiation to be detected directly, except if they are young, close, and/or accreting. This dark world may also be populated by small DM structures, such as (ultra-compact) minihalos [1][2][3][4][5] or more exotic objects such as primordial black holes (PBHs) [6], boson stars [7][8][9][10][11], and other composite DM objects [10,[12][13][14][15]. These clumps and structures may be invisible to us, but their presence can occasionally be revealed indirectly through gravitational waves [16,17], direct gravitational effects on visible stars [18][19][20][21], pulsar timing arrays [22][23][24][25][26][27][28][29][30], and gravitational lensing of light [31][32][33][34][35] and of gravitational waves [36][37][38][39].
Time-domain, astrometric, weak gravitational lensing of light has emerged as one of the most promising probes of compact objects in the Milky Way (MW) [40][41][42][43][44][45][46][47][48][49][50][51]. Following the foundational works proposing the astrometric and photometric observables of transient gravitational lensing [52][53][54][55], most observational efforts have relied primarily on photometric signatures in the strong lensing regime [56][57][58], e.g. managing to exclude PBHs comprising the totality of the DM abundance over a wide mass range [59][60][61][62][63][64][65][66]. The power of astrometric signatures has recently received an enormous boost from simultaneous advances in catalog size, observational cadence frequency, and positional precision of astrometric surveys, most notably that of the Gaia satellite [67,68], with great prospects for astrometric microlensing [69] and already many photometrically detected events [70], in addition to interesting candidate events from other surveys [71]. The (weak) astrometric gravitational lensing deflection signature decouples more slowly with increasing impact parameter than the (strong) photometric magnification, so it has parametric advantages in searches for rare dark objects [41]. Ref. [40] proposed a host of observables for time-domain astrometric weak lensing by dark objects: matched filters [72], and correlation functions or power spectra [73] of lensing-induced, correlated proper motion and acceleration corrections for many stars, or transient astrometric deflections of single (or multiple) stars.
In this paper, we present a robust data analysis pipeline to extract significant events of transient astrometric lensing on single stars, along with associated software tools and a procedure to generate faithful mock catalogs of compact objects in the MW. Our pipeline is developed with Gaia's fourth data release (DR4) in mind, but is applicable with minor modifications to other astrometric data sets (e.g. HSTPROMO [74] and PHAT [75]). Our robust and near-optimal data analysis pipeline is projected to detect several isolated astrophysical BHs in the MW (and perhaps other compact remnants such as neutron stars and white dwarfs), while having leading sensitivity to compact DM objects with masses between 1 M and 10 3 M .
In related work, ref. [76] expanded on the sensitivity estimates of ref. [40] by projecting the sensitivity of Gaia time series data to PBHs using a probabilistic model, and ref. [77] JCAP07(2023)037 reported the possible presence of a dark point-like lens in the observation of a single Gaia Data Release 3 (DR3) source based on a poor astrometric model fit. In this work, we faithfully produce mock data sets mimicking Gaia DR4 time series data which include not just statistical noise, but also backgrounds from astrophysical BHs and from binary systems. Additionally, we create an analysis pipeline that can be applied to Gaia DR4.
In section 2, we review the basics of astrometric observations and data products in Gaia, and how they can be affected by lensing dynamics. Section 3 details the generation of our realistic mock catalogs, while section 4 describes the steps in our data analysis. The results of data analyses on our mock catalogs are presented in section 5, and we conclude in section 6. Supporting materials such as derivations, extra plots, and minor results can be found in appendices A-F. The data and code are available on GitHub ( ), with links (ƭ) below each figure.
We primarily use two models of astrometric motion in our proposed search for dark compact objects in the Gaia DR4 data. We call the first the free model. It describes the apparent motion of a source moving across the sky without being subject to any gravitational effects, neither local nor along the line of sight (astrometric gravitational lensing). For trajectories across small patches of the sky, this motion can be modeled as entirely inertial. The second type of model, which we call the blip model, describes the apparent motion of a source subject to lensing due to a massive compact foreground object. By comparing the goodness-of-fit of these two models to any given source trajectory in the Gaia catalog, we may quantitatively probe various compact DM scenarios, as well as discover singular dark compact objects in the real Gaia data.
We analytically model the apparent astrometric motion of an unlensed or "free" source across the sky, as well as the motion of point-like lenses, as a function of five parameters. The model we employ is the angular component of the "standard model" of stellar motion described in refs. [78,79]. The angular barycentric coordinates of a free point-like celestial body θ θ θ free = (α * free , δ free ) in the standard barycentric celestial reference system (BCRS, [80]) at any given time t (with respect to some fixed reference time t 0 ) are given by θ θ θ free (t | θ θ θ 0 , µ µ µ, D) = θ θ θ 0 + µ µ µ(t -t 0 ) + (t | θ θ θ 0 , D), (2.1) where θ θ θ 0 is the BCRS parallax subtracted position of the body at reference time t 0 , µ µ µ = (µ α * , µ δ ) is the constant angular velocity of the body in the sky, and (t) is the parallax correction to the linear trajectory given by
sin(α) -cos(α) 0 cos(α) sin(δ) sin(α) sin(δ) -cos(δ)
x x x E,cart (t).
(2.2)
Here, D is the line of sight distance to the object, and x x x E,cart (t) are the Cartesian coordinates of Earth in the heliocentric frame, which we assume to follow a purely elliptical trajectory. We note that eq. (2.1) is equivalent to the 5-parameter astrometric model Gaia use to model each source trajectory that they measure. Hence the set of parameters (θ θ θ 0 , µ µ µ, D) are the same as reported by Gaia in all data releases thus far, except Gaia uses parallax as a parameter JCAP07(2023)037
. Astrometric lensing geometry. A point-like lens l at a line-of-sight distance D l from an observer o creates two displaced images s + and s -of a background source s at an angular impact parameter β and line-of-sight distance D s . The displaced images are separated from the true source location s by angles ∆θ + and ∆θ -as specified by eq. (2.6). We average the location of s + and s - weighted by their relative magnification to obtain a single lensed source location.
instead of distance D. The two are equivalent since D = (1 arcsec/ ) pc. 1 An example of a free trajectory can be seen in figure 2.
We additionally model source trajectories undergoing constant angular acceleration. We achieve this by adding two extra parameters to the free model
where γ γ γ = (γ α * , γ δ ) is the constant angular acceleration of the celestial body in the sky.
We use the acceleration model to discriminate between long period binaries and blips. See section 5.3 for more details.
We model the trajectory of a celestial body subject to detectable transient astrometric lensing caused by a point-like lens as a function of 11 parameters. We call this the "blip" model of celestial motion [40]. In addition to the 5 free motion parameters of eq. (2.1), there are 6 additional parameters: the position of the lens θ θ θ l,0 at reference time t 0 , the proper motion of the lens µ µ µ l , the distance to the lens D l , and the mass of the lens m l . In its most basic form, the blip model may be written as
We calculate the lensing deflection term ∆θ θ θ(t) assuming a point-like lens, and we employ the thin-lens approximation, in which we assume the lensing deflection takes place over a region that is very small compared to the line-of-sight distances involved. The point-like lens assumption allows us to construct a model that is valid in both the weak and strong lensing regimes. These approximations are valid as long as the Newtonian potential of the lens is small and the relative velocities of the observer, lens, and source are small compared to the speed of light, which is the case for all sources in the Gaia catalog. The Einstein radius θ E of a massive point-like object is given by
where D s and D l are the distance to the source and lens, respectively, m l is the mass of the lens, G is the gravitational constant, and c is the speed of light [81]. Using the Einstein radius, we then calculate the deflection of the two images created by the lens as
with relative (signed) magnification
where β β β(t) is the angular impact parameter pointing from the lens to the source, and the dimensionless impact parameter u(t) ≡ |β β β(t)/θ E |. The absolute value on the left hand side accounts for the fact that the inversion of the second image can be ignored since it is point-like. Given Gaia's point spread function (PSF) width of about 2 pixels or 100 mas [67], the two lensed source images are rarely resolved individually (especially if the lensing occurs inside the MW), meaning Gaia will usually only resolve the light centroid of the two images. Via eqs. (2.6) and (2.7), the light centroid deflection due to lensing is given by
which we insert into eq. (2.4) to obtain a complete expression for lensed trajectories in the Gaia catalog. We provide a schematic of the lensing geometry and notation in figure 1, and we show a realistic blip trajectory in figure 2. By contrast, adding the two image magnifications in (2.7), we obtain a total magnification of
which leads to an effective change in source magnitude of ∆Mag = -2.5 log 10 (1/µ). In figure 3, we show the astrometric deflection given by (2.8) and the photometric magnification given by (2.9) as a function of impact parameter, as well as the asymptotic behavior of each observable, for a gravitational lens with Einstein radius θ E = 10.0 mas.
In this section, we describe our method for creating mock catalogs that closely resemble the data products from the upcoming Gaia DR4. First, we discuss how to extrapolate the 5-parameter astrometric solution reported by Gaia EDR3 into the time-series data expected in Gaia DR4. Then, we describe the models we adopted for generating astrophysical BHs and compact DM. The mock catalog provides a way to understand the statistical background for event selection, detectable lensing events, and the projected compact DM constraints which are shown in section 5.
We take all the sources in Gaia EDR3 that have a 5-parameter astrometric solution and generate time-series data in the proposed format of Gaia DR4. By using astrometric parameters and stellar magnitudes directly from EDR3, we automatically capture extinction, crowding, and instrumental effects, which normally must be treated carefully in catalogs based on the injection of fully artificial stars. Some of the stars in EDR3 have negative parallaxes and large parallax uncertainty. To circumvent this issue, we take the median of the inferred distance posterior of each star with geometric and photometric priors prescribed in ref. [82]. With the unlensed mock catalog, we can test the false positive rate for lensing events and determine the distribution of our test statistics under the null hypothesis. We also inject lenses using astrophysically realistic priors on their phase space distribution to construct a lensed catalog. The epoch astrometry due to be released in DR4 will not provide timestamped twodimensional BCRS coordinates due to the scanning law of Gaia [67]. Instead, each epoch measurement will be reported as a one-dimensional displacement θ(t) with respect to a scan angle φ(t) in the so-called "Along Scan Direction" (AL) in the Gaia documentation. We convert the coordinates given by our model to this data format using the relation
where (α 0 , δ 0 ) are the BCRS coordinates of the source at a reference time t 0 provided by Gaia. Only the brightest stars will have a location offset in the perpendicular "Across Scan Direction" (AC). For simplicity, we will only use the AL location for all the stars in our mock catalog. The timestamp and scan angle for each epoch will be the same for all stars in the catalog.
Figure 3. Photometric and astrometric observables for a single lensing event, caused by a lens with a 10.0 mas Einstein radius passing by a luminous star. In solid red, magnification observable µ -1 as a function of impact parameter, with the asymptotic behavior (µ -1 ∼ (θ E /β) 4 for β θ E ) shown in dashed red. In solid blue, astrometric deflection observable ∆θ as a function of impact parameter, with the asymptotic behavior (∆θ ∼ θ E /β for β θ E ) shown in dashed blue. The lens Einstein radius is shown in green. ƭ two hours apart from the first set of 40 points (with the same set of scan angles) are added to the time series to mimic the scanning law described in ref. [67], for a combined total of 80 data points.
We note that our pipeline is also capable of using Gaia's Observation Forecast Tool (GOST) to obtain more accurate scan angles and observation timestamps for each source. 2In appendix C, we further discuss GOST and show limits on dark compact objects obtained using a mock catalog generated with GOST, analogous to the limits shown in section 5. We also discuss how GOST affects Gaia's ability to discover BHs. We emphasize that these limits and Gaia's discovery potential are only marginally different to the ones obtained in the simplified data scenario where all sources are observed exactly 80 times. Therefore, all subsequent sections assume this simplified scenario.
We inject isolated, electromagnetically quiet BHs and compact DM objects into the mock catalog. The priors for generating these two different populations are specified in the following.
MW stellar evolution simulations suggest that there should be of order 10 8 BHs in the MW [83], yet we have only observed a handful through the emission of electromagnetic waves from accretion and photometric microlensing. Gaia DR4 will provide an opportunity to discover isolated, non-accreting BHs via transient astrometric lensing.
BH Surface Number Density Left: BH surface number density obtained by integrating the BH volume number density inside a sphere of radius 5 kpc centered on the solar system, taken to be 8 kpc away from the Galactic Center. We assume there are a total of 10 8 BHs present in the entire Galactic thin disk [83], and 10 7 within 5 kpc of the Sun. Right: DM surface mass density obtained by integrating the BH volume mass density inside a sphere of radius 5 kpc centered on the solar system. We assume the MW DM density follows an NFW profile of scale radius 18 kpc with a value 10 -2 M /pcfoot_1 at the Sun's location. ƭ Since astrophysical BHs are remnants of stellar evolution, we assume that their distribution in the sky closely resembles the MW stellar distribution. The stellar population in the MW is commonly decomposed into the Galactic bulge, thin disk, thick disk, and the Galactic halo. The thin disk is of primary relevance for our purposes, due to its high stellar number density and its proximity to Earth. We model the Galactic thin disk with the exponential function
where n 0 is the central stellar number density, and z d and R d are the scale height and scale radius of the thin disk, respectively. Ref. [84] reports z d = 300 pc and R d = 2.6 kpc. Simply using stellar distributions to model the MW BH distribution does not account for BH natal kicks. These kicks -caused by the dynamics of supernova explosions -offset the BH velocity distribution from that of MW stars. These kicks explain the observed distribution of low mass X-ray binaries far away from the Galactic disk [85,86], because the BH velocity gain due to kicks will increase the scale height z d of the BH distribution relative to that of the stellar distribution, effectively "puffing up" the disk. We estimate in appendix D that the scale height will increase by a factor of about 10 due to this effect, so for astrophysical BHs, we use z d = 3 kpc and R d = 2.6 kpc. The surface number density of BHs across the sky is shown in figure 4. The probability density function (PDF) for BH distances D l at a given celestial location in galactic coordinates (l, b) is then
The combined PDF of BH proper motion and distance is
which we normalize such that P BH (µ µ µ l , D l |l, b)d 2 µ l dD l = 1. For a detailed derivation of the conditional PDF P BH (µ µ µ l |D l , l, b), see appendix D. We adopt the BH mass distribution reported by LIGO-Virgo [87] obtained from a combination of 47 binary BH merger observations. We thus assume -for now -that the BH mass distribution is similar for single BHs and for binary BHs. 3 We also assume that JCAP07(2023)037 the BH mass is independent of the position and the proper motion of the BH so that the two PDFs P BH (µ µ µ l , D l |l, b), P BH (M BH ) are separable. The model we use is the Power Law + Peak model reported by LIGO-Virgo, wherein the BH mass distribution follows a power law with a soft cutoff at the lower end and a hard cutoff at the upper end. A peak is added, motivated by a potential pile up of BHs just before the pair-instability gap of supernovae [88]. The resulting BH mass function is shown in figure 5.
Compact DM objects may comprise part or all of the DM abundance and thus produce transient astrometric lensing signals in Gaia DR4. A non-detection would set constraints on the fraction of DM composed of such compact objects (e.g. PBHs) as a function of their mass. Here, we only consider point-like sources, specifically lens objects with scale radii smaller than their Einstein radius
(In appendix A, we discuss the limitations on detecting lensing events from lenses with extended density profiles.) We assume that the DM distribution in the MW follows a Navarro-Frenk-White (NFW) profile [89] with a fiducial scale radius R s = 18 kpc and a local DM
and that the DM has a Gaussian velocity distribution The Gaia EDR3 error function [90], showing the median per transit astrometric error of a given source as a function of G magnitude. Right: The Gaia EDR3 G magnitude distribution [90]. Since we construct our mock catalogs based on EDR3, the astrometric errors are distributed exactly according to these two distributions. ƭ where σ DM = 166 km/s. The surface mass density of DM across the sky is shown in figure 4. With these parameters, there is roughly 7.0 × 10 10 M of DM mass within a 13 kpc radius around the solar system, corresponding to the 99th percentile of the stellar distances in our mock Gaia DR4 catalog (based on EDR3).
We perturb each astrometric positional data point generated via the free and blip models by subjecting the mock source trajectories to Gaussian noise. Since we base our mock catalogs on Gaia EDR3, we draw directly from the EDR3 error distribution. In practice, this is done by using the error function described in ref. [90] to convert each EDR3 source's reported photometric mean G magnitude into a Gaussian standard deviation quantifying the instrumental astrometric precision in the AL scan direction for a single transit. We then randomly shuffle each positional data point in every source trajectory by drawing from a normal distribution centered at each true source position and with standard deviation corresponding to the per transit error. The EDR3 error function and the EDR3 G magnitude distribution are shown in figure 6. Note that using the EDR3 error function is conservative, since errors are projected to decrease in future data releases across all G magnitudes [91]. For simplicity, the error function we use here is only a function of the stellar magnitude. In reality, the error function is position dependent as shown in ref. [90]. In crowded regions, such as inside the Galactic Bulge, uncertainties may be larger than for stars of similar magnitude located outside the bulge. However, since we directly apply the median Gaia EDR3 error function to our analysis, these error anisotropies should not significantly affect the final results reported in section 5.
In this section, we describe our construction of a data analysis pipeline to detect true blip events and set constraints on dark compact object populations in both the true Gaia DR4 catalog and the mock catalogs described in section 3. The pipeline systematically goes through an entire catalog and optimizes a set of test statistics for each source in order to discern the probability that any given source trajectory is a true blip event. By making cuts in the significance level of different test statistics, we can thus discriminate between blip and free stellar trajectories, and thus discover and flag true blip events effectively. We can also obtain limits on the compact object DM fraction in the MW using the Yellin method [92,93] applied on the distribution of these test statistics.
As pointed out in section 3.3, we assume the astrometric Gaia DR4 data to be subject to pure Gaussian noise, with the positional error of each source corresponding to its G magnitude. Hence, we use a Gaussian likelihood function to quantify the agreement between the astrometric data and our choice of model (either free or blip). Given a dataset θ θ θ obs = {θ n,obs } where the subscript n labels each data point in the source trajectory, as well as either 5 parameters y y y = y y y free (free model) or 11 parameters y y y = y y y blip (blip model), we may write the corresponding likelihood function as
where θ θ θ model (y y y) = {θ n,model } is the prediction given the model parameters y y y, and σ n is the error associated with the data point n. We then define our blip test statistic (TS) to be TS(θ θ θ obs ) ≡ -2 max namely, the test statistic for any given source trajectory is defined as the maximized log likelihood ratio between the free and blip model fits to the source trajectory data. We note that the negative log likelihood ratio is equivalent to the difference in χ 2 goodness of fit values between the two models. It should also be noted that under the assumption of trivial covariance between model parameters, the distribution of maximized test statistics follows a true χ 2 distribution in the asymptotic limit [94]. While eq. ( 4.2) provides a way to evaluate the quality of fit of our model to the data, the expression does not contain any prior information on the lens population being probed. To constrain our search, we therefore construct a second test statistic based on the posterior of a lensing event, rather than the likelihood. We define this constrained test statistic (TS * ) as
where max * indicates that rather than maximizing the blip likelihood directly, we are instead maximizing the log of the posterior probability associated with each source trajectory P post = log L(θ θ θ obs |y y y blip )P (y y y lens ) , (
where P (y y y lens ) is the prior probability density of the lens parameters, with the exact form of the prior depending on the lens population being probed, as described in section 3.2. Note
that the quantity inside the square brackets has nontrivial units, but these can be neglected since they amount to a constant offset in the test statistic and hence do not matter if eq. ( 4.4) is used as a loss function only. A further constraint implied by max * is the requirement blippiness(y y y blip ) ≡ t obs µ rel
where µ rel is the relative (linear) proper motion magnitude between the source and the lens, t obs is the total observation time, and β min is the minimum angular impact parameter between the lens and the source. We have coined the above quantity the blippiness of an event, as it is simply the ratio between the relative angular distance traversed by the source and lens over the full observation time τ , and the minimum angular impact parameter.
There are two reasons for imposing these extra constraints when maximizing the log likelihood ratio. First, maximizing the posterior rather than the likelihood means that we penalize choices of model parameters that are unphysical. Similarly, were we not to impose the blippiness constraint, we would be probing parts of parameter space which cannot produce a significant blip, simply because events that have a large minimal impact parameter are either too long or the lensing deflection is too weak to produce a signal. Second, imposing these constraints guides our choice of minimizer to a physical part of the blip parameter space, which reduces the amount of computational power needed to compute test statistics for all 2 × 10 9 events in the Gaia catalog.
We emphasize that constraining the maximization in eq. ( 4.3) only reduces the value of the test statistic compared to what would be obtained by calculating eq. (4.2), meaning the full test statistic distribution gets shifted to smaller (or even negative) values. However, for true blip events, the reduction in significance is minimal due to the distribution of true blip parameters coinciding with the prior probability distribution in eq. (4.4).
Finally, we note that our analysis pipeline is also capable of incorporating the prior on stellar distances reported in ref. [82]. This prior helps overcome Gaia's difficulty in determining stellar parallaxes for faint sources or sources located in the galactic bulge, where effects from blending and crowding can be significant. We do not include this prior in the analysis reported here; however, we tested how it affects the results reported in section 5 and found no significant difference. However, for the real DR4 data set, where bad parallax measurements have a more significant impact on the analysis, the Bailer-Jones prior will be beneficial. It is therefore enabled by default in the analysis software.
We employ the optimum interval method developed by Yellin [92,93] to determine (projected) limits on the DM fraction f l in compact DM objects. The Yellin method is suited to hypothesis testing of a known signal model in the presence of an unknown background distribution, in a fixed region of interest. For a one-dimensional distribution of events, it entails computing the integral of the signal distribution of all intervals of n events and assesses whether the largest interval significantly exceeds the expectation for the signal model, in which case the signal hypothesis is rejected.
In our analysis, the events are the constrained test statistics for all of the stars. We can compute the distribution of test statistics under the signal (lensing) hypothesis numerically by drawing compact DM objects from the distributions specified in section 3.2.2. For computational efficiency, we only consider stars in the distribution whenever a lens is present within a threshold impact parameter which causes a maximum deflection of at least ∆θ = 5 µas.
The background distribution is obtained by fitting the unlensed catalog; the background events are the large upwards statistical fluctuations in the constrained test statistics. Furthermore, we can consider a mock catalog contaminated with lensing by astrophysical BHs, and by binary systems with an undetected companion as astrophysical backgrounds.
The recipe of implementing the optimal interval method in this work is the following:
1. Generate a test statistic distribution only for stars that have a nearby lens. We call it the "signal distribution" S(TS).
2. Given the test statistics of the experiment, compute the maximum of expected number of events between all pairs of events e i , e i+(n+1) , which is the integral of S(TS) between e i , e i+(n+1) . We call this the maximum interval x n .
3. Generate many instances of Monte Carlo realizations of the signal events and perform step 2 on all of the realizations. The DM fraction f l is simply a scaling factor in the signal distribution S(TS). Following the steps outlined above, we find the limiting f l, * such that the DM fraction f l ≥ f l. * is excluded at 90% confidence level. For a more detailed discussion on the Yellin method, see refs. [92,93].
Gaia DR4 will contain time series data for about 2 billion sources. Scouring this vast catalog for blip events is a considerable computational challenge and requires a structured approach. We construct a modular analysis pipeline wherein key statistical assumptions, such as the lens priors, can be swapped to search for blips from different lens populations. Ancillary data from e.g. photometric surveys can also be incorporated via these priors.
A diagrammatic representation of the analysis pipeline's flow is shown in figure 7. We first fit the free model from section 2.1 to every source in the catalog. To do this, we employ SciPy's minimize function [95] to maximize the logarithm of eq. (4.1). This yields an optimized log likelihood -2 log Lfree value for each source trajectory, where the "hat" indicates that the likelihood has been maximized with respect to the source trajectory. After performing the initial fit, we impose our first cut. Any significant blip event should have a small optimal likelihood under the free trajectory hypothesis; we discard any events with an optimized negative log likelihood of -2 log Lfree < χ 2 5σ , where χ 2 5σ is the 5σ significance threshold of the -2 log Lfree distribution computed via Monte Carlo (MC). This distribution asymptotically matches a χ 2 distribution with m -5 degrees of freedom, where m is the number of data points in a given observation (see section 5.1), where χ 2 5σ = 152 for 75 degrees of freedom (all trajectories in the mock catalog consist of 80 data points) computed by matching the χ 2 distribution to the Gaussian 5σ p-value of 5.7 × 10 -7 . For any events that pass this cut, we rerun the free model fit, but this time using a nested sampling procedure using the Bayesian Gaia Catalog Figure 7. A flowchart representation of the analysis pipeline. The pipeline reads in astrometric data from the input catalog. Then, the free model is fitted to every source in the catalog, generating a set of optimized log likelihoods {-2 log Lfree }. Any events that are below a 5σ threshold in the free model χ 2 distribution is excluded. To ensure that all free fits have converged to global minima, we then rerun the same optimization procedure, except we use an nested sampling based optimizer, yielding a new set of log likelihoods {-2 log L free }. We then reimpose the 5σ cutoff on the new computed free log likelihoods. Any event that passes these cuts is then tested against the acceleration model, also using the nested sampling optimizer. Finally, we test the remaining events against the blip model by computing TS * . All events that pass the initial cut, as well as a 3σ cut in acceleration, and which satisfy TS * > 100 are flagged as blip events. The TS * distribution for events passing the initial 5σ free cut is also returned, which the pipeline uses to impose constraints on lens populations using the Yellin method.
inference tool PyMultinest [96][97][98][99]. This ensures that the global maximum of each free model log likelihood is found. Should any of the remaining sources fall under the 5σ threshold after this second fit, they also get discarded. This cut yields the most significant reduction in computational resources needed to search for blips, since it reduces the number of sources of interest by 6 orders of magnitude.
We then fit sources that pass the first two cuts against the acceleration model. This extra fit is primarily implemented to account for binaries (see section 5.3 for more details). Like with the free fit, we minimize -2 log L accel using first SciPy and then PyMultinest, yielding a set of optimized log likelihoods -2 log Laccel .
Finally, we compute the constrained test statistic TS * for each remaining source using again first SciPy and then PyMultinest to ensure convergence to global maxima. Any event that passes the 5σ free model cut, is above 3σ significance under the assumption of the acceleration model, and has a test statistic TS * > 100, is flagged as a blip event. Furthermore, for these events (and any other event that passed the initial 5σ free fit cut), the pipeline JCAP07(2023)037 The gray areas indicate regions in which an event meets our cut criteria and gets excluded from the list of blips found by the analysis pipeline. Note that all three of these true blip events cross the 5σ free log likelihood threshold and the TS * > 100 requirement at roughly the same value in the varied blip parameter. Also note that at very low blippiness or β min , the blip events enters the strongly lensed regime, which is the cause of the flatness of both the test statistic and likelihood in this range. ƭ outputs a list of test statistics (-2 log Lfree , -2 log Laccel , and TS * ), each model's best fit parameters and corresponding uncertainties, and nested sampling generated parameter space covariance data. See figure 8 for an example of the pipeline's sensitivity to changes in various lens parameters. The pipeline is finally also able to run a Yellin test on the computed TS * distribution and can generate 90% confidence limits on compact DM parameter space.
We run the data analysis pipeline of section 4.3 on the mock catalogs described in section 3 to test its ability to discover true blip events in quasi-realistic data, and to make projections for the discovery potential and expected constraints in Gaia DR4. We first apply the pipeline on a mock catalog unperturbed by lensing to quantify the distribution of test statistics generated by the analysis procedure, as well as to ensure that the pipeline is robust against random noise, misfitting errors, and other artifacts. We then run it blindly on the astrophysical BH catalog described in section 5.2 in order to probe its ability to detect this astrophysical signal In purple, the (upper) 5σ threshold of the analytic χ 2 /DOF distribution. The χ 2 distribution obtained by the analysis pipeline matches the theoretical distribution nearly perfectly, which is expected since this particular catalog only contains sources with free trajectories subject to Gaussian noise. The strong agreement demonstrates the robustness of the pipeline. ƭ that is guaranteed to be present in the data. Next, we run the pipeline on a series of mock binary events with dark companions to ensure that the pipeline will not flag binaries with dark companions as blips. Finally, we use the pipeline on the compact DM catalog described in section 3.2.2 to generate mock Yellin 90% limits on the compact DM fraction in the MW.
We first analyze the mock catalog consisting of 1,447,353,154 Gaia sources propagating freely across the sky, without undergoing any sort of lensing deflection. Figure 9 shows the -2 log Lfree distribution for these events. The log likelihood distribution closely follows an analytic χ 2 distribution, in line with expectations for Gaussian noise injection only, and highlighting that the free model's 5 parameters have minimal covariance. Some events in this catalog pass the initial 5σ cut in the free log likelihood distribution. This is expected from statistical noise and the sheer number of events in the catalog. Upon computing TS * , however, we see that none of the events in this catalog pass the TS * > 100 requirement for an event to be flagged as a blip. In fact, all of the events satisfy TS * < 60, meaning none of the events are even remotely close to being considered as a highly significant blip event. The stringent cuts in log likelihoods and in TS * effectively preclude statistical fluctuations from being classified as blips, at least under our assumption of high-quality data with Gaussian noise.
We analyze the mock catalog described in section 5.2 to test the pipeline's ability to search for isolated astrophysical BHs in the MW. We conduct this blip search blindly. A total of 6 events JCAP07(2023)037 The residuals from the fits, with the 1σ and 2σ bands being shown in green and yellow, respectively. ƭ pass the 5σ free model cut and our TS * > 100 requirement. Out of those six, two do not pass the 3σ cut after the acceleration fit. Upon comparing with truth information (unblinding), we learn that all six of these events are true blips, demonstrating that the pipeline is capable of flagging astrophysical BHs in the Gaia data and simultaneously not generating any false positives. These 6 events, their statistics, and their best fit parameters are shown in table 1. Furthermore, the raw AL scan fits and residuals for two of these events are shown in figures 10 and 11; the remaining four plots are available on GitHub ƭ. Finally, figure 12 shows the covariance between blip model parameters at the global maximum constrained log likelihood ratio (i.e. where TS * is computed) for one of the six events; the other five corner plots are available at this link ƭ.
We also compute the source brightness magnification due to photometric lensing for each of these 6 events. Of the 4 that pass the 3σ acceleration fit cut, only one event has a magnification above the Gaia photometric uncertainty. This demonstrates the advantage of looking for lensing signals with astrometric surveys. See appendix E for a more detailed discussion and supplemental photometry plots.
One of these events (top row, second column in table 1) has best fit values particularly close to the true lens parameters with narrow error bars. This is because the lens is rather close and has a high blippiness value. This event breaks much of the parameter degeneracy that plagues more distant and less significant blip trajectories. The lens distance degeneracy with mass and proper motion can also be seen in figure 12, and is much stronger for the other 5 blip events. These degeneracies explain why the parameters that maximize the constrained likelihoods do not necessarily coincide with the true lens parameters.
We conclude that we likely expect to see about four true blip events after both the acceleration and TS * cut, the closest and most blippy of which will have accurately determined lens parameters. The sources in question and candidate lens locations should then be followed up by other telescopes, providing exciting prospects for the study of phenomena associated with free-floating astrophysical BHs: e.g. accretion from the interstellar medium [100,101], and superradiance [102][103][104][105][106][107]. A free-floating astrophysical BH has only been claimed to have been detected once in the past [108][109][110] The residuals from the fits, with the 1σ and 2σ bands being shown in green and yellow, respectively. Note that despite the large uncertainties in the parameters of this source and the other 4, the blip model still provides an excellent fit. This is due to parameter degeneracy. ƭ
The exact fraction of stars in binary or higher order systems in the MW has not been accurately estimated, but surveys of Sun-like stars in the solar neighborhood suggest that it may be approximately half of all stars [111]. Binaries that are entirely or partially resolved have been studied extensively using Gaia data [112], and many of these sources are automatically flagged in Gaia's public data releases. Therefore, we may simply discard them from our analysis of the full astrometric DR4 catalog. Gaia's DR3 binary flagging procedure is described in ref. [113]. Gaia does not flag sources in binary orbits with a dark (or faint) companion, such as a neutron star, a brown dwarf, an exoplanet, or an astrophysical BH. It is known that the binary orbits of these sources induce a measurable correction to the free trajectory of the source [19,20,[114][115][116]. In particular, it is estimated that Gaia is capable of observing about 75 sources with BH companions [20]. Gaia's binary flagging system is also conservative since flawed binary flagging can hurt Gaia's science output -it is projected that some fully luminous binaries will bypass Gaia's flagging procedure. Ref. [117] describes Gaia BH1: a binary system consisting of a G-type star orbiting a BH. This event is well described by the Gaia DR3 binary orbit astrometric solution. It is estimated that there are sources like Gaia BH1 that go unflagged as binaries in the full Gaia database.
To test our pipeline's ability to distinguish between trajectories of sources with an unresolved binary companion and true blips, we follow a test procedure similar to the one carried out in ref. [19]; that is, we generate a mock catalog consisting of 10 3 luminous stars with masses of either 1 M or 10 M , each with dark companions with masses corresponding to brown dwarfs, white dwarfs, neutron stars, or BHs (0.05 M , 0.6 M , 1.4 M , and 10 M , respectively). Note that the light centroid of two luminous but unresolved stars will follow a trajectory similar to that of a star with a non-luminous companion, like the ones we sample over here. We place these companion objects at distances of 10 pc, 100 pc, and 1 kpc. For each of these combinations of masses and distances, we probe orbital periods of 10, 10 2 , 10 3 , and 10 4 days, with the binary eccentricity drawn from a uniform distribution ranging from 0 to 0.95 and orbital Euler angles drawn from a uniform distribution ranging from 0 to 2π. We then fit our free model, acceleration model, and blip model to the resultant stellar trajectories. We find that there are two classes of binaries, depending on which cuts are passed and which are failed. The first type of binary has an orbital period t bin longer than Gaia's observation time (t bin t obs ). These binaries can have significant free model log likelihoods, but their significance becomes much smaller when fit to the acceleration model due to their trajectory being well approximated by a star undergoing constant angular acceleration in a single direction. In our grid catalog, all of these binaries have acceleration fit log likelihoods below the 3σ interest threshold -2 log Laccel < χ 2 3σ . See figure 13 for an example of a fit of this type.
Figure 13. Residuals for free and acceleration models fitted to the trajectory of a source with a BH binary companion and orbital period t bin = 10 4 days. Left: In red, the residual from the free model fit to the source trajectory, with the 1σ and 2σ bands being shown in green and yellow, respectively. Right: In purple, the residual from the acceleration model fit to the source trajectory. Note that the dashed lines do not show the continuous trajectory of the model and data, but rather simply connect the data points since their order can otherwise be hard to gauge. The free fit exceeds the 5σ free log likelihood cutoff; however, it is far below 3σ significance in the acceleration fit, meaning it fails to qualify as a blip (even without accounting for its associated constrained test statistic TS * . The acceleration fit and cut effectively eliminates sources that are part of long-period binaries. ƭ
The second type of binary has a period comparable to or smaller than the observation time (t obs t bin ) and typically has a significant free log likelihood -2 log Lfree > χ 2 5σ , as well as a significant acceleration log likelihood -2 log Laccel > χ 2 3σ . However, because blips and short period binaries have very distinct trajectories, for most of the sources, TS * < 100. Binary trajectories are also disfavored by the priors we use to constrain the computation of TS * . However, for two sources in the catalog, even this cut is surpassed. To avoid accidentally flagging sources with dark companions as blips, we thus impose the cut TS * < TS 3σ , where the TS 3σ is the 3σ significance threshold for the blip model unconstrained TS distribution obtained via MC generated blip events. All six of the blip events in section 5.2 pass this cut. The actual number of dark companion that Gaia expects to see is much smaller than the number we have considered here, so it is likely that this extra cut is unnecessary. We nevertheless implement it into the analysis pipeline as a precautionary measure. Finally, we note that binaries that are completely dark, e.g. consisting of two black holes, are indistinguishable from isolated black holes in the sky when the binary angular separation is much smaller than the Einstein radius. These dark binaries have never been directly observed and are thus another interesting lens population to probe.
Figure 14 shows the projected constraining power of Gaia DR4 on compact DM, following the procedure of section 4.2. To arrive at this result, we inject 7.0 × 10 9 compact DM objects into the mock catalog (corresponding to 10% of total DM mass for 1 M compact objects). The blue curve shows the resulting 90%-CL limits on mock simulations with delta-function compact DM object mass functions over the range 10 -1 -10 5 M . For compact objects lighter than 0.3 M , there is no event in the signal region, so we are only able to quote an upper bound on f l as shown by the blue arrow. The sensitivity peaks at compact object masses between 1 M and 100 M . At smaller masses, the sensitivity sharply decreases due to the JCAP07(2023)037 saturation of astrometric deflection at the Einstein radius, while for larger masses it decreases more gradually due to the smaller expected number of compact objects with a large blippiness. Existing constraints from photometric microlensing [65,118], dwarf galaxy heating [119], and CMB spectral distortions (from X-ray accretion onto PBHs, not applicable for non-PBH compact objects) [120] are shown in gray.
We also show in figure 14 the initial analytic estimate from ref. [40] for the potentially accessible parameter space of compact DM objects (red dot-dashed curve). At the low-mass end, their estimate is a contour for which the local signal-to-noise ratio equals unity. Without any additional input from other surveys to identify potential astrometric lensing candidates, the look-elsewhere effect and the requirement of setting a 90%-CL limit drastically reduces the projected constraints on the DM fraction at low masses, equivalent to setting SNR = 15 in the language of ref. [40]. The requirement of such a high threshold for a blind search furthermore means that the weak lensing approximation no longer holds, further suppressing the sensitivity of a blind search purely based on astrometry alone. In appendix F, we recalculate the analytic estimate following the same procedure in ref. [40] with the above-mentioned effects and arrive at the updated analytic estimate shown in the red solid curve, which is much closer to the mock catalog simulation. We also show the projected reach of a futuristic 10-year mission with astrometric uncertainties 10 times better than the Gaia EDR3 uncertainties in solid orange. The contour of 2.3 detectable events (corresponding to a 90% constraint) from ref. [76] is shown by the dashed green curve. The difference between our work and that of ref. [76] can also be partially ascribed to differences in treatment of the look-elsewhere effect. The scaling difference at large compact object masses is because that we conservatively discard events that have an acceptable (within 3σ) 7-parameter acceleration fit, necessary to eliminate backgrounds from long-period binary systems.
Precision astrometric measurements from Gaia enable a new way to probe the MW for transient astrometric lensing caused by massive non-luminous objects of either astrophysical or primordial origin, with potential for discovering several free-floating BHs and searching for compact objects down to a very small fraction of DM. We construct an analysis pipeline ( ) capable of systematically and exhaustively searching for transient astrometric lensing events (or "blips") in the upcoming Gaia DR4 catalog. This pipeline works by first fitting a simple free (unlensed) model of stellar motion to more than a billion stars in DR4 using a combination of traditional optimization and bayesian inference. It then discards all events that are not more than 5σ outliers under the free stellar motion hypothesis. To account for binaries, the pipeline then fits a model of stellar motion in which the source being studied undergoes constant angular acceleration. Events that are not more than 3σ outliers under this constant angular acceleration hypothesis are similarly discarded. Finally, the pipeline fits a blip model, weighted by priors on lens proper motion, distance, and mass, to the remaining events. Any events that pass the free fit and acceleration fit cuts and that have blip test statistics TS * > 100 are flagged as blip candidates. Using the Yellin method, the pipeline furthermore infers constraints on dark compact object populations based on the test statistic distribution.
To test the pipeline, we create three types of mock DR4 catalogs based on the currently available EDR3 catalog. The first contains no dark lenses, meaning all sources undergo free stellar motion. In this catalog, the pipeline flags no events as being blips, and the log JCAP07(2023)037 The blue curve is the projected 90% constraint of DM fraction f l in the form of compact objects from the analysis in this work, assuming no other astrophysical backgrounds. Our sensitivity is peaked around 10-100 M and sharply evaporates below 1 M because there the Einstein radius is smaller than the astrometric precision of Gaia. At larger masses, the sensitivity to f l decreases linearly due to the decrease in lens number density at fixed f l . We overlay the analytic SNR = 1 estimate (dot-dashed red curve) of ref. [40] and our updated analytic estimate for the 90%-CL exclusion limit for a blind astrometry-only analysis (solid red). The reach of a hypothetical 10-year future mission with 10× better astrometric uncertainties is shown in the orange solid curve using the same analytic estimates. The 90%-CL exclusion curve using the probabilistic model of ref. [76] is depicted as the green dashed curve. Existing constraints from MW photometric microlensing [65,118], dwarf galaxy heating [119], and CMB spectral distortions from PBH accretion [120] are shown by gray shaded regions. ƭ likelihood distribution follows the χ 2 expectation (see figure 9). The second mock catalog is identical to the first, except we inject astrophysical BHs based on current priors on the BH number density and proper motion distribution across the MW. In this catalog, we find 4 lensing events that pass all of our cuts; namely, they are above 5σ significance under the free model expectation and above 3σ significance under the acceleration model fit, which separates the events from long-period binary systems with a dark companion, and they have a constrained test statistic TS * > 100. This gives us a benchmark of the total number of astrometric lensing events by isolated astrophysical BHs we expect to discover in Gaia DR4.
We inject the third mock catalog with compact objects of a single mass spanning the range 10 -1 -10 5 M to constrain their fraction of DM in the MW using the Yellin method. Our projected constraint indicates that Gaia has leading reach on the compact DM fraction JCAP07(2023)037 in the mass range of 1-10 3 M . We find that Gaia loses sensitivity for point-like DM lenses lighter than 0.3 M , is most sensitive between 10-100 M (projected exclusion fraction of f l ∼ 4 × 10 -3 ), and runs out of observable blip events for higher masses as the number of lenses and thus transient lensing events decreases. Our full Gaia DR4 mock catalog enables us to properly assess the statistical background of the large data set to obtain faithful projections of discovery potential and constraints.
We make a few assumptions and simplifications in creating the mock Gaia DR4 catalog which will be different from the actual Gaia DR4. Here we outline those points and the potential effect on the actual data analysis with real Gaia data.
• We assume all sources in Gaia will be observed exactly 80 times, roughly the skyaveraged expected number of observations. This is not the case for the real data. Each source will be observed roughly 60-140 times depending on the source's ecliptic latitude.
If the high-cadence region has a larger/smaller overlap with the region of higher stellar density (e.g. Galactic plane), then we would expect more/fewer lensing events discovered compared to the mock catalog.
• We assume Gaia only records the one-dimensional offset along the AL direction for all stars. This is not true for the brightest stars. They will have the full two-dimensional trajectory in the AL and AC direction recorded. However, the uncertainty in the AC direction is orders of magnitude worse than that of the AL direction due to design of the telescope. This will only improve sensitivity of the brightest stars by a small margin.
• The AL uncertainties we adopt in the mock catalog are the projected optimal uncertainties of DR4 reported in Gaia EDR3. If the actual uncertainties are different, the sensitivity projections in this work will be affected accordingly.
• We only inject astrophysical BHs for our search for compact objects. In reality, there will be other compact objects, such as neutron stars, white dwarfs, brown dwarfs, and faint main sequence stars. These objects could affect our projection, although we argue that their effect will be marginal (see appendix B). Potential contamination due to these other sources must nevertheless be carefully accounted for when DR4 is released and real data is available.
• We only use the effects of astrometric lensing for finding compact lens in this work. Gaia DR4 will also release time-series photometric measurements of the stars. Although Gaia's photometric capabilities are not optimal for lensing searches, a combination of its photometric and astrometric measurements will likely lead to more precise lens parameters and potentially stronger discovery potential, especially for low-mass lenses for which strong lensing events are more common.
Beyond the single-source blip search outlined here, it is also interesting to consider events in which a non-luminous lens affects the astrometric trajectory of multiple sources in a short time interval. Such events may not be detectable by probing for solitary blips, since the lensing deflection of any given source might be too small to be statistically significant. Furthermore, observing two or more sources undergoing gravitational lensing due to the same lens would likely yield a much better determination of the physical parameters of the lens. Conventional likelihood optimization, as used in this work, is likely not computationally feasible for carrying out a "multi-blip" search due to the number of free parameters in such a JCAP07(2023)037 model. Machine learning tools will likely accelerate the pattern recognition of those correlated lensing deflections -an avenue we will explore in future work.
Our analysis pipeline and mock catalog are not just applicable to Gaia DR4. The tools we provide in this work can be used on past astrometry legacy archives (e.g. HSTPROMO [74], PHAT [75]), as well as future astrometric surveys (e.g. the Nancy Grace Roman Space Telescope (formerly known as WFIRST) [121], GaiaNIR [122,123], THEIA [124] with minor adjustments, and of course Gaia DR5, which is projected to contain all collected Gaia data [122]). Charting out several isolated, electromagnetically quiet BHs will be a major milestone in astrophysics, and help in the understanding of their formation mechanisms. Finally, isolated BHs are also pristine laboratories for Beyond the Standard Model Physics searches. The extreme gravity near a BH can give rise to BSM signals, most notably through superradiance [102][103][104][105][106][107].
Transient astrometric weak lensing is a powerful probe of the distribution and properties of known compact remnants, such as BHs and neutron stars, as well as extreme overdensities in the DM distribution. We look forward to the application of our tools to these studies.
For simplicity of calculation, we assume the DM subhalo has a Gaussian density profile given as
where M l is the mass of the lens and r l is the scale radius of the lens. We define the mean lens density as ρ l,0 = M l /r 3 l . The resultant relation between r l and ρ l,0 is shown in figure 15. The solid line is the contour of total lens mass. The turning point near large scale radius is where the scale radius is equal to the Roche radius of the MW at 8 kpc. The dashed-dotted line is the contour of the maximum deflection an extended lens can induce. We require the deflection to be larger than 10 µas to be detected by Gaia, so extended lenses in the red-shaded region are not detectable. Another criterion for detection is eq. (4.5). Given that the maximum deflection of a extended lens occurs at the scale radius, the blip criterion is v rel τ ≥ r l . In the most conservative case where v rel = 1000 km/s, which corresponds to 2 objects moving back-to-back both at the galactic escape velocity, the requirement on r l is shown as the horizontal blue-dashed line.
The density of a subhalo that collapses at matter-radiation equality assuming a pure ΛCDM cosmology is shown as the vertical dotted line, marking the maximum density of a gravitationally collapsed subhalo in a standard cosmology. Since it is not within the range of the detectable parameter space, we only include lensing from point sources in this work.
The analysis pipeline we present in the main text is also suitable for carrying out a blip search on compact lens populations other than BHs. Here, we provide preliminary estimates of the blip discovery potential of neutron stars, white dwarfs, brown dwarfs, and faint main sequence stars in Gaia DR4.
The MW is estimated to contain 10 8 -10 9 neutron stars [132], which is 1-10 times the total number of BHs we inject in our mock catalog. The mass distribution of neutron stars is believed to lie within 1.0-2.2 M , peaking at 1.4 M [133]. Observations of neutron stars suggest that they, like BHs, receive natal kicks from supernovae, explaining their high velocities and large fractional abundance in the stellar halo [132]. Therefore, we assume that the spatial and velocity distribution of neutron stars follow that of BHs.
We may thus use our mock analysis of astrophysical BH lensing from section 5 to extrapolate the expected number of neutron star lensing events we will see in Gaia DR4. Astrophysical BHs typically have a mass of about 10 M and neutron stars typically have a mass of around 1 M . From figure 14, we see that the sensitivity from 10 M to 1 M drops by a factor of ∼ 5. We assume there are 10 9 neutron stars in the MW. Extrapolating, this means that the number of neutron star lensing events in DR4 with a significance level above 5σ is approximately one.
The MW is estimated to contain approximately 10 10 white dwarfs [134], which is 100 times the total number of BHs contained in our mock catalog. The mass distribution of white dwarfs covers a range of 0.4-1.4 M , peaking at 0.7 M [135]. Observations of white dwarfs
This figure shows the detectable region for an extended lens following a Gaussian density profile eq. (A.1). The gray solid line shows the contour of constant lens mass with the turning point being the Roche radius at 8 kpc. The red dashed-dotted line shows the maximum deflection an extended lens can cause (which equates to the minimum impact parameter being the scale radius). We exclude the parameter space of δ max < 10 µas which is our fiducial value of the Gaia DR4 sensitivity. The blue dashed line and the shaded region show the blip requirement for a source-lens relative velocity of 10 3 km/s, which corresponds to the lens and source both moving at the galactic escape velocity back-to-back. The vertical black dotted line shows the density of a subhalo that gravitationally collapsed at matter-radiation equality. ƭ combined with simulations suggest that white dwarfs can be categorized into three families based on their kinematics: thin disk, thick disk, and halo [134]. Each of the three families consist of O(1) of the total number of white dwarfs in the MW. To understand the discovery potential of isolated, faint white dwarfs in DR4, we follow the same procedure as in section B.1. Namely, we use the results of our BH mock analysis to extrapolate. Figure 14 shows that the sensitivity drops sharply for lenses with a mass less than 1 M . Thus, we are only sensitive to white dwarfs with a mass greater than 1 M . Ref. [135] suggests that roughly 10% of white dwarfs fall beneath this mass cutoff. So with 10 10 white dwarfs in the MW, the number of potential observable white dwarfs is 10 9 . This suggests that the number of white dwarf blip events in DR4 with a significance level greater than 5σ is approximately one.
Brown dwarfs are stellar objects with masses in the range 13 M J -80 M J (1.2-7.6 × 10 -2 M ), where M J is the mass of Jupiter. This is the mass range in which a star burns deuterium and hydrogen. Using the projected compact DM constraint shown in figure 14, we can see that the mass of a typical brown dwarf lies below Gaia's detectable range. This suggests that we will not see blip events caused by any isolated brown dwarfs in Gaia DR4. One can also see this by using the analytic SNR estimate described by eq. (F.5). For a brown dwarf with a mass of 5 × 10 -2 M located 10 (100, 1000) pc from the Sun, the maximum SNR one can get from astrometric lensing is 7 (4, 2), which is smaller than the SNR = 15 threshold. Thus, photometric microlensing is more suitable for the detection of brown dwarfs, cfr. the shaded gray region of figure 14.
Main sequence (MS) stars are another possible lens population. A MS star passing in front of a background star can cause a blip event. Refs. [136,137] discuss signals of star-star lensing and how to detect them in Gaia data when both the lens and the background star are above Gaia's photometric threshold. They propose that star-star lensing can be used to determine the mass of luminous foreground stars. Here, we discuss blip events caused by faint MS stars dimmer than the Gaia photometric threshold (G ≈ 21).
As discussed in section 5.2, all 5σ stellar BH events are within 1 kpc. This is due to the blippiness requirement described by eq. (4.5), which gives preference to lenses with large proper motions typically located at small line-of-sight distances. A MS star 1 kpc away from the Sun with an apparent magnitude of 20 will have an absolute magnitude of 10. Using the mass-luminosity relation
we estimate that the mass of such a star is roughly ∼ 0.3 M . Any MS star closer than 1 kpc that is too faint for Gaia to detect must be lighter than this, which places the star outside Gaia's projected blip sensitivity shown in figure 14. Therefore, our preliminary analysis using mock catalogs suggests that Gaia is not capable of discovering blips caused by faint MS stars.
Here, we discuss how using Gaia's Observation Forecast Tool (GOST) affects the results presented in section 5. To obtain accurate time-series data points for each Gaia source, we compute the average angular location of each local batch of sources (with each batch corresponding to one of the 3386 Gaia EDR3 files) by taking an average of their HealPIX location. Inputting this into GOST, we obtain the scanning law associated with each source. Using this method, the number of observations per source ranges from 43 to 249. The location dependence of Gaia's observation cadence is shown in figure 16. We rerun the analysis described in section 4 on mock catalogs generated using GOST. Limits obtained from the GOST DM mock catalog are shown in figure 17. We note that these limits are marginally weaker than those obtained using the 80 data points scenario. Furthermore, we also conduct a BH search on a BH GOST mock catalog and find exactly 3 (6) highly significant sources with (without) the acceleration test statistic cut. This result corresponds almost exactly to the one obtained using the 80 data points method presented in the main text.
Galactic GOST Number of Observations
Starting with the thin disk stellar distribution in eq. (3.2), we can estimate the increase in z d by considering the following. We assume all stars start at exactly z = 0 with some known velocity dispersion σ vz . The probability distribution function (PDF) of stars at z = 0 is
From energy conservation, the PDF of stars at z is
where φ(z) is the gravitational potential at z. Marginalizing over velocities gives:
where the second equals sign comes from eq. (3.2). Here we can see that if the background gravitational potential stays the same, the scale height z d ∝ σ 2 vz . BH X-ray binaries (figure 7 in ref. [138]) suggest a bimodal distribution of natal kick velocities. In Gaia DR2, the vertical velocity dispersion around the solar neighborhood is reported to be around σ vz ≈ 20 km/s [139]. Combining the stellar velocity dispersion and natal kick, the final velocity dispersion is approximately σ vz ≈ 70 km/s. In terms of the scale height of the thin disk distribution, this implies that the scale height of BH distribution is around 10 times that of the scale height of stellar distribution. Therefore, we use z d = 3 kpc for the BH distribution in the sky. (D.4)
The distance prior P (D l |l, b) is given by eq. (3.3). The conditional probability P (µ µ µ l |D l , l, b) can be calculated via the following process: we start with the conditional probability
where v C is the linear velocity vector in a cylindrical coordinate centered at the galactic center and with φ = π pointing towards the solar system. Σ Σ Σ is the velocity dispersion of the lens and we assume it is diagonal in this coordinate system. Next, we can rotate this into a Cartesian coordinate (U, V, W ) commonly used in astronomy where the galactic center sits at (0, 0, 0), the solar system sits at (-8, 0, 0) kpc, the V axis points towards the direction of the Sun's orbit around the galactic center, and the W axis points towards the galactic north pole.
And, shift into a frame where the Sun is stationary. Then, the joint PDF in the Cartesian coordinate is Here v R is the linear velocity relative to the Sun in the Cartesian coordinate, v R is the Sun's velocity, and φ is the angle in the cylindrical coordinate. Then, we can rotate from the Cartesian coordinate to galactic coordinate (r, l, b) Here v G is the linear velocity in galactic coordinate. One more rotation brings the velocity into equatorial coordinate (r, α, δ) [139] and the bimodal distribution of BH natal kicks from ref. [138]. ƭ
We can calculate the total magnification of a point-like background star due to a point-like foreground lens by summing up the magnification of the two images in eq. (2.7) when the two lensed image are not resolved independently using (2.9). Using this equation, we calculate the light curves of the 6 candidate lensing events found in our mock catalog (see section 5.2). These light curves are shown in figure 19. The horizontal black dashed line is the photometric uncertainty per transit for each of the background stars taken from ref. [67]. The Einstein radius and the minimum dimensionless impact parameter for each event in shown in the plot as well. We can see that the maximum brightening is below the Gaia sensitivity for all but 2 (5727504125199235456, 6262458554071571712) sources. And for the 4 sources that pass the 3σ acceleration fit, only one (5727504125199235456) has magnification larger than the Gaia photometric uncertainty.
The Gaia Photometric Alert System [140] will likely discover lensing events similar to 5727504125199235456. In fact, there is already a successful detection of a lensing event in Gaia with ID Gaia16aye [141]. This event was flagged using the Gaia Photometric Alert System and later confirmed to be a foreground lens consisting of a binary star system via Gaia astrometry coupled with ground-based photometry follow-up. However, this work demonstrates that many significant lensing events will evade a photometric alert system and only be detectable via astrometry.
Suppose that stars in the Gaia catalog are distributed evenly and are stationary at infinity. A lens with velocity v will sweep through an area of 2vτ b min . Thus, the expected minimum impact parameter of all lens is
where N * is the number of stars in the Gaia catalog. For this event to be a blip we require that b min < vτ . Plugging in N * = 1.4 × 10 9 , ρ DM = 10 -2 M pc -3 , D l = 10 kpc, we arrive at the rightmost branch of the analytic estimate:
On the other end, the ∆χ 2 used for the event selection is a proxy of SNR 2 , which can be parameterized by For the SNR to reach some threshold, we then arrive at the expression:
Accounting for the look-elsewhere effect and the average lens distance for significant events, we use SNR = 15 and D l = 1 kpc, which yields the left branch of the red dashed-dotted analytic estimate in figure 14, closer to the simulation done in this work. For strong lensing that saturates the astrometric deflection, eq. (F.3) is modified as The expected distance to the closest lens D l can be expressed as
Plug this back into eq. (F.5) to get the sharp cutoff in the left branch of the red-solid curve in figure 14.
Another thing we discovered is that eq. (F.3) only applies when the blippiness is large ( 10) because of the definition of δ max , which should be δ max -δ min for calculating ∆χ 2 . For events with large blippiness, δ min ≈ 0 so eq. (F.3) is valid. However, as figure 20 shows, events with small blippiness ( 10) do not follow this relation, becoming almost independent of blippiness, which we parametrize as the following: where 0.4 is the peak of the blue curve in figure 20. This gives the constraint:
which is the horizontal branch of the red-solid curve in figure 14.
JCAP07(2023)037 The relationship between SNR to blippiness. The red curve shows the relation adopted in ref. [40] as shown in the left branch of the red dashed-dotted curve in figure 14 and in eq. (F.4), which uses the maximum deflection δ max for calculating SNR. The blue curve shows the relation using the difference of maximum deflection and minimum deflection throughout the mission time δ max -δ min for calculating SNR, as is the relation used in the horizontal branch of the red solid curve in figure 14 and in eq. (F.7). Here we can see that the scaling changes for blippiness 10 and the SNR remains approximately constant in thie regime. ƭ
The Gaia's Observation Forecast Tool (https://gaia.esac.esa.int/gost/) provides a forecast of Gaia observations and scan angles.
One of the derived end products of our data analyses on Gaia DR4 and other data sets will be to pin down the mass function for isolated astrophysical BHs.
http://healpix.sourceforge.net.