# Operator lifetime and the force-free electrodynamic limit
of magnetised holographic plasma

Napat Poovuttikul<EMAIL_ADDRESS>Department of Mathematical
Sciences, Durham University, South Road, Durham DH1 3LE, UK  University of
Iceland, Science Institute, Dunhaga 3, IS-107, Reykjavik, Iceland  Aruna
Rajagopal<EMAIL_ADDRESS>University of Iceland, Science Institute, Dunhaga 3,
IS-107, Reykjavik, Iceland

###### Abstract

Using the framework of higher-form global symmetries, we examine the regime of
validity of force-free electrodynamics by evaluating the lifetime of the
electric field operator, which is non-conserved due to screening effects. We
focus on a holographic model which has the same global symmetry as that of low
energy plasma and obtain the lifetime of (non-conserved) electric flux in a
strong magnetic field regime. The lifetime is inversely correlated to the
magnetic field strength and thus suppressed in the strong field regime.

###### Contents

  1. I Introduction
  2. II The Holographic Model
    1. II.1 Linearised solutions in $\omega/T\ll 1$ limit and matching procedure
      1. II.1.1 Perturbation parallel to equilibrium magnetic field
      2. II.1.2 Perturbation perpendicular to equilibrium magnetic field
    2. II.2 Checking $T\gtrsim 0$ limit in $\omega/\sqrt{|\bf{B|}}\ll 1$ regime
      1. II.2.1 Zero temperature
      2. II.2.2 $T\lesssim\omega\ll\sqrt{|{\bf B}|}$ limit
  3. III Conclusion
  4. A Numerical solution and evaluation of operators lifetime
  5. B Frobenius analysis in $AdS_{3}\times\mathbb{R}^{2}$ region

##  I Introduction

Hydrodynamics LLfluid  is a well-established theoretical framework which
universally describes the long wavelength, low frequency behaviour of
interacting systems at finite temperature. Essentially, hydrodynamic theory is
a description of conserved quantities and the manifestation of the
corresponding symmetries in a system in thermal equilibrium. Theories with
widely varying microscopics can have the same macroscopic hydrodynamic
description. One possible explanation why such a universal description is
possible is that all operators except conserved charges have parametrically
short lifetimes compared to the scale of interest and, once the longest-lived
non-conserved operator111While this operator language is more familiar in the
context of quantum systems, it is also applicable to classical systems via
e.g. memory matrix formalism zwanzig1995 ; forster1995 . A more modern
introduction may be found in Hartnoll:2012rj . has decayed away, the
hydrodynamic description becomes viable (see Fig. 1).

Figure 1:  A cartoon illustration of the lifetime of operators of a theory
that exhibit hydrodynamic behaviour at late time. Here, there is a
parametrically large gap between conserved charges $\rho$ and the rest. The
life time $\tau_{1}$ of the longest-lived operator, denoted by
$\mathcal{O}_{1}$ set the time scale in which hydrodynamics becomes
applicable.

The hydrodynamic framework may be generalised to systems where the conserved
charges are those of a higher-form symmetry Gaiotto:2014kfa  which counts the
number density of extended objects. A recent exploration of this idea
Grozdanov_2017  (see also Schubring:2014iwa ; Hernandez:2017mch ;
Armas:2018atq ; Armas:2018zbe ) shows that the resulting hydrodynamics of a
one-form $U(1)$ charge reproduces the theory of magnetohydrodynamics
(MHD)222For the formulation of MHD that closely resembles higher-form symmetry
formulation, see e.g. dixon1982special ; anile2005relativistic ;
komissarov1999 .. This should not come as a big surprise. MHD is, after all, a
low energy effective theory of plasma where the (dynamical) electric field is
screened – the one-form U(1) symmetry associated to electric flux is
explicitly broken. This implies that in, for example, a plasma at zero
magnetic field (where the Ohm’s law ${\bf j}=\sigma{\bf E}$ is a good
approximation) the electric field has a finite lifetime,

$\delta{\bf E}\propto\exp(-t/\tau_{E})\qquad\Longleftrightarrow\qquad\langle
E^{i}(-\omega)E^{j}(\omega)\rangle\sim\frac{\delta^{ij}}{\omega+i/\tau_{E}}\,,$
(1)

with the lifetime of the electric field $\tau_{E}=1/\sigma$. The conductivity
$\sigma$ can be computed from first principles. For instance, in quantum
electrodynamics, it can be written as Arnold_2000 333The fact that this
quantity has only been computed at the beginning of this century indicates the
difficulty of the required computations.

$\sigma\propto\frac{T}{e^{2}\log e^{-1}}\,,$ (2)

where $e$ is the electromagnetic coupling. The lifetime of electric field
$\tau\sim 1/T$, is then much shorter than the scale $t\gg 1/T$ (or
$\omega/T\ll 1$ in Fourier space) where hydrodynamic behaviour is expected.
If, in this late-time limit, all other operators except energy density
$T^{tt}$ and the momentum $T^{ti}$ have already decayed away, one can expect
the hydrodynamic description of a plasma to be governed by

$\partial_{\mu}T^{\mu\nu}=0\,,\qquad\partial_{\mu}J^{\mu\nu}=0\,.$ (3)

The conserved currents $T^{\mu\nu}$ and $J^{\mu\nu}$ are expressed in terms of
energy, momentum, magnetic flux $J^{ti}\equiv B^{i}$ and their conjugates,
organised order by order in the gradient expansion. This formulation of MHD
only requires macroscopic consistency and does not require the introduction of
the gauge field $\star J=F=dA$ which, due to screening effect, is not a long-
lived degree of freedom.

This brings us to the central question of the present paper: Is a hydrodynamic
description of the form (3) applicable in the limit low temperature compared
to magnetic flux density $T^{2}/|{\bf B}|\ll 1$ ? This question is important
if one wants to apply the MHD description to astrophysical plasmas where the
magnetic field is many orders of magnitude larger than the scale set by the
temperature. %ูA more precise way to phrase this question is whether or not
non-conserved operators, decayed away at the following time scale If one were
to naively extrapolate (1)-(2), the lifetime of the electric field appears to
become arbitrarily long as the temperature decreases. However, there exists a
macroscopic description of plasma in this regime that has been successfully
applied. This theory is called force-free electrodynamics or FFE, and has been
used extensively in astrophysical setups such as in the magnetosphere of black
holes Blandford:1977ds ; Komissarov_2004 , neutron star 1969ApJ…157..869G  and
solar corona Wiegelmann_2012  just to name a few. In its conventional form,
this theory is applied to a system which is magnetically dominated (i.e.
$|{\bf B}^{2}|>|{\bf E}^{2}|$ or, covariantly $F_{\mu\nu}F^{\mu\nu}>0$) and
whose dynamics is governed by

$\displaystyle\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}$
$\displaystyle=0\,,$ (4a) $\displaystyle
F^{\mu\nu}\nabla_{\lambda}F^{\lambda}_{\;\;\nu}$ $\displaystyle=0\,.$ (4b)

Here, the first relation implies that ${\bf E}\cdot{\bf B}=0$ while the second
relation implies that the force $j_{el}^{\mu}F_{\mu\nu}$, with
$j^{\mu}_{el}:=\nabla_{\nu}F^{\nu\mu}$ via Maxwell’s equations, acting on
plasma vanishes (hence the name force-free electrodynamics). More details on
the geometric and effective action view point of FFE can be found in e.g.
Gralla_2014 ; Compere:2016xwa  and Uchida:1997 ; Thompson:1998ss ; Gralla_2019
; glorioso2018effective . One should emphasise that the system of equations in
Eq.(4) is independent of the microscopic details of the cold plasma, which
then strongly resembles hydrodynamic descriptions. In fact, it turns out that
(4) can arise in a special limit where $T\ll\sqrt{|{\bf B}|}$ of a
hydrodynamics description with one-form $U(1)$ symmetry in (3), see
Grozdanov_2017 ; Gralla_2019 ; glorioso2018effective 444Recasting of force-
free electrodynamics in the hydrodynamic language also allows the systematic
gradient expansions Gralla_2019 ; Benenowski:2019ule . This could serve to
classify correction to FFE in order to account for phenomena such as pulsar
radio emission where ${\bf E}\cdot{\bf B}\neq 0$. .

The existence of FFE is usually justified by saying that the cold plasma is,
on one hand, dense enough to screen the electric field (4a) but, on the other
hand, dilute enough so that force-free condition (4b) is applicable. This
statement can be made more precise in the light of relations between the
equations of FFE and hydrodynamics. Thus, we propose a criterion for testing
the validity of FFE using the lifetime of non-conserved operators – FFE, or
equivalently, hydrodynamic description of cold plasma in the $T\ll\sqrt{|{\bf
B}|}$ limit, is valid when the lifetime of all non-conserved operators is
parametrically shorter than the time scale of interest. A key advantage of
this approach is that the operator lifetime can be, in principle, computed
explicitly from microscopic description and therefore allows one to find the
‘cutoff’ scale where FFE description should break down.

Computing the operator lifetime from microscopic description is, however, not
always an easy task. In fact, we are not aware of a genuine computation
directly from quantum electrodynamics (in the sense of Arnold_2000 ) when both
$T$ and ${\bf B}$ are turned on. To simplify the computations, we shall
demonstrate the validity of FFE in the strongly interacting magnetised plasma
with a holographic dual as proposed in Grozdanov_2019b ; Hofman_2018  where
the one-form $U(1)$ global symmetry is taken into account via a two-form gauge
field in the gravity dual. This provides two key advantages. First, the
computation of correlation functions boils down to solving simple linearised
differential equations (see e.g. Son:2002sd ). Second, there is strong
evidences that charge neutral operators, apart from energy and momentum, have
a parametrically short lifetime in this class of theories 555To be more
precise, it has been shown in $\mathcal{N}=4$ supersymmetric Yang-Mills
theory, which constitutes the matter sector of the holographic model
Grozdanov_2019b ; Hofman_2018 , that there is no long-lived mode besides
hydrodynamic modes at any $T\neq 0$ and $|{\bf B}|=0$ Kovtun:2003wp . A
similar conclusion was reached for the same theory in the charge neutral
sector at finite non-dynamical magnetic field Fuini:2015hba ;
Janiszewski:2015ura . . Therefore, we shall focus on non-conserved operators
in the electromagnetic sector of the theory: the electric flux operators,
whose lifetime can be extracted via two-point correlation function as in (1).
This will provide strong evidence for the validity of FFE limit in a strongly
interacting holographic plasma.

On the technical side, the computations presented in this note show that there
are no quasinormal modes present in the vicinity of the hydrodynamic regime
$\omega/T\ll 1$ (and $\omega/\sqrt{|{\bf B}|}\ll 1$). The pole in the electric
flux correlation function in this regime then implies that the operator has a
parametrically long lifetime which could interfere with the hydrodynamic
modes. The presence of such long-lived mode can be determined analytically in
the usual hydrodynamic regime of $\omega/T\ll 1$ for a large class of
theories. It is usually difficult to go beyond this regime towards the limit
$\omega/T\sim 1,\omega/\sqrt{|{\bf B}|}\ll 1$. Such computation can, however,
be done analytically in the simple model of Grozdanov_2019b  thanks to the
presence of the BTZ$\times\mathbb{R}^{2}$ bulk geometry in the deep IR
DHoker:2009mmn . We should also note that the treatment of a (long-lived) non-
hydrodynamic modes has been extensively used to determine the breakdown of
hydrodynamic descriptions in the context of QFTs with holographic duals, see
e.g. Grozdanov_2019 ; Davison:2014lua ; Chen:2017dsy ; Davison:2018nxm  .

The remainder of this paper is organised as the follows. In section II, we
summarise the procedure involved in the computation of the two-point
correlation function in the holographic dual to one-form global symmetry. In
section II.1, we outline the method for exploring the existence of decaying
modes in the vicinity of the usual hydrodynamic limit $\omega/T\ll 1$ at
$T/\sqrt{{\bf B}|}\ll 1$. Due to the simplicity of the bulk geometry, we are
able to further extend the analysis to arbitrary value of $\omega/T$ with
$\omega/\sqrt{|{\bf B}|}\ll 1$ and $T/\sqrt{|{\bf B}|}\ll 1$. This is
described in section II.2. Further open questions and future directions are
discussed in section III.

##  II The Holographic Model

A simple holographic dual to a strongly interacting field theory of matter
charged under dynamical $U(1)$ electromagnetism (that is, the dynamical plasma
described by low energy MHD) and formulated in the language of higher-form
symmetry was constructed in Grozdanov_2019b ; Hofman_2018 . We present a brief
review here for completeness. The five-dimensional bulk theory is comprised of
Einstein gravity coupled to a two-form bulk gauge field, $B_{\mu\nu}$, and a
negative cosmological constant,

