NewReferences

# The black hole low mass X-ray binary V404 Cygni is part of a wide
hierarchical triple, and formed without a kick

Kevin B. Burdge1,2 ∗ Kareem El-Badry3 Erin Kara1,2 Claude Canizares1,2 Deepto
Chakrabarty1,2 Anna Frebel1,2 Sarah C. Millholland1,2 Saul Rappaport1,2 Rob
Simcoe1,2 Andrew Vanderburg1,2

###### Abstract

Evidence suggests that when compact objects such as black holes and neutron
stars form, they may receive a “natal kick,” where the stellar remnant gains
momentum. Observational evidence for neutron star kicks is substantial[1, 2],
yet limited for black hole natal kicks, and some proposed black hole formation
scenarios result in very small kicks [3, 4]. Here, we report the discovery
that the canonical black hole low-mass X-ray binary V404 Cygni is part of a
wide hierarchical triple with a tertiary companion at least 3500 astronomical
units away from the inner binary. Given the orbital configuration, the black
hole likely received a sub-5 kilometer per second kick to have avoided
unbinding the tertiary. This discovery reveals that at least some black holes
form with nearly no natal kick. Furthermore, the tertiary in this system lends
credence to evolutionary models of low-mass X-ray binaries involving a
hierarchical triple structure[5]. Remarkably, the tertiary is evolved,
indicating that the system formed 3-5 billion years ago, and that the black
hole has removed at least half a solar mass of matter from its evolved
secondary companion. During the event in which the black hole formed, it is
likely that at least half of the mass of the black hole progenitor collapsed
into the black hole; it may even have undergone a complete implosion, enabling
the tertiary to remain loosely bound.

Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA

Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, Cambridge, MA 02139, USA

Division of Physics, Mathematics and Astronomy, California Institute of
Technology, Pasadena, CA, USA

Black holes (BHs) are the stellar remnants of the most massive stars.
Gravitational wave astronomy has raised key questions regarding the formation
of these objects, including the role of natal kicks[6] and dynamical
interactions[5, 7, 8]. There are several proposed formation pathways for these
stellar remnants, including some channels in which there is an implosion with
no associated natal kicks or supernovae [3, 9, 10]. The best constraints on BH
kick physics come from X-ray binaries’ Galactic orbits[11], astrometric
microlensing[12], and orbital dynamics[13]. These constraints have largely
ruled out the need for kicks of $>100\,\rm km\,s^{-1}$ except for in a handful
of systems with exceptional orbits[14, 15].

We report the discovery that the black hole X-ray binary V404 Cygni, the first
low-mass X-ray binary (LMXB) widely accepted to host a black hole (with a BH
mass of $9^{+0.2}_{-0.6}\rm\,M_{\odot}$[16]), is part of a wide hierarchical
triple. Our serendipitous discovery stemmed from examining an optical image of
V404 Cygni on the Aladin Lite tool, illustrated in panel a of Figure 1, and
taking note that there is a nearby star just $1.43$ arcseconds from V404 Cygni
with _Gaia_ proper motions matching those of V404 Cygni, as seen in panel b of
Figure 1. We investigated the proper motions of nearby sources, and in panel a
of Figure 2, show that such a chance agreement in proper motions is unlikely.
In panel c of Figure 2, we quantify the probability of such an alignment
occurring by chance, which we find is about $10^{-7}$ (Methods). While
searching the literature on V404 Cygni, we found several works that commented
on the nearby star, with most simply assuming that it was an interloper[17,
18]. One notable exception was Maitra et al. 2017 [19], which speculated “in
passing whether the blended star is truly unrelated to the V404 system” due to
the similar estimated distance and extinction of the blended star. V404 Cygni
was noted as having a peculiar velocity in previous work[20], and this was
attributed to a kick. As part of our astrometric analysis, we investigated
nearby stars and found that the velocity of V404 Cygni is typical of stars in
the vicinity (Methods).

Given the 1.43 arcseconds separation and estimated distance of $2.39\pm
0.15\,\rm kpc$[20], we find that the wide tertiary companion is at least 3500
astronomical units (AU) away from the inner binary, which corresponds to a
Keplerian velocity of just a few kilometers per second. As illustrated in
panels a and c of Figure 1, this separation is approximately 90 times larger
than the distance of Pluto to the Sun, and 25000 times larger than that of the
inner V404 Cygni binary, which has a separation of just 0.14 AU, less than
half the distance between the Sun and Mercury.

To further confirm this association, we analyzed archival VLT X-shooter
spectroscopic observations (which targeted V404 Cygni, but contained a
spatially resolved trace of the tertiary) and obtained additional follow-up
GMOS spectroscopy of the source (Methods). We conducted a radial velocity
analysis and found excellent agreement with the reported systemic velocity of
V404 Cygni. As seen in the histogram of radial velocities of nearby stars
shown in panel c of Figure 2, this is unlikely to have occurred by chance,
further solidifying the association of the two components.

We fit the spectroscopic observations with model atmospheres, focusing on the
region around the Hydrogen Alpha and Beta absorption lines in the GMOS spectra
(Methods), and the region around the Calcium triplet lines in the X-shooter
observations, shown in panel d of Figure 3. In addition to fitting the
spectroscopic observations with model atmospheres, we also fit the broadband
spectral energy distribution (SED), carefully extracting photometry by
modeling the point-spread function in epochal Pan-STARRS images of the source
due to the blending of the tertiary and V404 Cygni (Methods). Additionally, we
use archival Hubble Space Telescope observations to obtain a measurement in
the near ultraviolet at 330 nm, and archival observations with the Keck
observatory’s NIRC2 instrument to measure the near-infrared flux at $\sim 2$
microns. By jointly fitting the spectra and SED, we obtained constraints on
the temperature, metallicity, and radius, with values reported in panel c of
Figure 3.

One significant result in our modeling of the SED and spectra is that we find
the tertiary in V404 Cygni has started to evolve off of the main sequence and
is about twice its initial radius. By fitting the tertiary with MIST
isochrones shown in panel a of Figure 3, we find that this constrains the
system’s age to about 3-5 gigayears, and the mass of the tertiary to around
1.2 solar masses (our full parameter estimates can be found in panel c of
Figure 3).

V404 Cygni’s secondary, like the donors seen in some other BH LMXBs, exhibits
enhanced lithium abundance, and this has been attributed to formation as a
result of accretion processes onto the black hole, or in the supernova that
formed it[21]. Given our constraint on the age of the system, we can rule out
that this lithium abundance is a result of recent formation; however, we
inspected our X-shooter and GMOS spectra and concluded that at this time there
is insufficient signal-to-noise and resolution to determine whether the
tertiary has enhanced lithium.

To learn about the physical constraints imposed on the formation of the black
hole, we simulated the dynamics of the triple with a range of configurations,
accounting for the BH kick, mass lost during the BH formation, and the initial
orbital periods of the secondary and tertiary, to investigate which BH
formation scenarios could retain the loosely bound tertiary.

As seen in panel a of Figure 4, the only way the inner binary could have
experienced a large kick and retained the tertiary, is if the tertiary started
at a short orbital period, and was kicked into a highly eccentric orbit
reaching $>3500$ AU at apastron. This scenario is unlikely, as one would need
to fine-tune the kick to the inner system to be large, but just barely below
the escape velocity of the system (Methods). We find that scenarios in which
the tertiary started in a wide orbit, and the inner binary received a small
kick of just a few kilometers per second are thus strongly favored.

When considering the inner binary, we simulated two scenarios. In one case, we
allowed the secondary to start in a wide orbit between 100 AU and 300 AU. This
scenario is viable for reproducing the current system, as we find that von
Zeipel-Lidov-Kozai cycles could readily cause the secondary to migrate into
its current 6.4-day orbit (Methods). As illustrated by the red histograms in
panel b of Figure 4, we find that in this scenario, the BH kick was likely
smaller than 3 kilometers per second and that the mass lost in the BH
formation could have been up to about 10 solar masses, or about half the mass
of the inner binary.

Alternatively, we consider a scenario where the secondary starts with an
initial orbital period between 1 and 6 days (its current orbit is 6.4 days).
These simulations are illustrated as the blue histograms in panel b of Figure
4. In this case, the inner binary’s barycenter receives a significant Blaauw
kick[22] as a result of mass loss, even if the BH itself does not receive a
kick. This results from matter being ejected in the BH progenitor’s rest
frame, causing it to be jettisoned out of the binary at the orbital velocity
of the BH in the barycentric frame. At orbital periods of days, the BH orbits
with a velocity of a few 10s of kilometers per second, and thus, as seen in
panel b of Figure 4, this scenario strongly constrains any possible mass loss
during the BH formation, because even if the BH does not receive a kick, the
inner binary does as a consequence of the mass loss, ejecting the tertiary.
One curiosity of this scenario is that in fine-tuned cases, one can achieve
slightly larger BH kicks while keeping the tertiary bound, because the
velocity imparted on the barycenter of the inner binary by mass loss can
absorb the effect of the kick oriented in the opposite direction, resulting in
a relatively low net kick to the inner binary, allowing the retention of the
tertiary. Overall, this scenario still favors small kick velocities of less
than 5 kilometers per second and does not allow for more than about a solar
mass to be ejected. Thus, if the secondary started in a tight orbit of a few
days, the most likely BH formation scenario is one in which there was a near-
complete implosion of the progenitor star, with negligible mass loss.

