# Multi-Messenger Constraints on Magnetic Fields in Merging Black Hole-Neutron
Star Binaries

Daniel J. D’Orazio Niels Bohr International Academy, Niels Bohr Institute,
Blegdamsvej 17, 2100 Copenhagen, Denmark<EMAIL_ADDRESS>Zoltán
Haiman Department of Astronomy, Columbia University, New York, NY, 10027, USA
Janna Levin Department of Physics and Astronomy, Barnard College of Columbia
University, New York, New York 10027, USA  Johan Samsing Niels Bohr
International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen,
Denmark  Alejandro Vigna-Gómez Niels Bohr International Academy, Niels Bohr
Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark  DARK, Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, 2200, Copenhagen, Denmark

###### Abstract

The LIGO-Virgo-KAGRA Collaboration recently detected gravitational waves (GWs)
from the merger of black-hole-neutron-star (BHNS) binary systems GW200105 and
GW200115. No coincident electromagnetic (EM) counterparts were detected. While
the mass ratio and BH spin in both systems were not sufficient to tidally
disrupt the NS outside of the BH event horizon, other, magnetospheric
mechanisms for EM emission exist in this regime and depend sensitively on the
NS magnetic field strength. Combining GW measurements with EM flux upper
limits, we place upper limits on the NS surface magnetic field strength above
which magnetospheric emission models would have generated an observable EM
counterpart. We consider fireball models powered by the black-hole battery
mechanism, where energy is output in gamma-rays over $\lesssim 1$ second.
Consistency with no detection by Fermi-GBM or INTEGRAL SPI-ACS constrains the
NS surface magnetic field to $\lesssim 10^{15}$ G. Hence, joint GW detection
and EM upper limits rule out the theoretical possibility that the NSs in
GW200105 and GW200115, and the putative NS in GW190814, retain $\gtrsim
10^{15}$ G dipolar magnetic fields until merger. They also rule out formation
scenarios where strongly magnetized magnetars quickly merge with BHs. We
alternatively rule out operation of the BH-battery powered fireball mechanism
in these systems. This is the first multi-messenger constraint on NS magnetic
fields in BHNS systems and a novel approach to probe fields at this point in
NS evolution. This demonstrates the constraining power that multi-messenger
analyses of BHNS mergers have on BHNS formation scenarios, the magnetic-field
evolution in NSs, and the physics of BHNS magnetospheric interactions.

Gravitational wave astronomy, Gamma-ray sources, Neutron star magnetic fields,
Black-hole-neutron-star binaries

##  1 Introduction

The LIGO-Virgo-KAGRA (LVK) collaboration recently announced the first
definitive detections of gravitational waves (GWs) from neutron-stars merging
with black-holes in the black-hole-neutron-star (BHNS) binaries GW200105 and
GW200115 (Abbott et al., 2021). These two events join the suspected BHNS
system GW190814 with inferred secondary mass on the cusp between
classification as a large NS and a small BH (Abbott et al., 2020). BHNS
systems have long attracted interest for their ability to generate louder GWs
than binary NSs (thanks to the BH), while still promising a bright
electromagnetic (EM) counterpart (thanks to the NS). Historically, the
disruption of the NS has been the focus of possible EM counterparts to BHNS
mergers (e.g., Narayan et al., 1992). However, disruption occurs only for a
small region of BHNS parameter space, when the BH is highly spinning and
small, $\lesssim 5-10{M_{\odot}}$ (see Foucart et al., 2018, and §3.1 below).
Furthermore, given the LVK population of compact object mergers discovered so
far, it is likely that most BHNS mergers will not disrupt the NS (Fragione,
2021).

However, as if often overlooked, the lack of a NS disruption does not imply
the lack of a bright EM counterpart, and hope is not lost for BHNS multi-
messenger astronomy. Bright EM emission is possible and even expected from
non-disrupting systems. The reason lies in the magnetic energy locked up in
the NS magnetosphere, which can be tapped without breaking apart the NS just
as it is tapped to power the many observational manifestations of known NSs,
e.g., the pulsars, anomalous x-ray pulses, soft-gamma repeaters, and others
(see, e.g., Kaspi, 2010).

A promising mechanism for brightening non-disrupting BHNS systems arises
through a coupling of the orbital and magnetic field energy of the merging
binary. Magnetic fields anchored in the NS pass over the BH horizon generating
an electromotive force that can supply power to accelerate charges in the
magnetosphere analogously to a unipolar inductor (homopolar generator). Such a
mechanism has been studied extensively in the literature, originally applied
to planetary systems (Goldreich & Lynden-Bell, 1969) and more recently applied
to exoplanetary systems (Laine & Lin, 2012), compact object mergers (Vietri,
1996; Hansen & Lyutikov, 2001; Lai, 2012; Piro, 2012) and systems including a
BH (McWilliams & Levin, 2011; Lyutikov, 2011; D’Orazio & Levin, 2013; D’Orazio
et al., 2016), in which case we refer to this energy source as the black-hole
battery. Encouragingly, numerical investigations of the BHNS magnetosphere
have found magnetospheric currents and Poynting flux measurements that agree
well with, and possibly exceed, estimates in the above, analytical works
(Paschalidis et al., 2013; Carrasco et al., 2021).

A key feature of BH-battery magnetospheric emission is its strong dependence
on the strength and configuration of the large-scale magnetic fields anchored
in the NS. The total power available for liberation via the BH-battery scales
with the square of the magnetic field strength at the BH horizon. Hence, if
the binary parameters are known from GW observations and consistent with a
non-disrupting system, then an EM observation consistent with a BH-battery
powered event provides information on the NS magnetic field as well as insight
into the BHNS magnetospheric physics. Similarly, a GW-confirmed non-disrupting
BHNS plus an EM upper limit places an upper limit on the magnetic field
anchored in the NS at merger.

Here we use this feature of the predicted BHNS magnetospheric emission, along
with GW observations and EM upper limits for LVK events GW200105, GW200115,
and GW190814 to make the first multi-messenger constraints of NS magnetic
field strengths in BHNS binaries. To do so we employ a model for the EM
signature of a BH-battery powered event developed in D’Orazio et al. (2016)
(hereafter DL16) which predicts a short-hard burst of gamma rays within
approximately a second after the merger. For each system we find that the
strength of a dipolar field at the NS surface is constrained to be $\lesssim
10^{15}$ G, within the astrophysically interesting realm of magnetar field
strengths. Future GW detections plus gamma-ray flux upper limits could offer
constraints as low as $10^{13-14}$ G depending on source distance, binary
parameters, and EM flux upper limits. Given the BHNS event rate, a few events
with meaningful flux upper limits are expected per year. This method offers
not only a new observational constraint on NS magnetic field evolution but
also a handle on BHNS formation channels.

This work is organized as follows. In §2 we provide an overview of the GW
observations and EM upper limits. In §3 we present a summary of the BH-battery
powered fireball model for non-disrupting magnetospheric emission in BHNS
systems. In §4 we apply this model to obtain magnetic field strength upper
limits for GW200105, GW200115, and GW190814. We discuss astrophysical
implications and prospects for the future in §5 and conclude in §6.

##  2 Overview of Observations

###  2.1 Gravitational Waves

The GW observations for GW200105 and GW200115 are detailed in Abbott et al.
(2021). We assume, as concluded in Abbott et al. (2021), that each binary
consists of a more massive primary BH component and a less massive secondary
NS component. For completeness we also consider GW190814 for which association
with a BHNS merger is likely but less certain given the large mass of the
secondary (Abbott et al., 2020; Antoniadis et al., 2021). Table 1 lists the
relevant binary parameters that we use in this work: BH mass
$M_{\mathrm{BH}}$, NS mass $M_{\mathrm{NS}}$, BH dimensionless spin
$S_{\mathrm{BH}}$ (ranging from -1,1), and luminosity distance $D_{L}$ 111To
convert between luminosity distance and redshift we assume the same cosmology
as Abbott et al. (2021), with $h=0.679$, $\Omega_{\mathrm{M}}=0.3065$.. For
all parameters we quote the median values with $90\%$ credible intervals. For
GW2001015 and GW200115 we assume the low-NS-spin prior of Abbott et al.
(2021).

Event |  $M_{\mathrm{BH}}$ [${M_{\odot}}$] |  $M_{\mathrm{NS}}$ [${M_{\odot}}$] | $S_{\mathrm{BH}}$ |  $D_{L}$ [Mpc]
---|---|---|---|---
GW200105 (Abbott et al., 2021) | $8.9^{+1.1}_{-1.3}$ | $1.9^{+0.2}_{-0.2}$ | $0.09^{+0.18}_{-0.08}$ |  $280^{+110}_{-110}$ ($z=0.06^{+0.02}_{-0.02}$)
GW200115 (Abbott et al., 2021) | $5.9^{+1.4}_{-2.1}$ | $1.4^{+0.6}_{-0.2}$ | $0.31^{+0.52}_{-0.29}$ |  $310^{+150}_{-110}$ ($z=0.07^{+0.03}_{-0.02}$)
GW190814 (Abbott et al., 2020) | $23.1^{+1.0}_{-1.1}$ | $2.59^{+0.09}_{-0.08}$ | $\lesssim 0.07$ |  $241^{+45}_{-41}$ ($z=0.053^{+0.009}_{-0.01}$)
Table 1: Double compact object parameters determined from gravitational wave
observations.