$S=\int
d^{5}X\sqrt{-G}\left(R-2\Lambda-\frac{L^{2}}{3}H_{abc}H^{abc}\right)+S_{bnd}-\frac{1}{\kappa(\Lambda)}\int_{r=\Lambda}d^{4}x\sqrt{-\gamma}(n^{a}H_{a\mu\nu})(n_{b}H^{b\mu\nu}),$
(5)

where $H=dB$ and $B_{ab}$ is the bulk 2-form gauge field, $\Lambda$ is the UV-
cutoff, $n^{a}$ is the unit normal to the boundary, and $S_{bnd}$ denotes the
Gibbons-Hawking and gravitational counter term. Roughly speaking, the two bulk
fields $G_{ab}$ and $B_{ab}$, asymptote to $g_{\mu\nu}$ and $b_{\mu\nu}$
respectively, which then source the currents, $T^{\mu\nu}$ and $J^{\mu\nu}$.

$\langle T_{\mu\nu}\rangle\equiv\frac{2}{\sqrt{-g}}\frac{\delta S}{\delta
g_{\mu\nu}}\,,\qquad\langle
J_{\mu\nu}\rangle\equiv\frac{1}{\sqrt{-g}}\frac{\delta S}{\delta b_{\mu\nu}}$
(6)

The generating functional takes the form,

$Z[g_{\mu\nu},b_{\mu\nu}]=\left\langle\exp\left[i\int
d^{4}x\sqrt{-g}\left(\frac{1}{2}T^{\mu\nu}g_{\mu\nu}+J^{\mu\nu}b_{\mu\nu}\right)\right]\right\rangle$
(7)

and diffeomorphism invariance and gauge symmetry lead to the following
equations,

$\nabla_{\mu}\langle T^{\mu\nu}\rangle=(db)^{\nu}_{\rho\sigma}\langle
J^{\rho\sigma}\rangle\,,\qquad\nabla_{\mu}\langle J^{\mu\nu}\rangle=0.$ (8)

$H=db$ is the three-form field strength of the two-form external source. The
equilibrium solution of this holographic model is a domain wall interpolating
between an asymptotic $AdS_{5}$ geometry in the UV ($r\to\infty$ in our
convention), and $BTZ\times\mathbb{R}^{2}$ in the near-horizon IR ($r=r_{h}$).
It is described by the following metric and gauge field

$\displaystyle ds^{2}$
$\displaystyle=G_{ab}dX^{a}dX^{b}=-r^{2}f(r)dt^{2}+\frac{dr^{2}}{r^{2}f(r)}+e^{2V(r)}(dx^{2}+dy^{2})+e^{2W(r)}dz^{2}\,,$
(9) $\displaystyle B$ $\displaystyle=h(r)dt\wedge
dz\qquad\text{with}\qquad\star_{5}H=\mathcal{B}dx\wedge dy$

Modulo the subtleties due to the mixed boundary conditions, this is nothing
but the hodge dual of the magnetised black brane solution of DHoker:2009mmn .
The radial coordinate is chosen such that $r\to\infty$ corresponds to the
usual asymptotic $AdS_{5}$ with

$f(r)=1\,,\qquad e^{2V(r)}=e^{2W(r)}=r^{2}$ (10)

in the $r\to\infty$ limit. The $BTZ\times\mathbb{R}^{2}$ solution near the
horizon can be written as

$f(r)=3\left(1-\frac{r_{h}^{2}}{r^{2}}\right)\,,\qquad
e^{2V}=\frac{\mathcal{B}}{\sqrt{3}}\,,\qquad e^{2W}=3r^{2}\,.$ (11)

The temperature is set by the horizon radius via $4\pi
T=r_{h}^{2}|f^{\prime}(r_{h})|=6r_{h}/L^{2}$. We set $L=1$ for simplicity.
Note also that $\mathcal{B}$ is related to the $z-$component of the ‘physical’
magnetic field ${\bf B}$ which differs by a prefactor $L$ or the 2-form gauge
field coupling in the bulk (e.g. if one were to define the action with
$S\sim\int(1/g^{2})H^{2}$). We will keep using $\mathcal{B}$ to emphasise its
holographic origin but there is no harm in thinking of it as simply ${\bf B}$.

One interesting feature of this model is that the leading divergence of
$B_{\mu\nu}$ in the Fefferman-Graham expansion is logarithmic. Thus, the
definition of the source $b_{\mu\nu}$ requires mixed boundary condition

$b_{\mu\nu}=B_{\mu\nu}(\Lambda)-\frac{1}{\kappa(\Lambda)}\langle
J_{\mu\nu}\rangle\,,\qquad\text{with}\qquad\langle
J^{\mu\nu}\rangle=-\sqrt{-G}n_{\alpha}H^{\alpha\mu\nu}$ (12)

Requiring the source $b_{\mu\nu}$ to be independent of the UV cutoff fixes the
form of the ‘coupling constant’ $1/\kappa(\Lambda)$ which turns out to be
logarithmically running. This is a common feature for fields with this type of
near-boundary behaviour where the counterterm also plays the role of the
double-trace deformation Witten:2001ua ; Berkooz:2002ug , see also Hofman_2018
; Grozdanov_2019b  for a discussion in the present context. Mapping
$J^{\mu\nu}$ in to a more familiar dynamical field strength via
$J^{\mu\nu}=\frac{1}{2}\mathcal{\epsilon}^{\mu\nu\rho\sigma}F_{\rho\sigma}$,
one can see that the double-trace deformation plays a role similar to the
Maxwell term for the dynamical gauge field in the dual QFT with
$1/\kappa(\Lambda)$ as a (logarithmically running) electromagnetic coupling.

The finite part of $1/\kappa(\Lambda)$ plays a crucial role in this setup.
While the finite counterterm in the ordinary bulk Maxwell theory simply
results in a contact term in the correlation function, the mixed boundary
condition for $B_{ab}$ implies the existence of the purely decaying mode
$\omega=-i/\tau_{E}$ that can interfere with the gapless hydrodynamic
excitation. This is nothing but the life-time of the electric flux operator
$Q_{E}\sim\int dS_{ij}J^{ij}$ which appears in the following correlation
function Hofman_2018 ; Grozdanov_2019

$\langle J^{ij}(t)J^{kl}(0)\rangle\sim\exp\left(-it/\tau_{E}\right)\,.$ (13)

Note that, due to the anisotropy introduced by finite equilibrium magnetic
field, the value of $\tau_{E}$ depends on which direction of the electric
field in consideration. The limit where $\tau_{E}$ is small, but finite,
compared to the length scale of interest (set by temperature or magnetic flux
density) is of particular interest as it allows one to extract $\tau_{E}$
analytically, via a matching procedure that we outline below. As argued in the
introduction, the lifetime of the electric flux determines the validity of MHD
and FFE description.

###  II.1 Linearised solutions in $\omega/T\ll 1$ limit and matching procedure

In this section, we outline the computation required to obtain the relaxation
time of the electric field. We focus on the hydrodynamic regime where
$\omega/T\ll 1$, and the low temperature limit 666Similar computation for the
holographic theory dual to a system with ordinary(zero-form) $U(1)$ symmetry
can be found in e.g. Davison:2013bxa ; Moitra:2020dal .  $T/\sqrt{|{\bf
B}|}\ll 1$. This allows us to solve the bulk equation of motion analytically
via a matching method similar to that was employed in Kovtun:2003wp  (see also
Grozdanov_2019  for a recent review). We consider the decay rate of the
electric field both along and perpendicular to the equilibrium magnetic field
denoted by $E^{\parallel}=J^{xy}$ and $E^{\perp}=J^{xz},J^{yz}$ respectively.

Before proceeding, let us summarise the matching procedure for the
$\omega/T\ll 1$ expansion. It involves separating the bulk into three suitably
defined pieces: inner region, intermediate region and outer region. The inner
region is a suitably defined region close to the horizon while the outer
region is defined to be the range of $r$ such that $\omega/r\ll 1$ so that one
can drop terms quadratic in $(\omega/r)^{2}$, which includes the near boundary
region. The integration constants of the solution in the outer region are
determined by matching the form of inner region solution for intermediate
value of $r$ that connect the two regions together. In our case, this is the
region of $r$ close to $r_{h}$ but

$\frac{\omega}{T}\log f(r)\ll 1$ (14)

This intermediate region defined above is also consistent with the outer
region assumption where $\omega/r\ll 1$ and thus we are able to match the two
solution together. Note that, while this procedure is applicable to any bulk
solution with event horizon, the limit $\omega/T\ll 1$ is crucial.

We now present the key equations and resulting lifetime of the electric flux.

####  II.1.1 Perturbation parallel to equilibrium magnetic field

As the magnetic field in equilibrium points along the $z-$direction, we are
interested in $E^{\parallel}=\frac{1}{2}\varepsilon^{zxy}\langle
J_{xy}\rangle$. The corresponding bulk perturbation is $\delta B_{xy}$ which
decouples from the metric perturbation in the zero wave vector limit. The bulk
equation of motion can be written as

$\left(r^{2}fe^{W-2V}\delta
B_{xy}^{\prime}\right)^{\prime}+\frac{\omega^{2}}{r^{2}f}e^{W-2V}\delta
B_{xy}=0$ (15)

where $(...)^{\prime}$ denotes a derivative w.r.t. the radial coordinate $r$.
The inner region solution for $\delta B_{xy}$, where we substitute the
$BTZ\times\mathbb{R}^{2}$ solution for $f,V,W$, with the ingoing boundary
condition can be written as

$\delta B_{xy}^{inner}=c^{H}\exp\left(-\frac{i\omega}{4\pi T}\log f(r)\right)$
(16)

The outer region solution can be obtained by considering the solution at
linear order in $\omega/r$ and one obtains,

$\displaystyle\delta B_{xy}^{outer}(r)$
$\displaystyle=c_{1}-c_{2}\left(\log\Lambda-\int^{\Lambda}_{\mathtt{r}=r}d\mathtt{r}\,\frac{e^{2V(\mathtt{r})-W(\mathtt{r})}}{\mathtt{r}^{2}f(\mathtt{r})}\right)$
(17) $\displaystyle=c_{1}-c_{2}\left(\log
r-\phi(r)+\frac{e^{2V-W}}{r_{h}^{2}f^{\prime}}\Big{|}_{r=r_{h}}\log
f\right)\,,$

where $\phi(r)$ is a function regular everywhere in the bulk defined as

$\phi(r)=\int^{\Lambda}_{\mathtt{r}=r}d\mathtt{r}\,\left[\frac{e^{2V(\mathtt{r})-W(\mathtt{r})}}{\mathtt{r}^{2}f(\mathtt{r})}-\left(\frac{e^{2V(\mathtt{r})-W(\mathtt{r})}}{r_{h}^{2}f^{\prime}(\mathtt{r})}\right)_{\mathtt{r}=r_{h}}\frac{f^{\prime}(\mathtt{r})}{f(\mathtt{r})}-\frac{1}{\mathtt{r}}\right]\,.$

This parametrisation allows us to single out leading contributions that
dominate when considering the solution near $r=\Lambda$, where $\phi(r)$ and
$\log(e^{-2V}r^{2}f)$ vanish, as well as near $r\approx r_{h}$ where the $\log
f$ term dominates. The integration constants $c_{1},c_{2}$ in (17) are related
to the source $b_{\mu\nu}$ and the 2-form current $\langle J^{xy}\rangle$. The
precise relations can be obtained via Eq. (12) to be

$\langle J^{xy}\rangle=c_{2}\,,\qquad
b_{xy}=c_{1}-\left(\log\,\Lambda+\frac{1}{\kappa(\Lambda)}\right)c_{2}$ (18)

Note that, for the source to be independent of the UV cutoff, one requires
$\kappa(\Lambda)^{-1}=\text{finite term}-\log\Lambda$. This is the
logarithmically running coupling usually found in a double-trace deformed
theory and resembles the running of electromagnetic coupling as pointed out in
Grozdanov_2019b ; Hofman_2018 ; Grozdanov_2019 .

For the outer and inner region solutions to match, we consider both solutions
in the intermediate region where we can write the inner solution as

$\exp\left(-\frac{i\omega}{4\pi T}\log f\right)\approx 1-i\frac{\omega}{4\pi
T}\log f+\mathcal{O}\left(\frac{\omega}{T}\right)^{2}$ (19)

The matching condition $\delta B_{xy}^{inner}=\delta B_{xy}^{outer}$ in this
region prompts yield the following algebraic relations between the boundary
quantities $b_{xy},\langle J_{xy}\rangle$:

