# Classifying High-cadence Microlensing Light Curves I: Defining Features

Somayeh Khakpash Department of Physics, Lehigh University, 16 Memorial Drive
East, Bethlehem, PA 18015, USA  Joshua Pepper Department of Physics, Lehigh
University, 16 Memorial Drive East, Bethlehem, PA 18015, USA  Matthew Penny
Department of Physics & Astronomy, 261-B Nicholson Hall, Tower Dr., Baton
Rouge, LA 70803, USA  B. Scott Gaudi Department of Astronomy, The Ohio State
University, 140 W. 18th Ave., Columbus, OH 43210  R. A. Street  LCOGT, 6740
Cortona Dr, Goleta, CA 93117, USA

###### Abstract

Microlensing is a powerful tool for discovering cold exoplanets, and the The
Roman Space Telescope microlensing survey will discover over 1000 such
planets. Rapid, automated classification of Roman’s microlensing events can be
used to prioritize follow-up observations of the most interesting events.
Machine learning is now often used for classification problems in astronomy,
but the success of such algorithms can rely on the definition of appropriate
features that capture essential elements of the observations that can map to
parameters of interest. In this paper, we introduce tools that we have
developed to capture features in simulated Roman light curves of different
types of microlensing events, and evaluate their effectiveness in classifying
microlensing light curves. These features are quantified as parameters that
can be used to decide the likelihood that a given light curve is due to a
specific type of microlensing event. This method leaves us with a list of
parameters that describe features like the smoothness of the peak, symmetry,
the number of peaks, and width and height of small deviations from the main
peak. This will allow us to quickly analyze a set of microlensing light curves
and later use the resulting parameters as input to machine learning algorithms
to classify the events.

##  1 Introduction

Microlensing is a phenomenon that happens when light emitted from a distant
object (the source) is lensed by a closer, massive object (the lens), and as a
result, multiple images of the source are formed. These images are typically
not resolved because their angular separation is much smaller than the angular
resolution of both ground- and space-based telescopes, and consequently, we
observe a brightening of the source.

Microlensing is a powerful tool for detecting small and dim objects that are
otherwise very hard to detect by their emitted light, and in particular, it is
so far the only method capable of investigating planetary systems with
terrestrial-mass planets orbiting beyond the snow line (Gaudi, 2010). For
several decades, ground-based surveys have searched the Galactic Bulge for
microlensing events. The Optical Gravitational Lensing Experiment (OGLE)
(Udalski et al., 2015), the Microlensing Observations in Astrophysics (MOA)
(Bond et al., 2001) and the Korea Microlensing Telescope Network (KMTNet) (Kim
et al., 2016) surveys detect thousands of microlensing events. These
observations have led to the successful discovery of over 100
exoplanets111https://exoplanetarchive.ipac.caltech.edu and some compact
substellar objects (Griest et al., 1995; Bennett et al., 2002; Mao et al.,
2002).

Identifying microlensing events from the huge numbers of light curves obtained
in these surveys is a challenge. It is useful to be able to distinguish
microlensing events either in real time during survey operations as the events
start, or after an observing campaign has finished and the data of being fully
analyzed. Microlensing surveys have generally relied on circulating real-time
alerts of potential microlensing events when there is an increase in the
brightness of an observed source. This method can lead to alerts for any
sudden rise in a light curve, which can include variability other than
microlensing, such as cataclysmic variables (CVs). There has been recent work
to increase the accuracy of these alerts, limiting them to genuine
microlensing events, using machine learning (ML) algorithms that can
distinguish microlensing light curves from other variability in real time
(Kessler et al., 2019; Godines et al., 2019). Real-time classification can
increase the accuracy of alerts. It can also identify high-value events for
which only partial coverage will be obtained by the survey, and allow
observers to schedule supplemental follow-up observations to increase the
coverage of the event.

There is also a need for classification of microlensing survey data after the
conclusion of observations. At that point, it is necessary to first
distinguish microlensing events from other types of variability, and also to
separate different types of microlensing events, such as those that include
planetary signals. Surveys like the KMTNet survey have employed an automated
approach to detect microlensing-like variability in complete light curves by
fitting simple functions (Kim et al., 2018). Belokurov et al. (2003) also
advocates the use of neural networks to distinguish well-sampled microlensing
light curves from other types of variability. In this paper, we wish to
explore the efficiency and efficacy of after-the-fact classification of
microlensing signals.

ML has now become a popular method for classifying astronomical time series.
Some ML classifiers like the Random Forest (Liaw & Wiener, 2002) and k-mean
(Lloyd, 1982) classifiers take light curve features as input and try to find a
connection between those features and the classification labels. In this
scenario, “features” are quantitative statistical or morphological
measurements of the time series. Some other methods like neural networks use
the light curves themselves as input, and then find common patterns or higher-
order correlated properties to classify them. As a specific example, Random
Forest is a widely-used algorithm that classifies time series by making
decisions based on features in the time series (Bluck et al., 2020; Pawlak,
2019). In order to use this algorithm efficiently, observable features in the
light curves that are most closely related to the canonical model parameters
must first be identified.

Here, we focus on analyzing complete high-cadence microlensing light curves.
These are simulated light curves for the Roman Galactic Bulge Exoplanet
Survey. The Nancy Grace Roman Space Telescope (Roman), formally known as the
Wide Field Infrared Survey Telescope (WFIRST) is a NASA future space mission
that is expected to be launched by mid 2020’s. One of its primary goals is to
detect exoplanets using the microlensing method (Spergel et al., 2015). It is
estimated that Roman will find about $54,000$ microlensing events and will
detect $1400$ planets (Penny et al., 2019). The $0.01$ AU distance of Roman
from the Earth enables measuring parallaxes for free-floating planets and
Earth-mass bound planets to help constrain their masses (Zhu & Gould, 2016;
Street et al., 2018). To assess the value of ML techniques to these surveys we
will apply ML features and classifiers to simulated Roman microlensing data
produced for the Roman microlensing data challenge222https://microlensing-
source.org/data-challenge/ (Street et al. in prep.).

Binary-lens light curves are often difficult to model because of the
complicated features in their light curves and the large parameter space that
needs to be fully searched (for a more thorough review refer to the
introduction of Penny (2014) and Khakpash et al. (2019)). Those challenges
highlight the importance of developing fast automated algorithms to quickly
analyze light curves and estimate the microlensing model parameters. Our
primary goal in this work is to identify a list of features specifically
defined for microlensing light curves that can be generated by fast and
efficient algorithms, and show that these features can help either in
differentiating microlensing light curves among other types of variability or
in classifying microlensing light curves into different types. This would
enable fast detection of planetary system lenses and other interesting lensing
cases in the released datasets of large surveys like Roman.

We present a collection of algorithms including various functional fits that
are applied to the light curves, and from these fits, we extract parameters
that quantify features of the light curve like smoothness of the peak,
symmetry, number of peaks, similarity to microlensing single-lens curves,
number of deviations from the peak, and width and height of the deviations
from the main peak. We then show how effective each of these functions are in
distinguishing between different types of microlensing events like identifying
single-lens versus multiple-lens systems, or planetary system lenses versus
stellar binary lens systems, and in some cases in detecting microlensing
events among other types of variability. A similar work was first done by Mao
& Di Stefano (1994) where they used features such as an estimate of the
asymmetry about the peak to detect binary-lens signatures in the light curves.
In Section 2, we introduce the different models of microlensing events that
need to be parameterized for the classification, and in Section 3, we discuss
the properties of our test dataset. In Section 4, we introduce our algorithm
package including the different functional fits and their respective output
parameters. We also evaluate the effectiveness of each function in capturing
the specific set of features in the light curve. In Section 6, we show
preliminary tests of using these features as input to machine learning
classifiers. Finally, in Section 7, we discuss our results and applicability
of our method to other datasets.

##  2 Microlensing Models

Microlensing occurs when the light coming from a distant source is lensed by a
closer object along the line of sight. As a result, the source appears
brighter as the angular separation of the two objects decreases. This
phenomenon then results in a peak in the light curve of the source star at the
time of closest angular approach. A simple single-lens microlensing light
curve has a single symmetric peak as shown in panel (a) of Figure 1. The
different panels in Figure 1, represent different shapes of microlensing
events that might be present in a large dataset of a galactic bulge survey.
Note that it is practically impossible to include all possible types of light
curve morphologies and this plot is only a small subset of all possibilities.

The light curves in Figure 1 have characteristics representative of the types
of features we are seeking to automatically detect in this work. When there
are multiple lenses or sources, and when there are second-order effects like
finite source effect and parallax, the shape of the single-lens peak will
deviate from a symmetric peak as in panel (a). In this section, we will
discuss the different physical phenomena that can cause these deviations.

Panels (b) and (c) in Figure 1 show two single-lens models affected by the
finite source effect (Witt & Mao, 1994; Nemiroff & Wickramasinghe, 1994; Gould
& Gaucherel, 1996). Panel (c) contains an event caused by a free-floating
planet and therefore the light curve is more affected by the finite source
effect. Panels (d) and (e) include examples light curves of binary-lens events
caused by binary stars, and panel (f) has an example of a planetary binary-
lens microlensing event with a caustic crossing. Panels (g) and (h) have
examples of planetary binary-lens events with no caustic crossings. The
example in panel (g) is an example of a major image perturbation and the event
in panel (h) is a minor image perturbation (for a discussion of major/minor
images refer to Section 2.4.2). In the following subsections, we introduce
different physical models that give rise to the light curve features seen in
Figure 1.

Figure 1: Light curves containing different types of major features in
microlensing light curves. The top row includes light curves created by single
lenses. Panel (a) is a single stellar lens with no significant finite source
effect. Panel (b) shows the effects of the finite source effect on a single-
lens stellar event. Panel (c) shows an event due to a free-floating planet.
Panels (d) and (e) show two examples of events caused by binary star lens
systems. Panels (g) and (h) are due to planetary systems and have no caustic
crossings, whereas Panel (f) is a planetary event that contains a caustic-
crossing event.

###  2.1 Point-source Point-lens Microlensing Light Curves

The Point-source Point-lens (PSPL) model of a lensing event when there is a
single lens and a single source, and the source and the lens can both be
approximated as point-like objects with zero angular size. When there is no
blending of the source with the neighboring stars or light from the lens or
companions to the lens or source, this model can be described by a simple
Paczyński curve as in Equation 1 and 2. $A$ is the magnification of the
source, $t_{0}$ is the time of the maximum magnification, $u_{0}$ is the
impact parameter, and $t_{E}$ is the Einstein crossing time. Figure 2 shows
five different PSPL models with the same Einstein timescale, and different
values of $u_{0}$. In the presence of blending, another parameter is added to
Equation 1 and will yield Equation 3 where $F(t)$ is the differential flux and
$f_{s}$ is the blending parameter and determines what fraction of the total
flux in the aperture is due only to the source. When there are no higher order
effects, these curves remain symmetric, and the sharpness of the peak in high-
magnification events remains intact. Detecting any deviation from this model
can help identify the presence of higher order effects.

