# A Contact Binary Mis-Classified as an Ellipsoidal Variable: Complications
for Detached Black Hole Searches

Tyrone N. O’Doherty International Centre for Radio Astronomy Research – Curtin
University, GPO Box U1987, Perth, WA 6845, Australia  Arash Bahramian
International Centre for Radio Astronomy Research – Curtin University, GPO Box
U1987, Perth, WA 6845, Australia  Adelle J. Goodwin International Centre for
Radio Astronomy Research – Curtin University, GPO Box U1987, Perth, WA 6845,
Australia  James C. A. Miller-Jones International Centre for Radio Astronomy
Research – Curtin University, GPO Box U1987, Perth, WA 6845, Australia  Jerome
A. Orosz Villanova University, Department of Astrophysics and Planetary
Sciences, 800 East Lancaster Avenue, Villanova, PA 19085, US  Jay Strader
Center for Data Intensive and Time Domain Astronomy, Department of Physics and
Astronomy, Michigan State University, East Lansing, MI 48824, USA

###### Abstract

Identifying sources exhibiting ellipsoidal variability in large photometric
surveys is becoming a promising method to search for candidate detached black
holes in binaries. This technique aims to exploit the orbital-phase dependent
modulation in optical photometry caused by the black hole distorting the shape
of the luminous star to constrain the mass ratio of the binary. Without
understanding if, or how much, contamination is present in the candidate black
hole samples produced by this new technique it is hard to leverage them for
black hole discovery. Here, we follow up one of the best candidates identified
from Gaia Data Release 3, Gaia DR3 4042390512917208960, with a radial velocity
campaign. Combined photometric and radial velocity modelling, along with
spectral disentangling, suggests that the true mass ratio (mass of the unseen
object divided by the mass of the luminous star) is an order of magnitude
smaller than that inferred assuming the modulations arise from ellipsoidal
variability. We therefore infer that this system is likely a contact binary,
or on the boundary of both stars nearly filling their Roche lobes, however,
further observations are required to confidently detect the secondary. We find
that the well-known problem of discriminating between ellipsoidal and contact
binary light curves results in a larger contamination from contact binaries
than previously suggested. Until ellipsoidal variables can be reliably
distinguished from contact binaries, samples of black hole candidates selected
based on ellipsoidal variability are likely to be highly contaminated by
contact binaries or similar systems.

††facilities: Gaia, ATT, OGLE, SOAR††software: aladin (Bonnarel et al., 2000;
Boch & Fernique, 2014), Astropy (Astropy Collaboration et al., 2022), emcee
(Foreman-Mackey et al., 2013), exoplanet (Foreman-Mackey et al., 2021), ELC
(Orosz & Hauschildt, 2000), iPython (Perez & Granger, 2007), IRAF/fxcor (Tody,
1986, 1993), Matplotlib (Hunter, 2007), Numpy (Harris et al., 2020), Pandas
(McKinney, 2010), PHOEBE (Prša & Zwitter, 2005), PyMC3 (Salvatier et al.,
2016), SAO DS9 (Joye & Mandel, 2003), The Joker (Price-Whelan et al., 2017),
Scipy (Virtanen et al., 2020).

##  1 Introduction

The Milky Way Galaxy is expected to contain $10^{8}-10^{9}$ stellar-mass black
holes (BHs) based on current stellar evolution theories (e.g., van den Heuvel
1992; Brown & Bethe 1994; Timmes et al. 1996; Samland 1998; Wiktorowicz et al.
2019; Olejak et al. 2020). However, it is challenging to constrain this
population observationally. So far, the only method that has proven successful
in identifying isolated BHs is through gravitational microlensing, the first
of which was confirmed only recently (Lam et al., 2022a, b; Sahu et al.,
2022). The occurrence of microlensing by isolated BHs is unpredictable, and
these events cannot be followed up once the event is over. As a result,
studies often focus on looking for BHs in binary systems. Until very recently,
the vast majority of confirmed and candidate Galactic BHs came from X-ray
binaries (XRBs) (e.g., Tetarenko et al. 2016; Corral-Santana et al. 2016).
XRBs are binary systems in which a BH accretes from a stellar companion, and
in the process radiates liberated gravitational potential energy. These XRBs
are primarily identified through X-ray emission arising from accretion.
However, broadly, this means the discovery of BHs has been limited to systems
in a configuration favourable to accretion. Furthermore, the discovery of
transient accreting BHs, which make up the majority, is limited to systems
with short recurrence time scales (X-ray all-sky surveys have only been
running for several decades). Whilst accreting BHs make up the majority of the
known population of Galactic BHs, it is believed that they are only a small
component of the total Galactic BH population (e.g., van den Heuvel 1992;
Corral-Santana et al. 2016).

Identifying detached BHs in binary systems in the Galaxy has been the focus of
numerous studies over the last decade. Of the many initially promising
candidates, very few have been found to actually host a BH. While there are
$\sim$ 20 accreting Galactic BHs with dynamical mass estimates (Tetarenko et
al., 2016; Corral-Santana et al., 2016), there are only 5 confirmed Galactic
detached systems with dynamical mass estimates that cannot accommodate a
neutron star (NS) in configurations where a non-degenerate companion is
infeasible. These are the two BHs discovered in the globular cluster NGC 3201
by a spectral survey (Giesers et al., 2018, 2019), two BHs in the Galactic
field (Shahaf et al., 2023; El-Badry et al., 2023a; Tanikawa et al., 2022; El-
Badry et al., 2023b) initially identified using binary astrometric solutions
from Gaia Data Release 3 (DR3; Gaia Collaboration et al. 2016, 2022), and one
system discovered through a spectroscopic survey of Galactic O-type stars
(Mahy et al., 2022). A sixth possible detached BH is 2MASS J05215658+4359220
(Thompson et al., 2019, 2020), although while a BH seems likely, the lower
limit of the mass estimate ($M=3.3^{+2.8}_{-0.7}$ M⊙; Thompson et al. 2020)
does not conclusively rule out a NS due to uncertainties in the NS equation of
state. A possible seventh detached BH is NGC 1850 BH1 (Saracino et al., 2022;
El-Badry & Burdge, 2022; Stevance et al., 2022; Saracino et al., 2023),
however, work is still ongoing. This means that in the Galactic field there
are only 3, possibly 4, known detached BHs thus far. Developing methods of
reliably identifying new BH candidates is critical for discovering more of
this elusive population.

The all-sky astrometric, photometric, and spectroscopic Gaia mission is
providing a wealth of data in which to search for detached black holes.
Indeed, two of the confirmed detached BHs in the Galactic field come from an
analysis of the astrometric orbits provided with DR3. However, while there are
$\sim 2$ billion sources in DR3, there are only $\sim 1$ million binary
orbital solutions. There are likely many more sources with detached BH
companions that have been observed by Gaia that are yet to be identified.
While novel techniques have been used to identify promising candidates from
Gaia (e.g. Andrew et al. 2022; Gomel et al. 2023) it is not understood if, or
how much, contamination from non-BH systems is present in these samples.
Without understanding this possible contamination, it is hard to exploit these
new techniques to their full potential. Following up candidates identified
using Gaia with spectroscopic radial velocity (RV) studies has the potential
to substantially increase the known sample of detached BH systems, as well as
helping to understand the contamination present in these samples and thereby
refine the selection techniques.

Accompanying the release of DR3 were papers detailing the analysis of numerous
kinds of variable sources and their classification (see e.g. Eyer et al. 2023,
and references therein; Rimoldini et al. 2023, and references therein; Gavras
et al. 2023, and references therein). This includes ellipsoidal variables,
which were used to identify candidate binaries with compact companions (Gomel
et al., 2023). In binary systems, the gravity of each component acts on the
other, distorting their shape. In BH-star binaries the BH can significantly
distort the star into a teardrop shape, resulting in sinusoidal modulations of
the observed luminosity of the star on the orbital period. Gomel et al.
(2021a) determined a relationship for a ‘modified minimum mass ratio’ (mMMR)
that can be computed from the light curve modulation, and does not depend on
the mass and radius of the visible star. This mMMR was found to always be less
than the minimum mass ratio, which is in turn less than the true mass ratio
($q=M_{\mathrm{BH}}/M_{\mathrm{star}}$). The main physical assumptions
underlying the derived relationship are that the primary star fills its Roche
Lobe and that the binary has an inclination of
90$\mathrm{\SIUnitSymbolDegree}$ (deviations from these assumptions imply the
true mass ratio is larger than the mMMR). The other main assumption is that
the light is coming from the primary star only. This method was applied to the
variable sources identified as ellipsoidal variables by the Gaia Variability
Pipeline (Eyer et al., 2023) classifier (Rimoldini et al., 2023), producing
6306 detached or weakly-accreting BH and NS candidates.

The 6306 candidates identified by Gomel et al. (2023) are a large sample, and
if all are real, would represent a substantial increase to the known
populations of BHs and NSs. However, as there is likely some level of
contamination from systems like contact binaries (Gomel et al., 2023), follow
up is required. This work builds on that carried out by Nagarajan et al.
(2023), who followed up 14 systems from Gomel et al. (2023) and found they
were unlikely to contain BHs, but did not conclusively characterise the
systems. In this paper, we investigate one of the most promising BH candidates
from the Gomel et al. (2023) sample and discuss the results of our
spectroscopic and photometric modelling that reveal it is likely a contact
binary. In Section 2.1 we discuss the selection of this source for targeted
follow-up and the details of the spectroscopic campaign. Details of the RV
extraction and modelling, joint photometric and RV modelling, spectral
disentangling, and the associated results of the analyses are in Section 3.
The results are discussed in Section 4.

