Microlensing: planets that bend starlight

In the southern winter of 2003 two survey telescopes were watching the same star in the Galactic bulge brighten. The Optical Gravitational Lensing Experiment (OGLE) had flagged it on 22 June as a microlensing event: a star in front of it, too faint to see, was acting as a lens and slowly magnifying it. On 14 July the light curve stopped following the smooth, symmetric shape that a single lens produces. For about a week it ran high, and on 21 July the 0.6-m telescope of the Microlensing Observations in Astrophysics (MOA) collaboration at Mount John, New Zealand, caught a sharp spike and then a sudden drop. Three groups modelled the 1,275 measurements independently, with different codes, and reached the same answer. The lens was a star with a companion 0.0039 times its mass: a planet of roughly one and a half Jupiter masses, several kiloparsecs away, which neither group could see and never would (Bond et al. 2004).

Nothing about that planet was luminous. It was detected because its gravity changed the paths of photons from a third body, a background star, that happened to pass behind its own host star. This is the only exoplanet method in this part of the book that does not depend on light from the planet or its host. It finds cold planets near and beyond the snow line, planets around faint stars far across the Galaxy, and planets with no host at all. It is also the method in which a discovery is least repeatable: the alignment happens once, the anomaly lasts hours or days, and the light curve is a single one-dimensional cut through a complicated two-dimensional pattern of magnification.

Those facts shape the computing. Microlensing has three separate search problems, and it pays to keep them apart. Event finding screens hundreds of millions of stellar light curves for the brightening caused by any lens. Anomaly finding searches the residuals of known events for the extra bump or dip a companion adds. Modelling explores a strongly multimodal parameter space to decide whether that companion is a planet, and whether some other configuration explains the data equally well. A fourth problem, completeness, repeats the second and third thousands of times on injected signals to turn a list of discoveries into an occurrence rate. This chapter follows those problems in order.

From a lensing experiment to a planet experiment

The idea

Shude Mao and Bohdan Paczyński’s 1991 Letter made the step from lensing to planets. Binary lenses, they argued, should be recognisable in monitoring of the Galactic bulge, and the perturbation stays strong even when the second mass is a planet. They proposed discovering planetary systems this way before any microlensing event had been confirmed. Their examples used a foreground stellar lens and a background bulge source, computed the probability that the source would cross a caustic, and stressed that short spikes would need real-time processing and advance warning. Their planetary disturbances lasted 2.5 to 10 hours for the event timescales they adopted. These were model predictions, but they set the cadence target for the next three decades (Mao & Paczyński 1991).

Andrew Gould and Abraham Loeb then worked out the observing strategy (Gould & Loeb 1992): find ordinary events with a wide-field survey, alert other observers while the events are still rising, then monitor selected stars rapidly from telescopes spread in longitude so that the bulge is never out of sight. OGLE and MACHO ran early alert systems. Follow-up networks, PLANET, GMAN, MPS and later μFUN, supplied the dense monitoring (Kim, D.-J. et al. 2018; Kim, H.-W. et al. 2018). Kim Griest and Neda Safizadeh added the observation that made follow-up efficient. In a high-magnification event the source passes very close to the host lens, where every planet, whatever its separation, creates a small central caustic. Sensitivity to planets is therefore concentrated in a short, predictable interval around peak, which is exactly when follow-up telescopes should be pointed at the event (Griest & Safizadeh 1998).

The surveys built the data stream. OGLE began in 1992, and its dedicated 1.3-m Warsaw telescope at Las Campanas enabled the OGLE-II phase from 1997. Difference imaging became central with OGLE-III in 2001, because the crowded bulge defeats aperture photometry of isolated stars (Udalski et al. 2015). The technique itself belongs to Chapter 10. What matters here is that it turned the bulge into a stable, dense set of light curves in which a planetary anomaly could be seen.

The first planet and the first warning

That last number is the thread through this chapter. The first secure microlensing planet already came with an alternative interpretation that the data did not decisively exclude, and with a physical mass that depended on a model of the Galaxy. Neither problem was introduced by later automation. Both are built into the method.

Low masses, survey networks and buried planets

OGLE-2005-BLG-390 showed the method could reach cold, low-mass planets. OGLE alerted the event on 11 July 2005. After peak magnification on 31 July, the PLANET network detected a deviation on 9 August. The combined data gave $q = (7.6 \pm 0.7) \times 10^{-5}$ and an inferred companion of $5.5^{+5.5}_{-2.7}\,M_\oplus$. The often quoted orbital separation of about 2.6 AU is a Bayesian interpretation, not the measured dimensionless separation. The source was a giant star, and its size smoothed the anomaly so much that the planet sat near the detection floor for that source (Beaulieu et al. 2006).

In the 2010s the observing model changed. OGLE-IV’s large mosaic, MOA-II’s high cadence and the Korea Microlensing Telescope Network (KMTNet) replaced scarce follow-up time with continuous survey coverage. KMTNet has three 1.6-m telescopes, each with a 4-deg² camera, at Cerro Tololo in Chile, the South African Astronomical Observatory and Siding Spring in Australia, so that the bulge can be watched around the clock. Commissioning observations began in 2015. The event-finding software for that season was developed only from mid-2016, after the data were already accumulating, and found 660 clear and 182 possible events (Kim, D.-J. et al. 2018). Telescope commissioning and discovery-software readiness are different milestones.

With continuous coverage, many planets could be found after the fact, in archival data, by software. Weicheng Zang and colleagues introduced the AnomalyFinder in 2021. Its first demonstration uncovered OGLE-2019-BLG-1053Lb, a low-mass-ratio planet whose signal was buried in a peripheral part of the light curve (Zang et al. 2021). Kyu-Ha Hwang and colleagues then found six new planets with $q < 2 \times 10^{-4}$ in the 2018–2019 data, five of them as subtle dips that visual searches had missed (Hwang et al. 2022). The systematic search has since been completed for the first four KMTNet seasons, 2016 to 2019, with 112 unambiguous planets, including previously known ones (Gui et al. 2025), and extended to the 2021 prime fields, where it found seven hidden planetary systems and three candidates (Shin et al. 2026). These planets were in the archive all along. Better software, not better telescopes, found them.

The physics of a lens

Einstein scales

A point mass $M$ at distance $D_L$ deflects light from a source at distance $D_S$. When the three are aligned, the source appears as a ring of angular radius