$A(t)=\frac{u(t)^{2}+2}{u(t)\sqrt{u^{2}+4}}$ (1)
$u(t)=\sqrt{{u_{0}}^{2}+{\left(\frac{t-t_{0}}{t_{E}}\right)}^{2}}$ (2)
$F(t)=f_{s}\times A(t)+(1-f_{s})$ (3)

Figure 2: Five PSPL models with different values of $u_{0}$. The figure shows
how the shape of the curve changes for different values of the impact
parameter.

###  2.2 Finite Source Effects

When the source is not point-like, and has a finite size, the magnification
function will be different, especially when the impact parameter is comparable
to the source size (Witt & Mao, 1994). In Figure 1, an example of the changes
due to finite source effects can be seen in panel (b). For small impact
parameters, this effect becomes stronger, and it smooths out the peak of the
event. This effect can be used to estimate $\theta_{E}$, the angular size of
the Einstein ring and therefore provide constraints on the lens mass (Yoo et
al., 2004; Jiang et al., 2004).

###  2.3 Free Floating Planets

Isolated, short duration microlensing events can be due to free floating
planets (FFPs) or planets in very wide orbits. Compred to most other exoplanet
detection methods, microlensing has the capability to detect these planets
since the method only depends on the mass of the lens and not its luminosity.
Microlensing events caused by FFPs have short timescales of $t_{E}<2$ days
because the size of the Einstein ring depends on the lens mass, and they can
be detected in high-cadence observations (Sumi et al., 2011; Mróz et al.,
2017, 2018, 2020). Roman’s unique combination of high cadence and small
photometric noise makes it particularly favorable for detecting the free-
floating planets (Johnson et al., 2020).

Microlensing due to FFPs is very likely to be affected by finite source
effects. The finite source effect (Witt & Mao, 1994) parameter, $\rho$, is
defined as ${\theta_{\ast}}/{\theta_{E}}$ where $\theta_{\ast}$ is the
apparent size of the source, and $\theta_{E}$ is the angular size of the
Einstein ring. $\theta_{E}$ depends on the lens mass, and it becomes smaller
for single planetary lenses, and therefore, $\rho$ becomes larger, and the
finite source effect is stronger. As a result, the top of the short-timescale
peak becomes flattened, and the overall shape of the event begins to resemble
a tophat function. An example of such an event can be seen in panel (c) of
Figure 1.

###  2.4 Binary-lens Microlensing Events

When there is more than one lens, the range of possible light curve
morphologies increases dramatically. In binary-lens systems, various shapes of
the light curves depend on parameters such as $q$, the mass ratio, $s$, the
projected separation, and $\alpha$, the angle the path of the source makes
with the line connecting the two lenses. Since $\alpha$ can be any value
between zero and $360^{\circ}$, this quantity along with $q$, $s$, and $t_{E}$
allow the time, height, width, and number of the features in the light curve
vary significantly. Based on OGLE (Udalski, 2004) and MOA (Bond et al., 2001)
observations, around 10% of microlensing light curves show lens binarity
signatures (Penny et al., 2011). Fully modeling each of these light curves is
challenging because there are multiple minima in the $\chi^{2}$ surface as a
function of their parameters, and therefore a wide range of parameter space
must be searched. Furthermore, there is no direct mapping between the
observable features and the canonical parameters
($q,s,\alpha,t_{0},t_{E},u_{0}$).

####  2.4.1 Stellar Binary-lens Events

Stellar binary lenses are likely as common or more common than planetary
system lenses, but the probability of seeing deviations from a PSPL curve is
much higher in nearly equal mass ($q\gtrsim 0.1$) binary lens light curves
than light curves due to planetary mass ratio ($q\lesssim 0.01$) binaries
because the size of the caustics (the loci of infinite magnification for a
point source) increase with increasing $q$. Therefore, they are typically easy
to distinguish from planetary system lenses since the deviations from the
standard PSPL form last much longer, and in particular can be comparable to
$t_{E}$. They usually have larger perturbations in their light curves, but can
be still misinterpreted as planetary binary-lens events (Han et al., 2016).
Two examples of such light curve can be seen in panels (d) and (e) of Figure
1.

####  2.4.2 Planetary Binary-lens Events

A lens system containing a star and a single planet is a binary-lens system
with a very small mass ratio ($q$). Note that there is no strict delineation
(in terms of $q$) between a stellar binary-lens event and a planetary binary-
lens event, and this is one of the reasons that make this a challenging
classification problem. In these cases, the light curve is dominated by the
host star, such that there is a main peak which can be described by the PSPL
model with one or more small deviations caused by the planet (Gaudi, 2012). As
in a single-lens microlensing event, two images of the source will be formed
as it passes close by the lens star; an image outside of the Einstein ring of
the lens star called the major image, and an image inside the Einstein ring of
the lens star called the minor image. Depending on whether or not the
projected separation of the planet and the lens star ($s$) is larger or
smaller than the Einstein ring, the planet will have to perturb one of these
images to leave a signature on the star light curve, and the perturbation will
have different characteristics. The effect of the planet on the light curve
can also be explained by the positions of caustics. In real cases, since the
source star is not a point source, we observe sharp high-magnification peaks
when there are caustic crossings. The numbers and sizes of these caustics
depend on the mass ratio ($q$) and the projected separation ($s$). The
structure of the caustics leads to the planetary features in the light curves
of these events (Gaudi, 2012).

Caustics have three main topologies. Depending on the values of the projected
separation and the mass ratio, there is either one (intermediate topology),
two (wide topology), or three (close topology) caustics (for a thorough
description of planetary caustics please refer to Gaudi (2012)). These
caustics are also classified into two classes, called the central caustic and
the planetary caustic, and depending on which one is closer to the path of the
source, the shapes of the light curve are different.

If the source crosses a caustic curve, two sharp deviations will be seen in
the light curve assuming sufficient sampling corresponding to the entrance and
exit of the source through the caustic. The shape of these deviations depends
on the path of the source which is determined by the source path angle
($\alpha$). Panel (f) of Figure 1 includes an example of a caustic-crossing
planetary microlensing event. In this work, we analyze these perturbations
based on whether or not the source has crossed the caustics, and we
investigate the features created in each of these cases. In cases with no
caustic crossing, the source only passes close to the caustics, and that
results in features like short bumps or troughs, or both, in the light curve.
These deviations are typically small in magnitude and short in duration which
makes them sometimes very difficult to detect, especially in low-cadence
surveys, and they can be created by either the central or the planetary
caustics. In panels (g) and (h) of Figure 1, we show examples of non-caustic-
crossing events caused by the planetary caustic. These perturbations happen
when the source passes closely by the planetary caustic, and finding the time
of the perturbation can help determine the location of the caustic, and
therefore estimate the projected separation of the planetary system. The event
in panel (g) is a perturbation in the major image and the event in panel (h)
is a perturbation in the minor image.

In cases with caustic crossing, we need to identify the times of the two sharp
deviations to find an estimate of the time the source spends inside the
caustic. This allows us to estimate the size and location of the caustics, and
therefore find an estimation of the mass ratio and the projected separation of
the planetary system (e.g., Poleski et al., 2014).

##  3 Scope of the Analysis

As mentioned in Section 1, classifying astronomical time series using ML
algorithms relies on the data sets used for training and testing. Features
computed for a data set cannot always be used to analyze another data set;
therefore, the connection between the dataset and the defined features should
be discussed. In this paper, we focus on computing features for simulated
Roman light curves. These light curves are simulated based on the Roman Cycle
7 mission design and have six 72-days seasons with 15 minute cadence. These
light curves are similar to the light curves generated by Penny et al. (2013)
and Penny et al. (2019), and the full details of the simulated data can be
found in these papers. This data set was designed in particular to be used for
the Roman microlensing data challenge333https://microlensing-source.org/data-
challenge/ and has four classes of variability types, consisting of different
types of microlensing light curves along with a particular type of stellar
variability that can be misinterpreted as microlensing. These classes include
cataclysmic variables, single-lens systems, stellar binary lenses and
planetary system lenses. An example of each class of light curves used in this
analysis can be seen in Figure 3.

Figure 3: Four classes of light curves in the dataset used for this analysis.
From left: is an example of the light curve of a cataclysmic variable. The
second panel shows a single-lens model which contains a single symmetric peak.
The third panel represents lens systems containing a stellar binary, and the
last panel is an example of an event with a planetary system as the lens. The
blue data points are from the original simulated light curves, and the black
data points are smoothed light curves using the low-pass filter.

The dataset contains $4181$ light curves, including $429$ cataclysmic
variables, $1626$ single-lens systems, $1386$ stellar binary lenses, and $688$
planetary system lenses. The CV class contains CV variabilities with either
zero, one or multiple instances of outbursts in the baseline. The single-lens
systems include both stars and planets as lenses, therefore their shape varies
in terms of the duration, amplitude, and morphology of the peak. As mentioned
before, the stellar binaries and the planetary systems can result in very
similar light curves, and the transition between these two classes is smooth
in terms of the observables, and distinguishing between them is one of our
main concerns.

The lightcurves selected for this project have a wide range of signal to noise
ratio. To enter our sample, events had to pass a $\Delta{\chi}^{2}>500$ cut of
the true model to a flat line, and binary and planetary events had to pass a
${\Delta\chi}^{2}>160$ cut of the true model relative to a single lens
lightcurve model (see Penny et al. (2019) for full details). While these
${\Delta\chi}^{2}$ values suggest quite high formal detection significance,
the cuts can be passed by lightcurves with signal to noise per data point less
than unity, e.g., a $t_{\rm E}=25$ day single lens event will have 4800 data
points within one $t_{\rm E}$ of the event peak. Similarly, a Jupiter-mass
planet with a planetary deviation lasting ${\sim}1$ day (i.e., ${\sim}96$ data
points) can pass the detection cut with an average contribution of
${\Delta\chi}^{2}$ per data point of less than 2. This said, our selection
cuts do not allow us to test our method’s sensitivity to the lowest signal to
noise ratio events that could potentially be formally detected.

##  4 Feature Detection and Parametrization

In Section 2, we described major types of microlensing variability with
different morphologies and amplitudes. In this section, we introduce a package
of statistical tests and model fits we have assembled into an algorithmic
approach to detect the features of a set of light curves. These features
individually can then be used to distinguish microlensing events from other
variability, and to classify the different types of microlensing systems.
Ultimately, in a subsequent paper, we will use the combined power of all of
the features as input to a ML classifier. The package includes the following
tools:

  * •