###  2.2 Electromagnetic Upper Limits

Here and in Table 2 we summarize the significant gamma-ray upper limits and
excesses reported for each event as reported in the GCN archives for GW200105
(GCN archive for S200105ae, 2020), GW200115 (GCN archive for S200115j, 2020),
and GW190814 GCN archive for LIGO/Virgo S190814bv (2020). For all three
events, no signal was detected by the Fermi-GBM above the background providing
a 3-sigma upper limit weighted by GW localization probability. Upper limits on
the gamma-ray flux were placed using the analysis pipeline of Blackburn et al.
(2015), updated by Goldstein et al. (2016). The upper-limit depends on the
expected signal duration and spectral template. The standard analysis pipeline
considers three spectral templates: ‘soft’, ‘normal’, and ‘hard’, and three
different durations: ‘short’ (0.128 s), ‘medium’ (1.024 s), and ‘long’ (8.192
s). We ignore the long-and soft combinations as these best describe long
gamma-ray bursts (GRBs) while the BHNS EM counterparts detailed below are
closer to hard, short GRBs. The ‘hard’ template assumes a Band-function (Band
et al., 1993) with a peak energy at $1$ MeV, and power law slopes $\alpha=0$,
$\beta=1.5$, while the ‘normal’ spectrum assumes a peak energy at $0.23$ MeV,
and power law slopes $\alpha=1$, $\beta=2.3$ (Blackburn et al., 2015). For
comparison, the energies of emission for the BHNS fireball (detailed below)
peak at $2.4\ \mathrm{MeV}\
\left(B_{\mathrm{NS}}/10^{14}\mathrm{G}\right)^{1/2}$ with approximate
durations of
$\sim$$0.1\mathrm{s}\left(B_{\mathrm{NS}}/10^{14}\mathrm{G}\right)$. Hence, we
choose to consider a range of upper limits given source templates that range
from medium-normal to short-hard, with the short-hard spectral template
providing the most conservative upper-limits. The resulting Fermi-GBM
(3-sigma) upper-limit ranges are listed in Table 2.

Table 2 also lists a flux upper limit for GW200105 from INTEGRAL SPI-
ACS222INTEGRAL did report one possibly associated excess with a $\sim 10^{49}$
erg/s/cm2 inferred luminosity, occurring $1.1$ sec before the merger with
signal-to-noise ratio of $3.3$ and a false-alarm probability of $0.01$.. This
limit is placed using a template in-between the Fermi medium-normal and short-
hard templates, with an exponentially cutoff $\alpha=-0.5$ power law and a
peak energy at 0.6 Mev, assuming a burst lasting $\leq 1$ s. This is
consistent with the range inferred from Fermi.

| Obs. Loc. Prob. | Time Window | Energy Range | Upper Limit (med.-norm. $-$ short-hard)
---|---|---|---|---
Detector | [$\%$] | [sec] | [MeV] | $\left[10^{-7}\mathrm{erg/s/cm}^{2}\right]$
| GW200105
Fermi-GBM | $83.2$ | $t_{\mathrm{mrg}}\pm 30$ | $10^{-2}-1$ | $1.8-11$
" | " | " | $10^{-3}-10$ | $2.5-25$
INTEGRAL SPI-ACS | $50.0$ | $t_{\mathrm{mrg}}\pm 300$ | $7.5\times 10^{-3}-2$ | 2.2
| GW200115
Fermi-GBM | $84.4$ | $t_{\mathrm{mrg}}\pm 30$ | $10^{-2}-1$ | $2.8-17$
" | " | " | $10^{-3}-10$ | $4.0-40$
| GW190814
Fermi-GBM | $100$ | $t_{\mathrm{mrg}}\pm 30$ | $10^{-2}-1$ | $0.9-5.8$
" | " | " | $10^{-3}-10$ | $1.3-13$
Table 2: Gamma-ray upper limits from GCN archives.

##  3 Non-Disruption and Emission Model

Given the binary properties inferred from GW observations, we will use the EM
upper limits and a model for EM emission to place constraints on the NS
magnetic field at merger. Hence, we next demonstrate that the NSs in GW200105,
GW200115, and GW190814 were very likely not disrupted by their respective BHs
and review the model of DL16 for BH-battery powered magnetospheric emission in
non-disrupting BHNS systems.

###  3.1 Disruption or not?

$\begin{array}[]{c}\includegraphics[scale={0.56}]{MinxSBH_to_disrupt_MBHvsMNSFouc}\end{array}$

Figure 1:  White-to-black contours, with specific values marked by the dashed
curves as labeled, represent the maximum aligned component of the BH spin for
which a BHNS system can be swallowed whole, given BH mass (y-axis), NS mass
(x-axis), and NS radius (fixed at $R_{\mathrm{NS}}=10$ km). The teal (orange)
dashed lines represent the 90$\%$ confidence upper limits on the BH spin from
GW200105 (GW200115) from Table 1. For example, systems lying on the white line
labeled $0.9$ will be swallowed whole if the orbit-aligned BH spin component
is below $0.9$ – they will at least partially disrupt for aligned BH spin
component above $0.9$. Solid red contours represent the required aligned spin
component above which a full disruption occurs. Teal (Orange) data points with
errorbars represent the values quoted in Table 1. GW190814 would lie to the
far upper right, beyond the plotted domain.

What portion of allowed binary parameter space results in NS disruption
outside of the BH horizon? We focus on GW200105 and GW2100115 as the
$23{M_{\odot}}$, non-spinning BH in GW190814 precludes it from any likely NS
disruption. In Figure 1, we use the criteria of Foucart et al. (2018),
calibrated using General Relativistic simulations of BHNS mergers on circular
orbits, to estimate the BH and NS parameters for which NS disruption occurs.
Given BH and NS masses and a choice of $R_{\mathrm{NS}}=10$ km, we solve for a
maximum value of the orbit-aligned component of the BH spin,
$S_{\mathrm{BH}}$, below which no remnant NS mass is liberated and escapes
beyond the BH event horizon. These represent maximum aligned spin components
above which at least a partial disruption would occur and are represented in
Figure 1 by the white-to-black contours, with selected values marked by the
dashed curves, as labeled. The two data points and associated errorbars
represent the median and $90\%$ confidence uncertainties on the BH and NS
masses in GW200105 (teal) and GW200115 (orange; see Table 1). In addition to
the BH spins below which the NS is swallowed whole, we also solve for the
aligned component of the black hole spin above which a full disruption would
occur. We draw red contours for these spins at values that pass through the
median masses for GW200105 and GW200115.

The dashed teal (orange) contours correspond to the $90\%$ confidence upper
bound on the BH spins of GW200105 (GW200115) as inferred from the GW
measurements (Table 1). For GW200105, the dashed teal contour demonstrates
that the maximum allowed spin magnitude at $90\%$ confidence of $0.27$ is well
below the value of $\sim$$0.9$ that is required to disrupt the NS. Hence, no
allowed combination of the BH spin and component masses results in even a
partial disruption for GW200105. For GW200115, no disruption will occur for
the median mass and spin values, but partial disruption can occur for the
largest allowed BH spin values (the dashed orange contour at
$S_{\mathrm{BH}}=0.83$) and all but the upper limits on the component masses,
assuming spin-orbit alignment. Any misalignment makes the NS even more
difficult to disrupt. The red contours show that a full disruption occurs only
if $S_{\mathrm{BH}}\gtrsim 0.885$ for GW200115, or $S_{\mathrm{BH}}\gtrsim
0.997$ for GW200105. Hence, neither system could have been fully disrupted
given the $90\%$ confidence system parameters. Hence, a short-GRB-like event
due to NS disruption is unlikely.

###  3.2 BH-Battery Powered Fireball

We next review the fireball model of DL16. In this model, a relative velocity
between the NS-sourced magnetic fields and the BH horizon generates an
electromotive force in analogy to a unipolar inductor (Goldreich & Lynden-
Bell, 1969). The voltage drop corresponding to this BH-battery can be written
in terms of binary parameters and the NS magnetic field strength as
(McWilliams & Levin, 2011)

$\displaystyle
V_{\mathcal{H}}=2R_{\mathcal{H}}\left[\frac{r\left(\Omega_{\rm{orb}}-\Omega_{\mathrm{NS}}\right)}{c}+\frac{S_{\mathrm{BH}}}{4\sqrt{2}}\right]B_{\rm
NS}\left(\frac{R_{\mathrm{NS}}}{r}\right)^{3},$ (1)

