Confirming candidate self-lensing binary detections with radial velocity measurements

Introduction

There are expected to exist vast numbers of compact objects in binary systems in the Milky Way (e.g. Shao & Li 2019; Wiktorowicz et al. 2020; Weller & Johnson 2023). The majority of these are not close enough to interact with their companions, making them challenging to detect. This is especially true of black holes, which would produce negligible electromagnetic emissions in such systems.

The detection of gravitational microlensing is a proven technique (Paczynski 1986; Paczyński 1996; Mao 2012), in which the gravity of a foreground object (the lens, which is often not visible) amplifies flux from a background object (the source). This has been used to detect and study various classes of object, including exoplanets (Bond et al. 2004; Gaudi 2012) and black holes (Sahu et al. 2022; Lam et al. 2022). Self-lensing is a special case of microlensing where the source and the lens are in a binary system together, in an orbit which is close to edge-on, producing a periodic microlensing amplification of the source once per orbit (Maeder 1973; Agol 2002). Self-lensing is only possible when the lens is a compact object. This means a detection of self-lensing (i.e. periodic microlensing) is de facto a detection of a compact object.

A self-lensing search is an intriguing way to search for compact objects in binaries which would be undetectable by other means. To date, five self-lensing binaries have been confirmed, each containing a white dwarf lens, all of which have relatively close, short-period orbits (Kruse & Agol 2014; Kawahara et al. 2018; Masuda et al. 2019). However, it is expected that stellar-compact object binaries should span the full range of compact object masses, as well as the full range of orbital configurations.

Photometry remains the best way to detect self-lensing binaries, and previous works have assessed the detectability of self-lensing binaries in various photometric surveys (Wiktorowicz et al. 2021). These results look promising, with modern all-sky surveys expected to imminently uncover a population of compact objects in self-lensing binaries, including the currently under-studied population of quiescent black holes.

While detection of self-lensing systems with photometry alone has been demonstrated, modelling of such systems and characterisation of the compact component can be subject to several degeneracies between the physical parameters which can make it challenging to accurately determine the mass of the lens in cases where the amplification due to self-lensing, the crossing time, or the stellar companion properties are not well known (see Figure 1 for an illustration of these relationships). For short-period binaries it should be possible to derive a photometric radial velocity (RV) by modelling Doppler boosting [37] (Loeb & Gaudi 2003; Zucker et al. 2007; Faigler & Mazeh 2011), but for systems with no detectable Doppler boosting effect, spectroscopic observations are the best way to solve for the binary’s physical parameters (Kawahara et al. 2018; Yahalomi et al. 2019). This is especially true for systems with the longest orbital periods, for which we may only be able to observe a single lensing event in a human lifetime (i.e. those with periods lasting decades, centuries, or longer).

In this work, we consider the question of how many spectroscopic RV measurements are required, and to what level of precision, in order to accurately retrieve the self-lensing system parameters as a function of those parameters. This includes systems with orbital periods on the order of centuries for which the orbital period can not be determined photometrically. This will inform observing strategies for follow-up of the large number of candidates expected to come from new surveys, allowing the most efficient use of resources. We also consider the usefulness of photometric-only fitting to guide RV observing strategy, developing a method of computing an RV observing schedule to maximise the separation between plausible solutions.

We will show that it is possible to accurately model very long period self-lensing binaries with a handful of RV measurements, even when periodicity can not be observed. In many cases, photometric-only fits can come close to determining the properties, removing the need for many RV measurements to retrieve the parameters accurately.

We begin by describing our methods for generating and fitting our simulated data in § 2, present and discuss the fitting results in § 3, and conclude in § 4.

Methods

We investigate the relationship between the number and precision of RV measurements and the relative error in recovered physical parameters as a function of lens mass, compact object type, and orbital period. We evaluate parameter accuracy using the median relative error, and separately consider recovery of broad compact-object class (i.e. distinguishing black holes from other compact objects even if the fitted mass remains imprecise).

When considering the recovery of compact object class, we adopt the following mass bins and object classes: $M_{\mathrm{lens}} < 3\,M_{\odot}$ (white dwarfs or neutron stars); $3 \le M_{\mathrm{lens}}/M_{\odot} < 5$ (lower mass gap); $5 \le M_{\mathrm{lens}}/M_{\odot} < 20$ (lower mass black holes); $20 \le M_{\mathrm{lens}}/M_{\odot} < 50$ (higher mass black holes). We separate the black hole mass ranges to assess how lens mass affects class recovery. A class is counted as recovered when the fitted lens mass falls in the same bin as the true lens mass.

Photometric detection of self-lensing signals is not the focus of the present work. Here we assume all of the simulated photometric lensing events have been detected and their precise timing has been determined. A robust method for carrying out such detections will be presented in a future work, including evaluation against self-lensing signals injected into real light curves. For now we focus only on the recovery of the physical properties of self-lensing systems which have already been detected, and how the RV observing strategy affects this.

As the basis of this work, we have generated a library of self-lensing photometric simulations covering a range of physical parameters, filtered to include only the simulated light curves which contain visible lensing, where the peak amplitude is larger than the photometric noise. We also exclude systems overflowing their Roche lobe and any systems where lensing is not geometrically possible. These simulations are described in § 2.1.

An MCMC fit is then performed against the simulated photometry, to sample the posterior distribution of the self-lensing physical parameters and find distinct physical solutions for each light curve. We assume the stellar properties, mass, radius, effective temperature, and distance from Earth, are known. The set of distinct solutions is then used to compute an RV observing schedule using a variable number of observations (between three and ten observations) over a fixed five-year campaign, such that each schedule is optimised to distinguish between the photometric solutions.

