# The impact of the free-floating planet (FFP) mass function on the event rate
for the accurate microlensing parallax determination: application to Euclid
and Roman parallax observation

Ban, M. University of Warsaw Astronomical Observatory, Warsaw, Poland

(Accepted 19th October 2023)

###### Abstract

A microlensing event is mainly used to search for free-floating planets
(FFPs). To estimate the FFP mass and distance via the microlensing effect, a
microlensing parallax is one of the key parameters. A short duration of FFP
microlensing is difficult to yield a parallax by the observer’s motion at a
recognisable level, so the FFP microlensing parallax is expected on the
simultaneous observation by multiple telescopes. Here, we approach the FFP
detection by considering a variation of FFP mass functions and the event rate
of accurately measured microlensing parallax. We used our FFP microlensing
simulator assuming a parallax observation between upcoming space-based
missions (Euclid and Roman) with full kinematics. As a result, we confirmed
that the event rate of accurately measured microlensing parallax (i.e. within
a factor of 2 uncertainty) does not simply follow the number of FFPs at a
given mass but the ratio of the FFP population per star. This is because the
population ratio determines the optical depth for a given mass and potential
sources. In addition, we found that the probability of the event that can
estimate the FFP mass and distance within a factor of 2 is not so high:
$\sim$40% of Earth-mass, $\sim$16% of Neptune-mass, and $\sim$4% of Jupiter-
mass FFP events under our criteria. The probability can be improved by some
technical approach such as using high cadence and observation in parallax more
than two observers.

Gravitational microlensing – Planet – Free-floating – Parallax – Simulation

††journal: Apj

##  1 Introduction

Gravitational microlensing observation is a method to find isolated planetary-
mass objects such as free-floating planets (FFPs). To find the FFP properties
via the microlensing effect, it is essential to measure a microlensing
parallax ($\pi_{E}$) which is a value of the stereo-visional observation
angle. The accuracy of the microlensing parallax measurement directly affects
the accuracy of the lens property estimation. However, the number of FFP event
observations is still not enough. In recent decades, microlensing surveys have
detected tens of thousands of microlensing events, of which there are only a
few published candidates (Mróz et al., 2018, 2019a, 2020; Kim et al., 2021;
Ryu et al., 2021; Gould et al., 2022). FFP microlensing events are difficult
to observe because of the short event duration and the small event rate. The
event rate issue has several views such as the telescope sensitivity limits,
FFP population, and the difficulty of the determination whether it is isolated
or bound in a wide orbit, and whether it is categorised into planets or low-
mass brown dwarfs (Ban et al., 2016; Bayo et al., 2017; Barclay et al., 2017;
Delorme et al., 2017; Hamolli et al., 2013, 2016; Henderson & Shvartzvald,
2016; Mróz et al., 2017; Gould, 2013; Gould & Horne, 2013; Pacucci et al.,
2013; Penny et al., 2017; Sumi et al., 2011; Yee, 2013).

The FFP population has been discussed in recent decades. Sumi et al. (2011)
first argued that the population of Jupiter-mass FFP is $\sim$1.8 times larger
than that of the main-sequence (MS) stars. Later, Clanton & Gaudi (2017) and
Mróz et al. (2017) suggested the overestimation of the Jupiter-mass FFP
population by calculation and observation, respectively. They provided the
population of $\sim$1.4 and $<$0.25 times larger than the MS star population,
respectively. Thus, the prediction of the FFP population has been drastically
changed. The overall FFP population is still not settled, but some recent
research (Gould et al., 2022; Johnson et al., 2020; Sumi et al., 2023)
partially agree with each other; especially for Jupiter-mass and Earth-mass
FFPs. Johnson et al. (2020) simulated the upcomming space-based mission Nancy
Gracy Roman Space Telescope (Roman) and derived the expected FFP population,
and Gould et al. (2022) and Sumi et al. (2023) evaluated the latest FFP
population from KMTNet data and MOA data, respectively. Aside from the
microlensing approach, Scholz et al. (2022) predicted the FFP population from
the star-forming clusters observed by James Webb Space Telescope that up to 5%
of the population in a cluster could be massive FFPs ($>5M_{Jupiter}$).

Parallax detection helps to estimate the lens properties such as mass and
distance (Refsdal, 1966). The observer’s motion yields a parallax effect, and
it appears as the deformation or asymmetry of the light curve. However, the
parallax effect by an observer’s motion is likely unrecognizable for FFP
microlensing events due to its short duration. Koshimoto et al. (2023) also
concluded the difficulty of recognising parallax from the FFP microlensing
light curve candidates observed by MOA. Hence, simultaneous parallax
observation by multi-observers is more effective for FFP microlensing events.
The upcoming space-based missions (Euclid and Roman) can yield a simultaneous
parallax observation between them and between any ongoing ground-based surveys
such as MOA, OGLE, and KMTNet. Some papers, including our previous research,
discuss the FFP event configurations including the event rate and the parallax
sensitivity (Bachelet et al., 2018, 2022; Ban et al., 2016; Ban, 2020; Hamolli
et al., 2013, 2016; Johnson et al., 2020; Penny, 2014; Spergel et al., 2013).
Particularly, Bachelet & Penny (2019); Bachelet et al. (2022) simulated the
capability of parallax and mass measurements between Euclid and Roman assuming
the correspondence of the FFP and bound planet population through detectable
microlensing event and two well-known predictions of FFP population: two Dirac
delta functions with two peaks at $10M_{\oplus}$ (Mróz et al., 2019a) and
$2M_{Jupiter}$ (Sumi et al., 2011) per main sequence star. In this paper, we
assume the parallax observation between Euclid and Roman with variations of
the FFP mass function referred from Gould et al. (2022); Johnson et al.
(2020); Sumi et al. (2023) and try to approach the relationship between the
FFP population, event sensitivity, and accuracy of the microlensing parallax
measurement. As a similar research, Hamolli et al. (2017) simulated the FFP
microlensing detectability between Kepler and Spitzer by varying the FFP mass
function, but the combination of these two telescopes seems to be not suitable
for research on the low-mass FFPs because of the wide separation and low
observation cadence.

The brief process of our research is as follows. First, we generate the sample
artificial events using our microlensing event simulator. Second, MulensModel,
one of the open-source microlensing light curve analysis tools (Poleski & Yee,
2019), is applied for all sample events to fit the light curve. Third, the
microlensing parallax value is found using the result of fitted parameters of
the light curve. The event rate of accurately measured microlensing parallax
is derived for different FFP mass functions. Hence, the structure of the paper
is also following the process. §2 describes our microlensing event simulator.
We summarise the different FFP mass functions used in our research in §3 and
the generalised light curve fitting to approach the event rate of the accurate
estimation of the microlensing parallax in §4. The result is in §5, and we
analyse the relationship among the estimation accuracy of the microlensing
parallax, FFP event rate, and FFP mass functions. In §6, we discuss some
additional aspects of our results. Finally, in §7, we summarise our results
and future research objectives.