where $r$ is the binary separation, $R_{\mathcal{H}}$ is the spin and mass
dependent BH horizon radius,
$\Omega_{\rm{orb}}=\sqrt{G(M_{\rm{BH}}+M_{\rm{NS}})/r^{3}}$,
$\Omega_{\mathrm{NS}}$ is the NS spin angular frequency, $S_{\mathrm{BH}}$ is
the dimensionless BH spin, and we have assumed a dipolar magnetic field with
NS surface field strength $B_{\rm NS}$. The power available to the BH-battery
is given by Ohm’s law, the horizon resistivity
$\mathcal{R}_{\mathcal{H}}=4\pi/c$ and the impedance-matching condition
$\mathcal{R}_{\mathcal{H}}=\mathcal{R}_{\rm NS}$ (see, e.g., Thorne et al.,
1986),

$\mathcal{P}(r)=\frac{2V^{2}_{\mathcal{H}}(r)}{(\mathcal{R}_{\mathcal{H}}+\mathcal{R}_{\rm
NS})^{2}}\mathcal{R}_{\rm NS}\rightarrow\frac{c}{8\pi}V^{2}_{\mathcal{H}}(r),$
(2)

where the prefactor counts emission from each hemisphere over which a full
voltage drop occurs.

DL16 show that at $\sim$$10\left(B_{\rm NS}/10^{14}\mathrm{G}\right)$ seconds
before merger, the magnetosphere becomes opaque to pair-production, primarily
from photon-magnetic-field ($\gamma-B$) interactions. Hence, after this point,
the BH-battery energy is injected into a hot pair-plasma. Because of the steep
$\sim$$r^{-4}$ scaling of BH-battery power, the injected power can be
evaluated as $\mathcal{P}(r_{\mathrm{mrg}})$, where we take the separation at
merger to be $r_{\mathrm{mrg}}=R_{\mathrm{NS}}+R_{\mathcal{H}}$ (see §5.3.2).
The optically thick pair plasma expands under its own pressure until it
reaches a radius $R_{\mathrm{ph}}$ a few
$\times\sqrt{B/(10^{14}~{}\mathrm{G})}$ milliseconds after merger where it
becomes optically thin to Thomson electron scattering and emits as a black
body with an effective temperature given approximately by,

$T_{\mathrm{eff}}=\left(\frac{\mathcal{P}(r_{\mathrm{mrg}})}{4\pi
r^{2}_{0}\sigma_{\mathrm{SB}}}\right)^{1/4},$ (3)

where $\sigma_{\mathrm{SB}}$ is the Stefan-Boltzmann constant, and the initial
injection radius is taken to be $r_{0}=r_{\mathrm{mrg}}$ to approximate the
size of the magnetosphere region surrounding the binary at merger (there is
very little dependence on this choice). This effective temperature takes into
account the adiabatic cooling of the fluid as it expands to $R_{\mathrm{ph}}$
and the counteracting effect (in the observer’s frame) of the Doppler boost
from the relativistically expanding fireball (see below and DL16).

The resulting spectrum is found from integrating the multi-component black-
body emission over a sphere of radius $R_{\mathrm{ph}}$ that is expanding with
Lorentz factor $\gamma=R_{\mathrm{ph}}/r_{0}$ and with rest frame photosphere
temperature $T_{0}/\gamma$, for initial injection temperature $T_{0}$. This
calculation is carried out in full in DL16, however, for simplicity here we
remove the integral over the sphere and approximate the observed flux with the
most significantly contributing portion of the emitting photosphere, the
portion within half opening angle $\gamma$ and a constant observed temperature
at redshift $z$ of $\gamma
T_{\mathrm{ph}}(1+z)^{-1}=\gamma\left(T_{0}/\gamma\right)(1+z)^{-1}=T_{0}(1+z)^{-1}$
(see, e.g., Paczynski, 1986, and Eq. 3),

$F_{\mathrm{FB},\rm{obs}}\approx
2\pi\mathcal{K}\left(1-\cos\gamma^{-1}\right)\left(\frac{R_{\mathrm{ph}}}{D_{L}}\right)^{2}\int^{\nu_{\rm{max}}}_{\nu_{\rm{min}}}{\frac{2h\nu^{3}}{c^{2}}\frac{d\nu}{\rm{exp\left[\frac{h\nu}{kT_{0}}\right]-1}}},$
(4)

which is the flux from a blackbody with temperature $T_{0}$ at redshift $z$,
corresponding to luminosity distance $D_{L}$333Written in terms of angular
diameter distance and $T_{0}(1+z)^{-1}$ in DL16.. Here $\nu$ is the frequency
in the observer’s frame, $2\pi\left(1-\cos\gamma^{-1}\right)$ is a geometrical
factor accounting for the significantly boosted region of the sphere, and
$\mathcal{K}$ is a numerical factor of order unity that corrects the single-
temperature approximation that allowed removal of the integral over the
hemisphere 444Under this single-temperature approximation the spectrum is
blackbody but in the full calculation, deviations from blackbody arise because
the observed temperature varies across the surface of the photosphere due to
the varying line of sight expansion velocity (DL16).. For magnetic field
strengths between $10^{12}$ and $10^{16}$, and binary parameters considered
here, $\mathcal{K}\approx 1.6$ to within a few $\%$.

$\begin{array}[]{cc}\includegraphics[scale={0.5}]{Constriant_GW200105_Norm-
mid_MNSeq1p9Msun_spin0p27}&\includegraphics[scale={0.5}]{Constriant_GW200105_Hard-
Short_MNSeq1p9Msun_spin0}\vspace{-23pt}\\\
\includegraphics[scale={0.5}]{Constriant_GW200115_Norm-
Mid_MNSeq1p4Msun_spin0p83}&\includegraphics[scale={0.5}]{Constriant_GW200115_Hard-
short_MNSeq1p4Msun_spin0}\vspace{-23pt}\\\
\includegraphics[scale={0.5}]{Constriant_GW190814_Norm-
Mid_MNSeq2p59Msun_spin0p07}&\includegraphics[scale={0.5}]{Constriant_GW190814_Hard-
Short_MNSeq2p59Msun_spin0}\end{array}$

Figure 2:  NS surface magnetic-field-strength constraints for GW200105 (top),
GW200115 (middle), and GW190814 (bottom). Grey-dashed contours represent the
magnetic field strength required to generate a signal at the Fermi-GBM upper
limit in the $10^{-3}-10$ MeV range. The dotted lines connecting
$M_{\rm{BH}}-D_{L}$ error bars represent the $90\%$ confidence contours for a
bivariate Gaussian distribution approximating the joint $M_{\rm{BH}}-D_{L}$
posteriors. The fiducial NS parameter values are the median $M_{\mathrm{NS}}$
from Table 1, $R_{\mathrm{NS}}=10$ km, and $\Omega_{\mathrm{NS}}=0$. The
thick-purple contour is the corresponding $90\%$ confidence upper limit on the
magnetic field strength. In each panel, the arrows point in the direction of
allowed magnetic field strength. The left column represents an optimistic
estimate and uses the medium-normal template spectrum and the maximum value of
BH spin allowed from the GW observations. The right column represents a more
conservative estimate and uses the short-hard template spectrum and zero BH
spin. The two scenarios bracket the range of magnetic field strength upper
limits, $B_{\rm{max}}$ in Table 3.

##  4 Results

The total flux predicted within a specified spectral range is determined by
the parameters
$(B_{\mathrm{NS}},M_{\mathrm{NS}},R_{\mathrm{NS}},\Omega_{\mathrm{NS}},M_{\mathrm{BH}},S_{\mathrm{BH}},D_{L})$.
For the NS radius we choose $R_{\mathrm{NS}}=10$ km. For the unconstrained NS
spin we conservatively choose zero spin as this is likely consistent with
astrophysical channels discussed in §5.1, which may bring a NS with strong
fields to merger. The results are not strongly dependent on the NS spin unless
it matches the orbital frequency of the binary (see Eq. 1). If this occurs in
the retrograde sense, then the available power increases by a factor of four.
In the prograde, i.e. tidally locked case, the BH power drops significantly,
however, tidal locking is not expected for these systems (Bildsten & Cutler,
1992; Lai, 2012).

Provided these NS radius and spin parameters, a flux upper limit
$F_{\mathrm{UL}}$ over frequency range
$\left[\nu_{\rm{min}},\nu_{\rm{max}}\right]$, and with the GW-measured ranges
on $(M_{\mathrm{NS}},M_{\mathrm{BH}},S_{\mathrm{BH}},D_{L})$, we solve
$F_{\mathrm{UL}}\leq F_{\mathrm{FB},\rm{obs}}(B_{\rm{max}})$ for the maximum
allowed NS surface magnetic field strengths. The most significant constraints
come from the $10^{-3}-10$ MeV flux upper limits as these extend to the high
energies where the fireball luminosity peaks. Hence, we use the corresponding
Fermi-GBM flux upper limits in Table 2 and choose $h\nu_{\rm{min}}=10^{-3}$
MeV and $h\nu_{\rm{max}}=10$ MeV in Eq. (4).