In either scenario, a BH formation event in which there is a complete
implosion and no kick always results in the survival of the system. This would
challenge current models for such a scenario, as this black hole has a mass of
just $\sim 9\rm\,M_{\odot}$, which is a smaller BH mass than predicted for
such an event[19]. We cannot rule out mass loss in the system, but we find it
is improbable that more than half the BH progenitor’s mass was suddenly
ejected during its formation.

We consider it unlikely that the system formed dynamically in a dense
environment such as a globular cluster or the Galactic center and was ejected,
as the escape velocities of these environments exceed the orbital velocity of
the tertiary by an order of magnitude (the young age of the tertiary also
disfavors an origin in a globular cluster). We also find it unlikely that the
system captured the tertiary in the field as a result of the low cross-section
of such an interaction occurring.

The presence of a tertiary companion in one of the most well-known LMXBs
supports theoretical work which has suggested that hierarchical triples may be
key to forming BH LMXBs. Forming BH LMXBs via purely binary evolution has been
theoretically challenging because of the large mass ratios involved, resulting
in a common envelope event that proceeds to a merger rather than a successful
ejection of the envelope. Thus, theoretical modeling of the formation of such
systems has explored the possibility that wide tertiary companions helped the
donor migrate into a tight orbit after the formation of the BH[5]. This is
achieved through a gravitational interaction of the tertiary with the inner
orbit, known as a von Zeipel-Lidov-Kozai cycle, in which the inner orbit
cycles between an inclined and eccentric orbit. We note that models such as
those presented in Naoz et al. 2016[5] predict tertiary companions in orbits
at $\sim 10^{4}\rm AU$, which is consistent with what we observe in V404
Cygni. One plausible evolutionary scenario for V404 Cygni is that the inner
binary began its life with a separation of $\sim 10^{2}\rm AU$, and over time
Kozai-Lidov cycles drove this to a shorter orbit as a result of tidal
dissipation and magnetic breaking draining orbital angular momentum from the
inner binary during the highly eccentric phase. Finally, the inner binary
hardened to an orbit of $<6.5$ days, too long for Roche-lobe overflow while
the secondary was on the main-sequence, but as the secondary evolved, it
overflowed its Roche-lobe. This is essentially the giant sub-channel described
in Naoz et al. 2016 [5].

We note that searches for wide binary companions are, in general, highly
incomplete. To be detected by Gaia, a companion must (a) be brighter than
$G\sim 20.7$, and (b) must be separated from the inner binary by at least
$\sim 1$ arcsec (depending somewhat on flux ratio). For V404 Cygni, this
corresponds to detection limits of $M\gtrsim 0.9\,M_{\odot}$ for main-sequence
stars, and separations $s\gtrsim 2500$ AU. The separation distribution of
solar-type tertiaries peaks at 10-100 AU, with companions at 100-2500 AU
outnumbering those at $>2500$ AU by a factor of two[23, 24]. While the
separation distribution of tertiaries to massive stars is quite uncertain,
this suggests that companions too close to be detected by Gaia may be common
and could have evaded detection thus far. Indeed, several BH LMXBs are thought
to have unresolved companions, which have thus far been interpreted as chance
alignments[25, 26, 27]. Proper motions for these companions have not yet been
measured.

V404 Cygni is nearer and brighter in the optical than most other BH LMXBs. If
the system were $\sim$50% more distant, the companion would be blended with
the inner binary, and Gaia would not have been able to measure its proper
motion. If the companion were $\gtrsim 10\%$ more massive, it would already be
a faint white dwarf below Gaia’s detection limit. If it were $\gtrsim 40\%$
less massive, it would be a main sequence star below the Gaia detection limit;
if it were $\gtrsim 20\%$ less massive it would be too faint for Gaia to have
measured a precise proper motion an establish the association with high
confidence. These considerations all suggest that harder-to-detect tertiaries
may well be hiding around other known BHs. It is quite possible that most BH
LMXBs formed through triple evolution, and deeper searches around other BHs
hold promise to detect them.

Evidence suggests that the spin axis of the BH in V404 Cygni is misaligned
from the orbital plane due to the rapidly changing orientation of the jets in
the system, and this has been attributed to a natal kick[28]. However, the
tertiary’s presence largely rules out a natal kick. An alternative possibility
is that the von Zeipel-Lidov-Kozai cycles induced this misalignment, as the
inner orbit would have evolved through a range of inclinations and may have
hardened at an inclination misaligned with the original orbit[29, 30]. If the
current spin of the black hole in the system was primarily inherited from the
progenitor star, this could naturally lead to a misalignment of this spin axis
with the current orbital plane. However, the evolved 1.2 solar mass tertiary
implies the black hole has removed at least 0.5 solar masses from the 0.7
solar mass secondary—which was originally more massive than the tertiary, as
it evolved first. If the black hole conservatively accreted this much mass, it
could account for the large spin, but one would not expect a misalignment from
the orbital plane.

The tertiary companion of V404 Cygni has provided favorable evidence for the
formation of at least some low-mass x-ray binaries in hierarchical triples.
Moreover, it has provided one of the strongest empirical constraints on natal
kicks in the formation of a black hole by indicating that the BH in V404 Cygni
likely formed with a kick of less than five kilometers per second,
demonstrating that at least some stellar mass black holes form without
substantial natal kicks. We conclude by noting that our simulations strongly
suggest that V404 Cygni’s secondary either started in a wider orbit, and
migrated in as a result of von Zeipel-Lidov-Kozai interactions, or if the
secondary originated in a tight orbit, the $9\,M_{\odot}$ BH formed without
ejecting more than a solar mass of matter–a near complete implosion.

Figure 1: a): A Pan-STARRS image of V404 Cygni and its companion, displayed
using the Aladin interface, with V404 Cygni and the tertiary labeled in red
and the separation of 1.43 arcseconds indicated in blue. We discovered the
tertiary while viewing the source on Aladin and inspecting the Gaia astrometry
on the interface, finding remarkable agreement in the measured proper motions
of these two sources. b): A plot of the positions and proper motion vectors of
all stars in the field, with V404 Cygni and its tertiary indicated in red. c):
A zoom-in on the inner binary of V404 Cygni, illustrating the $2.5\times
10^{4}$ ratio of the semi-major axes of the inner and outer orbits in the
triple.  Figure 2: a): A plot illustrating the proper motions of all Gaia
stars brighter than 18th magnitude within 5 arcminutes of V404 Cygni, and the
remarkable agreement between V404 Cygni and its tertiary (shown in red). b): A
histogram of Gaia-measured radial velocities of stars in the vicinity of V404
Cygni, with the systemic velocity of V404 Cygni shown as a black vertical
line, and the measured radial velocity of the tertiary shown as the grey
vertical line. c): A plot illustrating the chance alignment probability as a
function of separation for stars brighter than 18th magnitude (shown in red),
and 20th magnitude (shown in blue). This probability accounts for both the
agreement in proper motions, as well as radial velocity.  Figure 3: a): MIST
stellar isochrones, with the color bar indicating the age of stars evolving
according to these isochrones. The black star with error bars indicates the
position of the tertiary, which overlaps with tracks in the range of 1.2-1.3
solar masses at ages of around 3-5 billion years. b): The spectral energy
distribution of the tertiary in V404 Cygni (black diamonds), with the bluest
point coming from Hubble Space Telescope ACS observations, the reddest point
from Keck NIRC2 observations, and the remaining photometric measurements from
Pan-STARRS images. The red curve shows our best-fit spectrum to this data, and
the blue crosses indicate the filter-averaged flux values from this spectrum.
c): Our derived parameters as a result of the models illustrated in this
figure. d): The X-shooter spectrum of the Calcium triplet absorption features
and other nearby absorption lines in the tertiary (black), and our best fit
spectral model to it (red).  Figure 4: a): An illustration of the range of
possible kick velocities to the inner binary that could result in a bound
tertiary at a separation greater than 3500 AU. In general, large kicks are
only allowed if the tertiary originated at a very short orbital period, a
scenario that is unlikely to retain the tertiary (Methods). b): Histograms
representing our simulations of black hole kicks and mass loss in the system.
The red histograms represent a scenario in which the secondary originates at a
wide range of orbital periods, whereas the blue histograms represent a
scenario in which the secondary started at very short orbital periods.

We dedicate this work to the memory of our dear friend, Tom Marsh. K.B.B. is a
Pappalardo Postdoctoral Fellow in Physics at MIT and thanks the Pappalardo
fellowship program for supporting his research. This research was supported by
NSF grant AST-2307232.

The authors declare that they have no competing financial interests.

Correspondence and requests for materials should be addressed to K.B.B.
(email: kburdge@mit.edu).

##  0.1 Astrometric association

We investigated the robustness of the astrometric association of the tertiary
with the inner binary by analyzing a sample of all _Gaia_ sources within $30$
arcminutes of V404 Cygni. In _Gaia_ DR3, there are 80739 sources in this
region of the sky. This corresponds to a source density of $0.0079$ sources
per square arcsecond. As an initial check, we parsed these 80739 sources for
objects with proper motion values that fall within 1 sigma of V404 Cygni’s
values and found only 21 such sources. When we account for the error bars of
not just V404 Cygni’s measured proper motions, but also of all other sources,
we find that 1256 sources have 1 sigma confidence intervals that overlap the 1
sigma confidence interval of V404 Cygni. However, the mean magnitude in this
sample is _Gaia_ G of 20.5, and thus it is dominated by sources with large
uncertainties in proper motions, meaning that V404 Cygni’s well-measured
proper motion is consistent with over 1000 poorly measured proper motions in
this region.

