# Constraining Black Hole Natal Kicks with Astrometric Microlensing

Jeff J. Andrews Center for Interdisciplinary Exploration and Research in
Astrophysics (CIERA), 1800 Sherman Ave., Evanston, IL, 60201, USA  Physics
Department, University of Florida, 2001 Museum Rd., Gainesville, FL 32611
<EMAIL_ADDRESS>Vicky Kalogera Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA), 1800 Sherman Ave., Evanston,
IL, 60201, USA  Department of Physics and Astronomy, Northwestern University,
2145 Sheridan Rd., Evanston, IL 60208, USA

###### Abstract

Multiple pieces of evidence suggest that neutron stars receive large kicks
when formed from the remnant of a collapsing star. However, the evidence for
whether black holes (BH) receive natal kicks is less clear, reliant on weak
constraints from the analysis of BH X-ray binaries and massive runaway and
walkaway stars. Here we show for the first time that recent microlensing
detections offer a new method for measuring the kicks BHs receive at birth.
When a BH is identified through both photometric and astrometric microlensing
and when the lensed star has a known distance and proper motion, the mass,
distance and proper motion of the BH can be determined. We study the runaway
velocities for components of eccentric binaries disrupted during a supernova,
finding the peculiar velocity correlates strongly with the kick a BH received
a birth, typically within 20%, even when the natal kick is smaller than the
orbital velocity. Therefore, by measuring the peculiar velocity of a BH or
other compact object that formed from a binary that disrupted during core
collapse, we are in effect measuring the natal kick that object received. We
focus on MOA-2011-BLG-191/OGLE-2011-BLG-0462, an isolated, single BH detected
by microlensing and consider a range of possible formation scenarios,
including its formation from the disruption of a binary during a supernova
event. We determine that MOA-2011-BLG-191/OGLE-2011-BLG-0462 has a Milky Way
orbit consistent with a thick disk population, but if it were formed within
the kinematic thin disk it received a natal kick $\lesssim$100 km s-1.

black hole physics—binaries: general—stars: black holes

††journal: ApJ††software: astropy (Astropy Collaboration et al., 2018), gala
(Price-Whelan et al., 2021), matplotlib (Hunter, 2007), NumPy (van der Walt et
al., 2011), SciPy (Jones et al., 2001)

##  1 Introduction

As some of the most energetic phenomena in the Universe, supernovae (SN)
release $\sim$1051 erg s-1 into the surrounding environment in seconds
(Woosley & Weaver, 1986); even slight asymmetries can have significant
consequences for the kinematics of the neutron star (NS) and black hole (BH)
remnants left behind (Shklovskii, 1970; Janka, 2013; Wongwathanarat et al.,
2013). By measuring the proper motions of a sample of pulsars, Lyne & Lorimer
(1994) first and then later Hobbs et al. (2005) determined that NSs receive
kicks at birth of several hundred km s-1. Careful analysis of pulsar binary
systems largely confirm the requirement for NSs to be born with large kick
velocities (see Wong et al., 2010; Tauris et al., 2017, and references
therein).

It is not clear whether BHs are born with similarly large kick velocities.
Current constraints on BH natal kicks rely on a variety of data sets. The most
stringent limits are placed when analyzing individual X-ray binaries (XRBs)
known to host BH accretors (Mirabel et al., 2001; Jonker & Nelemans, 2004;
Willems et al., 2005; Fragos et al., 2009; Repetto et al., 2012; Wong et al.,
2012, 2014; Repetto & Nelemans, 2015; Repetto et al., 2017; Mandel, 2016). At
different levels of sophistication, these studies use a combination of binary
stellar evolution models and Milky Way orbital modelling to uncover the
formation of individual systems, thereby placing limits on the velocity of BH
natal kicks. Although the firmness of the constraints differ for each modelled
system, it appears that at least one XRB provides strong evidence that some
BHs receive natal kicks (e.g., $>$80 km s-1 for XTE J1118$+$480; Fragos et
al., 2009).

At the same time, at least two “failed” SN candidates have been identified
(Kochanek et al., 2008; Gerke et al., 2015; Reynolds et al., 2015) in which a
massive star disappeared, implying these stars collapsed into a BH but
emitting either very little or none of the radiation characteristic of a SN
(Smartt, 2009). To produce a failed SN, a massive star must undergo a prompt
collapse, likely without developing the asymmetries required to impart large
kicks, a possibility supported by some SN models (Lovegrove & Woosley, 2013;
Sukhbold et al., 2016; Patton & Sukhbold, 2020). Indeed, Mirabel & Rodrigues
(2003) used a combination of orbital dynamics, proper motion measurements, and
a lack of detected supernova remnant to suggest that the well-studied XRB
Cygnus X-1 may have formed through such a mechanism.

Recently, Callister et al. (2021) and Stevenson (2022) have used the catalog
of merging binary BHs detected by gravitational-wave observatories (Abbott et
al., 2021; The LIGO Scientific Collaboration et al., 2021) to show that, if
the population is formed through isolated binary evolution then BHs must
receive significant natal kick velocities ($\sim$100 km s-1) in order to
explain the observed spin and spin-orbit misalignment angle distributions.

Finally, a number of studies have focused on a subset of massive stars that
move with abnormally large velocities, exploring the possibility that these
runaway and walkaway stars are the former binary companions to collapsing
stars (Blaauw, 1961; van Rensbergen et al., 1996; De Donder et al., 1997; Dray
et al., 2005; Eldridge et al., 2011; Boubert & Evans, 2018; Renzo et al.,
2019; Aghakhanloo et al., 2022). In this scenario, the mass loss from the
collapsing star and possible natal kick imparted to the newly formed compact
object disrupted the binary, leaving behind the unbound companion. Renzo et
al. (2019) in particular focus on the effects that BH kicks have on the masses
and kinematics of such unbound OB stars. Although the effects of BH kicks on
massive star populations can be significant, placing constraints using
observed populations is challenging; observational selection effects are
significant and any constraint is population dependent, requiring an accurate
knowledge of massive binary star population characteristics at birth.

Recently, a new observational data set affords the opportunity to measure the
kicks that BHs receive at birth: BHs detected through microlensing. Several
free-floating BH candidates have been identified through the long-timescale
magnification they induce on background stars (Bennett et al., 2002; Mao et
al., 2002; Wyrzykowski et al., 2016; Wyrzykowski & Mandel, 2020); however
these have all been measured photometrically only. Recently, Sahu et al.
(2022) and Lam et al. (2022) have both separately presented six years of data,
both photometric and astrometric, of MOA-2011-BLG-191/OGLE-2011-BLG-0462
(hereafter MOA-11-191/OGLE-11-0462). The astrometric data in particular breaks
degeneracies in the model parameters and allows for the measurement of the
lens’s mass, distance, and transverse velocity. Using the same data, the
groups find different solutions for the lens properties. Sahu et al. (2022)
report a 7.1$\pm$1.3 ${{M_{\odot}}}$ BH at a distance of 1.58$\pm$ 0.18 kpc,
moving with a tangential velocity of $\simeq$45 km s-1. Lam et al. (2022) find
MOA-11-191/OGLE-11-0462 is either a 3.7 ${{M_{\odot}}}$ BH or a 2.1
${{M_{\odot}}}$ NS or BH, depending on the relative weighting between the
photometric and astrometric data. Both Sahu et al. (2022) and Lam et al.
(2022) hint that the observed tangential velocity could be due to a kick the
BH received at birth. Here we aim to explore this possibility, with the goal
of inferring new constraints on BH natal kicks. In Section 2 we provide our
method for calculating the runaway velocities of stars in eccentric binaries,
while we describe the results of those calculations in Section 3. We consider
possible formation scenarios of MOA-11-191/OGLE-11-0462 in Section 5, and
provide some concluding thoughts in Section 6.

##  2 The Runaway Velocities of Compact Objects in Unbound Binaries

Several studies have considered the effect of instantaneous mass loss on a
binary’s orbit (e.g., Blaauw, 1961; Hills, 1970), finding the landmark result
that an initially circular binary will disrupt if it instantaneously loses at
least half the system’s total mass. Hills (1983) expanded upon this result by
considering an initially eccentric binary and including the effects of a SN
natal kick (a separate derivation is also provided by Hurley et al., 2002).
Hills (1983) showed that if mass loss occurs near pericenter, then an
eccentric binary could dissolve with substantially less mass lost, while near
apocenter, a binary could withstand more mass lost. However, Hills (1983) and
Hurley et al. (2002) only described the properties of the binaries that
remained bound. Tauris & Takens (1998) calculated the runaway velocities of
binary components, but for initially circular binaries (see also Gott et al.,
1970). The complete solution for the runaway velocities of stars disrupted
from an initially eccentric binary, including the effect of mass loss and
natal kick within the SN was separately provided by Belczynski et al. (2008)
and Kiel & Hurley (2009) (see also Pijloo et al., 2012). For completeness, we
reproduce that procedure, but in its non-dimensional form.

We consider a binary in which a helium star with mass $M_{\rm He}$ undergoes
core collapse, evolving into a BH with mass $M_{\rm BH}$, while the secondary
star has a constant mass $M_{2}$. Although it is not important that the BH
progenitor is a helium star for the sake of this derivation, the massive star
progenitors of BHs are expected to be massive, helium-rich Wolf-Rayet stars
(for a review, see Langer, 2012). It is useful to define two parameters
dependent on the system’s mass:

$\displaystyle\mu$ $\displaystyle=$ $\displaystyle\frac{M_{\rm
BH}M_{2}}{M_{\rm BH}+M_{2}}$ (1) $\displaystyle\beta$ $\displaystyle=$
$\displaystyle\frac{M_{\rm BH}+M_{2}}{M_{\rm He}+M_{2}}.$ (2)

At the instant of core collapse, the distance between the two components of a
binary, $\bm{r}_{0}$, and its relative orbital velocity, $\bm{v}_{0}$,
expressed in terms of the binary’s orbital separation, $a$, and circular
orbital velocity, $V_{\rm c}=\left[\mathcal{G}(M_{\rm
He}+M_{2})/a\right]^{1/2}$, respectively, are:

$\displaystyle\bm{r}_{0}$ $\displaystyle=$ $\displaystyle\left[\cos
E-e,\sqrt{1-e^{2}}\sin E,0\right]$ (3) $\displaystyle\bm{v}_{0}$
$\displaystyle=$ $\displaystyle\frac{1}{1-e\cos E}\left[-\sin
E,\sqrt{1-e^{2}}\cos E,0\right],$ (4)

where $e$ is the pre-SN orbital eccentricity and $E$ is the eccentric anomaly
of the pre-SN orbit at core collapse. We have assumed the orbital plane lies
in the $X-Y$ plane with pericenter aligned along the $X$ axis. Belczynski et
al. (2008) provide a schematic describing the orbital orientation (see their
Figure 4). The individual stellar velocities around the center of mass are
therefore:

$\displaystyle\bm{v}_{1}$ $\displaystyle=$ $\displaystyle\frac{M_{2}}{M_{\rm
BH}+M_{2}}\bm{v_{0}}$ (5) $\displaystyle\bm{v}_{2}$ $\displaystyle=$
$\displaystyle\frac{M_{\rm BH}}{M_{\rm BH}+M_{2}}\bm{v_{0}}.$ (6)

The helium star (subscript 1) receives a kick $\bm{v_{k}}$ (also expressed in
terms of $V_{\rm c}$).

The velocity of the post-SN center of mass is:

$\bm{v_{\rm CM}}=\frac{M_{\rm
BH}\left(\bm{v_{1}}+\bm{v_{k}}\right)+M_{2}\bm{v_{2}}}{M_{\rm BH}+M_{2}},$ (7)

and the post-SN relative velocity of the two stars is:

$\bm{v}=\bm{v_{1}}-\bm{v_{2}}+\bm{v_{k}}.$ (8)

The specific angular momentum, $j=\left|\bm{j}\right|$, (in terms of $a$,
$V_{\rm c}$, and $\mu$) is:

$\bm{j}=\bm{r_{0}}\times\bm{v}.$ (9)