Figure 2 illustrates the results of this calculation in $M_{\mathrm{BH}}$ vs
$D_{L}$ space as these parameters most greatly affect the result given their
uncertainties and dependence in the flux calculation. Each row in Figure 2 is
for a different GW source, and the columns bracket the range of optimistic and
pessimistic choices for flux upper limit and BH spin. The higher the prograde
spin, the closer together the pair gets at merger, and the higher the power
input (Eq. 1). Hence, the highest spins allowed by the GW observation are used
in the optimistic cases, while we choose zero spin in the pessimistic case
(consistent with the lowest observationally allowed BH spins). In each panel
of Figure 2 we plot (dashed-grey) contours of $\log_{10}B_{\mathrm{NS}}$
required to generate an EM signal with the indicated choice of flux upper
limit and BH spin. For each panel we use the median value of $M_{\mathrm{NS}}$
since the allowed range quoted in Table 1 only affects our result for
$B_{\rm{max}}$ at the percent level.

To place the upper limit $B_{\rm{max}}$, we approximate the 2D posterior on
$M_{\mathrm{BH}}$ and $D_{L}$ as a bivariate Gaussian. The relevant $90\%$
confidence contours are plotted as the dotted quarter-ellipses in Figure 2. We
then find the values of $(M_{\mathrm{BH}},D_{L})$ which maximise
$B_{\rm{max}}$ along this ellipse giving the $90\%$ confidence upper limits
for $B_{\rm{max}}$. This value is denoted in each panel of Figure 2 by the
thick purple contour. From the pessimistic and optimistic values for each GW
source we generate a range of inferred maximum magnetic field values and list
these in Table 3. The majority of the difference between the optimistic (opt.)
and pessimistic (pes.) results in Table 3 is due to the change in the flux
upper limit (factor of $\sim$$\sqrt{10}$) as opposed to the change in spin
which contributes at most a factor of $\sim$$1.6$ in the case of GW200115,
where the spin ranges from $0.0-0.83$.

|  $B_{\rm{max}}$ [G] (opt. $-$ pes.)
---|---
GW200105 | $6.3\times 10^{14}-2.2\times 10^{15}$
GW200115 | $4.4\times 10^{14}-2.3\times 10^{15}$
GW190814 | $1.3\times 10^{15}-4.2\times 10^{15}$
Table 3: Range of magnetic field upper limits $B_{\rm{max}}$ using Fermi-GBM
upper limits over $10^{-3}-10$ MeV. Optimistic (opt.) and pessimistic (pes.)
values bracket the range as described in Fig. 2.

$\begin{array}[]{cc}\includegraphics[scale={0.53}]{Bmax_vs_rmrg_src3.pdf}\end{array}$

Figure 3:  Variation in NS magnetic field upper limit as a function of model
uncertainties. The x-axis varies the value of the separation where we assume
that the BH-battery power shuts off and the fireball begins expanding (see Eq.
3). The dashed lines correspond to the pessimistic case where
$S_{\mathrm{BH}}=0$, and the short-hard flux upper limit is used. The solid
lines correspond to the optimistic case, which employs $90\%$ confidence upper
limits on the BH spin and the medium-normal flux upper limits. The dots
represent the fiducial cases as plotted in Figure 2 and quoted in Table 3.

##  5 Discussion

###  5.1 Astrophysical Implications

We have shown that GW observations of BHNS binaries plus EM upper limits can
constrain the strength of the magnetic field anchored to the NS as it takes
its final plunge across the BH horizon. For the first two definitive BHNS
binary detections, an upper limit on the maximum NS field strength can be
conservatively placed at $\sim$$10^{15}$ G, and slightly higher for GW190814.
These values are at the high end of observed magnetar magnetic field strengths
and within theoretically allowed values. We now assess the astrophysical
implications of these and plausible future upper limits.

####  5.1.1 Long-Lived or Rejuvenated Fields

Theoretical arguments depending on hydro-magnetic instabilities that grow the
field during the proto-NS stage predict maximum initial fields strengths up to
$\sim 5\times 10^{16}$ G (Miralles et al., 2002). Theoretical work on the
coupled crust-core magnetic-field evolution in NSs is still under development,
but current uncertainty in our understanding leaves open the possibility that
large-scale magnetic fields anchored in the superconducting core could persist
indefinitely (Lyutikov & McKinney, 2011; Elfritz et al., 2016; Bransgrove et
al., 2018; Gusakov et al., 2020). Furthermore, while elegantly motivated by
the picture of NS unification (Kaspi, 2010), there is yet no direct evidence
for large-scale dipole-field decay in NSs (see Igoshev et al., 2021, for an
overview of existing indirect evidence). Hence, the theoretical possibility of
long-lived, strong magnetic fields anchored in NS cores combined with a lack
of direct evidence for fast dipolar field decay suggests the possibility that
NSs could be born with and retain $\lesssim 10^{16}$ G poloidal magnetic field
components that they take with them to merger. Indeed the only way to release
the energy in these core-anchored fields and so provide evidence for their
existence would be for a resulting NS collision with another compact object,
either through disruption, or through tapping these fields with the black-hole
battery.

Hence, the multi-messenger constraints considered here rule out the existence
of long-lived $\gtrsim 10^{15}$ G dipolar magnetic-fields retained on the
merging NSs in GW200105 and GW200115. An estimate of the fraction of NSs
expected to be born with fields above $10^{15}$ G, however, is beyond the
scope of this work. Similarly, we rule out any mechanism that pumps the large
scale poloidal magnetic-field component to high values during inspiral555See
for example the mechanism in Zrake & MacFadyen (2013), which, however,
generates a turbulent not and an ordered field, and likely requires full NS
disruption..

####  5.1.2 Magnetar Plus Black Hole Merger Scenarios

Taking a more conservative direction, we can employ the existence of the
observed population of high-magnetic-field NSs, also known as magnetars.
Magnetars have field strengths up to a few times $10^{15}$ G, inferred from
spin down arguments (see Igoshev et al., 2021, and references therein) and
also more direct cyclotron measurements (Tiengo et al., 2013). From a
combination of arguments involving kinematic association with a supernova
remnant or young star cluster, NS galactic position, or NS spin, one can infer
typical magnetar lifetimes of $\lesssim 10^{5}$ yrs (see Igoshev et al., 2021,
and references therein). Hence, under the assumption that the magnetic field
components responsible for magnetar behavior disappear at the end of the
magnetar lifetime, and that these are the large scale fields that power the
black-hole battery, then the upper limits derived here rule out GW200105 or
GW200115 as products of a (strong-field)-magnetar-BH merger. A necessary
condition for such a merger is that the time between magnetar formation and
merger is shorter than the magnetar lifetime of $10^{4}-10^{5}$ yrs. It is
then natural to ask what formation channels, if any, can act on this timescale
and can so be ruled out for these and similar future events.

##### Dynamical Pairing

If a magnetar is formed in a dense cluster, such as a nuclear star cluster, it
can be brought to merge with a BH through various dynamical pathways. One is
through direct GW capture, where the initially unbound NS becomes bound to a
BH through the emission of GWs at its first pericenter passage. Considering
now a system with constant velocity dispersion, $v$, and uniform density of
BHs, $n$, the characteristic time for a NS to be captured by a BH is
$t_{c}\sim 1/(n\sigma_{c}v)$, where $\sigma_{c}$ is the capture cross section
(e.g., Lee, 1993). For a population of BHs with mass $M=10M_{\odot}$, and a NS
with mass $m=1.4M_{\odot}$ one then finds $t_{c}\approx 10^{12}\
\mathrm{yrs}\times\ (v/(10\
\mathrm{kms}^{-1}))^{11/7}(n/10^{5}\mathrm{pc}^{-3})^{-1}$. Therefore, pairing
up magnetars through GW capture chance encounters in dense clusters is highly
unlikely. The same conclusion, but with different scalings, is found in the
case the merger is mediated through a 3-body interaction between a binary BH
and the NS, which suggests that BH-magnetar mergers are not easily assembled
through standard chaotic dynamics, unless more exotic pathways are considered
(see, e.g., the active galactic nuclei scenario McKernan et al. 2020; Yang et
al. 2020a; Tagawa et al. 2021).

That the merger rate involving NSs in stellar clusters generally is low
(Mandel & Broekgaarden, 2021) has also recently been shown using state-of-the-
art Monte Carlo cluster codes (Ye et al., 2020).

##### Binary Stellar Evolution

The merger rates of BHNS binaries are predicted to be dominated by massive
binaries that evolved in isolation in contrast to those that evolve in
dynamical environments (see Mandel & Broekgaarden, 2021, for a review).
Therefore, we use a synthetic population created with the rapid population
synthesis element of the COMPAS (Team COMPAS et al., 2021; Stevenson et al.,
2017; Vigna-Gómez et al., 2018) suite v02.19.01 (available via Vigna-Gómez et
al., 2021) to estimate the merger time distribution of BHNSs from isolated
binary evolution. The population consists of $10^{7}$ massive binaries
($M>5M_{\odot}$) at representative metallicities
$Z=\\{0.02,0.0142,0.01,0.001,0.0001\\}$ following the setup from Vigna-Gómez
et al. (2020). We only consider BHNS binaries in which the NS is the second
compact object to form.

