# Constraining Electromagnetic Signals from Black Holes with Hair

Nicole R. Crumpler*  —<EMAIL_ADDRESS>William H. Miller III Department of
Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21210, USA

###### Abstract

We constrain a broad class of “hairy” black hole models capable of directly
sourcing electromagnetic radiation during a binary black hole merger. This
signal is generic and model-independent since it is characterized by the black
hole mass ($M$) and the fraction of that mass released as radiation
($\epsilon$). For field energy densities surpassing the Schwinger limit, this
mechanism triggers pair-production to produce a gamma-ray burst. By cross-
referencing gravitational wave events with gamma-ray observations, we place
upper bounds of $\epsilon<10^{-5}-10^{-4}$ for $10-50$ $M_{\odot}$ black holes
depending on the black hole mass. We discuss the weak detection of a gamma-ray
burst following GW150914 and show that this event is consistent with rapid
electromagnetic emission directly from a “hairy” black hole with $\epsilon\sim
10^{-7}-10^{-6}$. Below the Schwinger limit, ambient charged particles are
rapidly accelerated to nearly the speed of light by the strong electromagnetic
field. For 1-50 $M_{\odot}$ black holes and $\epsilon$ ranging from $10^{-20}$
to $10^{-7}$, the typical proton energies are $\sim 20$ GeV-20 TeV and
electron energies are $\sim 0.01-10$ GeV. At these energies, cosmic ray
protons and electrons quickly diffuse into the Milky Way’s background magnetic
field, making it difficult to identify a point source producing them. Overall,
constraining $\epsilon$ in this less energetic regime becomes difficult and
future constraints may need to consider specific models of “hairy” black
holes.

††preprint: APS/123-QED

##  I Introduction

There is an important distinction between astrophysical and mathematical black
holes (BHs). Astrophysically, BHs are observed as compact regions of spacetime
in which gravity is so strong that even light cannot escape. These objects
have been detected merging with each other [1], emitting electromagnetically-
bright jets [2], consuming stars [3], and more. Mathematically, BHs are vacuum
solutions of the Einstein equations of general relativity describing spacetime
external to a compact mass distribution within an event horizon. The theory
underlying mathematical BHs has been used to characterize astrophysical BHs,
although there remain clear disconnects between these theoretical models and
physical reality.

One such disconnect stems from the so-called “no-hair” theorem. Canonically,
the “no-hair” theorem proposes that mathematical BHs are completely
characterized by the BH’s mass, charge, and spin as seen by an external
observer. This theorem has been tested astrophysically using a variety of
probes such as radio observations of the shadow of a BH event horizon [4, 5,
6], gravitational wave signals of binary BH (BBH) mergers [7, 8], and stellar
orbits around the galactic center [9, 10]. No evidence for its violation has
yet been discovered in astrophysical BHs. However, the “no-hair” theorem leads
directly to the BH information paradox [11]. In this paradox, different
configurations of matter, radiation, etc. that have fallen into a BH can be
described by the same mathematical BH solution, losing information about the
initial quantum state of the system. This violates a core tenet of quantum
mechanics in a regime in which both the theories of general relativity and
quantum mechanics are valid. There have been many attempts to resolve this
paradox; one can refer to Ref. [12] for a recent review. Despite these
attempts, no consensus has yet been reached.

The BH information paradox supports a possibility of richer physics underlying
BHs. “Hairy” BHs are novel solutions to the Einstein field equations which are
characterized by more than the three parameters of a canonical BH. Some of
these BH models have been proposed as explicit solutions to the BH information
paradox. One interesting possibility, proposed by Ref. [13], is the firewall
BH model. This mathematical BH has a singular shell (otherwise known as a
firewall) just outside the horizon, causing general relativity to break down
outside the BH horizon. Since general relativity no longer holds in this
regime, the BH information problem no longer applies. Such an exotic object
appears as a Schwarzschild BH to a distant observer, raising the question of
how a firewall BH (and, more generally, other “hairy” BHs) might be
distinguished from a canonical BH in astrophysical observations.

Electromagnetic (EM) radiation from astrophysical BHs in baryon-poor
environments would be a beacon of new fundamental physics and support the
existence of non-canonical BHs. Canonical BHs do not directly source EM
radiation, except the weak emission of thermal Hawking radiation [14], which
is not observable for BHs of astrophysically relevant masses. However, there
are convincing qualitative arguments indicating that “hairy” BH models could
radiate an appreciable proportion of their mass as EM radiation. As a
motivating example, consider the firewall BH discussed earlier. Ref. [13]
suggests several ways in which this model could produce EM radiation including
the explosion of an unstable firewall, a BH phase transition from a canonical
to a firewall BH, and BBH mergers involving a firewall BH. We emphasize that
no quantitative model for this effect yet exists, since working out the
details of such a model would likely be a challenging exercise in quantum
gravity. We hope this paper will convince those working in the field that such
a pursuit is worthwhile, given the fact that these models generically produce
observable signals and thus can be constrained via astrophysical observations.
Overall, our understanding of BHs is inconsistent, motivating searches for
generic signals of deviations from canonical BH models.

Given that only non-canonical BHs can radiate appreciably, there are two
important considerations in order to distinguish this radiation from typical
astrophysical sources. Firstly, we must be confident that there is a BH in the
region sourcing the radiation. Secondly, the BH must be in a sufficiently
baryon-poor environment and the emitted radiation must be sufficiently
energetic that we can be confident the radiation is not produced by standard
processes such as relativistic jets. The best observable available satisfying
these considerations is concurrent observations of BBH mergers with
gravitational wave detectors and EM radiation instruments. Gravitational wave
detectors such as the Advanced Laser Interferometer Gravitational-Wave
Observatory (LIGO, [15]) and Virgo [16] regularly observe the mergers of
stellar mass BHs. These events have only been observed extragalactically, with
local BBH merger rates measured to be $\sim 10\text{ Gpc}^{-3}\text{ yr}^{-1}$
[17]. Thus, any observable EM radiation from such events must be extremely
energetic, on the order of supernova energies.

There is typically insufficient baryonic matter surrounding BBHs to produce EM
radiation observable at these extragalactic distances and energy scales. Some
models have been proposed to produce a gamma-ray burst (GRB) during a BBH
merger. Most of these models require rapid accretion during the merger [18,
19, 20], necessitating a baryon-rich environment, and all of these models have
been contested [21, 22, 23]. Another class of model involves charged BHs [24],
but the required charge is unreasonably large [25]. Given the shortcomings of
these models, no extragalactically-observable EM signal is expected from a
stellar-mass BBH merger. Thus, in this paper, we investigate the observational
signatures of a stellar-mass BH directly releasing some of its mass as EM
radiation during a BBH merger as a novel indicator of the existence of “hairy”
BHs.

Since BBHs mergers are catastrophic events characterized by short timescales,
a large amount of energy could be emitted by the BH in a burst of EM
radiation. We parameterize the amount of energy released by $\epsilon$ such
that $\epsilon M$ is the energy emitted directly by the BH in EM radiation
into the ambient environment. The characteristic frequency of the EM signal is
independent of the “hairy” BH model sourcing the radiation because any such
emission directly from the BH must tunnel out of the gravitational well in the
same manner as Hawking radiation [26]. If the Schwarzschild radius of a BH is
$r_{s}$, then the characteristic frequency of the radiation emitted by the BH
is

$f=\frac{1}{2r_{s}}=3.3\times
10^{-17}\left(\frac{M}{M_{\odot}}\right)^{-1}\text{ MeV}.$ (1)

Throughout this paper, we work in natural units where
$\hbar=c=k_{B}=\epsilon_{0}=1$ unless otherwise stated. The frequencies
emitted by stellar mass BHs are very low-frequency radio waves. At these
frequencies, all of the radiation is absorbed in the interstellar medium,
predominantly via free-free absorption by the warm ionized medium [27]. Thus,
this radiation is not directly observable, but is absorbed and re-emitted as a
secondary signal. We calculate the range of $M$ and $\epsilon$ for which this
secondary signal could be detected. In future work, model-dependent effects
will need to be included to augment this generic parameterization. Again we
emphasize that, although there are “hairy” BH models capable of producing EM
radiation, no complete model able to make quantitative predictions for this
effect exists currently. We hope this paper will motivate others in the field
to work through the details of such models.

In this paper, we constrain a broad class of “hairy” BH models using a generic
and model-independent EM signal that is characterized by the BH mass ($M$) and
the fraction of that mass that is lost to EM radiation ($\epsilon$). In
Section II, we discuss the two phenomenologically distinct cases in which
radiation is emitted and derive the critical value of $\epsilon$ that
separates them. This division is set by the Schwinger limit, above which the
BH radiation triggers pair production resulting in a GRB and below which the
EM field accelerates ambient charged particles to create an overdensity of
cosmic rays. In Section III we characterize the extragalactic observability of
a GRB created by the BH radiation to constrain $\epsilon$ given the non-
detection of GRBs from BBH mergers. In Section IV we describe the electron and
proton cosmic ray energy spectrum created below the Schwinger limit and
discuss the difficulties of observationally constraining $\epsilon$ in this
less-energetic regime. In Section V we summarize our results.

##  II The Schwinger Limit

The Schwinger limit dictates the critical value of $\epsilon$ separating the
two phenomenologically distinct cases in which radiation is emitted for a
particular BH mass. This limit, derived from quantum electrodynamics, sets the
field strength at which an electric field becomes nonlinear due to the
spontaneous production of electron-positron pairs [28, 29]. Quantitatively,
the Schwinger limit occurs at an electric field strength of
$\mathcal{E}_{C}=m_{e}^{2}/e=0.86\text{ MeV}^{2}$. This corresponds to a field
energy density of $u_{C}=\mathcal{E}_{C}^{2}=0.74\text{ MeV}^{4}$.

This field energy density can be related to $\epsilon$ as follows. We assume
the radiation from the BH is spread over a volume one wavelength ($\lambda$)
in thickness outside of the BH Schwarzschild radius. Then, the energy density
of the BH radiation is

$\displaystyle u$ $\displaystyle=\frac{\epsilon
M}{\frac{4}{3}\pi((r_{s}+\lambda)^{3}-r_{s}^{3})}$ (2)
$\displaystyle=3.1\times
10^{9}\epsilon\left(\frac{M}{M_{\odot}}\right)^{-2}\text{ MeV}^{4}.$

Setting $u=u_{C}$ and solving for the critical value of $\epsilon$ gives

$\epsilon_{C}=2.4\times 10^{-10}\left(\frac{M}{M_{\odot}}\right)^{2}.$ (3)

Therefore, $\epsilon_{C}\sim 10^{-10}-10^{-6}$ for BH masses ranging from 1 to
50 $M_{\odot}$. For $\epsilon>\epsilon_{C}$, pair-production dominates and
results in a GRB as discussed in Section III. For $\epsilon<\epsilon_{C}$, the
EM field accelerates ambient charged particles, creating cosmic ray electrons
and protons as discussed in Section IV.

##  III Gamma-Ray Emission Above the Schwinger Limit

Above the Schwinger limit, the energy density of the field is large enough to
result in electron-positron pair production [30]. Qualitatively, the electron-
positron-photon gas thermalizes due to Thompson scattering and expands
relativistically as an ideal fluid, creating an object known as a fireball in
GRB literature [31, 32]. We denote the lab frame of an Earth observer as S and
the comoving frame of the fluid as S’.

When the EM radiation is first emitted by the BH, the lab and fluid frames
coincide. The initial temperature of the fireball in both S and S’ is given by

$T_{0}=\left(\frac{E}{V_{0}g_{0}a}\right)^{1/4}$ (4)

where $a=\pi^{2}/15$, $E$ is the energy dumped into a region of volume
$V_{0}$, and $g_{0}=2.75=11/4$ is half of the effective degrees of freedom for
a plasma consisting of photons, electrons, and positrons in thermal
equilibrium [31]. Again assuming that the radiation from the BH is spread over
a volume one wavelength in thickness outside of the BH Schwarzschild radius,
the initial temperature is

$\displaystyle T_{0}$ $\displaystyle=\left(\frac{\epsilon
M}{\frac{4}{3}\pi((r_{s}+\lambda)^{3}-r_{s}^{3})\frac{11}{4}a}\right)^{1/4}$
(5)
$\displaystyle=200\epsilon^{1/4}\left(\frac{M}{M_{\odot}}\right)^{-1/2}\text{
MeV}.$

Assuming that any remnants of stellar ejecta and envelopes have long
dispersed, we neglect any external baryon contributions to the dynamics of the
fireball. Thus, the fireball is a relativistic radiation-dominated fluid,
which rapidly accelerates to $\gamma\gg 1$ under its own super-Eddington
radiation pressure. Because the fireball is created outside the Schwarzschild
radius of the BH, the system originates in a region of small curvature.
Consequently, general relativistic and gravitational redshift effects can be
neglected. Employing the usual relativistic conservation equations of baryon
number and energy-momentum from Ref. [33] in the limit where $\gamma\gg 1$,
yields the following scaling relations for each fluid shell [34, 35]