##  2 Observations

###  2.1 Source Selection

We began by filtering the 6306 candidates from Gomel et al. (2023), retaining
sources for which we could estimate a reliable distance and with no evidence
of a luminous companion in Gaia. These restrictions manifested as a cut on
parallax significance greater than 3 and requiring the Gaia parameter
ipd_multi_frac = 0 (the fraction of successful windows in which their fitting
algorithm identified a double peak, and thus potentially a resolved double
star, be it a true binary or a visual double). From here, we began exploration
by restricting mMMR to be at least $2$, leaving a sample of 13 sources. From
this sample of thirteen, there is one particularly interesting source, Gaia
DR3 4042390512917208960.

Gaia DR3 4042390512917208960 has a mMMR of $2.54$ ($1.78$ at one sigma). The
mass estimated by the Gaia team for the luminous star by leveraging
astrometry, photometry, and spectroscopy is $1.809^{+0.058}_{-0.054}$ M⊙
(Creevey et al., 2023; Fouesneau et al., 2023). Assuming this mass is correct,
the inferred mass of the unseen object is $>4.6$ M⊙ ($>3$ M⊙ at one sigma),
making it an excellent BH candidate. This source is one of 262 sources with
mMMR significantly higher than 1, which were identified by Gomel et al. (2023)
as being the most promising BH candidates. Furthermore, it is bright
($G=13.8$), making spectroscopic follow up feasible with a wide range of
instruments.

Gaia DR3 4042390512917208960 has excellent astrometry (Lindegren et al.,
2021), with a parallax significance greater than 20. The renormalised unit
weight error (RUWE) for this source is 0.707, suggestive of over-fitting.
Other astrometric goodness-of-fit statistics (e.g. astrometic excess noise,
and all image parameter determination statistics) are as expected for a well-
behaved system. A RUWE of 0.707 is therefore likely not a cause for concern.
The blue and red photometric excess factor is 1.235, however, this excess is
fully consistent with the excess seen in standard sources (figure 18 in Riello
et al. 2021). The Gaia $G$, $G_{\mathrm{RP}}$, and $G_{\mathrm{BP}}$ phase
folded photometry can be found in Figure 1. The photometry was extracted from
the Gaia archive111https://gea.esac.esa.int/archive/. Sinusoidal modulation
can clearly be seen in all colours. However, the modulation is not completely
smooth, showing some scatter. We think this is likely a result of
contamination, which is not unexpected given the source’s location in the
bulge. Approximately 20% of the source’s transits of the Gaia field of view
were marked as blended, ((phot_bp_n_blended_transits
+phot_rp_n_blended_transits)/(phot_rp_n_obs \+ phot_bp_n_obs)), where a
blended transit means that there was at least one more source in the observing
window (Riello et al., 2021). It is important to note that a blended transit
occurring is dependent on the scan direction, and thus different transits of
the same source can be blended in different ways or not at all (Riello et al.,
2021). However, the scatter does not significantly alter the general shape of
the sinusoidal modulation. In summary, the Gaia astrometric and photometric
solution is believable and raises no major concerns.

Figure 1: Phase folded Gaia $G$, $G_{\mathrm{RP}}$, and $G_{\mathrm{BP}}$
photometry of Gaia DR3 4042390512917208960. The smooth sinusoidal modulation
is evident in all colours.

It is worth noting that this source lies slightly above the main sequence (see
figure 6 in Gomel et al. 2023), which could arise from the unseen companion
being a luminous object. It is possible to hide a more massive stellar
secondary in systems with evolved stars (e.g., El-Badry et al. 2022). However,
the luminous star in Gaia DR3 4042390512917208960 shows no evidence of being
significantly evolved off the main sequence (Fouesneau et al., 2023).
Furthermore, the estimated mass of the unseen object is large enough such that
its luminosity contribution (if it is a main sequence star) would be
significant. Cumulatively, as also identified by Gomel et al. (2023), we
identified this source as an interesting candidate and worthy of further
study.

###  2.2 Archival Observations of Gaia DR3 4042390512917208960

The location and brightness of Gaia DR3 4042390512917208960 means it has been
observed by both the Optical Gravitational Lensing Experiment (OGLE; Udalski
et al. 2015) and the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee
et al. 2014; Kochanek et al. 2017). The OGLE I-band photometry was extracted
from the OGLE Collection of Variable Stars online data
base222https://ogledb.astrouw.edu.pl/~ogle/OCVS/ecl_query.php, and is
presented in Figure 2. It was classified as an ellipsoidal variable by OGLE
(OGLE BLG-ELL-012306; Soszyński et al. 2016). In their study of OGLE
ellipsoidal variables in the Galactic Bulge, Gomel et al. (2021b) derived a
similar mMMR, also noting it as an interesting candidate. However, the ASAS-SN
V team’s machine learning classifier classed it as an eclipsing contact binary
(ASASSN-V J175613.02-335233.3; Jayasinghe et al. 2020). Note that ellipsoidal
variable was not a classification category in Jayasinghe et al. (2020). Later,
Rowan et al. (2021) searched for ellipsoidal variables within ASAS-SN data
using a combined $\chi^{2}$ ratio test followed by visual inspection,
identifying 369 candidates. However, ASASSN-V J175613.02-335233.3 was not
identified as an ellipsoidal candidate in this study. Gomel et al. (2023)
suggested that it is an ellipsoidal variable based on the Gaia and OGLE
classification. They also suggested that the comparatively low amplitude of
the modulation observed by ASAS-SN arises from contamination as the source is
in the Galactic bulge.

Figure 2: Phase folded OGLE I-band photometry. The same shape can be seen as
in Figure 1.

###  2.3 Dedicated Spectroscopic follow-up

We obtained spectroscopic observations of Gaia DR3 4042390512917208960 using
the Wide Field Spectrograph (WiFeS; Dopita et al. 2007, 2010), an integral
field unit (IFU) spectrograph on the Australian National University 2.3 m
optical telescope. The spectrograph has both a ‘blue’ and a ‘red’ camera, with
which we used the ‘B7000’ and ‘R7000’ gratings ($R\sim 7000$) offering
simultaneous wavelength coverage of $4184-5580$
$\mathrm{\SIUnitSymbolAngstrom}$ and $5294-7060$
$\mathrm{\SIUnitSymbolAngstrom}$, respectively. Ne-Ar arc lamp exposures were
taken immediately following each science observation for wavelength
calibration. Data were reduced using the
PyWiFeS333http://pywifes.github.io/pipeline/ pipeline (Childress et al.,
2014). Data were taken between 2022 June 23 and 2023 April 19 for a total of
24 observations. Individual observations were taken with two 700 $\mathrm{s}$
exposures (separated temporally only by readout time) so that the duration of
each exposure was less than 1% of the orbit of the binary
($P_{\mathrm{orb}}=0.8952134$ d; Soszyński et al. 2016). Resultant spectra
typically have a signal-to-noise of $>250$ per resolution element in the
continuum.

We obtained a single 300 s spectrum of the source with the Goodman
Spectrograph (Clemens et al., 2004) on the SOAR telescope on 2022 October 30.
This spectrum used the 2100 l mm-1 grating and a 1.2″ longslit, giving a full
width at half maximum resolution of about 0.9 Å over the wavelength range
6090–6660 Å. The spectrum was optimally extracted and wavelength calibrated
with a time-adjacent arc lamp exposure using standard tools in IRAF (Tody,
1986), with a signal-to-noise of $\sim 85$ per resolution element in the
continuum.

##  3 Analysis and Results

###  3.1 Radial Velocity Estimation

The output of the PyWiFeS pipeline is a calibrated data-cube with a spectrum
(frequency and brightness) for each pixel in the IFU. We extracted spectra of
Gaia DR3 4042390512917208960 from pixels of the WiFeS IFU centred on the
source whose total flux was greater than 5 times the root mean square (RMS)
flux of background pixels, followed by a subtraction of an averaged background
spectrum. To improve signal, for each observation, the two exposures were
averaged before proceeding to radial velocity (RV) estimation, leaving us with
24 RVs spread over 10 months.

RVs were estimated via cross-correlation using the IRAF task fxcor (Tody,
1986, 1993) and a synthetic template spectrum from the PHOENIX
library444http://phoenix.astro.physik.uni-goettingen.de (Husser et al., 2013),
which was convolved down to a comparable spectral resolution to the WiFeS
data. The library used was the PHOENIX medium resolution grid with uniform
grid spacing of 1 $\mathrm{\SIUnitSymbolAngstrom}$. As this spectrum is mainly
for determining RVs, and not detailed modelling of the star, we only briefly
experimented with metallicity and alpha element abundances, finding no obvious
benefit from the use of non-zero values. The effective temperature
($T_{\mathrm{eff}}$) and surface gravity ($\log g$) of Gaia DR3
4042390512917208960 reported in DR3 are $6524^{+31}_{-16}$ K and
$3.804^{+0.002}_{-0.003}$, respectively (Creevey et al., 2023; Fouesneau et
al., 2023), corresponding to a spectral type of F4V or F4IV. We matched
PHOENIX spectra to the WifeS spectra of the target by eye, finding good
agreement between the PHOENIX spectrum with $T_{\mathrm{eff}}=6500$
$\mathrm{K}$ and $\log g=4.0$ and the observed spectrum. We found the PHOENIX
spectrum with $T_{\mathrm{eff}}=6700$ $\mathrm{K}$ and $\log g=4.0$ to be the
best match when considering the depths of the absorption lines. The match
between the $6700$ $\mathrm{K}$ template and observation line depth was almost
identical for H$\beta$ and H$\gamma$. However, the observed H$\alpha$ line
depth is significantly shallower. Comparisons of the template and observed
spectra are shown in Figure 3.