Radial velocity measurements are then simulated for the computed observing schedule, using the same set of simulation parameters as the photometry (§ 2.4). Gaussian noise is added to simulate a range of spectrograph accuracies.

A second MCMC fit is then performed to sample the posterior distribution constrained jointly by the photometry and RV (§ 2.6). We then aggregate the median relative errors across observing parameters and self-lensing physical parameters to investigate how each observing schedule affects parameter retrieval.

Finally, we map the fitted outcomes onto previous population synthesis results to translate the estimated recovery error rates onto the expected real population of detectable self-lensing systems (§ 2.7).

Self-lensing and photometric binary effects

The photometric model used here is a modified version of COBALT-BLUE (Dixon et al. (2026)), created while development of COBALT-BLUE was still in progress. Since the main focus of the present work is radial velocity scheduling and fitting, not the development of a full photometric forward model, we adopt simplified approximations to some aspects of the full model. Compared with COBALT-BLUE, our model adopts a local approximation to the eccentric self-lensing geometry and lower-order treatments of Doppler boosting and ellipsoidal variability, and we only consider linear limb darkening (Dixon et al. (2026) shows that linear limb darkening is sufficient for the present purposes). This work also adds a stellar noise model (§ 2.2) and the spectroscopic radial-velocity component (§ 2.4) which are not included in COBALT-BLUE. A full mathematical description of the modified photometric model is given in Appendix § 7.

Simulations comprise a sample of self-lensing system parameters which is physically complete but not physically realistic (i.e. parameters are not sampled according to the frequency they are expected to occur in nature). This allows us to determine the ranges of parameters which produce self-lensing effects, independent of how common or uncommon particular parameter ranges are expected to be.

The photometric simulations comprise realistic simulated light curves with the same length and observing cadence as TESS light curves for each self-lensing system. When generating the simulations, the parameter values are sampled uniformly between the limits given in Table 1, with the exception of semi-major axis ($a$), eccentricity ($e$), and lens mass ($M_{\mathrm{lens}}$), which are sampled on a log scale to increase the number of samples at smaller values (the light curve morphology changes roughly in proportion to the relative change in these values; e.g. a 5 AU change in semi-major axis has a much larger effect going from 1 AU to 6 AU than for 100 AU to 105 AU).

Parameter	Description	Range	Spacing
$a$	Semi-major axis	0.1–200 AU	Logarithmic
$i$	Orbital inclination	89.99–90.00°	Linear
$M_{\rm lens}$	Lens mass	0.5–50 $M_\odot$	Logarithmic
$M_\star$	Stellar mass	0.5–50 $M_\odot$	Linear
$e$	Eccentricity	0–0.999999	Logarithmic
$\omega$	Periapsis phase angle	0–360°	Linear
$d$	Distance from Earth	0.1–100 kpc	Linear

Three photometric effects are included when generating the light curves: self-lensing, ellipsoidal variations, and Doppler boosting. At small values of $a$ the light curve is dominated by ellipsoidal variations and Doppler boosting (though lensing is often visible too); at large $a$ lensing is dominant, with ellipsoidal variation and Doppler boosting tending to zero as $a$ increases ((Sorabella et al. 2022)).

The peak amplitude and effective crossing time loosely correlate with $M_{\rm lens}$, $P_{\rm orb}$, and $i$ (Figure 1). Since the amplitude and self-lensing event duration are directly measurable from the light curve, this can in principle allow the physical properties of the lens to be estimated. However, there is scatter in the relationship between these properties, since multiple lens configurations can produce similar amplifications and crossing times. Care must therefore be taken to account for this degeneracy when using these relationships (see discussion of MCMC chain initialisation in § 2.3).

These effects all depend on the orbital phase at a given time during the simulation. The zero point for this is defined as the moment of maximum lensing magnification (i.e. the time of minimum projected separation between source and lens, $t_0$). The periapsis phase angle, $\omega$, is defined in terms of this zero point; i.e. the phase angle at which periastron passage occurs, relative to the position at $t_0$. Note that this is different from the convention used in COBALT-BLUE, which measures orbital phase relative to periastron passage.

All simulated light curves are 27 days long (the duration of a TESS sector). The peak time, $t_0$, is randomly selected to fall somewhere in the 27 day time window. An example light curve for a simulated self-lensing event can be seen in the top-left panel of Figure 2.

Stellar noise

Gaussian processes (GPs) are designed to model smoothly varying correlated noise ((Rasmussen & Williams 2005)). In the exoplanet community, this has been highly successful in allowing exoplanet signals to be disentangled from stellar noise, though it does create significant computational expense which makes it impractical to scale GP-based detection to the millions of light curves being produced by current surveys.

In this work, we utilise a GP to generate realistic light curves with stellar noise, sampled according to the TESS observing cadence, instrumental sensitivity, and bandpass. We combine several kernels to model granulation and asteroseismic modes, based on the simulated stellar properties. We use tin ygp (Foreman-Mackey et al. 2026) with quasi-separable simple harmonic oscillator (SHO) kernels.