Table 1: Survey settings of Euclid and Roman for our simulator. The Euclid configuration is taken from Table 2 of Penny et al. (2013) and Agency (2011), and the Roman configuration is from Spergel et al. (2013, 2015); Webster & Folta (2017). We treat sample source flux in the Johnson-Cousins photometric system throughout our simulation. Therefore, a proxy band (in the parenthesis of the Filter row) is assumed. The threshold amplitude ($A_{t}$) to make an event alert corresponds to the threshold impact parameter of $u_{t}\sim 3$ for the point source case. | Euclid | Roman
---|---|---
Location | Lagrangian Point 2 | Lagrangian Point 2
Orbital diameter [km] | $\sim 0.97\times 10^{6}$ | $\sim 1.5\times 10^{6}$
Orbital period [year] | 0.5 | 0.5
Filter | NIPS $H$ ($H$) | W149 ($H$)
$A_{\rm t}$ | 1.02 | 1.02
$m_{\rm sky}$ [mag$/$arcsec2] | 21.5 | 21.5
$\theta_{\rm psf}$ [arcsec] | 0.4 | 0.4
$m_{\rm zp}$ | 24.9 | 27.6
$t_{\rm exp}$ [sec] | 54 | 52
Cadense [min] | 20 | 15
SAA [deg] | 90-120 | 53-126

##  2 Microlensing event simulator

###  2.1 Event model

The basic structure of our microlensing simulator is almost the same as Ban et
al. (2016); Ban (2020). We assumed the noise parameters from the Besançon
Galactic model version m1612 (Marshall et al., 2006; Robin et al., 2003,
2012a, 2012b, 2017) and used the catalogues to draw sample source and lens
properties except for replacing the lens mass to FFPs. The range of FFP mass
is $10^{-6}M_{\odot}-10^{-2}M_{\odot}$ which roughly corresponds to 0.3 Earth-
mass up to 10 Jupiter-mass FFPs. The target field is
$(l,b)=(1.^{\circ}00,-1.^{\circ}75)$ where is expected to yield more events
(Penny et al., 2013; Ban et al., 2016).

The full kinematics of the source, lens, and observers are considered in our
simulator. Therefore, the parallax and xallarap are always included during the
simulation process. In this research, we assume both source and FFP lens have
no companions, so xallarap does not occur. Even though the parallax due to the
observer’s motion is simulated, we expect the parallax effect to be negligible
for the FFP event duration.

Refsdal (1966) discussed the two solutions of microlensing parallax and Gould
(2020) generalised the probability of two solutions. In the projected lens
frame, there are two possible paths of the source that have the same minimum
impact parameter (See also Figure 1 and 2 from Gould (1994)). In our
simulator, the two solutions of microlensing parallax are attributed to the
net motion of the source, lens, and observer. Instead, either of the two paths
of the source in the lens frame is randomly selected. We assume the
probabilities of two paths are equal whilst the impact parameter value is
quadratic from the centre of the lens.

###  2.2 Telescopes

Table 1 summarises the conditions and configurations of Euclid and Roman in
our microlensing event simulator. The Euclid orbit is assumed to take a
maximum insertion angle ($\sim 33^{\circ}$) so that Euclid can fully take
advantage of the solar aspect angle (SAA). As a result of these
configurations, the cooperation period of Euclid and Roman lasts for $\sim$40
days around equinoxes. Euclid has optical filters (VIS) and infrared filters
(NISP) (Agency, 2011), and the infrared filter is much more suitable for FFP
microlensing observation because it can cover the dimmer source regime and
lower magnification events (Bachelet et al., 2022). Roman is expected to use
infrared filters for the exoplanet search via microlensing events (Johnson et
al., 2020; Spergel et al., 2013, 2015). Therefore, we assume the infrared
filter for both Euclid and Roman. Since our stellar catalogues from Besançon
Galactic model have the stellar magnitudes in Johnson-Cousin filters, we
regard the $H$-band filter as a proxy of the Euclid’s and the Roman’s infrared
filters.

In the simulator, we manage the Euclid and Roman kinematics by phasing the
orbital motion from $-\pi$ to $\pi$ in the sky frame centring at the
Lagrangian Point 2 (L2). When initialising the event properties, both Euclid
and Roman phases are randomly selected. The same phases between them are still
acceptable because of the different orbital diameters. These selected phases
are for the reference time ($t$=0) at which the impact parameter is also set
up. We use Euclid as a reference observer because of the smaller SAA and
brighter zero-point magnitude than those of Roman. Thus, the Euclid
sensitivity likely determines the simultaneous parallax detectability between
them.

###  2.3 Event criteria

We set the criteria of event alerting that the magnified flux results in the
signal-to-noise ratio ($SN$) $>50$. The noise is parameterised in four ways:
the background noise from unresolved stars ($m_{\rm stars}$), the blending
noise from nearby stars ($m_{\rm blend}$), the detector noise by a point
spread function (PSF) and sky brightness ($m_{\rm sky}$), and the photometric
noise generalised by a square-root of source flux. Hence, the $SN$ equation
becomes

$S/N(t)=\frac{10^{0.2m_{\rm zp}}\,t^{1/2}_{\rm
exp}\,A(t)\,10^{-0.4m_{*}}}{\sqrt{10^{-0.4m_{\rm stars}}+10^{-0.4m_{\rm
blend}}+\Omega_{\rm psf}10^{-0.4m_{\rm sky}}+A(t)\,10^{-0.4m_{*}}}},$ (1)

where $A(t)$ is the amplitude for a given time. $m_{\rm zp}$, $t_{\rm exp}$,
and $\Omega_{\rm psf}$ are the zero-point magnitude, exposure time, and the
solid angle of the survey point spread function (PSF). All parameters except
for the source baseline magnitude ($m_{*}$) are unique values depending on the
telescope sensitivity (see Table 1). The unresolved background noise ($m_{\rm
stars}$) and nearby bright star blending ($m_{\rm blend}$) are found by the
same method as Ban (2020) using Besançon catalogue. The brief review of the
calculation of $m_{\rm stars}$ and $m_{\rm blend}$ is as follows.

For a given source star ($j$), a boundary between resolved and unresolved
magnitude is found. Supposing if the boundary magnitude is $m_{j_{lim}}$, then
the total flux of stars fainter than the boundary magnitude ($i$) is $m_{\rm
stars}$:

$B_{res,lim}=\textstyle\Omega_{\rm
psf}\sum\limits_{m_{i}>m_{j_{lim}}}\frac{10^{-0.4m_{i}}}{\textstyle\Omega_{\rm
cat,i}}=10^{-0.4m_{\rm stars}},$ (2)

where $\Omega_{\rm cat}$ is the solid angle of the Besançon data to adjust the
population for different magnitudes. To find $j_{lim}$, we set the criteria
that the given source star flux ($F_{j}$) over the noise from the total flux
of unresolved stars ($\sqrt{B_{res}}$) is $>3$:

$\frac{F_{j}}{\sqrt{B_{res}}}=\frac{10^{0.2m_{\rm zp}}\,t^{1/2}_{\rm
exp}\,10^{-0.4m_{j}}}{\sqrt{\textstyle\Omega_{\rm
psf}\sum\limits_{m_{i}>m_{j}}\frac{10^{-0.4m_{i}}}{\textstyle\Omega_{\rm
cat,i}}}}>3.$ (3)

The blending noise by the nearby resolved star is also found using Eq.(3). For
the blending, the total flux of the resolved stars within the PSF of the given
star ($j$) becomes $m_{blend}$:

$B_{blend,j}=\textstyle\Omega_{\rm
psf}\sum\limits_{nearby}\frac{10^{-0.4m_{i}}}{\textstyle\Omega_{\rm
cat,i}}=10^{-0.4m_{\rm blend}}.$ (4)

Note that the blending noise is a unique value for each source star in the
catalogues per observer whilst the unresolved background noise is a unique
value per observer.

In addition to the alerting criteria of $SN>50$, we apply the least event
duration limit of 1 hour. This means that the observed epochs that keep
$SN>50$ must be lasting at least 1 hour in total. Note that the event duration
is different from the Einstein timescale ($t_{E}=\theta_{E}/\mu_{rel}$). Once
the event satisfies these criteria and is determined as an observable event,
the weight of the event is calculated based on the event occurrence potential.
We set a formula for the weight of the event as

$W_{SL}=u_{{\rm t},SL}D_{L}^{2}\mu_{SL}\theta_{{\rm E},SL},$ (5)

where $SL$ stands for a given source-lens pair, $u_{t}$ is the threshold
impact parameter derived from the threshold amplitude ($A_{t}$, see Table 1),
$D_{L}$ is the lens distance, $\mu_{SL}$ is the source-lens relative proper
motion, and $\theta_{E}$ is the Einstein radius (Wambsganss, 1998). As our
previous research shows (Ban et al., 2016), $u_{t}$ varies between a point
source and a finite source. For the point source event, there is a simple
mutual relation between $A_{t}$ and $u_{t}$,

$A(t)=\frac{u(t)^{2}+2}{u(t)\sqrt{u(t)^{2}+4}}.$ (6)

For the finite source event, Eq.(6) is no longer available. Our simulator uses
the integration method over the apparent surface of the source for the finite
source $A_{t}-u_{t}$ relation:

$A(u(t))=\frac{1}{\pi\rho^{2}}\int^{2\pi}_{0}\int^{\rho}_{0}A[s(u(t),r,\phi)]r\dot{d}rd\phi,$
(7)

where $\rho$ is the angular radius of the source in units of $\theta_{E}$. It
is also necessary to set up the criteria for parallax detectability in
addition to event detectability. The parallax detectability is based on the
differential light curve between two observers. the differential light curve
is recognisable from any noises and uncertainties when either of the following
conditions is satisfied:

$\left[\frac{dA(t)}{\sigma_{dA(t)}}\right]_{max}>5,$ (8)
$\frac{dt_{0}}{\sigma_{dt_{0}}}>5,$ (9)

where $dA(t)$ is the differential amplitude observed in two telescopes at a
given time $t$, $dt_{0}$ is the time difference of the peak amplitude between
two observers, and $\sigma$ is the uncertainty of those values. Figure 1 shows
one of the sample parallax events generated from our microlensing simulator.

For the FFP events, the finite source effect is never ignorable because of the
smaller $\theta_{E}$ and larger $\rho$ than those of the stellar lens case
(Ban et al., 2016; Ban, 2020). As Eq.(7) shows, the amplitude with the finite
source effect depends on the source surface brightness. The source surface
brightness is usually not uniform, so we test the effect of limb darkening on
the integration of the amplitude for the finite source effect. Figure 2 shows
a comparison of the $\chi^{2}$ values between the assumption of uniform
brightness (UNI) and the limb-darkening coefficients (LDC). The colour bar
indicates the ratio of $\chi^{2}$ derived as

$\frac{|\Delta\chi^{2}|_{m}}{\chi^{2}_{v}}=\frac{|\chi^{2}_{m,LDC}-\chi^{2}_{m,UNI}|}{\chi^{2}_{v,UNI}},$
(10)

where the lower-case of $m$ stands for the discrepancy of data from a model
whilst $v$ stands for the weighted sum of the squared variation. That is to
say, $\chi^{2}_{m}$ presents how the theoretical light curve differs from the
point-source-uniform-brightness source case, and $\chi^{2}_{v}$ presents the
uncertainty of the observed photometry. As we can see in Figure 2, the light
curve difference of the uniform brightness and limb-darkening cases is much
smaller than the photometric uncertainty in the majority of events for FFP
lenses. We also found that the $T_{eff}<3300$ [K] regime did not yield the
finite source effect under the test conditions as Figure 2. Such a cool source
likely be a low-mass dwarf having a small physical radius, so $\rho$ got small
enough to be a point source. From this limb-darkening test, we finally
determined that limb darkening is ignorable for our research of the accuracy
of the microlensing parallax measurement. In Figure 2, we added the cross
marker to briefly show the distribution of our sample events for the parallax
measurement research that were assumed a uniform brightness case. The centre
of the marker is at
$(T_{eff},u_{0})\sim(4127^{+1044}_{-497}K,1.70^{+0.88}_{-1.31})$ with the
$\chi^{2}$ ratio of $\sim 5.75\times 10^{-7}$.

Figure 1: An example of the photometric data observed by Euclid and Roman in
simultaneous parallax generated by our microlensing simulator.
Figure 2: The map of the ratio of the chi-square values for the different
effective temperatures of a source and an impact parameter (see Eq.(10)). The
limb-darkening coefficients are referred from Claret et al. (2012); Claret
(2017). The finite source effects are simulated with the mean angular radii of
the source derived from the Besançon Galactic model catalogues assuming FFP
lenses. The cross marker represents the weighted mean and $1\sigma$
distribution of $T_{eff}$ and $u_{0}$ of our sample events. Note that the
uniform brightness was assumed for the sample drawing, and these sample events
survived our detectability criteria (see §2.3).

###  2.4 Microlensing parallax

The theoretical formula of the microlensing parallax ($\pi_{E}$) (Refsdal,
1966; Gould, 1992) is

$\pi_{E}=\frac{AU}{\tilde{r_{E}}}=\sqrt{\kappa M_{L}D_{rel}},$ (11)

where

$\kappa=\frac{4G}{c^{2}},$ (12)

and

$\frac{1}{D_{rel}}=\frac{1}{D_{L}}-\frac{1}{D_{S}}.$ (13)