Peak Finder

  * •

Gould Two-parameter Point-source Point-lens Fit (G-PSPL)

  * •

Symmetry Check

  * •

Trapezoidal Function Fit

  * •

Cauchy Distribution Fit

  * •

Planetary Parameters Finder

  * •

Chebyshev Polynomials Fit

In the sections below we introduce each of these tools, and explain how they
are applied to the simulated data set and evaluated. Not all tools are used
independently, and some are used in combination with other tools.

Before applying some of these tools, we employ a low-pass filter to the light
curve which allows low frequency signals to pass and therefore reduces the
high-frequency noise in the data. This low-pass filter has a smoothing window
of 10 subsequent data points. The width of the smoothing window was optimized
for the simulated Roman data set analyzed in this work, but should be
separately optimized for each data set so as to remove non-astrophysical noise
as much as possible.

###  4.1 Peak Finder

The primary idea behind this tool is to identify individual peaks in a light
curve, which are often characteristic of microlensing events. After employing
the low-pass filter, we bin the light curve in time, and then, in each bin, we
count the number of data points with positive deviations larger than some
number of standard deviations in that bin, which we define as the peak
threshold. If a bin has more than one data point beyond that threshold, that
bin is defined as a peak.

When using this algorithm, two parameters are important in determining the
success rate; the bin size and the peak threshold. Depending on the particular
implementation, we can use this method to search for a single large peak, such
as a single-lens microlensing event. Or after identifying a single-lens event,
we can subtract a model and then search for additional peaks in the light
curve residuals that would indicate a binary lens event. We can also search
simultaneously for multiple large peaks at once.

To examine the performance of this method in a particular case, we consider
the case of searching for a single peak within a light curve. We test
different values of the bin size and threshold parameters to determine what
combination is more successful at detecting peaks. We define a successful
detection as one where exactly one bin identifies a peak, and the light curve
is known to be due to a microlensing event. We define a failed detection as
one where a light curve labeled as containing a microlensing event did not
have any peaks identified, or multiple peaks identified. With this range of
bin sizes, this algorithm will not detect narrow deviations in microlensing
light curves, and will only look for the main event. Note that the algorithm
is looking for single-peaked variations, so it might also find single-peaked
transients. We expect that next stages of analysis will distinguish between
these cases. This method has a low false negative rate. For example, it is
possible that it would fail at detecting a microlensing light curve at low
thresholds. This happens when the event has a low magnification and it might
be considered similar to the noise in the light curve.

Figures 4, 5, and 6 show the performance of the algorithm across a range of
bin sizes and thresholds. They display the True Positive Rate (TPR), False
Positive Rate (FPR), and False Discovery Rate (FDR) for different options of
bin size and peak threshold. Figure 4 shows that bin sizes larger than about
50 days and a thresholds of $3\sigma$ to $5\sigma$ give the highest TPR of
about 80%.

Figure 5 shows the False Positive Rate which represents what fraction of non-
microlensing events are labeled as microlensing. Bin sizes do not consistently
affect the FPR. Thresholds larger than $5\sigma$ give unreasonably low FPR,
and that is because at these thresholds, most light curves will have zero
peaks identified. Thresholds of $3\sigma$ to $5\sigma$ give a FPR of 50%.
Since some fraction of non-microlensing events have a single peak, we expect
that a group of non-microlensing events will be classified as microlensing
even though the algorithm has successfully identified one peak in them.

Figure 6 show the False Discovery Rate which represents the fraction of events
labeled as microlensing that are are in fact non-microlensing. For the same
reason as mentioned above, thresholds larger than $5\sigma$ are not favorable,
and thresholds of $3\sigma$ to $5\sigma$ have reasonably low FDR of larger
than 7% for most bin sizes. Note that the extremely low rates of FDR is also
an indication of the unbalanced dataset (a dataset in which there is unequal
number of objects in each category.), and our purpose is to compare how
sensitive it is to different thresholds and bin sizes.

Using this peak finder algorithm with a bin size of $60$ days and a threshold
of $5\sigma$, we can identify which light curves show a single peak event. In
Section 5, we use this feature in conjunction with the other feature detection
tools.

Figure 4: Results of the Peak Finder algorithm is applied to all of the light
curves in the dataset for different values of bin size and threshold. The True
Positive Rate (TPR), the fraction of microlensing events that are labeled as
microlensing is plotted versus bin size, with various thresholds denoted by
symbol and color. The value of TPR is roughly independent of the bin sizes,
and that thresholds of $3\sigma$ to $5\sigma$ yield the highest TPR.

Figure 5:  As Figure 4, but showing the False Positive Rate (FPR), the
fraction of non-microlensing events that are labeled as microlensing. The
value of FPR is only weakly dependent of the bin sizes for most thresholds.

Figure 6:  As figure 4, but showing the False Discovery Rate (FDR), the
fraction of events labeled as microlensing that are in fact non-microlensing.
Overall, values of FDR are very low and also roughly independent of the bin
sizes for most thresholds.

###  4.2 Gould Two-parameter PSPL Fit

The PSPL function has been demonstrated to be a good fit to most microlensing
light curves, and to not match most other types of astrophysical variability
(Di Stefano & Perna, 1997). Therefore, the goodness of fit to the PSPL model
can be used as a feature of the light curve. In this section, we demonstrate
employing a two-parameter PSPL model first introduced by Gould (1996) to
detect curves similar to a Paczyński curve. The reason for not using the
regular PSPL model is that brute force searches over the three non-linear
parameters ($t_{E}$, $t_{0}$, $u_{0}$) of a PSPL model are slow, while the
two-parameter Gould approximate PSPL model (G-PSPL) captures the basic
observables of a full PSPL model but requires many fewer computations. For
more details refer to Gould (1996) and Woźniak & Paczyński (1997).

In order to employ the PSPL model as a light curve feature, we employ a
version of the Gould (1996) two-parameter PSPL model (G-PSPL) that has been
used by the KMTNet survey to detect microlensing events. For that survey, Kim
et al. (2018) fit a combination of two two-parameter G-PSPL models: one with
the assumption of $u_{0}=1$, and one with the assumption of $u_{0}$ being very
close to zero. This approach reduces the number of fitted parameters, and can
be used to detect both high and low-magnification events. This double two-
parameter G-PSPL fit is defined here as function $F(Q)$:

$F(Q)=f_{1}\times(A1(Q)+A_{2}(Q))+f_{0}$ (4)

where

$A_{1}(Q)=\frac{1}{\sqrt{Q}}$ (5)

and

$A_{2}(Q)=\frac{1}{\sqrt{1-\frac{1}{(1+\frac{Q}{2})^{2}}}}$ (6)

are functions that represent good fits to high and low magnification events,
respectively. The two parameters $f_{0}$ and $f_{1}$ do not have a physical
meaning in this equation, but in the limit of $u_{0}=1$ and $u_{0}$ close to
zero, they are related to the blending and source parameters, described in
detail in Kim et al. (2018). $A_{1}$ and $A_{2}$ are both functions of $Q$,

$Q(t)=1+\left(\frac{t-t_{0}}{t_{eff}}\right)^{2}.$ (7)

$Q$ is a function of time and depends on two parameters $t_{0}$ and $t_{\rm
eff}$. $t_{0}$ is the time of the maximum magnification, and $t_{\rm eff}$
approximates $u_{0}\times t_{E}$ in the limit of $u0\ll 0.5$, and helps
characterize the amplitude and duration of the microlensing event. For fitting
the function $F(Q)$, we set the initial values of the parameters $f_{0}$ and
$f_{1}$ to $0.5$, and $t_{0}$ is the time of the maximum magnification in the
light curve. For $t_{\rm eff}$, we choose a list of seven initial values
$0.01,0.1,0.5,2,5,10,$ and $20$. For high magnification events, low values of
$t_{\rm eff}$ usually provide a better fit, and for low magnification events,
larger values of $t_{\rm eff}$ are a better fit. The fit with the minimum
${\chi}^{2}$ value is used as the selected model. Note that for the remainder
of the paper, we specify whether we use A G-PSPL or a PSPL model fit.

###  4.3 Symmetry Check Algorithm

If we identify a peak in a light curve and fit a two-parameter G-PSPL model to
the data, we then employ a check on the symmetry of the event. In Section 4.2,
we showed that this method can be used to detect microlensing events. Here, we
calculate the reduced ${\chi}^{2}$ statistic separately for the right and left
side of the peak in the light curve. If the peak is symmetric, the ratio of
${{({\chi}^{2}_{left})}_{reduced}}/{{({\chi}^{2}_{right})}_{reduced}}$, which
we define as $\beta$, is very close to unity. Since the number of data points
in the right and left wings of the peak might not be exactly equal, we use the
reduced ${\chi}^{2}$. Deviations of $\beta$ from unity can be related to
additional physical details in the light curves, such as the existence of
planetary or binary star companion lensing signatures, parallax effects, or
the asymmetry in non-microlensing peaks like cataclysmic variables.

In Figure 7 we plot histograms showing the distributions of variability types
for a set of ranges of $\beta$ values. Each column represents events with a
specific range of $\beta$. We find that for CV light curves, the asymmetry is
extremely large, consistent with the morphology of CV light curves, displaying
steep rises and more gradual decreases. Single-lens microlensing events are
symmetric, with $\beta$ values close to unity. Binary lenses have a wide range
of $\beta$ values depending on whether they are caused by a planetary system
or a stellar binary system. These ratios can be used as a feature to obtain
information about the kinds of deviations from the PSPL model, and classify
the light curves.

Figure 7: $\beta$ is a measurement of asymmetry of the peaks in the light
curves that is obtained from the symmetry check algorithm. This plot shows the
distribution of $\beta$ for different classes of light curves in our dataset.

###  4.4 Trapezoidal Function Fit

As explained in Section 2.3, lensing events caused by free-floating planets
have a short duration and amplitude and are more often strongly affected by
the finite source effect. Therefore, their shapes approximately resemble a
trapezoidal function.

The parameters of the trapezoidal function are the baseline magnitude ($a$),
maximum magnitude ($b$), time of the first rise (${\tau}_{1}$), the time when
maximum magnitude is reached (${\tau}_{2}$), the time when maximum magnitude
ends (${\tau}_{3}$), and the time when the magnitude returns to the baseline
(${\tau}_{4}$). Figure 8 shows a diagram of the function and its six
parameters. We must first find initial guesses for the six parameters, then
fit the function and obtain more accurate estimates of the parameters. We also
calculate the duration of the trapezoidal portion $\Delta\tau_{full}$, and the
duration of the flat section of the trapezoid, $\Delta\tau_{top}$ as presented
in Equations 8 and 9 and shown on Figure 8.