$\displaystyle\frac{i\omega}{4\pi T}c^{H}$
$\displaystyle=\left(\frac{\mathcal{B}/r_{h}}{3r_{h}^{2}f^{\prime}(r_{h})}\right)\langle
J_{xy}\rangle\,$ (20) $\displaystyle c^{H}$
$\displaystyle=b_{xy}+\left[\frac{1}{\kappa(\Lambda)}+\log\left(\frac{\Lambda}{r_{h}}\right)+\phi(r_{h})\right]\langle
J_{xy}\rangle\,.$

Solving these equations at vanishing source $b_{xy}=0$ yields the spectrum of
the form $\omega=-i/\tau_{E^{\parallel}}$ where $\tau_{E^{\parallel}}$ is the
lifetime of the electric flux parallel to the equilibrium magnetic field. This
is the first key result that we advertised earlier, namely

$\tau_{E^{\parallel}}=\frac{2\pi
T}{\mathcal{B}}\left(e_{r}^{-2}+\phi(r_{h})\right)\,,$ (21)

where we write $e_{r}^{-2}=\log(\Lambda/r_{h})+\kappa(\Lambda)^{-1}$ which
plays the role of renormalised electromagnetic coupling. More details on the
$T/\sqrt{\mathcal{B}}$ dependence of $\phi(r_{h})$ can be found in Appendix A.

What does this result tell us about the lifetime of the electric flux
operator? While the integral $\phi(r_{h})$ can be a dimensionless function of
$T$ and $\mathcal{B}$, the renormalised electromagnetic coupling can be chosen
in such a way that $e_{r}^{-2}\gg\phi(r_{h})$ and
$e_{r}^{-2}T^{2}/\mathcal{B}\ll 1$ so that
$\omega\tau_{E^{\parallel}}\sim\omega/T\ll 1$. The second limit is essential
as the matching procedure assumes that $\omega/T\ll 1$ and the solution
outside this regime has to be discarded. Taking these factors into account,
one concludes that the temperature dependence of the electric flux is
different from the high temperature $T/\sqrt{\mathcal{B}}\gg 1$ limit where
$\tau_{E}\sim 1/T$ ( see Fig 2). Naively taking the limit $T\to 0$ in (21)
will result in the vanishing lifetime of the electric flux in contrast to the
result in (2). However, one has to carefully remove the limit $\omega/T\ll 1$
in order to access the lower temperature limit $\omega/T\sim
1,\omega/\sqrt{\mathcal{B}}\ll 1$.

Figure 2:  A sketch of the decay rate (inverse of the lifetime) of the
electric field as a function of $T/\sqrt{\mathcal{B}}$, measured in the unit
of $\sqrt{\mathcal{B}}$. The high temperature regime (red) depict the result
of decay rate at zero magnetic field found in Hofman_2018 ; Grozdanov_2019
which has the same temperature dependence as in (1)-(2). In the low
temperature regime (blue), however, the operator lifetime becomes those found
in (21).

####  II.1.2 Perturbation perpendicular to equilibrium magnetic field

Unlike the previous case, the perturbation $\delta B_{xz}$ that corresponds to
$E^{\perp}=\frac{1}{2}\epsilon^{yzx}\langle J_{zx}\rangle$ is coupled to the
metric perturbation. This is manifest in the equations of motion

$\displaystyle\frac{d}{dr}\Big{(}r^{2}fe^{-W}\delta
B_{xz}^{\prime}+\mathcal{B}(\delta
G^{x}_{t})\Big{)}+\frac{\omega^{2}e^{-W}}{r^{2}f}\delta B_{xz}$
$\displaystyle=0\,,$ (22)
$\displaystyle\frac{d}{dr}\left(e^{4V+W}(\delta{G^{x}_{t}})^{\prime}+4\mathcal{B}\delta
B_{xz}\right)$ $\displaystyle=0\,,$

where $\delta G_{\mu\nu}$ denotes the metric perturbations. Note that the
coupled perturbations $\\{\delta B_{xz},\delta G_{tx}\\}$ and
$\\{B_{yz},\delta G_{yz}\\}$ are equivalent due to $SO(2)$ symmetry in the
plane perpendicular to the equilibrium magnetic field. Also, the second
equation of motion in (22) can be written in a total derivative form
$d\pi_{tx}/dr=0$ with $\pi_{tx}$ is related to the momentum $\langle
T^{tx}\rangle$. Since we are working in the zero wavevector limit, the
conservation of momentum implies that $\pi_{tx}=0$ in Fourier space (which can
be shown explicitly using the $rx-$component of the Einstein equation).

The solution for $\delta B_{xz},\delta G_{tx}$ in the outer region can be
found by using the property of the background geometry combined with the
Wronskian method as in Grozdanov_2019 . To be more precise, one first notes
that the time-independent solution of the magnetised black brane can be
written in a total derivative form, which implies the existence of two
radially conserved currents.

$\displaystyle Q_{1}=r^{2}f(V^{\prime}-W^{\prime})e^{2V+W}+2\mathcal{B}h(r)$
$\displaystyle=0\,,$ (23a) $\displaystyle
Q_{2}=e^{4V+W}\frac{d}{dr}\left(e^{-2V}r^{2}f\right)-4\mathcal{B}h(r)$
$\displaystyle=sT\,,$ (23b) where we write the equilibrium ansatz for the
gauge field as $B=h(r)\,dt\wedge dz$ with gauge choice $h(r_{h})=0$, which,
together with the horizon regularity, sets $Q_{1}=0$. The relation between
$h(r)$ and the 3-form field strength is $e^{2V-W}h^{\prime}=\mathcal{B}\,.$
(23c)

More details on obtaining these radially conserved quantities can be found in
e.g. Gubser:2009cg . With this ansatz, we can compare (22) and (23) and find
that one of the solutions of (22) when $\omega/r\to 0$ are

$\delta B_{xz}=\Phi_{1}(r)=h(r)+\frac{sT}{4\mathcal{B}}\,,\qquad\delta
G^{x}_{t}=\Psi_{1}(r)=-e^{-2V}r^{2}f\,.$ (24)

One can use the Wronskian method to find find a pair of solution of (22) that
are linearly independent to $\\{\Phi_{1},\Psi_{1}\\}$. These solutions are

$\displaystyle\Phi_{2}(r)=\frac{1}{4\mathcal{B}}-\int^{\infty}_{r}d\mathtt{r}\left(\frac{\mathcal{B}e^{W(\mathtt{r})}\Psi_{2}(\mathtt{r})}{\mathtt{r}^{2}f(\mathtt{r})}\right)\,,\qquad\Psi_{2}(r)=\Psi_{1}(r)\int^{\infty}_{r}d\mathtt{r}\left(\frac{e^{-W(\mathtt{r})}}{\mathtt{r}^{4}f(\mathtt{r})^{2}}\right)$
(25)

As a result, the outer region solution can be written as

$\begin{pmatrix}\delta B_{xz}^{outer}\\\ (\delta
G^{x}_{t})^{outer}-\frac{1}{\mathcal{B}}\mathcal{J}_{xz}\end{pmatrix}=c_{1}\begin{pmatrix}\Phi_{1}\\\
\Psi_{1}\end{pmatrix}+c_{2}\begin{pmatrix}\Phi_{2}\\\ \Psi_{2}\end{pmatrix}$
(26)

where $\mathcal{J}_{xz}:=(r^{2}fe^{-W}\delta B_{xz}^{\prime}+\mathcal{B}\delta
G^{t}_{x})$ is an integration constant of (22) at $\omega=0$. One can
substitute the $BTZ\times\mathbb{R}^{2}$ ansatz into the solution in (26) to
check that $\Phi_{1},\Psi_{1,2}$ are finite at $r=r_{h}$ while $\Phi_{2}$ is
singular. It is convenient to separate out the singular part of $\Phi_{2}$ in
the following form

$\Phi_{2}(r)=\phi_{2}(r)-\left(\frac{\mathcal{B}e^{W}\Psi_{2}}{r^{2}f^{\prime}}\right)_{r=r_{h}}\log
f(r)$ (27)

where $\phi_{2}(r)$ is the integral in (25) with the logarithmic divergence
subtracted. The boundary condition where the source for both metric and 2-form
gauge field fluctuation vanishes corresponds to the following values of
$c_{1}$ and $c_{2}$

$c_{1}=\frac{\mathcal{J}_{xz}}{\mathcal{B}}\,,\qquad
c_{2}=-4\left(\frac{sT}{4\mathcal{B}}+h(\Lambda)+\frac{\mathcal{B}}{\hat{\kappa}(\Lambda)}\right)\mathcal{J}_{xz}$
(28)

One can also check that $\mathcal{J}_{xz}$ is identical to the one-point
function $\langle\delta J^{xz}\rangle$ via the definition (12). Note also that
the ratio $c_{2}/\mathcal{J}_{xz}$ is finite due to the cancellation of the
logarithmic divergence of $1/\kappa(\Lambda)$ and that of the near boundary
solution of $h(r)$, obtained via (23c).

Let us also pointed out another way to organise the equations of motion for
$\delta B_{xz}$. It turns out that (22) can be combined into a single equation
of motion that reduces to a total derivative form at $\omega=0$. Following the
procedure in e.g. Davison:2015taa  and some manipulation, we find

$\left([e^{4V+W}\left(e^{-2V}r^{2}f\right)^{\prime}]^{2}r^{2}fe^{-W}\delta\tilde{B}_{xz}^{\prime}\right)^{\prime}+\frac{\omega^{2}}{r^{2}fe^{W}}[e^{4V+W}(e^{-2V}r^{2}f)^{\prime}]^{2}\delta\tilde{B}_{xz}=0$
(29)

where $\delta\tilde{B}_{xz}=\delta B_{xz}/[e^{4V+W}(e^{-2V}r^{2}f)^{\prime}]$.
The outer region solution of (29) is easily obtained and can be shown to be
identical to those of (26).

We can now proceed to the inner region solution. This can be found by solving
Eq.(29) and one find

$\delta B^{inner}_{xz}=c^{H}\exp\left(-\frac{i\omega}{4\pi T}\log
f(r)\right)\,.$ (30)

In the intermediate region, we apply the expansion in (19). The coefficients
$c_{1},c_{2}$ are related to $c^{H}$ via

$\displaystyle\left(-\frac{i\omega}{4\pi T}\right)c^{H}$
$\displaystyle=-\left(\frac{\mathcal{B}e^{W}\Psi_{2}}{r^{2}f^{\prime}}\right)_{r=r_{h}}c_{2}\,,$
(31) $\displaystyle c^{H}$
$\displaystyle=\left(\frac{sT}{4\mathcal{B}}\right)c_{1}+\phi_{2}(r_{h})c_{2}\,,$

Substituting the form of $c_{1},c_{2}$ in terms of $\langle\delta
J_{xz}\rangle$, we can write the relations in a form similar to $\langle\delta
J_{xy}\rangle$, namely

$\left(-i\omega+\frac{1}{\tau_{E^{\perp}}}\right)\langle\delta
J_{xz}\rangle=0\,.$ (32)

In the case of vanishing sources, we can write
$\frac{c_{2}}{c_{1}}=-4\mathcal{B}\left(\frac{sT}{4\mathcal{B}}+h(\Lambda)+\frac{\mathcal{B}}{\kappa(\Lambda)}\right)$
and the relaxation time of the electric field perpendicular to the equilibrium
magnetic field is

$\tau_{E^{\perp}}=\frac{\sqrt{3}}{2\pi
T\mathcal{B}\Psi_{2}(r_{h})}\left[\frac{sT}{4\mathcal{B}}\frac{c_{1}}{c_{2}}+\phi_{2}(r_{h})\right]$
(33)

In contrast to the result at $e_{r}^{-2}\gg 1$ and zero equilibrium magnetic
field in Hofman_2018 ; Grozdanov_2019 , the lifetime at strong magnetic field
$\mathcal{B}/T^{2}$ has a very different form. To see this, it is useful to
examined that the combinations that enter the lifetime as follows

$\Psi_{2}(r_{h})\propto\frac{1}{\mathcal{B}T^{2}}\,,\qquad\phi_{2}(r_{h})\propto\frac{1}{\mathcal{B}}\,,\qquad\frac{c_{1}}{c_{2}}\propto\frac{1}{\mathcal{B}^{2}}\,\quad\text{for
large}1/\kappa(\Lambda)$ (34)

with proportionality constants given by some numbers of order
$\mathcal{O}(1)$. In the limit of large electromagnetic coupling
$1/\kappa(\Lambda)\gg 1$ and $\mathcal{B}/T^{2}\gg 1$, we find that this gives
a short lifetime of the form $\tau_{E^{\perp}}\propto T/\mathcal{B}$. However,
the location of this decaying mode $\omega=-i/\tau_{E^{\perp}}$ lies outside
the hydrodynamic regime $\omega/T\ll 1$. Thus, one conclude that there are no
modes with long lifetime in this regime777Note also that, if one were to
perform this analysis for a perturbation in the holographic dual to a theory
with zero-form $U(1)$ at $T>0,\mu=0$ (as in Kovtun:2003wp , see also
Grozdanov_2019 ), one would find a spectrum of the form $\omega\sim T$. This
solution is spurious as it lies outside the hydrodynamic regime $\omega/T\ll
1$ and, in fact, is not present in the genuine spectrum obtained numerically
at finite $\omega/T$ Kovtun:2005ev . .