To construct the curve shown in panel c of Figure 2, which we define as the
probability of finding a source within a given separation radius that has
proper motions in RA and Dec consistent within 1 sigma, as well as a radial
velocity consistent to within 1 sigma, we take the list of 80739 Gaia sources
within 30 arcminutes, and further downselect to 39577 sources with a Gaia G
magnitude greater than 20 (to construct the blue curve), and 10798 sources
with a Gaia G magnitude greater than 18 (to construct the red curve). We
follow the procedure outlined in El-Badry et al. 2018[31] to determine whether
two sources have consistent proper motions, accounting for possible orbital
motion influencing the astrometric solution. We find that for stars with
_Gaia_ $G>18$ (V404 Cygni’s tertiary is 17.9), there are only 3 sources within
30 arcminutes that have proper motions consistent with V404 Cygni (this number
ranges from 2-4 sources, as the number of sources consistent with the proper
motions of V404 Cygni depends on the assumed separation, as this is used to
compute possible variance in proper motions due to the orbital motion). In any
case, this leads to small chance alignment probabilities. Extended Data Table
1 lists the _Gaia_ astrometric solution of V404 Cygni and its tertiary.

We also consider the radial velocity of the sources in computing the chance
alignments. To do this, we selected all Gaia sources within 30 arcminutes of
V404 Cygni with a Gaia radial velocity error less than 5 kilometers per second
(this is comparable to our measured uncertainty on the RV of the tertiary). We
used these sources to construct the histogram seen in panel b of Figure 2.
From this distribution, we computed that there is approximately an 8.1 percent
probability of a source with a well-measured RV being consistent with that of
V404 Cygni to within 1 sigma, and we incorporated this information in our
chance alignment probabilities illustrated in panel c of Figure 2.

To translate these numbers into formal chance alignment probabilities, we ask,
what is the probability that such a source would fall within a circle of some
radius (represented by the x-axis in panel c of Figure 2) around V404 Cygni.
We use the vertical black line to denote the location of the actual tertiary,
which is 1.43 arcseconds away from V404 Cygni. A 1.43 arcsecond region
represents just 1/1584429 of the 30 arcminute region we queried for sources,
and thus the probability of any individual source falling by chance into such
a region is $6.3\times 10^{-7}$. Thus, if we consider the 3 sources with
_Gaia_ $G>18$ that have proper motions consistent with V404 Cygni, accounting
for possible orbital motion at a 1.43 arcsecond separation, the probability of
one of those sources falling within 1.43 arcseconds of V404 Cygni is
$1.26\times 10^{-6})$, and when we account for the probability of the radial
velocity agreeing to within 1 sigma, this becomes $\sim 10^{-7})$.

While analyzing the astrometry of V404 Cygni and its tertiary, we used the
astrometry of nearby stars to investigate the peculiar velocity reported in
Miller-Jones et al. 2009 [20]. Using stars within 1 degree of V404 Cygni, with
similar distances (parallax values between 1/3 and 1/2), and radial velocity
errors of less than $5\,\rm km\,s^{-1}$, we constructed the Toomre diagram
shown in Extended Data Figure 1, and found that the _Gaia_ astrometry and
radial velocities of nearby stars are largely consistent with that of V404
Cygni.

Table 1: Gaia astrometry of V404 Cygni and the Tertiary

Object | RA (epoch 2016) | Dec (epoch 2016) | $\omega$ | $\mu_{\rm RA}$ | $\mu_{\rm Dec}$
---|---|---|---|---|---
V404 Cygni | $306.0159085525$ | $33.8671768648$ |  $0.3024\pm 0.0783$11footnotemark: 1 | $-5.1775\pm 0.0785$ | $-7.7776\pm 0.0922$
Tertiary | $306.0159299085$ | $33.8675732739$ | $0.1423\pm 0.1161$ | $-5.1500\pm 0.1564$ | $-7.7647\pm 0.1316$

Table 1: 1 We note that the Gaia parallax measurement is much less precise
than the radio parallax of $\omega=0.418\pm 0.024$ reported in [20]. The
proper motions are consistent with what was measured in the radio. Extended
Data Figure 1: A Toomre diagram illustrating V404 Cygni’s Galactic orbital
velocity relative to stars within a degree of V404 Cygni. We find that V404
Cygni’s velocity is largely consistent with nearby stars.

##  0.2 Constraining kicks and mass loss

To constrain the possible natal kicks the black hole could have received while
still retaining the tertiary, we treated the system as a triple with component
masses of $M_{BH}=8.5\,M_{\odot}$ as the post supernova black hole mass,
$M_{2}=1.5\,M_{\odot}$ as the initial mass of the inner donor star, and
$M_{3}=1.2\,M_{\odot}$ as the initial mass of the tertiary.

We consider four free parameters in our simulations: the initial orbital
period of the tertiary, the initial orbital period of the secondary, the
amount of mass ejected in the supernova, and the kick that the black hole
received. We marginalized over the kick angle and the orientation of the inner
binary with respect to the outer tertiary’s orbit, and in all cases we assumed
initially circular orbits. In all cases, we assume initial orbits are
circular. To determine whether the post-kick orbits remain bound, we follow
the procedure outlined in Brandt & Podsiadlowski (1995)[32], which we outline
here. We define a non-dimensional characteristic mass:

$\tilde{m}=\frac{M_{1}+M_{2}}{M^{\prime}_{1}+M^{\prime}_{2}},$ (1)

where $M_{1}$ and $M_{2}$ are the initial masses in the system, and
$M^{\prime}_{1}$ is the post-supernova mass of component that underwent mass
loss. We also compute a dimensionless ratio of the magnitude of the kick
velocity to orbital velocity:

$\tilde{v}=\frac{v_{\textrm{kick}}}{v_{\textrm{orb}}},$ (2)

where $v_{\textrm{orb}}$ is the relative orbital velocity of the two
components computed using Kepler’s laws (in the case of the inner binary, the
two components are the BH progenitor and the secondary, and in the case of the
outer orbit, the components are the tertiary and the center of mass of the
inner binary).

Following Brandt & Podsiadlowski (1995)[32], this yields a post-supernova
energy in the new center of mass frame given by:

$E^{\prime}=-\frac{GM^{\prime}_{1}M_{2}}{2a}\left[2-\tilde{m}(1+2\tilde{v}\cos\phi\cos\theta+\tilde{v}^{2})\right],$
(3)

where $\phi$ and $\theta$ represent the polar and azimuthal kick angles,
respectively. As observed in Brandt & Podsiadlowski (1995)[32], for the final
system to remain bound, $E^{\prime}$ must be negative, and the final semi-
major axis is given by

$a^{\prime}=-\frac{GM^{\prime}_{1}M_{2}}{2E^{\prime}}.$ (4)

In each simulation, we first use the BH kick and level of mass loss and
consider the impact this has on the inner orbit’s Barycentric velocity,
$v_{\text{sys}}$. We then use this $v_{\text{sys}}$ as the input to a
calculation to check whether the outer tertiary will remain bound to the inner
binary. We compute $v_{\text{sys}}$ using the expression given in Brandt &
Podsiadlowski (1995)[32]:

$v_{\text{sys}}=\frac{v_{\text{orb}}}{M^{\prime}_{1}+M_{2}}\left[\left(\frac{\mu\Delta
M_{1}}{M_{1}}\right)^{2}-2\frac{\mu\Delta
M_{1}M^{\prime}_{1}}{M_{1}}\tilde{v}\cos\phi\cos\theta+(M^{\prime}_{1}\tilde{v})^{2}\right]^{1/2},$
(5)

where $\Delta M_{1}$ is the mass lost in the supernova, and $\mu$ is the
reduced mass of the pre-supernova binary, defined as
$\mu=\frac{M_{1}M_{2}}{M1+M2}$.

In constructing the histograms shown in panel b of Figure 4, we consider two
cases: one in which the inner companion started out in an orbit between 1 and
6 days (we use a uniform distribution of orbital periods between these two
values), and another case in which it started out at an orbital separation
between 100 and 300 AU (we also always enforce that the secondary must start
out at a shorter orbital period than the tertiary). We use 100 AU as a lower
bound to ensure that the von Zeipel-Lidov-Kozai timescales were short enough
that the secondary would have migrated into its current orbit over the few Gyr
age of V404 Cygni. We enforce the 300 AU upper bound because if the orbit were
significantly wider, it would be dynamically unstable. In both scenarios, we
drew from uniform distributions for the BH kick and mass loss during BH
formation, with the BH kick being drawn from a uniform distribution between
$0\,km\,s^{-1}$ and $300\,km\,s^{-1}$, and the mass loss being drawn from a
uniform distribution ranging from $0\,M_{\odot}$ and $20\,M_{\odot}$. We
consider “surviving” solutions to be ones where the secondary remains bound in
an orbit with a semi-major axis less than the final semi-major axis of the
tertiary and the tertiary remains bound, with a final semi-major axis greater
than $3500\,\rm AU$.

To construct panel a of Figure 4, which represents the range of possible kicks
the inner binary could have experienced while retaining the tertiary as a
function of the initial orbital period of the tertiary, we iterate over
initial orbital periods of the tertiary and simply apply a range of systemic
velocities to the inner binary, without assuming any mass loss, and
investigate which solutions retain the tertiary in a bound orbit greater than
$3500\,\rm AU$.