$\Delta\tau_{full}={\tau}_{4}-{\tau}_{1}$ (8)

$\Delta\tau_{top}={\tau}_{2}-{\tau}_{3}$ (9)

Figure 8: A trapezoidal function with its defined parameters.

In order to find initial guesses for the six parameters, we first subtract the
median baseline magnitude from the light curve, setting $a$ to zero. We then
select $20$, $2$, $0.2$, or $0.02$ days as initial guesses for the total
duration ($\Delta\tau_{full}$). We define $t_{0}$ as the time of maximum
brightness, and assume $\Delta\tau_{full}=2\Delta\tau_{top}$. We then define
the remaining parameters ${\tau}_{1}$, ${\tau}_{2}$, ${\tau}_{3}$, and
${\tau}_{4}$ as the following quantities, respectively:
$t_{0}-\frac{\Delta\tau_{full}}{2}$, $t_{0}-\frac{\Delta\tau_{full}}{4}$,
$t_{0}+\frac{\Delta\tau_{full}}{4}$, and $t_{0}+\frac{\Delta\tau_{full}}{2}$.
We fit the resulting models based on the initial guesses for
$\Delta\tau_{full}$, using a least-squared minimization approach. We finally
select the fit that minimizes the ${\chi}^{2}$ statistic.

We have applied this algorithm to our test dataset, with Figure 9 showing an
example of a best fit function for a simulated light curve. We find that half
of the total duration ($\Delta\tau_{full}$) is a good representation of the
Einstein timescale ($t_{E}$) of the event, and we compare it with true values
of $t_{E}$ of the events. Figure 10 shows the estimated versus true duration
of single-lens microlensing events with $t_{E}<2$ days that are candidates for
free-floating planets. We find that for events with $t_{E}<2$ days, there is a
median absolute deviation of $0.06$ days, while for events with $t_{E}\geq 2$
days, there is a median absolute deviation of $0.8$ days.

Figure 9: Trapezoidal function fitted to the light curve of a simulated event
caused by a free-floating planet. The flat top of the trapezoid describes the
duration of maximum brightness fairly well, and the duration of the event is
well represented by the full trapezoidal duration.

Figure 10: True versus estimated Einstein Timescale ($t_{E}$) of the events
with $t_{E}<1$ days found by the trapezoidal function fit. These events are
candidates for free-floating planets.

We define another parameter
$\kappa=\frac{\Delta\tau_{top}}{\Delta\tau_{full}}$ as the trapezoidal
timescale ratio, which is the ratio of the duration of the flat part of the
trapezoid to the total duration of the trapezoid. This quantity was initially
set to 0.5 for the initial guesses of the function fitting. This parameter can
be thought of as a measure of how square or peaked the event is, analogous to
the kurtosis of the event. We show that using $\kappa$ along with
$\Delta\tau_{full}$ allows us to flag the events that are strongly affected by
the finite source effect. We can characterize the approximate significance of
the finite source effect using $|u_{0}|/\rho$, the ratio of the impact
parameter to the finite source ratio. When $|u_{0}|/\rho$ is smaller than
unity, it implies the lens passes directly over the source and thus that can
be a sign of significant finite source effects.

The left panel of Figure 11 shows the trapezoidal timescale ratio ($\kappa$)
versus $\Delta\tau_{full}$ for the subset of all light curves that meet the
fitting criteria for single-lens models, which are described below in §5. Note
that both of the quantities on the axes in the left panel are found
experimentally using the trapezoidal function fit. In that panel, most of the
green data points have larger values of $\kappa$ which implies a more square
light curve shape caused by strong finite source effects. The right panel
shows the distributions of the true values of $\rho$ for all the data points
inside and outside of the black box in the left panel. That panel indicates
that light curves with small values of $\Delta\tau_{full}$ and large values of
$\kappa$ tend to have larger values of $\rho$ (grey), and light curves with
larger values of $\Delta\tau_{full}$ and smaller values of $\kappa$ tend to
have smaller values of $\rho$ (hatched white). This is because events with
short duration caused by free-floating planets are more affected by the finite
source effect, and therefore more resemble a trapezoidal functions, with
$\Delta\tau_{full}$ and $\kappa$ helping to indicate these events.

Figure 11:  Left Panel: Values of the trapezoidal timescale ratio ($\kappa$)
are plotted versus the total trapezoidal duration ($\Delta\tau_{full}$). The
purple data points show events with $|u_{0}|/\rho>1$, and the green data
points represent events with $|u_{0}|/\rho<1$, which indicates events with
free-floating planets or other PSPL events affected by the finite source
effect. $\kappa$ is larger for shorter events, and most of the events with
$|u_{0}|/\rho<1$ have large $\kappa$ and small $\Delta\tau_{full}$. Therefore,
the selected ranges of $\kappa$ and $\Delta\tau_{full}$ can be used to
identify likely events caused by free-floating planets. Right Panel: The
distribution of $|u_{0}|/\rho$ for all data points inside (grey) and outside
(hatched white) of the black box in the left panel. Light curves with small
values of $\Delta\tau_{full}$ and large values of $\kappa$ tend to have
smaller values of $|u_{0}|/\rho$, and light curves with larger values of
$\Delta\tau_{full}$ and smaller values of $\kappa$ tend to have larger values
of $|u_{0}|/\rho$.

###  4.5 Cauchy Distribution Fit

For microlensing events with a strong finite source effect in the absence of
limb darkening effects, the top of the peak becomes flatter and will not look
like a PSPL function. For detecting this phenomenon, we need to parameterize
the flatness of the top of the peak. Effectively, we would like to use a
functional fit that can apply to the range of morphologies between that of a
single-lens PSPL model, and a trapezoid that much more closely fits a strong
finite source effect. Our goal here is to find a correlation between the
fitted parameters and the finite source ratios. For that purpose, we use a
Cauchy distribution function to fit the light curves. The difference in the
minimum ${\chi}^{2}$ of the Cauchy fit and the PSPL fit (shown in Equations 3
and 1), along with one of the parameters of the Cauchy distribution, can be
used as features.

The Cauchy distribution shown in Equation 10 has four parameters, time of the
peak ($t_{0}$), the duration of the peak ($\sigma$), the flatness of the peak
($b$), and the amplitude of the peak ($a$).

$C(t)=\frac{a}{1+{\left|\frac{t-t_{0}}{\sigma}\right|}^{2b}}$ (10)

The functions $C(t)$ and PSPL have two common parameters, time and duration of
the peak. After fitting these two functions to a light curve, we have fitted
values for $t_{E}$ and $t_{0}$, which we refer to as $t_{E,PSPL}$,
$t_{E,Cauchy}$, $t_{0,PSPL}$, and $t_{0,Cauchy}$. Note that $t_{E,Cauchy}$ is
equivalent to $\sigma$ in Equation 10. We define the difference in the
${\chi}^{2}$ of the PSPL fit and the Cauchy fit as a feature characterizing
the flatness of the top of the curve represented in Equation 11. In order to
calculate $\psi$, we select a section of the curve at the peak, between
$t_{0}-(t_{E,PSPL}\times u_{0,PSPL})$ and $t_{0}+(t_{E,PSPL}\times
u_{0,PSPL})$.

$\psi={{\chi}^{2}}_{PSPL}-{{\chi}^{2}}_{Cauchy}$ (11)

The more the event is affected by the finite source effect, the more the top
of its peak deviates from the PSPL model and is more similar to the Cauchy
model. Thus for single-lens events with a flatter peak, $\psi$ is positive.
Even in cases where the Cauchy distribution fits much better than the PSPL
fit, for events that are not caused by free-floating planets that still
experience strong finite source effects, the Cauchy distribution will
typically be a better fit than the trapezoidal function. Figure 12 shows a
single-lens event highly affected by the finite source effect, along with the
best-fit PSPL and Cauchy models. The two curves deviate the most close to the
peak of the event.

Figure 12: Cauchy and PSPL functions fitted to a Roman simulated single-lens
microlensing light curve highly affected by the finite source effect. The blue
data points are from the simulated light curve, the red curve is the fitted
PSPL model, and the orange curve is the fitted Cauchy model.

Figure 13 shows values of the Cauchy feature ($\psi$) versus $\rho/u_{0}$ for
all single-lens events that have at least four observations within a time of
$t_{E,PSPL}\times u_{0,PSPL}$ from the maximum. The upper panel shows positive
values of $\psi$ and the lower panel shows the negative values. The events
with flatter peaks have a positive value of $\psi$, and on average, have a
shorter duration than those with negative $\psi$. Events with large positive
$\psi$ (flatter peaks) appear in the upper right of this figure and have
$\rho/u_{0}>1$. Events with a negative $\psi$ include both short and long
duration events. The black lines show the one-to-one relation for events with
positive and negative $\psi$, and have a median absolute deviations of 0.45
for positive $\psi$, and 0.64 for negative $\psi$. Events in the lower region
mostly have small negative values of $\psi$, indicating that neither the PSPL
and Cauchy functions are significantly better fits to the data. Those with
large negative values are those with sharp peaks that are well-described by
the PSPL function, which also have low $u_{0}$. This plot shows that the sign
and magnitude of $\psi$ can help identify the flattest and sharpest
microlensing peaks, and flag the ones that are more likely affected by the
finite source effect.

Figure 13:  The Cauchy feature ($\psi$) plotted versus $\rho/u_{0}$,
restricted to cases with at least four observations near the peak maximum. The
colors of the points represent the PSPL-estimated event duration $t_{E,PSPL}$,
while the sizes of the points represent the PSPL-estimated value of $u_{0}$
($u_{0,PSPL}$). The black trend lines show a one-to-one relation between
$\psi$ and $\rho/u_{0}$ (reversed in the bottom panel).

The other feature that we use for this purpose is parameter $b$ in Equation
10. We expect that if the event exhibits strong finite source effects which
means a flatter top, this feature would take a larger value, and if the finite
source effect is negligible, the feature would be very close to unity (in the
absence of limb darkening effects). In reality, $b$ does not correlate with
$\rho$ for most events, but once $\rho$ is greater than 0.1, we do see a
positive correlation with $b$, as seen in Figure 14. We therefore retain $b$
as a potentially useful feature.

Figure 14: Plot of the parameter $b$ in equation 10 versus true parameters of
$\rho$. It shows that the $b$ values are in general larger for events with
$\rho$ larger than $0.1$.

###  4.6 Planetary Parameters Finder

Our goal in this section is to measure the two physical parameters associated
with the planetary binary-lens events. These are $s$, the projected separation
between the planet and the lens star and $q$, the mass ratio of the planet and
the lens star. The approach that we take here is to first fit the main event
with a PSPL function, then find the deviations from the PSPL model by looking
at the residual, to identify peaks or troughs. We then seek to identify cases
where these residuals have one significant peak or trough, or if they have two
significant peaks.