$\displaystyle\gamma(r)$ $\displaystyle\sim\left(\frac{r}{R_{0}}\right)$ (6)
$\displaystyle T^{\prime}(r)$ $\displaystyle\sim
T_{0}\left(\frac{r}{R_{0}}\right)^{-1}\sim\frac{T_{0}}{\gamma(r)}$ (7)

where $r$ is the distance from the origin in the lab frame and $R_{0}$ is the
initial width of the fireball. These relations apply so long as the fireball
is ultra-relativistic, radiation-dominated, and opaque due to Thompson
scattering. So, as the fireball expands from the origin, the bulk Lorentz
factor continues to increase due to the acceleration from the radiation
pressure of the fluid. To first order in $\gamma$, the width of the fireball
in the lab frame is constant $R(r)=R_{0}$ [36]. This requires the width in the
comoving frame to increase as $R^{\prime}(r)=\gamma(r)R_{0}$, illustrating why
the fireball cools in its co-moving frame. In the lab frame, the temperature
is blue-shifted by

$T(r)=\gamma(r)T^{\prime}(r)=T_{0}$ (8)

since the fluid is moving relativistically towards the observer. Thus, so long
as the scaling relations apply, a lab observer sees each shell of the fireball
at the same constant temperature, $T_{0}$.

Within each fluid shell, the number density of electron-positron pairs in the
comoving frame decreases as the fireball cools. Eventually, the process of
pair creation and annihilation freezes out when the time for a positron to
annihilate with an electron is of the same order as the dynamical time. This
occurs at a comoving temperature of $T^{\prime}\sim 20\text{ keV}$ [37, 35].
At this temperature, the proportion of the initial energy from the BH
contained in the remaining electron-positron pairs is negligible [32]. So,
nearly all of the initial $\epsilon M$ is contained in photons that had been
trapped in the fluid by Thompson scattering. When pair production freezes out,
the Thompson opacity decreases dramatically and these photons escape. Since
the comoving temperature depends only on $r$, each shell experiences this
freeze out as it moves through the same radius in the lab frame. Thus, the
characteristic time delay between when photons free-stream from the inner and
outermost edges of the fireball is given by [31]

$\delta t\sim\frac{R_{0}}{c}=2.0\times
10^{-5}\left(\frac{M}{M_{\odot}}\right)\text{ s}$ (9)

since the initial radius of the fireball is set by the wavelength of the BH
radiation, $\lambda$. These timescales are short enough to be consistent with
a short GRB, $<O(1$ second).

When the photons are able to free stream from the fireball, an observer on
Earth sees a nearly thermal black body spectrum at temperature $T_{0}$
radiating from the BH [37]. The black body spectrum in photon number as a
function of photon frequency and the fireball temperature is given by

$B_{f}(T_{0})=\frac{2f^{2}}{e^{2\pi f/T_{0}}-1}\text{ MeV}^{2}\text{
s}^{-1}\text{ Hz}^{-1}\text{ sr}^{-1}.$ (10)

The peak photon energy for this spectrum is

$E_{peak}=1.6T_{0}=320\epsilon^{1/4}\left(\frac{M}{M_{\odot}}\right)^{-1/2}\text{
MeV}.$ (11)

These peak energies are plotted in Figure 1 as a function of BH mass and
$\epsilon$. These are the most likely photons to be emitted from the fireball,
and all peak energies are gamma rays (energy $>1$ MeV).

Figure 1: The peak photon energy, $E_{peak}$, as a function of BH mass and
$\epsilon$ for $\epsilon>\epsilon_{C}$. These peak energies are $>1$ MeV and
so are gamma rays.

Since nearly all the energy from the fireball is converted into photons at the
peak energy, the number of photons emitted is

$N_{\gamma}=\frac{\epsilon M}{E_{peak}}=3.4\times
10^{57}\epsilon^{3/4}\left(\frac{M}{M_{\odot}}\right)^{3/2}.$ (12)

Over the full range of photon energies shown in Figure 1, the most sensitive
telescope is the Fermi Gamma-ray Space Telescope. The two instruments onboard
Fermi are the Large Area Telescope (LAT) and the Gamma-ray Burst Monitor
(GBM). The LAT observes photon energies in the range $20$ MeV$-300$ GeV with a
sensitivity of $10^{-4}$ erg cm-2 [38]. The GBM observes photon energies in
the range $8$ keV$-40$ MeV with a sensitivity of $0.5$ ph cm-2 s-1 [39]. For
the gamma-ray energies emitted by 1-50 $M_{\odot}$ BHs over a timescale of
$\lesssim 1$ second, Fermi’s sensitivity limits requires a flux of $F\gtrsim
1$ ph cm-2. This minimum flux can be used to calculate the maximum distance,
$d$, to which EM emission by a BH above the Schwinger limit is observable.

$d=\sqrt{\frac{N_{\gamma}}{4\pi
F}}=5360\epsilon^{3/8}\left(\frac{M}{M_{\odot}}\right)^{3/4}\text{ Mpc}$ (13)

These distances are plotted in Figure 2. For a given distance, we can solve
for the minimum $\epsilon$ for which Fermi could observe such an event
extragalactically.

$\epsilon_{min}=2\times 10^{-5}\left(\frac{d}{100\text{
Mpc}}\right)^{8/3}\left(\frac{M}{M_{\odot}}\right)^{-2}$ (14)

Figure 2: The distance to which the gamma-ray photons emitted by 1-50
$M_{\odot}$ BHs are observable with Fermi, $d$, for various
$\epsilon>\epsilon_{C}$.

The Gravitational-wave (GW) Transient Catalog lists events detected by LIGO
[15], Virgo [16], and the Kamioka Gravitational Wave Detector (KAGRA, [40]).
Most BBH merger candidates listed in the catalog are posted on NASA’s General
Coordinates Network, which contains both concurrent and follow-up EM
observations of GW triggers by the Fermi GBM. The Fermi GBM has an 8 steradian
field-of-view [39], which is large enough to cover the full LIGO/Virgo 90%
confidence GW localization region if the telescope is well-aligned at the time
of the trigger. We cross-reference all BBH mergers in the GW Transient Catalog
with Fermi observations from the General Coordinates Network. For BH masses
ranging from 10 to 50 $M_{\odot}$ in intervals of 10 $M_{\odot}$, we identify
the nearest BBH merger for which Fermi observed at least 90% of the
localization region and recorded no GRB event. The GW events constraining each
mass interval are listed in Table 1. These non-detections are used to
constrain $\epsilon$ for each BH mass. To do this, we assume the furthest
distance (luminosity distance + error bar) measured by LIGO/Virgo and
calculate the minimum value of $\epsilon$ needed such that the event would be
observable with the Fermi GBM for the observed BH mass. All constrained values
of $\epsilon$ are above the Schwinger limit for the given BH mass and all
result in $\sim$ MeV photons which are in the energy range observable by the
GBM. These constraints are listed in Table 2. Assuming all BHs are “hairy” BHs
capable of producing this signal, the current upper bounds on $\epsilon$ are
$\epsilon<10^{-5}$ for 10, 30, 40 $M_{\odot}$ BHs and $\epsilon<10^{-4}$ for
20, 50 $M_{\odot}$ BHs since no high energy EM signal was observed from these
BBH mergers. These constraints will improve as more GW events with concurrent
Fermi observations are detected.

Table 1: The table lists the physical parameters for the GW events constraining $\epsilon$ for each BH mass interval. GW191216_213338-v1 was reported in Ref. [46], and all other events were reported in Ref. [47]. The table provides the event ID from the GW Transient Catalog, the measured constituent BH masses ($m_{1},$ $m_{2}$), the measured luminosity distance ($d_{L}$), the General Coordinates Network circular, and the percent of the localization region observed by the Fermi GBM (% LR). ID | $m_{1}/M_{\odot}$ | $m_{2}/M_{\odot}$ |  $d_{L}$/Mpc | GCN Reference | % LR
---|---|---|---|---|---
GW191216_213338 | $12.1_{-2.3}^{+4.6}$ | $7.7_{-1.9}^{+1.6}$ | $340_{-130}^{+120}$ | [41] | 99.8
GW190915_235702 | $32.6_{-4.9}^{+8.8}$ | $24.5_{-5.8}^{+4.9}$ | $1750_{-650}^{+710}$ | [42] | 98.0
GW190412_053044 | $27.7_{-6.0}^{+6.0}$ | $9.0_{-1.4}^{+2.0}$ | $720_{-220}^{+240}$ | [43] | 99.8
GW190521_074359 | $43.4_{-5.5}^{+5.8}$ | $33.4_{-6.8}^{+5.2}$ | $1080_{-530}^{+580}$ | [44] | 100.0
GW190701_203306 | $54.1_{-8.0}^{+12.6}$ | $40.5_{-12.1}^{8.7}$ | $2090_{-740}^{+770}$ | [45] | 100.0
Table 2:  This table lists the constraints on $\epsilon$, assuming all BHs are “hairy” BHs, for each mass range. The columns give the event ID from the GW Transient Catalog, the BH mass range ($O(m)$), the maximum luminosity distance ($d_{L,Max}$), the minimum value of $\epsilon$ such that the event would be observable ($\epsilon_{Min}$), and the peak photon energy ($E_{peak}$) for that $\epsilon_{min}$ and BH mass. ID | $O(m)/M_{\odot}$ |  $d_{L,Max}$/Mpc | $\epsilon_{Min}$ |  $E_{peak}$/MeV
---|---|---|---|---
GW191216_213338 | 10 | 460 | $1.4\times 10^{-5}$ | 6.3
GW190915_235702 | 20 | 2460 | $3.1\times 10^{-4}$ | 9.6
GW190412_053044 | 30 | 960 | $1.1\times 10^{-5}$ | 3.4
GW190521_074359 | 40 | 1660 | $2.7\times 10^{-5}$ | 3.7
GW190701_203306 | 50 | 2860 | $7.5\times 10^{-5}$ | 4.3

There has been one observation of a GRB and BBH merger occurring concurrently.
An offline search of Fermi GBM data following the detection of GW150914 [48]
revealed a 1 second short GRB occurring 0.4 seconds after the LIGO trigger
with a localization consistent with that of the GW signal [49]. The GRB
transient was near the detection threshold for the GBM (2.9$\sigma$
detection), was not detected by any other instrument, and the localization was
poorly constrained [49]. Thus, the GRB cannot be confidently associated with
GW150914. Assuming that the GRB did indeed originate from the BBH merger, we
can identify the range of $\epsilon$ that is consistent with the measured
properties of the merger and GRB. From the GW Transient Catalog, the
constituent BH masses range from $\sim$30-40 $M_{\odot}$ and the range of
luminosity distances is 270-590 Mpc. For these masses and distances, the
minimum value of $\epsilon$ such that the event is observable with the Fermi
GBM ranges from $\epsilon_{min}\sim 10^{-7}-10^{-6}$, resulting in peak photon
energies of $\sim 1-3$ MeV. This range of peak photon energies is consistent
with the properties of the observed GRB, which peaked near an MeV [49]. Since
this event was barely above the GBM’s detection threshold, we anticipate
$\epsilon\sim\epsilon_{min}$. Therefore, the observed GRB is consistent with a
GRB produced via rapid EM emission directly from a 30-40 $M_{\odot}$ “hairy”
BH for $\epsilon\sim 10^{-7}-10^{-6}$.

##  IV Galactic Cosmic Ray Signal Below the Schwinger Limit

Below the Schwinger limit the electric field emitted by the BH propagates
outwards, accelerating ambient charged particles. Mutual attraction between
the protons and electrons inhibits charge separation, requiring that the
particles have the same Lorentz factor on average. Since the protons are a
factor of $\sim 10^{3}$ heavier than the electrons and the particle energy
scales with mass, most of the BH field energy is absorbed by protons. Thus,
the dynamics of the system are set by the protons.

A rapid burst of isotropic EM radiation can be generically described by a
coherent single-wavelength pulse. The exact duration and coherence of emission
is model-dependent, but so long as the emission occurs on the short timescales
characteristic of BBH mergers, the signal is sufficiently similar to the
single-wavelength pulse approximation. First, we consider the acceleration of
a single charged particle due to this strong EM pulse. In the (+ - - -)
metric, the relativistic Lorentz force law is

$\frac{du^{\alpha}}{d\tau}=\frac{e}{m}F^{\alpha\beta}u_{\beta}$ (15)

where $e$ is the elementary charge, $m$ is the proton or electron mass,
$u^{\alpha}$ is the particle’s 4-velocity, and $\tau$ is the proper time in
the instantaneous rest frame of the particle. For a sinusoidal EM field, the
characteristic timescale of field variation is
$\mathcal{E}/\frac{d\mathcal{E}}{dt}\sim 1/f$ where $\mathcal{E}$ is the
electric field magnitude and $f$ is the frequency of the BH radiation. Since
$f\sim 10^{-17}$ MeV is very small, this characteristic timescale is very
long. So, the field magnitude is nearly constant in time and the system can be
approximated by a constant crossed field. For the ranges of $\epsilon$
considered here, the absorption length for the EM field is long enough that
the average charged particle experiences a sufficiently weak field to neglect
any special relativistic effects that transform the field but also a
sufficiently strong field that the protons are quickly accelerated to $v\sim
1$ and radiate negligibly. We choose a coordinate system such that the
electric field is directed along the x-axis and the magnetic field along the
y-axis, then the Faraday tensor is