The post-SN orbital energy can therefore be calculated:

$\mathcal{E}=\frac{v^{2}}{2}-\frac{\beta}{r_{0}},$ (10)

where $v=\left|\bm{v}\right|$ and $r_{0}=\left|\bm{r_{0}}\right|$. In this
work, we are only interested in systems that dissociate during core collapse,
i.e., $\mathcal{E}>0$.

The post-SN hyperbolic orbit is defined by an eccentricity, $e$, and semi-
latus rectum, $p$:

$\displaystyle e$ $\displaystyle=$
$\displaystyle\sqrt{1+2\mathcal{E}j^{2}/\beta^{2}}$ (11) $\displaystyle p$
$\displaystyle=$ $\displaystyle j^{2}/\beta.$ (12)

Depending on the orbital parameters, the orbital phase where core-collapse
occurs, and the kick velocity itself, the post-SN orbit will either be moving
toward pericenter ($\bm{v}\cdot\bm{r_{0}}<0$) or moving away from pericenter
($\bm{v}\cdot\bm{r_{0}}<0$). We can define an angle, $\phi_{0}$, which
measures how far along the hyperbola the two stars are immediately after core-
collapse, as measured from the focus:

$\phi_{0}=\begin{cases}\cos^{-1}\left[\left(p/r_{0}-1\right)/e\right],\text{for
}\bm{v}\cdot\bm{r_{0}}<0\\\
-\cos^{-1}\left[\left(p/r_{0}-1\right)/e\right],\text{for
}\bm{v}\cdot\bm{r_{0}}>0\end{cases}$ (13)

After some time, the binary’s separation becomes arbitrarily large, and we can
define the angle the BH leaves the system as

$\lim_{r_{0}\to\infty}\phi_{f}=\cos^{-1}\left(1/e\right).$ (14)

The runaway velocity of the two stars can be calculated by performing a
coordinate transformation, where we first apply a rotation $\mathcal{R}$ to
align the post-SN angular momentum vector to $\hat{z}$ then a clockwise
rotation $T$ to rotate the coordinate axis by $\phi_{f}-\phi_{0}$. In this
coordinate frame, the final orbital velocity is:

$\bm{v}^{\prime}=\sqrt{2\mathcal{E}}T(\phi_{f}-\phi_{0})\mathcal{R}\hat{r_{0}}.$
(15)

By applying the inverse transformation, we can now calculate the runaway
velocities of the two components in our initial inertial frame:

$\displaystyle\bm{v_{\rm BH,f}}$ $\displaystyle=$
$\displaystyle-\frac{M_{2}}{M_{\rm
BH}+M_{2}}\mathcal{R}^{-1}(\bm{v}^{\prime})+\bm{v_{\rm CM}}$ (16)
$\displaystyle\bm{v_{\rm 2,f}}$ $\displaystyle=$ $\displaystyle\frac{M_{\rm
BH}}{M_{\rm BH}+M_{2}}\mathcal{R}^{-1}(\bm{v}^{\prime})+\bm{v_{\rm CM}}.$ (17)

We have verified the validity of this procedure by comparing it with direct
orbital integration.

##  3 Results

Figure 1:  Distributions of runaway velocities for possible
MOA-11-191/OGLE-11-0462 progenitors for binaries with an initial eccentricity
of 0.1, adopting a BH mass of 7.1 $M_{\odot}$ from Sahu et al. (2022). Each
panel shows a different combination of possible He-star progenitor masses and
companion masses (left column shows low-mass companions and right column shows
high-mass companions), while the curves in each panel show the effects of
different kick velocities. In the top right of each panel we provide an inset
displaying the fraction of systems that disrupt during core collapse. For most
binary configurations, disruption velocities provide a close estimate of the
kick velocity the BH received at birth. Binaries with more massive companions
(right column of panels) have a somewhat broader distribution of runaway
velocities. However, even in these cases, the runaway velocities provides an
accurate estimate of the kick velocity.  Figure 2:  We provide the runaway
velocities for possible MOA-11-191/OGLE-11-0462 progenitors, for binaries with
an initial eccentricity of 0.9, again adopting a BH mass of 7.1 $M_{\odot}$.
Our results largely mimic those from Figure 1 except the distributions are
even narrower, suggesting the runaway velocity is an even closer
representation of the SN natal kick for highly eccentric pre-SN orbits. This
is because eccentric binaries spends most of their time near apocenter, where
the orbital velocities are small and the gravitational potential is shallow;
runaway velocities are therefore dominated by the natal kicks received.
Systems with high mass companions have distributions that peak toward somewhat
lower runaway velocities, suggesting that the runaway velocities of these BHs
may underestimate the natal kicks. Systems that are disrupted with lower kicks
produce a much larger range of possible runaway velocities; however the insets
in each panel suggest that these binaries rarely become disrupted during a SN.
Only in the bottom, left panel, representing low-mass companions and high-mass
He-star progenitors do low-kick systems tend to disrupt the binary with any
frequency. The black curve in this panel demonstrates that the runaway
velocity is an accurate representation of the BH’s natal kick.

We use a Monte Carlo method to determine the effects of different SN kicks and
mass loss on different possible pre-SN binary configurations. For each
combination of He-star mass, companion mass, pre-SN orbital eccentricity, and
kick velocity, we randomly and independently generate $10^{5}$ separate kick
directions (drawn from an isotropic distribution across all $4\pi$ steradians)
and eccentric anomalies (properly accounting for the fraction of an orbit an
eccentric binary spends at each $E$). After processing each of our randomly
drawn binaries using the method described in Section 2, we obtain a
distribution of runaway velocities for those binaries that are disrupted
during core-collapse.

Figure 3:  We provide the runaway velocities for possible
MOA-11-191/OGLE-11-0462 progenitors, for binaries with an initial eccentricity
of 0.1, again adopting a BH mass of 2.1 $M_{\odot}$ following Lam et al.
(2022). Our results largely mimic those from Figures 1 and 2; the BH runaway
velocities are an accurate proxy for the BH’s natal kick. The only exception
is in the bottom right panel, with large He-star masses ($\simeq$20
$M_{\odot}$) and large companion masses ($\simeq$20 $M_{\odot}$), in which
case the BH’s runaway velocity can be as much as twice the magnitude of its
natal kick. Nevertheless, the black bar in the bottom right panel’s inset
shows that BHs receiving such relatively low kick velocities only have a
$\simeq$40% of disrupting.

To explore the effects of different orbital parameters, we test somewhat
extreme orbital configurations to demonstrate how these affect our results and
with the idea that these will bracket the true progenitor to
MOA-11-191/OGLE-11-0462. We test both low-mass ($1\,{{M_{\odot}}}$) and high-
mass ($20\,{{M_{\odot}}}$) companions, as well as both low-mass
($8\,{{M_{\odot}}}$) and high-mass ($20\,{{M_{\odot}}}$) He-star progenitors
to the BH. To start, we adopt a mass of 7.1 $M_{\odot}$ for
MOA-11-191/OGLE-11-0462, following Sahu et al. (2022). Figure 1 shows the
resulting runaway velocities (as a fraction of the natal kick velocity) for
these combinations of masses for different non-dimensional kick velocities of
0.2, 1, 2, and 5 (in terms of $V_{c}$) for a binary with an initial
eccentricity of 0.1.

The left column of panels shows that a low-mass companion results in runaway
velocities within $\simeq$10% of the kick velocity. The lack of a black curve
in the top, left panel shows that when very little mass is lost from the
orbit, a kick velocity similar to or greater than the orbital velocity is
required to disrupt the binary. This is additionally expressed by the inset in
each panel of Figure 1, where we provide the fraction of systems disrupted
during SN for each SN kick velocity. The inset in the bottom, left panel shows
that when the He-star progenitor to the BH is sufficiently large
($\beta\lesssim$0.5), most binaries will disrupt, regardless of the natal kick
received by the BH. These binaries all have runaway velocities very similar to
the BH’s natal kick. Even in the case of $v_{k}=0.2$, the runaway velocities
are within 20% of the BH natal kick velocity.

In the right column of panels in Figure 1 we show the runaway velocities for
binaries with massive companions. These systems are more massive and therefore
require a larger kick magnitude ($v_{k}\gtrsim$1) to disrupt with any notable
frequency. The distribution of runaway velocities for these binaries peak near
the BH natal kick velocity, albeit with a somewhat wider distribution compared
with the left column of panels. The massive companion means that the He-star
progenitor to the BH is initially orbiting with a much larger velocity (with
respect to the system’s center of mass) which leads to a wider distribution of
possible runaway velocities, depending on the binary’s orbital phase at core
collapse and SN kick direction. Nevertheless, the runaway velocities tend to
be within 20% of the BH natal kick. For systems that receive moderate SN kicks
($v_{k}=$ 1 or 2), the distribution of possible runaway velocities extends to
somewhat lower values, implying the SN kick could be somewhat larger than the
observed runaway velocity.

Figure 4:  We provide the runaway velocities for possible
MOA-11-191/OGLE-11-0462 progenitors, for binaries with an initial eccentricity
of 0.9, again adopting a BH mass of 2.1 $M_{\odot}$ following Lam et al.
(2022). As was the case for Figures 1, 2, and 3, the BH runaway velocities are
an accurate proxy for the BH’s natal kick. The only exceptions (those with low
$v_{k}$ and large $M_{2}$) are unlikely to disrupt.

In Figure 2 we provide the same results as in Figure 1, but for binaries with
an initial eccentricity of 0.9. Comparison between the two figures shows
qualitatively similar results, with the largest difference being that the
highly eccentric binaries have distributions of runaway velocities even more
sharply peaked around the kick velocity. This is because highly eccentric
binaries spend most of their orbits far from pericenter, where the orbital
velocities and gravitational pull between the two components are low. In these
cases, the SN kick plays the significantly dominant role in the outcome of the
binary. A secondary effect can be seen in the $v_{k}=1,2$ curves in the right
column of panels in Figure 2; these distributions are peaked at runaway
velocities even smaller than the BH natal kick velocity, and extend to even
lower velocities than that seen in Figure 1 for binaries with initial
eccentricities of 0.1. If the true progenitor parameters of
MOA-11-191/OGLE-11-0462 are well represented by these orbital configurations,
then the runaway velocities in fact provide underestimates of the BH natal
kick velocity.

For comparison, we additionally provide the distributions of runaway
velocities adopting a BH mass of 2.1 $M_{\odot}$ for MOA-11-191/OGLE-11-0462,
following the Solution 2 from Lam et al. (2022). Figure 3 provides the
resulting distributions for an eccentricity of 0.1 while Figure 4 provides the
same distributions but for an eccentricity of 0.9. We have adjusted the helium
star mass, testing 6 $M_{\odot}$ and 20 $M_{\odot}$ models. In nearly all of
the curves shown in Figures 3 and 4 the runaway velocity is very similar to
the BH’s kick velocity, typically within 50%, often much closer. As was the
case with the distributions in Figures 1 and 2 for a 7.1 $M_{\odot}$ BH, the
only exceptions occur for our $v_{k}=0.2$ model with $M_{2}=20$ $M_{\odot}$,
in which case the runaway velocity can be as much as a factor of 2 different
from $v_{k}$. However, such cases ought to be rare as the black bars in the
insets of these panels indicate such binaries only have a $\lesssim$20%
probability of disrupting.

Figure 5:  The observed position, distance, and proper motion of
MOA-11-191/OGLE-11-0462 provide five of the six dimensions of phase space
describing its orbit. We calculate its orbit through the Milky Way adopting
different possible radial velocities. Columns show different 2-dimensional
projections of the resulting orbits, while rows show the solution from Sahu et
al. (2022) and the two solutions from Lam et al. (2022). The orbit is nearly
circular, but can extend to a large height above the Milky Way disk. Black
lines show Milky Way orbits of nearby field stars. MOA-11-191/OGLE-11-0462’s
orbit is inconsistent with the thin disk, but is representative of the
kinematic thick disk.