###  II.2 Checking $T\gtrsim 0$ limit in $\omega/\sqrt{|\bf{B|}}\ll 1$ regime

While the result in the previous section strongly indicated that the electric
flux lifetime becomes very short at extremely low temperature, the simplicity
of the holographic model also allows us to extend the analysis beyond the
usual hydrodynamic $\omega/T\ll 1$ regime. We will first show that the zero
temperature theory does not support the purely decaying mode of the form
$\omega=-i/\tau$ in the small $\omega/\sqrt{|{\bf B}|}$ regime. Next, we
further extend the regime of validity to that of $\omega/\sqrt{|{\bf B}|}\ll
1$ but for arbitrary $\omega/T$. The purpose of the latter is to show that
$\tau_{E}\propto T/\sqrt{|{\bf B}|}$ without relying on the $\omega/T\ll 1$
limit.

####  II.2.1 Zero temperature

A simple argument for the non-existence of such a slowly decaying mode, is the
presence of Lorentz symmetry at zero temperature on the $AdS_{3}$ submanifold
in the deep infrared. On the other hand, one can also show this, using
matching methods similar to those in DHoker:2010xwl ; DHoker:2010onp ;
Davison:2013bxa .

To obtain this result, one first realises that the geometry of the magnetised
black brane is that of an interpolation between IR
$AdS_{3}\times\mathbb{R}^{2}$ and UV $AdS_{5}$. Roughly speaking, the IR
geometry starts to becomes a good approximation as one starts to probe the
scale below the magnetic field i.e. $r\sim\sqrt{|\bf{B}|}$. The inner and
outer regions are defined such that they start off from the IR and UV geometry
respectively, and extend to cover the overlap region (see Figure 3). This is
achievable when $\omega/\sqrt{|{\bf B}|}\ll 1$.

Figure 3:  A sketch of the bulk geometry at zero temperature. The inner
region, whose solutions only depends on the ratio $\omega/r$ extended from the
near horizon limit $r\to 0$ to the one where $\omega/r\sim\omega/\sqrt{|{\bf
B}|}\to 0$ as we are working in the $\omega/\sqrt{|{\bf B}|}\ll 1$ limit. The
outer region is defined to be the region where the $\omega^{2}/r^{2}$ and
higher power in $\omega/r$ is suppressed, which can be extended toward
$r\gg\sqrt{|\bf{B}|}$ as long as the frequency is small.

For concreteness, let us demonstrate how this works in the $E^{\parallel}$
channel that involves the bulk field $\delta B_{xy}$ governed by Eq.(15). The
solution can be written in the same form as (17) evaluated at zero temperature
(i.e. $r_{h}=0$). It is worth noting that the singular behaviour near
$r/\sqrt{\mathcal{B}}\to 0$ is different from that in earlier section.
Instead, it can be written as

$\delta
B_{xy}^{outer}(r)=c_{1}-c_{2}\left(\log\Lambda-\bar{\phi}(r)+\frac{\mathcal{B}/3}{6r^{2}}\right)+\mathcal{O}\left(\frac{\omega^{2}}{r^{2}},\frac{\omega^{2}}{r^{2}}\log\left(\frac{\omega}{r}\right)\right)$
(35)

where the integration constants can be related to source and response via
(12). It is worth noting that the logarithmic divergence appears at order
$\omega^{2}$. This is can be confirmed via Frobenius analysis in $AdS_{5}$
region (see e.g. Kovtun:2006pf ) and $AdS_{3}\times\mathbb{R}^{2}$ region (see
appendix B). The prefactor of the $r^{-2}$ divergence is obtained by
evaluating $e^{2V(r)-W(r)}/f(r)$ at the horizon $r\to 0$. Here $\bar{\phi}(r)$
is the integral in (17) subtracted by the $r^{-2}$ divergent and logarithmic
divergent pieces. The resulting integral evalutated from
$r=r_{0}\sim\sqrt{\mathcal{B}}$ of the overlapping region to the UV cutoff
$r=\Lambda$ is finite and its number is not extremely relevant for us as long
as one keep $e_{r}^{-2}$ large.

Next, we consider the inner region solution, which can be obtained by solving
(15) in the $AdS_{3}\times\mathbb{R}^{2}$ region. Upon imposing horizon
regularity at $r\to 0$, we find that the inner region solution is

$\delta B_{xy}^{inner}=c^{H}\zeta
K_{1}(\zeta)\,,\qquad\zeta=\frac{3\omega}{r}$ (36)

For these two branches of solutions to match, we extend the inner region
solution to the regime where $\zeta=\omega/r\ll 1$. We find that the ‘near
boundary’ expansion takes the form

$\delta
B_{xy}^{inner}=c^{H}\left(1+\frac{1}{2}\gamma\zeta^{2}+\frac{1}{2}\zeta^{2}\log\zeta+...\right)$
(37)

Matching this solution to the outer region, we find that
$c_{2}\propto\omega^{2}$ unlike what happened in the previous section.
Carrying on the matching procedure, we find that the polynomial governing the
spectrum only depends on $\omega^{2}$ and thus rules out the purely imaginary
mode $\omega=-i/\tau$. The same argument can also be made for the $E^{\perp}$
channel involving $\delta B_{xz}$. This is because, the part that is relevant
to the matching procedure only depends on $\zeta^{2}$. See appendix B for more
details on the form of $\delta B_{xz}$ in the $AdS_{3}\times\mathbb{R}^{2}$
region.

####  II.2.2 $T\lesssim\omega\ll\sqrt{|{\bf B}|}$ limit

In this section, we show that the electric flux lifetime can also be obtained
regime where $\omega/T\gtrsim 1$ and $\omega/\sqrt{|\bf{B}|}\ll 1$ while
keeping $\sqrt{|\bf{B}|}/T\gg 1$. The calculations closely resembles that of
the zero temperature case except that the deep IR geometry is now
$BTZ\times\mathbb{R}^{2}$ instead of $AdS_{3}\times\mathbb{R}^{2}$. Figure 4
illustrates this geometry where the $AdS_{5}$ joined with the
$BTZ\times\mathbb{R}^{2}$ at the ‘boundary’ $AdS_{3}\times\mathbb{R}^{2}$ of
the IR geometry. We will only focus on the $E^{\parallel}$ fluctuations as it
is the only channel that contains the decaying modes in the $\omega/T\ll 1$
regime. Similar computation for this type of geometry can also be found in
DHoker:2010onp .

Figure 4:  A sketch of the bulk geometry at low temperature
$T\ll\sqrt{|\bf{B}|}$. The inner region, whose solutions only depends on the
ratio $\omega/r$ extends from the near horizon limit $r\to
r_{h}\ll\sqrt{\mathbf{B}}$ to the one where $\omega/r\sim\omega/\sqrt{|{\bf
B}|}\ll 1$, which corresponds to the near boundary region of
$BTZ\times\mathbb{R}^{2}$ geometry, described by
$AdS_{3}\times\mathbb{R}^{2}$. The outer region is defined to be the region
where $\omega^{2}/r^{2}$ (and higher powers) is negligible and, therefore, can
be extended toward $r\sim\sqrt{|\bf{B}|}$ in the $\omega/\sqrt{|\bf{B}|}\ll 1$
limit.

The outer region solution, which extends from the UV $AdS_{5}$ to the
intermediate $AdS_{3}\times\mathbb{R}^{2}$ region has the same form as in
(35). This is possible only in the limit where $\sqrt{\mathcal{B}}\gg T$ so
that $r/\sqrt{\mathcal{B}}$ is always much greater than
$T/\sqrt{\mathcal{B}}\sim r_{h}/\sqrt{\mathcal{B}}$ in this region.

The inner region solution in the $BTZ\times\mathbb{R}^{2}$ region can be
expressed in terms of a hypergeometric function (upon imposing ingoing
boundary condition)

$\delta
B_{xy}^{inner}=c^{H}\left(1-\frac{r_{h}^{2}}{r^{2}}\right)^{-i\mathfrak{w}/2}\,{{}_{2}F_{1}}\left(-\frac{i\mathfrak{w}}{2},-\frac{i\mathfrak{w}}{2},-\frac{i\mathfrak{w}}{2};1-\mathfrak{w};1-\frac{r_{h}^{2}}{r^{2}}\right)$
(38)

where $\mathfrak{w}=\omega/(2\pi T)=\omega/3r_{h}$. Extending this solution in
the $r\gg r_{h}$ limit (which is possible due to $r_{h}/r\to 0$ as we approach
the limit $\omega/r\to 0$) yields the following expansion abramowitz+stegun

$\delta B_{xy}^{inner}\propto c^{H}\left[1+\frac{i\omega
r_{h}}{6r^{2}}+\frac{1}{4}\left(\frac{\omega}{3r}\right)^{2}\left(2-2\gamma-2\psi(1-i\mathfrak{w}/2)-\log\left(\frac{r_{h}^{2}}{r^{2}}\right)\right)+\mathcal{O}(\omega^{3})\right]$
(39)

where $\psi(x)$ is the digamma function and the constants of proportionality
are combinations of gamma functions that can be absorbed in the definition of
$c^{H}$. The first two terms in $[...]$ are what important for us. By working
to leading order in $\omega/r\ll 1$ as one approaches the intermediate
$AdS_{3}\times\mathbb{R}^{2}$ region, we find the following matching solution

$c_{1}-c_{2}\log(\Lambda/\sqrt{\mathcal{B}})+\bar{\phi}=c^{H}\,,\qquad\left(\frac{\mathcal{B}}{3}\right)c_{2}=i\omega\left(\frac{2\pi
T}{3}\right)c^{H}$ (40)

We can convert $c_{1}$ to the source $b_{xy}$ and $c_{2}$ as done in the
previos sections. Upon taking $e_{r}^{-2}\gg\bar{\phi}$ (so that the solution
lies in the regime of validity $\omega/\sqrt{\mathcal{B}}\ll 1$), we find the
solution of the form $\omega=-i/\tau_{E^{\parallel}}$ where
$\tau_{E^{\parallel}}$ is the same as in (21). This indicates that the
lifetime indeed grows as $T/\sqrt{\mathcal{B}}$ increases regardless of the
ratio $\omega/T$.

##  III Conclusion

The higher-form symmetry viewpoint of magnetohydrodynamics and its low
temperature incarnation, the force-free electrodynamics, leads to new
insights. The central focus of the present work was to established the absence
of long-lived non-conserved operators. In turn, this indicates the validity of
a hydrodynamic description at low temperature and strong magnetic field. The
question of whether the only operators that govern the deep IR dynamics are
the conserved charges is important and ought to be asked before any
quantitative attempt is made to study hydrodynamic properties (such as shear
viscosity etc). All non-conserved operators must decay much faster than the
scale of interest if a hydrodynamic interpretation is to be meaningful.

We work with a holographic model which shares the same global symmetry as that
of the plasma, namely only the energy, momentum and magnetic flux commute with
the Hamiltonian. The model is simple enough for the lifetime of electric flux
to be determined by classical bulk dynamics and the precise question is
whether or not the electric flux is sufficiently long-lived to interfere with
hydrodynamic modes. Due to the anisotropy of the system in the presence of a
strong expectation value of magnetic field, the lifetime of the electric field
depends on its orientation. Our results can be summarised as follows

  * •

For electric flux $E_{\parallel}$ parallel to the magnetic field, the lifetime
has a strong dependence on the double-trace coupling $\kappa$ which plays a
role similar to the renormalised electromagnetic coupling. In the extreme
limit of $e^{-2}_{r}\gg|\bf{B}|/T^{2}$, the lifetime can be large enough to be
detectable by the analytic computation in both the ‘usual’ hydrodynamic regime
$\omega/T\ll 1$ and on even lower temperature regime where
$\omega/\sqrt{\bf{B}}\ll 1$ while $\omega/T$ may remains finite. We found that
the lifetime becomes shorter as one decreases the ratio of
$T/\sqrt{|\bf{B}|}$. The latter indicates that the lifetime will become
extremely short in the extremely strong magnetic field regime
$T/\sqrt{|\bf{B}|}\ll 1$ and cannot interfere with the low energy regime of
$\omega/\sqrt{|\bf{B}|}\ll 1$ where the FFE limit is thought to be applicable.

  * •