We used the ‘blue’ WiFeS spectra for RV estimation due to the abundance of
metal lines. Both H$\beta$ and H$\gamma$ were excluded when estimating the RVs
to minimise possible contamination arising from emission processes. We used
the following regions to estimate RVs with fxcor: 4200-4300 Å , 4385-4820 Å ,
4900-5400 Å (Figure 4). All velocities were corrected to the solar system
barycenter. A systematic uncertainty of 3
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ was added to the uncertainty of each
estimated RV (Kuruwita et al., 2018). The time of a measured RV was shifted
from local observation time to Barycentric Modified Julian Date (BMJD) in the
Barycentric Dynamical Time standard (TDB; Eastman et al. 2010). The RV from
the SOAR spectrum was estimated in the same way as the WiFeS RVs using the
same PHOENIX spectrum and fxcor. H$\alpha$ was excluded from the region that
was cross-correlated to minimise contamination from emission processes.
Furthermore, a correction of $1$ $\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ was
added based on the telluric oxygen lines in the observation. The estimated RVs
for Gaia DR3 4042390512917208960 are presented in Table 1, and the RV from the
SOAR spectrum is identified by an asterisk appended to the BMJD. Reported
times in Table 1 are the mid-observation time.

Figure 3: Comparison of PHOENIX template ($T_{\mathrm{eff}}=6700$ $\mathrm{K}$, $\log g=4.0$) spectrum with the shifted and averaged spectrum of all observations of Gaia DR3 4042390512917208960. The left, middle, and right panels are cutouts around H$\gamma$, H$\beta$, and H$\alpha$, respectively. Broadly there is good agreement between the synthetic PHOENIX spectrum and the averaged observed spectrum. Note that the depth of the H$\alpha$ line is much shallower in the observed spectrum. Figure 4: Comparison of PHOENIX template ($T_{\mathrm{eff}}=6700$ $\mathrm{K}$, $\log g=4.0$) spectrum with the shifted and averaged spectrum of all observations of Gaia DR3 4042390512917208960 in the regions used with fxcor. The shaded regions indicate sections of the spectra that were excluded from the cross-correlation. Note that while we exclude H$\gamma$ and H$\beta$, the region includes prominent lines from many atomic species (e.g., Fe, Ca, Na, Mg). BMJD | RV ($\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$)
---|---
59781.6833 | $19.3\pm 9.0$
59810.6610 | $-41\pm 10$
59753.4166 | $-46.4\pm 9.2$
59753.4775 | $-29.6\pm 7.2$
59753.5465 | $-12.5\pm 6.5$
59753.5947 | $-6.9\pm 6.5$
59753.6610 | $1.2\pm 6.6$
59753.6969 | $6.7\pm 7.1$
59781.4678 | $16.0\pm 7.4$
59781.5217 | $34.6\pm 8.2$
59781.5645 | $29.9\pm 8.2$
59781.6163 | $28.8\pm 8.6$
59810.4320 | $4.0\pm 6.4$
59810.4941 | $-7.3\pm 6.4$
59810.5569 | $-30.0\pm 7.7$
59822.4157 | $-32.8\pm 7.5$
59822.4776 | $-11.9\pm 8.7$
59822.5383 | $-4.3\pm 6.2$
59840.4369 | $-6.0\pm 5.9$
59840.4994 | $-4.6\pm 7.8$
59868.4355 | $34.4\pm 9.6$
59882.0102* | $-0.2\pm 9.0$
60032.7391 | $-43.4\pm 8.2$
60032.7580 | $-43.7\pm 7.9$
60053.7026 | $32.8\pm 9.5$
Table 1: RVs for Gaia DR3 4042390512917208960 measured using cross-
correlation. An asterisk appended to the BMJD indicates the RV was estimated
from the SOAR spectrum. The errors in the WiFeS RVs include the 3
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ systematic.

###  3.2 Binary Orbital Parameters

To fit Keplerian orbits to the estimated RVs we used The Joker (Price-Whelan
et al., 2017); a custom Monte Carlo (MC) sampler designed for this problem.
The parameters we consider are the binary orbital period ($P$), the orbital
eccentricity ($e$), the radial velocity semiamplitude of the luminous star
($K$), the systemic radial velocity of the binary ($\gamma$), a (constant)
jitter to account for underestimation of RV errors ($s$), the argument of
periastron ($\omega$), and the orbital phase at the reference time ($M_{0}$;
by default the reference time used in The Joker is the time of the earliest
observation). It is important to note that there may be a systematic offset
between the WiFeS and SOAR RVs, but with only 1 SOAR RV we cannot quantify it,
and thus we do not attempt to fit for the offset. The Joker samples from the
prior probability density functions (PDFs) for each parameter we are
considering and then performs rejection sampling, producing samples of the
posterior for each parameter.

The priors chosen have significant impact on the rejection sampling when the
data is sparse or uninformative and should be considered carefully. The
observations of Gaia DR3 4042390512917208960 are quite dense and thus the
choice of prior is unlikely to significantly affect the results. Nevertheless,
we trialled different priors. We used the default prior implemented in the The
Joker, passing values for the the standard deviation of the prior on systemic
velocity, $v_{0}$, and radial velocity semiamplitude, $K_{0}$, based on the
RVs. We trialled different values to asses the impact of these priors on the
final results. We used The Joker’s default log-normal prior on period, to
which we supplied $P_{min}=0.25$ d and $P_{max}=16384$ d as a blind search. We
also trialled a much narrower prior centred around the orbital period
determined from photometry.

When few samples survive and the posterior appears to be unimodal, the
posterior samples from the rejection sampling should be used in standard
Markov Chain Monte Carlo (MCMC). We used the No U-Turn Sampling (NUTS; Hoffman
& Gelman 2014) MCMC method to sample the posterior, as implemented in
pymc3_ext. pymc3_ext contains pymc3 (Salvatier et al., 2016) extensions
extracted from exoplanet (Foreman-Mackey et al., 2021). The priors used here
are the same as above, and the results of the rejection sampling are used to
initialise the MCMC. Convergence was verified using the Gelman-Rubin
diagnostic test (Gelman & Rubin, 1992), ensuring the Gelman-Rubin statistic
($\hat{R}$) was close to 1 for all best-fit parameters. We ran 16 chains with
2,000 tune iterations and 8,000 draws, producing 128,000 final draws.

Changing the standard deviation of the prior on $v_{0}$ and $K$ had no impact
on inferred values, and neither did the choice of period prior. The periods
independently inferred from OGLE photometry and WiFeS radial velocities are
consistent within error, supporting our belief that the spectra are not
affected by significant contamination from nearby stars in the Galactic Bulge.
Therefore, we only present one set of results. The binary orbital parameters
inferred from The Joker analysis of the RVs are listed in Table 2. While there
is some evidence in the posterior for non-zero eccentricity, $e=0$ cannot be
excluded. The 95% upper limit on the eccentricity is 0.18. A comparison
between the phase-folded RVs and the best fitting model, with residuals, is
shown in Figure 5.

There is a noticeable scatter in the RVs visible in Figure 5. Additionally,
the uncertainty in the RVs determined using fxcor is larger than we would have
expected, based on spectra taken with this instrument of other targets. A
plausible explanation for both of these is that the radial velocity from a
luminous secondary is causing the scatter. However, excluding spectral
disentangling (Section 3.4), we found no evidence for the presence of the
secondary in the spectra. It is believed that the WiFeS instrument has a
degree of inherent instability (see Kuruwita et al. 2018; standard deviation
of the instability is $\sim 3$ $\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$), this
could also be contributing to the scatter we see. However, the agreement in
orbital period separately determined from photometry and The Joker analysis
gave us confidence in there being no large systematic issue in the RV
modelling.

We then calculated the spectroscopic mass function of the secondary,
$\mathbf{m_{f}}$, using Equation 1

$\mathbf{m_{f}=}\frac{M^{3}_{2}\sin^{3}i}{(M_{1}+M_{2})^{2}}=\frac{PK_{1}^{3}}{2\pi
G}(1-e^{2})^{3/2}$ (1)

where $G$ is the Newtonian constant of gravitation, $i$ is the inclination,
$M_{1}$ is the mass of the luminous star, $M_{2}$ is the mass of the unseen
object, and the remaining parameters are as described in Table 2. We used the
posteriors from the binary parameter estimation to calculate the mass function
of the secondary star $\mathbf{m_{f}=}3.87^{+0.88}_{-0.77}\times 10^{-3}$ M⊙.

Parameter | Value
---|---
$P$ (d) | $0.89526^{+0.00013}_{-0.00012}$
$e$ | $0.043^{+0.057}_{-0.033}$
$K$ ($\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$) | $34.8^{+2.6}_{-2.5}$
$\gamma$ ($\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$) | $-5.3^{+1.7}_{-1.7}$
$s$ ($\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$) | $0.0088^{+0.4905}_{-0.0088}$
$\omega$ (rad) | $2.6^{+2.0}_{-1.4}$
$M_{0}$ (rad) | $-0.7^{+2.0}_{-1.4}$
Table 2: Orbital parameters determined through modelling the estimated RVs
with a Keplerian orbit using The Joker. Figure 5: Top panel: phase-folded RVs
overplotted with the best-fit Keplerian model. Bottom panel: residuals of
model minus data as a function of orbital phase. The blue data point is the RV
from the SOAR spectrum. The WiFeS RVs are in black.