The curves in Figures 3 and 4 closely match those of Figures 1 and 2,
calculated for a 7.1 $M_{\odot}$ BH. The close similarities allow us to infer
that our results can be applied to models with 3.7 $M_{\odot}$ BHs following
Solution 1 from Lam et al. (2022) and indeed any reasonable mass for
MOA-11-191/OGLE-11-0462. Taken together, the results presented in Figures 1,
2, 3, and 4 lead us to conclude that the runaway velocity of BHs formed from a
disrupting system provide an accurate proxy for the kicks these objects
received at formation. While the correlation between the natal kick velocity
and runaway velocity was pointed out by Tauris & Takens (1998) and Kiel &
Hurley (2009), we show it to be robust for eccentric pre-SN orbits and
highlight how it can actually be used to constrain natal kicks. By measuring
the peculiar velocity of a BH or other compact object that formed from a
binary that disrupted during core collapse, we are in effect constraining the
natal kick that object received. However, as we discuss in the following
sections, placing accurate constraints on the BH’s natal kick requires
knowledge of the pre-SN binary’s peculiar velocity.

##  4 MOA-11-191/OGLE-11-0462’s Peculiar Velocity

In addition to the BH masses, the combined astrometric and photometric
microlensing signals analyzed by Sahu et al. (2022) and Lam et al. (2022)
provide distances and proper motions for MOA-11-191/OGLE-11-0462. Sahu et al.
(2022) report a distance of $1.58\pm 0.18$ kpc and a proper motion of
$\mu_{\alpha}=-4.36\pm$0.22 mas yr-1 and $\mu_{\delta}=3.06\pm$0.66 mas yr-1.
Lam et al. (2022) report two solutions. For Solution 1, Lam et al. (2022) find
a parallax of 0.65${}^{+0.08}_{-0.08}$ mas and a proper motion
($\mu_{\alpha}$, $\mu_{\delta}$) of ($-2.36^{+0.12}_{-0.13}$,
1.49${}^{+0.73}_{-0.67}$) mas yr-1. For Solution 2, Lam et al. (2022) find a
parallax of 1.09${}^{+0.36}_{-0.33}$ mas and a proper motion of
($-0.72^{+0.96}_{-0.94}$, 1.54${}^{+1.29}_{-1.24}$) mas yr-1.

Since Sahu et al. (2022) and Lam et al. (2022) use the same photometric and
astrometric data, the different astrometric solutions reported by these
authors have their origin in the details of the lensing models they fit. In
particular, both sets of authors report discrepancies between the photometric
and astrometric lensing signal, with the photometric data preferring a
heavier, faster moving compact object. Indeed, the two different solutions
reported by Lam et al. (2022) depend on how these authors comparatively weigh
the photometric and astrometric data. Throughout this work, we consider all
three solutions for MOA-11-191/OGLE-11-0462.

Figure 6:  The eccentricity of MOA-11-191/OGLE-11-0462’s orbit around the
Milky Way with the maximum height it reaches above the Galactic Plane as a
function of its radial velocity. For comparison we also show nearby field
stars (black points) and globular clusters (red squares). Regardless of its
radial velocity MOA-11-191/OGLE-11-0462 is consistent with being part of the
kinematic thick disk.

Along with the sky position, the derived distances (or parallaxes) and proper
motions provide five of the six dimensions of phase space for
MOA-11-191/OGLE-11-0462; the radial velocity is unknown. In Figure 5 we show
the possible orbits in the Milky Way for MOA-11-191/OGLE-11-0462, adopting
possible radial velocities of $-$100, 0, and $+$100 km s-1 for each of the
three possible solutions (different rows). The different columns show the
orbits from face-on (left column) and edge-on from either side (middle and
right columns). These orbits are produced by integrating
MOA-11-191/OGLE-11-0462 backwards for 200 Myr using gala (Price-Whelan et al.,
2021) within the multi-component Milky Way model from Bovy (2015), adopting
astropy v4.0 (Astropy Collaboration et al., 2013, 2018) constants for the
position and peculiar velocity of the Sun. For comparison in Figure 5 we
additionally show the orbits (black lines) for stars in the same vicinity of
MOA-11-191/OGLE-11-0462, identified by selecting all Gaia EDR3 (Gaia
Collaboration et al., 2016, 2021; Lindegren et al., 2021) stars within 2∘ of
MOA-11-191/OGLE-11-0462 that have a parallax divided by the parallax
uncertainty of at least 10, a parallax between 0.55 and 0.7 mas, and radial
velocities measured by Gaia DR2 (Cropper et al., 2018).

The first row of Figure 5, shows that variations in radial velocity by
$\pm$100 km s-1 largely effect the eccentricity of the Galactic orbit, but do
not strongly affect the maximum height above or below the Milky Way disk,
$Z_{\rm max}$. This is a result of the strategy that microlensing surveys
typically focus on the Galactic Bulge where stellar densities are high to
maximize the likelihood of observing an event (Mróz et al., 2019). This
perspective allows for a lens’s motion in the Galactic $Z$ direction to be
well measured. Fortunately, $Z_{\rm max}$ is also a reasonable indicator of
the Milky Way population to which a star belongs. All the orbits for
MOA-11-191/OGLE-11-0462 in Figure 5 are consistent with disk orbits,
regardless of which astrometric solution or radial velocity is adopted
(excluding extreme velocities with magnitudes $>>$100 km s-1). Comparison
against stars in the same vicinity shows that the $Z_{\rm max}$ value is
sensitive to which astrometric solution is adopted.

In Figure 6 we compare the orbital eccentricities and $Z_{\rm max}$ of nearby
stars against MOA-11-191/OGLE-11-0462 for each of the three astrometric
solutions, adopting a range of possible radial velocities. This figure
demonstrates that the kinematic thin disk is clearly separate from any of the
three solutions for MOA-11-191/OGLE-11-0462, although Solution 2 from Lam et
al. (2022) is close. However, MOA-11-191/OGLE-11-0462 is consistent with
orbits from the Milky Way thick disk. There are two possibilities: either
MOA-11-191/OGLE-11-0462 was born in the kinematic thick disk or it was formed
in the kinematic thin disk and was kicked into a thick-disk-like orbit. In the
former, constraints are difficult to place as thick disk stars span a wide
range of kinematics characteristics; the different radial velocity solutions
for MOA-11-191/OGLE-11-0462 in Figure 6 show that a thick disk star could
receive a kick as large as 100 km s-1 and still appear as a thick disk star.

Figure 7:  We compare orbital eccentricity (in the Milky Way) and maximum
height reached above the Galactic Plane for MOA-11-191/OGLE-11-0462 and for
populations of nearby stars after artificially applying different velocity
kicks to their trajectories. Contours enclose 90% of the samples, and are
therefore representative of the kinematic thin disk. If
MOA-11-191/OGLE-11-0462 were formed in the thin disk and subsequently kicked
out, it likely had to receive a kick velocity in excess of 50 km s-1.

Because the kinematic thin disk is dynamically colder, we can place
constraints on the potential kick MOA-11-191/OGLE-11-0462 received at birth.
To each star in our sample of stars in the vicinity of
MOA-11-191/OGLE-11-0462, we apply 103 random kicks, isotropically distributed
with magnitudes of 10, 25, 50, and 100 km s-1, then integrate their orbits to
find their orbital eccentricity and $Z_{\rm max}$. Based on our results
described in Section 3, these kicks correspond identically to the natal kicks
of single BHs as well as those of BHs formed out of the disruption of a
stellar binary. Contours in Figure 7 contain the 90% confidence interval of
the resulting distributions. Increasing kick velocities lead to both an
increase in the orbital eccentricity as well as $Z_{\rm max}$. Different
constraints can be placed on the peculiar velocity of MOA-11-191/OGLE-11-0462,
depending on which astrometric solution is adopted. For the solution from Sahu
et al. (2022), we find a kick velocity of 50-100 km s-1 is required.
Alternatively, Solutions 1 and 2 from Lam et al. (2022) allow kick velocities
of up to $\simeq$50 km s-1, with Solution 2 more consistent with lower kick
velocities. We stress that these constraints are statistical, based on a
qualitative comparison between MOA-11-191/OGLE-11-0462 and nearby stars in its
vicinity for only two parameters summarizing their orbit in the Galaxy.
Furthermore, they are under the assumption that MOA-11-191/OGLE-11-0462 formed
within the kinematic thin disk and has a radial velocity within $\simeq$100 km
s-1 of zero. In the case that MOA-11-191/OGLE-11-0462’s unknown radial
velocity is significantly non-zero (which is very unlikely based on viewing
angle statistics), the required kick can become arbitrarily large.

##  5 How did MOA-11-191/OGLE-11-0462 Form?

Here we consider multiple scenarios for how MOA-11-191/OGLE-11-0462 could have
formed.

###  5.1 Single Star Scenario

Although $>$70% of massive stars have binary companions (Sana et al., 2012)
and most isolated BHs have a binary origin (Vigna-Gómez & Ramirez-Ruiz, 2022),
we will first consider the possibility that MOA-11-191/OGLE-11-0462 was one of
the minority of massive stars that formed and evolved as a single star. Given
its Galactic orbit described in Figures 5 and 6, MOA-11-191/OGLE-11-0462
either formed in the kinematic thick disk and received a kick $\lesssim$100 km
s-1 or was formed in the thin disk and received a kick of at most a few tens
of km s-1 into its current orbit, with a value dependent upon which
astrometric solution is adopted. Since it is a BH, we have no constraints on
MOA-11-191/OGLE-11-0462’s age; BHs of a few $M_{\odot}$ can be formed at any
metallicity (e.g., Riley et al., 2022). The relatively likelihood of whether
MOA-11-191/OGLE-11-0462 was formed in the thick disk or kicked out of the thin
disk is determined only from an astrophysical prior based on the relative mass
of the two populations and probability that they produce a detectable
microlensing signal. However, we note that in the single star scenario, any
perturbation to MOA-11-191/OGLE-11-0462’s kinematics is due to the kick it
received during core collapse; no other process in the evolution of isolated
single stars can produce substantial kicks. We therefore conclude that, if
MOA-11-191/OGLE-11-0462 was formed as a single star, it likely received a
natal kick $\lesssim$100 km s-1.

Figure 8:  We show the initial masses and pre-SN masses of single stars that
produce 2.1, 3.7, and 7.1 $M_{\odot}$ BHs using COSMIC. The possible
progenitors for MOA-11-191/OGLE-11-0462 depend strongly on metallicity: The
ZAMS masses are within a broad range of 15 to 90 $M_{\odot}$, while the pre-SN
masses are constrained to be between 6 and 20 $M_{\odot}$. Binary interactions
could alter this range somewhat.

###  5.2 Ejection from a Globular Cluster

The dynamics at the centers of globular clusters, especially those that are
massive and core-collapsed, can easily expel BHs into the Milky Way field.
Typical dynamical encounters that expel stars will cause them to only barely
escape the system, resulting in orbits that are similar to its host globular
cluster. In Figure 6 we provide the orbital eccentricities and $Z_{\rm max}$
of known globular clusters from the catalog of Baumgardt & Vasiliev (2021) and
Vasiliev & Baumgardt (2021), using updated Gaia EDR3 astrometry (Lindegren et
al., 2021; Gaia Collaboration et al., 2021). Although some globular clusters
have orbits consistent with MOA-11-191/OGLE-11-0462, particularly the solution
from Sahu et al. (2022), most have halo-like orbits. We therefore consider a
globular cluster origin for MOA-11-191/OGLE-11-0462 unlikely.

###  5.3 Formation From a Disrupted Binary

We consider the possibility that MOA-11-191/OGLE-11-0462 formed initially in a
binary system that was disrupted during BH formation. One possibility is that
the core collapse event forming MOA-11-191/OGLE-11-0462 disrupted the binary
system. In Figure 8 we use the binary population synthesis code COSMIC with
default code options (see Breivik et al., 2020, and references therein) to
determine initial masses and pre-SN masses of stars that evolve into 2.1, 3.7,
and 7.1 $M_{\odot}$ BHs. Although metallicity has a strong effect on the stars
that could evolve into MOA-11-191/OGLE-11-0462, the pre-SN masses (bottom
panel) are constrained to be within 6 and 20 $M_{\odot}$. Even though binary
interactions could lower this mass range somewhat, Figure 8 suggests that the
runaway velocity distributions in Figures 1 and 2 are applicable to the
formation of MOA-11-191/OGLE-11-0462. Under this scenario,
MOA-11-191/OGLE-11-0462 was unlikely to have formed with a kick velocity
$\gtrsim$100 km s-1 regardless of whether it formed in the kinematic thin or
thick disk.