the Einstein radius. Here $\pi_{\rm rel} = {\rm AU}\,(D_L^{-1} - D_S^{-1})$ is the relative parallax of lens and source and

The physical Einstein radius at the lens is $R_E = D_L\theta_E$, and the event timescale is $t_E = \theta_E/\mu_{\rm rel}$, where $\mu_{\rm rel}$ is the relative proper motion of lens and source (Penny et al. 2019).

Two features of these numbers explain where microlensing looks. First, the Einstein radius of a typical bulge lens is a few AU, near the snow line of a low-mass star. A planet is most likely to perturb a magnified image when its projected separation $s$, in units of $\theta_E$, is close to 1, so microlensing is most sensitive to planets a few AU from their hosts: cold planets that transits and radial velocities find hard. Second, the planetary signal is short compared with the event. The survey has to sample a 50-day event densely enough to catch a 3-hour feature in it.

The single-lens light curve

Write $u$ for the lens–source separation in units of $\theta_E$. For a point source, a single point lens and straight-line relative motion,

where $u_0$ is the impact parameter and $t_0$ the time of closest approach. For small $u$, $A \approx 1/u$. This is the Paczyński curve, and it is the shape every event-finding algorithm looks for. Each telescope and filter $d$ records

$F_d(t) = F_{s,d}\,A(t;\boldsymbol\vartheta) + F_{b,d},$

where $F_s$ is the source flux and the blend flux $F_b$ includes unresolved neighbours and possibly the lens itself. In a difference image one measures approximately $F_s[A(t) - A(t_{\rm ref})]$, the change since the reference image, with an offset that depends on the dataset (Zang et al. 2021; Koshimoto et al. 2023).

The two fluxes enter linearly. For fixed nonlinear parameters $\boldsymbol\vartheta$, with design matrix $X_d = [A, 1]$, data $\mathbf y_d$ and inverse covariance $W_d$, the best fluxes are a weighted least-squares solution,

$\widehat{\boldsymbol\beta}_d = (X_d^T W_d X_d)^{-1} X_d^T W_d\, \mathbf y_d .$

Every serious code profiles them out this way, removing two search dimensions per dataset. That matters more than it sounds: a modern event may have data from six or more site–filter combinations, so a naive search would carry a dozen extra parameters. With independent Gaussian errors the objective is the familiar $\chi^2$. With correlated errors it is $\mathbf r^T C^{-1}\mathbf r$, plus $\log\det C$ in the likelihood. In practice, bulge photometry is not well described by its nominal errors. Blending, sky, seeing-dependent subtraction residuals and shared reference images make errors correlated and heavy-tailed. The usual fix is an empirical rescaling, $\sigma' = k\sqrt{\sigma^2 + e_{\min}^2}$ in magnitudes, which has to be converted consistently when fitting in flux (Zhang et al. 2021).

Two masses and their caustics

A planet makes the lens a binary. In complex coordinates $z$ (image plane) and $\zeta$ (source plane), normalised by the Einstein radius of the total mass, with lens components at $z_j$ carrying mass fractions $\epsilon_j$, the lens equation and its Jacobian are

The point-source magnification is the sum over images, $A_{\rm point} = \sum 1/|J|$. For two masses, eliminating $\bar z$ gives a complex polynomial of fifth degree in $z$. Away from caustics only three or five of its roots are real images, and the others must be rejected by substituting them back into the lens equation (Bozza et al. 2018). This fifth-degree polynomial, solved millions of times, is the inner loop of the whole field.

Where $J = 0$ the magnification of a point source is infinite. These points form closed critical curves in the image plane, and their images in the source plane are the caustics. When a point source crosses a fold caustic, a pair of images appears or disappears, and the magnification on the inside diverges as the inverse square root of the distance from the fold. A planet produces two kinds of caustic. A small central caustic sits near the host star, and is what a high-magnification event probes. One or two planetary caustics sit near the position $s - 1/s$ from the host along the star–planet axis. For a wide planet ($s > 1$) there is a single four-cusped planetary caustic on the planet’s side. For a close planet ($s < 1$) there are two small triangular caustics on the far side, with a region of demagnification between them. Near $s = 1$ the caustics merge into a single resonant caustic.

The standard static model of a planetary event, with one source, has seven nonlinear parameters,

where $\alpha$ is the angle of the source trajectory to the star–planet axis and $\rho$ the angular source radius in Einstein units. A convenient trajectory is $\zeta(t) = e^{i\alpha}[(t - t_0)/t_E + i u_0]$, with the coordinate origin chosen consistently with the code. Conventions differ between codes (origin at the primary, at the centre of mass, or at the “centre of magnification”, and Einstein radius of the total or of the primary mass), and copying parameters between codes without converting them is a common way to get a wrong curve.

Two scaling results guide intuition. For $q \ll 1$, away from the resonant regime, planetary caustics have a linear size proportional to $\sqrt q$, while the central caustic’s width along the binary axis is approximately

$w_c \simeq \frac{4q}{(s - s^{-1})^2}$

(Chung et al. 2005). This is asymptotic and fails near $s = 1$, where its divergence is not physical. For $q = 10^{-3}$ and $s = 1.3$ it gives $w_c \approx 0.014$. For an Earth-mass ratio, $q = 6 \times 10^{-6}$, it gives $8.5 \times 10^{-5}$. The anomaly duration is often estimated as $t_p \sim \sqrt{q}\,t_E$, but central-caustic geometry, grazing trajectories and the source crossing time $t_\star = \rho t_E$ can control the observed duration instead. One cannot read $q$ off the duration of an anomaly (Wang et al. 2025).

The figure shows the two basic signatures. A trajectory through the planetary caustic of a wide planet produces a sharp bump, often double, as the source enters and leaves the caustic. A trajectory through the gap between the two triangular caustics of a close planet produces a dip, because the planet sits near the minor image and suppresses it. Dips were the signals visual searches missed most often. Any anomaly search must allow the extra signal to be negative.

Finite sources

A real source is a disk. The observed magnification is the point-source field averaged over it,