For the component of electric flux $E_{\perp}$ perpendicular to the magnetic
field, we find that there is no pole in the vicinity of $\omega/T\ll 1$. The
dependence of the lifetime on the renormalised electromagnetic coupling
disappears as one approaches the strong magnetic field limit.

We also performed a consistency check at $T\to 0$ to ensure that there are no
modes in the deep IR limit of $\omega/\sqrt{|\bf{B}|}\ll 1$. In this regime,
the modes that indicate (potentially) long lifetime of $E_{\parallel}$
disappear from the low energy spectrum as anticipated.

These computations are basic checks on the validity of FFE description. In the
holographic context, it would be interesting to check if all the accessible
non-conserved operator truly have a parametrically short lifetime as well as
confirming the low energy spectrum predicted by force-free electrodynamics
(and its subsequent derivative corrections). Extraction of FFE effective
action from gravity akin to Nickel:2010pr ; Glorioso:2018mmw ; deBoer:2018qqm
or the full constitutive relation as in Bhattacharyya:2008jc ; Banerjee:2008th
; Erdmenger:2008rm  would be desirable as a definitive proof of FFE
description in the dynamically magnetised black brane geometry. Last but not
least, it would be very interesting to investigate operators lifetime in
(weakly coupled) quantum electrodynamics at finite $T$ and ${\bf B}$ to better
understand FFE and its limitations in a system more directly connected to
astrophysical plasma than the strongly coupled holographic model considered
here.

## Acknowledgements

We would like to thank Jay Armas, Sašo Grozdanov, Nabil Iqbal, Kieran
Macfarlane, Watse Sybesma and Lárus Thorlacius for helpful discussions and
comments. We are particularly grateful to S. Grozdanov, N. Iqbal and L.
Thorlacius for commenting on the manuscript. The work of N. P. was supported
by Icelandic Research Fund grant 163422-052 and STFC grant number
ST/T000708/1. The work of A.R was supported in part by the Icelandic Research
Fund under grant 195970-052 and by the University of Iceland Research Fund.

##  Appendix A Numerical solution and evaluation of operators lifetime

In this section, remarks on the evaluation of the electric flux are
elaborated. The numerical background solution for this geometry can be
constructed in the same way as DHoker:2009mmn  using shooting method. The
solution is a one-parameter family characterised by $\mathcal{B}/T^{2}$ which
allows us the freedom to choose $r_{h}=1,r_{h}^{2}f^{\prime}(r_{h})=1$ (or
equivalently $T=1/4\pi$). It is also convenient to set $V(r_{h})=W(r_{h})=0$
which results in the UV boundary metric of the form

$\lim_{r\to\infty}ds^{2}=r^{2}\left(-dt^{2}+v(dx^{2}+dy^{2})+wdz^{2}\right)+\frac{dr^{2}}{r^{2}}$
(41)

Upon rescaling of spatial coordinates
$\\{dx,dy,dz\\}\to\\{dx/\sqrt{v},dy/\sqrt{v},dz/\sqrt{w}\\}$, we recover the
desired background solutions. Note also that the physical magnetic flux is
related to the input parameter (that produced the metric in (41)) by
$\mathcal{B}_{\text{physical}}=\mathcal{B}_{\text{input}}/v$. A small caveat
of this method is that one cannot find a smooth solution beyond
$\mathcal{B}_{\text{input}}\gtrsim\sqrt{3}/2$ which corresponds to the
temperature
$T/\sqrt{\mathcal{B}}=(4\pi\sqrt{\mathcal{B}_{\text{input}}/v})^{-1}\approx
0.05$. This is most likely an artifact of the presented numerical method as
there exists a smooth solution in the zero temperature limit corresponding to
the $AdS_{3}\times\mathbb{R}^{2}$ geometry in the deep IR. We should also note
that this is a sufficiently low energy temperature as the entropy becomes
sufficiently close to $s\propto T$ obtained from $BTZ\times\mathbb{R}^{2}$
geometry (c.f. DHoker:2009mmn ; Grozdanov_2019b ). The background is generated
for $r$ from $[1+10^{-3},10^{6}]$ and varying the (numerical) cutoffs within
this order of magnitude does not change the obtained numerical results.

Let us also remark on the the numerical value of the renormalised
electromagnetic coupling
$e_{r}^{-2}=\log(\Lambda/r_{h})+\kappa(\Lambda)^{-1}$. This quantity strongly
influences both the thermodynamics and low energy spectrum Grozdanov_2019b ;
Hofman_2018 ; Grozdanov_2019  of the model. In particular a small value of
$e_{r}^{-2}$ would result in the speed of sound becoming imaginary
Grozdanov_2019b . Another way to see that this quantity should be large is to
write it in terms of a renormalisation group independent scale $M_{*}$ that
denotes the energy scale of a Landau pole Hofman_2018  i.e.
$e_{r}^{-2}\sim\log(M_{\star}/T)$ where $M_{\star}\gg T$. We take this to be
the largest scale in the problem–much larger than the accessible value of
$\sqrt{\mathcal{B}}/T$.

Figure 5:  Numerical evaluation of $\phi(r_{h})$ in (21) as a function of
$T/\sqrt{\mathcal{B}}$. The black dots denote the numerical evaluation while
the red line denotes the fitting function for small $T/\sqrt{\mathcal{B}}$ as
$\phi\approx-(0.008)\frac{\mathcal{B}}{T^{2}}\log(5.7\mathcal{B}/T^{2})$. For
high temperatures, the value of $\phi(r_{h})$ is approximately constant around
$0.69$. The value of $\phi(r_{h})$ at lowest achievable temperature is at
$\phi(r_{h})=-23.49$.

Numerical value of the integral for $\phi(r_{h})$ in (21) is shown in Figure
5. For a larger temperature (when $\phi(r_{h})\approx\mathcal{O}(1)$), the
lifetime can be sensibly approximated to be $\tau_{E^{\parallel}}\approx
2\pi(T/\mathcal{B})e_{r}^{-2}$. As $T/\sqrt{\mathcal{B}}$ decreases, the
lifetime becomes shorter and, if we are to extrapolate the fitting function
$\phi\sim\frac{\mathcal{B}}{T^{2}}\log\frac{\mathcal{B}}{T^{2}}$ to even lower
temperature where $e_{r}^{-2}\gtrsim\phi$, it will escape the regime of the
validity of small $\omega/T,\omega/\sqrt{\mathcal{B}}$ expansions. In this
scenario, one shall conclude that there are no long-lived modes that can
interfere with the low energy excitations.

##  Appendix B Frobenius analysis in $AdS_{3}\times\mathbb{R}^{2}$ region

Consider the equation of motion for $\delta B_{xy}$ in the intermediate
$AdS_{3}\times\mathbb{R}^{2}$ region:

$\delta B^{\prime\prime}_{xy}(r)+\frac{3}{r}\delta
B_{xy}^{\prime}(r)+\frac{\omega^{2}}{9r^{2}}\delta B_{xy}(r)=0$ (42)

The solution in this region can be obtained via Frobenius method. More
precisely, one can change the radial coordinate into $\zeta=3\omega/r$ and
redefine $\delta B_{xy}=\zeta c(\zeta)$. It follows that $c(\zeta)$ is the
solution of the Bessel equation of order $1$, which has a regular singular
point at $\zeta=0$. The near-boundary $r\to\infty$, or equivalently $\zeta\to
0$, akin to the Fefferman-Graham expansion in the usual holographic
renormalisation, can be written as

$\delta
B_{xy}(\zeta)=c_{1}^{M}\mathcal{P}_{1}(\zeta)+\Big{(}c_{2}^{M}+\mathfrak{h}\log\zeta\Big{)}\mathcal{P}_{2}(\zeta)$
(43a)

where $c^{M}_{1},c^{M}_{2}$ are integration constants and
$\mathcal{P}_{i}(\zeta)$ are regular polynomials of the following form

$\mathcal{P}_{1}=1+\sum_{n=1}^{\infty}p_{1}^{[n]}\zeta^{n}\,,\qquad\mathcal{P}_{2}=\zeta^{2}\left(1+\sum_{n=1}^{\infty}p_{2}^{[n]}\zeta^{n}\right)$
(43b)

Similar to the usual procedure in the holographic renormalisation
deHaro:2000vlm , all the coefficients $p_{1}^{[n]},p_{2}^{[n]},\mathfrak{h}$
except $p_{1}^{[2]}$, which can be set to zero without loss of generality
Kovtun:2006pf , can be obtained recursively. The important piece of
information here is the coefficient $\mathfrak{h}=1$ which can be obtained by
recursively solving the equation (42). Another easy way to see this is to
recast (42) as the Bessel equation of order 1 as pointed out earlier. Then,
using the fact that the Bessel functions $K_{1}(\zeta)$ and $I_{1}(\zeta)$ are
two independent solutions of such equation and, for small $\zeta$ they admit
the following asymptotic expansions (see e.g. §3.3 of bender2013advanced )

$I_{1}(\zeta)=\frac{\zeta}{2}+\frac{\zeta^{2}}{16}+\mathcal{O}(\zeta)^{3}\,,\qquad
K_{1}(\zeta)=\left(\gamma+\log\frac{\zeta}{2}\right)I_{1}(\zeta)+\frac{1}{\zeta}$
(44)

will result in the series expansions of the solution in
$AdS_{3}\times\mathbb{R}^{2}$ region in (43a).

A similar procedure can also be applied for $E^{\perp}$ using Eq.(29).
Substituting $\delta\tilde{B}_{xz}=\zeta^{2}c(\zeta)$, one finds that it obeys
the Bessel equation of order $2$ whose $\zeta\ll 1$ expansion only yields even
power in $\zeta$.