In the real observation, the value of microlensing parallax is estimated
through the light curve fitting process, and the result is shown in 2D vectors
in the projected lens frame. The vector components are perpendicular to each
other and the line of sight. When the parallax effect is attributed mostly to
the telescope separation, the vector components are set in the direction of
the relative proper motion ($t$-component) and the direction of the minimum
impact parameter ($u$-component). In this case, the parallax formula (Eq.(11))
is rewritten as

$\bm{\pi_{E}}=\frac{AU}{D_{T}}\left(\frac{\Delta t_{0}}{t_{E}},\Delta
u_{0}\right),$ (14)

where $D_{T}$ is the projected observer separation in units of AU
perpendicular to the reference observer’s line of sight, $t_{0}$ is the time
when the amplitude peaks, $t_{E}$ is the Einstein timescale, and $u_{0}$ is
the minimum impact parameter.

In this paper, we derive the theoretical $\pi_{E}$ with Eq.(11) using the
values used for the event samples. The estimated $\pi_{E}$ is derived from the
fitted parameter of the light curve and Eq.(14). The estimation accuracy of
the microlensing parallax is calculated by taking the discrepancies of the
estimated $\pi_{E}$ from the theoretical $\pi_{E}$ as

$\Delta\pi_{E}[\%]=100\times\frac{|\bm{\pi_{E,est}}|-\pi_{E,th}}{\pi_{E,th}}.$
(15)

The weight of each sample event and any analytical constraints provide the
fraction of events ($FoE$). It is the weighted ratio of the events that
satisfy any analytical constraints over all events.

$FoE(limit)=\frac{\sum_{SL,limit}^{N_{limit}}f(M_{L})\cdot
W_{SL}}{\sum_{sl}^{n}{f(M_{l})\cdot W_{sl}}},$ (16)

where $f(M_{L})$ is the population factor to adjust the sampling bias for
different FFP mass functions. The detail of the population factor is described
in the next section (§3). The “limit” indicates any analytical constraints
such as the timescale, FFP mass, finite source effect, telescope phase
difference, and so forth. The lower case of $s$, $l$, and $n$ in the
denominator depends on the fraction type. If the fraction type is for all
sample events, then $s=S$, $l=L$, $n=N$. If the fraction type is for all
events (i.e. theoretically detectable by Euclid without considering the
parallax and event duration criteria), then the lowercase iterations include
both detectable and undetectable events under our parallax criteria. The
fraction of events derived for any analytical constraints finally provides the
event rate with the constraints as

$\Gamma_{limit}=\Gamma_{0}\times FoE(limit),$ (17)

where $\Gamma_{0}$ is the event rate of the theoretical solo observation. In
this paper, the basic FFP population of the theoretical solo observation is
that the 1 FFP per star and $\Gamma_{0}$ is derived with the same event
detectability criteria (i.e. $S/N>50$ and event duration $>1$ hour) described
in §2.3.

##  3 FFP Mass functions

The variation of the FFP population is taken as the mass function slope
($\delta N/\delta lnM\propto M^{-\alpha}$). We use the population of 1 FFP per
star as the standard mass function in our parallax simulation and label it as
the “Uniform” mass function case. According to the stellar distribution of the
Besançon Galactic model catalogues, the Uniform function corresponds to
$\delta N/\delta lnM\sim 3.7$ for our FFP mass range. In addition to the
Uniform function, we consider four more mass functions as listed below.

  * •

Uniform : $\delta N/\delta lnM\sim 3.7$

  * •

Johnson et al. (2020), Eq.(8) :
$\delta N/\delta lnM=(0.24/dex)\times(M_{FFP}/95M_{\oplus})^{-0.73}$ for
$M_{FFP}/M_{\oplus}\geq 5.2$
$\delta N/\delta lnM=100$ for $M_{FFP}/M_{\oplus}<5.2$

  * •

Gould et al. (2022), Eq.(10) with $p=0.9$ :
$\delta N/\delta lnM=(0.4/dex)\times(M_{FFP}/38M_{\oplus})^{-0.9}$

  * •

Gould et al. (2022), Eq.(10) with $p=1.2$ :
$\delta N/\delta lnM=(0.4/dex)\times(M_{FFP}/38M_{\oplus})^{-1.2}$

  * •

Sumi et al. (2023), Eq.(17) and Table 4 :
$\delta N/\delta lnM\propto(1.85/dex)\times(M_{FFP}/8M_{\oplus})^{-1.14}$ for
$M_{br}<M_{FFP}/M_{\odot}<0.02$
$\delta N/\delta lnM\propto M_{FFP}^{-0.13}$ for $-7<M_{FFP}/M_{\odot}<M_{br}$

where $m=log_{10}(M_{FFP}/M_{\odot})$, and the range is from -6 to -2 since
our sample event FFP mass range is $10^{-6}M_{\odot}-10^{-2}M_{\odot}$.
Johnson et al. (2020) estimated an FFP mass function with the FFP mass
boundary at 5.2$M_{\oplus}$ based on the bound planet mass function by Cassan
et al. (2012) and simulated the detectable event rate of FFP microlensing
using Roman configuration. Sumi et al. (2023) estimated the FFP mass function
from the short events observed by MOA. As Sumi et al. (2023) mentioned about
the total FFP mass per star, we need to adjust the event sampling bias between
different mass functions for our research. Hence, we apply a population factor
($f(M_{FFP})$) against the Uniform function case assuming that the population
factor forms a power law:

$f(M_{FFP})=\frac{N(M_{FFP})}{N(M_{*})}=\left(\frac{M_{FFP}}{M_{\odot}}\right)^{-\beta}.$
(18)

For the Uniform function case, the population factor is 1 so that the slope of
the power law ($\beta$) is 0 throughout the FFP mass. Figure 3 visualises
three mass functions and the slope of the power law of the population factor.

Figure 3: Three different FFP mass functions in the solid curves and the slope
of the population factor as a power law ($\beta$) in Eq.(18) in the dashed
curves. The left y-axis corresponds to the mass function values, and the right
y-axis corresponds to the $\beta$ value.

##  4 Fitting the light curve of sample events using MulensModel

MulensModel is one of the open-source tools for the microlensing light curve
fitting (Poleski & Yee, 2019). We use this tool to find the estimated
$\pi_{E}$ value for each sample event photometry. For all sample events, we
consider (1) single-source-single-lens event, (2) finite-source effect, and
(3) parallax by the observer’s motion. The MCMC method is applied to find the
best-fit parameters. Although MulensModel has a function to include ephemeris
of the space-based surveys, it works for a ground-based reference and space-
based partner observers case. Therefore, we use another approach for the
parallax by the motion of Euclid and Roman in MulensModel. Since we already
expect that the light curve torsion due to the parallax is negligible for the
FFP events because of the short event duration, we assume the motion of Euclid
and Roman during the observation of one event can be approximated to the shift
of the L2 point; hence, the orbital motion of Earth. On the other hand, the
impact parameter and the peak time of the light curve could largely differ
between Euclid and Roman, so we use the fitting steps as follows.

  1. 1.