where $I$ is the surface brightness, including limb darkening. The convolution regularises the caustic divergence. It also sets a physical floor on detectable planets. A Sun-like star at 8 kpc subtends about 0.58 μas, so for the lens of the worked example $\rho \approx 8 \times 10^{-4}$, comparable to the $\rho = 0.00096$ fitted for the first planet. A giant of ten solar radii subtends about 5.8 μas, larger than the 1.7-μas Einstein radius of an Earth-mass planet. Against such a source an Earth’s caustic is washed out, which is why OGLE-2005-BLG-390Lb, with its giant source, sat close to the floor. Small sources keep fine caustic structure but are faint, so photometric precision and finite-source smoothing pull in opposite directions (Beaulieu et al. 2006). If the light curve changes within an exposure, the model should also average over the shutter interval (Miyazaki & Kawahara 2025).

What is measured, and what is inferred

Lens flux measured years later in high-resolution images, the separation of lens and source as they move apart, and astrometric microlensing can all add constraints (Zang et al. 2021; Penny et al. 2019). Astrometric lensing is treated in Chapter 34. For demographics, the honest coordinates are $(s, q)$, and most occurrence results from microlensing are stated in them.

Detection statistics

The basic test in an anomaly search is a matched filter applied to residuals. Let $\mathbf r$ be the residual vector after the single-lens fit, $\mathbf h$ a trial anomaly shape, and $W$ the inverse covariance, with the constant baseline projected out. The best signed amplitude is $\hat a = (\mathbf h^T W\mathbf r)/(\mathbf h^T W\mathbf h)$, and the improvement in fit is

This is the weighted reduction at the heart of the fast bump-and-dip searches, written here in analytic form; pipelines differ in normalisation. The sign of $\hat a$ must be kept: a bump and a dip are different physical signals. The test that matters for a discovery compares optimised physical models,

$\Delta\chi^2_p = \min\chi^2_{\rm 1L1S} - \min\chi^2_{\rm 2L1S},$

the improvement of a binary-lens, single-source model over a single-lens, single-source model. Confirmation also requires comparison with a stellar binary lens and with a binary source: a faint second source star magnified by a single lens can mimic a short positive anomaly.

None of these numbers is a significance in the textbook sense. Under the null hypothesis of no planet, $s$ and $\alpha$ are not identified, the search maximises over many correlated templates, and pipelines often trim outliers adaptively. The tail probability of $\Delta\chi^2$ is therefore not the $\chi^2$ distribution, and nominal $\sqrt{\Delta\chi^2}$ is not a survey-wide false-alarm rate (see Chapter 1 on trials factors). The first planet paper already documented non-Gaussian outliers in its photometry (Bond et al. 2004). The field’s working answer is empirical: high thresholds, real-data null tests, inspection of the images, agreement between independent observatories, and injection and recovery.

Methods

The pipeline has six stages, and a gain in any one is worth only what the stages around it can use.

From crowded images to searchable light curves

Everything starts with difference imaging: register each image to a sharp reference, match the point-spread function with a convolution kernel, subtract, and measure flux at catalogue positions in the difference image (Chapter 10). A search in the image domain also finds sources too faint to be in the reference catalogue. Quality metadata (seeing, background, saturation, detector position and subtraction diagnostics) must travel with every measurement, because they are what later rejects false anomalies (Udalski et al. 2015; Koshimoto et al. 2023).

The scale is large before any lens is modelled. OGLE-IV reported 100 to 200 GB of uncompressed images per night. From 2010 to 2014 its photometry ran on 48 CPU cores and finished 2 to 16 hours after the night. A 260-core server installed in 2015 brought reductions to within minutes of preprocessing (Udalski et al. 2015). MOA’s nine-year analysis used 2,489,362 CCD images, equal to 248,936 full-camera exposures over 3,146 days, from its 1.8-m telescope with a 2.18-deg² field. Six fields were observed every 15 minutes from 2008 to 2014. Its image-domain search kept difference-image peaks above S/N 2.7, then required repeated detections at a consistent position and further light-curve cuts. The 2.7 is a candidate threshold, not a detection threshold (Koshimoto et al. 2023).

The most instructive recent result is about photometry rather than models. Qiyue Qian and colleagues rebuilt KMTNet full-frame subtraction with improved pySIS-based processing and Gaia input positions, then searched 483,068 effective bright-star light curves. They found 36 microlensing events, including five new events and six known to other surveys but missed by earlier KMTNet searches. Registration and subtraction cost about 50,000 CPU-hours and light-curve extraction about 3,500 CPU-hours for the tested data. Their extrapolated cost for all eight seasons is 16.3 million CPU-hours, before review and modelling (Qian et al. 2025). That is a projection, not a completed search. Its lesson is that a signal never represented accurately in the light curve cannot be recovered by any downstream optimiser, however good.

Finding events: EventFinder

A direct search for Paczyński curves would grid $(t_0, t_E, u_0)$: a three-dimensional template bank for each of hundreds of millions of stars. Dong-Jin Kim and colleagues reduced it to two shapes. With $Q = 1 + [(t - t_0)/t_{\rm eff}]^2$,

$H_0$ is the high-magnification limit of the Paczyński curve, with $t_{\rm eff} \simeq u_0 t_E$, and $H_1$ is the curve for $u_0 = 1$. The fitted $f_0, f_1$ are search parameters, not physical fluxes. Collapsing $u_0$ to two shapes is the key approximation (Kim, D.-J. et al. 2018).

for each star:
    for t_eff in geometric grid (ratio 4/3):
        for t0 spaced by t_eff/3:
            window = epochs within t0 ± 5 t_eff   (need ≥ 50 points)
            for shape in {H0, H1}:
                weighted linear fit for (f0, f1)
                drop worst 10% of points; refit
                compare with piecewise-constant baseline
    keep best trial above an empirical threshold
group neighbouring stars with similar events; inspect representatives
confirm with other sites; classify

For the 2015 season, durations spanned about 0.56 to 99 days and fits used only data within $t_0 \pm 5t_{\rm eff}$. For a 250-day season with about 4,000 points per star, the paper estimates about 10,000 trial fits per star. That number is easy to check. With a ratio of 4/3, the duration range needs $\log(99/0.56)/\log(4/3) \approx 18$ steps. The number of centres at each duration is $250/(t_{\rm eff}/3)$, and summing the geometric series over durations gives about 5,300 centres, or about 10,700 fits for two shapes. The catalogue had 71,555,640 stars. The often quoted $10^{12}$ measurements across $3 \times 10^8$ light curves a year is the paper’s scaling target, not the size of its commissioning sample (Kim, D.-J. et al. 2018).