## References

  * (1) L. D. Landau and E. M. Lifshitz, Fluid Mechanics.  Butterworth-Heinemann, 2nd ed., 1987.
  * (2) R. W. Zwanzig, Statistical mechanics of irreversibility.  Lectures on Theoretical Physics Volume 3, 139 (Interscience, 1961), 1961\.
  * (3) D. Forster, Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Functions.  Perseus Books, 1995.
  * (4) S. A. Hartnoll and D. M. Hofman, “Locally Critical Resistivities from Umklapp Scattering,” Phys. Rev. Lett. 108 (2012) 241601, arXiv:1201.3917 [hep-th].
  * (5) D. Gaiotto, A. Kapustin, N. Seiberg, and B. Willett, “Generalized Global Symmetries,” JHEP 02 (2015) 172, arXiv:1412.5148 [hep-th].
  * (6) S. Grozdanov, D. M. Hofman, and N. Iqbal, “Generalized global symmetries and dissipative magnetohydrodynamics,” Phys. Rev. D 95 no. 9, (2017) 096003, arXiv:1610.07392 [hep-th].
  * (7) D. Schubring, “Dissipative String Fluids,” Phys. Rev. D 91 no. 4, (2015) 043518, arXiv:1412.3135 [hep-th].
  * (8) J. Hernandez and P. Kovtun, “Relativistic magnetohydrodynamics,” JHEP 05 (2017) 001, arXiv:1703.08757 [hep-th].
  * (9) J. Armas and A. Jain, “Magnetohydrodynamics as superfluidity,” Phys. Rev. Lett. 122 no. 14, (2019) 141603, arXiv:1808.01939 [hep-th].
  * (10) J. Armas and A. Jain, “One-form superfluids & magnetohydrodynamics,” JHEP 01 (2020) 041, arXiv:1811.04913 [hep-th].
  * (11) W. G. Dixon, Special relativity: the foundation of macroscopic physics.  CUP Archive, 1982.
  * (12) A. M. Anile, Relativistic fluids and magneto-fluids: With applications in astrophysics and plasma physics.  Cambridge University Press, 2005.
  * (13) S. S. Komissarov, “A Godunov-type scheme for relativistic magnetohydrodynamics,” Monthly Notices of the Royal Astronomical Society 303 no. 2, (02, 1999) 343–366.
  * (14) P. B. Arnold, G. D. Moore, and L. G. Yaffe, “Transport coefficients in high temperature gauge theories. 1. Leading log results,” JHEP 11 (2000) 001, arXiv:hep-ph/0010177.
  * (15) R. Blandford and R. Znajek, “Electromagnetic extractions of energy from Kerr black holes,” Mon. Not. Roy. Astron. Soc. 179 (1977) 433–456.
  * (16) S. Komissarov, “Electrodynamics of black hole magnetospheres,” Mon. Not. Roy. Astron. Soc. 350 (2004) 407, arXiv:astro-ph/0402403.
  * (17) P. Goldreich and W. H. Julian, “Pulsar Electrodynamics,” Astrophys. J.  157 (Aug., 1969) 869.
  * (18) T. Wiegelmann and T. Sakurai, “Solar Force-free Magnetic Fields,” Living Rev. Sol. Phys. 9 (2012) 5, arXiv:1208.4693 [astro-ph.SR].
  * (19) S. E. Gralla and T. Jacobson, “Spacetime approach to force-free magnetospheres,” Mon. Not. Roy. Astron. Soc. 445 no. 3, (2014) 2500–2534, arXiv:1401.6159 [astro-ph.HE].
  * (20) G. Compère, S. E. Gralla, and A. Lupsasca, “Force-Free Foliations,” Phys. Rev. D 94 no. 12, (2016) 124012, arXiv:1606.06727 [math-ph].
  * (21) T. Uchida, “Theory of force-free electromagnetic fields. i. general theory,” Phys. Rev. E 56 (Aug, 1997) 2181–2197. https://link.aps.org/doi/10.1103/PhysRevE.56.2181.
  * (22) C. Thompson and O. Blaes, “Magnetohydrodynamics in the extreme relativistic limit,” Phys. Rev. D 57 (1998) 3219–3234.
  * (23) S. E. Gralla and N. Iqbal, “Effective Field Theory of Force-Free Electrodynamics,” Phys. Rev. D 99 no. 10, (2019) 105004, arXiv:1811.07438 [hep-th].
  * (24) P. Glorioso and D. T. Son, “Effective field theory of magnetohydrodynamics from generalized global symmetries,” arXiv:1811.04879 [hep-th].
  * (25) B. Benenowski and N. Poovuttikul, “Classification of magnetohydrodynamic transport at strong magnetic field,” arXiv:1911.05554 [hep-th].
  * (26) S. Grozdanov and N. Poovuttikul, “Generalised global symmetries in holography: magnetohydrodynamic waves in a strongly interacting plasma,” JHEP 04 (2019) 141, arXiv:1707.04182 [hep-th].
  * (27) D. M. Hofman and N. Iqbal, “Generalized global symmetries and holography,” SciPost Phys. 4 no. 1, (2018) 005, arXiv:1707.08577 [hep-th].
  * (28) D. T. Son and A. O. Starinets, “Minkowski space correlators in AdS / CFT correspondence: Recipe and applications,” JHEP 09 (2002) 042, arXiv:hep-th/0205051.
  * (29) P. Kovtun, D. T. Son, and A. O. Starinets, “Holography and hydrodynamics: Diffusion on stretched horizons,” JHEP 10 (2003) 064, arXiv:hep-th/0309213.
  * (30) J. F. Fuini and L. G. Yaffe, “Far-from-equilibrium dynamics of a strongly coupled non-Abelian plasma with non-zero charge density or external magnetic field,” JHEP 07 (2015) 116, arXiv:1503.07148 [hep-th].
  * (31) S. Janiszewski and M. Kaminski, “Quasinormal modes of magnetic and electric black branes versus far from equilibrium anisotropic fluids,” Phys. Rev. D 93 no. 2, (2016) 025006, arXiv:1508.06993 [hep-th].
  * (32) E. D’Hoker and P. Kraus, “Magnetic Brane Solutions in AdS,” JHEP 10 (2009) 088, arXiv:0908.3875 [hep-th].
  * (33) S. Grozdanov, A. Lucas, and N. Poovuttikul, “Holography and hydrodynamics with weakly broken symmetries,” Phys. Rev. D 99 no. 8, (2019) 086012, arXiv:1810.10016 [hep-th].
  * (34) R. A. Davison and B. Goutéraux, “Momentum dissipation and effective theories of coherent and incoherent transport,” JHEP 01 (2015) 039, arXiv:1411.1062 [hep-th].
  * (35) C.-F. Chen and A. Lucas, “Origin of the Drude peak and of zero sound in probe brane holography,” Phys. Lett. B 774 (2017) 569–574, arXiv:1709.01520 [hep-th].
  * (36) R. A. Davison, S. A. Gentle, and B. Goutéraux, “Impact of irrelevant deformations on thermodynamics and transport in holographic quantum critical states,” Phys. Rev. D 100 no. 8, (2019) 086020, arXiv:1812.11060 [hep-th].
  * (37) E. Witten, “Multitrace operators, boundary conditions, and AdS / CFT correspondence,” arXiv:hep-th/0112258.
  * (38) M. Berkooz, A. Sever, and A. Shomer, “’Double trace’ deformations, boundary conditions and space-time singularities,” JHEP 05 (2002) 034, arXiv:hep-th/0112264.
  * (39) R. A. Davison and A. Parnachev, “Hydrodynamics of cold holographic matter,” JHEP 06 (2013) 100, arXiv:1303.6334 [hep-th].
  * (40) U. Moitra, S. K. Sake, and S. P. Trivedi, “Near-Extremal Fluid Mechanics,” arXiv:2005.00016 [hep-th].
  * (41) S. S. Gubser and A. Nellore, “Ground states of holographic superconductors,” Phys. Rev. D 80 (2009) 105007, arXiv:0908.1972 [hep-th].
  * (42) R. A. Davison, B. Goutéraux, and S. A. Hartnoll, “Incoherent transport in clean quantum critical metals,” JHEP 10 (2015) 112, arXiv:1507.07137 [hep-th].
  * (43) P. K. Kovtun and A. O. Starinets, “Quasinormal modes and holography,” Phys. Rev. D 72 (2005) 086009, arXiv:hep-th/0506184.
  * (44) E. D’Hoker, P. Kraus, and A. Shah, “RG Flow of Magnetic Brane Correlators,” JHEP 04 (2011) 039, arXiv:1012.5072 [hep-th].
  * (45) E. D’Hoker and P. Kraus, “Magnetic Field Induced Quantum Criticality via new Asymptotically AdS5 Solutions,” Class. Quant. Grav. 27 (2010) 215022, arXiv:1006.2573 [hep-th].
  * (46) P. Kovtun and A. Starinets, “Thermal spectral functions of strongly coupled N=4 supersymmetric Yang-Mills theory,” Phys. Rev. Lett. 96 (2006) 131601, arXiv:hep-th/0602059.
  * (47) M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.  Dover, New York, ninth dover printing, tenth gpo printing ed., 1964.
  * (48) D. Nickel and D. T. Son, “Deconstructing holographic liquids,” New J. Phys. 13 (2011) 075010, arXiv:1009.3094 [hep-th].
  * (49) P. Glorioso, M. Crossley, and H. Liu, “A prescription for holographic Schwinger-Keldysh contour in non-equilibrium systems,” arXiv:1812.08785 [hep-th].
  * (50) J. de Boer, M. P. Heller, and N. Pinzani-Fokeeva, “Holographic Schwinger-Keldysh effective field theories,” JHEP 05 (2019) 188, arXiv:1812.06093 [hep-th].
  * (51) S. Bhattacharyya, V. E. Hubeny, S. Minwalla, and M. Rangamani, “Nonlinear Fluid Dynamics from Gravity,” JHEP 02 (2008) 045, arXiv:0712.2456 [hep-th].
  * (52) N. Banerjee, J. Bhattacharya, S. Bhattacharyya, S. Dutta, R. Loganayagam, and P. Surowka, “Hydrodynamics from charged black branes,” JHEP 01 (2011) 094, arXiv:0809.2596 [hep-th].
  * (53) J. Erdmenger, M. Haack, M. Kaminski, and A. Yarom, “Fluid dynamics of R-charged black holes,” JHEP 01 (2009) 055, arXiv:0809.2488 [hep-th].
  * (54) S. de Haro, S. N. Solodukhin, and K. Skenderis, “Holographic reconstruction of space-time and renormalization in the AdS / CFT correspondence,” Commun. Math. Phys. 217 (2001) 595–622, arXiv:hep-th/0002230.
  * (55) C. Bender and S. Orszag, Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory.  Springer New York, 2013.


# Quantitative Evaluation of Explainable Graph Neural Networks for Molecular
Property Prediction

Jiahua Rao  Shuangjia Zheng  Yuedong Yang

###### Abstract

Advances in machine learning have led to graph neural network-based methods
for drug discovery, yielding promising results in molecular design, chemical
synthesis planning, and molecular property prediction. However, current graph
neural networks (GNNs) remain of limited acceptance in drug discovery is
limited due to their lack of interpretability. Although this major weakness
has been mitigated by the development of explainable artificial intelligence
(XAI) techniques, the “ground truth” assignment in most explainable tasks
ultimately rests with subjective judgments by humans so that the quality of
model interpretation is hard to evaluate in quantity. In this work, we first
build three levels of benchmark datasets to quantitatively assess the
interpretability of the state-of-the-art GNN models. Then we implemented
recent XAI methods in combination with different GNN algorithms to highlight
the benefits, limitations, and future opportunities for drug discovery. As a
result, GradInput and IG generally provide the best model interpretability for
GNNs, especially when combined with GraphNet and CMPNN. The integrated and
developed XAI package is fully open-sourced and can be used by practitioners
to train new models on other drug discovery tasks.

Machine Learning, ICML

##  1 Introduction

Various concepts of Graph Neural Network (GNN) have been successfully adopted
for drug discovery tasks in the past few years (Gilmer et al., 2017; Yang et
al., 2019; Liu et al., 2018; Jin et al., 2018; Yan et al., 2020) such as
molecular property prediction and de novo molecular generation. Several graph
neural network models have been shown to yield more promising results than the
existing machine learning and quantitative structure-activity relationship
(QSAR) methods for drug discovery (Xiong et al., 2019; Song et al., 2020).This
advance is mostly owed to the ability of graph neural networks to effectively
model molecular graph data. With the increasing high-quality labeled data,
graph neural networks are promising to learn better molecular representations
that finally replace decades-old hand-crafted molecular fingerprint
representations.

Figure 1: Schematic of DrugXAI framework, including setups of explainable
tasks and attribution metrics. We built three levels of benchmark datasets to
quantitatively assess the interpretability of GNN models. Based on the
benchmarks, we implemented state-of-the-art XAI and commonly used GNN models,
and provided a uniform and rigorous framework to evaluate their performances.

Despite their promise, graph neural networks remain of limited acceptance in
drug discovery partly due to their lack of interpretability, and these models
are often considered as “black-box” (Jiménez-Luna et al., 2020, 2021). While
there have been efforts to model interpretability based on simplifications of
models, feature sub-selection, or attention (Ribeiro et al., 2016;
Sundararajan et al., 2017; Shrikumar et al., 2017; Vaswani et al., 2017), this
problem is further exacerbated by the fact that these models often produce the
correct answers for wrong reasons (Lapuschkin et al., 2019). Given the current
pace of AI in drug discovery, there will be an increased demand for
interpretability methods that help us understand and interpret the GNN models.
Thus, a highly accurate and mechanically interpretable model may be the key to
accelerated drug discovery with graph neural networks.

Explainable artificial intelligence (XAI) aims to help scientists to know how
the model reached a particular answer and explain why the answer provided by
the model is acceptable (Guidotti et al., 2018; Ying et al., 2019; Lundberg et
al., 2020). In drug discovery-related applications, in particular for property
prediction tasks, XAI methods could help the development of graph neural
networks by quantifying the molecular substructures that are critical for a
given prediction and explaining how reliable a prediction is (Jiménez-Luna et
al., 2021; Yu et al., 2021). Feature attribution is one approach to
interpretability, which measures how important the feature is to the model’s
prediction of a target property. Attribution methods have been studied in the
domains of drug discovery, for example, McCloskey et al. (2019) developed an
attribution method to verify whether each model trained on protein–ligand
binding data learns its corresponding binding logic correctly. Jiménez-Luna et
al. (2021) established an Integrated Gradients (IG) feature attribution
technique to examine the interpretability of GNN models in molecular property
predictions. Furthermore, with the development of graph neural networks, it
does not come as surprise in the fields of drug discovery to explore the
interpretability of models trained with graph convolution architecture. For
instance, Jin et al. (2020) employed Monte Carlo tree search to extract
molecular substructures with the help of a property predictor, which are
likely responsible for each property of interest. Yu et al. (2021) proposed a
framework of Graph Information Bottleneck (GIB) for the subgraph recognition,
which could recognize a compressed subgraph with minimum information loss in
terms of predicting the molecular properties. These methods, which can be
considered as subgraph recognition methods, are interested in directly
identifying the substructures that mostly represent certain properties of
molecules.

As mentioned above, many efforts have been made to mitigate the major weakness
of deep learning approaches, that is the lack of causal understanding.
Unfortunately, the quality and evaluation of model interpretations is hard to
examine because acquiring ground truth sub-structures or attributions requires
expensive wet experiment and subjective expert judgment. Although Sanchez-
Lengeling et al. (2020) have built an open-source synthetic benchmarking suite
for attribution methods on GNNs, the synthetic tasks were designed for
recognize simple subgraphs such as benzene from molecules. In fact, effective
XAI should reveal more complex information to scientists and render the
decision-making process. For example, in toxicity prediction, there are often
dozens of fragments that cause a particular toxicity, and some even have to be
considered many-to-one scenarios. Moreover, there is a kind of molecular pairs
with similar structures but completely different properties, which is known as
property cliffs (Stumpfe et al., 2014). Such realistic and sophisticated
scenarios will bring more challenges and opportunities to XAI.