An alternative possibility is MOA-11-191/OGLE-11-0462’s progenitor was ejected
from a binary system when its more massive companion collapsed, and
MOA-11-191/OGLE-11-0462 itself collapsed into a BH somewhat later. From the
top panel of Figure 8 we can derive possible binary configurations. At
$Z_{\odot}$, the primary in the system would have a pre-SN mass $\gtrsim$10
$M_{\odot}$, while MOA-11-191/OGLE-11-0462’s progenitor would have been
$\gtrsim$40 $M_{\odot}$. In this scenario, MOA-11-191/OGLE-11-0462’s
progenitor is the more massive of the two stars, suggesting its orbital
velocity compared to the system’s center of mass (and therefore its runaway
velocity upon the binary’s disruption; Tauris & Takens, 1998) will be small.
MOA-11-191/OGLE-11-0462’s peculiar velocity is therefore likely to be
dominated by whatever kick it received upon its own core collapse. At lower
metallicity the primary star could have been more massive at its core
collapse, and MOA-11-191/OGLE-11-0462’s progenitor could have received a
significant kick were the binary in a relatively tight orbit upon disruption.
MOA-11-191/OGLE-11-0462 would then receive an additional kick when its own BH
collapsed somewhat later.

One last possibility is that MOA-11-191/OGLE-11-0462 was the secondary in a
binary that survived its first SN, but was disrupted in the second SN.
Although two SNe are required to form double compact objects, it is unlikely
that MOA-11-191/OGLE-11-0462 formed from the second SN in such a system. Renzo
et al. (2019) find $\simeq$86% of massive binaries disrupt during the
primary’s core collapse. The surviving systems must avoid a merger during
subsequent mass transfer sequences to remain bound by the time the secondary
undergoes core collapse (Gallegos-Garcia et al., 2021). We therefore focus on
scenarios in which the system formed from the disruption of the first SN in a
binary.

##  6 Discussion and Conclusions

We calculate the runaway velocities of compact objects ejected from disrupted
binaries. While our analysis focuses on BHs, as we are primarily interested in
its implications for BH microlenses, this approach can also be applied to
runaway neutron stars, a topic already the subject of numerous studies (Gott
et al., 1970; Radhakrishnan & Shukre, 1985; Tauris & Takens, 1998). While the
connection between a compact object’s natal kick and its runaway velocity was
made by both Tauris & Takens (1998) and Kiel & Hurley (2009), it was not
recognized as a way to constrain natal kicks. Tauris & Takens (1998) showed
that the runaway velocity is principally dependent upon the natal kick
velocity, while other parameters such as the component masses have little
effect. By using Monte Carlo random draws Kiel & Hurley (2009) show that the
distribution of resultant runaway velocities are similar to the input
Maxwellian distribution of natal kick velocities. In this work, we focus on
individual systems, showing that the runaway velocity is typically within
$\simeq$20% of the compact object’s natal kick.

Of course, for sufficiently tight orbits or small natal kicks such that the
kick velocity is a small fraction of the orbital velocity, the correlation
between runaway velocity and natal kick velocity decreases. Nevertheless, in
our Monte Carlo simulations, even small kicks ($v_{\rm k}=0.2$) demonstrate a
significant correlation. How is this possible? Kicks significantly smaller
than the orbital velocity typically do not disrupt the binary unless a
substantial amount of mass is lost from the binary. For eccentric binaries
that spend most of their orbits far from periastron, more than half of the
system’s mass must be lost to disrupt (Hills, 1983). At the same time, the
orbital velocities far from periastron are lower than their circular velocity
equivalent. Therefore, systems with initially more eccentric orbits will
produce the tightest correlation between runaway velocities and natal kicks.

To determine the eccentricity distribution of pre-SN orbits, binary population
synthesis including an accurate distribution of binaries at the zero-age main
sequence is required. For sufficiently close orbits, mass transfer and tidal
interactions ought to circularize binaries prior to core collapse (van den
Heuvel, 1976). Yet if we are interested in the distribution of runaway BHs,
orbits of all periods, including those with wide orbits that may never
circularize prior to core collapse, ought to be considered. The tests we run
with different mass combinations, kick velocities, and eccentricities suggest
that our main result holds over a wide range of orbits: measuring the peculiar
velocity of a compact object provides a tight constraint of the kick velocity
it received at birth.

We apply this result to the recently detected BH MOA-11-191/OGLE-11-0462,
which has three solutions for mass, distance, and proper motion, each of which
we separately analyze. Although the lens’s radial velocity is unmeasured, we
find that for reasonable values, MOA-11-191/OGLE-11-0462’s motion through the
Milky Way is consistent with being part of the kinematic thick disk.
Alternatively were its progenitor formed in the kinematic thin disk, our
analysis suggests MOA-11-191/OGLE-11-0462 received a kick $\lesssim$100 km
s-1, depending on which solution is adopted.

While the natal kick constraints we place on MOA-11-191/OGLE-11-0462 are,
perhaps, quite broad, the method highlights the potential for constraints
derived from future, individual well-observed systems or even populations.
Recently, Vigna-Gómez & Ramirez-Ruiz (2022) use binary population synthesis to
analyze the formation of MOA-11-191/OGLE-11-0462 and suggest it – and other BH
microlenses like it – can be used to constrain binary evolution parameters in
general. In this work, we focus only on BH natal kicks, arguing that they are
especially constrained by astrometric microlenses with well-measured proper
motions.

When considering future observations, two competing effects lead to detection
biases that ought to be accounted for. Faster moving BHs are more likely to be
detected, as the likelihood of lensing a background star scales with its
tangential velocity. At the same time, degeneracy in photometric microlensing
solutions implies that lenses with longer timescales (and are therefore likely
to be massive) are more likely to be followed up with astrometric
observations. Regardless, the addition of future astrometric BHs with precise
proper motion and parallax measurements will only improve the constraints on
BH natal kicks.

We thank the anonymous referee for a suggestions which improved the quality of
this manuscript. We thank Jessica Lu and Casey Lam for useful discussions.
J.J.A. acknowledges support from CIERA and Northwestern University through a
Postdoctoral Fellowship. V.K. was partially supported through a CIFAR Senior
Fellowship and a Guggenheim Fellowship. This work has made use of data from
the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, Physical Review X, 11, 021053
  * Aghakhanloo et al. (2022) Aghakhanloo, M., Smith, N., Andrews, J., et al. 2022, arXiv e-prints, arXiv:2202.06887
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Baumgardt & Vasiliev (2021) Baumgardt, H., & Vasiliev, E. 2021, MNRAS, 505, 5957
  * Belczynski et al. (2008) Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223
  * Bennett et al. (2002) Bennett, D. P., Becker, A. C., Quinn, J. L., et al. 2002, ApJ, 579, 639
  * Blaauw (1961) Blaauw, A. 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Boubert & Evans (2018) Boubert, D., & Evans, N. W. 2018, MNRAS, 477, 5261
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29
  * Breivik et al. (2020) Breivik, K., Coughlin, S., Zevin, M., et al. 2020, ApJ, 898, 71
  * Callister et al. (2021) Callister, T. A., Farr, W. M., & Renzo, M. 2021, ApJ, 920, 157
  * Cropper et al. (2018) Cropper, M., Katz, D., Sartoretti, P., et al. 2018, A&A, 616, A5
  * De Donder et al. (1997) De Donder, E., Vanbeveren, D., & van Bever, J. 1997, A&A, 318, 812
  * Dray et al. (2005) Dray, L. M., Dale, J. E., Beer, M. E., Napiwotzki, R., & King, A. R. 2005, MNRAS, 364, 59
  * Eldridge et al. (2011) Eldridge, J. J., Langer, N., & Tout, C. A. 2011, MNRAS, 414, 3501
  * Fragos et al. (2009) Fragos, T., Willems, B., Kalogera, V., et al. 2009, ApJ, 697, 1057
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gallegos-Garcia et al. (2021) Gallegos-Garcia, M., Berry, C. P. L., Marchant, P., & Kalogera, V. 2021, ApJ, 922, 110
  * Gerke et al. (2015) Gerke, J. R., Kochanek, C. S., & Stanek, K. Z. 2015, MNRAS, 450, 3289
  * Gott et al. (1970) Gott, J. Richard, I., Gunn, J. E., & Ostriker, J. P. 1970, ApJ, 160, L91
  * Hills (1970) Hills, J. G. 1970, Nature, 226, 730
  * Hills (1983) —. 1983, ApJ, 267, 322
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897
  * Janka (2013) Janka, H.-T. 2013, MNRAS, 434, 1355
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, , .  http://www.scipy.org/
  * Jonker & Nelemans (2004) Jonker, P. G., & Nelemans, G. 2004, MNRAS, 354, 355
  * Kiel & Hurley (2009) Kiel, P. D., & Hurley, J. R. 2009, MNRAS, 395, 2326
  * Kochanek et al. (2008) Kochanek, C. S., Beacom, J. F., Kistler, M. D., et al. 2008, ApJ, 684, 1336
  * Lam et al. (2022) Lam, C. Y., Lu, J. R., Udalski, A., et al. 2022, arXiv e-prints, arXiv:2202.01903
  * Langer (2012) Langer, N. 2012, ARA&A, 50, 107
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2
  * Lovegrove & Woosley (2013) Lovegrove, E., & Woosley, S. E. 2013, ApJ, 769, 109
  * Lyne & Lorimer (1994) Lyne, A. G., & Lorimer, D. R. 1994, Nature, 369, 127
  * Mandel (2016) Mandel, I. 2016, MNRAS, 456, 578
  * Mao et al. (2002) Mao, S., Smith, M. C., Woźniak, P., et al. 2002, MNRAS, 329, 349
  * Mirabel et al. (2001) Mirabel, I. F., Dhawan, V., Mignani, R. P., Rodrigues, I., & Guglielmetti, F. 2001, Nature, 413, 139
  * Mirabel & Rodrigues (2003) Mirabel, I. F., & Rodrigues, I. 2003, Science, 300, 1119
  * Mróz et al. (2019) Mróz, P., Udalski, A., Skowron, J., et al. 2019, ApJS, 244, 29
  * Patton & Sukhbold (2020) Patton, R. A., & Sukhbold, T. 2020, MNRAS, 499, 2803
  * Pijloo et al. (2012) Pijloo, J. T., Caputo, D. P., & Portegies Zwart, S. F. 2012, MNRAS, 424, 2914
  * Price-Whelan et al. (2021) Price-Whelan, A., Sipőcz, B., Lenz, D., et al. 2021, adrn/gala: v1.4.1, vv1.4.1, Zenodo, doi:10.5281/zenodo.5057630.  https://doi.org/10.5281/zenodo.5057630
  * Radhakrishnan & Shukre (1985) Radhakrishnan, V., & Shukre, C. S. 1985, in Supernovae, their Progenitors and Remnants, ed. G. Srinivasan & V. Radhakrishnan, 155
  * Renzo et al. (2019) Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, A&A, 624, A66
  * Repetto et al. (2012) Repetto, S., Davies, M. B., & Sigurdsson, S. 2012, MNRAS, 425, 2799
  * Repetto et al. (2017) Repetto, S., Igoshev, A. P., & Nelemans, G. 2017, MNRAS, 467, 298
  * Repetto & Nelemans (2015) Repetto, S., & Nelemans, G. 2015, MNRAS, 453, 3341
  * Reynolds et al. (2015) Reynolds, T. M., Fraser, M., & Gilmore, G. 2015, MNRAS, 453, 2885
  * Riley et al. (2022) Riley, J., Agrawal, P., Barrett, J. W., et al. 2022, ApJS, 258, 34
  * Sahu et al. (2022) Sahu, K. C., Anderson, J., Casertano, S., et al. 2022, arXiv e-prints, arXiv:2201.13296
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444
  * Shklovskii (1970) Shklovskii, I. S. 1970, Soviet Ast., 13, 562
  * Smartt (2009) Smartt, S. J. 2009, ARA&A, 47, 63
  * Stevenson (2022) Stevenson, S. 2022, arXiv e-prints, arXiv:2202.03584
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38
  * Tauris & Takens (1998) Tauris, T. M., & Takens, R. J. 1998, A&A, 330, 1047
  * Tauris et al. (2017) Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021, arXiv e-prints, arXiv:2111.03606
  * van den Heuvel (1976) van den Heuvel, E. P. J. 1976, in Structure and Evolution of Close Binary Systems, ed. P. Eggleton, S. Mitton, & J. Whelan, Vol. 73, 35
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22
  * van Rensbergen et al. (1996) van Rensbergen, W., Vanbeveren, D., & De Loore, C. 1996, A&A, 305, 825
  * Vasiliev & Baumgardt (2021) Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978
  * Vigna-Gómez & Ramirez-Ruiz (2022) Vigna-Gómez, A., & Ramirez-Ruiz, E. 2022, arXiv e-prints, arXiv:2203.08478
  * Willems et al. (2005) Willems, B., Henninger, M., Levin, T., et al. 2005, ApJ, 625, 324
  * Wong et al. (2014) Wong, T.-W., Valsecchi, F., Ansari, A., et al. 2014, ApJ, 790, 119
  * Wong et al. (2012) Wong, T.-W., Valsecchi, F., Fragos, T., & Kalogera, V. 2012, ApJ, 747, 111
  * Wong et al. (2010) Wong, T.-W., Willems, B., & Kalogera, V. 2010, ApJ, 721, 1689
  * Wongwathanarat et al. (2013) Wongwathanarat, A., Janka, H. T., & Müller, E. 2013, A&A, 552, A126
  * Woosley & Weaver (1986) Woosley, S. E., & Weaver, T. A. 1986, ARA&A, 24, 205
  * Wyrzykowski & Mandel (2020) Wyrzykowski, Ł., & Mandel, I. 2020, A&A, 636, A20
  * Wyrzykowski et al. (2016) Wyrzykowski, Ł., Kostrzewa-Rutkowska, Z., Skowron, J., et al. 2016, MNRAS, 458, 3012