Extended Data Figure 2: a): The probability of retaining a bound tertiary as a
function of the kick velocity in the inner binary. In general, this
probability diminishes at large kicks because these are only allowed when the
tertiary originates in a tight orbit, and must be fine-tuned to be large
enough to send it to a wide separation without unbinding it. b): The
probability of retaining a bound tertiary at greater than 3500 AU as a
function of the initial orbital period of the tertiary. These probabilities
are diminished for tertiaries that originate in tight orbits due to the large
degree of fine-tuning required to send the tertiary to such a wide orbit
without unbinding it.

##  0.3 Spectroscopic Observations

We searched the ESO archive for observations of V404 Cygni in which the
companion may have serendipitously fallen inside the slit. This yielded a
VLT/X-shooter spectrum obtained on 16 July 2015 (Program 295.D-5027;
PI:Rahoui), between the main outburst and mini-outburst [33], when the black
hole was in quiescence. To our knowledge, these data have not been published
elsewhere. As seen in panel A of Extended Data Figure 3, two separate traces
are clearly visible in the raw spectra, separated by about 1.4 arcsec, with
one 5-10 times brighter than the other. Since V404 Cygni has no other
comparably bright neighbors, there is little doubt that the fainter of the two
traces corresponds to the companion. The seeing was about 1.3 arcsec (FWHM),
meaning that the center of the companion’s trace is separated from V404 Cygni
by 2.5 times the “$\sigma$” of the seeing disk. Given that the companion is
5-10 times fainter than V404 Cygni in the VIS data, with its relative flux
contribution increasing toward redder wavelengths, the companion was likely
only partially in the slit. The data were taken with the 1.2 arcsec slit,
yielding a typical resolution of 6500 in the VIS band.

The mid-exposure time is HJD 2457219.692, when the ephemeris of Casares et al.
2019[33] predicts an RV of $-49.9\,\rm km\,s^{-1}$ for the donor. Twelve
separate 148s exposures were obtained sequentially in nod-and-shuffle mode,
for a total exposure time of 1776s in the VIS band. After inspecting the 12
exposures, we rejected one in which the companion was unusually faint,
presumably because it fell farther outside the slit. We reduced and combined
the other 11, for an effective exposure time of 1628s in the VIS band. We
reduced the data using the ESO Reflex pipeline[34] with standard calibrations.
This performs bias-subtraction, flat fielding, wavelength calibration using
afternoon ThAr arcs, and order merging. We set the extraction window for the
two sources manually using the localize-slit-position and localize-slit-height
parameters. To minimize contamination from V404 Cygni, we extracted only the
$\approx 50\%$ of the companion’s trace on the far side of V404 Cygni. The
fact that the companion’s extracted spectrum shows no emission in H$\alpha$,
where V404 Cygni shows a strong double-peaked emission line (see panels b and
c of Extended Data Figure 3), suggests there is little contamination.

Extended Data Figure 3: a): A region in the archival X-shooter spectrum of
V404 Cygni, with the traces associated with the inner binary of V404 Cygni,
and the tertiary both labeled. b): The extracted X-shooter spectra around the
calcium triplet region for V404 Cygni (shown in red), and the tertiary (shown
in blue), illustrating that we were able to extract a spectrum of the tertiary
with minimal contamination by V404 Cygni. c): The spectrum of V404 Cygni
(shown in red) and the tertiary (blue) around the H Alpha Hydrogen absorption
line.

We used the “A” and “B” telluric bands to perform a flexure correction in the
wavelength solution, obtaining a $-12\rm\,km\,s^{-1}$ correction with “A”
band, and a $-11.0\rm\,km\,s^{-1}$ with “B” band, and adopted a correction of
$-11.5\rm\,km\,s^{-1}$, with a systematic uncertainty of $-1\rm\,km\,s^{-1}$.
After applying these corrections, we converted the X-shooter spectrum to air
wavelengths, and applied the Barycentric correction of $8.6\rm\,km\,s^{-1}$ in
the file headers.

We also obtained an additional spectroscopic observation using Gemini Multi-
Object Spectrographs (GMOS) on the 8.1-m Gemini North telescope on Mauna Kea,
to determine whether there was any RV variability within the tertiary (program
GN-2023B-DD-105). We used the 0.5-arcsecond slit and the R831_G5302 grating,
which provides a resolution of approximately 4400. Our spectrum covered a
wavelength range from $4576\,\text{\AA}$ to $6925\,\text{\AA}$. While being
slightly lower resolution than the X-shooter data, the spectrum provided a
significantly higher signal-to-noise ratio (SNR) at bluer wavelengths.

Extended Data Figure 4: a): The Gemini GMOS-N spectrum of the H Beta
absorption feature in the tertiary (black), and our best fit spectral model to
it (red). b): The Gemini GMOS-N spectrum of the H Alpha absorption feature in
the tertiary (black), and our best fit spectral model to it (red).

We aligned the position angle of the slit perpendicular to the angle between
V404 Cygni and the tertiary, effectively masking V404 Cygni. The observation
was obtained in good seeing (FWHM $\approx 0.7$ arcsec), meaning that V404
Cygni was nearly 4 times the “$\sigma$” of the seeing disk from the nearest
edge of the slit, and contamination is expected to be negligible. We used a
900 second exposure, yielding SNR of $\approx 30$ per pixel at 6600 Å,
$\approx 10$ per pixel at 5600 Å, and $\approx 5$ per pixel at 5000 Å. We
obtained a CuAr arc on-sky immediately after the science exposure.

We reduced the GMOS data using PypeIt, which required manual construction of a
new template for the R831 grating. PypeIt performs bias and flat field
correction, cosmic ray removal, wavelength calibration, sky subtraction,
extraction of 1d spectra, and heliocentric RV corrections. Like the X-shooter
data, we checked the wavelength solution using the telluric bands, and in this
case found that no correction was needed.

##  0.4 Radial Velocity Measurements

We used the X-shooter spectra to measure the radial velocities (RVs) of both
components, as well as perform atmospheric fits in conjunction with modeling
the spectral energy distribution (SED). To measure RVs, we used the ultranest
kernel density estimator to fit an atmospheric model constructed using the
PHOENIX stellar atmosphere library[35]. We selected a region between
$8480\rm\,\text{\AA}$ and $8700\rm\,\text{\AA}$ for the fit, as this
wavelength range contains the calcium triplet lines, which provide a strong RV
signal in the x-shooter spectrum. We allowed the RV to be a free parameter,
and fixed the other parameters in this atmospheric model to an effective
temperature of $T_{eff}=6100\,\rm K$, $logg=3.95\,\rm cgs$, and $[Fe/H]=0.0$
based on our joint atmospheric+SED fit. Our fit yielded an RV measurement of
$-4.5\pm 2.5\pm 1\rm\,km\,s^{-1}$ for the tertiary, in excellent agreement
with the $-2.0\pm 0.4\rm\,km\,s^{-1}$ systemic velocity of the inner binary
reported in Casares et al. 2019[33].

We also conducted a radial velocity analysis of the GMOS spectra, fitting the
region between $6400\rm\,\text{\AA}$ and $6700\rm\,\text{\AA}$ to capture the
H alpha absorption line and nearby features, and the region between
$4750\rm\,\text{\AA}$ and $5000\rm\,\text{\AA}$ for the H beta absorption
feature. We find that this analysis yields a radial velocity of $-4.1\pm
1.9\rm\,km\,s^{-1}$, consistent with that of the X-shooter data, and of V404
Cygni’s systemic velocity. We note that the orbital motion of the tertiary is
on the order of a few kilometers per second, and may also introduce some
difference between the systemic velocities of the two components at this
level.

##  0.5 Spectral energy distribution

To model the spectral energy distribution of the tertiary in V404 Cygni, we
used a combination of observations obtained with the Panoramic Survey
Telescope and Rapid Response System (Pan-STARRS)[36], the Advanced Camera for
Surveys (ACS)[37] aboard the Hubble Space Telescope (HST), and near-infrared
camera (NIRC2) on the W.M. Keck Observatory[38].

Because of the overlapping point spread function of V404 Cygni’s inner binary
and the tertiary in the Pan-STARRS images, we performed our own point spread
function photometry of the source by constructing a routine using the
photutils package to iteratively run over the epochal Pan-STARRs images. We
used a small cutout region, illustrated in panel b) of Extended Data Figure 5
for our analysis, constructing a PSF model using the nearby bright stars near
the top and bottom of this cutout. After constructing the PSF model, we use
photutils to identify bright sources in each image and fit for the position of
the PSF, as seen in panel c) of Extended Data Figure 5, and we found that this
yielded good quality photometry, leaving only small residuals, as seen in
panel d) of Extended Data Figure 5. We note that initially when we performed
this analysis on the stacked Pan-STARRS images, we were unable to eliminate
substantial residuals, and discovered after inspecting the epochal images that
there was a nearby focal plane artifact, in the shape of a heart, which may
have been contaminating some of the stacked images. One example of this
artifact is illustrated in panel a) of Extended Data Figure 5, and we found
that its position regularly changed at different pointings, sometimes
contaminating V404 Cygni or our comparison stars, leading to a variable PSF
model over the stacked images. To quantify the uncertainties in the apparent
magnitude of the V404 Cygni tertiary, we computed the standard deviation of
the estimated apparent magnitude across all epochal Pan-STARRS images for
which we were able to perform photometry (in each filter). In addition to
using the photometric scatter across epochal images to estimate the
measurement error, we compared the inferred apparent magnitude using the star
just to the north of the triple, and the one just to the south of the triple,
to quantify the systematic error in each filter. We note that we used the Pan-
STARRS1 reference catalog of apparent magnitudes for our two comparison stars
as a basis for calibrating the photometry.