The reason we take this approach is that the single-deviation residuals can be
characterized much more robustly than the double-peaked ones. For events of
the both groups, we fit the “busy” function introduced by Westmeier et al.
(2014). This function is primarily used to describe double-peaked features in
spectra, and with some changes in its parameters, we can also use it to find
the single-peaked deviations. Khakpash et al. (2019) suggest fitting the
single-peaked residuals with a Gaussian, and show that this approach is useful
mostly for low-mass ratio events. Here, we take the same approach to find the
initial parameters in the residual, but instead of a Gaussian, we fit all
deviations with the busy function. After fitting the residual, in the next
round of fittings, we fit a PSPL plus a busy function, and we use parameters
obtained from the first two fits as initial parameters of this fit. Next, we
calculate $s$ and $q$ using the parameters found by the final fit.

The busy function is commonly used to describe double-horn profile of galaxy
spectra (Westmeier et al., 2014). The function is shown in Equation 12, and
has nine parameters that are sketched in Figure 15. This function comprises
two error functions and a polynomial of degree $n$. The parameters $x_{e}$ and
$x_{p}$ determine location of the middle of the two error functions and the
middle of the polynomial, ${\delta}_{1}$ and ${\delta}_{2}$ are the steepness
of the two error functions, $n$ is the degree of the polynomial and determines
the steepness of the middle trough, $c$ determines the depth of the
polynomial, $w$ determines the distance of the error function zero points from
$x_{e}$, $a$ determines the height of the error functions, and $\varepsilon$
is the horizontal scaling of the function.

$\begin{split}A(t)\;=\;(a/4)\;\times\;\;({\rm
erf}\;(\;{\delta}_{1}(\;w+\varepsilon\times t-x_{e}))+1)\;\\\ \times\;\;({\rm
erf}\;(\;{\delta}_{2}(\;w-\varepsilon\times t+x_{e}))+1)\;\\\
\times\;\;(c\;\times\;{|\varepsilon\times t-x_{p}|}^{n}+1)\end{split}$ (12)

Figure 15: A schematic plot of the busy function, with features associated
with eight function parameters. The parameter $\varepsilon$ (not shown on the
figure) represents a horizontal scaling of the function without disrupting its
shape.

Depending on the values of the nine parameters, the shape of the function can
differ significantly. We take advantage of this fact and use different shapes
of the function to fit single-peaked and double-peaked deviations from the
PSPL model. The two forms of the busy function used to fit the residuals are
shown in Figure 16. The left panel shows a double-horn shape of the busy
function, and the values of the parameters can be found on the plot. The three
curves represent the shapes of the function when only the steepness of the two
error functions are altered. This set of parameters results in a form that can
describe the caustic-crossing features of the planetary microlensing light
curves. In the right panel, the value of $c$ in Equation 12 is set to zero,
and therefore the equation only comprises two error functions. The shape of
the function resembles the Gaussian function used by Khakpash et al. (2019) to
fit single-peaked deviations with an additional ability to describe asymmetric
deviations. The three curves show how the shape of the function changes as the
steepness of the two error functions are altered.

Figure 16: Two different shapes of the busy function that are fitted to the
single- and double-peaked residuals. Left Panel: A double-horn shape of the
busy function, and the values of the parameters. The three curves represent
the shapes of the function when only the steepness of the two error functions
are altered. Right Panel: The value of $c$ in Equation 12 is set to zero, and
the equation only comprises the two error functions and no polynomial. Only
one peak is present, with the two error functions controlling the steepness of
the two sides of the peak. The three curves show how the shape of the function
changes as the steepness of the two error functions are altered.

We apply this approach to the set of simulated Roman planetary microlensing
light curves. We follow the procedure in Khakpash et al. (2019), and first fit
a PSPL function as shown in Equations 1 and 3 to the light curves. We find an
initial value for $t_{0}$ by finding the time of the maximum flux in the light
curve, and we set the initial value of $f_{s}$ to be $0.5$. For an initial
guess for the duration of the event $t_{E}$, we first interpolate between the
data points using a cubic interpolation to estimate a continuous version of
the data, then, for events with $u_{0,initial}<0.5$, we assume $t_{E}$ is the
interval between times when the magnification is $1.34$. For events with
$u_{0,initial}>0.5$, we set $t_{E}$ equal to the intervals between times when
magnification is $1.06$. We use Equations 1 and 3 to calculate an initial
guess for $u_{0}$ based on the above other assumed values at time $t=t_{0}$.
Then, we fit the PSPL function using the initial guesses, and then apply the
peak-finding algorithm to the residual.

Next, we use the peak finder (§4.1) to look for deviations in the residual,
and we force it to find one or two peaks. We start with a large bin size of
about half the lightcurve baseline and a threshold of $3\sigma$. With these
values, the algorithm finds the most significant deviation. Then, we start
decreasing the bin size by 50% until it either finds two peaks or the bin size
reaches 0.2 days. At any of these decreasing steps, if it finds two peaks that
are separated by less than 10 days, we end the search. This condition aims at
preventing the algorithm from selecting multiple peaks in a caustic-crossing
microlensing event. If the search concludes with the detection of one
deviation, we accept that and move on to the next step. If it finds zero
deviations, we simply select the maximum or minimum data point in the residual
and identify that as a single deviation.

At this point, we fit the modified busy function to peaks identified in the
residual. For this purpose, we identify the time of the single deviation, or
the middle point of the double-peaked deviation to be initial guesses for
$x_{e}$ and $x_{p}$. We set these two values equal to zero, the values of
${\delta}_{1}$ and ${\delta}_{2}$ are set equal to $1$, and the value of
$\varepsilon$ is set to $5$. Parameters $a$, $w$, $n$ and $c$ are set up
differently depending on the number of deviations found in the residual. If
there is one deviation, $n$ and $c$ are set equal to zero, and $a$ is equal to
the largest peak or trough in the residual, and $w$ is set equal to unity. If
there are two deviations, initial values of $n$ and $c$ are set to $10$ and
$0.02$, and the value of $a$ is set equal to the median of the residual values
between the two peaks, and $w$ is set equal to half of the distance between
the two peaks.

Parameters describing the deviations in the residual are obtained by the busy
function fit. Using these parameters along with the PSPL parameters as initial
values, we then fit a PSPL plus a busy function to the light curve, and we
find final values of $13$ parameters, four of them describing the main event,
and nine of them describing the deviation. Using these values, we calculate
the two physical parameters $s$ and $q$ using two different approaches.

When only one deviation is found, the busy function will turn into two error
functions with six parameters. After fitting the function to the deviation, we
determine the duration of the deviation by looking at where the fitted model
is not zero, and we assume half of that to be the duration $t_{Ep}$. $t_{p}$
is set to $x_{e}$, and then we find the two values of $s$ by solving Equation
13. If the deviation is positive, a secondary check is done to find if there
is a large trough close to this peak. If there is a large de-magnification
close to a single caustic crossing, the event is very likely to have $s<1$
which could also be due to a resonant caustic. In the secondary check, ratio
of the sum of negative data points in the residual of the busy function fit to
the sum of positive data points is calculated, and if the ratio is larger than
unity, the value of $s<1$ is chosen, otherwise $s>1$ is selected. If the
deviation is negative, the value of $s<1$ is chosen. The estimated value of
$q$ is then obtained by Equation 14.

$|s-\frac{1}{s}|=u,\;{\rm
where}\;u=\sqrt{{u_{0}}^{2}+{\left(\frac{t_{0}-t_{p}}{t_{E}}\right)}^{2}}$
(13)

${q}={\left(\frac{t_{Ep}}{t_{E}}\right)}^{2}$ (14)

When there are two deviations, we assume that the deviations are caused by
crossing the planetary caustic, and that these epochs correspond to when the
source approaches two cusps of the caustic. Using these epochs, we calculate
the distance of the source from the lens at these times using
$u_{1,2}=\sqrt{{u_{0}}^{2}+{\left(\frac{t_{0}-t_{1,2}}{t_{E}}\right)}^{2}}$.
Han (2006) calculates the size of the planetary caustics in terms of $s$ and
$q$, and we find the size of the caustic using simple geometry, and then, we
use their formulation to calculate $s$ and $q$.

Figure 17 shows an example of the geometry involved in these lensing events.
In this figure, planetary caustics of two systems with projected separation of
$1.6$ and $0.8$ Einstein radius, and mass ratio of $0.03$ is shown. The shape
and sizes of the caustics depend on $s$ and $q$. An example path of the source
through the system is shown on both panels. The distance between the lens and
the center of the caustics, $LX$, gives us the value of $|s-\frac{1}{s}|$. At
this point, we do another secondary check as described earlier and we
determine whether we have a system with $s>1$ or $s<1$. Next, for systems with
$s>1$, we find the vertical dimensions of the caustic (left panel of Figure
17), and for systems with $s<1$, we find the distance between the two caustics
by geometry (right panel of Figure 17). We refer to both of these values as
height of the caustics. Han (2006) finds that these values are roughly equal
to $\frac{2\sqrt{q}}{s\sqrt{s^{2}-1}}$ and
$\frac{s^{3}\sqrt{q}}{4}\sqrt{(\frac{2}{s})^{8}+27}$, respectively. We then
use these equations to find $q$. Table 1 summarizes the process we adopt to
extract the system parameters from the light curve residuals in the cases of
one and two deviations (Poleski et al., 2014; Bozza, 2000).

Figure 17: Left panel: Planetary caustics of a system with a projected
separation of $1.6$ Einstein radius and mass ratio of $0.03$ (red). The size
of the caustic varies with $s$ and $q$. Right panel: Planetary caustics of a
system with a projected separation of $0.8$ Einstein radius and mass ratio of
$0.03$ (red). At this range of $s$, there are two smaller planetary caustics.
The distance between the lens and the center of the planetary caustic is equal
to $|s-1/s|$. The solid black line with the arrow in each panel represents the
path of the source through the caustics. The left panel is a reproduction of a
similar figure in Poleski et al. (2014). Table 1: The process taken to
calculate $s$ and $q$ is shown in this table.

One Deviation: |  1. $t_{p}$ is determined. 2. $u=\sqrt{{u_{0}}^{2}+{(t_{p}-t_{0}/t_{E})}^{2}}$ 3. $s-\frac{1}{s}-u$ is solved for $s$. 4. Secondary check determines if $s>1$ or $s<1$. 5. $q={(\frac{t_{E,p}}{t_{E}})}^{2}$.

---|---