The window is what makes this cheap. Let $S$ be the number of stars, $N$ the epochs per star, $K$ the number of durations and $\delta$ the centre step in units of duration. Each epoch falls in $O(Z/\delta)$ windows at each duration, where $Z$ is the window half-width in units of duration (here 5). Direct fitting therefore costs $O(SNKZ/\delta)$, not $O(SN \cdot N_{\rm templates})$, and the duration grid adds only logarithmically to the dynamic range. This derivation assumes roughly uniform sampling. The score is $N_{\rm data}(1 - \chi^2_{\rm lens}/\chi^2_{\rm flat})$ in positive-improvement form, and the operational thresholds of $\Delta\chi^2 > 500$, later 1,000, refer to that definition, not to a likelihood-ratio tail. Friends-of-friends grouping stops a single bright variable, whose light leaks into neighbouring catalogue positions, from multiplying the review burden.

Alerts and follow-up

A search of completed light curves cannot request more photons once an anomaly has passed. KMTNet’s AlertFinder handles the real-time problem. It first requires enough new points high above baseline, then consecutive high points in combinations of observatories, and finally a broken-line fit,

$F(t) = a + b\,\max(0,\, t - t_{\rm rise}),\qquad b > 0 .$

The published settings use $(N_{\rm high}, \Delta\chi^2_{\rm thresh}) = (5, 250)$ or $(10, 400)$, with high points at least $3\sigma$ above baseline. Using different site combinations prevents one poorly behaved site from vetoing a real rise. Operators inspect what passes (Kim, H.-W. et al. 2018).

For a fixed break time the fit is linear, so prefix sums of $w$, $wy$, $wt$, $wt^2$ and $wty$ allow every break to be scanned in linear time after sorting. That is an implementation detail. The scientific optimisation is scheduling: predicting which events will reach a short, highly sensitive peak, and allocating follow-up before it. Latency, weather and visibility matter as much as model throughput here.

Finding anomalies: AnomalyFinder

AnomalyFinder fits a point-source single-lens model, aligns the datasets, and scans the residuals with the EventFinder shapes. The crucial change is that the linear amplitude may be negative, so the same machinery finds dips. A caustic crossing may give two bumps or a U shape. The screening template only has to locate the anomaly, because the physical fit comes later (Zang et al. 2021).

The first study selected 1,216 events with catalogue $I < 19$, removed problematic ones after preliminary fits, and searched the remaining 997. It used 13 duration levels from 0.30 to 9.47 days, centres spaced by $t_{\rm eff}/6$, and windows of $\pm 3t_{\rm eff}$. A window needed at least five points and three consecutive residuals at least $2\sigma$ from zero. Two statistics ranked the candidates,

the improvement over zero residual and over a flat offset. The second separates a localised feature from a level shift. The initial thresholds were $\Delta\chi^2_0 > 120$, or $\Delta\chi^2_0 > 75$ together with $\Delta\chi^2_{\rm flat} > 35$. They produced 6,320 signals in 422 events, which review reduced to 24 potential planetary events and 59 other anomalous events (Zang et al. 2021). The 6,320 are overlapping template detections, not thousands of independent false positives.

The four-season sample used different settings: durations from 0.05 to 6.65 days and, in the prime fields, $\Delta\chi^2_0 > 200$, or $\Delta\chi^2_0 > 120$ with $\Delta\chi^2_{\rm flat} > 60$. Subprime fields used 100, or 60 and 30. The three-consecutive-$2\sigma$ requirement stayed (Zang et al. 2023). An occurrence analysis must use exactly the thresholds and reductions that defined the sample. Reproducing it with the preliminary settings would give a different selection function.

Independent data are decisive at the candidate stage. In the first study an apparent bump in one KMTNet event was absent from simultaneous OGLE data and was rejected. Another bump appeared at all three KMTNet sites, but image inspection traced it to a nearby variable star (Zang et al. 2021). Agreement between sites shows that the light varied, not that it was the intended source that varied.

Free-floating planets

A planet with no host, or one so distant that the host’s effect is negligible, gives a short event of its own. A search of residuals of long stellar events misses this channel. Its signature also differs from the Paczyński curve. When $\rho$ is large, the event’s duration is set by the time to cross the source, and the magnification is diluted. For a uniform source much larger than the Einstein radius, with the lens well inside the disk, $A - 1 \simeq 2/\rho^2$. Lower lens mass can mean a smaller amplitude without a proportionally shorter event (Johnson et al. 2020).

Qian and colleagues therefore used a model-independent search for positive excess rather than a point-source template. Windows ran from 0.05 to 28 days. The statistic sums squared positive standardised residuals, with threshold $250 + 0.1\max(N_{\rm window} - 20, 0)$, and also requires the summed significance of the ten strongest positive points to exceed 32, plus consecutive high points. From 11,876 initial candidates, matching against variable stars and visual and image checks left the 36 events described above. Seven of the candidates that looked like microlensing were asteroids (Qian et al. 2025), the population that Chapter 12 exists to find and that here is contamination.

The population numbers have moved a long way. In 2011 the MOA-II analysis suggested about 1.8 unbound or wide-orbit Jupiter-mass planets per main-sequence star (Sumi et al. 2011). In 2017 OGLE searched 50 million light curves, injected 8.6 million artificial events into images to calibrate efficiency, found no excess of Jupiter-timescale events, and set a 95% upper limit of 0.25 such planets per star (Mróz et al. 2017). The earlier large population was not confirmed at that mass.

The shortest events point instead to small masses. OGLE-2016-BLG-1928 has $t_E = 41.5$ minutes and $\theta_E = 0.842 \pm 0.064$ μas. Its lens could be between a Mars and an Earth mass, depending on distance, and stellar companions are excluded only out to a projected 8 AU (Mróz et al. 2020). MOA’s nine-year survey kept 3,535 high-quality single-lens events from 6,111 candidates, six of them with $t_E < 0.5$ day, and calibrated efficiency with 64 million image-level injections that included finite-source effects. The shortest, MOA-9y-5919, has $t_E = 0.057 \pm 0.016$ day and $\theta_E = 0.90 \pm 0.14$ μas (Koshimoto et al. 2023). The accompanying mass-function analysis infers $21^{+23}_{-13}$ free-floating or very wide-orbit planets per star between 0.33 and 6,660 Earth masses (Sumi et al. 2023). That is a model-dependent integral, dominated by extrapolation to low masses, not a count.