We used an archival HST observation (Proposal ID: 9686, PI: Hynes) to extract
flux in the ACS/HRC F330W filter. These observations consist of two 600s
exposures. As seen in panel a of Extended Data Figure 6, the tertiary is
visible and well resolved from V404 Cygni in these images. We used the
ACSTools python package to compute the flux of the tertiary, using a 0.2
arcsecond aperture and applying the appropriate corrections to the flux to
account for the portion of the PSF not encompassed by this aperture.

We performed infrared photometry using a single archival NIRC2 image of V404
Cygni, shown in panel b of Extended Data Figure 6. We used a 0.5 arcsecond
aperture to extract flux for the tertiary, as well as a comparison star just
to the north of the tertiary. We used the 2MASS reference magnitude of this
reference star to calibrate our flux, though we elected to apply a systematic
error of ten percent to our photometry due to the difference in filters in the
2MASS system (2MASS $K_{s}$), and the $K_{p}$ filter used in the NIRC2
observation. We elected not to perform more in-depth photometry on the full
archive of hundreds of images because we found that this photometric
measurement had little overall impact on our results.

Extended Data Figure 5: a): An example of a heart-shaped artifact that
contaminates many of the Pan-STARRS images of V404 Cygni. We believe this
artifact degraded the quality of the stacked PS images by occasionally
occluding either V404 Cygni or comparison stars, resulting in a non-uniform
PSF across the stacked images. This is why we ultimately decided to use
epochal images as the basis of our photometry. b): An epochal Pan-STARRS
r-band image of V404 Cygni. We used epochal images such as this to extract
fluxes for our spectral energy distribution analysis. c): An image
illustrating the point spread function (PSF) model we applied to the image to
extract photometry for V404 Cygni’s tertiary. d): A cutout illustrating the
residuals after we extract PSF photometry, with the PSF model constructed
using the bottom star in the image (hence the zero residuals around this
star).

##  0.6 Joint SED, Spectroscopic, and Isochrone analysis

We performed an analysis in which we fit the SED with model atmospheres,
synthesizing fluxes using the pyphot tool. In addition to the SED measurement,
the only other measurement we included in this analysis was the estimated
radio parallax of $0.418\pm 0.024$ to constrain the distance (and
appropriately account for the uncertainty in distance in estimating the
radius), and the reddening value of $E(g-r)=1.27$ reported in Green et al.
2019[39]. As with the analysis of the spectra, we used the ultranest kernel
density estimate to perform this fit. We fit for four free parameters, the
distance, effective temperature, metallicity, and radius of the star. This
analysis of the SED alone yielded the following parameters: $T_{eff}=6080\pm
86\,K$, $r=1.882\pm 0.047\,R_{\odot}$, $d=2.42\pm 0.14\,kpc$, and
$[Fe/H]=-0.21\pm 0.28$. We note that these estimates depend strongly on the
assumed reddening, but are broadly consistent with what we estimate from the
spectra.

Extended Data Figure 6: a) A Hubble Space Telescope ACS/HRC F330W image of
V404 Cygni, with both objects labelled. We used this archival image to extract
flux in the F330W filter, which was included in our analysis of the spectral
energy distribution of the object. b) A Keck NIRC2 adaptive optics image
obtained in the $K_{p}$ filter, used in our analysis of the spectral energy
distribution of the object.

## References

  * [1] Hobbs, G., Lorimer, D. R., Lyne, A. G. & Kramer, M.  A statistical study of 233 pulsar proper motions.  _MNRAS_ 360, 974–992 (2005).  astro-ph/0504584.
  * [2] Tauris, T. M. _et al._ Formation of Double Neutron Star Systems.  _ApJ_ 846, 170 (2017).  1706.09438.
  * [3] Fryer, C. L. & Kalogera, V.  Theoretical Black Hole Mass Distributions.  _ApJ_ 554, 548–560 (2001).  astro-ph/9911312.
  * [4] Burrows, A., Wang, T., Vartanyan, D. & Coleman, M. S. B.  A Theory for Neutron Star and Black Hole Kicks and Induced Spins.  _arXiv e-prints_ arXiv:2311.12109 (2023).  2311.12109.
  * [5] Naoz, S., Fragos, T., Geller, A., Stephan, A. P. & Rasio, F. A.  Formation of Black Hole Low-mass X-Ray Binaries in Hierarchical Triple Systems.  _ApJ_ 822, L24 (2016).  1510.02093.
  * [6] Callister, T. A., Farr, W. M. & Renzo, M.  State of the Field: Binary Black Hole Natal Kicks and Prospects for Isolated Field Formation after GWTC-2.  _ApJ_ 920, 157 (2021).  2011.09570.
  * [7] O’Leary, R. M., Meiron, Y. & Kocsis, B.  Dynamical Formation Signatures of Black Hole Binaries in the First Detected Mergers by LIGO.  _ApJ_ 824, L12 (2016).  1602.02809.
  * [8] Rodriguez, C. L., Haster, C.-J., Chatterjee, S., Kalogera, V. & Rasio, F. A.  Dynamical Formation of the GW150914 Binary Black Hole.  _ApJ_ 824, L8 (2016).  1604.04254.
  * [9] Fryer, C. L.  Mass Limits For Black Hole Formation.  _ApJ_ 522, 413–418 (1999).  astro-ph/9902315.
  * [10] Mirabel, F.  The formation of stellar black holes.  _New A Rev._ 78, 1–15 (2017).
  * [11] Mandel, I.  Estimates of black hole natal kick velocities from observations of low-mass X-ray binaries.  _MNRAS_ 456, 578–581 (2016).  1510.03871.
  * [12] Andrews, J. J. & Kalogera, V.  Constraining Black Hole Natal Kicks with Astrometric Microlensing.  _ApJ_ 930, 159 (2022).  2203.15156.
  * [13] Shenar, T. _et al._ An X-ray-quiet black hole born with a negligible kick in a massive binary within the Large Magellanic Cloud.  _Nature Astronomy_ 6, 1085–1092 (2022).  2207.07675.
  * [14] Fragos, T. _et al._ Understanding Compact Object Formation and Natal Kicks. II. The Case of XTE J1118 + 480.  _ApJ_ 697, 1057–1070 (2009).  0809.1588.
  * [15] Dashwood Brown, C., Gandhi, P. & Zhao, Y.  On the natal kick of the black hole X-ray binary H 1705-250.  _MNRAS_ 527, L82–L87 (2024).
  * [16] Khargharia, J., Froning, C. S. & Robinson, E. L.  Near-infrared Spectroscopy of Low-mass X-ray Binaries: Accretion Disk Contamination and Compact Object Mass Determination in V404 Cyg and Cen X-4.  _ApJ_ 716, 1105–1117 (2010).  1004.5358.
  * [17] Udalski, A. & Kaluzny, J.  CCD Photometry of the X-ray Nova V404 Cygni after the 1989 Outburst.  _PASP_ 103, 198 (1991).
  * [18] Casares, J., Charles, P. A., Jones, D. H. P., Rutten, R. G. M. & Callanan, P. J.  Optical studies of V404 Cyg, the X-ray transient GS 2023+338 -I. The 1989outburst and decline.  _MNRAS_ 250, 712 (1991).
  * [19] Maitra, D. _et al._ Simultaneous Multiwavelength Observations of V404 Cygni during its 2015 June Outburst Decay Strengthen the Case for an Extremely Energetic Jet-base.  _ApJ_ 851, 148 (2017).  1712.06668.
  * [20] Miller-Jones, J. C. A. _et al._ The First Accurate Parallax Distance to a Black Hole.  _ApJ_ 706, L230–L234 (2009).  0910.5253.
  * [21] Martin, E. L., Rebolo, R., Casares, J. & Charles, P. A.  High lithium abundance in the secondary of the black-hole binary system V404 Cygni.  _Nature_ 358, 129–131 (1992).
  * [22] Blaauw, A.  On the origin of the O- and B-type stars with high velocities (the “run-away” stars), and some related problems.  _Bull. Astron. Inst. Netherlands_ 15, 265 (1961).
  * [23] Raghavan, D. _et al._ A Survey of Stellar Families: Multiplicity of Solar-type Stars.  _ApJS_ 190, 1–42 (2010).  1007.0414.
  * [24] Tokovinin, A.  From Binaries to Multiples. II. Hierarchical Multiplicity of F and G Dwarfs.  _AJ_ 147, 87 (2014).  1401.6827.
  * [25] Armas Padilla, M. _et al._ Multiwavelength spectroscopy of the black hole candidate MAXI J1813-095 during its discovery outburst.  _MNRAS_ 485, 5235–5243 (2019).  1903.04498.
  * [26] Mata Sánchez, D. _et al._ Dynamical confirmation of a stellar mass black hole in the transient X-ray dipping binary MAXI J1305-704.  _MNRAS_ 506, 581–594 (2021).  2104.07042.
  * [27] Hynes, R. I., Haswell, C. A., Chaty, S., Shrader, C. R. & Cui, W.  The evolving accretion disc in the black hole X-ray transient XTE J1859+226.  _MNRAS_ 331, 169–179 (2002).  astro-ph/0111333.
  * [28] Miller-Jones, J. C. A. _et al._ A rapidly changing jet orientation in the stellar-mass black-hole system V404 Cygni.  _Nature_ 569, 374–377 (2019).  1906.05400.
  * [29] Liu, B. & Lai, D.  Spin-Orbit Misalignment of Merging Black Hole Binaries with Tertiary Companions.  _ApJ_ 846, L11 (2017).  1706.02309.
  * [30] Su, Y., Lai, D. & Liu, B.  Spin-orbit misalignments in tertiary-induced binary black-hole mergers: Theoretical analysis.  _Phys. Rev. D_ 103, 063040 (2021).  2010.11951.
  * [31] El-Badry, K. & Rix, H.-W.  Imprints of white dwarf recoil in the separation distribution of Gaia wide binaries.  _MNRAS_ 480, 4884–4902 (2018).  1807.06011.
  * [32] Brandt, N. & Podsiadlowski, P.  The effects of high-velocity supernova kicks on the orbital properties and sky distributions of neutron-star binaries.  _MNRAS_ 274, 461–484 (1995).
  * [33] Casares, J. _et al._ Accretion and outflow in V404 Cyg.  _MNRAS_ 488, 1356–1365 (2019).  1907.00005.
  * [34] Freudling, W. _et al._ Automated data reduction workflows for astronomy. The ESO Reflex environment.  _A &A_ 559, A96 (2013).  1311.5411.
  * [35] Husser, T. O. _et al._ A new extensive library of PHOENIX stellar atmospheres and synthetic spectra.  _A &A_ 553, A6 (2013).  1303.5632.
  * [36] Chambers, K. C. _et al._ The Pan-STARRS1 Surveys.  _arXiv e-prints_ arXiv:1612.05560 (2016).  1612.05560.
  * [37] Sirianni, M. _et al._ The Photometric Performance and Calibration of the Hubble Space Telescope Advanced Camera for Surveys.  _PASP_ 117, 1049–1112 (2005).  astro-ph/0507614.
  * [38] Matthews, K. & Soifer, B. T.  The Near Infrared Camera on the W. M. Keck Telescope.  In McLean, I. S. (ed.) _Astronomy with Arrays, The Next Generation_ , vol. 190 of _Astrophysics and Space Science Library_ , 239 (1994).
  * [39] Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S. & Finkbeiner, D.  A 3D Dust Map Based on Gaia, Pan-STARRS 1, and 2MASS.  _ApJ_ 887, 93 (2019).  1905.02734.