Granulation is a form of aperiodic correlated noise due to convection cells at the stellar surface. Granulation appears as low-amplitude correlated noise; aperiodic in nature, it prefers certain timescales dictated by the size and age of the star (Kjeldsen & Bedding 1995). As derived in Foreman-Mackey et al. (2021), we consider this a characteristic timescale of a damped stochastic process. We model such a process using a damped harmonic oscillator kernel with a fixed quality factor ($Q = \frac{1}{\sqrt{2}}$) parameter to create a damped random walk (DRW). This creates frequency-dependent broadband noise with a characteristic timescale but no dominant frequencies. The characteristic timescale in the model $\tau$ then reflects the average lifetime of a single convective cell, while the amplitude scales with the number of granules contributing to the observable flux and the source-star surface gravity, $g_\star$ (Kjeldsen & Bedding 1995; Pereira et al. 2019; Foreman-Mackey et al. 2021). The granulation timescale, $\tau_{\rm gran}$, for a source star of mass $M_\star$ and radius $R_\star$, goes as

This model captures the temporal decorrelation of the granulation component light curve over time, appropriate for granulation, where surface features evolve on a finite timescale and exhibit only short-range temporal correlations (Aigrain & Foreman-Mackey 2023). This model has also been tested and applied in previous GP transit studies ((Brewer & Stello 2009; Nicholson & Aigrain 2022)). Overall, the DRW kernel effectively models granulation noise without requiring high model complexity or prior knowledge of the star’s convective dynamics.

Asteroseismic oscillations are caused by the propagation of sound waves through the star’s outer convective zone (Kjeldsen & Bedding 1995; Chaplin & Miglio 2013). These occur on the timescales of minutes to tens of hours, placing them directly into the observing window of light curves from instruments such as TESS. Again, we choose a damped harmonic oscillator kernel, but this time with a free quality factor to allow for the presence of a periodic component (Brewer & Stello 2009; Foreman-Mackey et al. 2021). The frequency of the primary mode, given by the $\nu_{\rm seis}$ hyper-parameter, is proportional to the source-star surface gravity, $g_\star$, and the effective temperature, $T_{\rm eff}$, and therefore depends on the stellar properties of the simulated system. We draw the oscillation quality factor independently for each simulated star from a broad log-uniform prior, $Q \in[100, 10^4]$, chosen to represent coherent p modes.

where $M_\star$, $R_\star$, and $T_{\rm eff}$ are the mass, radius, and effective temperature of the source star.

In the GP implementation, these frequency scalings are evaluated from internally consistent stellar relations rather than fitted independently for every light curve. For a simulated stellar mass $M_\star$, we estimate the stellar radius and effective temperature from main-sequence scaling relations. Let $\mathcal{R}$ be the dimensionless mass-radius relation, $L(M_\star)$ the corresponding mass-luminosity relation, and $\sigma_{\rm SB}$ the Stefan–Boltzmann constant; $M_\odot$ and $R_\odot$ are the solar mass and radius. Then

The asteroseismic frequency is then initialised using the usual $\nu_{\max}$ scaling relation. Let $\nu_{\max,\odot}$ be the solar value of $\nu_{\max}$, and let $T_\odot$ be the solar effective temperature. Then

The granulation characteristic frequency is taken to be proportional to this value, adopting values around an empirical proportionality constant calibrated by Kallinger et al. 2014,

The kernel hyperparameters are stored as $\log_{10}\nu$ in Hz and converted to angular frequency internally when constructing the SHO kernels.

To calculate the source brightness, we approximate the stellar spectrum as a blackbody with the radius and effective temperature derived above. The spectral flux density at Earth is

where $B_\lambda$ is the Planck function, $\lambda$ is wavelength, and $d$ is the source distance. We integrate the corresponding photon flux over the TESS bandpass, weighted by the instrumental response and detector quantum efficiency. Dividing this integral by the same integral for an AB zero-magnitude reference spectrum ($f_\nu= 3631$ Jy) gives the ratio used to calculate the apparent TESS-band magnitude, $m_{\rm TESS}$, via $-2.5\log_{10}$ of that ratio. We then convert this magnitude to the expected detected source counts per cadence,

using the instrumental magnitude zeropoint of 20.44 for a count rate of one electron per second and the simulated cadence $\Delta t = 600$ s. Kernel amplitudes are set in counts per cadence by scaling the fractional stellar-variability amplitudes with $N_{\rm TESS}$. The standard deviation of GP component $q$ is initialised as

where $f_q$ is a fractional variability amplitude. We draw these amplitudes independently for each component and each simulated star from uniform distributions spanning 50–300 ppm for granulation and 100–1000 ppm for asteroseismic variability.

The same apparent magnitude is also used to set the TESS white-noise level. We linearly interpolate an empirical TESS noise model at $m_{\rm TESS}$ to obtain a fractional per-cadence uncertainty $\eta(m_{\rm TESS})$. For the GP covariance matrix $\Sigma_{ij}$ and the Kronecker delta $\delta_{ij}$, this gives a TESS white-noise standard deviation $\sigma_{\rm TESS}$ and diagonal covariance contribution

This term is fixed by the simulated source brightness rather than fitted as an additional free GP jitter parameter.

The variability caused by magnetic effects has also been included in prior studies ((Haywood et al. 2014; Luger et al. 2021; Foreman-Mackey et al. 2021)), as has the movement of star spots about the rotational axis of the star. However, such variations are often overpowered by more global variations, such as asteroseismic modes for distant stars, where surface structure cannot be resolved. We therefore assume that the jitter term is sufficient to account for magnetic contributions to the underlying stellar noise.

For a detailed review on modelling stellar variability through GPs, see Aigrain & Foreman-Mackey (2023).

Photometric MCMC fit