Two Deviations: |  1. $t_{1}$ and $t_{2}$ are determined. 2. $u_{1,2}=\sqrt{{u_{0}}^{2}+{(t_{1,2}-t_{0}/t_{E})}^{2}}$ 3. $LX$ is found by geometry. 4. $|s-\frac{1}{s}|-LX$ is solved for $s$. 5. Secondary check determines if $s>1$ or $s<1$. 6. For $s>1$: Height of the caustic $\sim\frac{2\sqrt{q}}{s\sqrt{s^{2}-1}}\implies q$ is found. 7. For $s<1$: Height of the caustic $\sim\frac{s^{3}\sqrt{q}}{4}\sqrt{(\frac{2}{s})^{8}+27}\implies q$ is found.

---|---

Figures 18, 19, and 20, show three examples of the PSPL plus the busy function
fitted to planetary light curves. Figure 18 shows the fit for a system with
one deviation in its PSPL residual where the true projected separation of the
system is smaller than unity. The algorithm first detects sharp peaks or
troughs, and after fitting the busy function to the light curve residuals from
the PSPL fit, we then obtain the residual of the busy function fit. At the
times where a single peak is found, we still cannot conclude that the event
includes a major image perturbation. At this point, if the sum of the flux of
negative points in the busy function fit residual is larger than the sum of
the flux of the positive points, the system will be considered to be affected
by a minor image perturbation, and $s$ will be chosen to be less than one.
Otherwise, $s$ will be larger than one. The example in Figure 18, shows an
event with a sharp positive deviation that appears to be a major image
perturbation; however, the preceding trough indicates that this a minor image
perturbation, and the algorithm correctly decides that $s$ is smaller than
one. The true parameters for this event are $s=0.63$ and $q=0.0003$ and its
fitted values are $s=0.64$ and $q=0.0002$.

Figure 18: An example of the PSPL plus the busy function to the light curve of
an event with $s=0.63$ and $q=0.0003$. The algorithm detects a sharp peak and
would simply decide that the system is caused by a major image if no secondary
checks are done. After a secondary check, it detects the negative deviations
from the model, and decides that $s$ should be less than unity. The estimated
parameters are $s=0.64$ and $q=0.0002$.

Figure 19 shows the PSPL plus the busy function fitted to a double-peaked
deviation with $s<1$. The busy function fit may appear to be a poor fit to the
deviations but this fact is not considered a failure, and in fact helps with
determining whether it is a major or a minor image perturbation. After the
secondary check, the algorithm decides to correctly choose $s$ smaller than
one. The true parameters for this event are $s=0.62$ and $q=0.0001$ and the
fitted values are $s=0.59$ and $q=0.0003$. Figure 20 shows another double-
peaked event that has $s>1$. The true parameters for this event are $s=1.11$
and $q=0.0008$ and the fitted values are $s=1.02$ and $q=0.007$.

Figure 19: An example of the PSPL plus the busy function to the light curve of
a double-peaked event with $s=0.62$ and $q=0.0001$. The busy function gives us
a good estimation of the location and duration of the event, and by doing a
secondary check, the large negative deviation from the model is detected, and
the algorithm decides that $s$ should be less than unity. The estimated
parameters are $s=0.59$ and $q=0.0003$.

Figure 20: An example of the PSPL plus the busy function to the light curve of
a double-peaked event with $s=1.11$ and $q=0.0008$. The busy function provides
a good estimate of the location and duration of the event, with the estimated
system parameters of $s=1.02$ and $q=0.007$.

In order to evaluate the performance of the algorithm, we compare the
calculated values of $s$ and $q$ with their true values as demonstrated by
Khakpash et al. (2019). Figure 21 shows the estimated values of $s$ and $q$
plotted versus their true values. Most values of $s$ are well estimated by the
method. The fitted values close to one arise from events caused by the central
caustics. While the estimated Values of $q$ show a lot of scatter compared to
the true value, they are generally estimated to within about one order of
magnitude of the true values. Note that this algorithm is slower than the
algorithm presented in Khakpash et al. (2019), but it fits a broader range of
events.

Figure 21: Left panel shows plot of fitted values of $q$ versus the true
values. Fitted values of $q$ are mostly within one order of magnitude of the
true values for the PSPL plus busy function fits

. The right panel shows fitted $s$ versus the true values. Most values of $s$
are well estimated by the method. The fitted values close to one occur for
events caused by the central caustics for which our formalism breaks down.
Blue dots are events without caustic crossings, and black dots are events with
caustic crossings.

###  4.7 Chebyshev Polynomials Fit

Di Stefano & Perna (1997) suggests that unlike binary-lens microlensing light
curves that are significantly perturbed, smooth binary-lens events can be
easily misclassified, and therefore, it is useful to develop methods that can
distinguish between this type of microlensing events and other types of
similar variabilities. The method should work fast and be effective at marking
light curves with smooth features. One of the approaches to approximate a
function is to expand it in form of a series of polynomials. For this purpose,
Di Stefano & Perna (1997) suggests to use the Chebyshev approximation on
microlensing light curves. The idea is that features of Chebyshev polynomials
are useful for capturing the smooth features of binary-lens microlensing light
curves without caustic crossings.

We have implemented the work of Di Stefano & Perna (1997), and applied it to
our simulated data set. For this purpose, we first detect the peaks in the
light curve, then for each peak, we choose the interval around the peak that
includes the wings of the event up to a magnification of 1.06 that corresponds
to an impact parameter of $2R_{E}$ (Di Stefano & Perna, 1997).

We then find the coefficients of the Chebyshev polynomials using the
formulation described by Press et al. (1992) and expand the selected parts of
the light curves by the expression in Equation 15. The polynomial $T_{k}$ has
$k+1$ number of extrema in the interval $[-1,1]$ where the Chebyshev
polynomials are defined. To use this method, the time coordinates of the
selected potion of our light curves should be within this interval, and so we
convert the interval of time to be between $-1$ and $1$ before fitting to the
Chebyshev polynomials.

$f(x)\approx\left[\sum_{k=0}^{m}c_{k}T_{k}(x)\right]-\frac{1}{2}c_{0}$ (15)

In this work, we choose the first 50 polynomials of this series ($m=50$) to
fit Equation 15 to the light curves, and then using the coefficients $c_{k}$,
we calculate
$\Lambda=-{log}_{10}((\sum_{k=0}^{m}{(\frac{c_{k}}{c_{0}})}^{2})-1)$ as a
feature parameter that can be used to distinguish between different light
curve types in our simulated data set.

In Figure 22, an example of a binary-lens event approximated by the Chebyshev
polynomials of degree 50 can be seen. It is important to note that the whole
light curve is not approximated in this method, and only the segment
containing the event excluding gaps is selected. This is because the Chebyshev
approximation is good for approximating functions that have finite number of
extrema in the interval of $[-1,1]$, and a non-varying light curve is not
usually well fitted by this approximation.

Figure 22: The stellar binary-lens light curve is approximated with the
Chebyshev polynomials of degree 50.

With this method, we now need to define Chebyshev parameters that are
different for various types of light curves. Here, we use $\Lambda$ and the
first six even coefficients as the parameters used to distinguish between
different types of light curves. The reason for that is that the coefficients
would decrease rapidly when going to higher order terms in the polynomials,
and therefore, the first coefficients would be dominant. Also, since the
microlensing events are closer to being an even function rather than an odd
function, the expansion only involves the even terms in the Chebyshev
polynomials.

We then check how well different values of $\Lambda$ match to different light
curve types. Figure 23 shows that the value of $\Lambda$ correlate fairly well
with four different light curve types. Note that the fractions are calculated
for a given range of $\Lambda$, across each light curve types, horizontally in
the figure. Although the fractions of different variability types in the
simulated Roman data set are not equal, there are large numbers of each
variability type, so the distinctions seen in Figure 23 are still quite
robust. We therefore expect that the continuous parameter $\Lambda$, can be a
useful input for a ML algorithm (such as a Random Forest) to classify event
types.

Figure 23:  This figure shows the breakdown by variability type for each range
of values of the Chebyshev feature $\Lambda$. $\Lambda$ is equal to
$-{log}_{10}((\sum_{k=0}^{m}{(\frac{c_{k}}{c_{0}})}^{2})-1)$, where all
coefficients are normalized by the first coefficient $c_{0}$.

In addition to the feature $\Lambda$, we also follow the suggestion of Di
Stefano & Perna (1997), to include the individual even-numbered coefficients
of the Chebyshev polynomial as classification features.

##  5 An Algorithmic Approach

In this section, we introduce a suggested algorithmic approach to use the
produced features as input to ML classifiers. The goal of this section is to
provide an efficient procedure to use these features to train classifiers on a
high-cadence data set and detect different types of microlensing light curves.
In Table 2, we have a summary of the tools in our package and the associated
features. Each of these tools can be either applied to the whole dataset or a
subset of it, and there may be more than a single algorithmic approach to use
these features.

Table 2: Summary of the introduced algorithms and their resulting features. Algorithms | Produced Features
---|---
Peak Finder | Number of the Peaks
G-PSPL Fit |  ${\chi}^{2}$: Goodness of the G-PSPL fit
Symmetry Check |  $\beta$: A measure of asymmetry of the peak
Trapezoidal Function Fit |  $\kappa$: Ratio of the duration of the flat top of the trapezoidal function to total duration
|  $t_{E,trap}$: Duration of the trapezoidal function
Cauchy Distribution Fit |  $\psi$: The difference between the goodness of PSPL and Cauchy fits
|  $b$: A measure of the flatness of the top
Planetary Parameter Finder |  $s$: Projected mass ratio in units of Einstein radius
|  $q$: Mass ratio
Chebyshev Polynomials Fit |  $\Lambda$: Sum of the square roots of the Chebyshev coefficients
|  $a_{2}$, $a_{4}$, $a_{6}$, $a_{8}$, $a_{10}$: First five even coefficients
of Chebyshev polynomials

In a paper under development (Khakpash et al, in preparation), we are
implementing an algorithmic approach to employing the features discussed in
this paper to comprehensively search for microlensing events, classify them by
type, and derive preliminary system parameters. For any set of lightcurve
features, there are numerous ways to conduct classification and our approach
is by no means certain to be optimal in every way. We encourage other
astronomers to make use of the features described here to develop independent
classification methods.

The particular structure or sequence of a classification approach can vary.
One might use all the features to detect all the categories at once, or
alternatively, find particular classes by doing a step-by-step classification.
In an ideal case where there are thousands examples of each class across a
full range of empirical properties of each feature, it is likely better to use
all the features to find all of the classes in a one-step classification.
Since in this case we have very limited examples of some of the classes, we
recommend a step-by-step classification approach.

The first step of classification is to distinguish microlensing light curves
from other types of variability (e.g., CV in our dataset). We refer to this
step as classification step I. Although our goal in this work is not focused
on this set, we believe some of the features identified in this work can be
used to improve the current existing classifiers focused on this task. We have
tested this type of classification using the features such as number of the
peaks found by the Peak Finder, $\Lambda$, $a_{2}$, $a_{4}$, $a_{6}$ produced
by the Chebyshev fit, and ${{\chi}^{2}}_{PSPL}$, ${t_{E,PSPL}}$, $\beta$
generated by the G-PSPL fit and the Symmetry Check.