# Physics-Consistent Data-driven Waveform Inversion with Adaptive Data
Augmentation

Renán Rojas-Gómez†,⋄, Jihyun Yang†,#, Youzuo Lin†,⋆, James Theiler†, and
Brendt Wohlberg† $\dagger$: Los Alamos National Laboratory
$\diamond$: Department of Electrical and Computer Engineering, University of
Illinois at Urbana-Champaign
$\\#$: Department of Geophysics, Colorado School of Mines
$\star$ Correspondence to: Y. Lin<EMAIL_ADDRESS>
###### Abstract

Seismic full-waveform inversion (FWI) is a nonlinear computational imaging
technique that can provide detailed estimates of subsurface geophysical
properties. Solving the FWI problem can be challenging due to its ill-
posedness and high computational cost. In this work, we develop a new hybrid
computational approach to solve FWI that combines physics-based models with
data-driven methodologies. In particular, we develop a data augmentation
strategy that can not only improve the representativity of the training set,
but also incorporate important governing physics into the training process and
therefore improve the inversion accuracy. To validate the performance, we
apply our method to synthetic elastic seismic waveform data generated from a
subsurface geologic model built on a carbon sequestration site at Kimberlina,
California. We compare our physics-consistent data-driven inversion method to
both purely physics-based and purely data-driven approaches and observe that
our method yields higher accuracy and greater generalization ability.

###### Index Terms:

Computational Imaging, Full-waveform Inversion, Convolutional Neural Networks,
Physics-consistent Machine Learning, Data Augmentation

##  I Introduction

In solid earth geosciences, characterizing the subsurface geology is crucial
for energy exploration, civil infrastructure, groundwater contamination and
remediation, etc. However, nearly all of the earth’s interior is inaccessible
to direct observation. Inference of unknown subsurface properties therefore
relies on indirect and limited geophysical measurements taken at or near the
surface. Seismic inversion attempts to reconstruct an image of subsurface
structures from measurements of natural or artificially produced seismic waves
that have travelled through the subsurface. A forward model describes how the
observations depend on the subsurface map, while the inverse problem involves
inferring that map from the observations. The forward model of wave
propagation is nonlinear. Travel-time inversion methods [14] are based on a
linear approximation of the forward model, while seismic full-waveform
inversion (FWI) addresses the full non-linear problem, leading to superior
inversion accuracy and resolution [15].

The FWI problem is challenging due to the non-linearity of the forward model
and its under-determined nature. Conventional computational methods for
solving FWI are based on optimization techniques and generic regularization
[15]. For simplicity of description, we call these approaches “physics-based
FWI methods” to distinguish them from data-driven methods, and from our
proposed hybrid approach. The major advantage of these methods is their
robustness to out-of-distribution data due to noise, change of station, and
other external factors, while the main disadvantage is their computational
expense. The primary computational cost involved in FWI is associated with the
solution of the wave equation, and is affected by the details of the finite
difference solver, the velocity model dimension, sources and receivers.

Most existing regularization techniques used for solving FWI employ generic
functions, such as $\ell_{1}$ or $\ell_{2}$ penalties on the gradient of the
solution [9].

Recently, a new class of algorithm has been developed, based on machine
learning applied to large datasets that are produced from many runs of the
forward physics model. In direct end-to-end learning [17, 1, 4], a large
number of velocity maps and corresponding seismic waveforms (usually
constructed through extensive simulation) are used as training data in
learning the mapping from seismic waveform to velocity map. In low-wave number
learning [12, 6], this type of learning approach is used to predict an initial
velocity map with low-frequency component, which is then used as the initial
guess for traditional physics-based optimization.

Here, we describe an approach for seismic FWI that incorporates the physics
model into the learning procedure. Specifically, this physics-consistent data-
driven full waveform inversion consists of a carefully designed encoder-
decoder-structured neural network and an adaptive data augmentation technique.
This augmentation employs the forward model to produce new training data that
are more representative of the solution we seek. To validate its performance,
we applied our inversion method to detect carbon sequestration leakage using
synthetic seismic data sets generated using a subsurface model for a potential
CO2 storage site at Kimberlina, California [3].

##  II Background

###  II-A Governing Physics: the Forward Model

Mathematically, the forward model can be expressed in terms of the seismic
elastic-wave partial differential equation [15]:

$\displaystyle\rho({r})\frac{\partial^{2}u({r},t)}{\partial t^{2}}=$
$\displaystyle(\lambda({r})+\mu({r}))\nabla(\nabla\cdot u({r},t))$
$\displaystyle+\mu({r})\nabla^{2}u({r},t)+s({r},\,t),$ (1)

where $\rho({r})$ is the density at spatial location ${r}$, $\lambda({r})$ and
$\mu({r})$ are the Lamé parameters, $s({r},\,t)$ is the source term,
$u(\mathbf{r},t)$ is the displacement wavefield, $t$ represents time, and
$\cdot$ is the divergence operator. When fluid such as supercritical CO2 leaks
into the subsurface formation, the geophysical parameters of P-wave and S-wave
velocities will be changed correspondingly.

Instead of inverting for $\rho({r})$, $\lambda({r})$ and $\mu({r})$, it is
customary to invert for a velocity map $m\in\mathbb{R}^{M\times N}$, where $M$
and $N$ are its vertical and lateral dimensions, respectively; here $m$ refers
to either P-wave or S-wave velocity can be expressed as a function of
$\rho({r})$, $\lambda({r})$ and $\mu({r})$. Similarly, we denote a seismic
data observation $d_{\text{obs}}\in\mathbb{R}^{T\times S\times R}$, where $T$
corresponds to the number of samples in the temporal domain, $S$ to the number
of sources and $R$ to the number of receivers used in the data acquisition
process. The seismic data can be expressed in terms of a highly nonlinear
forward mapping $f$:

$d_{\text{obs}}=f({m}).$ (2)

###  II-B Physics-Based Full-waveform Inversion

Various explicit regularization techniques have been developed to stabilize
the computation of seismic inversion, including $\ell_{1}$-norm [9] or
$\ell_{2}$-norm [5] methods. Given the forward model in Eq. (2), the
regularized seismic FWI can be posed as

$m=\underset{{m}}{\operatorname{argmin}}\left\\{\left\|d-f({m})\right\|_{2}^{2}+\lambda\,R({m})\right\\},$
(3)

where ${d}$ represents a recorded/field waveform dataset, $f({m})$ is the
corresponding forward modeling result, $\left\|d-f({m})\right\|_{2}^{2}$ is
the data misfit, $||\cdot||_{2}$ stands for the $\ell_{2}$ norm, $\lambda$ is
a regularization parameter and $R({m})$ is the regularization term. Note that
Eq. (3) is a general formulation for regularized FWI. More effective
regularization techniques have been developed for time-lapse monitoring [10].

##  III Physics-Consistent Data-Driven Full-waveform Inversion

###  III-A Data-Driven Inversion and Network Structure

A data-driven FWI structure based on an encoder-decoder architecture [16],
denoted $\mathcal{G}$ and characterized by hyperparameters
$\boldsymbol{\theta}$, is proposed to approximate the inverse mapping $f^{-1}$
and obtain accurate velocity map predictions
$\hat{m}(\boldsymbol{\theta})\triangleq\mathcal{G}(\boldsymbol{\theta},d_{\text{obs}})$
under a supervised learning scheme. Optimal parameters
$\boldsymbol{\theta}^{*}$ are obtained by adapting the architecture to a
representative training set with $L$ samples
$\\{d_{\text{obs},\ell},m_{\ell}^{*}\\},\ell\in\\{0,L-1\\}$.