Finite sources make the efficiency of a free-floating-planet search at least two-dimensional, $\epsilon(t_E, \theta_E)$. A one-dimensional efficiency in $t_E$ alone is an average over an assumed population,

so an efficiency calibrated under one mass function cannot be carried into a fit of a different one. MOA showed explicitly that the difference is large at the shortest durations (Koshimoto et al. 2023).

Modelling: rays, contours and global searches

Once an anomaly is found, the model has to be fitted, and every likelihood evaluation needs the finite-source magnification at every epoch. There are two families of solver.

Inverse ray shooting samples the image plane, maps each ray through the lens equation (which is explicit in that direction), and counts how many land on the source. For $R$ rays and $L$ lens masses the mapping costs $O(RL)$. A static magnification map over $P$ pixels costs $O(P)$ storage and can be reused for many trajectories. Convolving it with a source profile costs $O(P\log P)$ by FFT, after which each trajectory costs $O(N)$ interpolations. These costs exclude map construction, padding, resolution convergence and interpolation error. Maps suit global $(s, q)$ searches, but orbital motion of the lens changes the map, tiny caustics demand high resolution, and a uniform map wastes work far from where the images are. Image-centred ray shooting shoots rays only around the images of the source, which suits moving lenses, limb darkening and more than two masses (Miyazaki & Kawahara 2025).

Contour integration uses the fact that lensing conserves surface brightness. For a uniform source the magnification is the total image area divided by the source area, and Green’s theorem turns area into a boundary integral,

where $p_k$ is the parity of image $k$ and orientations must be consistent. The algorithm samples points on the source limb, solves the fifth-degree polynomial at each, connects corresponding roots into image boundaries, and integrates. Valerio Bozza’s method adds a parabolic correction and a local error estimate that place extra samples where they are needed, special rules for image creation and destruction, and handling of a cusp hidden between two limb samples. Limb darkening is integrated over concentric annuli (Bozza 2010). For $B$ limb samples and $I$ root-solver iterations the cost per epoch is about $O(BI)$, times the number of annuli. Near caustics $B$ varies wildly. Most epochs are far from caustics, and cheap tests decide when the point-source value suffices. VBBinaryLensing added quadrupole, ghost-image and planetary-caustic tests for exactly this switch (Bozza et al. 2018). Root solving itself is delicate: polynomial roots can be catastrophically inaccurate at small $q$, and the field has specialised solvers for it (Skowron & Gould 2012).

Global search is where the compute goes. The standard strategy grids $(\log s, \log q, \alpha)$, optimises $(t_0, u_0, t_E, \rho)$ and the linear fluxes at each grid point, collects every local minimum, then releases the grid parameters and refines. With $G$ geometric cells, $A$ angle seeds, $K$ optimiser evaluations, $N$ epochs and magnification cost $C_A$ per epoch, the cost is $O(GAKNC_A)$. A better sampler does not help if it starts in only one basin. The search must cover stellar-binary $q$ as well as planetary $q$, and must fit a binary-source model explicitly (Shang et al. 2026).

Yuxin Shang and colleagues supplied a concrete workload. For simulated KMTNet-like events they searched 533 geometric grid points and eight angle seeds, with 8,640 epochs per light curve. Over 198,720 generated light curves, of which 58,081 had a detectable anomaly, they found 2,705,397 local minima and refined them. Precomputed maps cut the grid-search cost from about 400 to about 4 CPU-hours per light curve. The total reported cost was $3.34 \times 10^5$ CPU-hours (Shang et al. 2026).

GPUs, gradients and neural posteriors

Three lines of work since 2021 attack the solver cost directly.

Twinkle, by Suwei Wang, Lile Wang and Subo Dong, ports contour integration to GPUs (CUDA for NVIDIA, HIP for AMD). It reformulates the polynomial coefficients to avoid catastrophic cancellation, makes ghost-image rejection more robust, and schedules work so that expensive unfinished cases are compacted and redispatched rather than stalling completed ones. On an RTX 4090 against one core of an AMD EPYC 7H12, it reports about a hundredfold speedup in its tests, and recommends batches of more than about 200 source positions to fill the device (Wang et al. 2025). That is GPU against one CPU core, not a hundredfold algorithmic gain on matched hardware, and not a speedup of complete inference.

microlux, by Haibin Ren and Wei Zhu, implements adaptive contour integration in JAX, with an error criterion for the convergence of gradients as well as of magnification. Its magnification is about 1.5 to 5 times slower than VBBinaryLensing, but gradients make Hamiltonian Monte Carlo possible. For KMT-2019-BLG-0371 the No-U-Turn sampler needed about 33 likelihood calls per effective sample against 158 for an ensemble sampler, a roughly fivefold wall-time advantage. The tests ran on an M2 MacBook Pro, and the gain is local: it assumes the mode is already known (Ren & Zhu 2025).

microJAX, by Shota Miyazaki and Hajime Kawahara, implements differentiable image-centred ray shooting for binary and triple lenses. On a 40-GB A100 against single-threaded VBBinaryLensing on an EPYC 7452, it is 5 to 6 times faster for small limb-darkened sources but 4 to 5 times slower for large uniform sources. Memory scales as $O(N_r N_\theta N_{\rm batch})$, for radial and angular grid sizes and batch size, so epochs must be batched. An example injection-and-recovery run with 500 warm-up steps and 5,000 retained samples took about 19 hours (Miyazaki & Kawahara 2025). Fast forward evaluation and cheap complete inference are different things.

Automatic differentiation through a root solver also needs care. For a simple root, implicit differentiation gives $dz/d\lambda = -(\partial p/\partial\lambda)/(\partial p/\partial z)$, which is ill-conditioned where roots merge, at exactly the caustics where the information is. microJAX uses custom derivatives, including smooth sigmoid approximations to step functions, so its gradients are derivatives of a numerical model and should be validated as such.

Neural posterior estimation takes a different route. Keming Zhang, Joshua Bloom, Scott Gaudi and colleagues trained a ResNet–GRU encoder and a 20-block masked autoregressive flow on simulated Roman binary-lens light curves. From one million simulations they kept 291,012 distinguishable binaries, split 276,461 for training and 14,551 for testing. Training took about 25 hours on four 2080 Ti GPUs, and sampling runs at about $10^5$ draws per second on one GPU (Zhang et al. 2021). These are draws from a learned approximation, not exact-likelihood samples. The principled use is as a proposal: suggest the plausible modes quickly, then evaluate the true likelihood and refine each one. The simulation prior, the cadence, the noise model and the higher-order effects included in training define where the network can be trusted. It cannot certify that a mode it never learned is absent.