###  3.3 Joint Radial Velocity and Photometry Modelling

Now that we have RVs for Gaia DR3 4042390512917208960 we can model them in
conjunction with the excellent OGLE photometry. We started jointly modelling
the photometry and RVs using the Eclipsing Light Curve (ELC) code described in
Orosz & Hauschildt (2000). Initial results from ELC suggested the binary was
composed of two stars, and that both the stars may be filling their Roche
lobes, prompting us to switch to PHOEBE (Prša & Zwitter, 2005; Conroy et al.,
2020). We trialled different binary configurations using PHOEBE, including
detached binaries, semi-detached binaries, binaries on the contact binary
boundary (both stars are almost filling their Roche lobe), and contact
binaries. We used the period and reference time inferred from the OGLE-IV
photometry due to the long time baseline of observations. We calculated the
flux from magnitudes using the reference wavelength ($\lambda_{\mathrm{eff}}$)
for the OGLE-IV I-band filter taken from the Spanish Virtual Observatory (SVO;
Rodrigo et al. 2012; Rodrigo & Solano 2020)555http://svo2.cab.inta-
csic.es/svo/theory/fps3/index.php assuming OGLE magnitudes are in the AB
magnitude system. We selected the logarithmic limb-darkening coefficients of
(0.5761,0.2210) from Claret (2017) using their radial least squares fit to the
log law model for $T_{\mathrm{eff}}=6500$ $\mathrm{K}$ and $\log g=4.0$.

PHOEBE modelling of this system supports the absence of a BH. The large
amplitude of the photometric modulation and the low radial velocity
semiamplitude of the luminous star are incompatible with the large mass of the
unseen object inferred from the photometry using the mMMR. They are also
incompatible with a compact object in general. Modelling the system with a
neutron star companion (maintaining orbital and stellar parameters constrained
from observations) produces a lightcurve with a much lower level of
photometric modulation, very different lightcurve minima, and a much larger RV
semiamplitude. These cannot be resolved through inclination effects. A larger
amplitude of photometric modulation requires a higher inclination, however,
decreasing the RV semiamplitude and difference in photometric minima requires
a lower inclination. Furthermore, reconciling the observed similarity in
photometric minima depths with the high level of photometric modulation is
challenging with only a single luminous star due to the physical distortion of
the star that is required to recover the large photometric modulation.

When also considering the shape of the lightcurve, both stars must have a high
Roche lobe filling factor. This places the system on the boundary between a
contact binary and a detached binary with high filling factors. We are not
able to discern which is the correct interpretation for this system. The
similar depth minima in the lightcurve suggest that each star has a similar
temperature. In a contact binary each star typically has a similar effective
temperature and in turn likely has a similar spectrum (e.g., Pribulla et al.
2003; Mitnyan et al. 2020). Assuming a contact binary configuration, we model
the system as a contact binary in PHOEBE to constrain system parameters. An
inclination and mass ratio of around $60$$\mathrm{\SIUnitSymbolDegree}$ and
0.2, respectively, seem likely. However, in this configuration, there is an
apparent contradiction with previous conclusions; the secondary is
contributing a significant amount of flux. This would suggest the observed
spectra should be double lined and not single lined as they appear.

###  3.4 Spectral Disentangling

The joint modelling of the photometry and RVs with PHOEBE led us to infer the
spectrum of the unseen companion was hidden in the observed spectrum. To test
this we explored spectral disentangling, using the shift-and-add approach
(e.g., González & Levato 2006; Shenar et al. 2020, 2022). When applied to
binaries this technique assumes the observed spectra are composites of two
individual spectra and exploits that spectral features belonging to each star
will move in anti-phase, arising from Doppler shifts. If $K_{1}$ and/or
$K_{2}$ are unknown and the other orbital parameters are known, a grid-based
exploration can be carried out to infer the $K_{1}$, $K_{2}$ values that best
reproduce the data. This involves evaluating $\chi^{2}(K_{1},K_{2})$, the
formalism of which is presented in Shenar et al. (2022).

We set the orbital parameters using the period and reference time from OGLE
and the systemic velocity from The Joker. Eccentricity was set to 0. This
choice is motivated by the indications that the system is a contact binary or
on the boundary, and that the RVs determined using the method in Section 3.1
are unreliable due to possible smearing (Sections 3.3 and later in this
Section). We trialled two different explorations of $K_{1}$ and $K_{2}$. The
first was holding $K_{1}$ constant as determined with The Joker and varying
$K_{2}$. The second was varying both $K_{1}$ and $K_{2}$ to investigate the
possibility of smearing impacting the RVs inferred using the method described
in Section 3.1. The initial guess for this second method was $K_{1,i}=50$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, $K_{2,i}=160$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, and the grid was split into 100
linearly spaced components between $0.1K_{i}$ and $2K_{i}$. We disentangled
the regions around H$\alpha$, H$\beta$, H$\gamma$, $4270-4330$
$\mathrm{\SIUnitSymbolAngstrom}$, and $5150-5200$
$\mathrm{\SIUnitSymbolAngstrom}$ as these are the deepest absorption features
in the spectrum. We disentangled these regions independently as well as
jointly, and also used the entire spectrum as the disentangling region.

Spectral disentangling reveals this system is likely a double-lined
spectroscopic binary (SB2), and not a single-lined spectroscopic binary (SB1)
as initially thought. Spectral disentangling while holding $K_{1}$ constant
suggested the presence of a luminous companion. However, the features in the
disentangled spectrum did not correspond to a physically believable stellar
spectrum. Therefore, we only present the results of the grid disentangling
where both $K_{1}$ and $K_{2}$ were varied. Disentangling the regions as
described above return consistent values of $K_{1}$ and $K_{2}$ in all tests,
and the resulting disentangled spectrum is consistent with a stellar spectrum.
Disentangling all regions excluding H$\alpha$ (it is not in the ‘blue’ WiFeS
spectra) produced $K_{1}=54.7\pm 7.2$ $\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$
and $K_{2}=177\pm 21$ $\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. The $\chi^{2}$
map is shown in Figure 6. The value for $K_{1}$ inferred from spectral
disentangling is discrepant at the $2\sigma$ level with the $K_{1}$ inferred
assuming the binary is single-lined. The disentangled spectrum for each
component is presented in Figure 7 along with the PHOENIX template spectrum.
They both appear to have a spectral type of either F4V or F4IV, sharing the
same prominent lines and species as the spectra discussed in Section 3.1.

Figure 6: $\chi^{2}$ surface explored in the $K1$, $K2$ grid search, with the
green dashed contour indicating 1$\sigma$. The semiamplitudes inferred from
this exploration were $K_{1}=54.7\pm 7.2$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ and $K_{2}=177\pm 21$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. The $K_{1}$ inferred assuming the
binary was single lined ($34.8^{+2.6}_{-2.5}$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$) is discrepant with the value inferred
using this technique at the $2\sigma$ level. Figure 7: Comparison between the
disentangled spectra of the components and the PHOENIX template model
identified as the best match to the observed data. The secondary has been
scaled for a light ratio such that $l_{2}/l_{\mathrm{tot}}=0.2$. Noticeably
the disentangled spectra are very similar. The lower left and right panels are
cutouts around H$\gamma$ and H$\beta$, respectively. In all panels the
normalised template has been shifted vertically by 0.3 for clarity.

##  4 Discussion

The results inferred from our dedicated spectroscopic campaign and existing
photometry imply there is no BH in this system. Considering the RVs in
isolation, the system would need to be almost face-on to accommodate a BH.
Jointly modelling both RVs and photometry demonstrate that there is no BH in
the system, and that the system is more in line with a contact binary or a
stellar binary where both stars are almost filling their Roche lobes.

###  4.1 Radial Velocity Analysis

Based on the spectroscopic mass function alone, the presence of a BH seems
unlikely. If we assume the mass of the luminous star in the binary to be
$M=1.7\pm 0.05$ M⊙ – the midpoint between the mass estimated by Gaia single
star evolutionary models (Creevey et al., 2023; Fouesneau et al., 2023) and
the mass from the PHOENIX template spectrum – the minimum mass of the
secondary is $0.25^{+0.02}_{-0.02}$ M⊙. This implies that for the secondary to
have a mass greater than 3 M⊙ the inclination must be less than
$5$$\mathrm{\SIUnitSymbolDegree}$. However, the large photometric variation
rules out low inclination angles (Section 3.3). For a conservative lower limit
on inclination from the photometry ($i>50$$\mathrm{\SIUnitSymbolDegree}$), the
mass of the primary would have to be greater than 53 M⊙ for the mass of the
unseen secondary to be $>3$ M⊙. The presence of a BH is therefore ruled out.

However, it is important to consider the results of spectral disentangling
discussed in Section 4.2. The RVs inferred as described in Section 3.1 are
likely smeared due to the presence of a luminous companion in the spectra.
This implies that $K$ and $e$ as inferred through this analysis are likely not
reliable. In particular, the eccentricity inferred here should not be used as
the basis for an argument against the contact binary scenario. If the
eccentricity is truly non-zero this argues strongly for a detached system with
both stars nearly, but not quite, filling their Roche lobes.

###  4.2 The Single- or Double-lined Nature of the Spectra

Based on the configuration determined using PHOEBE, the secondary is
contributing a significant fraction of the light of the system ($>10$%).
However, when visually inspecting the observed spectra we see no evidence for
the lines of the secondary. These two conclusions seem mutually exclusive
based on the quality of our spectra.