# Negative Flux Aggregation to Estimate Feature Attributions

Author Name Affiliation<EMAIL_ADDRESS>Xin Li1 Deng Pan1 Chengyin Li1 Yao
Qiang1&Dongxiao Zhu1 1Department of Computer Science, Wayne State University,
USA
{xinlee,pan.deng,cli,yao<EMAIL_ADDRESS>
###### Abstract

There are increasing demands for understanding deep neural networks’ (DNNs)
behavior spurred by growing security and/or transparency concerns. Due to
multi-layer nonlinearity of the deep neural network architectures, explaining
DNN predictions still remains as an open problem, preventing us from gaining a
deeper understanding of the mechanisms. To enhance the explainability of DNNs,
we estimate the input feature’s attributions to the prediction task using
divergence and flux. Inspired by the divergence theorem in vector analysis, we
develop a novel Negative Flux Aggregation (NeFLAG) formulation and an
efficient approximation algorithm to estimate attribution map. Unlike the
previous techniques, ours doesn’t rely on fitting a surrogate model nor need
any path integration of gradients. Both qualitative and quantitative
experiments demonstrate a superior performance of NeFLAG in generating more
faithful attribution maps than the competing methods.

##  1 Introduction

The growing demand for trustworthy AI in security- and safety-critic domains
has motivated developing new methods to explain DNN predicitons using image,
text and tabular data Chefer et al. (2021); Li et al. (2021); Pan et al.
(2020); Qiang et al. (2022). As noted in some pioneering works, e.g., Hooker
et al. (2019); Smilkov et al. (2017); Kapishnikov et al. (2021), faithful
explanation is indispensable for a DNN model to be trustworthy. However, it
remains to be challenging for human to understand a DNN’s predictions in terms
of its input features due to its black-box nature.

As such, the field of explainable machine learning has emerged and seen a wide
array of model explanation approaches. Among others, local approximation and
gradient based methods are the two major categories that are more intensively
researched. Local approximation methods mimic the local behavior of a black-
box model within a certain neighborhood using some simple interpretable
models, such as linear models and decision trees. However, they either require
an additional model training processes (e.g., LIMERibeiro et al. (2016),
SHAPLundberg and Lee (2017)), or rely on some customized propagation rules
(e.g., LRPBach et al. (2015), Deep Taylor DecompositionMontavon et al. (2017),
DeepLIFTShrikumar et al. (2017), DeepSHAPChen et al. (2021)). Gradient based
methods such as Saliency Map Simonyan et al. (2013), SmoothGrad Smilkov et al.
(2017), FullGrad Srinivas and Fleuret (2019), Integrated Gradient (IG) and its
variants Sundararajan et al. (2017); Hesse et al. (2021); Erion et al. (2021);
Pan et al. (2021); Kapishnikov et al. (2019, 2021) require neither surrogates
nor customized rules but must tackle unstable estimates of gradients w.r.t.
the given inputs. IG type of path integration based methods mitigate this
issue via a path integration for gradient smoothing, however, this also
introduces another degree of instability and noise sourced from arbitrary
selections of baselines or integration paths. Others (e.g., Saliency Map,
FullGrad) relax this requirement nevertheless can be vulnerable to the small
perturbations of the inputs due to its locality.

Ideally, we hope to avoid surrogates, special rules, arbitrary baselines
and/or path integrals in interpreting DNN’s prediction. In addition, since DNN
interpretation often works on per sample basis, efficient algorithms are also
crucial for a scalable DNN explanation technique. Here we examine the
prediction behavior of DNNs through the lens of divergence and flux in vector
analysis. We propose a novel Negative Flux Aggregation (NeFLAG) approach,
which reformulates gradient accumulation as divergence. By converting
divergence to gradient fluxes according to divergence theorem, NeFLAG
interprets the DNN prediction with an attribution map obtained by efficient
aggregation of the negative fluxes (see divergence theorem, flux, and
divergence in Preliminaries).

We summarize our contributions as follows: 1) To the best of our knowledge,
this is the first attempt to tackle the problem of explaining DNN model
prediction leveraging the concepts of gradient divergence and fluxes. 2) Our
NeFLAG technique eliminates the need for path integration via converting
divergence to flux estimation, opening a new avenue for gradient smoothing
techniques. 3) We propose an efficient approximation algorithm to enable a
scalable DNN model explanation technique. 4) We investigate the relationship
between flux estimation and Taylor approximation, bridging our method with
other local approximation methods. 5) Both qualitative and quantitative
experiment results demonstrate NeFLAG’s superior performance to the competing
methods in terms of interpretation quality.

##  2 Related Work

We categorize our approach as one of the local approximation methods. Other
examples include LIME Ribeiro et al. (2016), which fits a simple model (linear
model or decision trees etc.) in the neighborhood of a given input point, and
use this simple model as a local interpretation. Similarly, SHAP Lundberg and
Lee (2017) generalizes LIME via a Shapley value reweighing scheme that has
been proved to yield more consistent results. These methods, in general, enjoy
their own merit when the underlying model is completely black-box, i.e., no
gradient and model architecture information are available. However, in the
case of DNN model interpretation, we usually know both model architecture and
gradient information, so it is often beneficial to utilize this extra
information as we are essentially interpreting the DNN model itself rather
than a surrogate model. Saliency Map Simonyan et al. (2013) and FullGrad
Srinivas and Fleuret (2019) exploit the gradient information, yet remain
sensitive to small perturbations of the inputs. SmoothGrad Smilkov et al.
(2017) improves it by averaging attributions over multiple perturbed samples
of the input by adding Gaussian noise to the original input. As we will see in
Section of Preliminaries below, our NeFLAG method not only exploits the
smoothing gradient information to achieve robustness against perturbations,
but also eliminates the need for fitting an extra surrogate model.

Another line of research explains DNN predictions via accumulating/integrating
gradients along specific paths from baselines to the given input. Here we
denote them as path integration based methods, represented by Integrated
Gradients (IG) Sundararajan et al. (2017). Recent variants includes fast IG
Hesse et al. (2021), and IG with alternative paths and baselines (e.g.,
Adversarial Gradient Integration (AGI) Pan et al. (2021), Expected Gradient
(EG) Erion et al. (2021)), Guided Integrated Gradients (GIG)Kapishnikov et al.
(2021), Attribution through Regions (XRAI) Kapishnikov et al. (2019)). IG
chooses a baseline (usually a black image) as the reference to calculate the
attribution by accumulating gradients along a straight-line path from the
baseline to the given input in the input space. AGI, on the other hand,
relaxes the baseline and straight-line assumption by utilizing adversarial
attacks to automatically generate both baselines and paths, thereafter
accumulating gradients along these paths. Both would require paths and
baselines for gradient smoothing via path integration regardless of manually
picked or generated, which could introduce attribution noises along the path
(as noted by Sundararajan et al. (2017), different paths could result in
completely different attribution maps). On the contrary, our NeFLAG method
doesn’t need a path nor a baseline, in which the gradient smoothing is
controlled by a single radius parameter $\epsilon$, opening up a new direction
for gradient smoothing techniques.

##  3 Preliminaries: Divergence and Flux

We start explaining our approach by introducing the concept of divergence and
flux in vector analysis. Lets first consider a general scenario: to interpret
a DNN’s prediction, we hope to characterize its local behavior. Let’s define a
DNN model $f:\boldsymbol{X}\to\boldsymbol{Y}$, which takes inputs
$\boldsymbol{x}\in\boldsymbol{X}$, and outputs
$f(\boldsymbol{x})\in\boldsymbol{Y}$. For simplicity, we also assume that the
model is locally continuously differentiable. When we query the interpretation
for a given input $\boldsymbol{x}$, we are interested in how and why a
decision is made, i.e., what is the underlying decision boundary. In fact,
this is also the idea behind many other interpretation methods such as
LIMERibeiro et al. (2016) and SHAPLundberg and Lee (2017). In these methods,
an interpretable linear model (or other simple models) is fitted via sampling
additional points around the neighborhood of input of interest. Clearly,
taking advantage of a certain kind of neighborhood aggregation is a promising
route for interpretation of local approximation.

When the gradients $\nabla_{\boldsymbol{x}}f$ are available, it is already a
decent indicator of DNN’s local behavior. If we only calculate the gradients
at $\boldsymbol{x}$, the resulting attribution map is called the Saliency
MapSimonyan et al. (2013). However, without aggregation, these gradients are
usually unstable due to adversarial effect, where a small perturbation of the
inputs can lead to large variation of gradient values. On the other hand,
these gradients may be vanishing due to the so-called saturation effect
Miglani et al. (2020). To overcome the instability, lets denote a neighborhood
around input $\boldsymbol{x}$ to be $V_{\boldsymbol{x}}$, and estimate the
average gradients over it. Intuitively, the resulting vector
$\int_{V_{\boldsymbol{x}}}\nabla f\cdot\mathrm{d}V_{\boldsymbol{x}}$ can be
viewed as a local approximation for the underlying model. However, neither a
single gradient evaluation nor neighborhood gradient integration exploits the
fact that the function value $f(\boldsymbol{x})$ can be viewed as a result of
gradient field flow accumulation. The gradient field flow therefore inspires
us to accumulate the gradients along the flow direction, giving rise to the
new idea of gradient accumulation.

###  3.1 Gradient Accumulation