We choose the mean-absolute error (MAE) as our optimality criterion:

$\displaystyle\boldsymbol{\theta}^{*}=\underset{\boldsymbol{\theta}}{\text{argmin}}\
\frac{1}{L}\sum_{\ell=0}^{L-1}\|m_{\ell}^{*}-\hat{m}_{\ell}(\boldsymbol{\theta})\|_{1}.$
(4)

For a more detailed discussion of loss function selection, please refer to our
earlier work [16].

Our data-driven inversion network structure consists of an encoder network and
a decoder network [16]. Full details of our model are provided in the
Supplementary Material.

###  III-B Data Description

We apply our method to detect CO2 leakage in subsurface. To the best of our
knowledge, there is no real seismic data related to our problem of interest.
We use the simulated Kimberlina dataset from Lawrence Livermore National
Laboratory (and that we refer to here as CO2leak). The aim of the Kimberlina
dataset is to understand and assess the effectiveness of various geophysical
monitoring techniques in detecting CO2 shallow leakage in the wellbore [7]. A
portion of the of Kimberlina Simulations can be be downloaded from the DOE-EDX
platform [11]. The Kimberlina dataset is generated based on a hypothetical
numerical model built on the geologic structure of a commercial-scale geologic
carbon sequestration (GCS) reservoir at the Kimberlina site in the southern
San Joaquin Basin, 30 km northwest of Bakersfield, CA, USA. The P-wave and
S-wave velocity maps used in this work belong to the geophysical model, which
is created based on the realistic geologic-layer properties from the GCS site
[3].

The CO2leak dataset contains 991 CO2 leakage scenarios, each simulated over a
duration of 200 years, with 20 leakage maps provided (ie, at every ten years)
for each scenario. We obtain synthetic seismograms from elastic forward
modeling on CO2leak velocity maps. First, one-second traces with a time
interval of $0.5$ms using 7 sources and 114 receivers are generated. We then
down-sample each trace by a factor of $2$, resulting in a temporal dimension
of $1000$ time steps. The sources and receivers are evenly distributed along
the top of the model, with depths of 5m and 20m, respectively. The source
interval is 125m, and the receiver interval is 15m. We use a Ricker wavelet
with a central frequency of $25$Hz as the source to generate simulated seismic
waves due to its empirical success in processing seismic field data [5]. The
synthetic data is the staggered grid solution of the elastic wave equation
using a finite-difference scheme with a perfectly matched layered absorbing
boundary condition.

###  III-C Data Augmentation: Incorporation of Physics Knowledge

The CO2leak dataset includes $19,600$ velocity maps of $141\times 341$ grid
points describing CO2 and brine leakage plumes evolving with time. It is of
practical interest to detect plumes of leaking CO2 while they are still small.
This is particularly challenging when the available training data is dominated
by large plumes. Thus, CO2leak presents the opportunity to evaluate the
generalization of data-driven inversion with respect to different plume sizes.

Along with the data pairs included in the dataset, the ground-truth CO2 and
brine mass information for each sample, is provided. Based on this, the full
dataset is split into four parts, according to their CO2 leak mass plus brine
leak mass: tiny plumes (from $3.53\times 10^{2}$ to $9.10\times 10^{6}$Kg),
small plumes (from $9.10\times 10^{6}$ to $2.67\times 10^{7}$Kg), medium
plumes (from $2.67\times 10^{7}$ to $8.05\times 10^{7}$Kg), and large plumes
(from $8.05\times 10^{7}$ to $1.62\times 10^{9}$ Kg). These cover $20\%$,
$20\%$, $20\%$, and $40\%$ of the data samples, respectively. Figure 5 in the
Supplementary Material shows representative labels from each subset.

While conventional data augmentation techniques (such as rotation, flip, scale
etc.) have proved to be effective for image processing applications, it is not
clear that they have a useful role to play in our application. Our adaptive
data augmentation scheme provides additional training data that is not only
physically meaningful but also more closely related to the target unlabelled
data that we are trying to invert. We summarize our augmentation method as the
following four steps (illustrated in Fig. 1; a detailed description is
provided in the Supplementary Material):

  1. i.

Estimate approximate solver
$\mathcal{G}(\boldsymbol{\hat{\theta}},d_{\text{obs}})$;

  2. ii.

Generate approximate velocity maps from unlabeled data
$\hat{m}_{r}=\mathcal{G}(\boldsymbol{\hat{\theta}},d_{\text{obs},r})$;

  3. iii.

Create seismic data using forward model
$\hat{d}_{\text{obs},r}=f(\hat{m}_{r})$;

  4. iv.

Add new pairs to the original training set.

The augmented dataset plays a key role in model accuracy because it will not
only carry useful physics information, but also provides examples of velocity
maps that are consistent with the target geology feature of interests.
Furthermore, the full augmentation process can be applied in an iterative
fashion by re-training the approximate solver
$\mathcal{G}(\boldsymbol{\hat{\theta}},d_{\text{obs}})$ based on the extended
training set in order to generate new approximate velocity maps $\hat{m}_{r}$.
This approach allows further refinement of the mapping between velocity and
seismic subdomains, as empirically shown in Section IV-C.

Figure 1: Adaptive Data Augmentation: Approximate solver
$\mathcal{G}(\hat{\theta},d_{\text{obs}})$ is fully-trained over labeled set
$\\{m^{*}_{\ell},d_{\text{obs},\ell}\\}$, and applied to unlabeled seismic
data $d_{\text{obs},r}$ to generate new velocity maps $\hat{m}_{r}$.
Physically-coherent seismic data $\hat{d}_{\text{obs},r}$ is then generated
using the forward model $f$, producing a new labeled set
$\\{\hat{m}_{r},\hat{d}_{\text{obs},r}\\}$ which is added to the original
training set.

##  IV Numerical Experiments

###  IV-A Experimental Setup

For all evaluations, the training process is performed using a fixed stopping
criterion of $250$ epochs, random weight initialization, an initial learning
rate of $10^{-3}$ and its subsequent adaptive optimization via ADAM [8]. We
found that $250$ epochs is appropriate to guarantee diversity in the
reconstructed velocity maps while avoiding overfitting. The training is
performed using fixed batch sizes of $50$, where both features and labels are
normalized prior to their use in the training and inference processes. We
created $392$ batches in total. In consequence, reported loss values are
computed based on normalized training and testing pairs. All training and
testing routines are implemented using PyTorch code and executed on four
NVIDIA GeForce GTX 1080 GPUs. We mention that there are various inversion
strategies for elastic FWI [13], and we provide the inversion of P-wave
velocities in all tests considering its higher relative sensitivity to the
response of CO2 leakage. The S-wave velocities can be obtained through either
sequential inversion or empirical relations [13].

###  IV-B Generalization Performance without Data Augmentation

We assess network generalization performance by first training over one subset
of the training data (e.g., large plumes) and testing on a different subset
(e.g., tiny plumes). For each case, we quantify how informative is our data
augmentation approach by measuring the Mean Absolute Error (MAE), as expressed
in Eq. 4, attained by the network along epochs. We show the reconstruction
accuracy along epochs $\varepsilon(\boldsymbol{\theta}_{i})$ as the MAE
normalized with respect to the velocity map dimensions to clearly depict our
method’s effect in the predictions:

$\displaystyle\varepsilon(\boldsymbol{\theta}_{i})=$
$\displaystyle\frac{1}{\hat{L}MN}\sum_{\ell=0}^{\hat{L}-1}\|m_{\ell}^{*}-\hat{m}_{\ell}(\boldsymbol{\theta}_{i})\|_{1},$
(5)

where $\boldsymbol{\theta}_{i}$ corresponds to the set of hyperparameters at
the $i^{\text{th}}$ epoch, and $\hat{L}$ corresponds to the size of the
validation set.

Figure 2 shows the testing loss curve at each epoch. Although decreasing, the
curve shows that the network is not able to accurately reconstruct the
velocity maps under this scenario: given the reduced plume size in the testing
set, an $\ell_{1}$ loss value of $-8.7$ indicates either a large intensity
deviation among ground-truth and prediction, or a very low intersection over
the union between them. This is reflected in Figures 2 and 2, which show a
tiny ground-truth velocity map and its prediction by the network trained over
large plumes. In general, the network tries to generate plumes of large size,
which is consistent with the plumes on which it was trained.

(a) Testing loss curve.

(b) Tiny label.

(c) Tiny prediction.

Figure 2: Transferability: Training on the large subset, testing on the tiny
subset (0.205 (MAE) and 0.899 (SSIM)).

Figure 3: Reconstruction results for our data Augmentation approach applied to
labeled data: Testing MAE loss value for each training epoch.

Figure 4: Reconstruction results for our data Augmentation approach applied to
unlabeled data: Testing MAE loss value for each training epoch.

(a) (b) (c) (d) (e)

Figure 5: (a) Ground truth. Inversions and errors obtained using (b) physics-
based FWI (0.069, 0.972), data-driven model trained on (c) the large and
medium subsets (0.0620, 0.990), (d) large and augmented medium subsets
(augmented once) (0.0134, 0.993), and (e) large and augmented medium subsets
(augmented twice) (0.0122, 0.994). Error is in the format of (MAE, SSIM).

###  IV-C Generalization Performance with Data Augmentation