$F^{\alpha\beta}=\begin{pmatrix}0&-\mathcal{E}&0&0\\\ \mathcal{E}&0&0&0\\\
0&0&0&\mathcal{E}\\\ 0&0&-\mathcal{E}&0\end{pmatrix}$ (16)

where $\mathcal{E}$ is the magnitude of the electric field in MeV2. Since
$\frac{e}{m}F^{\alpha\beta}$ is constant, the Lorentz force equation has a
matrix exponential solution

$u^{\alpha}(\tau)=\exp{\left(\frac{e}{m}\tau F^{\alpha}_{\text{
}\beta}\right)}u^{\beta}(0).$ (17)

The particles accelerate from the thermal speed of the plasma ($v\sim 0$),
fixing $u^{\beta}(0)$. We compute $u^{\alpha}(\tau)$ for the given Faraday
tensor and use this to extract the Lorentz factor, $\gamma$, and the
components of the 3-velocity in the lab frame, v.

$\displaystyle\gamma(\tau)$
$\displaystyle=1+\frac{e^{2}\mathcal{E}^{2}\tau^{2}}{2m^{2}}$ (18)
$\displaystyle v_{x}(\tau)$
$\displaystyle=\frac{2e\mathcal{E}m\tau}{2m^{2}+e^{2}\mathcal{E}^{2}\tau^{2}}$
(19) $\displaystyle v_{y}(\tau)$ $\displaystyle=0$ (20) $\displaystyle
v_{z}(\tau)$
$\displaystyle=1-\frac{2m^{2}}{2m^{2}+e^{2}\mathcal{E}^{2}\tau^{2}}$ (21)

The final kinetic energy of the particle depends on the time at which the
particle exits the field in its instantaneous rest frame, $\tau_{f}$. By
definition, $\gamma(\tau)=\frac{dt}{d\tau}$. Thus,

$\int_{0}^{\tau_{f}}\gamma(\tau)d\tau=\tau_{f}+\frac{e^{2}\mathcal{E}^{2}\tau_{f}^{3}}{6m^{2}}=\int_{0}^{t_{f}}dt=\frac{1}{f}$
(22)

since in the lab frame the particle is in the field for $t_{f}=\frac{1}{f}$.
The magnitude of the EM field is sufficiently large, and we can neglect the
term that is linear in $\tau_{f}$ to find

$\displaystyle\tau_{f}$
$\displaystyle=\left(\frac{6m^{2}}{e^{2}\mathcal{E}^{2}f}\right)^{1/3}$ (23)
$\displaystyle=1.2\times
10^{8}\frac{1}{\mathcal{E}^{2/3}}\left(\frac{M}{M_{\odot}}\right)^{1/3}\text{
MeV}^{-1},$

where here, and in all following numerical expressions, we set $m=m_{p}$.

Now we characterize the absorption of energy from the BH EM field by ambient
protons. We partition the volume around the BH into spherical shells one
wavelength in thickness such that the distance from the BH is parameterized by
the dimensionless number $j$, the number of wavelengths from the Schwarzschild
radius of BH. Let $E$ be the total EM energy incident onto a shell $j$
wavelengths from the BH. Assuming $j$ is large, the initial energy density of
the shell is

$\displaystyle u$
$\displaystyle=\frac{E}{\frac{4}{3}\pi[(r_{s}+j\lambda)^{3}-(r_{s}+(j-1)\lambda)^{3}]}$
(24) $\displaystyle\sim\frac{E}{4\pi\lambda^{3}j^{2}}$
$\displaystyle=3.0\times
10^{-51}\frac{E}{j^{2}}\left(\frac{M}{M_{\odot}}\right)^{-3}\text{ MeV}^{4}.$

Then the magnitude of the electric field is

$\displaystyle\mathcal{E}$ $\displaystyle=\sqrt{u}$ (25)
$\displaystyle=5.4\times
10^{-26}\left(\frac{E}{j^{2}}\right)^{1/2}\left(\frac{M}{M_{\odot}}\right)^{-3/2}\text{
MeV}^{2}.$

The kinetic energy of each proton after interacting with the EM wave is

$\displaystyle K$ $\displaystyle\sim\gamma(\tau_{f})m$ (26)
$\displaystyle=1.0\times
10^{-5}\left(\frac{E}{j^{2}}\right)^{1/3}\left(\frac{M}{M_{\odot}}\right)^{-1/3}\text{
MeV}.$

This gives the kinetic energy of one proton in the $j$th shell. To derive the
total kinetic energy lost in the $j$th shell, we assume a homogeneously
distributed number density of charged particles with $n=n_{p}=n_{e}\sim
1\text{ cm}^{-3}\sim 10^{-32}\text{ MeV}^{3}$, consistent with the Milky Way
interstellar medium [50, 51]. The total number of protons in the $j$th shell
is

$N\sim 4\pi n\lambda^{3}j^{2}=2.6\times
10^{18}j^{2}\left(\frac{M}{M_{\odot}}\right)^{3}.$ (27)

So the total kinetic energy absorbed in the $j$th shell is

$K_{tot}=NK=2.6\times
10^{13}j^{4/3}E^{1/3}\left(\frac{M}{M_{\odot}}\right)^{8/3}\text{ MeV}.$ (28)

This gives a differential equation for the energy absorbed by the field

$\frac{dE}{dj}=-K_{tot}$ (29)

subject to the initial condition $E(0)=\epsilon M$. This equation can be
integrated to yield

$E(j)=\left(1.1\times
10^{40}\epsilon^{2/3}\left(\frac{M}{M_{\odot}}\right)^{2/3}-7.5\times
10^{12}\left(\frac{M}{M_{\odot}}\right)^{8/3}j^{7/3}\right)^{3/2}\text{ MeV}.$
(30)

Solving for the value of $j$ at which $E(j)=0$ gives the dimensionless
absorption length, $j_{abs}$,

$j_{abs}=4.4\times
10^{11}\epsilon^{2/7}\left(\frac{M}{M_{\odot}}\right)^{-6/7}.$ (31)

The absorption length, $j_{abs}\lambda$, gives the number of protons
accelerated by the field and the average kinetic energy and Lorentz factor of
each proton.

$\displaystyle N_{abs}$ $\displaystyle=\frac{4}{3}\pi
n(j_{abs}\lambda)^{3}=7.1\times
10^{52}\epsilon^{6/7}\left(\frac{M}{M_{\odot}}\right)^{3/7}$ (32)
$\displaystyle K_{avg}$ $\displaystyle=\frac{\epsilon M}{N_{abs}}=1.6\times
10^{7}\epsilon^{1/7}\left(\frac{M}{M_{\odot}}\right)^{4/7}\text{ MeV}$ (33)
$\displaystyle\gamma_{avg}$ $\displaystyle\sim\frac{KE_{avg}}{m}=1.7\times
10^{4}\epsilon^{1/7}\left(\frac{M}{M_{\odot}}\right)^{4/7}$ (34)

The number of electrons accelerated by the field is also $N_{abs}$ and the
average Lorentz factor of the electrons must be the same as that for the
protons.

We perform this calculation for values of $\epsilon$ ranging from
$\epsilon=10^{-20}$ to the Schwinger limit for 1-50 $M_{\odot}$ BHs. The
resulting average proton kinetic energies as a function of $\epsilon$ are
plotted in Figure 3. Overall, the average kinetic energy per proton ranges
from 20 GeV for a 1 $M_{\odot}$ BH with $\epsilon=10^{-20}$ to 20 TeV for a 50
$M_{\odot}$ BH at the Schwinger limit ($\epsilon\sim 10^{-6}$). The average
proton kinetic energies can be used to calculate the corresponding average
electron kinetic energies. Overall, the average kinetic energy per electron
ranges from 0.01 GeV for a 1 $M_{\odot}$ BH with $\epsilon=10^{-20}$ to 10 GeV
for a 50 $M_{\odot}$ BH at the Schwinger limit. These relativistic protons and
electrons are cosmic rays.

Figure 3: The average proton kinetic energy, $K_{avg}$, for 1-50 $M_{\odot}$
BHs as a function of $\epsilon$ for $\epsilon<\epsilon_{C}$. These
relativistic protons are at cosmic ray energies.

Cosmic rays of these energies are difficult to identify with a point source on
the sky. Ambient magnetic fields confine these cosmic rays to their host
galaxies and cause them to quickly diffuse (on timescales $\lesssim 0.01$
years) into the galactic background of cosmic rays. Thus, these cosmic rays
are indistinguishable from cosmic rays produced by supernova remnants [38],
stellar winds or flares [52], and other processes. Although these particles
create secondary signals as they lose energy to bremsshtrahlung, ionization,
synchrotron radiation and inverse Compton scattering for electrons and to
inelastic collisions for protons [53, 54], the timescales for these processes
is much longer than the diffusion time ($\gtrsim 10^{6}$ years). Therefore,
these secondary signals are also indistinguishable from those produced by
other cosmic ray processes.

Because these cosmic rays are mixed with and are indistinguishable from cosmic
rays produced via other mechanisms, it is not possible to confidently identify
cosmic rays due to EM radiation from BHs either in external galaxies or in the
Milky Way. Therefore, it becomes difficult to place strong constraints on
$\epsilon$ below the Schwinger limit. Although this avenue cannot constrain
$\epsilon$, the effect is still potentially observable. Should the BBH merger
occur on a sufficiently strong magnetic field background ($B>10$
MeV${}^{2}=0.05$ T), the ultrarelativistic electrons would produce synchrotron
radiation in the X-ray band, motivating X-ray observations of BBH mergers.

##  V Discussion and Conclusions

The theoretical understanding of BHs is inconsistent, motivating searches for
generic signals of deviations from canonical BH models. In this paper, we have
constrained a broad class of “hairy” BH models capable of emitting a fraction
of their mass as EM radiation. Since this radiation is sourced directly from
the BH, it must tunnel out of the BH’s gravitational well in the same manner
as Hawking radiation. Thus, the characteristic frequency of the radiation
depends only on the mass of the BH, resulting in a signal that is generic and
model-independent. We derive the critical value of $\epsilon$, the fraction of
the BH mass released as radiation, above which the field strength triggers a
GRB and below which ambient particles are accelerated to cosmic ray energies.
Because no extragalactically-observable EM signal is expected from a stellar-
mass BBH merger, we find that concurrent observations of BBH mergers with GW
detectors and EM radiation instruments offer the best data to detect such a
signal.

In the GRB regime, the BH mass and $\epsilon$ fix the initial volume and
temperature of the electron-positron fireball. The fireball expands
relativistically, maintaining constant temperature in the frame of an Earth
observer and cooling in its comoving frame. Once the fireball is sufficiently
cool in its frame, pair-production freezes out and the photons free stream. In
the frame of an Earth observer, these photons have energies described by a
black body spectrum at the initial temperature of the fireball. Thus, the
energy deposited by the BH is re-emitted as gamma-rays over a short timescale.
By cross-referencing GW events with concurrent Fermi GBM observations of the
localization region, we place upper bounds on $\epsilon$. These bounds are
$\epsilon<10^{-5}-10^{-4}$ for $10-50$ $M_{\odot}$ BHs depending on the BH
mass since no high energy EM signal was observed from these BBH mergers. These
constraints will improve as more GW events with concurrent Fermi observations
are detected. We also discuss the weak detection of a GRB following GW150914,
and find that this event is consistent with a GRB produced via rapid EM
emission directly from a “hairy” BH for $\epsilon\sim 10^{-7}-10^{-6}$.

Below the Schwinger limit, the EM radiation can be described by a constant-
crossed field. The dynamics of the system are fixed by the ambient protons,
which are rapidly accelerated to $v\sim 1$ by the field, absorbing energy as
the radiation propagates away from the BH. We solve the differential equation
describing the energy lost by the field to calculate the absorption length as
a function of BH mass and $\epsilon$. This absorption length fixes the average
energy of the ambient protons and electrons that interacted with the BH
radiation. For 1-50 $M_{\odot}$ BHs and $\epsilon$ ranging from $10^{-20}$ to
the Schwinger limit, the average kinetic energy per proton ranges from 20
GeV-20 TeV and the energy per electron ranges from 0.01-10 GeV. At these
energies, cosmic rays have a short diffusion length due to the galactic
magnetic field and are mixed in with other astrophysical cosmic rays.
Additionally, the secondary signals from cosmic rays of these energies are
produced on too long of a timescale to be attributed to EM radiation directly
from a BH. Overall, constraining $\epsilon$ in this less energetic regime is
difficult. Future work could investigate BBH mergers in strong background
magnetic fields. In this case, the ultrarelativistic electrons emit X-rays via
synchrotron radiation that may be observable.

Although this work benefits from employing a model-independent approach to
generically characterize radiation emitted directly from “hairy” BHs, model-
dependent effects will need to be included to augment this general
parameterization. Some “hairy” BH models, such as the firewall BH [13], are
capable of producing EM radiation. But currently, no complete model able to
quantitatively characterize this effect exists. We hope this paper will
motivate others in the field to work through the details of such models. The
firewall BH metric of Ref. [13] is particularly well-suited for this. Given
this metric and a parameterized charge distribution adhered to the firewall,
one could use numerical relativity to simulate the emission of gravitational
and EM radiation during a BBH merger involving a firewall BH. This approach
offers independent constraints on $\epsilon$ as a function of the charge
distribution parameter and may be able to constrain $\epsilon$ in cases below
the Schwinger limit.