Completeness, occurrence and Roman

To turn detections into rates, inject planets into each observed event over $(s, q, \alpha)$ and $\rho$, using that event’s actual sampling and noise, run the same detection and classification pipeline, and record whether the planet is recovered and uniquely interpreted. The sensitivity of event $j$ and of the survey are

For an occurrence density $f(s, q)$, a schematic population likelihood is

where $\mathcal L_k$ is the likelihood of planet $k$‘s parameters conditional on its having been retained; using an unconditional likelihood instead would count selection twice. Real analyses also have to handle normalisation and multi-planet events (Zang et al. 2025). The general machinery of selection functions is the subject of Chapter 16 and Chapter 35.

The KMTNet occurrence analysis used 63 planets in 60 events, down to $q = 6.79 \times 10^{-6}$ for OGLE-2016-BLG-0007. It infers $0.65^{+0.35}_{-0.15}$ planets per star per dex in $s$ over $-5.2 < \log q < -1.5$, of which $0.35^{+0.34}_{-0.14}$ fall in the super-Earth interval $-5.2 < \log q < -4.5$. A double-Gaussian mass-ratio distribution improved the population likelihood by $-2\Delta\ln\mathcal L = 22.6$, with a null-model probability of $1.60 \times 10^{-4}$ (Zang et al. 2025). That probability tests the shape of the population model, not the reality of any planet. The headline, that super-Earths are common in Jupiter-like orbits, is the kind of result only microlensing can currently provide.

Degeneracy enters the sensitivity directly. In Shang and colleagues’ study, a simulated planet counted as detected only if no competitive stellar-binary solution existed and competing planetary solutions did not have very different $q$. Of 58,081 anomalous simulations, 1,052 failed. After weighting by sensitivity, degeneracies reduce overall planet sensitivity by about 5 to 10%. The raw ratio and the weighted effect are different numbers. The authors project that a full realistic-sample calculation would need $O(10^8)$ CPU-hours, and suggest $O(10^5)$ GPU-hours as a hardware-based projection, not a benchmark (Shang et al. 2026).

The Nancy Grace Roman Space Telescope will run a Galactic Bulge Time-Domain Survey designed for this science. Pre-launch forecasts and the adopted design must be kept apart. Matthew Penny and colleagues’ forecast assumed six 72-day seasons, about 1.97 deg² and a 15-minute cadence, and predicted 1,428 bound planets between 0.1 and $10^4\,M_\oplus$ with a detection criterion of $\Delta\chi^2 > 160$ (Penny et al. 2019). Samson Johnson and colleagues’ isolated-planet forecast used $\Delta\chi^2 \ge 300$ plus six consecutive $3\sigma$ points and predicted 249 free-floating planets under its fiducial mass function (Johnson et al. 2020). Under the 2023 MOA mass function the isolated-planet prediction becomes $988^{+1848}_{-566}$, including $575^{+1733}_{-424}$ between 0.1 and 1 Earth mass (Sumi et al. 2023). That jump reflects uncertainty in the population, not a change in Roman’s hardware.

The implementation documentation updated on 24 August 2026 specifies six fields totalling 1.7 deg², a 12.1-minute cycle through them, six high-cadence seasons and four low-cadence seasons with visits about five days apart. High-cadence monitoring is generally 70.5 days per season, with an exception for the first (STScI 2026). The pre-launch strategy lists 66.9-s exposures (RGES-PIT 2026a). These differ from the forecasts’ assumptions, and the old totals have not been recomputed for this design. Roman’s planned efficiency pipeline injects planets into images or light curves, propagates them through photometry and modelling, and publishes the injected and recovered catalogues (STScI 2025). The Roman Microlensing Data Challenge 2026 supplies simulated light curves for classification and fitting, requires documented code and CPU-hours, and scores computational efficiency. Its experienced-tier deadline was 2 October 2026, with evaluation in October and final reports due in November, so no results were available when this chapter was written (RGES-PIT 2026b).

Landmark results

Software

A reproduction should pin the lens engine and wrapper versions; state whether Einstein units refer to total or host mass, the source profile, the time standard and the observatory ephemerides; fit fluxes analytically; record every mode seed and acceptance cut; keep image-quality masks; and publish both injected and recovered catalogues. For JAX and GPU codes, report cold-start compilation and warm throughput separately. A numerical tolerance on magnification is not the same thing as a fractional photometric error. Online photometry from survey event pages is optimised for rapid triage and is usually not the reduction used to claim a planet (Poleski & Yee 2019).

Limitations and open problems

Missing data cannot be recomputed. A short perturbation may fall between exposures, in daylight, behind clouds, or outside a space mission’s season. More computation can recover a marginally sampled signal. It cannot reconstruct an unobserved one without adding prior assumptions. The finite source sets a further physical limit, a source-dependent loss of contrast.

Degeneracy is more than close versus wide. Central caustics often give similar light curves at $s$ and $1/s$ with similar $q$, which is harmless for demographics. But inner versus outer trajectories, close versus wide planetary-caustic events, central versus resonant solutions, and accidental degeneracies created by gaps in the data can give different $q$ or $\rho$. A planetary model can resemble a stellar binary, and a binary source can mimic a planet. Colour information helps against binary sources, but sparse colour sampling and similar source colours can defeat it (Shang et al. 2026).

Higher-order effects interact. Annual parallax, lens orbital motion and source orbital motion (xallarap) distort the trajectory in partly similar ways. Satellite parallax has discrete sign alternatives. These are problems of identifiability, not of optimisation, and only new information breaks them: observatories with useful separation, late-time imaging that resolves lens and source, or astrometry.

An isolated lens is not necessarily unbound. The absence of a detectable host excludes a range of separations, orientations and host masses, not every host. A short $t_E$ also admits unusually fast or nearby lenses. Measuring $\theta_E$ is strong evidence for low mass but leaves the distance uncertain, and the large free-floating population integrals rest on extrapolation and Galactic priors (Mróz et al. 2020; Sumi et al. 2023).