Each simulated photometric light curve is fitted via MCMC using a random-walk Metropolis-Hastings sampler implemented with the BlackJAX package (128 chains, 1,000 total steps, including 200 burn-in steps). More modern MCMC samplers, such as Hamiltonian Monte Carlo, may be more effective or more efficient, but this requires a fully differentiable model. Our photometric model is differentiable in principle, however the current implementation around edges/parameter boundaries would result in discontinuous gradients in practice, potentially negatively affecting such samplers. This could be addressed in future work.

For these fits it is assumed we know the stellar physical parameters exactly, along with the relative lensing amplification (Appendix § 7) and effective crossing time, $\tau_{\rm eff}$ (Wiktorowicz et al. 2021). In reality there will be uncertainty in these values, but determining them is a separate problem which is out of scope for the present work. The photometric fit therefore represents an idealised best-case scenario for recovery of the physical properties from the photometry, eliminating the effects of the uncertainties in stellar models and allowing us to assess the RV observing schedule in isolation.

Light curves are binned according to their crossing time to minimise computation time while preserving the shape of the self-lensing peak (Table 2). A limited number of bin widths was chosen for computational efficiency, as fits are run as vectorised batches on a GPU, requiring all the light curves in a batch to have the same length.

Effective crossing time, $\tau_{\rm eff}$	Photometric bin width
$\tau_{\rm eff} \le4$ h	10 min (native cadence)
$4$ h $\le\tau_{\rm eff} \le1$ d	1 h
$\tau_{\rm eff} > 1$ d	6 h

We fit only the fundamental physical parameters of the system,

with the stellar mass, stellar radius, effective temperature, distance, and lensing-peak time held fixed at their simulated values. Orbital period, $P_{\rm orb}$, impact parameter, maximum possible impact parameter, and effective lensing-event duration are derived from the fitted physical parameters rather than fitted directly. The orbital-period calculation is given in Appendix § 7.

Sampler chains are initialised using an estimated parameter distribution calculated from $\tau_{\rm eff}$ and $A_{\rm SL}$ (Appendix § 7). Similar to Figure 1, parameters are calculated in the region around the observed values, and then these are clustered by the parameters in $\boldsymbol{\theta}$ using $k$-means clustering to find distinct groupings of estimated parameters. The central point of each cluster is used to initialise the sampler chains. This ensures the random walk thoroughly explores the area of the parameter space with the highest likelihood given the known amplification and crossing time, without constraining it from exploring beyond those regions.

We use a uniform prior within the corresponding parameter limits in Table 1, with negative infinity probability outside these limits. We also require trial systems to satisfy the self-lensing geometry condition $b \le b_{\max}$, where $b$ is the impact parameter in units of the Einstein radius and $b_{\max}$ is its maximum value for self-lensing. Their expressions are given in Appendix § 7.

At the end of the MCMC run, we analyse the joint posterior distribution to identify distinct candidate solutions, accounting for degeneracies in which multiple self-lensing configurations can produce similar light curves. Examining the marginal posterior of each parameter separately is insufficient because the fitted parameters are correlated and only certain combinations represent physically plausible solutions. We therefore apply a heuristic clustering procedure that combines farthest-first traversal (Gonzalez 1985) with Lloyd’s $k$-means refinement algorithm (Lloyd 1982). This is applied in a scaled representation of the joint parameter space, with the periodic phase represented using its sine and cosine. The number of clusters, $k$, is determined heuristically: candidate cluster centres are added using farthest-first traversal until no sample lies beyond a specified minimum separation from the existing centres, subject to a maximum of $k = 8$. The centres are then refined using $k$-means, and clusters containing less than 3 per cent of the posterior samples are merged

into the nearest retained cluster. The upper limit of eight clusters, the minimum centre separation, and the minimum cluster size of 3 per cent were chosen empirically from inspection of a random sample of fits, balancing the representation of distinct solutions against the computational cost of constructing the RV observing schedule. Because clustering is used only to construct the observing schedule and does not enter the subsequent parameter fitting or error calculation directly, we do not interpret the resulting clusters as formal estimates of the number of posterior modes.

Radial velocity calculation

The RV simulations use the same physical parameters and phase convention as the photometric simulations (§ 2.1). We define the orbital phase relative to $t_0$ as the true orbital longitude, using the true anomaly $\nu(t)$,

where $\nu(t)$ is measured relative to periapsis. The phase satisfies $\phi(t_0) = 0$, and the true anomaly is calculated from Kepler’s equation as described in Appendix § 7.

Let $\mu_g = G(M_{\mathrm{lens}} + M_\star)$ be the standard gravitational parameter, and $h = [\mu_g a(1 - e^2)]^{1/2}$ be the specific angular momentum of the relative orbit. The source-star barycentric radial and tangential orbital-velocity components are

The simulated line-of-sight radial velocity is then

This reduces to the circular-orbit convention $v_r(t) \propto\sin\phi(t)$ when $e = 0$, but for eccentric orbits the radial component means the RV at $t_0$ is not generally zero.

We do not attempt to simulate realistic spectropic observations, which would require accounting for various instrumental factors including spectrograph stability, spectral resolution, and signal-to-noise ratio. Instead we add varying amounts of instrumental noise ($\sigma_{\mathrm{ins}} \in\{0.1, 1, 10, 100, 1000\}$ m/s) to simulate the total uncertainty in the radial velocity measurements that would result from the different instrumental effects. The values chosen for $\sigma_{\mathrm{ins}}$ span the orders of magnitude of expected RV accuracy of spectrographs in service today, from ESPRESSO (Pepe et al. 2021, $\sigma_{\mathrm{ins}} = 0.1$ m/s) to APO ARCES (Gaulme et al. 2022, $\sigma_{\mathrm{ins}} = 0.5$ km/s). We also add an additional fixed Gaussian noise to represent stochastic stellar variability (Wright 2005, $\sigma_\star= 10$ m/s). For simplicity, we do not include proper motion or heliocentric motion, as these would be corrected before fitting real observations.