Fitting both Euclid and Roman photometry together

  2. 2.

Fitting each Euclid and Roman photometry with a fixed value of $t_{E}$ and
$\rho$ from the step 1

  3. 3.

Checking reduced chi square ($\chi^{2}/d.o.f.$) for each Euclid and Roman

  4. 4.

Recording $t_{0}$, $u_{0}$, and $\pi_{E}$ for each Euclid and Roman from the
step 2

  5. 5.

Going back the step 1 with average $t_{E}$ and $\rho$ from the step 2

  6. 6.

Repeating the step 1-5 until $\chi^{2}/d.o.f.<2$ for each Euclid and Roman

In the end, we gain a final record of $t_{0}$, $u_{0}$ and $\pi_{E}$ for each
Euclid and Roman individually and common $t_{E}$ and $\rho$ between them.

It is essential to derive both the microlensing parallax and the angular
source radii from the light curve to estimate the lens properties. We assume
the theoretical constraint of the recognisable finite source effect is
$\rho/u_{0}<$1.0, and the sensitivity-based constraint is $\rho/u_{0}<$0.5
(Johnson et al., 2020) for the reference observer Euclid. We consider both
constraints of the angular source radii to calculate the accuracy of the
microlensing parallax measurement from the fitted parameters.

The significance of the consideration of (2) finite-source effect and (3)
parallax by the observer’s motion are also confirmed by comparing
$\chi^{2}/d.o.f.$ of with and without these considerations. As a result,
$\chi^{2}/d.o.f.$ changes more significantly for the low-mass FFPs and small
impact parameters. The comparison of the fitting results with and without the
finite source consideration also takes part in the finite source constraints
mentioned above. On the other hand, the observer’s motion is ignorable as we
expected and does not affect $\chi^{2}/d.o.f.$ and fitted parameters.

Although we optimised the parallax due to the observer’s motion and did not
include ephemeris in the light curve fitting process, we assume that the
space-based surveys usually have a good record of the telescope coordinates in
the sky. Our simulator yields the coordinates information for each detection
epoch. The separation between Euclid and Roman on the reference line ($D_{T}$)
is derived from these coordinates.

##  5 Results

###  5.1 Light curve fitting uncertainty

The procedure of the light curve fitting shown in §4 yielded the best-fit
parameters with uncertainties no greater than 20%. We also calculated the
discrepancy of the fitted parameters from the theoretical parameters used in
the simulator. Those theoretical parameters were used as the initial
parameters for the fitting process using MulensModel. Table 2 shows the
fitting parameter discrepancies and their likelihood among our samples for
Earth-, Neptune-, and Jupiter-mass FFPs. The same concept as Eq.(15) were used
to calculate the discrepancies for $t_{E}$, $\rho$, and $|u_{0}|$ for both
Euclid and Roman. The likelihood was the same as the fraction of the event
which the fundamental formula is shown as Eq.(16). The discrepancy of $t_{0}$
is removed from Table 2 because it reaches zero due to the MulensModel fitting
configurations, and the likelihood becomes $\sim$1\. This is because we
required all sample events to successfully observe the peak magnitude, so
MulensModel was given a true value of $t_{0}$, which also is directly readable
from the input photometric data. Unlike the other event parameters such as
$u_{0}$ and $t_{E}$, $t_{0}$ did not fluctuate so much during the fitting
process as long as a given parameter and read data from the photometric data
matches quite well.

Although we applied the theoretical parameters as the initial parameters of
the light curve fitting, the MCMC approach yielded the best-fit parameters
with some discrepancy from the theoretical values. For example, the Jupiter-
mass FFP case likely yields a large discrepancy of the $\rho$ estimation. This
is because $\rho$ of Jupiter-mass FFP generally becomes smaller than that of
the low-mass FFPs so the discrepancy becomes relatively large. The comparison
of the light curve fitting with and without finite source consideration does
not show a large $\chi^{2}/d.o.f.$ difference for massive FFPs. Thus, there is
a calibration limit of the small $\rho$ fitting, and the possible range of
$\rho$ likely increases. Aside from the $\rho$ discrepancy, $t_{E}$ and
$|u_{0}|$ are fitted well. The $t_{E}$ value of more than 80% events and the
$|u_{0}|$ value of around 90% events are fitted within 100% discrepancy. Here,
the 100% discrepancy corresponds to a factor of 2 uncertainty. Some sample
events show $>$100% discrepancy for the fitted parameters, and this is because
of the number of epochs, magnification level, and lack of peak observation.
Especially for the former two causes, those sample events passed our
detectability criteria, but very barely.

Table 2: Likelihood of the fitted parameters of the Euclid-Roman simultaneous microlensing observation for Earth-, Neptune-, and Jupiter-mass FFPs cases. The different discrepancy level ($<$10%, $<$50%, and $<$100%) are extracted here. The discrepancy calculation follows the same concept as Eq.(15) for each parameter’s theoretical and estimated values. The likelihood is derived as the fraction of events from the sample event weight over all samples. $E$ and $R$ notes for $|u_{0}|$ stand for Euclid and Roman. FFP mass | Discrepancy | $t_{E}$ | $\rho$ | $|u_{0,E}|$ | $|u_{0,R}|$
---|---|---|---|---|---
Earth |  $<$10% | 0.25 | 0.10 | 0.15 | 0.18
$<$50% | 0.66 | 0.34 | 0.53 | 0.58
$<$100% | 0.84 | 0.89 | 0.90 | 0.87
Neptune |  $<$10% | 0.35 | 0.12 | 0.23 | 0.29
$<$50% | 0.68 | 0.36 | 0.56 | 0.65
$<$100% | 0.82 | 0.73 | 0.94 | 0.92
Jupiter |  $<$10% | 0.33 | 0.05 | 0.22 | 0.24
$<$50% | 0.69 | 0.22 | 0.57 | 0.63
$<$100% | 0.85 | 0.48 | 0.93 | 0.92

###  5.2 Event rate and the accuracy of the microlensing parallax measurement

For a given FFP mass, the fraction of events for the simultaneous parallax
observations does not change over the different FFP mass functions. We found
that the probability of simultaneous parallax observation with finite source
effect was $\sim$40%, 16%, and 4% for Earth-mass, Neptune-mass, and Jupiter-
mass FFPs under our criteria, respectively. Our parallax detectability
criteria affect more on massive FFPs than low-mass FFPs since it focuses on
the differential light curve. The criteria result in limiting the impact
parameter of the Jupiter-mass FFP events. For instance, the differential
magnitude of $u_{0,a}=0.5$ and $u_{0,b}=0.7$ is smaller than that of
$u_{0,a}=0.1$ and $u_{0,b}=0.3$ even though the differences of the impact
parameters are both 0.2. The latter case yields a larger magnitude difference
and more possible to survive under our criteria. This tendency also relates to
the occurrence of the finite source effect for Jupiter-mass FFPs. The Einstein
ring of Jupiter-mass FFPs is generally large enough so it requires a small
impact parameter to yield a finite source effect. On the other hand, Earth-
mass FFPs yield a distinctive differential light curve more easily due to the
small Einstein ring, and the finite source events are dominant (Ban, 2020).