Table 1: Statisitics of the datasets Tasks | Type | Dataset Name | Compounds | Train | Test | Subgraph ground truth
---|---|---|---|---|---|---
Single-rationale | Graph classification | 3MR | 3152 | 2521 | 631 | Three-membered ring
Graph classification | Benzene | 12000 | 9600 | 2400 | Benzene ring
Multiple-rationales | Graph classification | Mutagenicity | 6506 | 5204 | 1302 | Mutagenicity alerts
Graph classification | Liver | 587 | 469 | 118 | Hepatotoxic alerts
Property cliff | Graph regression | hERG | 6993 | 6483 | 510 | Structural motifs
Graph classification | CYP450 | 9122 | 9025 | 97 | Structural motifs

In this work, we established three levels of benchmarks, from easy to hard,
for quantitatively and comprehensively assessing the interpretability of
current graph neural networks and XAI modules, highlighting theirs benefits,
limitations and future opportunities for drug discovery. In addition, we
provide a uniform and rigorous framework to evaluate the performance of
commonly-used XAI methods on several GNN model types. We examine the quality
of those XAI methods through accuracy and consistency, which evaluate the
ability of XAI methods in identifying the subgraphs that mostly represent
certain properties of the molecules and perturbations that are most likely to
cause property cliffs. We offer three main contributions:

  * •

We build three levels of benchmark datasets to quantitatively assess the
interpretability of the state-of-the-art GNN models.

  * •

We provide a uniform and rigorous framework to evaluate the performance of
commonly-used XAI methods on several GNN model types.

  * •

We provide a comprehensive overview of recent XAI methods, highlighting their
benefits, limitations and future opportunities for drug discovery.

##  2 Method

In this section, We will first define the tasks and datasets used in DrugXAI
and provide necessary information analyzing its effectiveness in evaluating
the model interpretations. For the model interpretations, we will introduce
commonly-used XAI methods on several GNN model types. Finally, detailed
information about training and evaluating the XAI methods in drug discovery
applications will be presented in the end.

Table 2: Prediction Performance on benchmarks | Dataset | 3MR | Benzene | Liver | Mutagenicity | hERG | CYP450
---|---|---|---|---|---|---|---
| Metric | AUROC | AUROC | ACC | AUROC | RMSE | AUROC
GNN Model | CMPNN | 0.998 | 0.999 | 0.432 | 0.866 | 1.328 | 0.618
GraphSAGE | 0.989 | 0.999 | 0.415 | 0.843 | 1.230 | 0.586
GraphNet | 0.995 | 0.998 | 0.423 | 0.823 | 1.304 | 0.600
GAT | 0.995 | 0.989 | 0.381 | 0.852 | 1.361 | 0.564
Figure 2: XAI method accuracy across Single-rationale tasks and model
architectures.  Colors are used to distinguish two metric types — Attribution
AUROC and Attribution ACC. With the perfect predictive performances, the
existing popular XAI methods on the trained GNN models are consistently
superior to the random baseline and the improvements are significant.

###  2.1 Tasks and Benchmark Datasets

Our goal is to evaluate whether current XAI methods allow GNNs to perform
graph classification or regression tasks while identifying task-relevant
molecular substructures (i.e., rationales). The substructures can be a
particular subgraph, the conjunction of multiple subgraphs or a local
transformation between two molecular graphs (i.e, property cliff).

Correspondingly, we design three levels (single rationale, multiple
rationales, and property cliff) of benchmark datasets to quantitatively assess
the interpretability of the state-of-the-art GNN models. We summarized the
statistics of the datasets in Table 1.

####  2.1.1 Single-rationale tasks

Following the strategy of McCloskey et al. (2019), we first established two
single rationale benchmarks, whose goal is to identify if a molecular graph
contains particular subgraphs of interest. We defined two synthetic subgraph
logics: Benzene and 3MR, aiming to verify whether each model is able to
capture the corresponding subgraphs correctly. In particular, we collected 50K
unlabeled molecules from ZINC15 lead-like subset (Sterling & Irwin, 2015) and
used substructure matching (Benzene and 3MR) to identify positive molecules
that matched, and then randomly selected the same amount of unmatched
molecules as negative samples. Both datasets are graph classification tasks, 1
for those containing a specific substructure, 0 otherwise.

####  2.1.2 Multiple-rationales tasks

In real-world scenario, there are often dozens of rationales that cause a
particular property. To this end, we construct two multiple-rationales
benchmarks in which ground truths are the hand-crafted structural alerts that
raise the problems of hepatotoxicity and ames mutagenicity. In particular, for
hepatotoxicity, we collected the data from Liu et al. (2015), which contains
174 hepatotoxic, 230 possible hepatotoxic and 183 non-hepatotoxic compounds,
and 12 molecular rationales have been identified as structural alerts for
human liver injuries. For ames mutagenicity, we used a well-established
dataset compiled by Hansen et al. (2009), including 6506 compounds and
corresponding ames mutagenicity results. Sushko et al. (2012) summarizes 46
structural toxic alerts that raise the mutagenicity, which have been computed
as ground truth subgraphs.

####  2.1.3 Property cliff tasks

To further explore the interpretability of the commonly-used XAI methods in
extreme cases, we compiled two publicly available datasets, hERG inhibition
and Cytochrome P450 Inhibition, for property cliff tasks. We use the dataset
constructed by Jiménez-Luna et al. (2021), which collected by a literature
survey and publicly available data. Specifically, hERG inhibition is an
important anti-target that must be avoided during drug development. For this
endpoint, 6993 compounds with reported activity (IC50 values) in the nanomolar
range were collected and IC50 values were transformed into the pIC50 scale
during model training (Sato et al., 2018). Inhibition of cytochrome P450
(CYP450) enzymes is the most common mechanism leading to drug–drug
interactions. For this endpoint, 9120 CYP3A4 inhibitors and substrates with
binary activity information (active/inactive) were determined by Veith et al.
(2009).

Previously, Jiménez-Luna et al. (2021) established these endpoint datasets but
missed ground truths for interpretations. The “ground truth” in their
experiments ultimately rests with subjective human judgment. Differently, we
have constructed a quantitatively assessable scheme based on the
characteristics of the property cliff. Particularly, we first split these
datasets into training and testing, forcing those eligible property cliff
molecular pairs to appear only in the test set. Then, the property cliff pairs
of molecules with similar structure but completely different properties were
captured by MMPA technique (Hussain & Rea, 2010) in the RDKit package and the
ground truth attributions were computed according to the local differences
between the molecular pairs. For the property cliff dataset which is the graph
regression task, we define a local transformation to create a pair of active
cliff pairs when it causes more than a 100-fold change ($\bigtriangleup$ pIC50
$\textgreater$ 2) in properties.

###  2.2 Graph Neural Network

When viewing molecules as graphs with atoms as nodes and chemical bonds as
edges, Graph Neural Network (GNN) (Kipf & Welling, 2017; Gilmer et al., 2017)
is a neural network that takes a graph as input and outputs a graph with the
same topology, but with updated node, edge and/or graph-level information. One
key feature of the GNNs we study is the message passing function, which allows
nodes to update their states by aggregating feature information from
neighboring nodes and edges. Depending on the message passing strategy, the
message can contain information about the node, edge, or the global
information.

Our experiments use four existing popular GNN architectures distinguished by
their message-passing strategies. The first is a framework for inductive
representation learning (GraphSAGE)(Hamilton et al., 2017), which is used to
generate low-dimensional vector representations for nodes. The second model is
the Communicative Message Passing Neural Network (CMPNN) (Song et al., 2020)
which improve the molecular embedding by strengthening the message
interactions between nodes and edges through a communicative kernel. Another
model is the Graph Attention Network (GAT)(Veličković et al., 2018) which
aggregates node information via an attention mechanism. The fourth is Graph
Nets(Battaglia et al., 2018) in which the message mechanism utilizes a global
features vector in addition to node and edge features.

Figure 3: XAI method accuracy across Multiple-rationales tasks and model
architectures.  Colors are used to distinguish two metric types — Attribution
AUROC and Attribution ACC. GradInput and IG consistently well across these
tasks and GNN models. The attribution performances, for example, the
Attribution-ACC on Liver datasets ranging between 0.835 and 0.938, have shown
that the GNN models with XAI methods were able to detect the implicit
relationships compared to the random baselines.

###  2.3 XAI methods

An XAI method $A$ takes a model $M$ and a molecular graph $G$ as inputs to
generate an attribution score $G_{A}=(v_{A},e_{A})$ where $v_{A}$ and $e_{A}$
are node and edge weightings relevant for predicting property $y$. These
weightings can be visualized as a heatmap superimposed on a graph. Our ground-
truths for attributions are node-level, so we redistribute edge attributions
equally onto their endpoint nodes’ attributions. In our framework, we utilize
the following methods for molecular graphs: GradInput (Shrikumar et al.,
2017), GradCAM (Selvaraju et al., 2017), Integrated Gradients (Sundararajan et
al., 2017), CAM (Zhou et al., 2016), MCTS (Jin et al., 2020) and Attention
weights (Vaswani et al., 2017).

Class Activation Map (CAM) is a feature attribution method, which uses a
global average pooling (GAP) layer prior to class outputs to obtain
attributions. CAM attributions could be derived by multiplying the final
convolutional layer’s feature map activations with the output weights at the
last message passing layer. GradInput attributions are calculated by the
element-wise product of the input graph with the gradient of property with
respect to the input node and edge features while GradCAM extends GradInput by
using intermediate activations in which corresponds to the element-wise
product of the activations of intermediate message-passing layer with the
gradient of property. Integrated Gradients (IG) integrates the element-wise
product of an interpolated input with the gradient of property. Attention
mechanism method is specific to the GAT model, which will produce attention
scores on edges to adjacent nodes. We can use these attention scores as a
measure of importance for propagating information relevant to the predictive
task. Monte Carlo Tree Search (MCTS), a subgraph recognition method, is also
used to extract candidate rationales from molecules with the help of a
property predictor.

##  3 Experiments

In this section, We first present the predictive performance evaluation for
graph classification and regression tasks on different benchmarks. It’s
necessary as the interpretations generated by the model with poor predictive
capability bear little trust. Next, we evaluate the attribution performance of
commonly-used XAI methods on the trained GNN models. Finally, we analyze the
performances on the three levels of benchmarks, from easy to hard, from
recognizing particular subgraphs of interest to identifying property cliffs.
All the GNN models and XAI methods were implemented based on Rao et al. (2021)
and Sanchez-Lengeling et al. (2020). Particularly, no explicit hyper-parameter
optimization is performed on any model training to avoid unnecessary model
selection bias. All the implementation details can be found in
https://github.com/biomed-AI/MolRep.

###  3.1 Predictive Performance

The results of predictive performance on our benchmarks are presented in Table
2, where the root-mean-square error (RMSE) of regression tasks, the receiver-
operator characteristic area under the curve (AUROC) of classification tasks
and the accuracy (ACC) of multi-classification tasks are reported. We observed
that all trained models showed perfect predictive capabilities on the single-
rationale datasets, with AUROC ranging between 0.989 and 0.999, as shown in
Table 2. These values suggested that the meaningful molecular graph features
were identified by all trained models in the learning process. On the
multiple-rationales benchmarks, the results obtained were also markedly better
than random, with ACC ranging between 0.381 and 0.432 for the three-
classification dataset, and AUROC ranging between 0.823 and 0.866 for the
binary classification model. Finally, on the most difficult datasets, the
poorest performances of the existing GNN models were obtained, suggesting that
current GNN models have difficulties to capture the property-relevant features
during the learning process.

###  3.2 Single Rationale Recognition

On the single-rationale datasets, we used Attribution-AUC and Attribution-ACC,
designed by McCloskey et al. (2019), to assess attribution accuracy and
consistency in identifying the corresponding substructure correctly. If
multiple attributions are valid (e.g., a subgraph is present twice in a
graph), we take the maximum attribution value of all possible solutions.
Particularly, the subgraph recognition methods such as MCTS could not be
evaluated by Attribution-AUC since their outputs are directly subgraphs rather
than a node probability.

Jointly analyzing Table 2 and figure 2, we found that with the perfect
predictive performances, the existing popular XAI methods on the trained GNN
models are consistently superior to the random baseline and the improvements
are significant, suggesting that the trained models have identified the
meaningful features that are relevant to molecular property and the current
XAI approaches could effectively recognize subgraph logics corresponding. This
makes sense since the graph labels of single-rationale datasets are whether
the corresponding molecule contains the particular subgraphs of interest. The
direct relationship between target rationales and graph labels helps the
models to capture accurate information.