An example of the radial velocity measurements for a simulated self-lensing system can be seen in the bottom-left panel of Figure 2.

Targeted radial velocity scheduling

The simulated RV curves start at randomised times between 1 year and 10 years from the end of the photometry. This is to simulate a realistic delay in identifying candidates in the photometry and organising follow-up of those candidates. For this work we assume each observing campaign will last five years, and we trial varying numbers of RV observations, $3 \leq N_{\mathrm{obs}} \leq10$. We do not attempt to simulate fully realistic observing (with observing seasons, seeing, varying weather conditions, etc.) as we are mainly interested in measuring the effects of instrument selection and number of observations.

The distinct solutions derived in § 2.3 are used to compute a schedule for the radial velocity observations, in which observations are deliberately timed to distinguish between the solutions and maximise the information gained about the binary orbit.

First, in the event that the photometric fit only found one plausible solution, the orbital period of that solution is checked to determine whether a lensing event will occur during the simulated observing window. If so, three observations are scheduled around this event, one at $t_0$, and two at $t_0 \pm\delta t$. $\delta t$ is chosen by drawing random samples from the photometric MCMC posterior and choosing the smallest value of $\delta t$ where the change in RV is larger than $5\sigma_{\mathrm{ins}}$ for $> 90\%$ of samples, assuming $\omega= \pi$; i.e. the smallest time offset where a change in RV can be reliably measured given the specified value of $\sigma_{\mathrm{ins}}$. This triplet of observations provides information about $e$ and $\omega$: if $v_r(t_0) \ne0$, then $e > 0$; if $v_r(t_0) = 0$, then either $e = 0$ or $\omega= 0$ (or both) and the value of $v_r(t_0 \pm\delta t)$ provides information about these values.

If there are multiple plausible photometric MCMC solutions, the predicted RV values for each of them over the observing campaign are calculated, and the difference in RV values for each pair of solutions is calculated (each “solution” is a cluster of posterior samples, so the separation is actually the squared difference between the mean RVs of each cluster, normalised by the error). The final set of observation times is chosen by iteratively selecting the times which maximise the minimum separation between solutions, for a minimum of three observations and a maximum of $\lfloor N_{\mathrm{obs}}/2\rfloor$.

Any remaining observations are scheduled to minimise the gaps between observations. For a circular orbit, this is equivalent to limiting the gaps in orbital phase between observations, maximising the orbital phase coverage (within the constraint of the length of the observing campaign), and minimising the variance in fitted model parameters later derived from the observations (Koen 2024).

The bottom-left panel of Figure 2 shows the scheduled RV measurements for a simulated self-lensing system. There are ten RV observations in total, unevenly spaced so as to effectively distinguish between solutions.

Joint photometric-RV MCMC fit

Each simulated pair of photometric light curves and RV curves is then fitted via MCMC using the same setup as in § 2.3, using the final chain state from the photometric fit as the initial state for the joint fit. Photometry is binned as per the first fit, while RV data are not binned. The Metropolis–Hastings random-walk proposal scale is reduced based on the number of RV measurements, when RV data are present, or otherwise on the number of photometric points, to allow the sampler to efficiently explore the narrower posterior where more data are present.

For a trial parameter vector $\theta$, let $v_{\mathrm{mod}}(t_k;\theta)$ be the model RV at the $k$th RV observation time, and let $N_{\mathrm{RV}}$ be the number of RV measurements. The RV contribution to the fit statistic is

where $v_k$ and $\sigma_{\mathrm{RV},k}$ are the simulated RV measurement and uncertainty. For the light curve, let $f_j$ and $\sigma_{\mathrm{LC},j}$ be the binned photometric flux and uncertainty at the $j$th binned time, let $f_{\mathrm{mod}}(t_j;\theta)$ be the corresponding model flux, and let $N_{\mathrm{LC}}$ be the number of binned lightcurve bins. The light-curve contribution is

The quantity sampled by the MCMC is then the sum of these terms plus the parameter-prior term, $\ln\pi(\boldsymbol{\theta})$,

The reported fitted parameters are posterior summaries rather than single maximum-posterior samples: after burn-in, we retain the 256 samples with the highest finite posterior values for each fit and take the mean of the fitted parameters, using a circular mean for the periapsis phase angle.

An example joint photometric-RV fit is shown in Figure 2. The corner plot in the right hand panel shows that the true solution is well described by the posterior distribution with the highest posterior model being close to the true values for $a$, $M_{\mathrm{lens}}$, $e$, and $\omega$, with only $i$ fitting to an incorrect value.

Predicting recovery of real self-lensing systems

In order to assess how the RV observing strategy would affect the recovery of the physical parameters of a physically realistic sample of self-lensing systems, we seek to combine our estimates of the recoverability of the physical properties of self-lensing systems with the detectability estimates of Wiktorowicz et al. (2021). While the present work deliberately assesses a sample which is not physically realistic, Wiktorowicz et al. (2021) uses population synthesis to generate a physically plausible distribution of physical properties, as the basis for simulating the total number of detectable systems observed by each photometric survey.