Strengthening constraints on $\epsilon$ will be an important pursuit for
future work given the relatively loose bounds found in this paper. In general,
“hairy” BHs are well-motivated by the BH information paradox. However, such
BHs must appear nearly canonical to a distant observer due to observational
evidence in favor of the “no hair” theorem. Therefore, novel approaches to
constraining “hairy” models will continue to be vital in the search for new
fundamental physics.

###### Acknowledgements.

We thank Surjeet Rajendran and David Kaplan for useful discussions. We also
thank Nadia Zakamska and Erwin Tanin for their edits to this manuscript. This
work is supported by the National Science Foundation Graduate Research
Fellowship Program under Grant No. DGE2139757. Any opinions, findings, and
conclusions or recommendations expressed in this material are those of the
author and do not necessarily reflect the views of the National Science
Foundation.

## References

  * [1] B.P. Abbott, R. Abbott, T.D. Abbott, S. Abraham, F. Acernese, K. Ackley et al., GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Physical Review X 9 (2019) 031040.
  * [2] A. Vayner, N. Zakamska, S.A. Wright, L. Armus, N. Murray and G. Walth, Multiphase Outflows in High-redshift Quasar Host Galaxies, The Astrophysical Journal 923 (2021) 59.
  * [3] S. Komossa, Tidal disruption of stars by supermassive black holes: Status of observations, Journal of High Energy Astrophysics 7 (2015) 148.
  * [4] A.E. Broderick, T. Johannsen, A. Loeb and D. Psaltis, Testing the No-hair Theorem with Event Horizon Telescope Observations of Sagittarius A*, The Astrophysical Journal 784 (2014) 7 [1311.5564].
  * [5] D. Psaltis, N. Wex and M. Kramer, A Quantitative Test of the No-hair Theorem with Sgr A* Using Stars, Pulsars, and the Event Horizon Telescope, The Astrophysical Journal 818 (2016) 121 [1510.00394].
  * [6] D. Wang, Shaving the Hair of Black Hole with Sagittarius A∗ from Event Horizon Telescope, arXiv e-prints (2022) arXiv:2205.08026 [2205.08026].
  * [7] M. Isi, M. Giesler, W.M. Farr, M.A. Scheel and S.A. Teukolsky, Testing the No-Hair Theorem with GW150914, Phys. Rev. Lett. 123 (2019) 111102 [1905.00869].
  * [8] K. Wang, Retesting the no-hair theorem with GW150914, European Physical Journal C 82 (2022) 125 [2111.00953].
  * [9] L. Sadeghian and C.M. Will, Testing the black hole no-hair theorem at the galactic center: perturbing effects of stars in the surrounding cluster, Classical and Quantum Gravity 28 (2011) 225029 [1106.5056].
  * [10] H. Qi, R. O’Shaughnessy and P. Brady, Testing the black hole no-hair theorem with Galactic Center stellar orbits, Phys. Rev. D 103 (2021) 084006 [2011.02267].
  * [11] S.W. Hawking, Breakdown of predictability in gravitational collapse, Phys. Rev. D 14 (1976) 2460.
  * [12] S. Raju, Lessons from the information paradox, Physics Reports 943 (2022) 1.
  * [13] D.E. Kaplan and S. Rajendran, Firewalls in general relativity, Phys. Rev. D 99 (2019) 044033 [1812.00536].
  * [14] S.W. Hawking, Particle creation by black holes, Communications in Mathematical Physics 43 (1975) 199.
  * [15] J. Aasi, J. Abadie, B.P. Abbott, R. Abbott, T. Abbott, M.R. Abernathy et al., Characterization of the LIGO detectors during their sixth science run, Classical and Quantum Gravity 32 (2015) 115012 [1410.7764].
  * [16] F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca et al., Advanced Virgo: a second-generation interferometric gravitational wave detector, Classical and Quantum Gravity 32 (2015) 024001 [1408.3978].
  * [17] I. Mandel and S. Mink, Merging binary black holes formed through chemically homogeneous evolution in short-period stellar binaries, Monthly Notices of the Royal Astronomical Society 458 (2015) .
  * [18] A. Loeb, Electromagnetic Counterparts to Black Hole Mergers Detected by LIGO, The Astrophysical Journal 819 (2016) L21 [1602.04735].
  * [19] R. Perna, D. Lazzati and B. Giacomazzo, Short Gamma-Ray Bursts from the Merger of Two Black Holes, The Astrophysical Journal 821 (2016) L18 [1602.05140].
  * [20] S.E. Woosley, The Progenitor of GW150914, The Astrophysical Journal 824 (2016) L10 [1603.00511].
  * [21] L. Dai, J.C. McKinney and M.C. Miller, Energetic constraints on electromagnetic signals from double black hole mergers, Monthly Notices of the Royal Astronomical Society 470 (2017) L92 [1611.00764].
  * [22] J.M. Fedrow, C.D. Ott, U. Sperhake, J. Blackman, R. Haas, C. Reisswig and A. De Felice, Gravitational Waves from Binary Black Hole Mergers inside Stars, Phys. Rev. Lett. 119 (2017) 171103 [1704.07383].
  * [23] S.S. Kimura, S.Z. Takahashi and K. Toma, Evolution of an accretion disc in binary black hole systems, Monthly Notices of the Royal Astronomical Society 465 (2017) 4406 [1607.01964].
  * [24] B. Zhang, Mergers of Charged Black Holes: Gravitational-wave Events, Short Gamma-Ray Bursts, and Fast Radio Bursts, The Astrophysical Journal 827 (2016) L31 [1602.04542].
  * [25] M. Lyutikov, Fermi GBM signal contemporaneous with GW150914 - an unlikely association, arXiv e-prints (2016) arXiv:1602.07352 [1602.07352].
  * [26] M.K. Parikh and F. Wilczek, Hawking Radiation As Tunneling, Phys. Rev. Lett. 85 (2000) 5042 [hep-th/9907001].
  * [27] R.J. Reynolds, The Low Density Ionized Component of the Interstellar Medium and Free-Free Absorption at High Galactic Latitudes, in _Low Frequency Astrophysics from Space_ , N.E. Kassim and K.W. Weiler, eds., vol. 362, p. 121 (1990), DOI.
  * [28] W. Heisenberg and H. Euler, Consequences of Dirac’s theory of positrons, Z. Phys. 98 (1936) 714 [physics/0605038].
  * [29] J. Schwinger, On gauge invariance and vacuum polarization, Phys. Rev. 82 (1951) 664.
  * [30] R. Lieu, Y. Takahashi and T.W.B. Kibble, Gamma Ray Burst as Vacuum Discharge of Super-Schwinger Electric Fields, arXiv e-prints (1998) astro [astro-ph/9803072].
  * [31] J. Goodman, Are gamma-ray bursts optically thick?, The Astrophysical Journal 308 (1986) L47.
  * [32] B. Paczynski, Gamma-ray bursters at cosmological distances, The Astrophysical Journal 308 (1986) L43.
  * [33] S. Weinberg, _Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity_ (1972).
  * [34] T. Piran, A. Shemi and R. Narayan, Hydrodynamics of Relativistic Fireballs, Monthly Notices of the Royal Astronomical Society 263 (1993) 861.
  * [35] P. Kumar and B. Zhang, The physics of gamma-ray bursts & relativistic jets, Physics Reports 561 (2015) 1 [1410.0679].
  * [36] P. Mészáros, _Gamma-ray bursts_ , Reports on Progress in Physics 69 (2006) 2259 [astro-ph/0605208].
  * [37] T. Piran, Gamma-ray bursts and the fireball model, Physics Reports 314 (1999) 575 [astro-ph/9810256].
  * [38] C.D. Dermer, Sources of GeV Photons and the Fermi Results, Saas-Fee Advanced Course 40 (2013) 225 [1202.2814].
  * [39] The Fermi GBM Collaboration, Overview of the Fermi GBM, 2020.
  * [40] Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa, T. Sekiguchi, D. Tatsumi and H. Yamamoto, Interferometer design of the KAGRA gravitational wave detector, Phys. Rev. D 88 (2013) 043007 [1306.6747].
  * [41] LIGO and Virgo Collaborations, GCN 26454, .
  * [42] LIGO and Virgo Collaborations, GCN 25752, .
  * [43] LIGO and Virgo Collaborations, GCN 24098, .
  * [44] LIGO and Virgo Collaborations, GCN 24629, .
  * [45] LIGO and Virgo Collaborations, GCN 24948, .
  * [46] R. Abbott, The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, T.D. Abbott, F. Acernese et al., GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, arXiv e-prints (2021) arXiv:2111.03606 [2111.03606].
  * [47] R. Abbott, Abbott, The LIGO Scientific Collaboration, T.D. the Virgo Collaboration, F. Acernese, K. Ackley et al., GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, arXiv e-prints (2021) arXiv:2108.01045 [2108.01045].
  * [48] B.P. Abbott, R. Abbott, T.D. Abbott, M.R. Abernathy, F. Acernese, K. Ackley et al., Localization and Broadband Follow-up of the Gravitational-wave Transient GW150914, The Astrophysical Journal 826 (2016) L13 [1602.08492].
  * [49] V. Connaughton, E. Burns, A. Goldstein, L. Blackburn, M.S. Briggs, B.B. Zhang et al., Fermi GBM Observations of LIGO Gravitational-wave Event GW150914, The Astrophysical Journal 826 (2016) L6 [1602.03920].
  * [50] J. Canto, On the density and energy of supernova remnants., Astronomy and Astrophysics 61 (1977) 641.
  * [51] R.J. Reynolds, The warm ionized medium, AIP Conference Proceedings 278 (1992) 156 [https://aip.scitation.org/doi/pdf/10.1063/1.44005].
  * [52] L.I. Miroshnichenko, Solar Cosmic Rays, vol. 260 (2001), 10.1007/978-94-015-9646-6.
  * [53] M.S. Longair, _High energy astrophysics. Vol.1: Particles, photons and their detection_ (1992).
  * [54] M.S. Longair, _High energy astrophysics. Vol.2: Stars, the galaxy and the interstellar medium_ , vol. 2 (1994).


Email<EMAIL_ADDRESS>
#  Electronic transport in copper-graphene composites

Kashi N. Subedi  Department of Physics and Astronomy, Ohio University, Athens,
OH 45701, USA  Kishor Nepal  Department of Physics and Astronomy, Ohio
University, Athens, OH 45701, USA  Chinonso Ugwumadu  Department of Physics
and Astronomy, Ohio University, Athens, OH 45701, USA  Keerti Kappagantula
Pacific Northwest National Laboratory, Richland, WA 99352, USA  D. A. Drabold
[  Department of Physics and Astronomy, Ohio University, Athens, OH 45701, USA

###### Abstract

We investigate the electronic transport properties of copper-graphene
composites using a density-functional framework. Conduction in composites by
varying the interface distance of a copper/graphene/copper (Cu/G/Cu) interface
models was studied. The electronic density of states reveals increasing
contributions from both copper and carbon atoms near the Fermi level with
decreasing Cu-G interfacial distance. Electronic conductivity of the models
computed using the Kubo-Greenwood formula showed the conductivity increases
with decreasing Cu-G distance. We also find that the conductivity saturates
below a threshold Cu-G distance. By computing the space-projected conductivity
of the Cu/G/Cu models, we show that the graphene forms a bridge to the
electronic conduction at small copper-graphene distances, thereby enhancing
the conductivity.

Metal composites are key industrial materials that may be manufactured using
various novel techniques Heringhaus and Raabe (1996); Bouaziz, Kim, and Estrin
(2013); Bruschi  _et al._ (2021), leading to metastable materials exhibiting
desirable mechanical, thermal, electronic, and transport properties without
the constraints typically seen in alloys. Copper is used extensively for
energy transport and electromagnetic applications such as motors, generators,
cable, busbars and transformers. In the last two decades, we have witnessed
copper composites with improved transport properties Singh and Gautam (2018),
devised as a means to reduce energy loss in applications.

One material that has demonstrated improved electrical conductivity and
current density is a copper-graphene (Cu-G) composite. Various Cu-G composites
ranging from a simple interface laminated to bulk Cu with interconnected
graphene flakes have been reported in the literature, with some showing
enhanced conductivity Weiping  _et al._ (2015); Rho  _et al._ (2017); An  _et
al._ (2019); Yang  _et al._ (2021), and more commonly enhancement in
mechanical performance Hwang  _et al._ (2013); An  _et al._ (2019).
Conventional wisdom suggests that the introduction of additives to metal would
increase scattering that would lead to increased electrical resistivity Ellis
_et al._ (1993). This is the reason why addition of silver to copper leads to
reduced electrical conductivity of Cu, even though Ag is 8% higher in
conductivity than Cu Felicia, Rochiem, and Laia (2018). The additives, on the
other hand, also arrest dislocation movement in the metals due to several
mechanisms such as precipitate hardening, solid solution strengthening and
dispersion strengthening, resulting in enhanced mechanical performance
Kuhlmann-Wilsdorf (1989).

At this juncture, it is known that the electrical conductivity of the Cu-G
composites are dependent on the quantity of graphene used in the composites,
the type of graphene and their defect density, and the arrangement of Cu
matrix and graphene in the composite microstructure, which eventually dictates
the carrier transport in the material during conduction Kappagantula  _et al._
(2022); Pan  _et al._ (2022). Processing routes adopted to manufacture the
composites are also known to influence electrical performance. Processes that
can control atomic scale deposition of Cu and graphene, such as vapor
deposition and molecular beam epitaxy, have shown promise in manufacturing
nano-to-micron scale samples demonstrating enhanced electrical performance
Zheng  _et al._ (2018). At the bulk scale, solid phase processing techniques
have shown promise in enhancing electrical conductivity Kappagantula  _et al._
(2022).

In all this, however, while there is an evolving understanding of processing
approaches and material chemistries that result in Cu-G composites with
improved electrical performance, little is known about the nature of charge-
carrier transfer between Cu and graphene. Several papers discuss carrier
mobility in graphene that is deposited on Cu foils Orofeo  _et al._ (2012);
Banszerus  _et al._ (2015). However, limited experimental and computational
research is published on transport phenomena across metal/graphene interfaces.
There is minimal understanding on how the nature of the interface (such as
atom arrangement, interfacial distance) affects the electron density at the
interface. In Cu-G composites demonstrating enhanced conductivity compared to
the corresponding Cu substrates, two causes were hypothesized for the
enhancement in conductivity. In one case, where the manufacturing conditions
allow for it, the Cu grains are surmised to be templated on the adjacent
graphene flakes leading to the formation of predominantly 111 textures in the
composite during the manufacturing process. Such crystallographic orientations
in Cu show higher conductivity compared to other textures such as 200 or 110.
If there is a large concentration of 111 grains in the Cu microstructure, it
is reasonable to expect an enhanced conductivity overall. In this hypothesis,
graphene is thought to only template the grains but not actually participate
in the conduction process. In the second case, graphene is assumed to
participate during conduction through the exchange of carriers with the
surrounding Cu matrix. One key consideration for this hypothesis is the nature
of interface between the Cu and graphene that may be suitable for carrier
transport across the interface. Some papers indicate that the 111 Cu grains
offer the least lattice mismatch with graphene, which may also help with
carrier transport. However, it is currently not well known how much
conductivity the introduction of graphene brings about for Cu per this
hypothesis.

In this Letter, we explore the electrical transport properties of the Cu-G
composites as a function of the interfacial distance between Cu and graphene,
a parameter that can be modulated during Cu-G manufacturing and a simple proxy
for a possible strained local configuration. We quantify the electronic
transport by computing the electronic conductivity using the Kubo–Greenwood
formula (KGF) Kubo (1957); Greenwood (1958), ideally suited for density-
functional simulations of materials with its single particle Kohn–Sham
orbitals and energies Martin (2008). The KGF method has been utilized in
tight-binding and DFT computations of the conductivity of liquids and solids
Allen and Broughton (1987); Galli  _et al._ (1989). Besides quantifying the
conductivity, we determine the conduction-active sites in the Cu-G composites
by projecting the KGF conductivity onto real space using the space-projected
conductivity (SPC), described elsewhere Prasai  _et al._ (2018). The method
has been implemented to study conduction processes in crystalline systems with
defects, amorphous and semi-conducting systems at atomistic level Subedi  _et
al._ (2022); Thapa  _et al._ (2022); Subedi  _et al._ (2019); Subedi, Prasai,
and Drabold (2020).

We simulated the Cu-G composite by constructing an interface model with a
geometry of the form copper/graphene/copper (Cu/G/Cu) as shown in Figure 1a.
We adopted one of the low-indexed surfaces, namely 111, to construct the Cu
surface. We placed the graphene sheet with one sub-lattice site (A) above the
top of the first layer of Cu 111 and the other sub-lattice site (B) placed
above the third layer of Cu 111 as shown in Figure 1b. This arrangement is
also known as top-fcc, and has shown to represent the low-energy structure for
Cu-G surface models Yang  _et al._ (2017). We simulated the external pressure
on Cu-G composites by varying the copper-graphene distance ($d_{Cu-G}$) of the
Cu/G/Cu interface model as shown in Figure 1a. Besides Cu/G/Cu models, we also
constructed an orthorhombic supercell of 216 atoms and simulated the external
pressure by reducing the vertical dimension of the supercell. A similar
approach has been used to simulate copper and aluminum under external pressure
in earlier work Lanzillo  _et al._ (2014).

Figure 1:  (a) Geometry of copper/graphene/copper (Cu/G/Cu) model with dCu-G
representing copper-graphene distance. b) Top-fcc configuration of graphene
above the Cu 111 surface with A and B representing sub-lattice sites above the
Cu atom at first and third layer respectively. The pink colored spheres
represent Cu atoms and brown colored spheres represent C atoms.

Density-functional calculations were carried out using the Vienna Ab Initio
Simulation Package Kresse and Hafner (1993) (VASP) code. For atomic structure
relaxations of Cu/G/Cu models, we used a plane-wave basis set with a kinetic
energy cutoff of 400 eV. For static calculations, we used a larger cutoff of
520 eV. Projected augmented wave (PAW) Blöchl (1994) potentials was used to
account for ion-election interactions, and the generalized gradient
approximation of Perdew-Burke-Ernzerhof (PBE) Perdew, Burke, and Ernzerhof
(1997) as the exchange-correlation functional. All calculations were performed
at 4 irreducible k-points, and periodic boundary conditions were implemented
throughout. To determine the partial occupancy of electrons near the Fermi
level, Fermi-Dirac distribution function with a smearing temperature of 1000 K
was considered. $\delta$ function in the expression of KGF by was approximated
by a Gaussian distribution function of width 0.05 eV. Factors like thermal
broadening, k-point sampling, and finite size effects discussed in references
Belmonte  _et al._ (2015); Knyazev and levashov (2013); Calder’m, Karaseiv,
and Trickey (2017); Bulanchuk (2021) play a role in employing the KGF.

To obtain representative low-energy configurations, atomic structure
relaxation of the Cu/G/Cu models were performed. The initial and final Cu-G
distances after relaxation are tabulated in the first and the second column of
table 1 respectively. It is apparent that the extent of repulsion of Cu atoms
away from the graphene is higher for models with short Cu-G distances,
indicated by larger shift in Cu-G layers (third column of table 1). For long
Cu-G distances (up to $d_{Cu-G}$ = 2.61 Å ), most of the Cu-Cu bonds are
virtually unaffected and only few Cu-Cu bonds are formed up to $\approx$2.52
Å, whereas for short Cu-G distances, the vertical Cu-Cu bonds are also formed
at $\approx$2.49 and $\approx$2.46 Å.

Table 1:  Variation of Cu-G distance for different Cu/G/Cu models after atomic structure relaxation and corresponding transverse component of Stress Tensor.  $d_{Cu-G}$111$d_{Cu-G}$ is copper-graphene distance shown in Figure 1a)  (Å) | $d_{Cu-G}$ (Å) | Shift (Å) | Pz222Transverse component of Stress Tensor [GPa]
---|---|---|---
(Initial) | (Relaxed) |  |
2.88 | 2.94 | 0.06 | 2.7
2.71 | 2.78 | 0.07 | 3.2
2.54 | 2.61 | 0.07 | 3.1
2.37 | 2.46 | 0.09 | 4.2
2.20 | 2.37 | 0.17 | 6.5
2.02 | 2.29 | 0.27 | 10.3
1.86 | 2.23 | 0.37 | 14.9
1.69 | 2.19 | 0.50 | 20.2