The results of spectral disentangling suggest the system is an SB2 binary,
with the spectrum of the secondary buried in the observed spectra. However,
the spectra of each component are very similar (Figure 7). Small errors in
$K_{1}$ can lead to the secondary’s disentangled spectrum spuriously mimicking
that of the primary (Shenar et al., 2022). However, in our case the $K_{1}$
value determined from spectral disentangling is $\sim 60$% larger than the
$K_{1}$ value determined assuming the spectra are single lined. Furthermore,
the disentangled spectrum and the $K_{1}$ and $K_{2}$ estimates are consistent
regardless of the region of the spectra that is used (as stated in Section
3.4). Assuming the secondary contributes $15-30$% of the light (motivated by
PHOEBE modelling), we simulated synthetic SB2 spectra with similar resolution,
sampling, and signal to noise as the observed spectra using the PHOENIX
template as the input spectra for both components. Exploring light ratios from
0.0001 to 0.5 supports our results (Appendix A). Additionally, the magnitude
of $K_{1}$ identified from The Joker is less than the $K_{1}$ identified from
spectral disentangling. This is typical behaviour when the spectra of both
components are blended together (Bodensteiner et al., 2021; Banyard et al.,
2022; Shenar et al., 2022). As a result, we are inclined to trust the results
of the spectral disentangling analysis. However, we stress that further
observations with a higher resolution spectrograph and more stable instrument
are required to confirm the detection of the secondary.

That the disentangled spectra are similar is not surprising if the binary is a
contact binary, and in turn provides evidence that the system is a contact
binary. While this is a circular argument, it is self-consistent, and in line
with all other avenues of investigation for this system. For a detached binary
with such a mass ratio it would seem unlikely that both components have such a
similar spectrum.

With fitted values for $K_{1}$ and $K_{2}$ we can determine the mass ratio of
the binary. Assuming a circular Keplerian orbit,
$q=\frac{K_{1}}{K_{2}}=0.309^{+0.060}_{-0.051}$. Assuming $M_{1}=1.7\pm 0.05$
M⊙ (as discussed in Section 4.1), then $M_{2}=0.525^{+0.103}_{-0.087}$ M⊙.

###  4.3 Discrepant mass ratios

The mass ratio determined from spectral disentangling
($0.309^{+0.060}_{-0.051}$) is significantly different than the mMMR of 2.54
determined by assuming the system is an ellipsoidal variable. The mMMR derived
by Gomel et al. (2021a) assumes that all the light in the ellipsoidal system
is coming from one star, and thus the observed photometric modulation arises
due to the distorted star. In the scenario described, and under the
assumptions that the star fills its Roche lobe and the binary is observed
edge-on, the mMMR is always less than the actual mass ratio. The discrepancy
between the mMMR and the mass ratio we determine for Gaia DR3
4042390512917208960 can be explained the same way regardless of whether the
system is a contact binary or a binary with both stars being close to filling
their Roche lobes. The assumption that all the observed light is coming from a
single star breaks down for stellar binaries. In stellar binaries the change
in projected surface area as a function of orbital phase is significantly
larger than in a star-BH binary due to the presence of two stars. As a result,
large modulations can be produced in the light curve purely from the change in
projected surface area. As the Gomel et al. (2021a) method assumes the light
is coming from one tidally distorted star, and the mMMR depends primarily on
the amplitude of modulation, binaries with two stars will have inflated mMMRs.

Furthermore, this biases samples selected based on high mMMR towards stellar
binaries. A random sample of contact binaries will, on average, have higher
mMMRs than a random sample of star-BH binaries with the same broad system
parameters. If an initial sample of ellipsoidals are contaminated with stellar
binaries, then a sample of high mMMR systems will preferentially include these
stellar binaries. The problem is therefore one of classification. While the
Gomel et al. (2021a) method appears robust, it is only so if it is applied to
true ellipsoidal variables. Reducing the contamination should be the priority
of future studies applying the Gomel et al. (2021a) methodology to identify
BHs.

###  4.4 Classification of Gaia DR3 4042390512917208960

The photometry, spectral disentangling, and PHOEBE modelling we have carried
out indicate that Gaia DR3 4042390512917208960 is likely a contact binary, or
on the boundary with both stars nearly filling their Roche lobes.
Classification of the system is not straightforward from photometry alone. It
was incorrectly classified as an ellipsoidal variable by Gaia and OGLE, and
the initial ASAS-SN variable classification did not have ellipsoidal as a
category. The ellipsoidal catalogue from ASAS-SN produced later is not public,
preventing examination of their classification (Rowan et al., 2021). Gaia
(Eyer et al., 2023) and ASAS-SN (Jayasinghe et al., 2020) used machine
learning classifiers, whereas OGLE (Soszyński et al., 2016) primarily used
template fitting. The templates used by the OGLE team were constructed using
the OGLE-IV photometry of bright, previously classified ellipsoidal and
eclipsing systems in the best-sampled fields. The solution for better
discrimination between contact binaries and ellipsoidal variables from
photometry alone is also not straightforward. Whilst mis-classification is
often assumed to be due to sparse data, in the case of Gaia DR3
4042390512917208960, simply having more photometric data is not sufficient to
resolve the issue. The OGLE-IV I-band light curve for Gaia DR3
4042390512917208960 has 649 data points spanning a time baseline of more than
2000 d (see Figures 1 and 2), and they also mis-classify the system. However,
Gaia and OGLE use different methods of classification. Clearly, the
classification problem is not limited to one technique, or easily corrected by
having more data.

###  4.5 Contamination of Ellipsoidal Variable Samples

The source studied in this work is either a contact binary or on the boundary
with both stars nearly filling their Roche lobes. The 14 systems studied by
Nagarajan et al. (2023) are likely contact binaries, and unlikely to host BHs.
It seems likely that the contamination of the Gomel et al. (2023) sample with
non-ellipsoidal systems is very high. As discussed in Section 4.3, when there
is contamination, stellar binaries are likely to dominate in samples selected
for high mMMR.

Gomel et al. (2023) discussed the difficulty in discriminating between contact
binaries and true ellipsoidal systems. They attempted to minimise the
contamination by (a) limiting the systems to those with orbital periods
$>0.25$ d, and (b) requiring the light curve to have unequal minima. The
periods explored by Nagarajan et al. (2023) extend to 0.75 d, and the system
in this work has a period of 0.8952 d. This suggests that the $P>0.25$ d
restriction is not conservative enough to exclude contact binaries.
Furthermore, while unequal minima may help discriminate against other types of
variables, it does not exclusively exclude contact binaries and include
ellipsoidal variables. Examples of contact binaries with unequal minima
include the system studied in this work, systems in Nagarajan et al. (2023) if
they are truly contact binaries, and also classical contact binary systems
such as 44 Boötis, VW Cephei, and Y Sextantis, which can be seen in Transiting
Exoplanet Survey Satellite (TESS; Ricker et al. 2015) photometry.

Whilst the mMMR is an excellent technique for identifying candidate compact
objects when applied to ellipsoidal variables, the contamination by stellar
binaries in the Gomel et al. (2023) sample hampers the viability of using the
technique in practice. Further studies should prioritise minimising
contamination by stellar binaries as much as possible before turning to RV
studies. Again, we re-iterate that removing the contaminating binaries is not
a simple problem, as discussed in Section 4.4. Furthermore, while for these
studies reliable classification is critical, it is also important that the
classification is not computationally expensive. For example, Eyer et al.
(2023) analysed the photometric and spectroscopic time series of 1.8 billion
sources from Gaia DR3, identifying 10.5 million sources as variable. In the
future, the number of known variables will only grow, and so too will the
volume of data available (e.g., future Gaia data releases, data from new
missions such as the Vera C. Rubin Observatory (Ivezić et al., 2019) and the
Nancy Grace Roman Space Telescope (Spergel et al., 2015)).

##  5 Conclusions

We have presented a follow-up investigation of one of the most promising BH
candidates from the Gomel et al. (2023) sample of compact object candidates.
This sample was identified by leveraging the relationship between optical
light curve modulation amplitude and mass ratio for ellipsoidal variables
(Gomel et al., 2021a) and applying it to sources from Gaia DR3. The system we
identified as being one of the most interesting candidates from this sample,
Gaia DR3 4042390512917208960, had a mMMR of 2.5, implying the unseen secondary
was significantly more massive than the luminous star. With detailed multi-
epoch spectral observations, we characterised the system, concluding that the
binary is unlikely to host a BH. We find the system is most likely a contact
binary, with a mass ratio of $0.309^{+0.060}_{-0.051}$, an order of magnitude
less than the mMMR. We suggest that this discrepancy arises from the
assumption underlying the Gomel et al. (2021a) method that all the light comes
from one star and is due to ellipsoidal modulation. In contact binaries there
are two stars, which gives rise to a significantly larger change in projected
surface area over the binary orbit than for ellipsoidal variables with similar
binary system parameters. This causes a larger modulation in the optical light
curve and results in an artificially large mMMR. This highlights the main
issue with leveraging the Gomel et al. (2021a) technique to find BHs: without
better filtering to remove contact binaries (or other stellar binaries),
samples selected based on high mMMR are likely to be heavily contaminated.
However, while the contact binary interpretation requires a future study to
confirm the detection of the secondary’s lines, the rationale for the
contamination of the Gomel et al. (2021a) sample is still valid for a detached
binary in which both stars are close to filling their Roche lobe. Finally,
this work highlights the need for dynamical studies to both confirm new BHs
and evaluate exciting new techniques for identifying BH candidates.