Assuming we have identified all the microlensing light curves, the next type
of classification is to classify them into single-lens versus binary-lens
events which we refer to as step II. Note that in a real dataset, multi-lens
events also exist which can either be added as a category or can be included
in one category along with the binary-lens events. We use the same feature as
in step I excluding the number of peaks.

Once we have single-lens and binary-lens events, we classify the binary-lens
events into stellar binary-lens and planetary binary-lens systems
(classification step III) using ${{\chi}^{2}}_{PSPL}$, ${t_{E,PSPL}}$ produced
by the PSPL fit and $s$, and $q$ produced by the PSPL plus busy fit.

Furthermore, single-lens events can be classified into classes of isolated
short timescale events likely caused by free-floating planets or planets on
very wide orbits, stellar PSPL, and stellar FSPL events (Classification step
IV). We suggest using $\kappa$, $\Delta\tau_{full}$, $\psi$, and $b$ produced
by the Trapezoidal Fit and the Cauchy/PSPL Fit to obtain better results. This
step will be investigated in the future work.

Figure 24 displays a diagram showing our suggested algorithmic approach to use
the different features produced by our package. We have selected this approach
to optimize the classification results considering our small dataset. We have
tested classification steps I, II, and III with four ML classifiers including
a k-nearest neighbors classifier, a decision tree classifier, a random forest
classifier and a neural network classifier. The preliminary results for that
are presented in the next section.

Figure 24:  This diagram shows an algorithmic approach for using the features
described in this paper

as input to ML classifiers.

##  6 Preliminary Testing of ML Algorithms

As stated above, there are many ways to implement an approach to identify and
characterize microlensing events with classification based on light curve
features. We are developing a thorough investigation into that question
(Khakpash, et al. in preparation), but for now we present a preliminary
analysis using a few simple analysis steps based on the classification
approach described in §5. We present the results of training four ML
classifiers including a k-nearest neighbor classifier (KNN), a decision tree
classifier (DT), a random forest classifier (RF), and a neural network
classifier (NN) using the features we introduced in this paper. In order to
test these algorithms, we follow the step-by-step scheme of Figure 24. It is
important to note again that the step IV of the classification in Figure 24 is
not tested here and will be pursued in the future. At each step, we set the
size of our test set to be 20% of the whole dataset. The test set is then
randomly chosen in a 5-fold cross-validation process, and the average scores
are reported at the end.

In order to compare results if these classifiers, we show confusion matrices
made with the test set at each step along with their Receiver operating
characteristic (ROC) curves. As mentioned in Section 5, we understand that the
one-step classification is a more ideal approach, and we tested this approach
with our current data set. However, the limited number of object instances in
each class was insufficient for achieving an overall acceptable accuracy.
Nevertheless, the isolated confusion matrices of the step-by-step
classifications are valuable to evaluate the utility of the features presented
in this paper. We are planning to thoroughly investigate the one-step
classification in a future paper.

At each step, we use the same dataset and features for all of the four
classifiers, but we optimize their hyperparameters separately. The RF, KNN,
and DT classifiers are implemented using the scikit-learn package in python
(Pedregosa et al., 2011). The NN classifier is implemented using the
Tensorflow and Keras packages in python (Abadi et al., 2015; Chollet et al.,
2015). A summary of the features and hyperparameters of each classifier at
each step is given in Table 3.

###  6.1 Classification Step I

As shown in Figure 24, the first step of the algorithmic approach is to detect
microlensing light curves among other stellar variability. We should first
note that our current light curve data set is not completely representative of
what we expect from the Roman mission, since the non-microlensing type of
variability in the light curves in our current dataset only includes
cataclysmic variables. Most (but not all) forms of non-microlensing
variability are expected to be periodic or quasi-periodic, and thus simply
including CVs is a good starting point. Our dataset for this set contains 4181
light curves among which there are 3752 microlensing light curves (labeled as
1) and 429 non-microlensing light curves (labeled as 0).

We find that the test set and training set accuracy for all of the four
classifiers in this step are very close. A common way to evaluate the results
of a classification model is to plot a confusion matrix for it. A confusion
matrix is a table containing the percentages of both correctly and incorrectly
classified objects for each class in the dataset. According to the confusion
matrices shown in Figure 25, most of ML tools can find the microlensing light
curves with very small classification error, whereas, about 40% of the CV
light curves are misclassified. This could be a result of having a small
training set, or incomplete hyperparameter tuning, or might be an indication
of the need to include more features.

A more robust method of comparing different ML classifiers is to plot their
ROC curves, and calculate their values of the Area Under the Curve (AUC). The
closer the AUC is to unity, the better the performance of that classifier is.
This includes plotting True Positive Rate (TPR) versus False Positive Rate
(FPR) for different decision thresholds, and finding the area under that.
Figure 26 shows the ROC curves of the four classifiers trained in step I. The
diagonal dashed line represents the ROC curve of a random classification. The
ROC curves show that the NN and RF classifiers have a better performance and
can achieve a higher TPR without lowering the FPR.

(a) Random Forest Classifier

(b) Neural Networks Classifier

(c) Decision Tree Classifier

(d) K-nearest Neighbors

Figure 25: Confusion matrices of the four trained classifiers of
classification Step I. In this step, we aim at classifying all of the light
curves into two classes of CV and microlensing.

Figure 26: ROC curves of the four trained classifiers of classification Step I
along with their AUC values. In this step, we aim at classifying all of the
light curves into two classes of CV and microlensing. RF and NN have the
largest AUC implying that they are able to achieve higher TPR while the FPR is
also low.

###  6.2 Classification Step II

The second classification step is distinguishing between single-lens and
binary-lens microlensing light curves as shown in Figure 24. Our dataset for
this step contains 4181 light curves among which there are 2143 binary-lens
microlensing light curves (labeled as 1) and 1626 single-lens microlensing
light curves (labeled as 0). For this step, we find that DT and RF have
similar training and test accuracy and seem to work better than NN and KNN,
although NN seems to work much better than KNN.

Figure 27 shows confusion matrices of the four classifiers trained for this
step. The confusion matrices suggest that RF works better at predicting both
labels, whereas NN works best at finding a larger fraction of the binary-lens
light curves. Figure 28 shows the ROC curves and AUC values of this step. The
classifiers are overall working better in this step mainly because the dataset
is more balanced here.

(a) Random Forest Classifier

(b) Neural Networks Classifier

(c) Decision Tree Classifier

(d) K-nearest Neighbors

Figure 27: Confusion matrices of the four trained classifiers of
classification Step II. At this stage, we assume that the microlensing light
curves are already detected, and our goal is to classify them into groups of
single-lens and binary-lens microlensing light curves.

Figure 28: ROC curves of the four trained classifiers of classification Step
II along with their AUC values. In this step, we are classifying the
microlensing light curves into groups of single-lens and binary-lens
microlensing light curves. The dataset is more balanced in this step and all
of the classifiers seem to work well. RF and NN still show a better
performance.

###  6.3 Classification Step III

After finding the binary lens light curves, the next classification step is
distinguishing between stellar binary-lens and planetary binary-lens
microlensing light curves as shown in Figure 24. Our dataset for this step
contains 2074 light curves among which there are 688 planetary binary-lens
microlensing light curves (labeled as 1) and 1386 stellar binary-lens
microlensing light curves (labeled as 0).

This classification step is particularly important in this context since
planetary binary-lens systems are the ones that astronomers would like to
distinguish from the rest of the dataset. In our tested example the number of
light curve is lower than the previous classification steps and this decreases
the accuracy of the classifiers. Because of this, we find that the overall
accuracy values of this step are smaller than the previous steps.
Additionally, stellar and planetary binary-lens systems are much less
distinguishable from each other compared to the previous tasks, and for this
reason the simpler algorithms of DT and KNN appear to have lower accuracy. RF
and NN have higher overall accuracy, but NN has a higher test accuracy which
results in a larger fraction of the test set being correctly labeled. It seems
that an algorithm like NN is more capable of distinguishing between these two
categories which is expected as NN is theoretically more complex and is
designed to find complicated patterns in a data set.

Figure 29 shows confusion matrices of the four classifiers. The confusion
matrix of the NN shows a larger value of TPR compared to all the other
confusion matrices, and this is more favorable since our ultimate goal is to
detect planetary microlensing light curves. Figure 30 shows the ROC curves and
AUC values of the four classifiers in step III. The dataset in this step is
smaller compared to the other two steps and is not well balanced. Therefore,
the performance of the different classifiers are not as well-differentiated as
in the other steps. However, NN shows significantly better performance
compared to the other classifiers.

(a) Random Forest Classifier

(b) Neural Networks Classifier

(c) Decision Tree Classifier

(d) K-nearest Neighbors

Figure 29: Confusion matrices of the four trained classifiers of
classification Step III. Once the binary-lens microlensing light curves are
found, we use this classification step to classify them into stellar and
planetary binary-lens microlensing light curves.

Figure 30: ROC curves of the four trained classifiers of classification Step
III along with their AUC values. In this step, we aim at classifying binary-
lens light curves into two classes of stellar and planetary binary-lens
microlensing. Since the dataset is smaller and not well balanced, the overall
performance of the classifiers in this step is lower than the other two steps.
However, NN shows significantly better performance than the other classifiers.

Table 3: Classifiers Hyperparameters and Features.

| Features | KNN | DT | RF | NN
---|---|---|---|---|---
Step I |  |  $\Lambda$, $a_{2}$, $a_{4}$, $a_{6}$, $\beta$,
---
Number of peaks
(binsize = 60,
threshold = 4 & 6),
${\chi}^{2}_{PSPL}$ , $t_{E,PSPL}$
| metric = “manhattan”,
---
n_neighbors = 4
| criterion = “gini”,
---
max_depth = 5,
splitter = “best”
| max_features = 7,
---
min_samples_leaf = 1
n_estimators = 100,
criterion = “entropy”,
bootstrap = False
| 4 Dense Layers,
---
3 Dropout Layers,
activation = “relu”
batch = 32,
epochs = 300
Step II |  |  $\Lambda$, $a_{2}$, $a_{4}$, $a_{6}$, $\beta$,
---
${\chi}^{2}_{PSPL}$ , $t_{E,PSPL}$
| metric = “euclidean”,
---
n_neighbors = 20
| criterion = “entropy”,
---
max_depth = 5,
splitter = “best”
| max_features = 7,
---
min_samples_leaf = 1
n_estimators = 100,
max_depth = 5
| 3 Dense Layers,
---
2 Dropout Layers,
activation = “relu”
batch = 16,
epochs = 300
Step III |  |  $\Lambda$, $a_{2}$, $a_{4}$, $a_{6}$,
---
${\chi}^{2}_{PSPL}$, $t_{E,PSPL}$, $s$, $q$
| metric = “minkowski”,
---
n_neighbors = 3, p = 15
| criterion = “gini”,
---
max_depth = 7,
splitter = “best”
| max_features = 8,
---
min_samples_leaf = 1
n_estimators = 300,
max_depth = 8
| 4 Dense Layers,
---
2 Dropout Layers,
activation = “tanh”
batch = 32,
epochs = 100