KGF conductivity for the Cu/G/Cu models in the direction normal to the
graphene layer (i.e., along c-axis in Figure 1a) was computed. Figure 2a
displays conductivities for the models with different Cu-G distances
normalized by the KGF conductivity of copper crystal computed at 300 K
corresponding to $d_{Cu-Cu}$ = 2.56 Å. From Figure 2a, we see that the
conductivity increases with decreasing Cu-G distance in an almost exponential
manner up to $d_{Cu-G}$ = 2.37 Å. Beyond 2.37 Å, we find a sharp rise in
conductivity followed by a saturation below $\approx$2.23 Å. The increase in
conductivity may be attributed to increased constructive overlapping of copper
and graphene orbitals with decreasing Cu-G distance. For the short Cu-G
distances with $d_{Cu-G}$ = 2.23 Å and beyond, the conductivity is obtained to
be $\approx$10 % higher compared to that of copper computed at 300 K.

Figure 2:  a) KGF conductivities for relaxed Cu/G/Cu models with different
Cu-G distance normalized by the KGF conductivity ($\sigma_{0}$) of a copper
crystal computed at 300K with $d_{Cu-Cu}$ = 2.56 Å. b) Relative conductivity
($\sigma/\sigma_{0}$) of orthorhombic Cu model computed as a function of
vertical Cu-Cu bond length. Figure 3: Electronic density of states for
different Cu/G/Cu models. Top subplot corresponds to the total EDOS and the
bottom subplot corresponds to the projected density of states (PDOS) onto
carbon atoms. The compression of models resulted in an enhanced EDOS near the
Fermi level, see inset for details. Figure 4: Space-projected conductivity
plotted as a RGB colormap along a plane normal to the graphene sheet for
different Cu/G/Cu models corresponding to $d_{Cu-G}$ = 3.25 Å, 2.94 Å, 2.29 Å
and 2.23 Å respectively. The colorbar on the left represents the magnitude of
SPC values that increases from red towards blue. The SPC values are scaled
from the minimum value for $d_{Cu-G}$ = 3.25 Å. Cu and C atoms are represented
by pink and brown-colored spheres respectively.

To understand such an increasing trend in the conductivity with decreasing
Cu-G distance, total electronic density of states (EDOS) and the projected
density of states (PDOS) for selected models were computed and are shown in
Figure 3. From Figure 3, we find relatively higher EDOS near the Fermi level
with decreasing Cu-G distance. Such increase in the EDOS is contributed by
both copper and carbon atoms of the Cu-G composites, which are more pronounced
for short Cu-G distances. According to Mott and Davis Mott and Davis (1979),
the DC conductivity, $\sigma_{dc}\propto(N(\epsilon_{f}))^{2}$, where $N$
represents the density of states. So, one expects increasing conductivity with
decreasing Cu-G distance for these composites as a combined effect of enhanced
EDOS near the Fermi level. We note that it was not at all obvious that the
graphene would enhance the DOS at the Fermi level, and we conjecture that this
is what underlies experiments showing improved conductivity of the composites.

Further exploration was made by investigating a spatial description of
conductivity as a function of Cu-G distance. We computed the SPC for selected
Cu/G/Cu models and are shown in Figure 4. Figure 4 displays the SPC projected
on 100 plane as a RGB colormap for four such models with different Cu-G
distance. The colorbar on the left depicts the magnitude of SPC values, where
red color represents low SPC values and blue color represents high SPC values.
From Figure 4, it is apparent that the contribution from both Cu and graphene
to the conduction increases with decreasing Cu-G distance, and is more
significant at short distances. The SPC plots at short Cu-G distances show
that the graphene directly participates in conduction and forms a bridge for
conduction between Cu atoms on opposite layers. In addition to the graphene,
we also see significant contributions from Cu atoms to conduction at short
Cu-G distances. This contribution from Cu atoms to conductivity can be
associated with different Cu bonding environments in the Cu-G composites.

KGF conductivity for the orthorhombic Cu model for different vertical Cu-Cu
distances ranging from 2.43 to 2.55 Åwere calculated. Figure 2b shows the
relative conductivity of Cu models for different Cu-Cu vertical distances. One
observes that the conductivity of Cu increases with reducing the Cu-Cu
distance from $\approx$2.55 to $\approx$2.46 Å after which it drops at
$\approx$2.44 Å. We observed that for short Cu-G distances ($d_{Cu-G}\leq$
2.23 Å), the Cu-Cu bond-lengths are $\approx$2.49 and $\approx$2.46 Å. So, the
contribution to the conductivity from the Cu atoms with nearest neighbors
close to 2.46 and 2.49 Å is higher compared with Cu atoms at other distances
and, in agreement with SPC plots in Figure 4.

Saturation in conductivity of Cu/G/Cu model for compression below $d_{Cu-G}$ =
2.23 Å (refer Figure 2a) were observed. Such a characteristic is attributed to
the contribution of the Cu atoms forming different nearest-neighbor distances.
At $d_{Cu-G}$ = 2.23 Å, the repulsion of Cu layers on both sides of graphene
leads to formation of Cu-Cu bonds at $\approx$2.49 and $\approx$2.46 Å. With
further reduction of Cu-G distance, here $d_{Cu-G}$ = 2.19 Å, the Cu-Cu bonds
are formed at $\approx$ 2.46 Å and also at $\approx$2.44 Å. So, the additional
disorder in Cu leads to slightly less contribution from the Cu atoms to the
conduction (refer Figure 2b) even as the graphene still contributes.