Systematics dominate the tails. Nearby variables, reference-image artefacts, seeing correlations, saturation and moving objects can all mimic bumps and dips. The KMTNet full-frame result showed that improved subtraction and better source positions can matter more than a new classifier, and a catalogue-only search misses sources absent from the reference list (Qian et al. 2025). Only injection into images tests the whole chain.

Numerical and statistical convergence are different. Roots can be inaccurate at small $q$. Root continuation can connect images wrongly across a caustic. Magnification convergence does not ensure gradient convergence. And a perfectly evaluated likelihood, sampled by a long, well-mixed chain in one minimum, can still miss another minimum that matters. Global mode recovery, effective sample size within a mode, and relative weights of modes each need their own validation (Skowron & Gould 2012; Ren & Zhu 2025).

The field’s own statement of its unsolved problem is complete, scalable sensitivity calculations that include degenerate models and realistic gaps and systematics. Shang and colleagues state that their corrections from representative events are not sufficient for direct application to every real event, and that their grid still misses closely spaced resonant minima (Shang et al. 2026). Put more sharply: no method today gives a calibrated upper bound on the probability of having missed a viable alternative interpretation. That probability matters for single discoveries and for occurrence rates alike.

Where the room is

It helps to sort the limits first. Alignment probability, photon statistics, finite-source averaging, missing time coverage and genuine degeneracy in the observables are physical: no algorithm removes them. Inaccurate photometry, slow or incomplete search, late alerts, numerical failures and impractically expensive completeness experiments are engineering. A good new method attacks the second list without disguising the first as smaller error bars. The broader argument for where order-of-magnitude gains are plausible is in Chapter 42. For microlensing the opportunities are specific.

Make the hierarchy preserve rare modes. A sensible cascade is image-quality filtering, then a signed anomaly scan, a geometric proposal bank, a map-based coarse likelihood, an accurate solver, multi-mode refinement, and finally physical interpretation. Neural posterior proposals can replace much of a blind grid, but an independent coarse bank or explicit exploratory proposals should be kept to catch what the network omits. Caustic morphology suggests separate parameterisations for central, planetary and resonant events, because one uniform logarithmic grid is inefficient in all three. Early rejection is possible: partial sums of nonnegative squared residuals bound the final $\chi^2$ when the model and fluxes are fixed. But when fluxes are refitted as data are added, the running residual is not automatically a lower bound on the final profiled score, and a safe implementation needs bounds that account for that. A rigorous branch-and-bound search near caustic topology changes would be new. It is a research proposal, in the spirit of the safe pruning discussed in Chapter 38.

Batch irregular work explicitly. Twinkle already compacts unfinished tasks. The idea extends naturally: separate point-source epochs, mild finite-source corrections and difficult caustic crossings into queues bucketed by estimated contour complexity, and batch many parameter vectors and many events together rather than the epochs of one light curve, since a single event rarely has enough work to fill a GPU. Reuse polynomial coefficients, nearby root seeds and ephemerides where accuracy tests allow (Chapter 40). Mixed precision suits map interpolation and regular reductions, but small-$q$ coefficients and nearly coincident roots are numerically sensitive, and a final FP64 solve or independently verified residual may be needed. Tensor cores fit neural encoders and large dense contractions. An adaptive complex root solver is not a tensor-core problem. An order of magnitude over a scalar CPU baseline is already observed. Another order of magnitude over a tuned GPU code, at matched error and with complete mode recovery, is unproven.

Optimise the population experiment, not the single fit. Completeness is the larger repeated workload: one observed event may need many injected planets, each with its own search over alternative models. Importance sampling toward the boundaries where efficiency changes, stratification by magnification, source size and cadence, and sharing of calculations between injections with the same geometry could all reduce the cost. The target is the uncertainty of the final occurrence distribution per unit compute, with unbiased weights. An injection scheme that skips rare ambiguous configurations gives precise-looking but biased answers. For free-floating planets, keep the two-dimensional $\epsilon(t_E, \theta_E)$ table and integrate it jointly with the population model, rather than regenerating one-dimensional corrections. For bound planets, the measured 5 to 10% degeneracy effect is a scale for the scientific error: compute savings are worth most if they make realistic, degeneracy-aware calibration affordable for every event, instead of a universal correction from idealised ones.

Transforms help only where the sampling allows. At fixed duration, an anomaly-search numerator is a correlation of residuals with a translated template, and on a regular time grid FFTs or multiresolution filter banks could share work across centres (Chapter 39). The denominator needs its own correlation with weights and masks. Irregular times, weather gaps, clipping and per-site calibration erode the advantage, and naive interpolation can invent correlations or dilute a short dip. The window-limited direct scan is already efficient, so a transform bank must beat it end to end. And there is no folding shortcut of the kind a transit search uses (Chapter 30): a microlensing event happens once.

Build an equal-sensitivity benchmark. Claims of speed are cheap when sensitivity is not held fixed. A credible comparison uses the same raw images or frozen light curves, the same timestamps, masks, source population, injected lens parameters and non-planetary contaminants. It compares against both a pinned historical code and a strong current one. It measures weak-anomaly completeness at a fixed candidate-review budget or false-positive rate, recovery of alternative modes, bias and coverage in $q$, $\rho$ and $\pi_E$, and end-to-end wall time, energy, peak memory, compilation and training cost. It includes short free-floating events, negative dips, incomplete caustic crossings, stellar-binary alternatives and lens orbital motion. The Roman data challenge is a starting point for the inference half. Image-level selection needs tests of its own.

The summary is that microlensing’s hardest computational problem is not making one likelihood fast. It is proving, at survey scale and for every injected planet, that nothing else fits. A method that delivers that proof cheaply would change both the discovery papers and the occurrence rates that rest on them.

References