In the present work, we do not simulate a physically representative distribution of parameters, because our primary aim is to explore the ability to recover those parameters regardless of how common a given configuration will be. In other words, we want to simulate very rare configurations as thoroughly as more common ones. We are also concerned only with characterisation of lenses/recovery of physical parameters via spectroscopic RV measurements after detection in photometry, while Wiktorowicz et al. (2021) investigates photometric detectability.

Still, it is useful to derive from our simulations an estimate of how well we can recover the properties of a more realistic sample of self-lensing systems. To this end we build a nearest-neighbour emulator ($k=25$) to map the population synthesis results onto the fitting/recovery results for our radial velocity simulations.

We select population-synthesis entries marked as detectable and draw a sample of 10,000 entries with replacement, with selection probabilities proportional to their population-synthesis weights. This incorporates the relative frequencies represented by those weights into the sampled population, so each sampled entry contributes equally to the final statistics.

The emulator estimates the recovery probability for each population-synthesis system from the fitted outcomes of nearby systems in our simulation library. Here $i$ indexes our simulations and $j$ indexes the population-synthesis systems. We assign one scaled feature vector to each individual system in either sample, using its true physical parameters:

For each of the four coordinates, the centring constant $\mu$ is the median and the scale factor $s$ is the difference between the 84th and 16th percentiles, calculated from the combined simulation library and detected population-synthesis sample. For mass and period, these statistics are calculated after taking the logarithm. The same constants are applied to every system in both samples. This transformation shifts the origin and rescales each axis, placing both samples in a common dimensionless coordinate system and preventing parameters with larger numerical spreads from dominating the neighbour search.

For a population-synthesis system $j$, we calculate the Euclidean distance to each simulation $i$,

Thus, a vector describes an individual system, while a distance describes a pair of systems; the common centring constants cancel when taking the difference. We select the $k=25$ simulations with the smallest distances. The indices of the closest simulations form the set $\mathcal{N}_k(j)$, which can differ for each population-synthesis system. We assign these neighbours inverse-distance weights,

where $\epsilon=10^{-12}$ avoids division by zero, and the weights sum to unity over the selected neighbours. For each observing strategy (specified by observing baseline, number of RV measurements, and RV accuracy), let $q_i$ denote the compact-object class recovery outcome of simulation $i$. We set $q_i=1$ if the fitted compact-object class matches the true class and $q_i=0$ otherwise.

The emulated class recovery probability for population-synthesis system $j$ is the weighted mean of its neighbours’ class recovery outcomes,

The final class recovery fraction is the arithmetic mean of these emulated probabilities over the sampled population,

where $N_{\mathrm{pop}}$ is the number of sampled population-synthesis entries.

We also emulate the relative error in lens mass. For simulation $i$ and the specified observing strategy, we define the relative error in fitted lens mass as $r_i=\lvert M_{\mathrm{lens},i}^{\mathrm{fit}}-M_{\mathrm{lens},i}^{\mathrm{true}}\rvert/M_{\mathrm{lens},i}^{\mathrm{true}}$. We then interpolate these errors using the same neighbour weights,

The reported median relative error is the median of the $\hat{r}_j$ values for sampled entries within each true lens-mass bin.

We note that the survey-specific detection yields we use here come from Table 2 in Wiktorowicz et al. (2021), which is an optimistic upper limit including events for which only one photometric detection is made (i.e. one photometric point in the light curve falls within the lensing event). In reality, detections would require multiple photometric points to estimate $\tau_0$ and $\tau_{\mathrm{eff}}$.

Results and discussion

We evaluate the fits based on their ability to reproduce the true simulated parameters, with particular interest in lens mass and orbital period. The following subsections evaluate the fitting results in terms of of the median relative error in lens mass and orbital period. We also evaluate recovery of the correct compact-object class regardless of the accuracy of the fitted mass, considering if we can distinguish black holes from white dwarfs and neutron stars, as this distinction affects planning for further follow-up observations. We then evaluate the efficacy of fitting photometry alone versus fitting RV alone or jointly fitting photometry and RV, assessing how the number of RV measurements affects the relative error. Finally we use these results to emulate recovery rates for the synthesised population from Wiktorowicz et al. (2021), to estimate the relative errors in parameters for the real population of self-lensing systems.

It is worth noting that these results represent an ideal best case scenario for obtaining follow-up observations. We aim to assess a representative range of observing schedules, rather than rigorously simulating realistic observations. We do not account for observing seasons, we assume every target is always observable, we do not simulate weather or bad seeing, and we assume only a simple model of noise in the RV measurements (varying amounts of Gaussian noise to simulate spectrograph accuracy, ranging from sub-metre per second accuracy up to kilometre per second accuracy, to represent the range of spectrograph facilities currently available), plus a fixed 10 m/s Gaussian stellar noise. We also assume the exact time of the lensing peak ($t_0$) and the properties of the stellar companion (mass, radius, effective temperature, and distance from Earth) are known with no uncertainty; in practice these would all be subject to some amount of uncertainty.

Lens mass relative error

Figure 3 shows that after performing joint photometric/RV fitting, the results for higher lens masses generally have lower median relative errors. The number of observations has only a modest effect on the recovery of lens mass. For lower lens masses, the RV accuracy ($\sigma_{\rm ins}$) is an important factor, with $\sigma_{\rm ins} = 100$ m s$^{-1}$ and 1 km s$^{-1}$ accuracies producing substantially worse results. Interestingly, better RV accuracy does not always result in better recovery of lens mass; for lens masses $> 5\,M_\odot$, $\sigma_{\rm ins} = 100$ m s$^{-1}$ is better, while for the lowest mass lenses 10 m s$^{-1}$ is best. The presumptive explanation for this is that a $\sigma_{\rm ins}$ smaller than the simulated stellar noise can confound the fitting procedure. Future work could incorporate a GP-based noise model into the RV likelihood function to improve recovery further.