The authors thank Adela Kawka, Chris Lidman, Fiona H. Panther, Ian Price,
Kathryn Ross, Katie Auchettl, Michael S. Bessel, and Michael Ireland for
assistance with WiFeS, Michael Abdul-Masih, Kyle E. Conroy, Sara Saracino, and
Sebastian Kamann for assitance with PHOEBE, Robin Humble for assistance with
the OzSTAR Supercomputer, and Tsevi Mazeh and Thomas J. Maccarone for helpful
discussions. The authors thank the reviewer for their constructive comments
that helped improve this work. T.N.O’D was supported by a Forrest Research
Foundation Scholarship, and an Australian Government Research Training Program
(RTP) Stipend and RTP Fee-Offset Scholarship. This work was supported by the
Australian government through the Australian Research Council’s Discovery
Projects funding scheme (DP200102471). We acknowledge extensive use of the
SIMBAD database (Wenger et al., 2000), NASA’s Astrophysics Data System, and
arXiv. 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. We acknowledge the
traditional owners of the land on which the ANU 2.3 m telescope stands, the
Gamilaraay people, and pay our respects to elders, past and present. Based in
part on observations obtained at the Southern Astrophysical Research (SOAR)
telescope, which is a joint project of the Ministério da Ciência, Tecnologia e
Inovações (MCTI/LNA) do Brasil, the US National Science Foundation’s NOIRLab,
the University of North Carolina at Chapel Hill (UNC), and Michigan State
University (MSU). Part of this work was performed on the OzSTAR national
facility at Swinburne University of Technology. The OzSTAR program receives
funding in part from the Astronomy National Collaborative Research
Infrastructure Strategy (NCRIS) allocation provided by the Australian
Government, and from the Victorian Higher Education State Investment Fund
(VHESIF) provided by the Victorian Government. This research has made use of
the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded
by MCIN/AEI/10.13039/501100011033/ through grant PID2020-112949GB-I00.

##  Appendix A Spectral Disentangling

Disentangling was performed as described in Section 3.4. The SB2 synthetic
data were generated using the PHOENIX template spectrum as the model for each
star, reflecting the similarity between the disentangled spectra for Gaia DR3
4042390512917208960. 24 synthetic observations were created to reflect the
number of observations of the source, and the model spectra were convolved to
a resolution and sampling factor similar to the WiFeS data. The signal-to-
noise of the synthetic spectra were set so such that the synthetic data were
similar to the WiFeS spectra. The same orbital period as Gaia DR3
4042390512917208960 was used, an eccentricity of 0, and $K_{1}=50$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ and $K_{2}=160$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$. The six light contribution fractions
($l_{2}/l_{\mathrm{tot}}$) examined were 0.5, 0.3, 0.2, 0.1, 0.05, and 0.0001
(reflective of no luminous companion). The $K_{1}$, $K_{2}$ grid that was
explored was quite coarse; split into 15 linearly spaced components between
$0.1K_{i}$ and $2K_{i}$. For this test, $K_{1,i}=60$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, $K_{2,i}=192$
$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ to shift the centre of the explored
grid away from the correct values.

The results presented in Figure 8 suggest we could reliably constrain $K_{2}$
and the general shape of the secondary’s spectrum when it is contributing
$>20$% of the total light when both components have the same spectrum. At
light contribution fractions of 0.1, 0.05, and 0.0001, the disentangled
$K_{2}$ is spurious and the disentangled spectrum of the secondary can mimic
the primary. The $\chi^{2}$ surface for light contribution fractions of 0.1,
0.05, and 0.0001 is not smooth and has no clear minima, whereas for light
contribution fractions of 0.2, 0.3, and 0.5 the surface is smooth with a clear
minimum. The smooth $\chi^{2}$ map for Gaia DR3 4042390512917208960 (Figure 6)
resembles an accurate disentanglement far more than any of the spurious
results described in this section. While not conclusive in isolation, we
believe it is another piece of supporting evidence that our disentangling
results are accurate.

Figure 8: Comparison of disentangled spectra with the input spectrum around
H$\gamma$ and H$\beta$ for six different light contribution fractions. All
secondary spectra are scaled by their light contribution except the bottom row
which is scaled assuming $l_{2}/l_{\mathrm{tot}}=0.05$. These results suggest
that we could reliably detect the presence of a luminous companion if it was
contributing at least 20% of the total light of the system.