The event rate of accurately measured microlensing parallax varies with the
FFP mass function because the population of a given FFP mass changes in the
total FFP mass per star. Figure 4 shows the event rate for three different FFP
mass functions listed in §3 along the FFP mass from $10^{-6}M_{\odot}$ to
$10^{-2}M_{\odot}$. Two limits of finite source effect are plotted in the same
diagram, and the discrepancy of $<100\%$ is assumed as the least acceptable
level. Table 3 summarises the event rate of Earth-, Neptune-, and Jupiter-mass
FFPs for three different mass functions. In this table, we assumed that the
calculation of all events was done by a typical case (i.e. $\Delta
u_{0}=u_{0,a}-u_{0,b}$ where $a$ and $b$ indicate separated observers
(Refsdal, 1966)) since it is usually difficult to distinguish between two
solutions of the microlensing parallax from the light curve. In case the
solutions are distinguished will be discussed in the later section §6.

According to Figure 3, the FFP population density at $M_{FFP}\sim 10^{-3.8\pm
0.1}$ is almost agreed for all different FFP mass functions except for Gould
et al.’s function with $p=1.2$. However, the event rates for the mass in
Figure 4 do not agree with each other. This disagreement reflects the
different population ratios of the mass to the total FFP mass per star. In
other words, the optical depth for a given FFP mass lens and potential sources
varies and determines the event rate. The Gould et al.’s functions and Sumi et
al.’s function cross over the Johnson et al.’s function at some FFP mass in
Figure 3, and the FFP population densities roughly agree with each other. In
this case, we can regard that the optical depth for the FFP mass lens and
potential sources is almost equal between those mass functions.

Figure 4: Event rate [event/year/deg2] for accurately measured microlensing parallax for different FFP mass functions. The FFP mass on the x-axis is binned by $dlog_{10}M_{FFP}\sim 0.1$. The least acceptable discrepancy is assumed as $<100\%$, which corresponds to a factor of 2 uncertainty. The limitations of the finite source effect are $\rho/u_{0}<0.5$ in solid and $<1.0$ in dashed curves. Two solutions of the microlensing parallax are not considered in this plot and the $\delta u_{0}$ is taken from a typical subtraction formula. Table 3: Event rate [event/year/deg2] for accurately measured microlensing parallax for different FFP mass functions. The numerical values are extracted for the conditions of Earth-, Neptune-, and Jupiter-mass FFPs, and $<$10%, $<$50%, and $<$100% discrepancies. For each column of the mass function, the left sub-columns are for the limit of the finite source effect $\rho/u_{0}<$0.5, and the right sub-columns are for $<$1.0. Two solutions of the microlensing parallax are not considered in this plot and the $\delta u_{0}$ is taken from a typical subtraction formula. FFP mass | Discrepancy | Uniform | Johnson et al. |  Gould et al., $p=0.9$ |  Gould et al., $p=1.2$ | Sumi et al.
---|---|---|---|---|---|---
$<$0.5 |  $<$1.0 |  $<$0.5 |  $<$1.0 |  $<$0.5 |  $<$1.0 |  $<$0.5 |  $<$1.0 |  $<$0.5 |  $<$1.0
Earth |  $<$10% | 0.71 | 1.31 | 22.4 | 44.3 | 25.3 | 52.0 | 9.12 | 18.8 | 15.6 | 32.0
$<$50% | 2.68 | 6.24 | 87.0 | 208 | 101 | 242 | 36.7 | 87.2 | 62.4 | 149
$<$100% | 3.74 | 9.53 | 126 | 325 | 147 | 379 | 53.2 | 137 | 90.7 | 233
Neptune |  $<$10% | 12.2 | 21.1 | 15.1 | 21.3 | 17.0 | 23.9 | 3.56 | 4.45 | 10.5 | 14.7
$<$50% | 25.9 | 77.9 | 28.8 | 85.1 | 32.2 | 94.4 | 6.08 | 17.7 | 19.8 | 58.1
$<$100% | 30.2 | 95.5 | 32.2 | 103 | 36.0 | 114 | 6.72 | 21.3 | 22.1 | 70.1
Jupiter |  $<$10% | 24.4 | 69.5 | 4.86 | 13.7 | 1.36 | 3.15 | 0.15 | 0.30 | 0.84 | 1.94
$<$50% | 76.0 | 164 | 12.7 | 26.2 | 3.48 | 6.42 | 0.39 | 0.66 | 2.14 | 3.95
$<$100% | 94.4 | 198 | 15.8 | 31.9 | 4.19 | 7.84 | 0.47 | 0.82 | 2.58 | 4.83

##  6 Discussion

###  6.1 Comparison with the other research

Bachelet et al. (2022) considered the probability of FFP lens mass and
distance estimation assuming the mass function of 10 Earth-mass FFPs per Main-
sequence star (Mróz et al., 2019b) and 2 Jupiter-mass FFPs per Main-sequence
stars (Sumi et al., 2011). They showed that 110 Earth-mass and 450 Jupiter-
mass FFP microlensing events per year can be observed in parallax and $\sim$30
and $\sim$19 of them can determine the lens mass and distance within a factor
of 2 uncertainty. Under the similar FFP population and total field of view
(2.0 [deg2]) as Bachelet et al. (2022), our simulation indicates that 23
Earth-mass and 11 Jupiter-mass FFPs microlensing events per year are
observable in parallax and $\sim$12 and $\sim$2 of them can determine the mass
and distance within a factor of 2 uncertainty. Here, we assumed that the
accuracy of the fitted parameters (see §4) is improved, and 90% of events
resulted in all fitted parameters being within 100% discrepancy for any FFP
mass.

Our event rates are smaller than the values of Bachelet et al. (2022) because
our observation condition is more limited, and detectability criteria are more
strict than their configurations. Hence, to compare the results, we can
consider the ratio of the detectable number of events that yields the event
rate values above. Table 3 of Bachelet et al. (2022) shows the fraction of
events (SN only):(SN+parallax):(SN+finite) $\sim$1:0.85:0.24 for Earth-mass
and $\sim$1:0.92:0.04 for Jupiter-mass FFPs whilst our results are
$\sim$1:0.69:0.58 and $\sim$1:0.20:0.18, respectively. For the Earth-mass FFP
case, our criteria of parallax detectability are more strict than that of
Bachelet et al. (2022) whilst a more finite source case is allowed. This is
because our criteria take advantage of the extension of the threshold impact
parameter with large $\rho$ (see Ban et al. (2016) for the details). As a
total probability, Bachelet et al. (2022) shows $\sim$20% of the Earth-mass
FFP events and $\sim$4% of the Jupiter-mass FFP events are observed in
parallax, and the lens mass and distance can be determined in a factor of 2
uncertainty. Our corresponding results are $\sim$40% of the Earth-mass and
$\sim$4% Jupiter-mass FFP events (see the discussion in §5.2). The larger
probability of our result for the Earth-mass FFPs is due to the consideration
of the extension of the threshold impact parameter which results in more
detectable events satisfying the event duration limit of $>1$ hour. We also
confirmed that the average expanded threshold impact parameter for Earth-mass
FFPs is about twice as large as the point source case. Hence, it is plausible
that our probability is doubled from Bachelet et al. (2022). On the other
hand, the Jupiter-mass FFP generally has less advantage of such extension of
the threshold impact parameter, so the similar probability indicates that we
reached a similar result to Bachelet et al. (2022) using the different
approach of the event criteria.