The number of BHNS binaries per total mass evolved, the yield, is
approximately $10^{-5}M_{\odot}^{-1}$ for all metallicities. The fraction of
BHNSs with merger times $<0.1$ Myr is $\lesssim 2/1000$. In this population,
short-lived BHNS binaries are the product of exotic formation channels and/or
fortuitous natal kicks. This fraction rises to $\lesssim 1/100$ for merger
times below a few tens of Myr, and $\lesssim 1/10$ for for merger times below
a few hundreds of Myr.

Hence, if we assume that magnetars can be formed in binary systems (Xu et al.,
2021) and that the NS can only hold onto large magnetic fields for the
magnetar lifetime, then $\sim$$10^{3}$ BHNS mergers formed through isolated
binary evolution are required within a detectable EM horizon before the signal
described here is expected. On the other hand, we can use this calculation to
estimate that after observing $\sim$$10$ BHNS mergers via GWs within a
detectable EM horizon and with no EM counterpart, we begin to constrain large
magnetic field retention timescales to less than $\sim$$10^{8}$ yrs. Below we
estimate how many such events within a detectable EM horizon are expected with
present-day instruments.

Note that both formation scenarios above assume that the magnetar lifetime
spans a continuous period starting from NS formation, but there is still much
to learn about magnetar life-cycles. There is evidence for evolution into and
out of the magnetar stage suggested by some radio pulsars for which a growing
poloidal magnetic field is inferred via measurements of the breaking index.
For example, at its current rate of inferred field growth, PSR 1734-3333 would
be reclassified as a magnetar in a few tens of thousands of years (Espinoza et
al., 2011), suggesting the possibility of a bidirectional evolution between
the two classes. Additionally, magnetic fields can be buried by gas accretion
and re-emerge later, suggesting the possibility of “hidden magnetars" Geppert
et al. (1999). These possibilities could increase the chances for magnetar-BH
mergers by delaying magnetar field strengths to a later time after magnetar
formation.

###  5.2 Prospects for future observations

In the future, when we have $N$ total constraints similar to those considered
here, we will be able to establish that $\lesssim 1/N$ BHNS mergers do not,
for example, contain NSs with long-lived dipole fields above a given strength.
Above we estimated what this fraction might be for some viable astrophysical
scenarios; here we estimate what $N$ may grow to given current merger rates
and EM horizons.

##### How many do we expect to see in coming years?

Assuming that these two events are representative of the population, the BHNS
merger rate is estimated at $45^{+75}_{-33}$ Gpc-3yr-1, or higher
$130^{+112}_{-69}$ Gpc-3yr-1 when assuming a broader mass distribution (Abbott
et al., 2021). Considering that the majority of the BHNS parameter space
represented in Figure 1 does not result in NS disruption, we assume that most
of the coming BHNS events will not result in a NS disruption (Fragione, 2021)
and provide additional upper limits on the NS magnetic fields strength (if not
an actual EM counterpart). Assuming that useful limits can be made for systems
out to 300 Mpc, we expect $1.2^{+2.0}_{-0.9}$ or $3.5^{+3.0}_{-1.9}$ events
per year of observing time. For even more useful, closer events within 100
Mpc, we expect $0.05^{+0.08}_{-0.03}$ or $0.13^{+0.11}_{-0.07}$ events after
one year of observations. Hence a few events similar to those studied here are
expected in a coming year of observations and ten years may be needed to
detect a sub-100 Mpc system.

##### What field strengths can we expect to probe?

From Figure 2 we see that moving GW2001105 or GW2001115 closer, to within 100
Mpc can decrease the upper limit to $\sim$$10^{14}$ G. Because the magnetic
field strength upper limit scales as the square root of the flux upper limit,
a future BHNS merger with similar properties as GW2001105 could result in a
$10^{13}$ G constraint for a $100\times$ better flux upper limit than used
here. More robust upper limits may be possible in the future not only with
more sensitive EM detectors but with GW pre-warnings from low-frequency
detectors such as LISA (Sesana, 2016; Amaro-Seoane & et al., 2017), TianQin
(Luo et al., 2016; Liu et al., 2020), or DECIGO (Isoyama et al., 2018).

###  5.3 Caveats

####  5.3.1 Observational Constraints

Recently Mandel & Smith (2021) carried out a re-analysis of the GW200115
event, specifically on the recovery of mass and spins using a more
astrophysically informed prior on the BH spin. Mandel & Smith (2021) recover a
BH spin that is consistent with zero and a tighter constraint on the BH mass
of $7.0^{+0.4}_{-0.4}{M_{\odot}}$. This lower high-end of the BH mass estimate
results in a decreased $B_{\rm{max}}$ in the pessimistic zero spin case (see
Figure 2). However, in the optimistic maximum spin case, reducing the BH spin
from $0.83$ to $0.0$ has the overall effect of increasing $B_{\rm{max}}$ for
GW200115. This results in a $B_{\rm{max}}$ range for GW200115 of $8.1\times
10^{14}~{}\mathrm{G}\lesssim B_{\rm{max}}\lesssim 2.6\times
10^{15}~{}\mathrm{G}$ when using the binary component values inferred in
Mandel & Smith (2021).

Throughout, we have used the low-NS-spin prior for the GW parameter values. If
we instead use the high-spin priors of Abbott et al. (2021), the largest
change in the recovered system parameters is at the $\sim$0.1${M_{\odot}}$
level for the BH mass, resulting in an insignificant change to our magnetic
field upper limits. A larger, though still minimal effect comes from including
a non-zero $\Omega_{\mathrm{NS}}$, as discussed in §4.

As discussed in §2.2, the gamma-ray flux upper limits are derived from
temporal-spectral models of the expected emission, based on gamma-ray bursts.
We used upper limits that derived from choices of the spectrum and signal
duration which most closely match the fireball emission models. These are,
however, not tailored to the fireball models. Given uncertainties in the
fireball model that we discuss below, however, we expect that the chosen upper
limits are sufficient for the level of approximation here.

####  5.3.2 Fireball Model

Throughout we have assumed that the total energy injected into the fireball is
the full BH-battery power at a separation
$r_{\mathrm{mrg}}=R_{\mathrm{NS}}+R_{\mathcal{H}}$. This is motivated by the
final binary separation where magnetic fields from the NS can be dragged
across the BH horizon and generate the voltage of Eq. (1). Figure 3 shows how
our results for the magnetic field upper limits varies for a range of choices
for $r_{\mathrm{mrg}}$. For the largest $r_{\mathrm{mrg}}$, we conservatively
choose the ISCO of a test mass orbiting a non-spinning BH. For the smallest
$r_{\mathrm{mrg}}$, we choose the natural scale of the problem, the horizon
radius. At such small separations one might expect that the BH-battery
operation as we describe it here breaks down or that power is swallowed by the
BH. Meanwhile, the larger, non-spinning ISCO separations may be too
conservative as the BH-battery may not yet be out of power at this earlier
point in the merger. We make our fiducial choice (denoted by the dots in Fig.
3) because even if the BHNS is no longer on adiabatically decaying circular
orbits below the ISCO (due also to tidal and higher order relativistic
effects), it still has relativistic motion that continues to generate an
electromotive force and power the BH-battery (e.g., D’Orazio & Levin, 2013).
For highly spinning BHs, the ISCO approaches $R_{\mathcal{H}}$ and is smaller
than our fiducial $r_{\mathrm{mrg}}$. Further understanding of the BH-battery
powered emission may eliminate the need to include this parameter or help to
choose its value more accurately.

We have not considered mis-alignment of the BH spin, NS magnetic field, and
orbital angular momentum (however, see Carrasco et al., 2021) nor have we
considered non-dipolar magnetic field configurations. We have not considered
beaming of the signal, which would increase the volume to which events could
be observed, while decreasing the fraction of BHNS systems for which
magnetospheric emission would be observed. Thus far, numerical force-free
General-Relativistic simulations of BHNS magnetospheres predict only weak
emission anisotropy – Paschalidis et al. (2013) find that the Poynting flux
carries away energy in a broad beam that could result in a lighthouse-like
effect.

####  5.3.3 Non-Fireball Emission