In conclusion, copper-graphene composites were modeled by considering the
Cu/G/Cu interface model with varying Cu-G distance. We showed that the DC
conductivity of the Cu/G/Cu interface models increases with decreasing copper-
graphene distance, and such a increase is a combined contributions from both
Cu atoms and the graphene sheet. We attribute this to an enhanced EDOS near
the Fermi-level. We also provided the spatial view of the electronic
conductivity by computing the SPC for varying copper-graphene distance. The
SPC calculations showed that graphene forms the bridge to electronic
conduction between Cu atoms that lie on the opposite layers. We also showed an
interesting characteristic of saturation of the conductivity at a short
copper-graphene distance (below $d_{Cu-G}$ = 2.23 Å), and such a
characteristic is mostly attributed to slightly reduced contribution to the
conduction from Cu atoms that form Cu-Cu bonds at $\approx$2.44 Å even when
the graphene showing enhanced contribution.

This paper provides one of the first demonstrations of graphene participating
in the conduction process along with the copper atoms in the composite,
especially at low copper-graphene distances. An implication of this finding is
that manufacturing process conditions may influence the electrical performance
metal-graphene composites by engineering the interfaces to encourage
graphene’s participation in the conduction processes, rather than leaving it
as a scattering site. It is important to note that in addition to the
interfacial distance, the arrangement of copper-carbon atoms (geometry) can
also influence electrical performance of the composite; however, those effects
were not explored in this work. Additional factors of importance, such as the
number of graphene layers, their dimensions and the defect density of the
graphene additives will be explored in future work.

Acknowledgements: We acknowledge Extreme Science and Engineering Discovery
Environment (XSEDE), supported by NSF Grant No. ACI-1548562 at the Pittsburgh
Supercomputing Center, for providing computational resources under allocation
DMR-190008P,Department of Energy Advanced Manufacturing Office, and Pacific
Northwest National Laboratory operated by the Battelle Memorial Institute for
the U.S. Department of Energy under Contract No. DE-AC06-76LO1830.

## References

  * Heringhaus and Raabe (1996) F. Heringhaus and D. Raabe, “Recent advances in the manufacturing of copper-base composites,” Journal of Materials Processing Technology 59, 367–372 (1996).
  * Bouaziz, Kim, and Estrin (2013) O. Bouaziz, H. S. Kim, and Y. Estrin, “Architecturing of metal-based composites with concurrent nanostructuring: A new paradigm of materials design,” Advanced Engineering Materials 15, 336–340 (2013).
  * Bruschi  _et al._ (2021) S. Bruschi, J. Cao, M. Merklein, and J. Yanagimoto, “Forming of metal-based composite parts,” CIRP Annals 70, 567–588 (2021).
  * Singh and Gautam (2018) M. K. Singh and R. K. Gautam, “Mechanical and tribological properties of plastically deformed copper metal matrix nano composite,” Materials Today: Proceedings 5, 5727–5736 (2018).
  * Weiping  _et al._ (2015) L. Weiping, L. Delong, F. Qiang, and P. Chunxu, “Conductive enhancement of copper/graphene composites based on a high-quality graphene,” RSC Advances 5, 80428–80433 (2015).
  * Rho  _et al._ (2017) H. Rho, Y. S. Jang, S. Kim, S. Bae, T.-W. Kim, D. S. Lee, J.-S. Ha, and S. H. Lee, “Porous copper–graphene heterostructures for cooling of electronic devices,” Nanoscale 9, 7565–7569 (2017).
  * An  _et al._ (2019) Z. An, J. Li, A. Kikuchi, Z. Wang, Y. Jiang, and T. Ono, “Mechanically strengthened graphene-cu composite with reduced thermal expansion towards interconnect applications,” Microsystems & Nanoengineering 5, 20 (2019).
  * Yang  _et al._ (2021) J. Yang, Y. He, X. Zhang, W. Yang, Y. Li, X. Li, Q. Chen, X. Chen, K. Du, and Y. Yan, “Improving the electrical conductivity of copper/graphene composites by reducing the interfacial impurities using spark plasma sintering diffusion bonding,” Journal of Materials Research and Technology 15, 3005–3015 (2021).
  * Hwang  _et al._ (2013) J. Hwang, T. Yoon, S. H. Jin, J. Lee, T.-S. Kim, S. H. Hong, and S. Jeon, “Enhanced mechanical properties of graphene/copper nanocomposites using a molecular-level mixing process,” Advanced Materials 25, 6724–6729 (2013).
  * Ellis  _et al._ (1993) T. Ellis, I. Anderson, H. Downing, and J. Verhoeven, “Deformation-processed wire prepared from gas-atomized cu-nb alloy powders,” Metallurgical Transactions A 24, 21–26 (1993).
  * Felicia, Rochiem, and Laia (2018) D. M. Felicia, R. Rochiem, and S. M. Laia, “The effect of silver (ag) addition to mechanical and electrical properties of copper alloy (cu) casting product,” AIP Conference Proceedings 1945, 020075 (2018).
  * Kuhlmann-Wilsdorf (1989) D. Kuhlmann-Wilsdorf, “Theory of plastic deformation: - properties of low energy dislocation structures,” Materials Science and Engineering: A 113, 1–41 (1989).
  * Kappagantula  _et al._ (2022) K. S. Kappagantula, J. A. Smith, A. K. Nittala, and F. F. Kraft, “Macro copper-graphene composites with enhanced electrical conductivity,” Journal of Alloys and Compounds 894, 162477 (2022).
  * Pan  _et al._ (2022) C. Pan, A. P. S. Gaur, M. Lynn, M. P. Olson, G. Ouyang, and J. Cui, “Enhanced electrical conductivity in graphene–copper multilayer composite,” AIP Advances 12, 015310 (2022).
  * Zheng  _et al._ (2018) L. Zheng, H. Zheng, D. Huo, F. Wu, L. Shao, P. Zheng, Y. Jiang, X. Zheng, X. Qiu, Y. Liu, _et al._ , “N-doped graphene-based copper nanocomposite with ultralow electrical resistivity and high thermal conductivity,” Scientific reports 8, 1–7 (2018).
  * Orofeo  _et al._ (2012) C. M. Orofeo, H. Hibino, K. Kawahara, Y. Ogawa, M. Tsuji, K.-i. Ikeda, S. Mizuno, and H. Ago, “Influence of cu metal on the domain structure and carrier mobility in single-layer graphene,” Carbon 50, 2189–2196 (2012).
  * Banszerus  _et al._ (2015) L. Banszerus, M. Schmitz, S. Engels, J. Dauber, M. Oellers, F. Haupt, K. Watanabe, T. Taniguchi, B. Beschoten, and C. Stampfer, “Ultrahigh-mobility graphene devices from chemical vapor deposition on reusable copper,” Science advances 1, e1500222 (2015).
  * Kubo (1957) R. Kubo, “Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems,” J.Phys.Soc.Jpn 12, 570–586 (1957).
  * Greenwood (1958) D. A. Greenwood, “The boltzmann equation in the theory of electrical conduction in metals,” Proc.Phys.Soc. 71, 585–596 (1958).
  * Martin (2008) R. M. Martin, Electronic Structure, Cambridge University Press, Cambridge (2008).
  * Allen and Broughton (1987) P. B. Allen and J. Q. Broughton, J. Phys. Chem. 91, 4964 (1987).
  * Galli  _et al._ (1989) G. Galli, R. M. Martin, R. Car, and M. Parrinello, Phys. Rev. Lett. 63, 988 (1989).
  * Prasai  _et al._ (2018) K. Prasai, K. N. Subedi, K. Ferris, P. Biswas, and D. A. Drabold, “Spatial projection of electronic conductivity: The example of conducting bridge memory materials,” Phys Status solidi:Rapid Res.Lett. 12 (2018).
  * Subedi  _et al._ (2022) K. N. Subedi, K. Kappagantula, F. Kraft, A. Nittala, and D. A. Drabold, “Electrical conduction processes in aluminum: Defects and phonons,” Phys. Rev. B 105, 104114 (2022).
  * Thapa  _et al._ (2022) R. Thapa, C. Ugwumadu, K. Nepal, J. Trembly, and D. A. Drabold, “Ab initio simulation of amorphous graphite,” Phys. Rev. Lett. 128, 236402 (2022).
  * Subedi  _et al._ (2019) K. N. Subedi, K. Prasai, M. N. Kozicki, and D. A. Drabold, “Structural origins of electronic conduction in amorphous copper-doped alumina,” Phys. Rev. Materials 3, 065605 (2019).
  * Subedi, Prasai, and Drabold (2020) K. N. Subedi, K. Prasai, and D. A. Drabold, Phys. Status Solidi B 258, 2000438 (2020).
  * Yang  _et al._ (2017) T. Yang, L. Yang, H. Liu, H. Zhou, S. Peng, X. Zhou, F. Gao, and X. Zu, “Ab initio study of stability and migration of point defects in copper-graphene layered composite,” Journal of Alloys and Compounds 692, 49–58 (2017).
  * Lanzillo  _et al._ (2014) N. A. Lanzillo, J. B. Thomas, B. Watson, M. Washington, and S. K. Nayak, “Pressure-enabled phonon engineering in metals,” Proceedings of the National Academy of Sciences 111, 8712–8716 (2014).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, “Ab initio molecular dynamics for liquid metals,” Phys. Rev. B 47, 558–561 (1993).
  * Blöchl (1994) P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B 50, 17953–17979 (1994).
  * Perdew, Burke, and Ernzerhof (1997) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple [phys. rev. lett. 77, 3865 (1996)],” Phys. Rev. Lett. 78, 1396–1396 (1997).
  * Belmonte  _et al._ (2015) A. Belmonte, U. Celano, A. Redolfi, A. Fantini, R. Muller, W. Vandervorst, M. Houssa, M. Jurczak, and L. Goux, IEEE Trans.Elec.Devices 62, 4349 (2015).
  * Knyazev and levashov (2013) D. Knyazev and P. levashov, Computational Materials Science 79, 817 (2013).
  * Calder’m, Karaseiv, and Trickey (2017) L. Calder’m, V. Karaseiv, and S. Trickey, Computer Physics Communications 221, 118 (2017).
  * Bulanchuk (2021) P. Bulanchuk, Computational Materials Science 261, 107714 (2021).
  * Mott and Davis (1979) N. Mott and E. Davis, “Electronic processes in non-crystalline materials, 2nd edn. clarendon,” Oxford (1979).


# On Shock Profiles in Four-Field Formulations
of Dissipative Relativistic Fluid Dynamics

Heinrich Freistühler1,a and Blake Temple2,b

(March 17, 2023)

###### Abstract

This paper shows that in second-order hyperbolic systems of partial
differential equations proposed in authors’ earlier paper (J. Math. Phys. 59
(2018)) for modelling the relativistic dynamics of barotropic fluids in the
presence of viscosity and heat conduction, shock waves of arbitrary strength
have smooth, monotone dissipation profiles. The results and arguments extend
classical considerations of Weyl (Comm. Pure Appl. Math. 2 (1949)) and Gilbarg
(Amer. J. Math. 73 (1951)) to the relativistic setting.

Affiliations:
1) Department of Mathematics and Statistics, University of Konstanz, 78457
Konstanz, Germany
2) Department of Mathematics, University of California at Davis, Davis CA
95616, USA

a) E-mail<EMAIL_ADDRESS>
b) E-mail<EMAIL_ADDRESS>
##  1 Introduction

In the theory of relativity, the state of a barotropic fluid can be described
by a 4-vector $\psi^{\alpha},\alpha=0,1,2,3,$ which, as a function of the
space-time coordinates $x^{\beta},\beta=0,1,2,3,$ is governed by a system of
partial differential equations,

$\displaystyle{\partial\over\partial x^{\beta}}\left(T^{\alpha\beta}+\Delta
T^{\alpha\beta}\right)=0,\quad\alpha=0,1,2,3,$ (1.1)

where in case viscosity and / or heat conduction are active,

$\displaystyle\Delta
T^{\alpha\beta}=-B^{\alpha\beta\gamma\delta}\frac{\partial\psi_{\gamma}}{\partial
x\delta}.$ (1.2)

The tensors

$\displaystyle T^{\alpha\beta}\text{ and
}B^{\alpha\beta\gamma\delta},\quad\alpha,\beta,\gamma,\delta=0,1,2,3,$ (1.3)

are given functions of the four fields, i.e., the components of
$\psi^{\alpha}$. In the absence of viscosity and heat conduction, the
equations of motion reduce to the relativistic Euler equations

$\displaystyle{\partial\over\partial
x^{\beta}}T^{\alpha\beta}=0,\quad\alpha=0,1,2,3.$ (1.4)

The present paper focusses on shock waves, whose ideal version is given by
discontinuous solutions to the latter, (1.4), of the (prototypical) form

$\displaystyle\psi_{\alpha}(x)=\begin{cases}\psi_{\alpha}^{-},&x^{\beta}\xi_{\beta}<0,\\\
\psi_{\alpha}^{+},&x^{\beta}\xi_{\beta}>0,\end{cases}$ (1.5)

and asks whether they can be properly represented in the dissipative setting.
A standard way to achieve such representation is a ‘dissipation profile’,
i.e., a regular solution of (1.1) that depends also only on
$x^{\beta}\xi_{\beta}$ and connects the two states forming the shock, in other
words, a solution $\hat{\psi}$ of the ODE