Recovering object class

While the ultimate goal of our observing campaign will be to accurately measure the lens mass and other physical properties of the population of self-lensing systems, we anticipate having to grapple with many orders of magnitude more false positives than true self-lensing systems. This situation will therefore require an informed, systematic way of triaging and prioritising targets for follow-up. The volume of false positives is anticipated to be so large that we will want to obtain as few observations of each target as possible in order to characterise what type of compact object the lens is (or, indeed, whether each target is even a self-lensing system to begin with). In the present work we assume every target is some kind of self-lensing system, but even with that assumption, some classes of object will be more astrophysically interesting than others, and some will require more or fewer observations in order to confirm. It is therefore informative to know whether we can accurately predict the class of the lens, even in cases where we can not accurately predict the mass. Being able to do so will allow us to allocate resources accordingly when we plan follow-up observations.

Figure 4 shows how well the object class is recovered as a function of the number of observations when performing a joint photometric/RV fit. A system with a fitted mass over 5 $M_\odot$ can often be identified as a black hole even when its fitted mass has a large error. This can inform our priorities when planning observations; e.g. we can choose to obtain additional observations of targets with fitted masses that fall in the 3–5 $M_\odot$ mass gap, or we can treat higher mass systems as confirmed with fewer observations. We find that the number of observations has only a small effect on the recovery of the object class, while RV accuracy is more important for lower mass systems but not for high mass systems. This suggests a staged observing strategy would be appropriate, where lower accuracy observations are first obtained for all targets, with higher accuracy observations reserved for those which seem to be lower mass.

Orbital period relative error

The orbital period strongly influences the shape of the RV curve. In the circular-orbit limit, the RV is zero at the lensing event and reaches an extremum half an orbit later; for eccentric systems, eccentricity and periapsis phase shift the detailed RV shape and extrema. When we can fix the orbital phase at a specific time (as we can when lensing is observed), obtaining an RV curve, even over a partial orbit, can tightly constrain the orbital period, eccentricity, and periapsis phase angle even for very long orbits.

Figure 5 shows the median relative error in orbital period as a function of number of observations, RV accuracy, and orbital period. Shorter orbital periods are recovered with smaller errors than longer ones; this is especially true for the shortest bin, where many systems will complete multiple orbits during the observing window.

Photometry-only and RV-only vs joint fits and $N_{\rm RV}$

The simulations are generated on the assumption that we photometrically detect a candidate self-lensing event and then follow up with RV measurements some time later. The joint photometry and RV fits use the photometric peak time and the targeted schedules described above. For comparison, we also carried out RV-only fits using independently generated, evenly spaced observations and fit the RV data in isolation, with a free orbital phase and no photometric information.

Figure 6 shows the median relative error in the recovered values of $M_{\rm lens}$, $a$, and $P_{\rm orb}$ as a function of the number of RV measurements. In all three cases, an RV-only fit performs worse than a photometry-only fit, and a joint fit has the lowest relative error. The median relative error in RV-only and joint fit converges after 5-6 observations.

Figure 7 shows the median improvement of fitted parameters broken down by the lens mass. The relative error is comparable for each $M_{\rm lens}$ bin when fitting $a$ and $P_{\rm orb}$ (with a modest improvement for larger masses). However, when fitting $M_{\rm lens}$ we find that higher lens masses can be fitted more accurately than lower masses. The improvement in the relative error when doing a joint fit is smaller for higher lenses, and in fact the RV-only fit is worse than the photometry-only fit at higher masses, suggesting RV measurements are the most useful for lower mass lenses. For the highest lens masses, fitting photometry alone results in recovered lens masses almost as accurate as a joint fit.

Figure 8 shows the percentage of systems which fall into the correct compact object mass bin regardless of the accuracy of the fitted mass. There is a clear trend where higher mass objects can be more accurately determined to be such, even if their exact mass can not be constrained. A higher mass black hole creates such a prominent lensing curve that it can be identified more easily, whereas neutron stars and white dwarfs are more ambiguous. It may be surprising that the photometry-only fits look so much worse for low mass lenses in Figure 8 than in Figure 7; this is because Figure 7 shows the relative error, with a median error close to 500% for the smallest mass bin. An error of that scale is more likely to result in lower mass systems falling in the wrong mass bin than larger mass systems (which have wider bins in this analysis). This is important astrophysically, as it makes it harder to distinguish neutron stars from mass gap black holes with photometry alone; for such candidates it is therefore more important to obtain RV measurements in order to characterise the nature of the lens.

Predicted recovery rates

We used the nearest-neighbour emulator to estimate the recovery rates for parameters based on the population synthesis-based parameter distribution. Recovery rates are the predicted percentages of detectable systems we could recover, calculated by emulating the results of our RV simulations across the physically-motivated population-synthesis distribution used in Wiktorowicz et al. (2021). These recovery rates are independent of photometric survey in this model, but still depend on the RV observing strategy. Figure 9 shows the recovery rates for compact-object class. Class recovery spans roughly 60–100% depending on number of RV observations, with higher mass lenses being recovered more successfully than lower mass.

Figure 10 shows the median relative errors for lens mass. Larger lens masses are recovered with smaller errors than smaller masses, and errors are significantly worse for $N_{\mathrm{obs}} < 5$, especially for lower lens masses.