We have focused on magnetospheric emission via the DL16 fireball model because
it relies on relatively simple, robust energetics arguments and predicts a
frequency dependent flux with which we can compare with observed upper limits.
However, just as the escaping radiation can take on many forms in the
analogous pulsar magnetosphere so could it for the BHNS binary magnetosphere.
For example, Mingarelli et al. (2017) argue that the BH-battery could power
radio emission reminiscent of fast radio bursts. Carrasco et al. (2021) carry
out global force-free simulations of a NS with a pure dipole magnetic field
orbiting in a Kerr spacetime finding $\sim$$10\times$ BH-battery luminosities
and finding that flux tubes transport electric currents similarly to what is
predicted in the analytical BH-battery models (see also Paschalidis et al.
2013), finding additionally that energy can be transported outwards by re-
connection and large scale Alfven waves. They consider EM signatures at radio
frequencies, but other manifestations of this energy transport are possible
and require further study. A number of recent works have also begun to develop
possible signatures due to charging and discharging of the BH via the Wald
(1974) mechanism (Chen & Dai, 2021; Levin et al., 2018; Dai, 2019; Zhang,
2019; Pan & Yang, 2019; Zhong et al., 2019; Yang et al., 2020b). Interestingly
many of these mechanisms still have total power that scales with the square of
the magnetic field strength, but expected energies of emission will likely
vary. Crustal cracking in the NS may allow for an EM counterpart which is also
strongly dependent on magnetic field strength and could be considered
alongside magnetospheric emission in future work (see Neill et al., 2021, and
references therein). In the future, as models for EM emission in non-
disrupting BHNSs mature, a compendium of signatures at different observing
bands should be generated and used to place a range of model dependent
magnetic field constraints.

##  6 Conclusion

We have presented the first multi-messenger constraints on neutron-star
magnetic field strengths, providing a rare observational handle on the
magnetic fields of neutron stars at these likely late stages of their lives,
before merger with a black hole. The present observations require magnetic
field strengths below $\sim$$10^{15}$ G, in the realm of magnetar field
strengths. This rules out the theoretically possible, high-strength long-lived
or rejuvenated fields at merger, and is consistent with the non-expectation of
mergers of black holes and young magnetars. Future observations with present
technology may enhance these constraints by adding a few similar multi-
messenger BHNS constraints per year, with the expectation of a nearby
($\lesssim 100$ Mpc) event, which could constrain field strengths below
$10^{13-14}$ G, in the next ten years. The limiting factor in achieving better
constraints is the gamma-ray flux sensitivity – more sensitive EM upper-limits
would enhance these constraints with the square root of the flux sensitivity.

Observations with the power to constrain the neutron-star magnetic field at
astrophysically interesting scales warrants further understanding of the BHNS
magnetospheric physics. This will lead to more accurate predictions for the
expected frequency-dependent flux and timing of radiation from these systems,
which will contribute to searches for these signatures and to tighter
constraints from non-detection. It will also hone predictions for multi-
wavelength signatures of these events, which can be leveraged in addition to
multi-messenger constraints as a more powerful probe of the BHNS magnetosphere
and neutron star magnetic field evolution.

We thank Yuri Levin for useful discussions. DJD received funding from the
European Union’s Horizon 2020 research and innovation programme under Marie
Sklodowska-Curie grant agreement No. 101029157. ZH acknowledges support by
NASA grant NNX17AL82G and NSF grants 1715661 and AST-2006176. JL is supported
in part by the Tow Foundation. AVG acknowledges support by the Danish National
Research Foundation (DNRF132). J.S. is supported by the European Union’s
Horizon 2020 research and innovation programme under Marie Sklodowska-Curie
grant agreement No. 844629 and through Villum Fonden grant No. 29466.

## References

  * Abbott et al. (2020) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJL, 896, L44
  * Abbott et al. (2021) —. 2021, ApJL, 915, L5
  * Amaro-Seoane & et al. (2017) Amaro-Seoane, P., & et al. 2017, ArXiv e-prints, arXiv:1702.00786
  * Antoniadis et al. (2021) Antoniadis, J., Aguilera-Dena, D. R., Vigna-Gómez, A., et al. 2021, arXiv e-prints, arXiv:2110.01393
  * Band et al. (1993) Band, D., Matteson, J., Ford, L., et al. 1993, ApJ, 413, 281
  * Bildsten & Cutler (1992) Bildsten, L., & Cutler, C. 1992, ApJ, 400, 175
  * Blackburn et al. (2015) Blackburn, L., Briggs, M. S., Camp, J., et al. 2015, ApJS, 217, 8
  * Bransgrove et al. (2018) Bransgrove, A., Levin, Y., & Beloborodov, A. 2018, MNRAS, 473, 2771
  * Carrasco et al. (2021) Carrasco, F., Shibata, M., & Reula, O. 2021, PRD, 104, 063004
  * Chen & Dai (2021) Chen, K., & Dai, Z. G. 2021, ApJ, 909, 4
  * Dai (2019) Dai, Z. G. 2019, ApJL, 873, L13
  * D’Orazio & Levin (2013) D’Orazio, D. J., & Levin, J. 2013, PRD, 88, 064059
  * D’Orazio et al. (2016) D’Orazio, D. J., Levin, J., Murray, N. W., & Price, L. 2016, PRD, 94, 023001
  * Elfritz et al. (2016) Elfritz, J. G., Pons, J. A., Rea, N., Glampedakis, K., & Viganò, D. 2016, MNRAS, 456, 4461
  * Espinoza et al. (2011) Espinoza, C. M., Lyne, A. G., Kramer, M., Manchester, R. N., & Kaspi, V. M. 2011, ApJL, 741, L13
  * Foucart et al. (2018) Foucart, F., Hinderer, T., & Nissanke, S. 2018, PRD, 98, 081501
  * Fragione (2021) Fragione, G. 2021, arXiv e-prints, arXiv:2110.09604
  * GCN archive for LIGO/Virgo S190814bv (2020) GCN archive for LIGO/Virgo S190814bv. 2020, https://gcn.gsfc.nasa.gov/other/GW190814bv.gcn3, accessed: 2021-10-20
  * GCN archive for S200105ae (2020) GCN archive for S200105ae. 2020, https://gcn.gsfc.nasa.gov/other/S200105ae.gcn3, accessed: 2021-09-21
  * GCN archive for S200115j (2020) GCN archive for S200115j. 2020, https://gcn.gsfc.nasa.gov/other/S200115j.gcn3, accessed: 2021-09-21
  * Geppert et al. (1999) Geppert, U., Page, D., & Zannias, T. 1999, A&A, 345, 847
  * Goldreich & Lynden-Bell (1969) Goldreich, P., & Lynden-Bell, D. 1969, ApJ, 156, 59
  * Goldstein et al. (2016) Goldstein, A., Burns, E., Hamburg, R., et al. 2016, arXiv e-prints, arXiv:1612.02395
  * Gusakov et al. (2020) Gusakov, M. E., Kantor, E. M., & Ofengeim, D. D. 2020, MNRAS, 499, 4561
  * Hansen & Lyutikov (2001) Hansen, B. M. S., & Lyutikov, M. 2001, MNRAS, 322, 695
  * Igoshev et al. (2021) Igoshev, A. P., Hollerbach, R., Wood, T., & Gourgouliatos, K. N. 2021, Nature Astronomy, 5, 145
  * Isoyama et al. (2018) Isoyama, S., Nakano, H., & Nakamura, T. 2018, Progress of Theoretical and Experimental Physics, 2018, 073E01
  * Kaspi (2010) Kaspi, V. M. 2010, Proceedings of the National Academy of Science, 107, 7147
  * Lai (2012) Lai, D. 2012, ApJL, 757, L3
  * Laine & Lin (2012) Laine, R. O., & Lin, D. N. C. 2012, ApJ, 745, 2
  * Lee (1993) Lee, M. H. 1993, ApJ, 418, 147
  * Levin et al. (2018) Levin, J., D’Orazio, D. J., & Garcia-Saenz, S. 2018, PRD, 98, 123002
  * Liu et al. (2020) Liu, S., Hu, Y.-M., Zhang, J.-d., & Mei, J. 2020, PRD, 101, 103027
  * Luo et al. (2016) Luo, J., Chen, L.-S., Duan, H.-Z., et al. 2016, Classical and Quantum Gravity, 33, 035010
  * Lyutikov (2011) Lyutikov, M. 2011, PRD, 83, 064001
  * Lyutikov & McKinney (2011) Lyutikov, M., & McKinney, J. C. 2011, PRD, 84, 084019
  * Mandel & Broekgaarden (2021) Mandel, I., & Broekgaarden, F. S. 2021, arXiv e-prints, arXiv:2107.14239
  * Mandel & Smith (2021) Mandel, I., & Smith, R. J. E. 2021, ApJL, 922, L14
  * McKernan et al. (2020) McKernan, B., Ford, K. E. S., & O’Shaughnessy, R. 2020, MNRAS, 498, 4088
  * McWilliams & Levin (2011) McWilliams, S. T., & Levin, J. 2011, ApJ, 742, 90
  * Mingarelli et al. (2017) Mingarelli, C. M. F., Lazio, T. J. W., Sesana, A., et al. 2017, Nature Astronomy, 1, 886
  * Miralles et al. (2002) Miralles, J. A., Pons, J. A., & Urpin, V. A. 2002, ApJ, 574, 356
  * Narayan et al. (1992) Narayan, R., Paczynski, B., & Piran, T. 1992, ApJL, 395, L83
  * Neill et al. (2021) Neill, D., Tsang, D., van Eerten, H., Ryan, G., & Newton, W. G. 2021, arXiv e-prints, arXiv:2111.03686
  * Paczynski (1986) Paczynski, B. 1986, ApJL, 308, L43
  * Pan & Yang (2019) Pan, Z., & Yang, H. 2019, PRD, 100, 043025
  * Paschalidis et al. (2013) Paschalidis, V., Etienne, Z. B., & Shapiro, S. L. 2013, PRD, 88, 021504
  * Piro (2012) Piro, A. L. 2012, ApJ, 755, 80
  * Sesana (2016) Sesana, A. 2016, Phys. Rev. Lett., 116, 231102
  * Stevenson et al. (2017) Stevenson, S., Vigna-Gómez, A., Mandel, I., et al. 2017, Nature Communications, 8, 14906
  * Tagawa et al. (2021) Tagawa, H., Kocsis, B., Haiman, Z., et al. 2021, ApJ, 908, 194
  * Team COMPAS et al. (2021) Team COMPAS, Riley, J., Agrawal, P., et al. 2021, arXiv e-prints, arXiv:2109.10352
  * Thorne et al. (1986) Thorne, K. S., Price, R. H., & MacDonald, D. A. 1986, Black holes: The membrane paradigm (Yale University Press)
  * Tiengo et al. (2013) Tiengo, A., Esposito, P., Mereghetti, S., et al. 2013, Nature, 500, 312
  * Vietri (1996) Vietri, M. 1996, ApJL, 471, L95
  * Vigna-Gómez et al. (2018) Vigna-Gómez, A., Neijssel, C. J., Stevenson, S., et al. 2018, MNRAS, 481, 4009
  * Vigna-Gómez et al. (2020) Vigna-Gómez, A., MacLeod, M., Neijssel, C. J., et al. 2020, PASA, 37, e038
  * Vigna-Gómez et al. (2021) Vigna-Gómez, A., Aguilera-Dena, D. R., & Willcox, R. 2021, Dataset from: Fallback Supernova Assembly of Heavy Binary Neutron Stars and Light Black Hole-Neutron Star Pairs and the Common Stellar Ancestry of GW190425 and GW200115, doi:10.5281/zenodo.4682798
  * Wald (1974) Wald, R. M. 1974, PRD, 10, 1680
  * Xu et al. (2021) Xu, K., Li, X.-D., Cui, Z., et al. 2021, arXiv e-prints, arXiv:2110.10438
  * Yang et al. (2020a) Yang, Y., Gayathri, V., Bartos, I., et al. 2020a, ApJL, 901, L34
  * Yang et al. (2020b) Yang, Y.-S., Zhong, S.-Q., Zhang, B.-B., et al. 2020b, ApJ, 899, 60
  * Ye et al. (2020) Ye, C. S., Fong, W.-f., Kremer, K., et al. 2020, ApJL, 888, L10
  * Zhang (2019) Zhang, B. 2019, ApJL, 873, L9
  * Zhong et al. (2019) Zhong, S.-Q., Dai, Z.-G., & Deng, C.-M. 2019, ApJL, 883, L19
  * Zrake & MacFadyen (2013) Zrake, J., & MacFadyen, A. I. 2013, ApJL, 769, L29