$\displaystyle\xi_{\beta}\xi_{\delta}B^{\alpha\beta\gamma\delta}(\hat{\psi})\hat{\psi}_{\gamma}^{\prime}=\xi_{\beta}T^{\alpha\beta}(\hat{\psi})-q^{\alpha},\quad
q^{\alpha}:=\xi_{\beta}T^{\alpha\beta}(\psi\pm),$ (1.6)

on $\mathbb{R}$ which is heteroclinic to them,

$\displaystyle\hat{\psi}_{\alpha}(-\infty)=\psi_{\alpha}^{-},\quad\hat{\psi}_{\alpha}(+\infty)=\psi_{\alpha}^{+}.$
(1.7)

Concretely, the state variables of a barotropic fluid are given by

$\psi^{\alpha}=\frac{U^{\alpha}}{\theta},$

where $U^{\alpha}$ and $\theta$ are the 4-velocity and the temperature,

$\theta=\left(-\psi_{\alpha}\psi^{\alpha}\right)^{-1/2},$

the fluid is specified by prescribing its pressure as a function of the
temperature,

$p=\tilde{p}(\theta),$

and the ideal part of the energy-momentum tensor is given by

$\displaystyle
T^{\alpha\beta}=\frac{\partial(\tilde{p}(\theta)\psi^{\beta})}{\partial\psi_{\alpha}}=\theta^{3}\frac{d\tilde{p}(\theta)}{d\theta}\psi^{\alpha}\psi^{\beta}+\tilde{p}(\theta)g^{\alpha\beta}.$
(1.8)

We assume strict causality in the sense that

$\displaystyle\left(\frac{\partial^{2}(\tilde{p}(\theta)\psi^{\beta})}{\partial\psi_{\alpha}\partial\psi_{\gamma}}T_{\beta}\right)_{\alpha,\gamma=0,1,2,3}\text{is
negative definite, all future non-spacelike directions }T_{\beta}.$ (1.9)

As regards the dissipative part we specify, following [8, 9, 10], $-\Delta
T^{\alpha\beta}$ as111We use the metric $g^{\alpha\beta}=(-+++)$ and the
projector $\Pi^{\alpha\beta}=g^{\alpha\beta}+U^{\alpha}U^{\beta}$.

$\displaystyle\begin{aligned} -\Delta
T^{\alpha\beta}_{\Box}&=&\eta\Pi^{\alpha\gamma}\Pi^{\beta\delta}\left[\frac{\partial
U_{\gamma}}{\partial x^{\delta}}+\frac{\partial U_{\delta}}{\partial
x^{\gamma}}-\frac{2}{3}g_{\gamma\delta}\frac{\partial U^{\epsilon}}{\partial
x^{\epsilon}}\right]+\check{\zeta}\Pi^{\alpha\beta}\frac{\partial
U^{\gamma}}{\partial x^{\gamma}}\\\
&&+\sigma\left[U^{\alpha}U^{\beta}\frac{\partial U^{\gamma}}{\partial
x^{\gamma}}-\left(\Pi^{\alpha\gamma}U^{\beta}+\Pi^{\beta\gamma}U^{\alpha}\right)U^{\delta}\frac{\partial
U_{\gamma}}{\partial x^{\delta}}\right]\\\
&&+\chi\bigg{[}\bigg{(}U^{\alpha}\frac{\partial\theta}{\partial
x_{\beta}}+U^{\beta}\frac{\partial\theta}{\partial
x_{\alpha}}\bigg{)}-g^{\alpha\beta}U^{\gamma}\frac{\partial\theta}{\partial
x^{\gamma}}\bigg{]},\end{aligned}$ (1.10)

where

$\displaystyle\sigma=((4/3)\eta+\zeta)/(1-c_{s}^{2})-c_{s}^{2}\chi\theta\quad\text{and}\quad\check{\zeta}=\zeta+c_{s}^{2}\sigma-
c_{s}^{2}(1-c_{s}^{2})\chi\theta$ (1.11)

with $\eta,\zeta,\chi$ the coefficients of shear viscosity, bulk viscosity,
and thermal conductivity, and $0<c_{s}<1$ the speed of sound.

The following are our main results.

######  Theorem 1.

Consider a barotropic fluid with viscosity and without heat conduction. Assume
that the acoustic mode is genuinely nonlinear. Then any Lax shock has a
dissipation profile with respect to $\Delta T_{\Box}$.

######  Theorem 2.

Consider a barotropic fluid with viscosity and with heat conduction. Assume
that the acoustic mode is genuinely nonlinear. Then any Lax shock has a
dissipation profile with respect to $\Delta T_{\Box}$.

Sections 2 and 3 are devoted to the proofs of Theorems 1 and 2, respectively.
In Section 4, we contrast these findings with recently established properties
of a different formulation that was proposed by Bemfica, Disconzi, and Noronha
in [2].

##  2 Shock profiles in viscous barotropic fluids without heat conduction

To demonstrate Theorem 1, we write the pressure also as a function

$p=\hat{p}(\rho)$

of the energy

$\rho=\theta\tilde{p}^{\prime}(\theta)-\tilde{p}(\theta),$

to which the sound speed is of course connected as
$c_{s}^{2}=\hat{p}^{\prime}(\rho)$, and recall from [4, 3, 15] that genuine
nonlinearity of the acoustic mode is characterized by the condition

$\displaystyle(\rho+\hat{p}(\rho))\hat{p}^{\prime\prime}(\rho)+2(1-\hat{p}^{\prime}(\rho))\hat{p}^{\prime}(\rho)>0.$
(2.1)

Consider an ideal shock wave (1.5) and assume w. l. o. g. that the
spatiotemporal direction of propagation is $(0,1,0,0)$, i.e.,
$\xi_{\beta}=\delta_{\beta 1}$, and $\psi_{2}=\psi_{3}=0$ (as can always be
achieved by a Lorentz transformation). The profile ODE system,

$\displaystyle-(\Delta T^{\alpha 1}_{\Box})^{\prime}=T^{\alpha 1}-q^{\alpha},$
(2.2)

then has two active equations, $\alpha=0,1$.

We first characterize the rest points in terms of their dependence on the free
constant $q^{\alpha}$.

######  Lemma 1.

For every $q^{1}>0$ there exists a unique $Q(q^{1})>0$ such that the following
holds: The algebraic system

$T^{\alpha 1}(\cdot)=q^{\alpha},\quad\alpha=0,1,$

has more than one solution if and only if

$\displaystyle q_{1}^{2}<q_{0}^{2}<q_{1}^{2}+Q(q^{1}),$ (2.3)

in which case it has precisely two solutions.

###### Proof.

To see this, note first that the equation $T^{11}(.)=q^{1}$ is equivalent to

$p<q^{1}\quad\text{and}\quad u_{1}^{2}=\frac{q^{1}-p}{\rho+p}.$

Under this condition, the equation $T^{01}(.)=q^{0}$ is equivalent to the
combination of

$\displaystyle(\rho+p)^{2}\frac{q^{1}-p}{\rho+p}\left(1+\frac{q^{1}-p}{\rho+p}\right)=q_{0}^{2}$
(2.4)

and

$\displaystyle u^{1}q^{0}>0.$ (2.5)

We write (2.4) as

$\displaystyle
g(\rho)\equiv-\rho\hat{p}(\rho)+q^{1}(-\hat{p}(\rho)+\rho)=q_{0}^{2}-q_{1}^{2}$
(2.6)

with $g$ defined on the interval

$I\equiv[0,\bar{\rho}]\quad\text{with }\hat{p}(\bar{\rho})=q^{1}.$

Note now that any stationary point of $g$ is a nondegenerate maximum. This
follows as assuming $0=g^{\prime}(\rho)$ at some point $\rho\in I$ implies

$q^{1}=\frac{\hat{p}(\rho)+\rho\hat{p}^{\prime}(\rho)}{1-\hat{p}^{\prime}(\rho)}$

and thus, using (2.1),

$\displaystyle
g^{\prime\prime}(\rho)=-(\rho+q_{1})p^{\prime\prime}(\rho)-2\hat{p}^{\prime}(\rho)=-\frac{\rho+\hat{p}(\rho)}{1-\hat{p}^{\prime}(\rho)}\hat{p}^{\prime\prime}(\rho)-2\hat{p}^{\prime}(\rho)<0.$
(2.7)

As $g^{\prime}(0)>0$, this means that along $I$, $g$ increases from $g(0)=0$
to $Q\equiv\max_{I}g>0$ and then decays to $g(\bar{\rho})=-q_{1}^{2}<0$, the
monotonicity being strict in both parts. Equation (2.6) thus has more than
one, namely two, positive solutions if and only (2.3) holds. ∎

Returning to our shock wave (1.5), we note that it must correspond to certain
parameter values $q^{0},q^{1}$ with the properties recorded in Lemma 1;
regarding (2.5), we fix signs222The opposite case $u^{1},q^{0}<0$ differs only
by a transformation $x^{1}\mapsto-x^{1}$. as

$u^{1},q^{0}>0.$

Since

$-(\Delta T^{\alpha 1}_{\Box})^{\prime}=\sigma(u^{\alpha})^{\prime},$

we can rewrite the traveling wave system (2.2) equivalently as

$\displaystyle\sigma u_{\alpha}(u^{\alpha})^{\prime}$ $\displaystyle=$
$\displaystyle u_{\alpha}(T^{\alpha 1}-q^{\alpha})$ (2.8) $\displaystyle\sigma
v_{\alpha}(u^{\alpha})^{\prime}$ $\displaystyle=$ $\displaystyle
v_{\alpha}(T^{\alpha 1}-q^{\alpha})$ (2.9)

with $(v^{0},v^{1})=(u^{1},u^{0})$ orthogonal to $(u^{0},u^{1})$. As
$u_{\alpha}(u^{\alpha})^{\prime}=0$, equation (2.8) is an algebraic
constraint,

$\displaystyle q^{0}u^{0}-(\rho+q^{1})u^{1}=0,$ (2.10)

or

$\displaystyle
u_{1}=U(\rho):=\left(\left(\frac{\rho+q^{1}}{q^{0}}\right)^{2}-1\right)^{-1/2}$
(2.11)

By virtue of (2.10),

$\displaystyle
u^{0}v_{\alpha}(u^{\alpha})^{\prime}=u^{0}(v^{0}u_{0}^{\prime}+v^{1}u_{1}^{\prime})=u_{1}^{\prime},$
(2.12)

and

$\displaystyle u^{0}v_{\alpha}(T^{\alpha
1}-q^{\alpha})=u^{0}v_{\alpha}((\rho+p)u^{\alpha}u^{1}+pg^{\alpha
1}-q^{\alpha})=u^{0}(-v_{0}q^{0}+v_{1}(p-q^{1})),$ (2.13)

equation (2.9) reduces to

$\sigma u_{1}^{\prime}=(\rho+p)u_{1}^{2}+(p-q^{1})$

or, using (2.10) again,

$\sigma
U^{\prime}(\rho)\rho^{\prime}=R(\rho):=(\rho+p)\left(\left(\frac{\rho+q^{1}}{q^{0}}\right)^{2}-1\right)^{-1}+(p-q^{1})=(\rho+p)\left(\frac{\rho+q^{1}}{q^{0}}\right)^{2}-(\rho+q^{1}).$

As $R>0$ between its two zeros, the heteroclinic solution connects them, and
$\rho$ increases in the direction in which the fluid moves. The latter
correctly fits the fact that Lax shocks are compressive.

_Remarks._ (i) The argument, notably inequality (2.7), reveals the geometric
meaning of the genuine nonlinearity condition (2.1) for the Rankine-Hugoniot
relations.
(ii) For the special case of pure radiation, $p=\rho/3$, the result updates
considerations in [8] to the dissipation tensor $\Delta T_{\Box}$ as derived
in [10], which (cf. a remark towards the end of Sec. 1 in [10]) is not exactly
identical with the version originally proposed in [8].

##  3 Shock profiles in viscous barotropic fluids with heat conduction

To include both viscosity and heat conduction, we have to use the full
dissipation tensor (1.10). In this case, we can work directly with (1.6); in
other words, we express (2.2) via (1.2) as

$\displaystyle B^{\alpha 1\gamma
1}(\hat{\psi})\hat{\psi}_{\gamma}^{\prime}=T^{\alpha
1}(\hat{\psi})-q^{\alpha},\quad q^{\alpha}:=T^{\alpha 1}(\psi^{\pm})$ (3.1)

and consider this system on its natural domain of definition,

$\Psi=\\{(\psi^{0},\psi^{1})\in\mathbb{R}^{2}:\psi^{0}>|\psi^{1}|\\},$

with $q^{\alpha}$ corresponding to a fixed Lax shock $\psi^{-}\to\psi^{+}$ as
in Sec. 2, the characteristic speeds $\lambda_{1,2}$ being both positive at
$\psi^{-}$ and having different signs at $\psi^{+}$.

Since

$T^{\alpha\beta}=\theta^{3}p^{\prime}(\theta)\psi^{\alpha}\psi^{\beta}+p(\theta)g^{\alpha\beta}=\frac{\partial
L^{\beta}(\psi)}{\partial\psi_{\alpha}}$

with