Figure 4: XAI method accuracy across Property-Cliff tasks and model
architectures.  Colors are used to distinguish two metric types — Attribution
F1 and Attribution ACC. Most of the XAI methods failed to recognize property
cliffs on the Property-Cliff tasks. Among them, IG and GradInput performed
relatively well.

###  3.3 Multiple Rationales Recognition

Furthermore, in real-world scenario, the relationships between the molecular
structures and their related properties are often implicit. Therefore, whether
the current GNN model with the popular XAI methods is capable of capturing the
implicit relationships is the motivation for us to conduct the attribution
experiments on multiple-rationales datasets.

As shown in Figure 3 and Table 2, owing to the prediction performances have
demonstrated that the useful graph features were identified during the
learning process, the attribution performances, for example, the Attribution-
ACC on Liver datasets ranging between 0.835 and 0.938, have shown that the GNN
models with XAI methods were able to detect the implicit relationships
compared to the random baselines. More specifically, we could observe that
GraphNet and CMPNN with high predictive performance have better attribution
performance than other baselines across the popular XAI methods. And GradInput
and IG across tasks and model architectures have better attribution
performances. We attribute the improvements to the utilization of the
gradients of the predicted property, which directly provide detailed insight
into the relationship between the input features and the property.

###  3.4 Property Cliff Identification

Finally, to explore the interpretability of GNN models on the more complex
relationships between substructures and graph properties, we evaluated the
capabilities of the models to recognize property cliffs. Since the ground
truths (local transformations) on these datasets have been classified to
positive and negative based on the change of property, we only report
Attribution-ACC and Attribution-F1 to assess the attribution accuracy and
verify whether each model learns the corresponding key substructure
transformations relevant for predicting property correctly.

Jointly Analyzing Figure 4 and Table 2, we observed that most of the XAI
methods failed to recognize property cliffs effectively under this condition,
with low Attribution-F1 ranged from 0.104 to 0.294 on hERG and CYP datasets.
There are two possible reasons for the results. Firstly, the poor predictive
performances of the existing GNN models were obtained, suggesting that current
GNN models were unable to capture the property-relevant features on the
property-cliff tasks. Secondly, the relationships between property cliffs and
their corresponding property are extremely diverse and complex. But considered
from another perspective, on property cliff benchmarks, GNN methods are ten or
even tens of times better than the random baseline, proving that although the
problem is very difficult, the GNNs have captured some key features from it.

##  4 Discussion and conclusions

In this work we have established a series of XAI benchmarks in the context of
drug discovery and created a framework to quantitatively evaluate attribution
methods in GNNs. We expect that our benchmark and framework will aid in
developing better methodology for interpretability in molecular graph tasks.
Overall, our experiments show existing popular XAI methods have good
performance with GNN architectures in recognizing regular subgraphs in
molecules but performed poorly in identifying the property cliff on realistic
datasets. As a final caveat, we find that XAI techniques that are directly
related to predicted labels like GradInput and IG also tend to have better
attribution performance, across many models. This also hints that current XAI
methods cannot yet be used as a recipe for drug discovery, requiring
significant human expertise for correct interpretation. The reasons for these
phenomena remain a topic of further study.

There also remains much room for improvement in interpretation performance for
GNNs. Nonetheless, further development of XAI applications would greatly
benefit from this meaningful benchmarking. We hope that our quantitative
benchmark and framework of interpretations in GNNs are useful for developing
new XAI methods. Our framework and dataset are available for further study,
broadening the development of explainable artificial intelligence in drug
discovery-related applications.

## References

  * Battaglia et al. (2018) Battaglia, P. W., Hamrick, J. B., Bapst, V., Sanchez-Gonzalez, A., Zambaldi, V., Malinowski, M., Tacchetti, A., Raposo, D., Santoro, A., Faulkner, R., et al.  Relational inductive biases, deep learning, and graph networks.  _arXiv preprint arXiv:1806.01261_ , 2018.
  * Gilmer et al. (2017) Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E.  Neural message passing for quantum chemistry.  In _International Conference on Machine Learning_ , pp. 1263–1272. PMLR, 2017.
  * Guidotti et al. (2018) Guidotti, R., Monreale, A., Ruggieri, S., Turini, F., Giannotti, F., and Pedreschi, D.  A survey of methods for explaining black box models.  _ACM computing surveys (CSUR)_ , 51(5):1–42, 2018\.
  * Hamilton et al. (2017) Hamilton, W. L., Ying, R., and Leskovec, J.  Inductive representation learning on large graphs.  _arXiv preprint arXiv:1706.02216_ , 2017.
  * Hansen et al. (2009) Hansen, K., Mika, S., Schroeter, T., Sutter, A., Ter Laak, A., Steger-Hartmann, T., Heinrich, N., and Müller, K.-R.  Benchmark data set for in silico prediction of ames mutagenicity.  _Journal of chemical information and modeling_ , 49(9):2077–2081, 2009.
  * Hussain & Rea (2010) Hussain, J. and Rea, C.  Computationally efficient algorithm to identify matched molecular pairs (mmps) in large data sets.  _Journal of chemical information and modeling_ , 50(3):339–348, 2010.
  * Jiménez-Luna et al. (2020) Jiménez-Luna, J., Grisoni, F., and Schneider, G.  Drug discovery with explainable artificial intelligence.  _Nature Machine Intelligence_ , 2(10):573–584, 2020.
  * Jiménez-Luna et al. (2021) Jiménez-Luna, J., Skalic, M., Weskamp, N., and Schneider, G.  Coloring molecules with explainable artificial intelligence for preclinical relevance assessment.  _Journal of Chemical Information and Modeling_ , 61(3):1083–1094, 2021.
  * Jin et al. (2018) Jin, W., Barzilay, R., and Jaakkola, T.  Junction tree variational autoencoder for molecular graph generation.  In _International Conference on Machine Learning_ , pp. 2323–2332. PMLR, 2018.
  * Jin et al. (2020) Jin, W., Barzilay, R., and Jaakkola, T.  Multi-objective molecule generation using interpretable substructures.  In _International Conference on Machine Learning_ , pp. 4849–4859. PMLR, 2020.
  * Kipf & Welling (2017) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  In _International Conference on Learning Representations (ICLR)_ , 2017.
  * Langley (2000) Langley, P.  Crafting papers on machine learning.  In Langley, P. (ed.), _Proceedings of the 17th International Conference on Machine Learning (ICML 2000)_ , pp. 1207–1216, Stanford, CA, 2000\. Morgan Kaufmann.
  * Lapuschkin et al. (2019) Lapuschkin, S., Wäldchen, S., Binder, A., Montavon, G., Samek, W., and Müller, K.-R.  Unmasking clever hans predictors and assessing what machines really learn.  _Nature communications_ , 10(1):1–8, 2019.
  * Liu et al. (2018) Liu, Q., Allamanis, M., Brockschmidt, M., and Gaunt, A. L.  Constrained graph variational autoencoders for molecule design.  _arXiv preprint arXiv:1805.09076_ , 2018.
  * Liu et al. (2015) Liu, R., Yu, X., and Wallqvist, A.  Data-driven identification of structural alerts for mitigating the risk of drug-induced human liver injuries.  _Journal of cheminformatics_ , 7(1):1–8, 2015\.
  * Lundberg et al. (2020) Lundberg, S. M., Erion, G., Chen, H., DeGrave, A., Prutkin, J. M., Nair, B., Katz, R., Himmelfarb, J., Bansal, N., and Lee, S.-I.  From local explanations to global understanding with explainable ai for trees.  _Nature machine intelligence_ , 2(1):56–67, 2020\.
  * McCloskey et al. (2019) McCloskey, K., Taly, A., Monti, F., Brenner, M. P., and Colwell, L. J.  Using attribution to decode binding mechanism in neural network models for chemistry.  _Proceedings of the National Academy of Sciences_ , 116(24):11624–11629, 2019.
  * Rao et al. (2021) Rao, J., Zheng, S., Song, Y., Chen, J., Li, C., Xie, J., Yang, H., Chen, H., and Yang, Y.  Molrep: A deep representation learning library for molecular property prediction.  _bioRxiv_ , 2021.
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., and Guestrin, C.  ” why should i trust you?” explaining the predictions of any classifier.  In _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_ , pp. 1135–1144, 2016.
  * Sanchez-Lengeling et al. (2020) Sanchez-Lengeling, B., Wei, J., Lee, B., Reif, E., Wang, P., Qian, W. W., McCloskey, K., Colwell, L., and Wiltschko, A.  Evaluating attribution for graph neural networks.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Sato et al. (2018) Sato, T., Yuki, H., Ogura, K., and Honma, T.  Construction of an integrated database for herg blocking small molecules.  _PLoS One_ , 13(7):e0199348, 2018.
  * Selvaraju et al. (2017) Selvaraju, R. R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., and Batra, D.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In _Proceedings of the IEEE international conference on computer vision_ , pp. 618–626, 2017.
  * Shrikumar et al. (2017) Shrikumar, A., Greenside, P., and Kundaje, A.  Learning important features through propagating activation differences.  In _International Conference on Machine Learning_ , pp. 3145–3153. PMLR, 2017.
  * Song et al. (2020) Song, Y., Zheng, S., Niu, Z., Fu, Z.-H., Lu, Y., and Yang, Y.  Communicative representation learning on attributed molecular graphs.  In _IJCAI_ , 2020.
  * Sterling & Irwin (2015) Sterling, T. and Irwin, J. J.  Zinc 15–ligand discovery for everyone.  _Journal of chemical information and modeling_ , 55(11):2324–2337, 2015.
  * Stumpfe et al. (2014) Stumpfe, D., Hu, Y., Dimova, D., and Bajorath, J.  Recent progress in understanding activity cliffs and their utility in medicinal chemistry: miniperspective.  _Journal of medicinal chemistry_ , 57(1):18–28, 2014.
  * Sundararajan et al. (2017) Sundararajan, M., Taly, A., and Yan, Q.  Axiomatic attribution for deep networks.  In _International Conference on Machine Learning_ , pp. 3319–3328. PMLR, 2017.
  * Sushko et al. (2012) Sushko, I., Salmina, E., Potemkin, V. A., Poda, G., and Tetko, I. V.  Toxalerts: a web server of structural alerts for toxic chemicals and compounds with potential adverse reactions, 2012.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I.  Attention is all you need.  _arXiv preprint arXiv:1706.03762_ , 2017.
  * Veith et al. (2009) Veith, H., Southall, N., Huang, R., James, T., Fayne, D., Artemenko, N., Shen, M., Inglese, J., Austin, C. P., Lloyd, D. G., et al.  Comprehensive characterization of cytochrome p450 isozyme selectivity across chemical libraries.  _Nature biotechnology_ , 27(11):1050–1055, 2009\.
  * Veličković et al. (2018) Veličković, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., and Bengio, Y.  Graph Attention Networks.  _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=rJXMpikCZ.  accepted as poster.
  * Xiong et al. (2019) Xiong, Z., Wang, D., Liu, X., Zhong, F., Wan, X., Li, X., Li, Z., Luo, X., Chen, K., Jiang, H., et al.  Pushing the boundaries of molecular representation for drug discovery with the graph attention mechanism.  _Journal of medicinal chemistry_ , 63(16):8749–8760, 2019.
  * Yan et al. (2020) Yan, C., Ding, Q., Zhao, P., Zheng, S., Yang, J., Yu, Y., and Huang, J.  Retroxpert: Decompose retrosynthesis prediction like a chemist.  _arXiv preprint arXiv:2011.02893_ , 2020.
  * Yang et al. (2019) Yang, K., Swanson, K., Jin, W., Coley, C., Eiden, P., Gao, H., Guzman-Perez, A., Hopper, T., Kelley, B., Mathea, M., et al.  Analyzing learned molecular representations for property prediction.  _Journal of chemical information and modeling_ , 59(8):3370–3388, 2019.
  * Ying et al. (2019) Ying, R., Bourgeois, D., You, J., Zitnik, M., and Leskovec, J.  Gnnexplainer: Generating explanations for graph neural networks.  _Advances in neural information processing systems_ , 32:9240, 2019.
  * Yu et al. (2021) Yu, J., Xu, T., Rong, Y., Bian, Y., Huang, J., and He, R.  Graph information bottleneck for subgraph recognition.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=bM4Iqfg8M2k.
  * Zhou et al. (2016) Zhou, B., Khosla, A., Lapedriza, A., Oliva, A., and Torralba, A.  Learning deep features for discriminative localization.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 2921–2929, 2016.