# Causal Discovery with Missing Data
in a Multicentric Clinical Study

Alessio Zanga  Department of Informatics, Systems and Communication (DISCo),
University of Milano - Bicocca, Milan, Italy  Data Science and Advanced
Analytics, F. Hoffmann - La Roche Ltd, Basel, Switzerland  Alice Bernasconi
Department of Informatics, Systems and Communication (DISCo), University of
Milano - Bicocca, Milan, Italy  Evaluative Epidemiology Unit, Department of
Research, Fondazione IRCCS Istituto Nazionale dei Tumori, Milan, Italy  Peter
J.F. Lucas  University of Twente, Enschede, The Netherlands  Hanny Pijnenborg
RadboudUMC, Nijmegen, The Netherlands  Casper Reijnen  RadboudUMC, Nijmegen,
The Netherlands  Marco Scutari  Istituto Dalle Molle di Studi
sull’Intelligenza Artificiale (IDSIA), Lugano, Switzerland  Fabio Stella
Department of Informatics, Systems and Communication (DISCo), University of
Milano - Bicocca, Milan, Italy

###### Abstract

Causal inference for testing clinical hypotheses from observational data
presents many difficulties because the underlying data-generating model and
the associated causal graph are not usually available. Furthermore,
observational data may contain missing values, which impact the recovery of
the causal graph by causal discovery algorithms: a crucial issue often ignored
in clinical studies. In this work, we use data from a multi-centric study on
endometrial cancer to analyze the impact of different missingness mechanisms
on the recovered causal graph. This is achieved by extending state-of-the-art
causal discovery algorithms to exploit expert knowledge without sacrificing
theoretical soundness. We validate the recovered graph with expert physicians,
showing that our approach finds clinically-relevant solutions. Finally, we
discuss the goodness of fit of our graph and its consistency from a clinical
decision-making perspective using graphical separation to validate causal
pathways.

##  1 Introduction

Much of the data collected in clinical research is observational, collected as
part of daily clinical practice. Correctly interpreting them requires a good
understanding of their characteristics and of possible sources of bias. A
common one is missing values, which may arise in three different ways [3]:
_data missing completely at random_ (MCAR), _data missing at random_ (MAR),
and _data missing not at random_ (MNAR) that are neither MCAR nor MAR. MNAR is
common in clinical observational data and thus interesting to study, as it is
often possible to unravel the reason for the missingness: for instance, a
laboratory test may be skipped in favour of a more precise ones available at a
later stage.

Missing values in clinical data are commonly imputed with heuristics or with
single/multiple imputation. Such techniques assume that the data are MCAR or
MAR; we cannot test whether these assumptions are valid without knowing the
missingness mechanism but, at the same time, if these assumptions do not hold
our clinical conclusions are likely to be biased [4]. A possible approach to
this problem is _causal discovery_ : modelling the missingness mechanism is to
recover the underlying causal graph $\mathcal{G}^{*}$, given the data
$\mathcal{D}$ and the prior knowledge $\mathcal{K}$ [6]. In our previous work
on endometrial cancer (EC) [5], we proposed a new causal discovery approach
based on bootstrapping for clinical data with low sample size and high
missingness assuming MAR (_Bootstrap SEM_). Algorithms assuming MNAR were not
available until recently when _HC-aIPW_ [1] was introduced.

Our aim is to showcase how modern causal discovery techniques can model the
biases in observational data, in particular for MNAR. For this purpose, we
applied different causal discovery algorithms with different assumptions to
data from a multicenter study on EC, highlighting the clinical implications of
their biases on recovering the causal mechanisms behind the prognosis of EC.

##  2 Background

A causal graph $\mathcal{G}=(\mathbf{V},\mathbf{E})$ is a directed acyclic
graph (DAG) where for each directed edge $(X,Y)\in\mathbf{E}$, $X$ is a direct
cause of $Y$ and $Y$ is a direct effect of $X$. The vertex set $\mathbf{V}$ is
usually split into two disjoint subsets $\mathbf{V}=\mathbf{O}\cup\mathbf{U}$,
where $\mathbf{O}$ is the set of the _fully observed_ variables (with no
missing values), while $\mathbf{U}$ is the set of _fully unobserved_ variables
(the _latent_ variables). A missingness graph
$\mathcal{M}=(\mathbf{V}^{*},\mathbf{E}^{*})$ [2] is a causal graph where the
vertices in $\mathbf{V}^{*}$ are partitioned into five disjoint subsets:
$\mathbf{V}^{*}=\mathbf{O}\cup\mathbf{U}\cup\mathbf{M}\cup\mathbf{S}\cup\mathbf{R}$,
where $\mathbf{M}$ is the set of the _partially observed_ variables, that is,
the variables with at least one missing value; $\mathbf{S}$ is the set of the
proxy variables, that is, the variables that are actually observed;
$\mathbf{R}$ is the set of the _missingness indicators_. Missingness graphs
can be queried for independencies using _d-separation_. The set of variables
$\mathbf{Z}$ d-separates $X$ from $Y$, denoted by $X\perp\\!\\!\\!\perp
Y\mid\mathbf{Z}$, if it _blocks_ every path $\pi$ between $X$ and $Y$. A path
$\pi$ is blocked by $\mathbf{Z}$ if and only if $\pi$ contains: a fork
$A\leftarrow B\rightarrow C$ or a chain $A\rightarrow B\rightarrow C$ so that
$B$ is in $\mathbf{Z}$, or, a collider $A\rightarrow B\leftarrow C$ so that
$B$, or any descendant of it, is not in $\mathbf{Z}$. MCAR, MAR and MNAR
result in different independence statements [3] which are linked to the
independency statements implied by the missingness graph: MCAR implies
$\mathbf{O}\cup\mathbf{U}\cup\mathbf{M}\perp\\!\\!\\!\perp\mathbf{R}$, the
missingness is random and independent from the fully observed and the
partially observed variables; MAR implies
$\mathbf{U}\cup\mathbf{M}\perp\\!\\!\\!\perp\mathbf{R}\mid\mathbf{O}$,
missingness is random only conditionally on the fully observed variables; MNAR
if neither MCAR nor MAR. Since MCAR implies MAR, any method assuming MAR can
be used on MCAR.