Plenty of methods have already embraced the idea of accumulation without
explicitly defining it. They typically accumulate the gradients along a
certain path, for example, Integrated Gradients (IG)Sundararajan et al. (2017)
integrates the gradients along a straightline in the input space from a
baseline to the given input. Specifically, they study not a single input
point, but the gradient flows from a baseline point towards the given input
(Figure 1). The baseline here is selected assuming no information is contained
for decision. Another example is Adversarial Gradient Integration (AGI) Pan et
al. (2021), which instead integrates along multiple paths generated by an
adversarial attack algorithm, and aggregates all of them. This method also
accumulates gradient flows, and it differs from IG only in accumulation paths
and selection of baselines.

Figure 1: In a non-uniform field, IG accumulates gradient flow from a baseline
to the given input. Here black solid arrows denote the gradient field, the red
dashed arrow denotes the accumulation direction, and the length of green solid
arrows denote the gradient magnitude on the accumulation direction.

Although both IG and AGI exploit the idea of gradient accumulation, they only
tackle one-dimensional accumulation, i.e., accumulating over a path/paths. In
fact, Pan et al. (2021) point out that neither an accumulation path is unique,
nor a single path is necessarily sufficient. Ideally, an accumulation method
should take into account of every possible baseline and path. However, it is
nearly impossible to accumulate all possible paths. Then how could we overcome
such an obstacle? A promising direction is to examine the problem through the
lens of divergence and flux.

###  3.2 Divergence and Flux

For studying accumulation, the definition of divergence perfectly fits our
needs. Let $\mathbf{F}=\nabla_{x}f$ be the gradient field, the divergence
$\operatorname{div}\mathbf{F}$ is defined by

$\displaystyle\operatorname{div}\mathbf{F}=\nabla\cdot\mathbf{F}=\left(\frac{\partial}{\partial
x_{1}},\frac{\partial}{\partial
x_{2}},...\right)\cdot\left(F_{x_{1}},F_{x_{2}},...\right),$ (1)

where $F_{x_{i}}$ and $x_{i}$ are the gradient vector’s and input’s $i$th
entry, respectively. The above definition is convenient for computation,
however, difficult to interpret. A more intuitive definition of divergence can
be described as follows.

######  Definition 1.

Given that $\mathbf{F}$ is a vector field, the divergence at position
$\boldsymbol{x}$ is defined by the total vector flux through an infinitesimal
surface enclosure $S$ that encloses $\boldsymbol{x}$, i.e.,

$\displaystyle\left.\operatorname{div}\mathbf{F}\right|_{\mathbf{x}}=\lim_{V\rightarrow
0}\frac{1}{|V|}\oiint_{S(V)}\mathbf{F}\cdot\hat{\mathbf{n}}dS,$ (2)

where $V$ is an infinitesimal volume around $\boldsymbol{x}$ that is enclosed
by $S$, and $\mathbf{F}\cdot\hat{\mathbf{n}}$ denotes the normal vector flow
(flux) through the surface $S$.

It is not hard to see that $\nabla\cdot\mathbf{F}$ is essentially the
(negatively) accumulated gradients at point $\boldsymbol{x}$ from Definition
1, as it defines the total gradient flow contained in the infinitesimal volume
$V$.

Therefore, given input $\boldsymbol{x}$, to interpret the model within its
neighborhood denoted by $V_{\boldsymbol{x}}$, we only need to calculate the
accumulated gradients by integrating the divergence over all points within
this neighborhood. The attribution heatmap can be obtained by simply replacing
the dot product in Eq. 1 by an element-wise product, and then integrate, i.e.,

$\displaystyle\text{Attribution}=\int_{V_{\boldsymbol{x}}}\left(\frac{\partial}{\partial
x_{1}},\frac{\partial}{\partial
x_{2}},...\right)\odot\left(F_{x_{1}},F_{x_{2}},...\right)dV.$ (3)

Eqs. 2 and 3 are equivalent because of the divergence theorem described below.
We use element-wise product because we want to disentangle the effects from
different input entries.

However, there is one major obstacle, that is, how to integrate the whole
volume $V_{\boldsymbol{x}}$ especially when the input dimension is high. It is
tempting to sample a few points inside the neighborhood surface enclosure, and
sum up the divergence. But the computational cost for divergence, which
involves second order gradient computation, is much higher than only
calculating the first order gradients.

Fortunately, we can gracefully convert the volume integration of divergence
into surface integration of gradient fluxes using divergence theorem, and thus
simplify the computation.

###  3.3 Divergence Theorem

In vector analysis, the divergence theorem states that the total divergence
within an enclosed surface enclosure is equal to the surface integral of the
vector field’s flux over such surface enclosure. The flux is defined as the
vector field’s normal component to the surface. Formally, let $\mathbf{F}$ be
a vector field in a space $U$, $S$ is defined as a surface enclosure in such
space, and we call the neighborhood volume that is enclosed in $S$ as $V$.
Assuming that $\mathbf{F}$ is continuously differentiable on $V$, the
_Divergence Theorem_ states that

$\displaystyle\iiint_{V}(\nabla\cdot\mathbf{F})\mathrm{d}V=\oiint_{S}(\mathbf{F}\cdot\hat{\mathbf{n}})\mathrm{d}S.$
(4)

Here the integrating symbols $\iiint$ and $\oiint$ don’t necessarily need to
be $3$-dimensional and $2$-dimensional, respectively. They can be any higher
dimension as long as the surface integration is one dimension lower than the
volume integration.

In DNNs, the gradient field $\mathbf{F}=\nabla_{\boldsymbol{x}}f$ w.r.t. the
inputs is exactly a vector field in the input space. Note that there could be
non-differentiable points (that could violate the continuously differentiable
condition due to some activation functions such as ReLU may not differentiable
at some point). Nevertheless we can safely assume that it is at least
continuously differentiable within a certain small neighborhood.

Figure 2: The fluxes on closed surfaces with different shapes. a. The total
flux of a closed sphere surface is zero when the field is uniform, and the
vector obtained by negative flux vector aggregation is coherent to the field
vector. Here red arrows and blue arrows denote positive and negative flux. The
green arrow bar denotes the direction of the aggregated vector. b. When the
closed surface is not symmetric with respect to the field, the resulting
aggregated vector is not in the same direction as the field.

##  4 Negative Flux Aggregation

In this Section, we describe our interpretation approach, Negative Flux
Aggregation (NeFLAG). We first explain the necessity of differentiating
positive and negative fluxes, and then describe an algorithm on how to find
negative fluxes on a $\epsilon$-sphere and aggregate them for interpretation.

###  4.1 Negative Fluxes Define a Linear Model

Interpreting via directly integrating all fluxes is an intriguing approach,
however, the positive and negative fluxes should be interpreted differently.
By convention, positive fluxes are pointed outward the enclosure surface, and
the negative fluxes are pointed inward (Figure 2a). Let $S$ be a surface
enclosure and $V$ is its corresponding volume, a positive flux can be viewed
as the gradient loss of $V$, and similarly a negative flux is the gradient
gain of $V$. It means the confidence gain or loss of the prediction of $x$.
For example, from $f(\boldsymbol{x})=0.99$ to $f(\boldsymbol{x})=0.8$
represents a confidence loss. The sum of negative fluxes here means a total
confidence loss if $\boldsymbol{x}$ moves out of the neighborhood. When
interpreting a DNN model prediction $f$, assuming $\boldsymbol{x}$ is a point
in the input space $U$, the prediction of it by $f$ is class $A$. To interpret
this prediction at a location in the input space (i.e., $\boldsymbol{x}$), we
need to find out that moving to which direction could make it less likely to
be $A$? Apparently, a good choice is to move in the direction of negative
gradient, i.e., with negative gradient fluxes. In contrast, the positive flux
can be interpreted as what makes it more likely to be class $A$. But since
$\boldsymbol{x}$ is already predicted as $A$, this information is less
informative than the negative fluxes in terms of interpretation. (An
additional intuition is that for a given input $\boldsymbol{x}$ we are more
interested in its direction towards decision boundary, instead of the other
side of it.)

Linear Case: It is not hard to see that in a linear vector field
$\mathbf{F}_{l}$ with underlying linear function $f_{l}$, the sum of fluxes
within any enclosure is 0 (Figure 2). Moreover, when the enclosure is an
$N$-dimensional sphere (where $N$ is the dimension of the space), we can sum
up over all vectors on the sphere surface that are contributing negative flux.
The obtained summed-up vector is then equivalent to the weight vector of the
underlying linear function $f_{l}$.

The observation of the linear case study inspires us to learn the weight
vector of a local linear function by summing up all negative fluxes on an
$N$-d sphere. Hence we say that the vector sum of all negative fluxes
interprets the local behavior of the original function via a local linear
approximation. Formally, we state Assumption 1 as follows.

######  Assumption 1.

Given model $f$ whose gradient field is $\mathbf{F}$, and sample
$\boldsymbol{x}$, let $S_{\boldsymbol{x}}$ be an $N$-d $\epsilon$-sphere that
centered at $\boldsymbol{x}$. The vector sum of all negative flux can be
obtained by

$\displaystyle\boldsymbol{w}=\oiint_{S_{\boldsymbol{x}}^{-}}\mathbf{F}\odot\hat{\mathbf{n}}dS_{\boldsymbol{x}},$
(5)

where $\odot$ denotes element-wise product, and $S_{\boldsymbol{x}}^{-}$ is a
set of point on the sphere where flux is negative. The linear model defined by
weight vector $\boldsymbol{w}$ is then a local linear interpretation to the
original model of input $\boldsymbol{x}$.

Note that this assumption is derived from the observation when
$y=a_{1}x_{1}+a_{2}x_{2}+...+a_{n}x_{n}$, so we interpret this linear model by
attributions $\mathbf{attr}=(a_{1}x_{1},a_{2}x_{2},...,a_{n}x_{n})$. The entry
$attr_{i}=a_{i}x_{i}$ represents the contribution of the $i$th attribute to
the final prediction. Note that the attribution can also be written as
$\mathbf{attr}=\boldsymbol{a}\odot\boldsymbol{x}$, hence the element-wise
product in the Assumption. It is also critical to set $S_{\boldsymbol{x}}$ as
a sphere in Assumption 1, any other shape introducing any degree of asymmetry
would cause disagreement between local linear approximation and the negative
flux aggregation, as shown in Figure 2b.

Interpretation: An interpretation can be achieved by plotting an attribution
map using the vector $\boldsymbol{w}$ from Eq. 5. The rationale is directly
from the correspondence between this vector and the underlying linear model.

###  4.2 An Approximation Algorithm

To calculate Eq. 5, we have to find out both the point
$\tilde{\boldsymbol{x}}$ with negative flux and its normal vector
$\hat{\mathbf{n}}$. For the latter, given a candidate point
$\tilde{\boldsymbol{x}}$ on the sphere surface $S_{\boldsymbol{x}}$, we can
replace $\hat{\mathbf{n}}$ by
$(\tilde{\boldsymbol{x}}-\boldsymbol{x})/|\tilde{\boldsymbol{x}}-\boldsymbol{x}|$
because $\boldsymbol{x}$ is the center of the sphere $S_{\boldsymbol{x}}$, and
the radius $\epsilon=|\tilde{\boldsymbol{x}}-\boldsymbol{x}|$. As for point
$\tilde{\boldsymbol{x}}$ with negative flux, a straightforward solution is to
randomly subsample a list of points on the $\epsilon$-sphere, then select
those with negative flux. However, this trick has no guarantee on how many
subsampling is sufficient, and also cannot guarantee that a negative flux
exists. Therefore, we need an approximation algorithm that can guide us to
those points with negative fluxes.

Since we are interested in interpreting the local behavior, meaning that the
radius of $\epsilon$-sphere should be sufficiently small. This provides us
with the possibility to approximate the gradient (flux). i.e.,

$\displaystyle\mathbf{F}(\tilde{\boldsymbol{x}})\cdot\hat{\mathbf{n}}\approx\frac{\left(f(\tilde{\boldsymbol{x}})-f(\boldsymbol{x})\right)}{\epsilon}.$
(6)