###  6.2 Two solutions of microlensing parallax

As Refsdal (1966) explained, the two solutions of the microlensing parallax
are attributed to the direction of the impact parameters. If the impact
parameter vector is opposite between observers, then $\Delta
u_{0}=u_{0,a}+u_{0,b}$ would be true; where $a$ and $b$ indicate separated
observers. This equation makes the microlensing parallax larger than the
typical subtraction of $\Delta u_{0}=u_{0,a}-u_{0,b}$. We checked all sample
events in different FFP mass functions and found that the fraction of the
addition case ($\Delta u_{0}(+)$) to the typical subtraction case ($\Delta
u_{0}(-)$) is no greater than $1/10$ events for the FFP mass rage of
$10^{-6}M_{\odot}-10^{-2}M_{\odot}$. The addition case is occurable for all
FFP mass ranges in our research, and the low-mass FFPs are more possible. We
also derived the event rate for accurately measured microlensing parallax if
the solution is identified as either subtraction or addition cases for each
sample event. If it is distinctive in either the subtraction or addition case,
the event rate slightly increases from the case that all events are treated as
the subtraction case (i.e. the results shown in §5). However, the event rate
becomes almost the same when the discrepancy of the microlensing parallax is
accepted up to 100%. Hence, it will be an advantage if either typical
subtraction or addition calculation is distinctive, but it is not a big issue
for FFPs under our criteria because of the low occurrence rate of an addition
case.

##  7 Conclusion

We explored the influence of the different FFP mass functions on the event
rate for accurately measured microlensing parallax and found that the event
rate is related to the FFP population based on the extent of the finite source
effect. For instance, the finite source event is dominant for the Earth-mass
FFP events (Ban, 2020), and $\sim$40% of them can measure the FFP mass and
distance in a factor of 2 uncertainty. The probability becomes $\sim$16% and
$\sim$4% for Neptune-mass and Jupiter-mass FFPs, respectively. However, the
event rate does not simply rely on the number of FFPs for a given mass because
the ratio of the FFP population over the total FFP mass per star determines
the optical depth and event rate.

Aside from the influence of the different FFP mass functions on the event
rate, the accuracy of the microlensing parallax measurement can be improved to
some extent by improving the accuracy of the light curve fitting. Other than
the technical issue of the fitting, the simple way of improvement is to use a
higher cadence than we are expecting in our simulation or to have a multi-
parallax of more than two telescopes. The ground-based surveys such as MOA,
OGLE, and KMTNet can also be partners for the parallax observation with any
space-based telescopes (Ban et al., 2016; Ban, 2020). Moreover, there are some
other telescopes planned to launch for L2 orbit such as Ariel and PLATO
(Gardner et al., 2006; Rauer et al., 2014; Gaia Collaboration et al., 2016;
Agency, 2020). Although the microlensing event is not a main target of these
future surveys, they still have the potential to detect and alert microlensing
events (Mislis et al., 2016; Yee et al., 2017; Nikolaus & Hundertmark, 2018).
By operating these space-based telescopes one after another, it would also be
possible to overtake the parallax observation by the space-based surveys after
the co-operation period of the Euclid and Roman. It is necessary to gather the
observation data to determine the FFP mass function, so a long-lasting mission
as a whole is worthwhile.

##  8 Acknowledgement

Numerical computations were carried out on a PC cluster at the Center for
Computational Astrophysics, National Astronomical Observatory of Japan. The
work was supported by the Polish National Agency for Academic Exchange grant
“Polish Returns 2019” to Radosław Poleski.