1. Mao, S. & Paczyński, B. (1991). Gravitational microlensing by double stars and planetary systems. ApJ Letters 374, L37. doi:10.1086/186066
2. Gould, A. & Loeb, A. (1992). Discovering planetary systems through gravitational microlenses. ApJ 396, 104. doi:10.1086/171700
3. Griest, K. & Safizadeh, N. (1998). The use of high-magnification microlensing events in discovering extrasolar planets. ApJ 500, 37. arXiv:astro-ph/9710342
4. Bond, I. A., Udalski, A. et al. (2004). OGLE 2003-BLG-235/MOA 2003-BLG-53: a planetary microlensing event. ApJ Letters. arXiv:astro-ph/0404309
5. Chung, S.-J., Han, C., Park, B.-G. et al. (2005). Properties of central caustics in planetary microlensing. ApJ. arXiv:astro-ph/0505363
6. Beaulieu, J.-P., Bennett, D. P. et al. (2006). Discovery of a cool planet of 5.5 Earth masses through gravitational microlensing. Nature. arXiv:astro-ph/0601563
7. Bozza, V. (2010). Microlensing with an advanced contour integration algorithm: Green's theorem to third order, error control, optimal sampling and limb darkening. MNRAS 408, 2188. arXiv:1004.2796
8. Sumi, T. et al. (2011). Unbound or distant planetary mass population detected by gravitational microlensing. Nature. arXiv:1105.3544
9. Skowron, J. & Gould, A. (2012). General complex polynomial root solver and its further optimization for binary microlenses. arXiv preprint. arXiv:1203.1034
10. Udalski, A., Szymański, M. K. & Szymański, G. (2015). OGLE-IV: fourth phase of the Optical Gravitational Lensing Experiment. Acta Astronomica 65, 1. arXiv:1504.05966
11. Mróz, P. et al. (2017). No large population of unbound or wide-orbit Jupiter-mass planets. Nature. arXiv:1707.07634
12. Kim, D.-J., Kim, H.-W., Hwang, K.-H. et al. (2018). Korea Microlensing Telescope Network microlensing events from 2015: event-finding algorithm, vetting, and photometry. AJ 155, 76. arXiv:1703.06883
13. Kim, H.-W., Hwang, K.-H. et al. (2018). The Korea Microlensing Telescope Network (KMTNet) alert algorithm and alert system. arXiv preprint. arXiv:1806.07545
14. Bozza, V., Bachelet, E., Bartolić, F. et al. (2018). VBBinaryLensing: a public package for microlensing light-curve computation. MNRAS 479, 5157. arXiv:1805.05653
15. Poleski, R. & Yee, J. C. (2019). Modeling microlensing events with MulensModel. Astronomy and Computing 26, 35. arXiv:1803.01003
16. Penny, M. T. et al. (2019). Predictions of the WFIRST microlensing survey. I. Bound planet detection rates. ApJS 241, 3. arXiv:1808.02490
17. Johnson, S. A., Penny, M., Gaudi, B. S. et al. (2020). Predictions of the Nancy Grace Roman Space Telescope Galactic Exoplanet Survey. II. Free-floating planet detection rates. AJ. arXiv:2006.10760
18. Mróz, P. et al. (2020). A terrestrial-mass rogue planet candidate detected in the shortest-timescale microlensing event. ApJ Letters. arXiv:2009.12377
19. Zang, W. et al. (2021). Systematic KMTNet planetary anomaly search. I. OGLE-2019-BLG-1053Lb, a buried terrestrial planet. AJ 162, 163. arXiv:2103.11880
20. Zhang, K., Bloom, J. S., Gaudi, B. S., Lanusse, F., Lam, C. & Lu, J. R. (2021). Real-time likelihood-free inference of Roman binary microlensing events with amortized neural posterior estimation. AJ 161, 262. arXiv:2102.05673
21. Hwang, K.-H. et al. (2022). Systematic KMTNet planetary anomaly search, Paper II: six new q < 2 × 10⁻⁴ mass-ratio planets. AJ 163, 43. arXiv:2106.06686
22. Zang, W. et al. (2023). Systematic KMTNet planetary anomaly search, Paper VII: complete sample of q < 10⁻⁴ planets from the first four-year survey. AJ. arXiv:2210.12344
23. Koshimoto, N., Sumi, T., Bennett, D. P. et al. (2023). Terrestrial- and Neptune-mass free-floating planet candidates from the MOA-II 9 yr Galactic Bulge survey. AJ. arXiv:2303.08279
24. Sumi, T., Koshimoto, N. et al. (2023). Free-floating planet mass function from MOA-II 9 yr survey toward the Galactic Bulge. AJ. arXiv:2303.08280
25. Bozza, V., Saggese, V., Covone, G., Rota, P. & Zhang, J. (2025). VBMicrolensing: three algorithms for multiple lensing with contour integration. A&A 694, A219. arXiv:2410.13660
26. Wang, S., Wang, L. & Dong, S. (2025). Twinkle: a GPU-based binary-lens microlensing code with the contour integration method. ApJS 276, 40. arXiv:2501.03322
27. Ren, H. & Zhu, W. (2025). A differentiable binary microlensing model using adaptive contour integration method. AJ. arXiv:2501.07268
28. Qian, Q., Yang, H. et al. (2025). Systematic search for FFPs in KMTNet full-frame images. I. Photometry pipeline. PASP 137, 064401. arXiv:2503.24097
29. Gui, Y., Zang, W. et al. (2025). Systematic KMTNet planetary anomaly search. XII. Complete sample of 2017 subprime field planets. arXiv preprint. arXiv:2504.20155
30. Zang, W., Jung, Y. K., Yee, J. C. et al. (2025). Microlensing events indicate that super-Earth exoplanets are common in Jupiter-like orbits. Science 388, 400. arXiv:2504.20158
31. Miyazaki, S. & Kawahara, H. (2025). microJAX: a differentiable framework for microlensing modeling with GPU-accelerated image-centered ray shooting. ApJ 994, 144. arXiv:2510.02639
32. Shang, Y., Yang, H., Zhang, J., Mao, S., Gould, A., Zang, W., Qian, Q. & Yee, J. C. (2026). Binary-lens microlensing degeneracy: impact on planetary sensitivity and mass-ratio function. arXiv preprint, v3. arXiv:2509.16191
33. Shin, I.-G., Yee, J. C., Zang, W. et al. (2026). Systematic KMTNet planetary anomaly search. XIII. Complete sample of 2021 prime field planets. arXiv preprint. arXiv:2605.17053
34. STScI (2025). MSOS detection efficiency pipeline. Roman User Documentation. roman-docs.stsci.edu
35. STScI (2026). Galactic Bulge Time-Domain Survey (updated 24 August 2026). Roman User Documentation. roman-docs.stsci.edu
36. RGES-PIT (2026a). RGES design. Roman Galactic Exoplanet Survey Project Infrastructure Team. rges-pit.org/design
37. RGES-PIT (2026b). Roman Microlensing Data Challenge 2026. Roman Galactic Exoplanet Survey Project Infrastructure Team. rges-pit.org/data-challenge