Since $\boldsymbol{x}$ is the center of $\epsilon$-sphere
$S_{\boldsymbol{x}}$, to make the flux (left hand side of Eq. 6) negative, we
only need to find $\tilde{\boldsymbol{x}}$ such that
$f(\tilde{\boldsymbol{x}})<f(\boldsymbol{x})$. Moreover, for interpreting a
classification task, it is usually the case that $f(\boldsymbol{x})$ is of a
high value (when a prediction is confident, the output probability would be
close to 1). Therefore, if we could find a random local minimum
$\tilde{\boldsymbol{x}}$ on the $\epsilon$-sphere $S_{\boldsymbol{x}}$, it
would most likely have $f(\tilde{\boldsymbol{x}})<f(\boldsymbol{x})$.

In order to obtain a local minimum starting with any arbitrary initialization,
we propose the following Lemma.

######  Lemma 1.

Let $\boldsymbol{x}$ be the center of sphere $S_{\boldsymbol{x}}$, and
$\boldsymbol{x}^{(0)}$ be a random point on the $\epsilon$-sphere
$S_{\boldsymbol{x}}$. We define the following recurrence

$\displaystyle\boldsymbol{x}^{(k)}=\boldsymbol{x}-\epsilon\cdot\frac{\mathbf{F}\left(\boldsymbol{x}^{(k-1)}\right)}{\left|\mathbf{F}\left(\boldsymbol{x}^{(k-1)}\right)\right|}.$
(7)

Let $V_{x}$ be the space enclosed by $S_{\boldsymbol{x}}$, assuming
$\mathbf{F}$ is continuous on $V_{x}$, then $f(\boldsymbol{x}^{(k)})$
converges to a local minima on $S_{\boldsymbol{x}}$ as $k\to\infty$.

###### Proof.

To find a local minima on the surface $S_{\boldsymbol{x}}$, given the current
candidate $\boldsymbol{x}^{(k-1)}$ and its corresponding gradient
$\mathbf{F}(\boldsymbol{x}^{(k-1)})$, the updating formula needs to take
$\boldsymbol{x}^{(k-1)}$ to $\boldsymbol{x}^{(k)}$ such that
$\boldsymbol{x}^{(k)}$ is on the opposite direction of the tangent component
of $\mathbf{F}(\boldsymbol{x}^{(k-1)})$ to the surface $S_{\boldsymbol{x}}$.
Assuming the updating step to be $\epsilon$, the updating rule is then

$\displaystyle\boldsymbol{x}^{(k)}=\boldsymbol{x}^{(k-1)}-\epsilon\cdot\left(\frac{\mathbf{F}\left(\boldsymbol{x}^{(k-1)}\right)}{\left|\mathbf{F}\left(\boldsymbol{x}^{(k-1)}\right)\right|}-\frac{\mathbf{F}\left(\boldsymbol{x}^{(k-1)}\right)\cdot\hat{\mathbf{n}}}{\left|\mathbf{F}\left(\boldsymbol{x}^{(k-1)}\right)\right|}\cdot\hat{\mathbf{n}}.\right).$
(8)

Note that the term in the parentheses is exactly the tangent direction of the
gradient (i.e., the direction of the gradient minus the normal component). We
can derive that

$\displaystyle\boldsymbol{x}$
$\displaystyle=\boldsymbol{x}^{(k-1)}+\epsilon\cdot\frac{\left(\boldsymbol{x}-\boldsymbol{x}^{(k-1)}\right)}{\epsilon}$
(9)
$\displaystyle\approx\boldsymbol{x}^{(k-1)}+\epsilon\cdot\frac{\mathbf{F}\left(\boldsymbol{x}^{(k-1)}\right)\cdot\hat{\mathbf{n}}}{\left|\mathbf{F}\left(\boldsymbol{x}^{(k-1)}\right)\right|}\cdot\hat{\mathbf{n}}.,$
(10)

where $\boldsymbol{x}$ is the center of the $\epsilon$-sphere
$S_{\boldsymbol{x}}$, substituting Eq. 10 into Eq 8, we then have the updating
rule in Eq. 7. ∎

Lemma 1 tells us that given any initial point $\boldsymbol{x}^{(0)}$ on the
sphere, we can always find a point with local minimal flux value. Ideally, by
seeding multiple initial points, following the recurrence in Eq. 7, we can
obtain a set of local minima points on $S_{\boldsymbol{x}}$. The integration
in Eq. 5 can then be approximated by summation of the gradient fluxes at these
local minima points. But this approximation is suboptimal because the updating
rule in Eq. 8 will cause bias towards those points with _minimal fluxes_ ,
instead of distributing evenly to points with _negative fluxes_. To overcome
this issue, we add a sign operation on the recurrence in Eq. 7, implemented in
Algorithm 1, to introduce additional stochasticity, hence making the negative
flux points distributed more evenly. Algorithm 1 describes this trick as well
as the step-by-step procedure of calculating NeFLAG.

Here we note Eq. 5 provides a formulation (or a framework) to calculate a
representing vector for the neighborhood near $\boldsymbol{x}$. We develop
Algorithm 1 as one way of approximation to efficiently perform the
experiments. There are other algorithms with better approximation accuracy
that we may explore the different approximation tricks in the future.

Algorithm 1 $\text{NeFLAG}(f,\boldsymbol{x},n,S_{\boldsymbol{x}},\epsilon,m$)

Input: $f$: Classifier, $\quad\boldsymbol{x}$: input, $\quad n$: number of
negative flux samples, $S_{\boldsymbol{x}}$: $\epsilon$-sphere,
$\quad\epsilon$: radius of $S_{\boldsymbol{x}}$, $\quad m$: max number of
backpropagation steps;
Output: Attribution map NeFLAG;

1: $\text{NeFLAG}\leftarrow 0$; $\quad j\leftarrow 0$; $\quad i\leftarrow 0$ ;

2: while $i=1:n$ do

3: while $j=1:m$ do

4: Randomly sample $\tilde{\boldsymbol{x}}$ on sphere $S_{\boldsymbol{x}}$;

5:
$\tilde{\boldsymbol{x}}=\boldsymbol{x}-\epsilon\cdot\text{sign}\left(\frac{\mathbf{F}\left(\tilde{\boldsymbol{x}}\right)}{\left|\mathbf{F}\left(\tilde{\boldsymbol{x}}\right)\right|}\right)$
;

6: end while

7:
$\text{NeFLAG}=\text{NeFLAG}+\mathbf{F}\left(\tilde{\boldsymbol{x}}\right)\odot(\boldsymbol{x}-\widetilde{\boldsymbol{x}})$;

8: end while

Figure 3: Examples of attribution maps obtained by NeFLAG, AGI and IG methods.
The underlying prediction model is InceptionV3 (additional examples for
ResNet152 and VGG19 can be found in the supplementary materials). Compared to
AGI and IG, we observe that NeFLAG’s attribution map has a clearer shape and
focuses more densely on the target object.

###  4.3 Connection to Taylor Approximation

NeFLAG can be viewed as a generalized version of Taylor approximation. Here we
show that NeFLAG can not only be viewed as a gradient accumulation method, but
also an extension to some local approximation methods. Given that
$S_{\boldsymbol{x}}$ is an $N$-d $\epsilon$-sphere centered at
$\boldsymbol{x}$, for any point $\tilde{\boldsymbol{x}}$ on the sphere
surface, we can represent the normal vector $\hat{\mathbf{n}}$ as
$(\tilde{\boldsymbol{x}}-\boldsymbol{x})/\epsilon$. This inspires us to
investigate its connection to the first order Taylor approximation. As pointed
out by Montavon et al. (2017), Taylor decomposition can be applied at a nearby
point $\tilde{x}$, then the first order decomposition can be written by

$\displaystyle\begin{aligned}
f(\boldsymbol{x})&=f(\widetilde{\boldsymbol{x}})+\left(\left.\nabla_{\boldsymbol{x}}f\right|_{\boldsymbol{x}=\widetilde{\boldsymbol{x}}}\right)^{\top}\cdot(\boldsymbol{x}-\widetilde{\boldsymbol{x}})+\eta,\end{aligned}$
(11)

where $\eta$ is the error term of second order and higher. We use
$R(\boldsymbol{x})=\left(\left.\nabla_{\boldsymbol{x}}f\right|_{\boldsymbol{x}=\widetilde{\boldsymbol{x}}}\right)\odot(\boldsymbol{x}-\widetilde{\boldsymbol{x}})$
to represent a heatmap generated by this approximation process. Note that the
heatmap $R(\boldsymbol{x})$ completely attributes $f(\boldsymbol{x})$ if
$f(\tilde{\boldsymbol{x}})=0$ and the error term $\eta$ can be omitted.

One may notice that if $\tilde{\boldsymbol{x}}$ is by chance on the N-d
$\epsilon$-sphere, the first order Taylor approxmation of $\boldsymbol{x}$
from $\tilde{\boldsymbol{x}}$ is then equivalent to the flux estimation at
$\tilde{\boldsymbol{x}}$. To explain it more clearly, we can simply view
$-(\boldsymbol{x}-\tilde{\boldsymbol{x}})/|\boldsymbol{x}-\tilde{\boldsymbol{x}}|$
as the normal vector $\hat{\mathbf{n}}$ in Eq. 5, and
$\left(\left.\nabla_{\boldsymbol{x}}f\right|_{\boldsymbol{x}=\widetilde{\boldsymbol{x}}}\right)$
is exactly the gradient vector $\mathbf{F}$ at $\tilde{\boldsymbol{x}}$. Hence
the first order Taylor decomposition and a single point flux estimation
differs only by a factor of $-\epsilon$ as
$\epsilon=|\boldsymbol{x}-\tilde{\boldsymbol{x}}|$.

Since we only care about the negative fluxes, then the second term of Eq 11,
i.e.,
$\left(\left.\nabla_{\boldsymbol{x}}f\right|_{\boldsymbol{x}=\widetilde{\boldsymbol{x}}}\right)^{\top}\cdot(\boldsymbol{x}-\widetilde{\boldsymbol{x}})$
must be positive. If we further omit the error term $\eta$, we must have
$f(\boldsymbol{x})=f(\tilde{\boldsymbol{x}})+\text{positive value}$. It tells
us that in the perspective of first order Taylor decomposition, the negative
flux indeed attributes to positive predictions, bridging the connection
between NeFLAG and other local approximation methods.

Metrics | Methods/Models | IG | SG | GIG | EG | AGI | NeFLAG-1 | NeFLAG-10 | NeFLAG-20
---|---|---|---|---|---|---|---|---|---
Deletion Score $\downarrow$ | VGG19 | 0.071 | 0.049 | 0.054 | 0.041 | 0.040 | 0.034 | 0.052 | 0.059
ResNet152 | 0.132 | 0.064 | 0.090 | 0.066 | 0.068 | 0.055 | 0.076 | 0.081
InceptionV3 | 0.122 | 0.057 | 0.096 | 0.049 | 0.059 | 0.048 | 0.066 | 0.068
Insertion Score $\uparrow$ | VGG19 | 0.223 | 0.269 | 0.304 | 0.338 | 0.401 | 0.416 | 0.521 | 0.535
ResNet152 | 0.332 | 0.382 | 0.437 | 0.447 | 0.480 | 0.485 | 0.568 | 0.578
InceptionV3 | 0.375 | 0.410 | 0.465 | 0.478 | 0.480 | 0.544 | 0.618 | 0.625
Table 1: Quantitative evaluation using deletion and insertion scores.
NeFLAG-1, NeFLAG-10 and NeFLAG-20 are NeFLAG methods with various number of
negative flux samples. Note NeFLAG-1 already outperforms all the competing
methods with the computational complexity grows linearly with the number of
backprogation steps just like the most efficient IG method.

##  5 Experiments

In this Section, we demonstrate and evaluate NeFLAG’s performance both
qualitatively and quantitatively using a diverse set of DNN architectures and
baseline interpretation methods.

###  5.1 Experiment Setup

Models: InceptionV3 Szegedy et al. (2015), ResNet152 He et al. (2015) and
VGG19 Simonyan and Zisserman (2014) are selected as the pre-trained DNN models
to be explained.