##  7 Conclusion and Future Work

Classifying light curves that manifest different types of stellar variability
is still a major challenge. Although using ML methods to classify astronomical
time series is a powerful tool, it is important to understand which method we
should choose and what are the steps we need to take in order to use these
methods effectively.

In this paper, we introduce a new approach towards classifying microlensing
light curves based on light curve morphologies. We introduce a package of
tools including several functional fits that can be applied to the light curve
in a fast and efficient way to extract information about the different
morphological features in the light curves. This information is quantified as
features that can be used to make decisions about the light curve types.

This approach will be useful when it comes to analyzing large microlensing
data sets from high-cadence surveys like Roman (Spergel et al., 2015) in the
future and ongoing surveys like KMTNet (Kim et al., 2010). Our preliminary
results in Section 6 show that the features produced by our tools can be used
as input to ML classifiers like Random Forest to distinguish between different
types of microlensing light curves, and help prioritize the ones that are more
likely to be caused by planetary systems.

The simulated data used in this paper included CV-like events as the only non-
microlensing instance. An ideal dataset would include a variety of stellar
variability, and our developed package should be modified to include all other
variability in its analysis. We believe that the same approach will be
successful in recovering microlensing events from a wide variety of non-
microlensing variability on its own.

There are currently a number of methods that attempt to categorize photometric
variability in large data sets, such as parametric statistical methods and ML.
ML methods for detecting different types of variabilities are becoming more
common (e.g. Richards et al., 2011; Pichara, 2013; Pashchenko, 2017;
Valenzuela, 2017). For example, as mentioned before, Godines et al. (2019) has
trained a Random Forest classifier to detect microlensing light curves in a
low-cadence survey dataset in real time. Our produced features can be a
complementary set of features for such algorithms and not only can improve
their results but also can give them the ability to detect possible planetary
binary-lenses as well.

In ML classifiers, increasing the number of object instances can significantly
improve the results. Our largest training set included 4181 light curves,
which is not a large number for most ML applications. Increasing this number
to about $\sim 10,000$ light curves would yield more robust and reliable
results. Some of the tools presented in this paper produce other parameters as
well, and including different subsets of those parameters could also improve
the results. This task needs to be done carefully, though, since adding more
features that are not important might result in overfitting. Additionally, a
great avenue to improve the results would be to test other classifiers like
the Support Vector Classifier and Naive Bayesian, and also investigate deeper
neural networks. The improvement of the data set and the ML algorithms will be
presented in a second paper of this series.

##  8 Acknowledgement

S.K thanks the LSSTC Data Science Fellowship Program, which is funded by
LSSTC, NSF Cybertraining Grant $\\#1829740$, the Brinson Foundation, and the
Moore Foundation; her participation in the program has benefited this work.
She also thanks the International Microlensing Conference that took place in
2018 in New York City, NY, where there were discussions and works that
initiated this project. R.S gratefully acknowledges support from NASA grant
80NSSC19K0291. B.S.G. was supported by a Thomas Jefferson Chair for Space
Exploration endowment at the Ohio State University.

## References

  * Abadi et al. (2015) Abadi, M., Agarwal, A., Barham, P., et al. 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, software available from tensorflow.org
  * Belokurov et al. (2003) Belokurov, V., Evans, N. W., & Du, Y. L. 2003, Monthly Notices of the Royal Astronomical Society, 341, 1373
  * Bennett et al. (2002) Bennett, D., Becker, A. C., Quinn, J., et al. 2002, The Astrophysical Journal, 579, 639
  * Bluck et al. (2020) Bluck, A. F., Maiolino, R., Sánchez, S. F., et al. 2020, Monthly Notices of the Royal Astronomical Society, 492, 96
  * Bond et al. (2001) Bond, I., Abe, F., Dodd, R., et al. 2001, Monthly Notices of the Royal Astronomical Society, 327, 868
  * Bozza (2000) Bozza, V. 2000, Journal of Mathematical Physics, 41, 6284
  * Chollet et al. (2015) Chollet, F., et al. 2015, Keras, https://keras.io
  * Di Stefano & Perna (1997) Di Stefano, R., & Perna, R. 1997, The Astrophysical Journal, 488, 55
  * Gaudi (2010) Gaudi, B. S. 2010, arXiv preprint arXiv:1002.0332
  * Gaudi (2012) —. 2012, Annual Review of Astronomy and Astrophysics, 50, 411
  * Godines et al. (2019) Godines, D., Bachelet, E., Narayan, G., & Street, R. 2019, Astronomy and Computing, 28, 100298
  * Gould (1996) Gould, A. 1996, The Astrophysical Journal, 470, 201
  * Gould & Gaucherel (1996) Gould, A., & Gaucherel, C. 1996, arXiv preprint astro-ph/9606105
  * Griest et al. (1995) Griest, K., Alcock, C., Allsman, R., et al. 1995, arXiv preprint astro-ph/9506016
  * Han (2006) Han, C. 2006, The Astrophysical Journal, 638, 1080
  * Han et al. (2016) Han, C., Bennett, D. P., Udalski, A., & Jung, Y. K. 2016, The Astrophysical Journal, 825, 8
  * Jiang et al. (2004) Jiang, G., DePoy, D., Gal-Yam, A., et al. 2004, The Astrophysical Journal, 617, 1307
  * Johnson et al. (2020) Johnson, S. A., Penny, M. T., Gaudi, B. S., et al. 2020, arXiv preprint arXiv:2006.10760
  * Kessler et al. (2019) Kessler, R., Narayan, G., Avelino, A., et al. 2019, arXiv preprint arXiv:1903.11756
  * Khakpash et al. (2019) Khakpash, S., Penny, M., & Pepper, J. 2019, The Astronomical Journal, 158, 9
  * Kim et al. (2018) Kim, D.-J., Kim, H.-W., Hwang, K.-H., et al. 2018, The Astronomical Journal, 155, 76
  * Kim et al. (2010) Kim, S.-L., Park, B.-G., Lee, C.-U., et al. 2010, in Ground-based and Airborne Telescopes III, Vol. 7733, International Society for Optics and Photonics, 77333F
  * Kim et al. (2016) Kim, S.-L., Lee, C.-U., Park, B.-G., et al. 2016, Journal of the Korean Astronomical Society, 49, 37
  * Liaw & Wiener (2002) Liaw, A., & Wiener, M. 2002, R News, 2, 18
  * Lloyd (1982) Lloyd, S. 1982, IEEE transactions on information theory, 28, 129
  * Mao & Di Stefano (1994) Mao, S., & Di Stefano, R. 1994, Interpretation of gravitational microlensing by binary systems, Tech. rep., SCAN-9411403
  * Mao et al. (2002) Mao, S., Smith, M. C., Woźniak, P., et al. 2002, Monthly Notices of the Royal Astronomical Society, 329, 349
  * Mróz et al. (2017) Mróz, P., Udalski, A., Skowron, J., et al. 2017, Nature, 548, 183
  * Mróz et al. (2018) Mróz, P., Ryu, Y.-H., Skowron, J., et al. 2018, The Astronomical Journal, 155, 121
  * Mróz et al. (2020) Mróz, P., Poleski, R., Han, C., et al. 2020, The Astronomical Journal, 159, 262
  * Nemiroff & Wickramasinghe (1994) Nemiroff, R. J., & Wickramasinghe, W. 1994, arXiv preprint astro-ph/9401005
  * Pashchenko (2017) Pashchenko. 2017, Monthly Notices of the Royal Astronomical Society, 475, 2326
  * Pawlak (2019) Pawlak, M. 2019, in Astronomical Data Analysis Software and Systems XXVI, Vol. 521, 319
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, the Journal of machine Learning research, 12, 2825
  * Penny et al. (2013) Penny, M., Kerins, E., Rattenbury, N., et al. 2013, Monthly Notices of the Royal Astronomical Society, 434, 2
  * Penny (2014) Penny, M. T. 2014, The Astrophysical Journal, 790, 142
  * Penny et al. (2019) Penny, M. T., Gaudi, B. S., Kerins, E., et al. 2019, The Astrophysical Journal Supplement Series, 241, 3
  * Penny et al. (2011) Penny, M. T., Mao, S., & Kerins, E. 2011, Monthly Notices of the Royal Astronomical Society, 412, 607
  * Pichara (2013) Pichara. 2013, The Astrophysical Journal, 777, 83
  * Poleski et al. (2014) Poleski, R., Udalski, A., Dong, S., et al. 2014, The Astrophysical Journal, 782, 47
  * Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, The art of scientific computing, 2, 1002
  * Richards et al. (2011) Richards, J. W., Starr, D. L., Butler, N. R., et al. 2011, The Astrophysical Journal, 733, 10
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv preprint arXiv:1503.03757
  * Street et al. (2018) Street, R., Lund, M., Donachie, M., et al. 2018, arXiv preprint arXiv:1812.04445
  * Sumi et al. (2011) Sumi, T., Kamiya, K., Bennett, D., et al. 2011, Nature, 473, 349
  * Udalski (2004) Udalski, A. 2004, arXiv preprint astro-ph/0401123
  * Udalski et al. (2015) Udalski, A., Szymański, M., & Szymański, G. 2015, arXiv preprint arXiv:1504.05966
  * Valenzuela (2017) Valenzuela. 2017, Monthly Notices of the Royal Astronomical Society, 474, 3259
  * Westmeier et al. (2014) Westmeier, T., Jurek, R., Obreschkow, D., Koribalski, B. S., & Staveley-Smith, L. 2014, Monthly Notices of the Royal Astronomical Society, 438, 1176
  * Witt & Mao (1994) Witt, H. J., & Mao, S. 1994, The Astrophysical Journal, 430, 505
  * Woźniak & Paczyński (1997) Woźniak, P., & Paczyński, B. 1997, The Astrophysical Journal, 487, 55
  * Yoo et al. (2004) Yoo, J., DePoy, D., Gal-Yam, A., et al. 2004, The Astrophysical Journal, 603, 139
  * Zhu & Gould (2016) Zhu, W., & Gould, A. 2016, arXiv preprint arXiv:1601.03043