These results follow the same trends as the unscaled results, showing that the unscaled results are informative about the potential for recovering the properties of the real population, even though the underlying parameter distribution is not physically representative.

Conclusions

Previous works have shown that self-lensing systems can be confirmed from photometry alone, or with a combination of photometry and radial velocity measurements, by the discovery of five relatively short period white dwarf lenses. We show it is possible to characterise lenses across the full compact object mass range, including for very long orbital periods where only one self-lensing event is detectable.

We have presented a proof-of-concept for a method of scheduling a small number of RV follow-up observations of self-lensing candidates with previously-observed photometric lensing signals. We show that with only a small number of RV measurements it is possible to determine the broad compact object class in many cases and to measure the lens mass to an acceptable relative error. This can then inform prioritisation of further observations, if needed. With this in mind, our recommendation for following up on photometrically-detected self-lensing candidates is a staged follow-up strategy, where all candidates are observed a small number of times to determine which ones are truly binaries and what the lens class is in those cases. A second stage of follow-up observations can then concentrate on the true binaries to ascertain the lens mass to a higher precision.

We show that the highest mass lenses tend to have the lowest relative errors, even with large orbital periods. We also show that the number of observations matters more than RV precision in many cases. This is helpful when it comes to planning follow-up observations, as it widens the pool of facilities which can provide useful RV measurements. It is also promising that we are able to constrain the simulation parameters even for very long period systems. This means it is worthwhile following up on single-event detections as well as those which are seen to be periodic.

We show that a targeted radial velocity observing campaign, carefully timing observations to distinguish between possible solutions, is effective in recovering object properties with a minimal number of observations. It is possible to compute such a targeted schedule by fitting photometry alone.

Combining parameter recovery rates with previous population synthesis results suggests that it is possible to confirm the existence of a substantial number of self-lensing systems in surveys such as TESS, ZTF, LSST, and PLATO. The estimates presented here are optimistic upper limits for confirmed detections under ideal conditions, but do well for our ability to find self-lensing systems in these and other surveys.

Acknowledgements

The authors acknowledge the use of the Codex coding agent to assist in creating and debugging the code underlying these simulations.

This work was supported by the Science & Technology Facilities Council [Grant Reference No ST/Y001699/1].

Data availability

The Python code and LaTeX source for this paper are available on GitHub (adammc̀master/self_lensing_radial_velocity), including scripts for generating all simulated data.

Photometric model

As with the full COBALT-BLUE model, we include three photometric effects: self-lensing, ellipsoidal variations, and Doppler boosting. We use the orbital phase convention described in § 2.1, with the phase and true anomaly defined in § 2.4.

The source-lens distance is calculated according to the phase angle at a given time, meaning the source-lens distance changes during the microlensing event according to the eccentricity. Strictly speaking, the Einstein ring radius should also change as a function of orbital phase, but this effect is expected to be small and so is not accounted for in the present work. Instead, the radius at $t_0$ is used for all calculations.

We use the Einstein radius and impact-parameter convention of Wiktorowicz et al. (2021), generalised from the circular-orbit case by replacing the orbital separation with the instantaneous Keplerian lens-source separation,

where $r(t_0)$ is the lens-source separation at the lensing peak. For a circular orbit, $r(t_0) = a$ and this reduces to the expression used by Wiktorowicz et al. (2021). The true anomaly, $\nu(t)$, is obtained from Kepler’s equation using the mean anomaly $\mathcal{M}(t)$ and eccentric anomaly $E(t)$,

with $\mathcal{M}_0$ chosen such that $\phi(t_0)=0$. $\arctan2$ is the standard quadrant-aware inverse tangent. The true anomaly at the lensing peak is $\nu_0=-\omega$, so the corresponding eccentric anomaly and mean anomaly are

For this orbital phase, the true orbital lens-source separation is

The projected lens-source separation in Einstein-radius units is

where $u(t)$ is dimensionless. The self-lensing magnification is calculated with the integrated Witt–Mao finite-source formalism (Witt & Mao 1994). The implementation evaluates the closed-form uniform-disc magnification,

where $K$, $E$, and $\Pi$ are complete elliptic integrals of the first, second, and third kind. The coefficients $C_E$, $C_K$, and $C_{\Pi}$, and the elliptic-integral arguments $m$ and $n$, are functions only of the projected lens-source separation and $\rho_{\star}=R_{\star}/R_{\mathrm{E}}$. The disc-integrated magnification used in the simulated light curves is then evaluated by quadrature over a linearly limb-darkened stellar disc,

where $\gamma=0.6$ is the fixed linear limb-darkening coefficient. This is the value corresponding to Eddington’s grey atmosphere (Eddington 1926). It is a fiducial approximation; in analyses of individual candidates, the coefficient should instead be interpolated from atmosphere models using the measured stellar parameters.

Calculation of ellipsoidal variations and Doppler boosting are adapted from Sorabella et al. (2022), with the orbital separation and radial velocity generalised to the eccentric case. The ellipsoidal-variation term is

where $\beta$ is the ellipsoidal-variation coefficient and $g=0.32$ is the gravity-darkening coefficient for a Sun-like star (Lucy 1967).

The Doppler-boosting term is

where $\nu_{\mathrm{TESS}}$ is the TESS-band effective frequency and $v_r(t)$ is the source radial velocity. The combined model flux is

which makes the simplifying assumption that the ellipsoidal variations and Doppler boosting are of small amplitude.

This paper has been typeset from a \TeX/\LaTeX file prepared by the author.