## References

  * Andrew et al. (2022) Andrew, S., Penoyre, Z., Belokurov, V., Evans, N. W., & Oh, S. 2022, MNRAS, 516, 3661, doi: 10.1093/mnras/stac2532
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Banyard et al. (2022) Banyard, G., Sana, H., Mahy, L., et al. 2022, A&A, 658, A69, doi: 10.1051/0004-6361/202141037
  * Boch & Fernique (2014) Boch, T., & Fernique, P. 2014, in Astronomical Society of the Pacific Conference Series, Vol. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay, 277
  * Bodensteiner et al. (2021) Bodensteiner, J., Sana, H., Wang, C., et al. 2021, A&A, 652, A70, doi: 10.1051/0004-6361/202140507
  * Bonnarel et al. (2000) Bonnarel, F., Fernique, P., Bienaymé, O., et al. 2000, A&AS, 143, 33, doi: 10.1051/aas:2000331
  * Brown & Bethe (1994) Brown, G. E., & Bethe, H. A. 1994, ApJ, 423, 659, doi: 10.1086/173844
  * Childress et al. (2014) Childress, M. J., Vogt, F. P. A., Nielsen, J., & Sharp, R. G. 2014, Ap&SS, 349, 617, doi: 10.1007/s10509-013-1682-0
  * Claret (2017) Claret, A. 2017, A&A, 600, A30, doi: 10.1051/0004-6361/201629705
  * Clemens et al. (2004) Clemens, J. C., Crain, J. A., & Anderson, R. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5492, Ground-based Instrumentation for Astronomy, ed. A. F. M. Moorwood & M. Iye, 331–340, doi: 10.1117/12.550069
  * Conroy et al. (2020) Conroy, K. E., Kochoska, A., Hey, D., et al. 2020, ApJS, 250, 34, doi: 10.3847/1538-4365/abb4e2
  * Corral-Santana et al. (2016) Corral-Santana, J. M., Casares, J., Muñoz-Darias, T., et al. 2016, A&A, 587, A61, doi: 10.1051/0004-6361/201527130
  * Creevey et al. (2023) Creevey, O. L., Sordo, R., Pailler, F., et al. 2023, A&A, 674, A26, doi: 10.1051/0004-6361/202243688
  * Dopita et al. (2007) Dopita, M., Hart, J., McGregor, P., et al. 2007, Ap&SS, 310, 255, doi: 10.1007/s10509-007-9510-z
  * Dopita et al. (2010) Dopita, M., Rhee, J., Farage, C., et al. 2010, Ap&SS, 327, 245, doi: 10.1007/s10509-010-0335-9
  * Eastman et al. (2010) Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935, doi: 10.1086/655938
  * El-Badry & Burdge (2022) El-Badry, K., & Burdge, K. B. 2022, MNRAS, 511, 24, doi: 10.1093/mnrasl/slab135
  * El-Badry et al. (2022) El-Badry, K., Seeburger, R., Jayasinghe, T., et al. 2022, MNRAS, 512, 5620, doi: 10.1093/mnras/stac815
  * El-Badry et al. (2023a) El-Badry, K., Rix, H.-W., Quataert, E., et al. 2023a, MNRAS, 518, 1057, doi: 10.1093/mnras/stac3140
  * El-Badry et al. (2023b) El-Badry, K., Rix, H.-W., Cendes, Y., et al. 2023b, arXiv e-prints, arXiv:2302.07880, doi: 10.48550/arXiv.2302.07880
  * Eyer et al. (2023) Eyer, L., Audard, M., Holl, B., et al. 2023, A&A, 674, A13, doi: 10.1051/0004-6361/202244242
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Foreman-Mackey et al. (2021) Foreman-Mackey, D., Luger, R., Agol, E., et al. 2021, exoplanet: Gradient-based probabilistic inference for exoplanet data & other astronomical time series, 0.5.1, doi: 10.5281/zenodo.1998447
  * Fouesneau et al. (2023) Fouesneau, M., Frémat, Y., Andrae, R., et al. 2023, A&A, 674, A28, doi: 10.1051/0004-6361/202243919
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211, doi: 10.48550/arXiv.2208.00211
  * Gavras et al. (2023) Gavras, P., Rimoldini, L., Nienartowicz, K., et al. 2023, A&A, 674, A22, doi: 10.1051/0004-6361/202244367
  * Gelman & Rubin (1992) Gelman, A., & Rubin, D. B. 1992, Statistical Science, 7, 457, doi: 10.1214/ss/1177011136
  * Giesers et al. (2018) Giesers, B., Dreizler, S., Husser, T.-O., et al. 2018, MNRAS, 475, L15, doi: 10.1093/mnrasl/slx203
  * Giesers et al. (2019) Giesers, B., Kamann, S., Dreizler, S., et al. 2019, A&A, 632, A3, doi: 10.1051/0004-6361/201936203
  * Gomel et al. (2021a) Gomel, R., Faigler, S., & Mazeh, T. 2021a, MNRAS, 504, 2115, doi: 10.1093/mnras/stab1047
  * Gomel et al. (2021b) Gomel, R., Faigler, S., Mazeh, T., & Pawlak, M. 2021b, MNRAS, 504, 5907, doi: 10.1093/mnras/stab1235
  * Gomel et al. (2023) Gomel, R., Mazeh, T., Faigler, S., et al. 2023, A&A, 674, A19, doi: 10.1051/0004-6361/202243626
  * González & Levato (2006) González, J. F., & Levato, H. 2006, A&A, 448, 283, doi: 10.1051/0004-6361:20053177
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362, doi: 10.1038/s41586-020-2649-2
  * Hoffman & Gelman (2014) Hoffman, M. D., & Gelman, A. 2014, J. Mach. Learn. Res., 15, 1593
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/mcse.2007.55
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6, doi: 10.1051/0004-6361/201219058
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jayasinghe et al. (2020) Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2020, MNRAS, 491, 13, doi: 10.1093/mnras/stz2711
  * Joye & Mandel (2003) Joye, W. A., & Mandel, E. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook, 489
  * Kochanek et al. (2017) Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502, doi: 10.1088/1538-3873/aa80d9
  * Kuruwita et al. (2018) Kuruwita, R. L., Ireland, M., Rizzuto, A., Bento, J., & Federrath, C. 2018, MNRAS, 480, 5099, doi: 10.1093/mnras/sty2108
  * Lam et al. (2022a) Lam, C. Y., Lu, J. R., Udalski, A., et al. 2022a, ApJS, 260, 55, doi: 10.3847/1538-4365/ac7441
  * Lam et al. (2022b) —. 2022b, ApJ, 933, L23, doi: 10.3847/2041-8213/ac7442
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2, doi: 10.1051/0004-6361/202039709
  * Mahy et al. (2022) Mahy, L., Sana, H., Shenar, T., et al. 2022, A&A, 664, A159, doi: 10.1051/0004-6361/202243147
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the Python in Science Conference (SciPy), doi: 10.25080/majora-92bf1922-00a
  * Mitnyan et al. (2020) Mitnyan, T., Szalai, T., Bódi, A., et al. 2020, A&A, 635, A89, doi: 10.1051/0004-6361/201937214
  * Nagarajan et al. (2023) Nagarajan, P., El-Badry, K., Rodriguez, A. C., van Roestel, J., & Roulston, B. 2023, arXiv e-prints, arXiv:2304.07324, doi: 10.48550/arXiv.2304.07324
  * Olejak et al. (2020) Olejak, A., Belczynski, K., Bulik, T., & Sobolewska, M. 2020, A&A, 638, A94, doi: 10.1051/0004-6361/201936557
  * Orosz & Hauschildt (2000) Orosz, J. A., & Hauschildt, P. H. 2000, A&A, 364, 265, doi: 10.48550/arXiv.astro-ph/0010114
  * Perez & Granger (2007) Perez, F., & Granger, B. E. 2007, Computing in Science & Engineering, 9, 21, doi: 10.1109/mcse.2007.53
  * Pribulla et al. (2003) Pribulla, T., Kreiner, J. M., & Tremko, J. 2003, Contributions of the Astronomical Observatory Skalnate Pleso, 33, 38
  * Price-Whelan et al. (2017) Price-Whelan, A. M., Hogg, D. W., Foreman-Mackey, D., & Rix, H.-W. 2017, ApJ, 837, 20, doi: 10.3847/1538-4357/aa5e50
  * Prša & Zwitter (2005) Prša, A., & Zwitter, T. 2005, ApJ, 628, 426, doi: 10.1086/430591
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003, doi: 10.1117/1.JATIS.1.1.014003
  * Riello et al. (2021) Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3, doi: 10.1051/0004-6361/202039587
  * Rimoldini et al. (2023) Rimoldini, L., Holl, B., Gavras, P., et al. 2023, A&A, 674, A14, doi: 10.1051/0004-6361/202245591
  * Rodrigo & Solano (2020) Rodrigo, C., & Solano, E. 2020, in XIV.0 Scientific Meeting (virtual) of the Spanish Astronomical Society, 182
  * Rodrigo et al. (2012) Rodrigo, C., Solano, E., & Bayo, A. 2012, SVO Filter Profile Service Version 1.0, IVOA Working Draft 15 October 2012, doi: 10.5479/ADS/bib/2012ivoa.rept.1015R
  * Rowan et al. (2021) Rowan, D. M., Stanek, K. Z., Jayasinghe, T., et al. 2021, MNRAS, 507, 104, doi: 10.1093/mnras/stab2126
  * Sahu et al. (2022) Sahu, K. C., Anderson, J., Casertano, S., et al. 2022, ApJ, 933, 83, doi: 10.3847/1538-4357/ac739e
  * Salvatier et al. (2016) Salvatier, J., Wieckiâ, T. V., & Fonnesbeck, C. 2016, PyMC3: Python probabilistic programming framework.  http://ascl.net/1610.016
  * Samland (1998) Samland, M. 1998, ApJ, 496, 155, doi: 10.1086/305368
  * Saracino et al. (2022) Saracino, S., Kamann, S., Guarcello, M. G., et al. 2022, MNRAS, 511, 2914, doi: 10.1093/mnras/stab3159
  * Saracino et al. (2023) Saracino, S., Shenar, T., Kamann, S., et al. 2023, MNRAS, 521, 3162, doi: 10.1093/mnras/stad764
  * Shahaf et al. (2023) Shahaf, S., Bashi, D., Mazeh, T., et al. 2023, MNRAS, 518, 2991, doi: 10.1093/mnras/stac3290
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48, doi: 10.1088/0004-637X/788/1/48
  * Shenar et al. (2020) Shenar, T., Bodensteiner, J., Abdul-Masih, M., et al. 2020, A&A, 639, L6, doi: 10.1051/0004-6361/202038275
  * Shenar et al. (2022) Shenar, T., Sana, H., Mahy, L., et al. 2022, A&A, 665, A148, doi: 10.1051/0004-6361/202244245
  * Soszyński et al. (2016) Soszyński, I., Pawlak, M., Pietrukowicz, P., et al. 2016, Acta Astron., 66, 405, doi: 10.48550/arXiv.1701.03105
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757, doi: 10.48550/arXiv.1503.03757
  * Stevance et al. (2022) Stevance, H. F., Parsons, S. G., & Eldridge, J. J. 2022, MNRAS, 511, L77, doi: 10.1093/mnrasl/slac001
  * Tanikawa et al. (2022) Tanikawa, A., Hattori, K., Kawanaka, N., et al. 2022, arXiv e-prints, arXiv:2209.05632, doi: 10.48550/arXiv.2209.05632
  * Tetarenko et al. (2016) Tetarenko, B. E., Sivakoff, G. R., Heinke, C. O., & Gladstone, J. C. 2016, ApJS, 222, 15, doi: 10.3847/0067-0049/222/2/15
  * Thompson et al. (2019) Thompson, T. A., Kochanek, C. S., Stanek, K. Z., et al. 2019, Science, 366, 637, doi: 10.1126/science.aau4005
  * Thompson et al. (2020) —. 2020, Science, 368, eaba4356, doi: 10.1126/science.aba4356
  * Timmes et al. (1996) Timmes, F. X., Woosley, S. E., & Weaver, T. A. 1996, ApJ, 457, 834, doi: 10.1086/176778
  * Tody (1986) Tody, D. 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, Instrumentation in astronomy VI, ed. D. L. Crawford, 733, doi: 10.1117/12.968154
  * Tody (1993) Tody, D. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 173
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, Acta Astron., 65, 1, doi: 10.48550/arXiv.1504.05966
  * van den Heuvel (1992) van den Heuvel, E. P. J. 1992, Endpoints of stellar evolution: the incidence of stellar mass black holes in the Galaxy., In ESA, Environment Observation and Climate Modelling Through International Space Projects. Space Sciences with Particular Emphasis on High-Energy Astrophysics p 29-36 (SEE N93-23878 08-88)
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9, doi: 10.1051/aas:2000332
  * Wiktorowicz et al. (2019) Wiktorowicz, G., Wyrzykowski, Ł., Chruslinska, M., et al. 2019, ApJ, 885, 1, doi: 10.3847/1538-4357/ab45e6


MIDLMedical Imaging with Deep Learning 2019 MIDL 2019 – Extended Abstract
Track Grzegorz Chlebus1,2<EMAIL_ADDRESS>
Nasreddin Abolmaali3<EMAIL_ADDRESS>
Andrea Schenk1<EMAIL_ADDRESS>
Hans Meine1,4<EMAIL_ADDRESS>
1 Fraunhofer Institute for Digital Medicine MEVIS, Bremen, Germany
2 Diagnostic Image Analysis Group, Department of Radiology and Nuclear
Medicine, Radboud University Medical Center, Nijmegen, The Netherlands
3 Department of Radiology, Städtisches Klinikum Dresden, Dresden, Germany
4 University of Bremen, Medical Image Computing Group, Bremen, Germany

# Relevance analysis of MRI sequences for automatic liver tumor segmentation

###### Abstract

Explainability of decisions made by deep neural networks is of high value as
it allows for validation and improvement of models. This work proposes an
approach to explain semantic segmentation networks by means of layer-wise
relevance propagation. As an exemplary application, we investigate which MRI
sequences are most relevant for liver tumor segmentation.

###### keywords:

explainability, deep learning, segmentation, MRI

##  1 Introduction