##  3 Multicentric Clinical Data on Endometrial Cancer

The observational data we explore in this paper comprise 763 patients with
endometrial cancer from 10 gynecological oncological clinics in Europe that
are part of the European Network for Individualized Treatment of Endometrial
Cancer (ENITEC). Clinical experts selected the variables that they considered
most important for predicting the presence of lymph node metastases (LNMs) and
survival [5]. The selected variables were: the cytology of the cervix uteri,
the preoperative tumour grade, the postoperative tumour grade (after
pathological examination of the tumour tissue obtained after surgical removal
of the uterus), treatment by chemotherapy or radiatiotherapy, lymphvascular
space invasion (that is, whether there is tumour growth into the lymph or
blood vessels), the levels of estrogen and progesterone in blood, the presence
of lymph node metastasis according to CT or MRI imaging, the CA125 tumour
marker, L1CAM (an intracellular protein that promotes tumour cell motility),
the p53 tumour suppressor gene, the number of platelets, presence of lymph
node metastases, recurrence of the tumour, and lastly survival before and
after 1, 3, and 5 years. The tumour markers (p53, CA125, L1CAM, estrogen and
progesterone levels) are thought to offer causal prognostic information about
tumour cell behaviour and thus tumour ingrowth, metastases, recurrence, and
survival.

##  4 Experiments

We performed numerical experiments to compare the graphs recovered under MAR
and MNAR by the Bootstrap SEM and HC-aIPW. For reference, we also reported the
results for HC on data completed with single imputation, denoted _HC-
complete_. Prior knowledge elicited from experts consists of forbidden and
required edges (Survival$1$yr $\to$ Survival$3$yr, Survival$3$yr $\to$
Survival$5$yr).

Firstly, we evaluated the goodness of fit of the recovered graphs by computing
the log-likelihood (LL) of both the data used to recover the graph (in-sample)
and those held aside for validation (out-of-sample). The former allows us to
see which algorithm fits a particular data set the best. The latter
approximates the Kullback-Leibler divergence between the recovered causal
graph and the unknown causal graph underlying the data, and allows us to see
how close the two are and how well the recovered graph generalises to new
data. We repeated causal discovery for 100 bootstrap replicates and computed
the mean and the standard deviation of both in-sample and out-of-sample LL.

Figure 1: In-sample and out-of-sample LL for each algorithm. Re-scaled by
sample size and absolute maximum value. Lower values are better.

We observe (Figure 1) that HC-aIPW, which assumes MNAR, dominates Bootstrap
SEM and HC-complete, which assume MAR, for both in-sample and out-of-sample
LL. In the case of in-sample LL, this may be attributed (at least in part) to
making the correct missingness assumption: MCAR and MAR are too restrictive
and limit how well the recovered graph fits the data; MNAR it is not strict
enough when MCAR or MAR hold and let causal discovery algorithms overfit. This
decreases the out-of-sample LL because overfitted models are too complex and
do not generalize well. The fact that HC-aIPW, which assumes MNAR, outperforms
Bootstrap SEM and HC-complete, which assume MAR, for both in-sample and out-
of-sample LL suggests that the MNAR assumption is correct for the data and
that it allows HC-aIPW to recover a causal graph that is close to the
underlying model and that generalises better to new data.

Goodness-of-fit measures allow for a _quantitative_ evaluation of the
recovered graphs, but they say little about the _qualitative_ information they
encode in terms of _independence statements_. We denote the graph recovered by
Bootstrap SEM as $\mathcal{M}_{\textrm{MAR}}$; Figure 2 shows that recovered
by HC-aIPW, $\mathcal{M}_{\textrm{MNAR}}$. For readability, we colored the
vertices depending on their _semantic_ interpretation: treatments
(Radiotherapy and Chemotherapy) are colored in blue, outcomes (Survival1yr,
etc.) in red, the event of interest (LNM) in orange, relevant biomarkers (ER,
PR, CA125, etc.) in lightblue and the _context_ variable Hospital in gray.

Focusing on the interactions of LNM, we observe using d-separation that
$\textrm{LNM}\perp\\!\\!\\!\perp\\{\textrm{CA125},\textrm{p53}\\}\mid\textrm{PostoperativeGrade}$
is true in $\mathcal{M}_{\textrm{MAR}}$, but false in
$\mathcal{M}_{\textrm{MNAR}}$, where CA125 and p53 are effects of LNM. This
makes $\mathcal{M}_{\textrm{MNAR}}$ close to the clinical practice where both
CA125 and p53 are considered relevant biomarkers linked to LNM, providing
additional information on LNM even if PostoperativeGrade is observed. Indeed,
a crucial difference between $\mathcal{M}_{\textrm{MAR}}$ and
$\mathcal{M}_{\textrm{MNAR}}$ is that the biomarkers and LNM are d-separated
from LNM in $\mathcal{M}_{\textrm{MAR}}$ if PostoperativeGrade is observed,
but are descendants of LNM in $\mathcal{M}_{\textrm{MNAR}}$. Hence, if our
goal is to detect the presence of LNM in EC patients, measuring CA125, p53 or
any of their descendants is coherent with $\mathcal{M}_{\textrm{MNAR}}$.

Figure 2: Causal graph $\mathcal{M}_{\textrm{MNAR}}$ recovered by HC-aIPW
under MNAR.

Shifting the focus to the treatment variables Chemotherapy and Radiotherapy,
$\textrm{LNM}\perp\\!\\!\\!\perp\textrm{Chemotherapy}$ does not hold in either
$\mathcal{M}_{\textrm{MAR}}$ or $\mathcal{M}_{\textrm{MNAR}}$, since
Chemotherapy is a direct cause of LNM in both. Therefore, Chemotherapy is
expected to influence the likelihood of LNM, which is exactly the reason why
it is prescribed by clinicians. On the other hand,
$\textrm{LNM}\perp\\!\\!\\!\perp\textrm{Radiotherapy}$ does hold in
$\mathcal{M}_{\textrm{MNAR}}$ but not in $\mathcal{M}_{\textrm{MAR}}$,
suggesting a spurious correlation induced by MAR. This is confirmed by the
clinical literature: since Radiotherapy is aimed at local treatment of the
tissue surrounding the uterus, and there is a clear dependence with
MyometrialInvasion of the tumour (in both models), Radiotherapy effects on LNM
are not expected.

##  5 Discussion & Conclusions

In this work we presented a systematic analysis of the impact of missingness
assumptions using state-of-the-art causal discovery algorithms. We applied
these methods to a real-world, observational multicentric study on EC
patients, extending them to include expert prior knowledge without sacrificing
theoretical soundness. Furthermore, we validated the obtained causal models
with experienced physicians and clinical literature. We evaluated the
goodness-of-fit of the recovered graphs with respect to the underlying data
distribution, showing that stricter assumptions are associated to models that
generalize poorly. Moreover, by leveraging the test for graphical separation,
we explained how the missingness mechanism affects the causal pathways
associated to the clinical decision-making perspective. Quantifying the bias
due to missingness in other case studies and its overlap with the effects of
hidden and selection variables are open problems.

### Acknowledgements

Alessio Zanga is funded by F. Hoffmann-La Roche Ltd.

## References

  * [1] Yang Liu and Anthony C. Constantinou.  Greedy Structure Learning from Data that Contain Systematic Missing Values.  Machine Learning, 111(10):3867–3896, 10 2022.
  * [2] Karthika Mohan and Judea Pearl.  Graphical Models for Processing Missing Data.  Journal of the American Statistical Association, 116(534):1023–1037, 1 2018.
  * [3] Donald B. Rubin.  Inference and Missing Data.  Biometrika, 63(3):581–592, 1976.
  * [4] Marianne Riksheim Stavseth, Thomas Clausen, and Jo Røislien.  How Handling Missing Data may Impact Conclusions: A Comparison of Six Different Imputation Methods for Categorical Questionnaire Data.  SAGE Open Medicine, 7:205, 1 2019.
  * [5] Alessio Zanga, Alice Bernasconi, Peter J F Lucas, Johanna M. A. Pijnenborg, Casper Reijnen, Marco Scutari, and Fabio Stella.  Risk Assessment of Lymph Node Metastases in Endometrial Cancer Patients: A Causal Approach.  In Proceedings of the 1st Workshop on Artificial Intelligence For Healthcare, 2022.
  * [6] Alessio Zanga, Elif Ozkirimli, and Fabio Stella.  A Survey on Causal Discovery: Theory and Practice.  International Journal of Approximate Reasoning, 151:101–129, 12 2022.