$L^{\beta}(\psi)\equiv p(-(\psi^{\gamma}\psi_{\gamma})^{-1/2})\psi^{\beta},$

and due to sharp causality $B^{\alpha 1\gamma 1}$ is positive definite, and
$L$ with

$L(\psi)=L^{1}(\psi)-q^{\gamma}\psi_{\gamma}$

is a strict Liapunov function for (3.1).

The Jacobian of the vector field $F^{\alpha}=T^{\alpha 1}-q^{\alpha}$ is the
Hessian of $L$,

$\frac{\partial
F^{\alpha}}{\partial\psi_{\beta}}=\frac{\partial^{2}L^{1}(\psi)}{\partial\psi_{\alpha}\partial\psi_{\beta}}=\frac{\partial^{2}L(\psi)}{\partial\psi_{\alpha}\partial\psi_{\beta}}=:H^{\alpha\beta}(\psi),$

and the eigenvalues of $H^{\alpha\beta}(\psi)$ relative to the positive
definite matrix

$\frac{\partial T^{\alpha
0}}{\partial\psi_{\beta}}=\frac{\partial^{2}L^{0}(\psi)}{\partial\psi_{\alpha}\partial\psi_{\beta}}>0$

are the characteristic speeds $\lambda_{1,2}$ at the fluid state $\psi$. If we
choose $q^{\alpha}$ and signs as in Sec. 2, these speeds are both positive at
the left hand state $\psi^{-}$, and of different signs at the right hand state
$\psi^{+}$ of the shock. This implies that

$H^{\alpha\beta}(\psi^{-})>0\quad\text{and}\quad\det
H^{\alpha\beta}(\psi^{+})<0;$

the exactly two critical points of the Liapunov function $L$ thus are a strict
local minimum at $\psi^{-}$ and a hyperbolic saddle at $\psi^{+}$.

For a shock of sufficiently small amplitude, i.e., for a parameter value
$q=(q^{0},q^{1})$ with $q^{1}>0$ while $q_{0}^{2}-q_{1}^{2}>0$ is small
enough, the shock profile exists, i.e., system (3.1) possesses a heteroclinic
orbit with $\alpha$-limit $\psi^{-}$ and $\omega$-limit $\psi^{+}$.
Consequently, for this value of $q$

$c^{+}\equiv L(\psi^{+})>L(\psi^{-})\equiv c^{-}.$

and the $c^{+}$ level line of $L$ contains a closed curve (with a corner at
$\psi^{+}$) whose interior

$\Omega\subset L^{-1}((-\infty,c^{+}))\cap W^{u}(\psi^{-})$

contains $\psi^{-}$. Now, on the one hand, as

$\frac{\partial
L}{\partial\psi_{1}}(\psi^{0},\psi^{1})=p(\theta)+\theta^{3}p^{\prime}(\theta)\psi_{1}^{2}-q^{1}>0$

for sufficiently small

$\theta^{-2}=-\psi^{\gamma}\psi_{\gamma}<\psi_{0}^{2},$

$L$ increases strictly on the line segment

$S\equiv\\{(\psi^{0},\psi^{1})\in\Psi:\psi^{0}=\sigma\\}$

if $\sigma>0$ is chosen small enough (in dependence on $q^{1}$). On the other
hand, $L$ tends to $\pm\infty$ near $\partial\Psi$,

$\lim_{\psi^{1}\searrow-\psi^{0}}L(\psi^{0},\psi^{1})=-\infty,\quad\lim_{\psi^{1}\nearrow+\psi^{0}}L(\psi^{0},\psi^{1})=+\infty.$

Therefore the situation, including the shock profile is robust against
perturbations of $q$ within the range given by (2.3). This means that every
Lax shock has a profile.

##  4 Non-existence of shock profiles in the Bemfica-Disconzi-Noronha model

In 1940, Eckart has proposed equations for dissipative relativistic fluid
dynamics [5]. For barotropic fluids, the Eckart theory reads

$\displaystyle{\partial\over\partial x^{\beta}}\left(T^{\alpha\beta}+\Delta
T^{\alpha\beta}_{E}\right)=0,$ (4.1)

with

$\displaystyle-\Delta
T^{\alpha\beta}_{E}=\eta\Pi^{\alpha\gamma}\Pi^{\beta\delta}\left[\frac{\partial
U_{\gamma}}{\partial x^{\delta}}+\frac{\partial U_{\delta}}{\partial
x^{\gamma}}-\frac{2}{3}g_{\gamma\delta}\frac{\partial U^{\epsilon}}{\partial
x^{\epsilon}}\right]+\zeta\Pi^{\alpha\beta}\frac{\partial U^{\gamma}}{\partial
x^{\gamma}}+\chi\left(\Pi^{\alpha\gamma}U^{\beta}+\Pi^{\beta\gamma}U^{\alpha}\right)\frac{\partial\theta}{\partial
x^{\gamma}}$ (4.2)

where, as before, $\eta,\zeta,\chi>0$ denote the fluid’s coefficients of shear
and bulk viscosity, and heat conduction. Sometimes referred to as
‘relativistic Navier-Stokes’, the Eckart equations (4.1) have been fundamental
to the development of fluid dynamics. However, like the classical Navier-
Stokes equations, they do not have the property of finite speed of
propagation. This violates causality, i.e., the principle that speeds of
propagation must not be larger than that of light.

Since physically, Navier-Stokes is a first order theory, in [10] we address
the causality issue by introducing a rigorous notion of when two different
dissipation tensors $\Delta T^{\alpha\beta}$ produce equations which are
equivalent in a first order sense, and show the following.

######  Theorem 3.

(i) Our hyperbolic Navier-Stokes system (1.1) with $\Delta
T^{\alpha\beta}=\Delta T^{\alpha\beta}_{\Box}$ and (1.11) is causal in the
sense that all signal speeds are not larger than that of light. (ii) The
system is first-order equivalent to the Eckart system (4.1).

Here, first-order equivalence is defined in terms of gradient-expansion
transformations between different four-field theories (1.1) of dissipative
fluid dynamics – two systems are equivalent if the change in $\Delta T$ occurs
only at second order in the magnitude of the dissipation coefficients. See
[10, 9] for a precise definition and a formal algebraic setup of first-order
equivalence transformations.

Our theory with $\Delta T^{\alpha\beta}=\Delta T^{\alpha\beta}_{\Box}$ is not
the only formulation (1.1) that is first-oder equivalent to Eckart’s. After
all, the prominent theory given in Landau and Lifshitz [12] also is, while it
is again not causal. In [2], Bemfica, Disconzi, and Noronha have proposed for
the dynamics of the pure radiation fluid, $p(\theta)=\theta^{4}/3$, another
four-field PDE formulation that is first-order equivalent to Eckart’s (with
$\zeta=\chi=0$), namely

$\displaystyle\partial_{\beta}(T^{\alpha\beta}+\Delta
T^{\alpha\beta}_{BDN})=0$ (4.3)

by setting

$\displaystyle-\Delta
T^{\alpha\beta}_{BDN}=B^{\alpha\beta\gamma\delta}_{BDN}\displaystyle{\frac{\partial\psi_{\gamma}}{\partial
x^{\delta}}},$ (4.4)

with

$\displaystyle B^{\alpha\beta\gamma\delta}_{BDN}=\eta
B^{\alpha\beta\gamma\delta}_{E}-\mu B^{\alpha\beta\gamma\delta}_{1}-\nu
B^{\alpha\beta\gamma\delta}_{2}$ (4.5)

where the classical Eckart viscosity tensor

$\displaystyle\begin{aligned}
B^{\alpha\beta\gamma\delta}_{E}&=\Pi^{\alpha\gamma}\Pi^{\beta\delta}+\Pi^{\alpha\delta}\Pi^{\beta\gamma}-\frac{2}{3}\Pi^{\alpha\beta}\Pi^{\gamma\delta}\end{aligned}$
(4.6)

is augmented via

$\displaystyle
B^{\alpha\beta\gamma\delta}_{1}=(3U^{\alpha}U^{\beta}+\Pi^{\alpha\beta})(3U^{\gamma}U^{\delta}+\Pi^{\gamma\delta}),\quad
B^{\alpha\beta\gamma\delta}_{2}=(U^{\alpha}{\Pi^{\beta}}_{\epsilon}+U^{\beta}{\Pi^{\alpha}}_{\epsilon})(U^{\gamma}\Pi^{\delta\epsilon}+U^{\delta}\Pi^{\gamma\epsilon}).$
(4.7)

It was shown in [2] that this formulation, which we will here briefly refer to
as ‘the BDN model’, is indeed causal if and only if, relative to the classical
coefficient $\eta$ of viscosity, the coefficients $\mu$ and $\nu$ of the
“regulators” $B^{\alpha\beta\gamma\delta}_{1},B^{\alpha\beta\gamma\delta}_{2}$
satisfy

$\displaystyle\mu\geq\frac{4}{3}\eta\quad\text{and}\quad\nu\leq\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1}.$
(4.8)

On the other hand the following theorem was proven in [6]:

######  Theorem 4.

If the dissipation coefficients $\eta,\mu,\nu>0$ satisfy the strict causality
condition

$\displaystyle\mu\geq\frac{4}{3}\eta\quad\text{and}\quad\nu<\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1},$
(4.9)

then the BDN model always possesses Lax shocks that do not admit any
dissipation profile.

This contrasts sharply with our Theorems 1 and 2 above.

Theorem 4 would not preclude the possibility that for some _sharply causal_
choice of $\mu$ and $\nu$, i.e., in the case

$\displaystyle\mu\geq\frac{4}{3}\eta\quad\text{and}\quad\nu=\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1}$
(4.10)

(cf. [6]), all Lax shocks do have profiles again. However, Pellhammer has
shown [14]:

######  Theorem 5.

Whatever values $\eta,\mu,\nu>0$ with (4.10) are assumed, the BDN model always
possesses a range of Lax shocks that have either no profile at all or an
oscillatory profile, i.e., a profile whose orbit infinitely spirals around one
of the endstates.

Compared with non-relativistic gas dynamics [11], this property seems exotic.
It would be interesting to know whether oscillatory shock profiles in the BDN
model are dynamically stable.

For positive stability results on shock profiles in hyperbolically regularized
systems of conservation laws see [1].

Conclusion. Our formulation (1.1), (1.10) of relativistic Navier-Stokes is
causal, takes a solvable hyperbolic form, is physically justified by its
leading-order equivalence with Eckart and Landau-Lifshitz, and, like classical
Navier-Stokes, appears phenomenologically correct notably in the sense that it
captures all shock waves consistently.

## References

  * [1] J. Bärlin: Spectral stability of shock profiles for hyperbolically regularized systems of conservation laws, arXiv:2208.12165
  * [2] F. S. Bemfica, M. M. Disconzi, J. Noronha: Causality and existence of solutions of relativistic viscous fluid dynamics with gravity, Phys. Rev. D 98 (2018), 104064.
  * [3] G. Boillat: Sur l’existence et la recherche d’équations de conservation supplémentaires pour les systèmes hyperboliques, C. R. Acad. Sci. Paris Sér. A 278 (1974), 909–912.
  * [4] Y. Choquet-Bruhat: Ondes asymptotiques et approchées pour des systèmes d’équations aux dérivées partielles non linéaires. J. Math. Pures Appl. 48 (1969), 117–158. Cf. p. 149.
  * [5] C. Eckart, The thermodynamics of irreversible processes. 3: Relativistic theory of the simple fluid, Phys. Rev. 58 (1940), 919-924.
  * [6] H. Freistühler: Non-existence and existence of shock profiles in the Bemfica-Disconzi-Noronha model. Phys. Rev. D 103 (2021), 124045.
  * [7] H. Freistühler, M. Reintjes, B. Temple: Decay and subluminality of modes of all wave numbers in the relativistic dynamics of viscous and heat conductive fluids, J. Math. Phys. 62 (2021), 053101.
  * [8] H. Freistühler, B. Temple: Causal dissipation and shock profiles in the relativistic fluid dynamics of pure radiation, Proc. R. Soc. A 470 (2014), 20140055.
  * [9] H. Freistühler, B. Temple: Causal dissipation for the relativistic dynamics of ideal gases, Proc. R. Soc. A 473 (2017), 20160729.
  * [10] H. Freistühler, B. Temple: Causal dissipation in the relativistic dynamics of barotropic fluids, J. Math. Phys. 59 (2018), 063101.
  * [11] D. Gilbarg: The existence and limit behavior of the one-dimensional shock layer, Amer. J. Math. 73 (1951), 256-274.
  * [12] L. D. Landau, E. M. Lifshitz: Fluid Mechanics. Pergamon Press, London 1959, Section 127.
  * [13] A. Majda, R. L. Pego: Stable viscosity matrices for systems of conservation laws, J. Differ. Eqs. 56 (1985), 229 –262.
  * [14] V. Pellhammer: Oscillating shock profiles in relativistic fluid dynamics, arXiv:2301.11743
  * [15] J. Smoller, B. Temple: Global solutions to the relativistic Euler equations, Comm. Math. Phys. 156 (1993), 67–99.
  * [16] H. Weyl: Shock waves in arbitrary fluids, Comm. Pure Appl. Math. 2 (1949), 103–122.