We evaluate the reconstruction accuracy of the augmentation method described
in Section III-C, and compare to traditional data-driven methods with no
augmentation strategies.

For the first scenario, considering the data partitioning described in Section
III-C, we compute the reconstruction accuracy along epochs for the following
setups:

  1. i.

Train on the large subset.

  2. ii.

Train on the large and medium subsets.

  3. iii.

Train on the large, medium and small subsets.

  4. iv.

Train on the large and augmented medium subset.

In every case, performance is evaluated on the tiny subset.

For setups (i), (ii) and (iii), the main differences are the number of
training samples and the size of the plumes included in each training set.
Based on the network transferability results shown in Section IV-C, the
network trained over large, medium and small plumes is expected to attain the
best reconstruction accuracy on tiny plumes. On the other hand, setup (iv),
corresponding to the network trained over large and augmented medium subsets,
utilizes the same information as setup (ii). However, setup (ii) also includes
medium samples. In other words, the difference between setup (ii) and (iv) is
that setup (iv) increases the number of medium samples using our augmentation
strategy.

Figure 3 shows the reconstruction accuracy along epochs for the four described
setups. Regarding the training scenarios with no augmentation strategies, the
Mean Absolute Error along epochs behaves as expected: after 250 epochs,
training over large, medium and small plumes attains the best reconstruction
value of $-13.5$ for the error metric in Eq. 5, followed by the network
trained over large and medium plumes at $-9.3$, and the network trained over
large plumes at $-8.7$. On the other hand, our augmented approach from setup
(iv) reaches an accuracy of approximately $-12$, which is smaller than the
error obtained by setup (ii). This implies that the augmented medium samples
generated using the forward model provide physically-coherent information to
the training process, allowing a better velocity map reconstruction.

It is important to remark that the augmentation approach from setup (iv)
generates new samples by augmenting labeled data from the medium subset.
Specifically, we train a surrogate network over the labeled large subset, use
it to generate new data pairs over the labeled medium dataset, and then fully
train the network using the large subset along with both original and
augmented medium datasets. Although this experimental setup and its comparison
against non-augmented strategies illustrates the advantages of our proposed
augmentation approach, a more interesting scenario is the generation of
physically-consistent pairs from unlabeled data.

With this in mind, a second scenario is considered. Let the medium unlabeled
subset correspond exclusively to the medium seismic data (the corresponding
velocity maps are not taken into account). Then, consider the following two
experiments:

  1. i.

Train on the large subset.

  2. ii.

Train on the large and augmented medium unlabeled subset.

  3. iii.

Test on the small subset.

The augmentation process follows the procedure described in Section III-C. In
contrast with our first experiment, however, in which both medium and
augmented medium datasets are added to the large dataset, only the augmented
medium dataset is added to the initial large dataset. Figure 4 shows
reconstruction results along $250$ epochs. The network trained over large
plumes attain a reconstruction accuracy of approximately $-10.5$. On the other
hand, the network trained over large and augmented medium unlabeled plumes
obtain a better reconstruction accuracy of approximately $-11.7$. These
results strengthen our observation regarding the information encapsulated in
the samples generated by the augmentation process: by including the forward
modeling operation in the process, physically-consistent data pairs are
generated, allowing a better domain adaptation. The use of unlabeled data
shows how our proposed method is not limited to labeled data, which can be
difficult to obtain in real applications.

Finally, we provide visualization of the inverted velocity maps using both
physics-based FWI and data-driven approaches in Fig. 5. Our network trained
over the large and medium subsets obtains a reasonably accurate estimate of
the tiny plume. Also, the estimate obtained by our augmentation further
refines the plume shape and location. This example shows how our augmentation
approach used iteratively improves the network output, which reflects the
potential of including the physics-based forward model in learning. An
additional test and detailed analysis of our inversion model robustness with
respect to variant levels of noisy seismic data is provided in the
Supplementary Material.

##  V Conclusion and Future Work

We develop a physics-consistent data-driven seismic full-waveform inversion
method. We design a novel data augmentation strategy that incorporates
critical physics information and improves the representability of the training
set. We validate its performance to detect small CO2 leakage. Compared with
purely physics-based and purely data-driven inversion methods, our physics-
consistent data-driven inversion yields higher accuracy and better
generalization. With respect to computational cost, the most expensive
component of our approach is in data generation (including modeling) and
training. The cost of the inference is very low. Considering a single
inversion, the overall modeling cost of our method can be more expensive (up
to 10 times) than that of physics-based FWI. However, for those applications
that require inverting multiple seismic surveys at the same location such as
time-lapse monitoring, our approach can be very cost-effective.

Nonrepeatability is an important issue for time-lapse seismic imaging [2]. Our
technique is general, and can be combined with existing time-lapse imaging
methods. We will study the robustness of our technique in nonrepeatable
scenarios. Different means have been proposed to incorporate physics in
solving inverse problems, cyclic consistency being one of those [18]. We will
compare the performance between our techniques and consistent generative
methods. The high cost in training may hinder the wide application of our
method to a broader class of problems. We believe, however, that there may be
room for further improving the efficiency of the training. Another future
direction would be to explore data augmentation in the low-data regime.

## Acknowledgements

This work was supported by the Center for Space and Earth Science at Los
Alamos National Laboratory (LANL) and by the Laboratory Directed Research and
Development program of LANL under project number 20200061DR. Renán A. Rojas-
Gómez would also like to thank Javier E. Santos and Manish Bhattarai for the
insightful discussions. We also thank two anonymous reviewers and the
Associate Editor, Dr. Luis Gómez, for their constructive suggestions and
comments that improved the quality of this work.

## References

  * Araya-Polo et al., [2018] Araya-Polo, M., J. Jennings, A. Adler, and T. Dahlke, 2018, Deep-learning tomography: The Leading Edge, 37, 58–66.
  * Asnaashari et al., [2015] Asnaashari, A., R. Brossier, S. Garambois, F. Audebert, P. Thore, and J. Virieux, 2015, Time-lapse seismic imaging using regularized full-waveform inversion with a prior model: which strategy?: Geophysical Prospecting, 63, 78–98.
  * Buscheck et al., [2019] Buscheck, T., K. Mansoor, X. Yang, H. Wainwright, and S. Carroll, 2019, Downhole pressure and chemical monitoring for CO2 and brine leak detection in aquifers above a CO2 storage reservoir: Int. J. Greenhouse Gas Control, 91.
  * Farris et al., [2018] Farris, S., M. Araya-Polo, J. Jennings, B. Clapp, and B. Biondi, 2018, Tomography: a deep learning vs full-waveform inversion comparison, in First EAGE Workshop on High Performance Computing for Upstream in Latin America, European Association of Geoscientists & Engineers.
  * Fichtner, [2010] Fichtner, A., 2010, Full seismic waveform modelling and inversion: Springer Science & Business Media.
  * Hu et al., [2019] Hu, W., Y. Jin, X. Wu, and J. Chen, 2019, Progressive transfer learning for low frequency data prediction in full waveform inversion: arXiv preprint arXiv:1912.09944.
  * Jordan and Wagoner, [2017] Jordan, P., and J. Wagoner, 2017, Characterizing construction of existing wells to a CO2 storage target: The Kimberlina site, California: Technical report, U.S. Department of Energy - Office of Fossil Energy.
  * Kingma and Ba, [2014] Kingma, D. P., and J. Ba, 2014, Adam: A method for stochastic optimization: arXiv preprint arXiv:1412.6980.
  * Lin and Huang, [2015] Lin, Y., and L. Huang, 2015, Acoustic- and elastic-waveform inversion using a modified total-variation regularization scheme: Geophysical Journal International, 200, 489–502.
  * Maharramov et al., [2016] Maharramov, M., B. Biondi, and M. Meadows, 2016, Time-lapse inverse theory with applications: Geophysics, 81, R485–R501.
  * NETL, [2018] NETL, 2018, LLNL Kimberlina 1.2 Simulations.  (https://edx.netl.doe.gov/dataset/llnl-kimberlina-1-2-nuft-simulations-june-2018-v2).
  * Ovcharenko et al., [2019] Ovcharenko, O., V. Kazei, M. Kalita, D. Peter, and T. Alkhalifah, 2019, Deep learning for low-frequency extrapolation from multioffset seismic data: Geophysics, 84, no. 6, 58–66.
  * Raknes and Arntsen, [2014] Raknes, E. B., and B. Arntsen, 2014, Strategies for elastic full waveform inversion, in SEG Technical Program Expanded Abstracts: Society of Exploration Geophysicists, 1222–1226.
  * Tarantola, [2005] Tarantola, A., 2005, Inverse problem theory and methods for model parameter estimation: SIAM.
  * Virieux et al., [2014] Virieux, J., A. Asnaashari, R. Brossier, L. Métivier, A. Ribodetti, and W. Zhou, 2014, Chapter 6: An introduction to full waveform inversion, in Encyclopedia of Exploration Geophysics: SEG.
  * Wu and Lin, [2019] Wu, Y., and Y. Lin, 2019, InversionNet: An efficient and accurate data-driven full waveform inversion: IEEE Transactions on Computational Imaging, 6, 419–433.
  * Yang and Ma, [2019] Yang, F., and J. Ma, 2019, Deep-learning inversion: A next-generation seismic velocity model building method: Geophysics, 84, no. 4, R583–R599.
  * Zhu et al., [2017] Zhu, J., T. Park, P. Isola, and A. Efros, 2017, Unpaired image-to-image translation using cycle-consistent adversarial networks: Presented at the IEEE International Conference on Computer Vision.