## References

  * Agency (2011) Agency, E. S. 2011, Euclid : Mapping the geometry of the dark Universe.  https://arxiv.org/pdf/1110.3193.pdf
  * Agency (2020) —. 2020, Ariel : ESA/SCI(2020)1 Enabling planetary science across light-years European Space Agency.  https://sci.esa.int/documents/34022/36216/Ariel_Definition_Study_Report_2020.pdf
  * Bachelet et al. (2018) Bachelet, E., Hinse, T. C., & Street, R. 2018, AJ, 155, 191, doi: 10.3847/1538-3881/aab3c8
  * Bachelet & Penny (2019) Bachelet, E., & Penny, M. 2019, ApJ, 880, L32, doi: 10.3847/2041-8213/ab2da5
  * Bachelet et al. (2022) Bachelet, E., Specht, D., Penny, M., et al. 2022, A&A, 664, A136, doi: 10.1051/0004-6361/202140351
  * Ban (2020) Ban, M. 2020, MNRAS, 494, 3235, doi: 10.1093/mnras/staa786
  * Ban et al. (2016) Ban, M., Kerins, E., & Robin, A. C. 2016, A&A, 595, A53, doi: 10.1051/0004-6361/201629166
  * Barclay et al. (2017) Barclay, T., Quintana, E. V., Raymond, S. N., & Penny, M. T. 2017, ApJ, 841, 86, doi: 10.3847/1538-4357/aa705b
  * Bayo et al. (2017) Bayo, A., Joergens, V., Liu, Y., et al. 2017, ApJ, 841, L11, doi: 10.3847/2041-8213/aa7046
  * Cassan et al. (2012) Cassan, A., Kubas, D., Beaulieu, J. P., et al. 2012, Nature, 481, 167, doi: 10.1038/nature10684
  * Clanton & Gaudi (2017) Clanton, C., & Gaudi, B. S. 2017, ApJ, 834, 46, doi: 10.3847/1538-4357/834/1/46
  * Claret (2017) Claret, A. 2017, A&A, 600, A30, doi: 10.1051/0004-6361/201629705
  * Claret et al. (2012) Claret, A., Hauschildt, P. H., & Witte, S. 2012, A&A, 546, A14, doi: 10.1051/0004-6361/201219849
  * Delorme et al. (2017) Delorme, P., Dupuy, T., Gagné, J., et al. 2017, A&A, 602, A82, doi: 10.1051/0004-6361/201629633
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485, doi: 10.1007/s11214-006-8315-7
  * Gould (1992) Gould, A. 1992, ApJ, 392, 442, doi: 10.1086/171443
  * Gould (1994) —. 1994, ApJ, 421, L75, doi: 10.1086/187191
  * Gould (2013) —. 2013, ApJ, 763, L35, doi: 10.1088/2041-8205/763/2/L35
  * Gould (2020) —. 2020, Journal of Korean Astronomical Society, 53, 99, doi: 10.5303/JKAS.2020.53.5.99
  * Gould & Horne (2013) Gould, A., & Horne, K. 2013, ApJ, 779, L28, doi: 10.1088/2041-8205/779/2/L28
  * Gould et al. (2022) Gould, A., Jung, Y. K., Hwang, K.-H., et al. 2022, Journal of Korean Astronomical Society, 55, 173, doi: 10.5303/JKAS.2022.55.5.173
  * Hamolli et al. (2017) Hamolli, L., De Paolis, F., Hafizi, M., & Nucita, A. A. 2017, Astrophysical Bulletin, 72, 73, doi: 10.1134/S1990341317030099
  * Hamolli et al. (2016) Hamolli, L., Hafizi, M., De Paolis, F., & Nucita, A. A. 2016, Ap&SS, 361, 274, doi: 10.1007/s10509-016-2860-7
  * Hamolli et al. (2013) Hamolli, L., Hafizi, M., & Nucita, A. A. 2013, International Journal of Modern Physics D, 22, 1350072, doi: 10.1142/S0218271813500727
  * Henderson & Shvartzvald (2016) Henderson, C. B., & Shvartzvald, Y. 2016, AJ, 152, 96, doi: 10.3847/0004-6256/152/4/96
  * Johnson et al. (2020) Johnson, S. A., Penny, M., Gaudi, B. S., et al. 2020, AJ, 160, 123, doi: 10.3847/1538-3881/aba75b
  * Kim et al. (2021) Kim, H.-W., Hwang, K.-H., Gould, A., et al. 2021, AJ, 162, 15, doi: 10.3847/1538-3881/abfc4a
  * Koshimoto et al. (2023) Koshimoto, N., Sumi, T., Bennett, D. P., et al. 2023, AJ, 166, 107, doi: 10.3847/1538-3881/ace689
  * Marshall et al. (2006) Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S. 2006, A&A, 453, 635, doi: 10.1051/0004-6361:20053842
  * Mislis et al. (2016) Mislis, D., Bachelet, E., Alsubai, K. A., Bramich, D. M., & Parley, N. 2016, MNRAS, 455, 626, doi: 10.1093/mnras/stv2333
  * Mróz et al. (2017) Mróz, P., Udalski, A., Skowron, J., et al. 2017, Nature, 548, 183, doi: 10.1038/nature23276
  * Mróz et al. (2018) Mróz, P., Ryu, Y. H., Skowron, J., et al. 2018, AJ, 155, 121, doi: 10.3847/1538-3881/aaaae9
  * Mróz et al. (2019a) Mróz, P., Udalski, A., Bennett, D. P., et al. 2019a, A&A, 622, A201, doi: 10.1051/0004-6361/201834557
  * Mróz et al. (2019b) Mróz, P., Udalski, A., Skowron, J., et al. 2019b, ApJS, 244, 29, doi: 10.3847/1538-4365/ab426b
  * Mróz et al. (2020) Mróz, P., Poleski, R., Han, C., et al. 2020, AJ, 159, 262, doi: 10.3847/1538-3881/ab8aeb
  * Nikolaus & Hundertmark (2018) Nikolaus, K. S., & Hundertmark, M. 2018, arXiv e-prints, arXiv:1804.10136, doi: 10.48550/arXiv.1804.10136
  * Pacucci et al. (2013) Pacucci, F., Ferrara, A., & D’Onghia, E. 2013, ApJ, 778, L42, doi: 10.1088/2041-8205/778/2/L42
  * Penny (2014) Penny, M. T. 2014, ApJ, 790, 142, doi: 10.1088/0004-637X/790/2/142
  * Penny et al. (2017) Penny, M. T., Rattenbury, N. J., Gaudi, B. S., & Kerins, E. 2017, AJ, 153, 161, doi: 10.3847/1538-3881/aa61a1
  * Penny et al. (2013) Penny, M. T., Kerins, E., Rattenbury, N., et al. 2013, MNRAS, 434, 2, doi: 10.1093/mnras/stt927
  * Poleski & Yee (2019) Poleski, R., & Yee, J. C. 2019, Astronomy and Computing, 26, 35, doi: 10.1016/j.ascom.2018.11.001
  * Rauer et al. (2014) Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38, 249, doi: 10.1007/s10686-014-9383-4
  * Refsdal (1966) Refsdal, S. 1966, MNRAS, 134, 315, doi: 10.1093/mnras/134.3.315
  * Robin et al. (2017) Robin, A. C., Bienaymé, O., Fernández-Trincado, J. G., & Reylé, C. 2017, A&A, 605, A1, doi: 10.1051/0004-6361/201630217
  * Robin et al. (2012a) Robin, A. C., Marshall, D. J., Schultheis, M., & Reylé, C. 2012a, A&A, 538, A106, doi: 10.1051/0004-6361/201116512
  * Robin et al. (2003) Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523, doi: 10.1051/0004-6361:20031117
  * Robin et al. (2012b) Robin, A. C., Luri, X., Reylé, C., et al. 2012b, A&A, 543, A100, doi: 10.1051/0004-6361/201118646
  * Ryu et al. (2021) Ryu, Y.-H., Mróz, P., Gould, A., et al. 2021, AJ, 161, 126, doi: 10.3847/1538-3881/abd55f
  * Scholz et al. (2022) Scholz, A., Muzic, K., Jayawardhana, R., Quinlan, L., & Wurster, J. 2022, PASP, 134, 104401, doi: 10.1088/1538-3873/ac9431
  * Spergel et al. (2013) Spergel, D., Gehrels, N., Breckinridge, J., et al. 2013, arXiv e-prints, arXiv:1305.5422, doi: 10.48550/arXiv.1305.5422
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report.  http://arxiv.org/abs/1503.03757
  * Sumi et al. (2011) Sumi, T., Kamiya, K., Bennett, D. P., et al. 2011, Nature, 473, 349, doi: 10.1038/nature10092
  * Sumi et al. (2023) Sumi, T., Koshimoto, N., Bennett, D. P., et al. 2023, AJ, 166, 108, doi: 10.3847/1538-3881/ace688
  * Wambsganss (1998) Wambsganss, J. 1998, Living Reviews in Relativity, 1, 12, doi: 10.12942/lrr-1998-12
  * Webster & Folta (2017) Webster, C. M., & Folta, D. C. 2017, 1–20
  * Yee (2013) Yee, J. C. 2013, ApJ, 770, L31, doi: 10.1088/2041-8205/770/2/L31
  * Yee et al. (2017) Yee, J. C., Fazio, G. G., Benjamin, R., et al. 2017, arXiv e-prints, arXiv:1710.04194, doi: 10.48550/arXiv.1710.04194