Dataset: The ImageNet Deng et al. (2009) dataset is used for all of our
experiments. ImageNet is a large and complex data set (compared with smaller
and simpler data sets such as Places, CUB200, or Flowers102) for us to better
demonstrate the key advantages of our approach.

Baseline Interpretation Methods: IG Sundararajan et al. (2017) and AGI Pan et
al. (2021) are selected as baselines for both qualitative and quantitative
interpretation performance comparison. We also include Expected Gradients
(EG)Erion et al. (2021), Guided Integrated Gradient Kapishnikov et al. (2021),
and SmoothGrad Smilkov et al. (2017) in quantitative comparison as they can be
viewed as smoothed versions of attribution methods. We use the default setting
provided by captum 111https://captum.ai/ for IG, EG and SmoothGrad methods.
For AGI, we adopt the default parameter settings reported in Pan et al.
(2021), i.e., step size $\epsilon=0.05$, number of false classes $n=20$. For
GIG, we also used the default settings provided in their latest official
implementation Kapishnikov et al. (2021). Here we focus on the comparison with
gradient based methods since non-gradient methods typically require a
surrogate and additional optimization processes.

Figure 4: Quantitative Performance Comparison in terms of Difference between
the Insertion and Deletion Scores.

###  5.2 Qualitative Evaluation

We first show some examples to demonstrate a better quality of the attribution
maps generated by NeFLAG than by other baselines. The NeFLAG is configured as
follows: $\epsilon$-sphere radius is set to $\epsilon=0.1$, and the number of
random negative flux point $\tilde{\boldsymbol{x}}$ is $n=20$. Figure 3 shows
the attribution maps generated by various interpretation methods for the
Inception V3 model (examples for ResNet152 and VGG19 models can be found in
the Appendix). Qualitatively speaking, we denote an attribution map is
faithful to model’s prediction if 1) it focuses on the target objects, and 2)
has clear and defined shapes w.r.t. class label instead of sparsely
distributed pixels. Based on this notion, it is clearly observed that NeFLAG
has a better attribution heatmap quality than AGI and IG.

A key observation from Figure 3 is that the NeFLAG attribution heatmap
provides a defined shape of the entity of the class whereas IG and AGI
attribution heatmaps do not reflect that ability. In fact, the latter heatmaps
are scattered and do not retain a definite shape. For example, in Figure 3, IG
and AGI’s heatmaps fail to delineate a definite and slender shape of the Rock
Python. The same notion happens with the Impala and Chimpanzee where the
attribution heatmaps of IG and AGI are scattered across the background and
class entity. We note that the results shown in Figure 3 are very typical (not
cheery-picked) in our experiments. We also provide more examples from each
class to demonstrate the superior quality of attribution maps in the Appendix.

We note that the striking performance disparities referenced above could
mainly be caused by the inconsistency of path accumulation methods. Since both
IG and AGI incorporate a certain kind of path integration, any gradient vector
on that path is taken into account of the final interpretation. However, the
gradient vectors on the path aren’t necessarily unique and consistent. The
rationale behind this argument is that neither IG nor AGI is capable of
guaranteeing that their choice of accumulation path is optimal. For example,
in IG, the path is a straight-line, which isn’t by any means ideal. Similarly
in AGI, the path of an adversarial attack is chosen. However, due to the
stochastic effect, a slightly different initialization could result in
completely different attack paths. Our method NeFLAG, on the other hand, takes
advantage of negative flux aggregation without the burden of path
accumulation.

In fact, when we set the number of negative flux points $n=1$, NeFLAG reduces
to the IG method with a random baseline. In this special case, the only
difference is that NeFLAG with $n=1$ uses single point gradient evaluation
instead of a path integration in IG. By comparing NeFLAG and IG in Figure 3,
we can see that by simply replacing a path accumulation in IG ($n=1$, step
size $m=100$) with a surface flux integration in NeFLAG ($n=20$, step size
$m=1$), the interpretation becomes more faithful.

###  5.3 Quantitative Evaluation

In the quantitative experiments, we compare the performance of IG, AGI, EG,
GIG, SmoothGrad and three NeFLAGs variants with different numbers of sampled
negative flux points. Since our NeFLAG is considered as a robust local
attribution method, the size the neighborhood (radius of $\epsilon$-sphere)
and number of negative fluxes are tuning parameters. Similarly, the baseline
local attribution method, SmoothGrad, also has two tuning parameters: the
noise level and the number of samples to average over. Contrastively, IG, AGI,
EG, and GIG are global attribution methods so there is no tuning parameter
involved. We extensively investigate the neighborhood size and experiment with
various number of negative fluxes (e.g., NeGLAG-1, NeGLAG-10, and NeGLAG-20)
to faithfully demonstrate the stability of our method in comparison with
others.

Here VGG19, InceptionV3, and ResNet152 are used as underlying DNN prediction
models. Since we focus on explaining DNN prediction on a per sample basis just
like other DNN explanation methods, we randomly select 5,000 samples from
ImageNet dataset with 5 samples from each class as a good representation of
the classes. We use insertion score and deletion score as our evaluation
metricsPetsiuk et al. (2018). We replace the top pixels with black pixels in
the first round Petsiuk et al. (2018) whereas with Gaussian blurred pixels in
the second round Sturmfels et al. (2020). We report the average performance of
two rounds of experiments. Experimental details can be found in the Appendix.

Table 1 demonstrates that NeFLAG outperforms other baselines even with only a
single negative flux point (NeFLAG-1), and becomes much better when we
incorporate sufficient amount of negative flux points (NeFLAG-10 and
NeFLAG-20). In Figure 4, we systematically compared the performance of our
NeFLAG with the competing methods on three DNN architectures. We used the
Difference between the Insertion and Deletion scores as the metric for
comparison Shah et al. (2021). It is because the deletion and insertion scores
can be influenced by distribution shift caused by removing/adding pixies
Hooker et al. (2019). Since the shift occurs during both pixel insertion and
deletion, focusing on their relative difference instead of their absolute
values helps in neutralizing the effects of distribution shift. It is clearly
observed from Figure 4 that our NeFLAG method outperforms all the competing
methods across the three DNN models in terms of the Difference between the
Insertion and Deletion scores. We note there still a room for more reliable
and comprehensive evaluation metrics for the attribution methods.

##  6 Conclusion

We develop a novel DNN model explanation technique NeFLAG built off the
concept of divergence, and point out a connection to other local approximation
methods. NeFLAG doesn’t require a baseline, nor an integration path. This is
achieved by converting a volume integration of the second order gradients to a
surface integration of the first order gradients using divergence theorem.
Both qualitative and quantitative experiments demonstrate a superior
performance of NeFLAG in explaining DNN predictions over the strong baselines.

## References

  * Bach et al. [2015] Sebastian Bach, Alexander Binder, Grégoire Montavon, Frederick Klauschen, Klaus-Robert Müller, and Wojciech Samek.  On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation.  PloS one, 10(7), 2015.
  * Chefer et al. [2021] Hila Chefer, Shir Gur, and Lior Wolf.  Transformer interpretability beyond attention visualization.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 782–791, 2021.
  * Chen et al. [2021] Hugh Chen, Scott M Lundberg, and Su-In Lee.  Explaining a series of models by propagating shapley values.  arXiv preprint arXiv:2105.00108, 2021.
  * Deng et al. [2009] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei.  ImageNet: A Large-Scale Hierarchical Image Database.  In CVPR09, 2009.
  * Erion et al. [2021] Gabriel Erion, Joseph D Janizek, Pascal Sturmfels, Scott M Lundberg, and Su-In Lee.  Improving performance of deep learning models with axiomatic attribution priors and expected gradients.  Nature Machine Intelligence, pages 1–12, 2021.
  * He et al. [2015] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition. corr abs/1512.03385 (2015), 2015.
  * Hesse et al. [2021] Robin Hesse, Simone Schaub-Meyer, and Stefan Roth.  Fast axiomatic attribution for neural networks.  Advances in Neural Information Processing Systems, 34, 2021.
  * Hooker et al. [2019] Sara Hooker, Dumitru Erhan, Pieter-Jan Kindermans, and Been Kim.  A benchmark for interpretability methods in deep neural networks.  Advances in neural information processing systems, 32, 2019.
  * Kapishnikov et al. [2019] Andrei Kapishnikov, Tolga Bolukbasi, Fernanda Viégas, and Michael Terry.  Xrai: Better attributions through regions.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 4948–4957, 2019.
  * Kapishnikov et al. [2021] Andrei Kapishnikov, Subhashini Venugopalan, Besim Avci, Ben Wedin, Michael Terry, and Tolga Bolukbasi.  Guided integrated gradients: An adaptive path method for removing noise.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 5050–5058, 2021.
  * Li et al. [2021] Xin Li, Xiangrui Li, Deng Pan, and Dongxiao Zhu.  Improving adversarial robustness via probabilistically compact loss with logit constraints.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 8482–8490, 2021.
  * Lundberg and Lee [2017] Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  In Advances in neural information processing systems, pages 4765–4774, 2017.
  * Miglani et al. [2020] Vivek Miglani, Narine Kokhlikyan, Bilal Alsallakh, Miguel Martin, and Orion Reblitz-Richardson.  Investigating saturation effects in integrated gradients.  arXiv preprint arXiv:2010.12697, 2020.
  * Montavon et al. [2017] Grégoire Montavon, Sebastian Lapuschkin, Alexander Binder, Wojciech Samek, and Klaus-Robert Müller.  Explaining nonlinear classification decisions with deep taylor decomposition.  Pattern Recognition, 65:211–222, 2017.
  * Pan et al. [2020] Deng Pan, Xiangrui Li, Xin Li, and Dongxiao Zhu.  Explainable recommendation via interpretable feature mapping and evaluation of explainability.  In Christian Bessiere, editor, Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence, IJCAI 2020, pages 2690–2696. ijcai.org, 2020.
  * Pan et al. [2021] Deng Pan, Xin Li, and Dongxiao Zhu.  Explaining deep neural network models with adversarial gradient integration.  In Thirtieth International Joint Conference on Artificial Intelligence (IJCAI), 2021.
  * Petsiuk et al. [2018] Vitali Petsiuk, Abir Das, and Kate Saenko.  Rise: Randomized input sampling for explanation of black-box models.  arXiv preprint arXiv:1806.07421, 2018.
  * Qiang et al. [2022] Yao Qiang, Deng Pan, Chengyin Li, Xin Li, Rhongho Jang, and Dongxiao Zhu.  Attcat: Explaining transformers via attentive class activation tokens.  In Advances in Neural Information Processing Systems, 2022.
  * Ribeiro et al. [2016] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  " why should i trust you?" explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 1135–1144, 2016.
  * Shah et al. [2021] Harshay Shah, Prateek Jain, and Praneeth Netrapalli.  Do input gradients highlight discriminative features?  Advances in Neural Information Processing Systems, 34:2046–2059, 2021.
  * Shrikumar et al. [2017] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje.  Learning important features through propagating activation differences.  In International Conference on Machine Learning, pages 3145–3153. PMLR, 2017.
  * Simonyan and Zisserman [2014] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * Simonyan et al. [2013] Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.  Deep inside convolutional networks: Visualising image classification models and saliency maps.  arXiv preprint arXiv:1312.6034, 2013.
  * Smilkov et al. [2017] Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda Viégas, and Martin Wattenberg.  Smoothgrad: removing noise by adding noise.  arXiv preprint arXiv:1706.03825, 2017.
  * Srinivas and Fleuret [2019] Suraj Srinivas and François Fleuret.  Full-gradient representation for neural network visualization.  Advances in neural information processing systems, 32, 2019.
  * Sturmfels et al. [2020] Pascal Sturmfels, Scott Lundberg, and Su-In Lee.  Visualizing the impact of feature attribution baselines.  Distill, 2020.  https://distill.pub/2020/attribution-baselines.
  * Sundararajan et al. [2017] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  In International Conference on Machine Learning, pages 3319–3328. PMLR, 2017.
  * Szegedy et al. [2015] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jonathon Shlens, and Zbigniew Wojna.  Rethinking the inception architecture for computer vision. corr abs/1512.00567 (2015), 2015.