Algorithms employing deep neural networks achieved state-of-the-art results
for many computer vision tasks [Hu et al.(2018)Hu, Shen, and Sun]. Thanks to
millions of parameters and nonlinear behavior deep models can learn very
complex input-output dependencies allowing them to surpass human expert
performance [Bejnordi et al.(2017)Bejnordi, Veta, Van Diest, Van Ginneken,
Karssemeijer, Litjens, Van Der Laak, Hermsen, Manson, Balkenhol, et al.].
Understanding how such models work is of big importance because it allows to
verify the reasoning of the system and possibly identify model or dataset
problems [Selvaraju et al.(2017)Selvaraju, Cogswell, Das, Vedantam, Parikh,
and Batra, Lapuschkin et al.(2019)Lapuschkin, Wäldchen, Binder, Montavon,
Samek, and Müller]. Several approaches for classifier explainability have been
proposed: guided backpropagation [Springenberg et al.(2014)Springenberg,
Dosovitskiy, Brox, and Riedmiller], layer-wise relevance propagation (LRP)
[Bach et al.(2015)Bach, Binder, Montavon, Klauschen, Müller, and Samek], and
gradient-weighted class activation mapping [Selvaraju et al.(2017)Selvaraju,
Cogswell, Das, Vedantam, Parikh, and Batra]. These methods aim at a
visualization of input regions influencing the model decision and were shown
to work well for image classification and reinforcement learning models
[Lapuschkin et al.(2019)Lapuschkin, Wäldchen, Binder, Montavon, Samek, and
Müller].

In this work, we propose a method to explain semantic segmentation models by
means of LRP. We apply the proposed method to analyze the importance of input
MRI sequences for the task of automatic liver tumor segmentation. Our
motivation is that models requiring fewer sequences would find a broader
clinical application, as not all hospitals employ the same acquisition
protocols. The importance analysis would allow for an informed selection of
most relevant sequences that could be used to train a segmentation model.

##  2 Materials and methods

###  2.1 Layer-wise relevance propagation

Layer-wise relevance propagation allows to obtain a pixel-wise decomposition
of a model decision and is applicable for most state-of-the-art architectures
for image classification and neural reinforcement learning [Bach et
al.(2015)Bach, Binder, Montavon, Klauschen, Müller, and Samek]. LRP uses a
notion of relevance $R$ which is equal to the model output for the output
neurons and can be propagated to lower layers.

####  2.1.1 LRP for image classification

LRP can be employed to explain classification decisions for a given class $i$
by relevance propagation from the corresponding model output $y^{i}$ according
to:

$y^{i}=R=\ldots=\sum_{d\in L^{l}}R_{d}^{(l)}=\ldots=\sum_{d\in
L^{1}}R_{d}^{(1)}=\sum M^{i}$ (1)

where $l$ refers to the layer index, $L^{l}$ to all neurons of layer $l$, and
$R_{d}^{(l)}$ to a relevance of neuron $d$ in layer $l$. Typically, the
relevances are propagated to the input layer ($l=1$) yielding a relevance map
$M^{i}$, which enables visualization of input regions influencing the model
decision.

####  2.1.2 LRP for semantic segmentation

In order to apply LRP to semantic segmentation models, we cast the
segmentation problem as a voxel-wise classification. This means that in order
to explain a decision of a segmentation model for a given output region $A$,
we propose to compute the input relevance maps according to Eq. 1 for each
considered output location $a\in A$. Then the relevance map $M^{i}$ explaining
the model decision for class $i$ in the region $A$ can be calculated as:

$M^{i}_{A}=\sum_{a\in A}\frac{M_{a}^{i}}{\sum M_{a}^{i}}$ (2)

We normalize $M_{a}^{i}$ by its sum to ensure that each output location $a$
equally contributes to the final relevance map $M^{i}$.

###  2.2 Liver tumor segmentation model

We train and evaluate a 3D u-net [Çiçek et al.(2016)Çiçek, Abdulkadir,
Lienkamp, Brox, and Ronneberger] model with a 6-channel input and 2-channel
output using MRI data of 69 patients (49 training, 20 evaluation) with primary
liver cancer and/or liver metastates acquired on a 3T MRI scanner (GE
Healthcare, USA) at Städtisches Klinikum Dresden, Germany. Imaging data of
each patient contains 6 MRI sequences: T2, non contrast enhanced T1
(plain-T1), and four dynamic contrast enhanced (DCE) T1 images acquired 20 s
(T1-20s), 60 s (T1-60s), 120 s (T1-120s), and 15 min (T1-15min) after contrast
agent administration (Gd-EOB-DTPA). Reference tumor segmentation was performed
manually by an experienced radiologist assistant on the last DCE phase. All
sequences were motion corrected using a non-rigid registration using the
T1-15min image as reference [Strehlow et al.(2018)Strehlow, Spahr, Rühaak,
Laue, Abolmaali, Preusser, and Schenk].

###  2.3 Relevance analysis of MRI sequences

We use Eq. 2 to compute a relevance map $M$ for both background ($i=0$) and
tumor ($i=1$) class: $M=M^{0}_{B}+M^{1}_{T}$. We choose regions $B$ and $T$
classified as background and tumor, respectively, to contain the same amount
of voxels to prevent bias towards the more frequent background class. The
importance of an input channel is corresponding to a global sum of a
corresponding channel in the relevance map $M$. As absolute values of $M$ have
no meaning, we normalize $M$ such that its values sum up to 1.

##  3 Results and discussion

The relevance distribution across input MRI sequences for 20 test patients is
shown in Fig. LABEL:fig:example. According to mean values, the most important
sequence was T1-120s and the least important T1-20s. The biggest differences
in attributed relevance were observed for T2 and T1-15min sequences. The
relevance attribution for other sequences follows almost the same pattern for
all test patients. The model learned to use information from all inputs as no
sequence received a zero relevance. The T1-15min sequence, which was used to
create reference segmentations, was not the most relevant for all test
patients, which was contrary to our expectations.

##  4 Conclusion

In this work, we proposed a method to explain semantic segmentation models by
means of layer-wise relevance propagation. We applied our method to analyze
which MRI sequences are the most important for the task of liver tumor
segmentation. Further applications of our method, for example to investigate
false predictions, are future work.

fig:example

Figure 1: Normalized relevance distribution across input MRI sequences for 20
test patients denoted by different colors.

We would like to thank Maximilian Alber, Sebatian Lapuschkin, and other
authors of the iNNvestigate toolbox [Alber et al.(2018)Alber, Lapuschkin,
Seegerer, Hägele, Schütt, Montavon, Samek, Müller, Dähne, and Kindermans].

## References

  * [Alber et al.(2018)Alber, Lapuschkin, Seegerer, Hägele, Schütt, Montavon, Samek, Müller, Dähne, and Kindermans] Maximilian Alber, Sebastian Lapuschkin, Philipp Seegerer, Miriam Hägele, Kristof T Schütt, Grégoire Montavon, Wojciech Samek, Klaus-Robert Müller, Sven Dähne, and Pieter-Jan Kindermans.  innvestigate neural networks!  _arXiv preprint arXiv:1808.04260_ , 2018.
  * [Bach et al.(2015)Bach, Binder, Montavon, Klauschen, Müller, and Samek] 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):e0130140, 2015.
  * [Bejnordi et al.(2017)Bejnordi, Veta, Van Diest, Van Ginneken, Karssemeijer, Litjens, Van Der Laak, Hermsen, Manson, Balkenhol, et al.] Babak Ehteshami Bejnordi, Mitko Veta, Paul Johannes Van Diest, Bram Van Ginneken, Nico Karssemeijer, Geert Litjens, Jeroen AWM Van Der Laak, Meyke Hermsen, Quirine F Manson, Maschenka Balkenhol, et al.  Diagnostic assessment of deep learning algorithms for detection of lymph node metastases in women with breast cancer.  _Jama_ , 318(22):2199–2210, 2017.
  * [Çiçek et al.(2016)Çiçek, Abdulkadir, Lienkamp, Brox, and Ronneberger] Özgün Çiçek, Ahmed Abdulkadir, Soeren S Lienkamp, Thomas Brox, and Olaf Ronneberger.  3d u-net: learning dense volumetric segmentation from sparse annotation.  In _International conference on medical image computing and computer-assisted intervention_ , pages 424–432. Springer, 2016.
  * [Hu et al.(2018)Hu, Shen, and Sun] Jie Hu, Li Shen, and Gang Sun.  Squeeze-and-excitation networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 7132–7141, 2018.
  * [Lapuschkin et al.(2019)Lapuschkin, Wäldchen, Binder, Montavon, Samek, and Müller] Sebastian Lapuschkin, Stephan Wäldchen, Alexander Binder, Grégoire Montavon, Wojciech Samek, and Klaus-Robert Müller.  Unmasking clever hans predictors and assessing what machines really learn.  _Nature communications_ , 10(1):1096, 2019.
  * [Selvaraju et al.(2017)Selvaraju, Cogswell, Das, Vedantam, Parikh, and Batra] Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pages 618–626, 2017.
  * [Springenberg et al.(2014)Springenberg, Dosovitskiy, Brox, and Riedmiller] Jost Tobias Springenberg, Alexey Dosovitskiy, Thomas Brox, and Martin Riedmiller.  Striving for simplicity: The all convolutional net.  _arXiv preprint arXiv:1412.6806_ , 2014.
  * [Strehlow et al.(2018)Strehlow, Spahr, Rühaak, Laue, Abolmaali, Preusser, and Schenk] Jan Strehlow, Nadine Spahr, Jan Rühaak, Hendrik Laue, Nasreddin Abolmaali, Tobias Preusser, and Andrea Schenk.  Landmark-based evaluation of a deformable motion correction for dce-mri of the liver.  _International journal of computer assisted radiology and surgery_ , pages 1–10, 2018.