# Axion Production in Pulsar Magnetosphere Gaps

Anirudh Prabhu<EMAIL_ADDRESS>Stanford Institute for Theoretical
Physics, Stanford University, Stanford, California 94305, USA

###### Abstract

Pulsar magnetospheres admit non-stationary vacuum gaps that are characterized
by non-vanishing ${\bf{E}}\cdot{\bf{B}}$. The vacuum gaps play an important
role in plasma production and electromagnetic wave emission. We show that
these gaps generate axions whose energy is set by the gap oscillation
frequency. The density of axions produced in a gap can be several orders of
magnitude greater than the ambient dark matter density. In the strong pulsar
magnetic field, a fraction of these axions may convert to photons, giving rise
to broadband radio signals. We show that dedicated observations of nearby
pulsars with radio telescopes (FAST) and interferometers (SKA) can probe
axion-photon couplings that are a few orders of magnitude lower than current
astrophysical bounds.

††preprint: APS/123-QED

##  I Introduction

Despite the many successes of the Standard Model (SM), a number of
observations point to the need for new physics. Limits on the neutron electric
dipole moment (EDM) [1] suggest a finely tuned cancellation of the topological
$\theta$ parameter in QCD and the overall phase of the quark mass matrix, a
fine-tuning constraint known as the Strong CP problem. In addition, a number
of cosmological observations point to an abundance of dark matter that does
not appear to consist of SM particles. One of the best-motivated solutions to
both of these problems is the QCD axion, which is the pseudo Nambu-Goldstone
boson of the spontaneously broken Peccei-Quinn symmetry [2, 3, 4, 5]. Non-
perturbative effects generate a potential for the QCD axion that dynamically
relaxes the $\theta$ parameter to its CP-conserving value, explaining the
smallness of the neutron (EDM). In addition, QCD axions may be produced in the
correct abundance to account for dark matter [6, 7, 8]. Axions also appear
more broadly in string theory in the compactification of antisymmetric tensor
fields on Calabi-Yau manifolds [9, 10]. The experimental discovery of axions
would be an important step towards understanding physics beyond the SM.

A number of direct detection experiments have been proposed to detect axions
by exploiting their coupling to photons, $\mathcal{L}\supset-
g_{a\gamma\gamma}a(x){{\bf{E}}\cdot{\bf{B}}}$. Some of these experiments make
use of high-$Q$ microwave cavities [11, 12, 13, 14, 15, 16, 17], gapped
toroidal magnets [18], and dielectric haloscopes [19, 20]. In the
aforementioned experiments, axions are assumed to be produced in the early
universe and make up a significant portion of dark matter. Other experiments
do not rely on axions being dark matter. For example, axions may be produced
by the Primakoff process in the Sun [21, 22] or by an oscillating background
${\bf{E}}\cdot{\bf{B}}$ field in the laboratory [23, 24, 25, 26, 27, 28, 29].
When produced in the laboratory by photons with incident ${\bf{E}}$ parallel
to an externally applied ${\bf{B}}$, axions may penetrate optical barriers and
convert back to photons giving the appearance of light shining through walls
(LSW). LSW experiments are not currently competitive with astrophysical bounds
on axions, due in large part to their reliance on two low-probability events,
photon-to-axion conversion and axion-to-photon conversion. Signal power is
additionally suppressed by the limitations on the strength of magnetic fields
that can be produced in the laboratory.

Beyond the laboratory, astrophysical settings are often ideal testbeds for
theories of physics beyond the SM. Neutron stars, in particular, provide
excellent laboratories for axion detection due to their enormous magnetic
fields and dense plasma magnetospheres. Infalling axion dark matter may
resonantly convert to photons in neutron star magnetospheres to produce thin
lines that may be detectable with current and planned radio missions [30, 31,
32, 33, 34, 35, 36]. Another interesting possibility is that neutron stars and
white dwarfs may be factories for axions. In particular, axions may be
thermally produced in the cores of extreme astrophysical objects and converted
to photons in their magnetospheres [37, 38]. In this paper we discuss a new
mechanism for axion production by pulsars. Pulsars, young, rapidly-rotating
neutron stars, are surrounded by a dense plasma of electrons, positrons, and
ions which are thought to be responsible for gamma ray, X-ray, and radio
emission. In response to plasma outflow along open magnetic field lines, the
magnetosphere forms non-stationary vacuum gaps in which
${\bf{E}}\cdot{\bf{B}}\neq 0$ [39, 40]. The oscillating
${\bf{E}}\cdot{\bf{B}}$ in a gap acts as a source for axions, similar to axion
generation in LSW experiments. The energy density of axions produced by a gap
can be several orders of magnitude higher than the corresponding dark matter
energy density, and hence can give rise to very bright radio signals.

The paper is organized as follows. In Section II, we review some basics of
pulsar magnetospheres, including the formation of non-stationary vacuum gaps
for plasma regeneration. In Section III, we present a concrete model for a
gap, based on which we compute the spectrum of the axions produced in that
gap. In Section IV we discuss the conversion of axions to radio photons in the
strong magnetic fields surrounding the pulsar. In Section V we discuss an
observation scheme using existing and planned radio telescopes, with a
discussion of instrumental and astrophysical backgrounds. We conclude this
section with an estimate of the reach of our scheme. We end with some
concluding remarks in Section VI.

##  II  Theory of Pulsar Magnetospheres

Pulsars are highly magnetized compact stars that emit electromagnetic
radiation over a wide range of frequencies. Many pulsars have been observed
due to their coherent, pulsed emission of radio waves. While the mechanism for
coherent radio emission remains a mystery, it is widely believed that the
dense electron-positron plasma that fills the magnetosphere plays an important
role. Early explanations include coherent radio emission by dense, charged
bunches [40] and maser-like two-stream plasma instabilities [41]. The
explanation based on emission by charged bunches has fallen out of favor on
theoretical grounds as a result of the difficulty in explaining the origin and
longevity of these bunches [42] as well as on observational grounds [43]. A
more recent proposal claims that non-uniform plasma generation in the
magnetosphere, to be discussed subsequently, gives rise to electromagnetic
radiation whose power and spectrum are in agreement with observations [44]. We
show that the processes responsible for plasma generation are also very
efficient axion factories.

The presence of plasma is essential not only to the emission of coherent radio
waves, but also to the stability of the magnetosphere [45]. Generation of
electron-positron pairs relies on the formation of charge-starved “gaps” in
the magnetosphere, which arise as a result of plasma outflow, along open field
lines, through the light cylinder. The gap possesses a large component
${\bf{E}}_{\parallel}$, parallel to the magnetic field, that accelerates
charged particles to high energy. Once sufficient charge generation has taken
place, ${\bf{E}}_{\parallel}$ is screened and the gap is fully discharged. As
a result, large amplitude fluctuations of ${\bf{E}}\cdot{\bf{B}}$ occur in the
gap. The formation of gaps and their subsequent discharge has been
demonstrated in numerical simulations [46].

A pulsar is described by its angular velocity ${\bm{\Omega}}$ and magnetic
moment ${\bm{\mu}}$. We limit our discussion to aligned rotators
(${\bm{\mu}}\parallel{\bm{\Omega}}$). While aligned rotators cannot describe
the observed pulsation characteristic of pulsars, they capture the important
details of pulsar dynamics. The surface of a NS is an excellent electric
conductor. Therefore the surface electric field measured in the co-rotating
frame of a NS is zero. In the laboratory frame, electrons redistribute
themselves to counteract the Lorentz force due to the internal magnetic field
yielding the following (equivalent) surface condition

$\displaystyle{\bf{E}}_{\text{int}}+{\bf v}\times{\bf{B}}_{\text{int}}=0$ (1)

where ${\bf v}={\bm{\Omega}}\times{\bf r}$. For a vacuum magnetosphere (1) may
be used to find the external electrostatic potential,

$\displaystyle\Phi(r,\theta)=-{B_{0}\Omega R^{5}\over
3r^{3}}P_{2}\left(\cos\theta\right)$ (2)

where $B_{0}$ is the magnitude of the surface magnetic field, $R$ is the NS
radius, and $P_{n}(x)$ are Legendre polynomials. An important consequence of
(2) is that there exists an external electric field component,
${\bf{E}}_{\parallel}$, parallel to the dipolar magnetic field. It was shown
by Goldreich and Julian [45] that the vacuum magnetosphere solution is
unstable. The electric force on charges due to ${\bf{E}}_{\parallel}$
dominates gravity and strips charges from the NS surface, accelerating them
along magnetic field lines and populating the magnetosphere. Eventually, the
plasma surrounding the NS screens ${\bf{E}}_{\parallel}$ and co-rotates with
the pulsar. The co-rotating plasma charge density is given by [45]

$\displaystyle\rho_{c}={2{\bm{\Omega}}\cdot{\bf{B}}\over
1-\Omega^{2}r^{2}\sin^{2}\theta}.$ (3)

Co-rotation cannot be sustained at a radial distance (defined in cylindrical
coordinates) $r>R_{LC}=c/\Omega$ where $R_{LC}$ is the radius of the light
cylinder. The light cylinder defines two categories of magnetic field lines:
_closed_ and _open_ , with the former being entirely contained within the
light cylinder and the latter exiting it. Particles moving along open field
lines may escape freely to infinity. As a result, the region containing the
open field lines lacks the charges necessary to screen ${\bf{E}}_{\parallel}$,
leading to the formation of gaps. While the existence of pulsar gaps has been
well-established, there is some uncertainty about the mechanism for generating
${\bf{E}}_{\parallel}$ and the location of gaps. Various gap models have been
proposed including the polar-cap [39, 40], slot gap [47, 48], and outer gap
models [49]. In this article we focus on the observable signatures of polar-
cap gap models.

###  II.1 Polar Cap Production

Along open field lines co-rotation cannot be sustained. As a result, the
electric potential deviates from that given by (2) and the plasma that lies on
open field lines decelerates relative to the NS. Additionally, plasma can
leave the light cylinder through open field lines, leading to a violation of
the longitudinal electric field screening condition in a gap right above the
polar cap (the region on the NS surface out of which open field lines flow).
The longitudinal voltage drop across the gap can exceed $10^{12}$ V [50]. The
strong longitudinal electric field in this gap is unstable to $e^{+}e^{-}$
plasma generation. Repeated formation and discharge of the vacuum gap is
responsible for pulsar activity and, as we show below, the production of
axions. In the ensuing discussion we follow the concrete model developed by
Ruderman and Sutherland [40] (hereafter referred to as the RS model)111The RS
model is restricted to anti-aligned pulsars where
${\bm{\Omega}}\cdot{\bf{B}}<0$ and the plasma above the polar cap is
positively charged. The mechanism has been extended to aligned pulsars in, for
example, [51].. In the RS model, positive charges flow out of the light
cylinder along open field lines and cannot be easily replenished by either the
charge-separated magnetosphere or the NS surface, which strongly binds
positively charged ions. Consequently, a gap of thickness $h\ll R$ forms above
the polar cap. The potential drop across the polar gap is $\Delta
V=\Omega^{*}Bh^{2}$ where $\Omega^{*}$ is the rotational velocity of plasma in
the open magnetosphere, which is the same as the rotational velocity of the NS
up to $\mathcal{O}(h^{2}/R^{2})$.

###  II.2 Dynamics of the Gap

The gap width $h$ grows at roughly the speed of light with the potential drop
across the gap growing as $h^{2}$. Gamma rays with energy $\omega\gg 2m_{e}$
may convert to $e^{+}e^{-}$ pairs in the strong gap fields; the $e^{+}e^{-}$
pairs travel along curved field lines and give off curvature radiation which
can further pair-produce in the electric and magnetic fields. Pair-production
proceeds until the gap width exceeds the mean free path of energetic gamma
rays with $\omega>2m_{e}$, at which point the gap breaks down into a series of
sparks. This condition occurs when [40]

$\displaystyle h\approx 360\text{ cm}\left(\rho\over 10\text{
km}\right)^{2/7}\left(\Omega\over\text{Hz}\right)^{-3/7}\left(B_{s}\over
10^{14}\text{ G}\right)^{-4/7},$ (4)

where $\rho$ is the radius of curvature of field lines exiting the polar cap
and $B_{s}$ is the magnitude of the surface magnetic field. The potential drop
across the gap vanishes once $e^{+}e^{-}$ density coincides with (3) and
${\bf{E}}_{\parallel}$ can be effectively screened. To capture these broad
features, we assume ${\bf{E}}\cdot{\bf{B}}$ is periodic in time (with period
$T$) with

$\displaystyle{\bf{E}}\cdot{\bf{B}}({\bf{x}},t)=$ $\displaystyle 2\Omega
B_{s}^{2}h\left\\{\begin{array}[]{ccc}t/h&0\leq z\leq t,&0\leq t\leq h\\\
1&0\leq z\leq h,&h<t\leq T-h\\\ \left(T-t\over h\right)&z\leq 0\leq
T-t,&T-h<t\leq T\end{array}\right.,$ (8)

and the radial distance from the rotation axis is less than the polar cap
radius, $r_{\text{pc}}$. Here $z$ represents the direction normal to the NS
surface at the polar cap. Before discussing axion production in this model, we
highlight a few effects that should be considered in a more sophisticated
model. Firstly, the mean free path of the aforementioned gamma rays depends
exponentially on the magnetic field and gamma ray energy. Therefore, small
fluctuations in the magnetic field can lead to inhomogeneous spark formation
in the gap and a complex sub-burst structure. These fluctuations, however, do
not change either the axion density or frequency spectrum by more than
$\mathcal{O}(1)$ fractional amounts and will thus be neglected in the
following discussion. Secondly, it is unlikely that gap formation and
discharge occurs periodically. More realistically, there will be variations in
the time between successive discharge events. We discuss the relaxation of the
periodicity assumption in the following section.

The periodic time oscillation of the gap produces axions with energy given by
the Fourier modes of the oscillation,

$({\bf{E}}\cdot{\bf{B}})({\bf{x}},t)=\displaystyle\sum_{n}({\bf{E}}\cdot{\bf{B}})_{n}({\bf{x}})e^{i\omega_{n}t}$
(9a) $({\bf{E}}\cdot{\bf{B}})_{n}({\bf{x}})={2\over
T}\displaystyle\int_{0}^{T}({\bf{E}}\cdot{\bf{B}})({\bf{x}},t)e^{i\omega_{n}t}.$
(9b)

where $\omega_{n}=2\pi n/T$. Equation (6) represents the axion spectral
content. For some standard pulsar parameters these frequencies are

$\displaystyle\frac{\omega_{1}}{2\pi}=42\text{
MHz}\left(\frac{B_{s}}{10^{14}\text{
G}}\right)^{4/7}\left(\frac{\Omega}{\text{
Hz}}\right)^{3/7}\left(\frac{cT}{2h}\right)^{-1}.$ (10)

The time-scale of oscillation here is set by the crossing time of the gap at
maximum thickness. Another time-scale which may contribute to axion production
is the plasma frequency near the pulsar surface.

Figure 1: Comparison of axion density produced by gap oscillation to that of
infalling axion dark matter at a distance of 100 km from the NS surface. The
comparison is made for pulsar J1856–3754 at angles $\theta=0^{\circ}$ (blue)
and $\theta=90^{\circ}$ (red) with respect to the rotation axis. The axion
dark matter density infalling onto the NS (black) is larger than the local DM
density as a result of Liouville’s theorem. The axion-photon coupling is taken
to be $g_{a\gamma\gamma}=10^{-11}$ GeV-1.

##  III Axion Production in Gaps

Axions couple to EM fields through the effective Lagrangian,

$\displaystyle\mathcal{L}_{\text{EM}+a}=$ (11) $\displaystyle-{1\over
4}F_{\mu\nu}F^{\mu\nu}+{1\over 2}(\partial_{\mu}a)^{2}-{1\over
2}m_{a}^{2}a^{2}-{g_{a\gamma\gamma}\over 4}aF_{\mu\nu}\tilde{F}^{\mu\nu},$

where $a$ is the axion field, $m_{a}$ is the axion mass, $g_{a\gamma\gamma}$
is the axion-photon coupling, and
$\tilde{F}^{\mu\nu}=\varepsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}/2$ is the
dual field tensor. The dynamics of axion-photon coupling are governed by the
following equations,

$\displaystyle\nabla\cdot{\bf{E}}$
$\displaystyle={\rho\over\epsilon_{0}}-g_{a\gamma\gamma}\nabla
a\cdot{\bf{B}},$ (12) $\displaystyle\nabla\times{\bf{B}}$
$\displaystyle={\partial{\bf{E}}\over\partial
t}+\mu_{0}{\bf{J}}+g_{a\gamma\gamma}\left(\dot{a}{\bf{B}}-\nabla
a\times{\bf{E}}\right)$ (13) $\displaystyle\left(\Box+m_{a}^{2}\right)a$
$\displaystyle=-g_{a\gamma\gamma}{\bf{E}}\cdot{\bf{B}},$ (14)

where $\rho$ and ${\bf{J}}$ are the charge and current densities. The source-
free Maxwell equations remain unaltered by the presence of an axion. According
to (14), a region where ${\bf{E}}\cdot{\bf{B}}\neq 0$, such as a pulsar gap,
can source axions. The production of axions using a background
${\bf{E}}\cdot{\bf{B}}$ component is employed in Light Shining through Walls
(LSW) and related experiments [23, 24, 25, 26, 27], Any Light Particle Search
(ALPS) [52, 53, 54], and the CERN Resonant Weakly Interacting sub-eV Particle
Search (CROWS) [28]. In experiments such as CROWS, axions are produced in a
cavity in which ${\bf{E}}\cdot{\bf{B}}\neq 0$ where ${\bf{E}}$ comes from an
excited cavity mode and ${\bf{B}}$ is an externally applied, static magnetic
field. A single mode of frequency $\omega$ sources an axion field,

$\displaystyle a_{\omega}({\bf{x}},t)$ $\displaystyle=-{g_{a\gamma\gamma}\over
4\pi}e^{i\omega t}\displaystyle\int
d^{3}{\bf{x}}^{\prime}{e^{ik_{a}|{\bf{x}}-{\bf{x}}^{\prime}|}\over|{\bf{x}}-{\bf{x}}^{\prime}|}\left({\bf{E}}\cdot{\bf{B}}\right)_{\omega}({\bf{x}}^{\prime}),$
(15)

where $k_{a}=\sqrt{\omega^{2}-m_{a}^{2}}$ is the axion momentum and
$\left({\bf{E}}\cdot{\bf{B}}\right)_{\omega}$ is the component of
${\bf{E}}\cdot{\bf{B}}$ with frequency $\omega$. For our purposes the integral
is performed over the pulsar gap. In pulsar gaps ${\bf{E}}\cdot{\bf{B}}$ is
time-dependent due to the outflow of plasma through the light cylinder and the
discharge of the gap induced by pair production. The axion density computed
using (14) and (8) is shown in Fig. 1 at two different magnetosphere locations
in pulsar J1856–3754. The appearance of resonance peaks at $m_{a}=\omega_{n}$
is a result of the assumption of gap oscillation periodicity, which is likely
not realized in realistic models. In calculations of the signal power that
follow we do not consider the effects of the resonance peaks shown in Fig. 1.
The total energy released in gap discharge is predominantly in the form of
ultra-relativistic $e^{+}e^{-}$ pairs. The fraction of the gap energy that is
released in the form of axions is $f\sim(g_{a\gamma\gamma}\Omega
B_{s}h^{2})^{2}\ll 1$ and hence the production of axions does not affect
standard pulsar processes. Axions produced in the vacuum gap are also highly
beamed in the direction normal to the polar cap, as demonstrated in Fig. 2.

Figure 2: Comparison of the axion density produced by the gap to the dark
matter density as a function of position in the magnetosphere. Axion emission
is beamed in the direction normal to the polar cap ($z$-direction). The
parameters are taken to be: $m_{a}=10^{-7}$ eV, $g_{a\gamma\gamma}=10^{-11}$
GeV-1, $B_{0}=10^{14}$ Gauss, $P=1$ sec.

##  IV Axion-Photon Conversion

Once produced, axions can penetrate the dense plasma until they reach a
distance at which photon conversion becomes kinematically allowed. In the
presence of a strongly magnetized plasma, axion electrodynamics (12–14) must
be supplemented by strong QED terms,

$\displaystyle\mathcal{L}_{\text{QED}}\supset{\alpha^{2}\over
90m_{e}^{4}}\left[(F_{\mu\nu}F^{\mu\nu})^{2}+{7\over
4}(F_{\mu\nu}\tilde{F}^{\mu\nu})^{2}\right],$ (16)

where $\alpha$ is the electromagnetic fine-structure constant and $m_{e}$ is
the electron mass. Axion-photon conversion results from mixing of the axion
with the (parity-odd) state of the photon, $A_{\parallel}$, whose polarization
is parallel to the external magnetic field, as shown below,

$\displaystyle\left[\partial_{z}^{2}+\omega^{2}+2\omega^{2}\begin{pmatrix}\Delta_{\parallel}(z)&\Delta_{B}(z)\\\
\Delta_{B}(z)&\Delta_{a}\end{pmatrix}\right]\begin{pmatrix}A_{\parallel}\\\
a\end{pmatrix}$ $\displaystyle=0$
$\displaystyle\Delta_{\parallel}=-{\omega_{p}^{2}(z)\sin^{2}\theta\over
2(\omega^{2}-\omega_{p}(z)^{2}\cos^{2}{\theta})}+{14\alpha^{2}\over
45}{B(z)^{2}\sin^{2}\theta\over m_{e}^{4}}$ (17)
$\displaystyle\Delta_{B}={g_{a\gamma\gamma}\omega B(z)\sin{\theta}\over
2(\omega^{2}-\omega_{p}(z)^{2}\cos^{2}{\theta})},\ \Delta_{a}=-{m_{a}^{2}\over
2\omega^{2}},$

where $B$ is the magnitude of the magnetic field, $\theta$ is the angle it
makes with the $z$ direction, $\omega_{p}(z)$ is the spatially varying plasma
frequency, $\omega$ is the frequency of the axion/photon fields. The
refractive index in the direction parallel to the transverse magnetic field,
$n_{\parallel}$ ($\Delta_{\parallel}$ \+ 1) receives a contribution from
strong-field QED. This term may be neglected for non-relativistic axions as it
does not modify the resonance condition even for the highest magnetic fields
observed in magnetars [30], however it must be restored when considering
possibly ultrarelativistic axions converting close to the pulsar surface. We
have ignored the Faraday effect, which can mix the two transverse EM modes
[55].

The conversion probability between axions and photons depends critically on
the relationship between the axion mass and various pulsar parameters. When
the maximal mixing condition,

$\displaystyle{\omega_{p}^{2}\sin^{2}{\theta}\over
1-{\omega_{p}^{2}\over\omega^{2}}\cos^{2}{\theta}}-{28\alpha^{2}\over
45}{B^{2}\sin^{2}\theta\over m_{e}^{4}}\omega^{2}=m_{a}^{2},$ (18)

is satisfied, conversion between axions and photons is resonantly enhanced.
Strong QED effects prevent maximal mixing from being achieved for axions with
frequency $\omega>(360\pi^{2}\Omega^{2}m_{e}^{2}/7e^{2})^{1/4}\approx
1.6\times 10^{-4}$ eV $(\Omega/\text{Hz})^{1/2}$ such as those produced
thermally in the cores of NSs [37] and white dwarfs [38]. At lesser
frequencies, strong QED effects do not contribute significantly to conversion.

The fundamental Fourier mode of gap oscillations is $\omega_{1}\sim 10^{-7}$
eV, which means axions produced by such oscillations can convert to photons
resonantly even if they are ultra-relativistic. We restrict our analysis to
$m_{a}\lesssim 10^{-5}$ eV since for heavier axions the axion dark matter
density dominates that of the axions produced by gap oscillations. For light
axions ($m_{a}<\omega_{p}(R_{LC})$) resonant conversion does not occur in the
magnetosphere. Intermediate mass axions
($\omega_{p}(R_{LC})<m_{a}\ll\omega_{1}$) require some care, as resonant
conversion occurs in the outer magnetosphere where magnetic fields are
considerably weaker than near the NS surface.

The length scale over which the magnetic field and plasma properties change is
generally much larger than the de Broglie wavelength of the axions and photons
generated by the gap. Therefore we may approximate
$\omega^{2}+\partial_{z}^{2}\approx(\omega+k)(\omega-i\partial_{z})$. As
mentioned above, for light axions the axion and photon dispersion relations do
not coincide within the light cylinder. However, since light axions are ultra-
relativistic, both the photon and axion momenta are well-approximated by
$k_{a,\gamma}\approx\omega$. This approximation breaks down when $m_{a}$
becomes comparable to the fundamental oscillation frequency of the gap.
However at this mass the conversion occurs resonantly and hence
$k_{a}=k_{\gamma}=\sqrt{\omega^{2}-m_{a}^{2}}$. In both regimes, we may reduce
(17) to a system of coupled first-order differential equations [55]. For
ultra-relativistic axions the simplified equation is

$\displaystyle\left[i\partial_{z}+\omega+\omega\begin{pmatrix}\Delta_{\parallel}(z)&\Delta_{B}(z)\\\
\Delta_{B}(z)&\Delta_{a}\end{pmatrix}\right]\begin{pmatrix}A_{\parallel}\\\
a\end{pmatrix}=0.$ (19)

In the weak-mixing limit ($\Delta_{B}\ll\Delta_{\parallel}-\Delta_{a}$), (19)
may be solved to first-order using time-dependent perturbation theory, giving
a conversion probability

$\displaystyle
P_{a\to\gamma}=\omega^{2}\left|\displaystyle\int_{z_{\text{min}}}^{\infty}dz^{\prime}\
\Delta_{B}(z^{\prime})e^{i\omega\left(\Delta_{a}z^{\prime}-\int_{z_{\text{min}}}^{z^{\prime}}dz^{\prime\prime}\Delta_{\parallel}(z^{\prime\prime})\right)}\right|^{2}.$
(20)

where $z_{\text{min}}$ is the radius at which the axion energy is equal to the
photon mass, making photon conversion kinematically possible. In the limit of
maximal mixing, conversion occurs in a narrow region surrounding a critical
radius, $z_{c}$ defined by $\omega_{p}(z_{c})=m_{a}$. In this limit the
dispersion relations of the axion and photon coincide, allowing us to make the
following ansatz $A_{\parallel}(z,t)=\bar{A}_{\parallel}(z)e^{-i\omega
t+ikz},\ a(z,t)=\bar{a}(z)e^{-i\omega t+ikz}$. This ansatz, along with the WKB
approximation lead to the following first-order equation [30],

$\displaystyle\left[i\partial_{z}+{\omega^{2}\over
k}\begin{pmatrix}\Delta_{\parallel}(z)-\Delta_{a}&\Delta_{B}(z)\\\
\Delta_{B}(z)&0\end{pmatrix}\right]\begin{pmatrix}A_{\parallel}\\\
a\end{pmatrix}=0.$ (21)

where $k=\sqrt{\omega^{2}-k_{a}^{2}}$. We emphasize that (21) is only valid in
a narrow region around the conversion radius, $z_{c}$. The amplitude of the
photon field once it has left the conversion region is modulated by the
varying plasma density and asymptotically
$A_{\parallel}(\infty)\approx\sqrt{k/\omega}A_{\parallel}(r_{c})$. The final
conversion probability, which includes the modulation of outgoing photons is

$\displaystyle P_{a\to\gamma,\text{res}}=$ $\displaystyle{\omega^{3}\over
k}\left|\displaystyle\int_{z_{\text{min}}}^{\infty}dz^{\prime}\
\Delta_{B}(z^{\prime})e^{i{\omega^{2}\over
k}\left(\Delta_{a}z^{\prime}-\int_{z_{\text{min}}}^{z^{\prime}}dz^{\prime\prime}\Delta_{\parallel}(z^{\prime\prime})\right)}\right|^{2}.$
(22)

Although equations (19) and (21) were derived in different limits, they are
almost identical. The total radio flux seen on Earth is related to the
conversion probability as,

$\displaystyle S={2\over D^{2}}\rho_{a}(z_{c})P_{a\to\gamma}v_{c}z_{c}^{2},$
(23)

where $\rho_{a}(z)$ is the axion density at distance $z$ from the NS surface,
$z_{c}$ is the radius at which conversion takes place, $v_{c}$ is the axion
velocity at the conversion radius, and $D$ is the distance to the pulsar.

##  V Observation and Backgrounds

Periodicity of gap formation and discharge would lead to a series of narrow
axion-induced radio emission lines at integer multiples of
$\omega_{1}=2\pi/T$. In reality, the gap breakdown is both inhomogeneous
within the gap and likely not perfectly periodic, as suggested by the complex
sub-burst structure of pulsar emissions. While the assumption of periodicity,
and hence the calculation above provides a good approximation for the axion
flux produced by gap oscillations, the breakdown of this assumption
necessarily changes the signal search strategy. We therefore adopt a broadband
strategy for signal detection. Most radio telescopes employ a spectrometer
that divides the passband into sub-bands of width
$\delta\nu=\Delta\nu/N_{\nu}$ to gather spectral information about the source.
Radio data will consist of a series power measurements in the $N_{\nu}$
frequency bins and $N_{t}$ time bins each of width $\delta
t=T_{\text{int}}/N_{t}$. The power of the axion-induced signal in frequency
bin $i$ and time bin $j$ is denoted by $S_{\text{ax}}(i,j)$. Under the null
hypothesis the power measurement will receive contributions from astrophysical
and receiver noise sources to be discussed below. To test the hypothesis of an
axion-induced signal we construct a $\chi^{2}$ test statistic,

$\displaystyle\chi^{2}=2\displaystyle\sum_{i=1}^{N_{\nu}}\displaystyle\sum_{j=1}^{N_{t}}{S_{\text{ax}}(i,j)^{2}\over\sigma_{S}(i)^{2}}$
(24)

where $\sigma_{S}$ is the flux density error in each bin, which is related to
the familiar system equivalent flux density (SEFD) by
$\sigma_{S}(\nu)=\text{SEFD}$. The factor of $2$ comes from the number of
linear polarization states being measured.

The relevant quantity to calculate the sensitivity of various radio telescopes
to broadband signals from axion-photon conversion is the flux density
$\Phi=1/(D^{2}\Delta\nu_{\text{sig}})dP/d\Omega$ where $dP/d\Omega$ is the
signal power emitted in the direction of Earth, $D$ is the distance of the
pulsar from Earth, and $\Delta\nu_{\text{sig}}$ is the bandwidth of the
signal. For a given radio telescope or array the minimum detectable flux
density at unit signal-to-noise is

$\displaystyle S_{\text{min}}$
$\displaystyle={\text{SEFD}\over\sqrt{n_{\text{pol}}\Delta\nu_{\text{rec}}T_{\text{int}}}}$
(25)

where SEFD is the system equivalent flux density, $n_{\text{pol}}$ is the
number of polarization states present in the image, $\Delta\nu_{\text{rec}}$
is the telescope bandwidth, and $T_{\text{int}}$ is the total integration
time.

The SEFD is defined as the flux density of a source that would deliver the
same amount of power to the receiver as conventional noise sources. It may be
expressed as SEFD$=2T_{\text{noise}}/A_{\text{eff}}$ where $T_{\text{noise}}$
is the noise temperature, which receives contributions from instrumental as
well as astrophysical sources, and $A_{\text{eff}}$ is the effective area of
the telescope or array. For a single-dish telescope the effective area is
related to the physical area by an $\mathcal{O}(1)$ aperature efficiency term
and for an array with $N\geq 2$ elements, the effective area is greater than
that of a single telescope by a factor of $\sqrt{N(N-1)}$. In the following
subsection we estimate the instrumental and astrophysical backgrounds that
contribute to SEFD.

###  V.1 Backgrounds

Conventional noise sources in radio astronomy include radiometer losses,
galactic and extragalactic radio sources, the CMB, atmospheric emission, and
spillover radiation from terrestrial sources. The total noise temperature at
frequency $\nu$ is then approximately
$T_{\text{noise}}(\nu)=T_{R}(\nu)+T_{\text{astro}}(\nu)$ where $T_{R}$ is the
receiver noise, $T_{\text{astro}}$ is the noise coming from astrophysical
sources unrelated to the NS. The receiver noise depends on the particular
radio telescope(s) being used and is discussed below. Astrophysical noise is
dominated by anisotropic galactic synchrotron emission as demonstrated by
Haslam et al. [56, 57, 58]. The radio sky also receives contributions from an
isotropic extragalactic component. While these sources differ in their spatial
and spectral properties, the distinction is not important for our purposes. A
recent measurement of the absolute sky temperature by the ARCADE 2
collaboration [59] shows a spectrum

$\displaystyle T_{\text{astro}}=\left(24.1\pm 2.1\right)\text{
K}\left(\nu\over 310\text{ MHz}\right)^{-2.599\pm 0.036}.$ (26)

####  V.1.1 Pulsar and Nebula Radio Emission

Ideal targets for observing the effects of axion conversion are nearby pulsars
with high magnetic fields and minimal pulsed radio emission. The absence of
coherent radio emission is often associated with inefficient plasma generation
in the polar cap, and hence inefficient axion production therein. However,
there are simple explanations, consistent with polar cap pair and axion
production, for the lack of observed radio emission from certain pulsars. For
example, radio-quiet pulsars may be beamed away from Earth, but are still
detectable through quasi-isotropic thermal surface emission of soft X-rays
[60]. Additionally, while axion emission is strongly beamed in the direction
normal to the polar cap, coherent radio emission can be generated at a
considerable distance and angle from the polar axis [44]. Other classes of
targets include soft gamma repeaters (SGRs) and anomalous X-ray pulsars
(AXPs), which are highly magnetic pulsars that are widely believed to be
magnetars. While few SGRs and AXPs have been observed to emit pulsed radio
waves [61], many of them are characterized by their X-ray and gamma-ray
emissions, with no observed radio counterpart. For the strong magnetic fields
in SGRs and AXPs, the lack of radio emission may not be incompatible with
efficient charge generation near the polar cap. Once primary pairs are
produced in the gaps of SGRs and AXPs, the third-order QED process
$\gamma\to\gamma\gamma$ may quench further pair production from curvature
radiation and hence coherent radio emission [62]. One potential pitfall of
SGRs and AXPs is that they possess magnetic fields large enough to produce
$e^{+}e^{-}$ pairs without the need for gaps. While the exact plasma
generation mechanism in magnetars remains an open problem, it is possible that
they possess vacuum gaps [63]. In the table below we include the relevant
properties of three candidate stars with ultra-strong magnetic fields and no
observed radio emission, perhaps due to one of the reasons mentioned in the
preceding discussion. SGR properties are taken from [64]. For the selected
pulsars we assume radio emission from these sources is dominated by standard
astrophysical noise and receiver noise of the antennas.

Name, Ref | $B_{s}\text{ (Gauss)}$ | $P_{0}\text{ (sec)}$ | Dist. (kpc)
---|---|---|---
RX J1856.5–3754 [65] | $1.5\times 10^{13}$ | 7.1 | 0.16
SGR 1900+14 [66] | $7.0\times 10^{14}$ | 5.2 | 12.5
SGR 1806–20 [67] | $2.1\times 10^{15}$ | 7.6 | 13

The fundamental gap oscillation frequency of RX J1856–3754 is $\nu_{1}=6$ MHz,
which is undetectable by terrestrial radio telescopes. Thus any detectable
signal is generated by higher order gap oscillation harmonics. The SGRs listed
in the table do not suffer from the same issue as $\nu_{1,1900+14}\sim 60$
MHz, $\nu_{1,1806-20}\sim 100$ MHz. For axion mass $m_{a}$, the highest signal
flux comes from the $n$th harmonic where
$n=\lceil\text{max}(m_{a}/2\pi,\nu_{\text{min}})/\nu_{1}\rceil$ and
$\nu_{\text{min}}$ is the lowest frequency observed by the telescope or
interferometer. We propose detecting the axion-induced signal over an
observing bandwidth equal to the fundamental gap oscillation frequency of the
pulsar being observed, that is $\delta\nu=\nu_{1}$. Thus, in contrast with
axion dark matter searches that look for thin radio lines, we assume the
signal generated at frequency $\nu_{n}$ ($n$ corresponding to the gap
oscillation harmonic generating the signal), is spread out over the separation
between lines. The sensitivity of this proposal, shown in Fig. 3, corresponds
to 100 hours of observation of RX J1856.5–3754 and SGR 1806–20 with the Five-
hundred-meter Aperture Spherical Telescope (FAST) [68].

Figure 3: Projected reach of 100 hours of observation of magnetar SGR 1806–20
(Blue) and pulsar J1856–3754 (Red) with the Five-hundred-meter Aperture
Spherical Telescope (FAST), assuming the line-of-sight is $0^{\circ}$ (solid)
and $90^{\circ}$ (dashed) from the rotation axis (here assumed to coincide
with the direction of the pulsar magnetic dipole moment). The search is
broadband with bandwidth set by the fundamental gap oscillation frequency
$\Delta\nu\approx 100$ MHz (SGR 1806–20), 6 MHz (J1856–3754). Relevant
parameters describing the pulsars are listed in the table.

##  VI Conclusions

In this work we presented a new mechanism for axion production in pulsar
magnetospheres. As shown, oscillating gap regions, where
${\bf{E}}\cdot{\bf{B}}\neq 0$, can source axions (as in Light Shining through
Walls experiments) with energies corresponding to those of the gap oscillation
modes and whose abundance greatly exceeds that of axion dark matter near the
neutron star. The conversion of the created axions to photons gives rise to
enhanced radio emission at the Fourier frequencies of the gap oscillations. We
found that 100 hours of observation of various Galactic pulsars with minimal
radio emission could display an indirect signal of axion-photon conversion
events. We found that radio missions operating at frequencies above 50 MHz,
such as FAST and SKA, can probe axion-photon couplings that are orders of
magnitude lower than astrophysical bounds over a wide range of masses. At high
axion-photon couplings that are still not excluded by current bounds, axion-
photon conversion events may be bright enough to explain the enigmatic fast
radio bursts (FRBs) observed by a number of collaborations [69, 70, 71, 72,
73, 74, 75, 76, 77, 78, 79, 80]. The mechanism we have proposed is consistent
with claims that FRBs have a magnetar origin [81, 82, 83]. A more detailed
analysis of the connection between gap oscillation-sourced axions and FRBs is
the subject of ongoing investigation. The projected reach of this proposal
also motivates the study of more sophisticated models of gap dynamics, such as
those found in recent particle-in-cell simulations [44].

###### Acknowledgements.

The author thanks Savas Dimopoulos for helpful conversations and Robert
Lasenby for useful comments on signal detection and on the manuscript. The
author also thanks Anatoly Spitkovsky and Roger Blandford for detailed
discussion of vacuum gaps and radio emission properties of pulsars and
magnetars. This work was supported by the National Science Foundation under
Grant No. PHYS- 1720397 and the Gordon and Betty Moore Foundation Grant
GBMF7946. The author acknowledges the support of the Fletcher Jones Foundation
and the National Science Foundation (NSF) Graduate Research Fellowship
Program.

## References

  * Baker  _et al._ [2006] C. A. Baker _et al._ , Improved experimental limit on the electric dipole moment of the neutron, Phys. Rev. Lett. 97, 131801 (2006).
  * Peccei and Quinn [1977a] R. D. Peccei and H. R. Quinn, $\mathrm{CP}$ conservation in the presence of pseudoparticles, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn [1977b] R. D. Peccei and H. R. Quinn, Constraints imposed by $\mathrm{CP}$ conservation in the presence of pseudoparticles, Phys. Rev. D 16, 1791 (1977b).
  * Weinberg [1978] S. Weinberg, A new light boson?, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Problem of strong $P$ and $T$ invariance in the presence of instantons, Phys. Rev. Lett. 40, 279 (1978).
  * Abbott and Sikivie [1983] L. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120, 137 (1983).
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the invisible axion, Physics Letters B 120, 127 (1983).
  * Svrcek and Witten [2006] P. Svrcek and E. Witten, Axions in string theory, Journal of High Energy Physics 2006, 051?051 (2006).
  * Arvanitaki  _et al._ [2010] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String axiverse, Phys. Rev. D 81, 123530 (2010).
  * Sikivie [1983] P. Sikivie, Experimental tests of the “invisible” axion, Phys. Rev. Lett. 51, 1415 (1983).
  * Bradley  _et al._ [2003] R. Bradley _et al._ , Microwave cavity searches for dark-matter axions, Rev. Mod. Phys. 75, 777 (2003).
  * Asztalos  _et al._ [2004] S. J. Asztalos _et al._ , Improved rf cavity search for halo axions, Phys. Rev. D 69, 011101 (2004).
  * Asztalos  _et al._ [2010] S. J. Asztalos _et al._ , Squid-based microwave cavity search for dark-matter axions, Physical Review Letters 104, 10.1103/physrevlett.104.041301 (2010).
  * Du  _et al._ [2018] N. Du _et al._ (ADMX Collaboration), Search for invisible axion dark matter with the axion dark matter experiment, Phys. Rev. Lett. 120, 151301 (2018).
  * Brubaker  _et al._ [2017] B. M. Brubaker _et al._ , First results from a microwave cavity axion search at $24\text{ }\text{ }\mu\mathrm{eV}$, Phys. Rev. Lett. 118, 061302 (2017).
  * McAllister  _et al._ [2017] B. T. McAllister, G. Flower, J. Kruger, E. N. Ivanov, M. Goryachev, J. Bourhill, and M. E. Tobar, The organ experiment: An axion haloscope above 15 ghz (2017), arXiv:1706.00209 [physics.ins-det] .
  * Kahn  _et al._ [2016] Y. Kahn, B. R. Safdi, and J. Thaler, Broadband and resonant approaches to axion dark matter detection, Phys. Rev. Lett. 117, 141801 (2016).
  * Caldwell  _et al._ [2017] A. Caldwell _et al._ (MADMAX Working Group), Dielectric haloscopes: A new way to detect axion dark matter, Phys. Rev. Lett. 118, 091801 (2017).
  * Baryakhtar  _et al._ [2018] M. Baryakhtar, J. Huang, and R. Lasenby, Axion and hidden photon dark matter detection with multilayer optical haloscopes, Phys. Rev. D 98, 035006 (2018).
  * Arik  _et al._ [2009] E. Arik _et al._ , Probing ev-scale axions with cast, Journal of Cosmology and Astroparticle Physics 2009 (02), 008?008.
  * Arik  _et al._ [2015] M. Arik _et al._ (CAST Collaboration), New solar axion search using the cern axion solar telescope with ${}^{4}\mathrm{He}$ filling, Phys. Rev. D 92, 021101 (2015).
  * Van Bibber  _et al._ [1987] K. Van Bibber _et al._ , Proposed experiment to produce and detect light pseudoscalars, Phys. Rev. Lett. 59, 759 (1987).
  * Hoogeveen [1992] F. Hoogeveen, Terrestrial axion production and detection using RF cavities, Phys. Lett. B 288, 195 (1992).
  * Sikivie  _et al._ [2007] P. Sikivie, D. B. Tanner, and K. van Bibber, Resonantly enhanced axion-photon regeneration, Phys. Rev. Lett. 98, 172002 (2007).
  * Jaeckel and Ringwald [2008] J. Jaeckel and A. Ringwald, A Cavity Experiment to Search for Hidden Sector Photons, Phys. Lett. B 659, 509 (2008), arXiv:0707.2063 [hep-ph] .
  * Caspers  _et al._ [2009] F. Caspers, J. Jaeckel, and A. Ringwald, Feasibility, engineering aspects and physics reach of microwave cavity experiments searching for hidden photons and axions, JINST 4, P11013, arXiv:0908.0759 [hep-ex] .
  * Betz  _et al._ [2013] M. Betz, F. Caspers, M. Gasior, M. Thumm, and S. W. Rieger, First results of the cern resonant weakly interacting sub-ev particle search (crows), Phys. Rev. D 88, 075014 (2013).
  * Ballou  _et al._ [2015] R. Ballou _et al._ (OSQAR Collaboration), New exclusion limits on scalar and pseudoscalar axionlike particles from light shining through a wall, Phys. Rev. D 92, 092002 (2015).
  * Hook  _et al._ [2018] A. Hook, Y. Kahn, B. R. Safdi, and Z. Sun, Radio signals from axion dark matter conversion in neutron star magnetospheres, Phys. Rev. Lett. 121, 241102 (2018).
  * Safdi  _et al._ [2019] B. R. Safdi, Z. Sun, and A. Y. Chen, Detecting axion dark matter with radio lines from neutron star populations, Phys. Rev. D 99, 123021 (2019).
  * Leroy  _et al._ [2020] M. Leroy, M. Chianese, T. D. P. Edwards, and C. Weniger, Radio Signal of Axion-Photon Conversion in Neutron Stars: A Ray Tracing Analysis, Phys. Rev. D 101, 123003 (2020), arXiv:1912.08815 [hep-ph] .
  * Battye  _et al._ [2020] R. A. Battye _et al._ , Dark matter axion detection in the radio/mm waveband, Phys. Rev. D 102, 023504 (2020).
  * Foster  _et al._ [2020] J. Foster _et al._ , Green bank and effelsberg radio telescope searches for axion dark matter conversion in neutron star magnetospheres, Phys. Rev. Lett. 125, 171301 (2020).
  * Buckley  _et al._ [2021] J. H. Buckley, P. S. B. Dev, F. Ferrer, and F. P. Huang, Fast radio bursts from axion stars moving through pulsar magnetospheres, Phys. Rev. D 103, 043015 (2021).
  * Prabhu and Rapidis [2020] A. Prabhu and N. M. Rapidis, Resonant Conversion of Dark Matter Oscillons in Pulsar Magnetospheres, JCAP 10, 054, arXiv:2005.03700 [astro-ph.CO] .
  * Buschmann  _et al._ [2021] M. Buschmann, R. T. Co, C. Dessert, and B. R. Safdi, X-ray Search for Axions from Nearby Isolated Neutron Stars, Phys. Rev. Lett. 126, 021102 (2021), arXiv:1910.04164 [hep-ph] .
  * Dessert  _et al._ [2019] C. Dessert, A. J. Long, and B. R. Safdi, X-ray signatures of axion conversion in magnetic white dwarf stars, Phys. Rev. Lett. 123, 061104 (2019).
  * Sturrock [1971] P. A. Sturrock, A Model of Pulsars, Astrophys. J. 164, 529 (1971).
  * Ruderman and Sutherland [1975] M. A. Ruderman and P. G. Sutherland, Theory of pulsars: polar gaps, sparks, and coherent microwave radiation., Astrophys. J. 196, 51 (1975).
  * Lyutikov  _et al._ [1999] M. Lyutikov, R. D. Blandford, and G. Machabeli, On the nature of pulsar radio emission, Monthly Notices of the Royal Astronomical Society 305, 338 (1999), https://academic.oup.com/mnras/article-pdf/305/2/338/18630642/305-2-338.pdf .
  * Melrose [1995] D. B. Melrose, The models for radio emission from pulsars - The outstanding issues., Journal of Astrophysics and Astronomy 16, 137 (1995).
  * Lesch  _et al._ [1998] H. Lesch, A. Jessner, M. Kramer, and T. Kunzl, On the possibility of curvature radiation from radio pulsars, A& A 332, L21 (1998).
  * Philippov  _et al._ [2020] A. Philippov, A. Timokhin, and A. Spitkovsky, Origin of Pulsar Radio Emission, Phys. Rev. Lett. 124, 245101 (2020), arXiv:2001.02236 [astro-ph.HE] .
  * Goldreich and Julian [1969] P. Goldreich and W. H. Julian, Pulsar Electrodynamics, Astrophys. J. 157, 869 (1969).
  * Timokhin and Arons [2013] A. N. Timokhin and J. Arons, Current Flow and Pair Creation at Low Altitude in Rotation Powered Pulsars’ Force-Free Magnetospheres: Space-Charge Limited Flow, Mon. Not. Roy. Astron. Soc. 429, 20 (2013), arXiv:1206.5819 [astro-ph.HE] .
  * Arons [1983] J. Arons, Pair creation above pulsar polar caps : geometrical structure and energetics of slot gaps., Astrophys. J. 266, 215 (1983).
  * Muslimov and Harding [2004] A. G. Muslimov and A. K. Harding, High-altitude particle acceleration and radiation in pulsar slot gaps, Astrophys. J. 606, 1143 (2004).
  * Cheng  _et al._ [1986] K. S. Cheng, C. Ho, and M. Ruderman, Energetic Radiation from Rapidly Spinning Pulsars. I. Outer Magnetosphere Gaps, Astrophys. J. 300, 500 (1986).
  * Beskin  _et al._ [1993] V. S. Beskin, A. V. Gurevich, and Y. N. Istomin,  _Physics of the pulsar magnetosphere_ (1993).
  * Gil, J.  _et al._ [2003] Gil, J., Melikidze, G. I., and Geppert, U., Drifting subpulses and inner acceleration regions in radio pulsars*, A&A 407, 315 (2003).
  * Ehret  _et al._ [2009] K. Ehret _et al._ (ALPS), Resonant laser power build-up in ALPS: A ’Light-shining-through-walls’ experiment, Nucl. Instrum. Meth. A 612, 83 (2009), arXiv:0905.4159 [physics.ins-det] .
  * Ehret  _et al._ [2010] K. Ehret _et al._ , New ALPS Results on Hidden-Sector Lightweights, Phys. Lett. B 689, 149 (2010), arXiv:1004.1313 [hep-ex] .
  * Bähre  _et al._ [2013] R. Bähre _et al._ , Any light particle search II —Technical Design Report, JINST 8, T09001, arXiv:1302.5647 [physics.ins-det] .
  * Raffelt and Stodolsky [1988] G. Raffelt and L. Stodolsky, Mixing of the photon with low-mass particles, Phys. Rev. D 37, 1237 (1988).
  * Haslam [1981] C. G. T. o. Haslam, A 408 MHz all-sky continuum survey. I - Observations at southern declinations and for the North Polar region., A&A 100, 209 (1981).
  * Haslam  _et al._ [1970] C. G. T. Haslam, M. J. S. Quigley, and C. J. Salter, A 408 MHz survey of the galactic anticentre region, MNRAS 147, 405 (1970).
  * Haslam  _et al._ [1982] C. G. T. Haslam, C. J. Salter, H. Stoffel, and W. E. Wilson, A 408 MHz all-sky continuum survey. II. The atlas of contour maps., A&A 47, 1 (1982).
  * Fixsen  _et al._ [2011] D. J. Fixsen _et al._ , ARCADE 2 MEASUREMENT OF THE ABSOLUTE SKY BRIGHTNESS AT 3-90 GHz, The Astrophysical Journal 734, 5 (2011).
  * Brazier and Johnston [1999] K. T. S. Brazier and S. Johnston, The implications of radio-quiet neutron stars, MNRAS 305, 671 (1999), arXiv:astro-ph/9803176 [astro-ph] .
  * Camilo  _et al._ [2006] F. Camilo _et al._ , Transient pulsed radio emission from a magnetar, Nature 442, 892?895 (2006).
  * Baring and Harding [1998] M. G. Baring and A. K. Harding, Radio-quiet pulsars with ultrastrong magnetic fields, The Astrophysical Journal 507, L55?L58 (1998).
  * Thompson [2008] C. Thompson, Electrodynamics of magnetars. iv. self?consistent model of the inner accelerator with implications for pulsed radio emission, The Astrophysical Journal 688, 499?526 (2008).
  * Olausen and Kaspi [2014] S. A. Olausen and V. M. Kaspi, The McGill Magnetar Catalog, ApJS 212, 6 (2014), arXiv:1309.4167 [astro-ph.HE] .
  * van Kerkwijk and Kaplan [2008] M. H. van Kerkwijk and D. L. Kaplan, Timing the nearby isolated neutron star rx j1856.5-3754, The Astrophysical Journal 673, L163?L166 (2008).
  * Mereghetti  _et al._ [2006] S. Mereghetti _et al._ , The FirstXMM-Newton Observations of the soft gamma-ray repeater SGR 1900+14, The Astrophysical Journal 653, 1423 (2006).
  * Woods  _et al._ [2007] P. Woods _et al._ , The Prelude to and Aftermath of the Giant Flare of 2004 December 27: Persistent and Pulsed X-Ray Properties of SGR 1806-20 from 1993 to 2005, Astrophys. J. 654, 470 (2007), arXiv:astro-ph/0602402 [astro-ph] .
  * Li and Pan [2016] D. Li and Z. Pan, The five-hundred-meter aperture spherical radio telescope project, Radio Science 51, 1060 (2016).
  * Lorimer  _et al._ [2007] D. R. Lorimer, M. Bailes, M. A. McLaughlin, D. J. Narkevic, and F. Crawford, A bright millisecond radio burst of extragalactic origin, Science 318, 777?780 (2007).
  * Thornton  _et al._ [2013] D. Thornton _et al._ , A population of fast radio bursts at cosmological distances, Science 341, 53?56 (2013).
  * Burke-Spolaor and Bannister [2014] S. Burke-Spolaor and K. W. Bannister, The Galactic Position Dependence of Fast Radio Bursts and the Discovery of FRB011025, Astrophys. J. 792, 19 (2014), arXiv:1407.0400 [astro-ph.HE] .
  * Spitler  _et al._ [2014] L. G. Spitler _et al._ , Fast radio burst discovered in the arecibo pulsar alfa survey, The Astrophysical Journal 790, 101 (2014).
  * Ravi  _et al._ [2015] V. Ravi, R. M. Shannon, and A. Jameson, A fast radio burst in the direction of the carina dwarf spheroidal galaxy, The Astrophysical Journal 799, L5 (2015).
  * Masui  _et al._ [2015] K. Masui _et al._ , Dense magnetized plasma associated with a fast radio burst, Nature 528, 523?525 (2015).
  * Champion [2016] D. J. o. Champion, Five new fast radio bursts from the htru high-latitude survey at parkes: first evidence for two-component bursts, Monthly Notices of the Royal Astronomical Society: Letters 460, L30?L34 (2016).
  * Caleb  _et al._ [2016] M. Caleb _et al._ , Fast radio transient searches with utmost at 843 mhz, Monthly Notices of the Royal Astronomical Society 458, 718?725 (2016).
  * Bannister  _et al._ [2017] K. W. Bannister _et al._ , The detection of an extremely bright fast radio burst in a phased array feed survey, The Astrophysical Journal 841, L12 (2017).
  * Shannon  _et al._ [2018] R. Shannon _et al._ , The dispersion-brightness relation for fast radio bursts from a wide-field survey, Nature 562, 386 (2018).
  * Amiri  _et al._ [2018] M. Amiri _et al._ (CHIME/FRB), The chime fast radio burst project: System overview, The Astrophysical Journal 863, 48 (2018).
  * Amiri  _et al._ [2019] M. Amiri _et al._ (CHIME/FRB), Observations of fast radio bursts at frequencies down to 400 megahertz, Nature 566, 230 (2019), arXiv:1901.04524 [astro-ph.HE] .
  * Bochenek [2020] C. D. o. Bochenek, A fast radio burst associated with a galactic magnetar, Nature 587, 59?62 (2020).
  * Bhandari  _et al._ [2020] S. Bhandari _et al._ , The host galaxies and progenitors of fast radio bursts localized with the australian square kilometre array pathfinder, The Astrophysical Journal 895, L37 (2020).
  * Zhang [2020] B. Zhang, The physical mechanisms of fast radio bursts, Nature 587, 45?53 (2020).


# Identification and prediction of time-varying parameters of COVID-19 model:
a data-driven deep learning approach

Jie Longa, A. Q. M. Khaliqb & K. M. Furatic Jie Long. Email:
<EMAIL_ADDRESS>Q. M. Khaliq. Email<EMAIL_ADDRESS>M. Furati.
Email<EMAIL_ADDRESS>a,bDepartment of Mathematical Sciences, Middle
Tennessee State University, Murfreesboro, TN 37132, USA; bDepartment of
Mathematics and Statistics, King Fahd University of Petroleum and Minerals,
Dhahran 31261, Saudi Arabia

###### Abstract

Data-driven deep learning provides efficient algorithms for parameter
identification of epidemiology models. Unlike the constant parameters, the
complexity of identifying time-varying parameters is largely increased. In
this paper, a variant of physics-informed neural network (PINN) is adopted to
identify the time-varying parameters of the Susceptible-Infectious-Recovered-
Deceased model for the spread of COVID-19 by fitting daily reported cases. The
learned parameters are verified by utilizing an ordinary differential equation
solver to compute the corresponding solutions of this compartmental model. The
effective reproduction number based on these parameters is calculated. Long
Short-Term Memory (LSTM) neural network is employed to predict the future
weekly time-varying parameters. The numerical simulations demonstrate that
PINN combined with LSTM yields accurate and effective results.

###### keywords:

PINN; LSTM; Time-Varying Parameters; SIRD; COVID-19; Deep Neural Network

##  1 Introduction

The novel SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) has
spread all over the world since its discovery at the end of 2019, with
millions of confirmed cases and more than a million deaths [29]. On March 11,
2020, COVID-19 was characterized as a pandemic by the World Health
Organization (WHO) [28]. Due to the huge negative effects of COVID-19,
precautionary measures have been aggressively carried out worldwide, such as
facial masking, contact tracing, social distancing, and some governmental
actions such as lockdowns. Hence, it is significant to analyze the dynamics of
COVID-19 so that the effectiveness of those implemented measures can be
verified.

Epidemiological models provide an efficient tool for determining and
explaining the dynamics of disease transmission. In the early stage, one of
the classical models is the Susceptible-Infectious-Recovery (SIR) model
presented by Kermack and McKendrick in 1927 [13]. This compartmental model
computes the theoretical number of susceptible people that became infected by
the disease and then recovered. Based on the SIR model, many models have been
proposed for different diseases, like SARS and COVID-19 [24, 8, 27]. The
Susceptible-Exposed-Infectious-Removed (SEIR) model is proposed to analyze the
effect of the precautionary measures on the dynamics of the epidemic [2].
Taking quarantine into account, a SEIQR model is investigated analytically and
numerically [19]. The SEIR model has been extended by two extra classes of
populations namely, ’C’, which is the number of cumulative cases, and ’D’ that
is the number of severe, critical, and deceased cases, respectively, to
understand the trends of COVID-19 outbreak in Wuhan, China [15].

The SIRD model is governed by three parameters $\beta$, $\gamma$, and $\mu$.
Commonly, these parameters are assumed as independent of time during the
epidemic. However, the interventions and measures implemented by governments
to control the spread of COVID-19 result in significantly varying parameters
with time. Thus, if parameters are still considered as constants, the
epidemiology model cannot discover the varied dynamic of this disease
transmission that is affected by those efficient precautions. Regarding
parameters as time-series enable the model to react to the changes in real
situations.

Recently, data-driven algorithms were developed to the determination of such
time-varying parameters. Raissi et al. [21] proposed PINN to approximate
parameters, which are constants, of the SIR model using a small data set.
Furthermore, Wang et al. [26] divided the period time into four phases so that
their parameters become piece-wise constants. Also, there are some studies
relating to the time-dependent parameter, like a time-dependent SIR model is
proposed to track their parameter time-series [3]. Similarly, the time-
dependent parameters of the SIR model were identified by utilizing a neural
network to analyze the COVID-19 spread in South Korea [11].

Predicting the trend of COVID-19 transmission is a significant and arduous
topic. A model of COVID-19 in China was developed to predict the cumulative
number of reported cases [16]. By averaging the slope of approximated
parameters of the SIRD model over the last several days, the extrapolated
trends of infected, recovered, deceased, and susceptible cases were obtained
via solving the model with estimated parameters [18]. An artificial
intelligence method was used to predict the pandemic by training on the 2003
SARS data [30]. Recently, several forecasting models including autoregressive
integrated moving average, support vector regression, LSTM, and bidirectional
LSTM (BiLSTM) were proposed to predict the confirmed cases, deaths, and
recoveries [23]. Similarly, Recurrent Neural Network (RNN), LSTM, BiLSTM,
Gated recurrent units (GRUs), and Variational AutoEncoder (VAE) algorithms
were employed to forecast the number of infected and recovered cases [31].

Based on data collected, the SIRD model is used in this study to simulate the
dynamics of the disease. Parameters of this model are regarded as daily and
weekly time-varying. As for daily time-varying parameters, it represents that
parameters would be varied from each day so that they are more similar to real
situations. However, the difficulty of parameter identification increases with
the smaller time interval by utilizing some general regression methods such as
least square methods. PINN is still able to approximate these time-series
parameters with some changes in architecture. Then, they are verified by being
substituted into the SIRD model, and then a numerical algorithm is used to
solve this ordinary differential equation (ODE) system.

When it comes to weekly time-series, we divide the collected data week by week
and employ the interpolation method for each week to gain more data so that
the neural network has enough information to learn. The neural network
implemented here is the PINN, which can estimate parameters to fit the
collected data [20]. When the identification of these weekly time-series
parameters is validated, LSTM is implemented to predict future parameters. In
this way, by examining weekly values instead of daily ones, LSTM can still
accurately predict future parameters. Another reason to predict the parameters
rather than forecast the infectious cases directly is that we could obtain
corresponding reproduction number $R_{0}$. It is an essential threshold to
reflect the disease transmission speed and verify the effectiveness of these
measures in curbing the outbreak of COVID-19.

Figure 1: Flow Chart of Data-Driven Algorithm

In Figure 1, we describe that PINN approximates the parameters of the SIRD
model by training collected data. Then, LSTM is proposed to predict the future
value of parameters. After that, obtained parameters are substituted into the
SIRD model so that we could predict the infectious cases.

Outline of the paper: In Section 2 we introduce the SIRD model and explain
model parameters ($\beta$, $\gamma$, and $\mu$). Following that, we introduce
the effective reproduction number and how to compute it. Deep learning neural
networks, including feedforward neural network (FNN), PINN, and LSTM, are
described briefly in Section 3. In Section 4, the collected data is described,
and the loss function of PINN is formulated. After that, we present LSTM
neural networks for the future prediction of parameters. The simulation
results are discussed in Section 5.

##  2 Epidemiology Model

In this section, we introduce the SIRD model employed in our study. One
promise of the SIRD model is that natural birth and death rates are neglected
or equivalent so that the total population is considered as constant. Then,
the population could be divided into four mutually exclusive groups, which are
susceptible, infected, recovered, and deceased, respectively. Also, the
reproduction number is introduced in this section.

###  2.1 SIRD Model

We consider the SIRD model described by the following system of ordinary
differential equations [5]:

$\displaystyle\begin{aligned} \dot{S}(t)&=-\frac{\beta(t)S(t)I(t)}{N}\\\
\dot{I}(t)&=\frac{\beta(t)S(t)I(t)}{N}-(\gamma(t)+\mu(t))I(t)\\\
\dot{R}(t)&=\gamma(t)I(t)\\\ \dot{D}(t)&=\mu(t)I(t),\end{aligned}$ (1)

where $S(t)$, $I(t)$, $R(t)$, and $D(t)$ are the numbers of susceptible,
infected, recovered, and deceased individuals, respectively, and $N$ is the
total population, $N=S(t)+I(t)+R(t)+D(t)$. Without loss of generality, $N$ is
a large number and difficult to compute in some optimization problems. Thus,
we used the fraction of the population in each, which means that $S(t)$,
$I(t)$, $R(t)$, and $D(t)$ are divided by $N$. In this way, our calculation is
simpler, and we can still adhere to the same dynamic system of the epidemic.

The parameter $\beta(t)$ represents the number of contacts each day for
infected individuals. Besides, there are two essential assumptions. First,
contacts between the infected and uninfected people are sufficient to spread
the disease. Second, the population is mixed homogeneously. Thereby,
$\beta(t)S(t)$ susceptible people are infected by a patient each day on
average, and $\beta(t)S(t)I(t)$ is the number of newly infected people,
$\gamma(t)$ is the rate at which infected cases become recovered. Then, there
are $\gamma I(t)$ infected individuals turned into the third compartment.
Similarly, $\mu(t)$ is the rate of mortality by disease, and $\mu(t)I(t)$
infected individuals are deceased.

###  2.2 Reproduction Number

There is an essential threshold used to get a better understanding of the
transmission rate. This threshold is characterized by the reproduction number,
which is an estimate of newly infected cases caused by an infected individual
[25]. Assuming there is just one person infected in the beginning, then in the
SIRD model, $I_{0}=1$, $S_{0}=N-1$, and no recovery at this moment. When we
utilized the fraction for each compartment, $I_{0}\approx 0$ and $S_{0}\approx
1$. By the second equation of eq. (1), we have:

$\dot{I}(t)=I(t)[\beta(t)S(t)-(\gamma(t)+\mu(t))],$ (2)

where $\beta(t)$, $\gamma(t)$, and $\mu(t)$ are greater than 0. In eq. (2), it
is simple to obtain $\dot{I}>0$ when $\beta(t)S(t)>(\gamma(t)+\mu(t))$, which
means the number of infected cases keeps increasing.

A key idea about the reproduction number $R_{0}$ is that it is the gauge of
the secondary infections at the beginning of this disease invasion [9]. After
the disease invasion, the contact rate $\beta$ decreases because of the
reduction of this susceptible group. Those who are infected or have already
recovered cannot be infected once again. Thereby, another similar threshold,
which is the effective reproduction number $R_{e}$ [22], is used to indicate
whether the disease keeps spreading, and $R_{0}$ is the upper bound of $R_{e}$
[12].

$R_{e}(t)=\frac{\beta(t)S(t)}{\gamma(t)+\mu(t)}<\frac{\beta(0)}{\gamma(0)+\mu(0)}=R_{0}.$
(3)

##  3 Deep Neural Network

In this section, we provide a brief description of the deep learning neural
network architectures utilized in our algorithm for parameter identification
and disease dynamics prediction.

###  3.1 Feedforward Neural Network

The feedforward neural network (FNN) is an artificial neural network with a
simple architecture in which the connections between the nodes do not form a
cycle. FNN is utilized to approximate the target function by applying several
activation functions to input recursively and then output a value or vector
that is close to the target values.

The input layer consists of input neurons, bringing the training data into the
network for further processing by hidden layers. Layers between the input
layer and output layer are hidden layers, which perform linear or nonlinear
transformations of training data to produce an intended output. These
transformations are implemented by activation functions. The first hidden
layer takes the weighted sums of input through an activation function to
produce output for each neuron. Then, these outputs plus a constant, which is
bias in the first hidden layer, are viewed as an input of the second hidden
layer. After that, as the last layer of the network, the output layer presents
desired results [6].

The forward propagation from one layer to the next layer’s nodes is given as
follows [7].

$z_{j}^{l+1}=\sum_{i}^{n_{l}}w^{l}_{ij}f_{l-1}(z_{i}^{l})+b^{l},$ (4)

where $n_{l}$ represents the number of neurons in the $l^{th}$ layer,
$z_{i}^{l}$, $i=1,\cdots,n_{l}$, denotes the output of the $i^{th}$ node in
$(l-1)^{th}$ layer, $f$ represents the activation function, $w_{ij}^{l}$ are
the weights between node $i^{th}$ and node $j^{th}$, $b^{l}$ are the bias in
the $l^{th}$ layer. The activation functions employed in this study are the
hyperbolic tangent function

$\tanh(x)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}},$ (5)

and sigmoid function

$\sigma(x)=\frac{1}{1+e^{-x}}.$ (6)

When output values are obtained from the output layer, they are used to build
a loss function, which is used to estimate the error of the model. An
optimizer, like Adam [14] or gradient descent method [1], is employed to
minimize the loss function so that the output is close to the observations.
According to different problems, the selections of optimizer and loss function
are different.

###  3.2 Physics-Informed Neural Network

The physics-inform neural network (PINN) is a data-driven algorithm to
approximate the solution of differential equations and identify parameters
[20]. It could utilize any type of neural network architectures, like the FNN
and Convolutional neural networks (CNN), as the main framework. The applied
activation functions and optimization methods in PINN are the same as the
usual deep learning techniques. The fascinating part of this algorithm is the
loss function, which is generally comprised of boundary conditions, initial
values, and physical constraints.

The outputs of the neural network are constrained to satisfy the system of
differential equations by penalizing the residuals of differential equations
into the loss function. In these residual equations, the derivatives of
outputs with respect to time are computed by a black-box, which is automatic
differentiation [17]. Automatic differentiation is employed to train these
models and allows PINN to take derivatives with respect to input coordinates
so that it could discover the latent physical laws. More details could be
found in [20].

###  3.3 Long Short-term Memory Neural Network

As an advanced tool to process and predict time-series, LSTM is a special type
of RNN. The difference between RNN and feedforward neural network is that it
is capable of storing past information by taking previous layers’ output as an
input for the next input. Also, all layers of RNN share the same parameters
because they perform the same tasks for the elements of the sequence, which is
a so-called recurrent. Based on this feature, RNN could predict sequential
information. However, on account of only involving the output of the last
layer, it cannot work out long-term memory problems [10]. Thus, LSTMs are
created for solving them.

Figure 2: Architecture of LSTM Cell [4]

LSTMs as shown in Figure 2 are chains of memory cells containing three gates,
which are input gate, forget gate, and output gate as explained in [4]. The
key part of these memory cells is the cell state that could maintain the
global information of sequence in each time-step. The content of the cell
state would be modified at different time-steps through the forget gates and
input gates. More specifically, forget gates decide what information would be
thrown away by utilizing the sigmoid function. Instead, input gates decide
what new information would be stored by combining sigmoid eq. (6) and tanh eq.
(5) functions. Then, the cell state is updated and moves forward to the next
state by implementing some operations. After that, the output is calculated in
the output gate. By analyzing the structure of a memory cell, we could see
that LSTMs are capable of storing long-term memory by putting new information
into cell states.

##  4 Data and Algorithms

In this section, we describe the collected data and introduce our algorithms
for parameter identification and making predictions. In our work, the
parameters of the SIRD model are viewed as daily-varying and weekly-varying
parameters to simulate the complex real situation of COVID-19. Based on this
perspective, two algorithms are employed to identify daily and weekly time-
series parameters. Having the identified parameters by PINN as an input, we
introduce the LSTM algorithm for making predictions.

###  4.1 Data

The data considered in our study is downloaded from
https://covidtracking.com/data/download for New York, New Mexico, and Texas
states. It is comprised of commutative infected, recovered, and deceased cases
for the period from March 30 to September 30. A plot of the data is shown in
Figure 3.

(a)

(b)

(c)

Figure 3: COVID-19 data for New York, New Mexico, and Texas

As can be observed, the three states have different dynamics between the
compartments. In New York, while the number of cases in I, R, and D
compartments is increasing, the number of recovered and deceased cases is
relatively small compared to infections. New Mexico data shows that the
infected and recovered cases kept almost the same increasing rate before
August but then the infections started decreasing slowly, while the number of
recovered individuals kept increasing at a high rate. As for Texas, the
infected population started growing fast from the second half of June, and
then went down after about 30 days.

###  4.2 Parameter Identification Algorithm

In this subsection, we describe how to identify parameters of the SIRD model
by utilizing PINN.

####  4.2.1 PINN Architecture for SIRD Model

The utilized PINN for learning the parameters is composed of an FNN
architecture. This architecture consists of 4 hidden layers and 60 neurons in
each layer. The activation functions of hidden layers and the output layer are
the hyperbolic tangent function, and sigmoid function, respectively.

Figure 4: Physics Informed Neural Network

Figure 4 shows the basic framework, which is FNN, and the loss function of
PINN in our study. The loss function of PINN has two parts, namely $GE_{loss}$
and $OB_{loss}$ in this figure. The $GE_{loss}$ represents the residuals that
are obtained from the SIRD model eq. (1) by subtracting the right side from
the left side:

$\displaystyle\begin{aligned} R_{1}&=\frac{dS}{dt}+\beta(t)S(t)I(t)\\\
R_{2}&=\frac{dI}{dt}-\beta(t)S(t)I(t)+(\gamma(t)+\mu(t))I(t)\\\
R_{3}&=\frac{dR}{dt}-\gamma(t)I(t)\\\
R_{4}&=\frac{dD}{dt}-\mu(t)I(t).\end{aligned}$ (7)

The $OB_{loss}$ is the mean squared error between the outputs of the neural
network and data. We employ Adam algorithm, which is the first-order gradient-
based optimization of stochastic objective functions [14], to update the
parameters of the neural network by minimizing the loss function.

####  4.2.2 Parameters Identification Algorithm

Algorithm 1 describes how to use PINN to identify the daily time-varying
parameters. The input is time $t$, and outputs include the three parameters
and the four compartments of the SIRD model. The weights $w$ and bias $b$ are
initialized randomly.

Data: $S$,$I$,$R$,$D$,$t$

_Initialize weights $w$, bias $b$_;

for  _number of epochs_ do

Each compartment of SIRD is obtained from forwarding propagation of neural
network  $\hat{S},\hat{I},\hat{R},\hat{D},\beta,\gamma,\mu=NN(t)$

Build the loss function components, where $N_{ob}$ is the number of collected
data of each compartment, and $N_{f}$ is the number of collocation points
$OB_{loss}=\frac{1}{N_{ob}}\sum_{i=1}^{N_{ob}}\Big{(}\big{|}\hat{S}^{i}-S^{i}\big{|}+\big{|}\hat{I}^{i}-I^{i}\big{|}+\big{|}\hat{R}^{i}-R^{i}\big{|}+\big{|}\hat{D}^{i}-D^{i}\big{|}\Big{)}\;$
$OB_{loss}$ denotes the difference between the output of the neural network
and observation data
$GE_{loss}=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}\Big{(}\Big{|}\frac{d\hat{S}^{i}}{dt^{i}}+\beta\hat{S}^{i}\hat{I}^{i}\Big{|}+\Big{|}\frac{d\hat{I}^{i}}{dt^{i}}-\beta\hat{S}^{i}\hat{I}^{i}+(\gamma+\mu\hat{I}^{i})\Big{|}\\\
+\Big{|}\frac{d\hat{R}^{i}}{dt^{i}}-\gamma\hat{I}^{i}\Big{|}+\Big{|}\frac{d\hat{D}^{i}}{dt^{i}}-\mu\hat{D}^{i}\Big{|}\Big{)}$
$GE_{loss}$ represents the SIRD model residual loss, then we have the total
loss  $loss=OB_{loss}+GE_{loss}$ Using the Adam optimizer to update the
weights and bias by minimizing the loss function

end for

Algorithm 1 PINN identify daily time-varying parameters

The approach of identifying the weekly-varying parameters is described in
algorithm 2. We divide the utilized data of algorithm 1 into weekly intervals,
and then they are trained to identify corresponding parameters. Cubic spline
interpolation is employed to obtain adequate training set for the neural
network training.

Data: $S$,$I$,$R$,$D$,$t$

_Form the weekly data_ ;

for  _each week_ do

$S,I,R,D,t=CubicSpline(S,I,R,D,t)$ ;

$N=length(t)$ ;

$\triangleright$Each dataset has the same length

_Initialize weights $w$, bias $b$,$\beta$,$\gamma$,$\mu$_;

for  _number of epochs_ do

Each compartment of SIRD is obtained from forwarding propagation of neural
network  $\hat{S},\hat{I},\hat{R},\hat{D}=NN(t)$

Build the loss function components.
$OB_{loss}=\frac{1}{N}\sum_{i=1}^{N}\Big{(}\big{|}\hat{S}^{i}-S^{i}\big{|}^{2}+\big{|}\hat{I}^{i}-I^{i}\big{|}^{2}+\big{|}\hat{R}^{i}-R^{i}\big{|}^{2}+\big{|}\hat{D}^{i}-D^{i}\big{|}^{2}\Big{)}\;$
$OB_{loss}$ denotes the difference between the output of the neural network
and observation data
$GE_{loss}=\frac{1}{N}\sum_{i=1}^{N}\Big{(}\Big{|}\frac{d\hat{S}^{i}}{dt^{i}}+\beta\hat{S}^{i}\hat{I}^{i}\Big{|}^{2}+\Big{|}\frac{d\hat{I}^{i}}{dt^{i}}-\beta\hat{S}^{i}\hat{I}^{i}+(\gamma+\mu\hat{I}^{i})\Big{|}^{2}\\\
+\Big{|}\frac{d\hat{R}^{i}}{dt^{i}}-\gamma\hat{I}^{i}|^{2}+\Big{|}\frac{d\hat{D}^{i}}{dt^{i}}-\mu\hat{D}^{i}\Big{|}^{2}\Big{)}$
$GE_{loss}$ represents the SIRD model residual loss
$loss=OB_{loss}+GE_{loss}$ Using the Adam optimizer to update the weights and
bias by minimizing the loss function

end for

 _save the value of $\beta$ $\gamma$ and $\mu$_

end for

Algorithm 2 PINN identify weekly time-varying parameters

###  4.3 LSTM for Prediction of Infectious Cases

Having the learned parameters by PINN, the LSTM is employed to predict the
future parameters by applying Keras as depicted in Figure 5. At first, we
normalize these parameters and create new datasets with multiple inputs and
one output. Inputs are previous three-time steps parameters, and output, which
is the prediction of these inputs, is the next one-time step parameter. Then,
the data is used to train the LSTM neural network. We utilize three hidden
layers of 80 nodes each to predict the parameters of the next four weeks.

Figure 5: Flow-Chart of LSTM

##  5 Simulation

In this section, we present and discuss our simulation results. The accuracy
of the learned parameters by PINN is validated using the relative error in the
numerical solution of System (1), obtained using the Fourth Order Runge-Kutta
(RK4) method, with respect to the exact solution. Furthermore, the relative
error in the learned solution by PINN with respect to the data is examined. We
start by presenting the results for the learned daily and weekly time-varying
parameters followed by the associated reproduction numbers. Then, predictions
of infectious cases are provided.

###  5.1 Time-Varying Parameters

In general, the daily-varying values of the parameters are more reflective of
the daily data than the single constant average values of these parameters.
However, identifying the daily-varying parameters is traditionally costly and
adds to the complexity of the problem. Fortunately, PINN provides an efficient
approach to overcome this difficulty. On the other hand, our simulations show
that the learned weekly-varying data produce precisely identified parameters
and more stable approximation results.

Algorithms 1 and 2 give the identified parameters and the learned values of
the SIRD model by PINN. We use the RK4 method to solve the SIRD model with
learned parameters. Then, the solutions of the ODE system (1) are compared
with the learned values.

In Figures 6, 7, and 8, we present the learned values by PINN (NN), the
solution of the ODE system (1) (ode), and the reported data (data). Moreover,
the relative errors with respect to the reported data are also computed.

(a)

(b)

Figure 6: Simulation Results of New York

(a)

(b)

Figure 7: Simulation Results of Texas

(a)

(b)

Figure 8: Simulation Results of New Mexico

There are some large relative errors at the early stage of the outbreak of the
disease. This could be due to the small size of the recovered and deceased
populations. But the relative errors of the solution of the SIRD model are in
a reasonable range, which means that we obtained ideal identified parameters.
For the weekly-varying parameters, the period is divided into weeks and the
learned results are presented in Figures 9, 10, and 11.

(a)

(b)

Figure 9: Simulation Results of New York

(a)

(b)

Figure 10: Simulation Results of Texas

(a)

(b)

Figure 11: Simulation Results of New Mexico

The relative errors in Figures 9, 10, and 11 are smaller compared to the
results in Figures 6, 7, and 8. Thus, we could infer that the approximation of
constant parameters is more accurate and stable.

###  5.2 Simulated Reproduction Numbers

By eq. (3), the effective reproduction number can be obtained using the
learned daily and weekly varying parameters. Figure 12 shows the effective
reproduction number in the three states. We could see the daily and weekly
time-varying reproduction numbers have the similar trends, which proves the
feasibility and effectiveness of our methods to some extent. In the result of
New York, we could see the $R_{e}$ is greater than 1 in the whole period,
which represents that the number of infectious cases keeps increasing. As for
the situation in New Mexico, $R_{e}$ decreases to 1 from August, then the
population of infected individuals decreased at that period. By observing the
change of $R_{e}$ in Texas, the value is too large in the beginning because
the parameter identification in that period is not ideal. But after that
period, the $R_{e}$ goes back to normal and oscillates between 0 and 4.
Besides, we could see that the growth rate of infected individuals in New York
is faster than the other two states by Figure 3.

(a) New York

(b) New Mexico

(c) Texas

Figure 12: Effective Reproduction Number

###  5.3 Prediction of Infectious Cases

The LSTM is employed to predict the future values of the parameters. Figure 13
presents the predictions of infected cases for the next four weeks. The
prediction of New York is close to the original trend of infectious cases.
Predictions in New Mexico show general accord with real situations, but it
cannot capture the small fluctuations. As for Texas, there is a sharp increase
and decrease in the infectious cases at the end of September, and our
prediction could not detect this situation.

(a) New York

(b) New Mexico

(c) Texas

Figure 13: Predictions

##  6 Conclusions

We introduced a data-driven deep learning approach based on physics informed
neural network to solve the epidemiological models and identify weekly and
daily time-varying parameters. The PINN was used for parameters identification
and solving the ODE system. LSTM was implemented to predict the weekly time-
varying parameters and then substituted predicted parameters into the SIRD
model so that we could obtain the future trend of COVID-19. The results and
errors have shown that the weekly time-varying parameters provided a good fit
to the real data. The algorithms could be further developed to achieve a more
accurate approximation for complex problems. We intend to explore other
architectures of PINN in future work.

## References

  * [1] S.i. Amari, _Backpropagation and stochastic gradient descent method_ , Neurocomputing 5 (1993), pp. 185–196.
  * [2] J.M. Carcione, J.E. Santos, C. Bagaini, and J. Ba, _A simulation of a COVID-19 epidemic based on a deterministic SEIR model_ , arXiv:2004.03575 (2020).
  * [3] Y.C. Chen, P.E. Lu, C.S. Chang, and T.H. Liu, _A Time-dependent SIR model for COVID-19 with undetectable infected persons_ , IEEE Transactions on Network Science and Engineering 7 (2020), pp. 3279–3294.
  * [4] H. Fan, M. Jiang, L. Xu, H. Zhu, J. Cheng, and J. Jiang, _Comparison of Long Short Term Memory Networks and the Hydrological Model in Runoff Simulation_ , Water 12 (2020), p. 175.
  * [5] L. Ferrari, G. Gerardi, G. Manzi, A. Micheletti, F. Nicolussi, and S. Salini, _Modelling provincial COVID-19 epidemic data in Italy using an adjusted time-dependent SIRD model_ , arXiv:2005.12170 (2020).
  * [6] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_ , MIT Press, 2016, http://www.deeplearningbook.org.
  * [7] D.L. Gray and A.N. Michel, _A training algorithm for binary feedforward neural networks_ , IEEE Transactions on Neural Networks 3 (1992), pp. 176–194.
  * [8] D. Guanghong, L. Chang, G. Jianqiu, W. Ling, C. Ke, and Z. Di, _SARS epidemical forecast research in mathematical model_ , Chinese Science Bulletin 49 (2004), pp. 2332–2338.
  * [9] R.J. Hatchett, C.E. Mecher, and M. Lipsitch, _Public health interventions and epidemic intensity during the 1918 influenza pandemic_ , Proceedings of the National Academy of Sciences 104 (2007), pp. 7582–7587.
  * [10] S. Hochreiter and J. Schmidhuber, _Long short-term memory_ , Neural Computation 9 (1997), pp. 1735–1780.
  * [11] H. Jo, H. Son, S.Y. Jung, and H.J. Hwang, _Analysis of COVID-19 spread in South Korea using the SIR model with time-dependent parameters and deep learning_ , medRxiv (2020).
  * [12] A. Kergassner, C. Burkhardt, D. Lippold, S. Nistler, M. Kergassner, P. Steinmann, D. Budday, and S. Budday, _Meso-scale modeling of COVID-19 spatio-temporal outbreak dynamics in Germany_ , medRxiv (2020).
  * [13] W.O. Kermack and A.G. McKendrick, _A contribution to the mathematical theory of epidemics_ , Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character 115 (1927), pp. 700–721.
  * [14] D.P. Kingma and J. Ba, _Adam: A method for stochastic optimization_ , arXiv:1412.6980v9 (2017).
  * [15] Q. Lin, S. Zhao, D. Gao, Y. Lou, S. Yang, S.S. Musa, M.H. Wang, Y. Cai, W. Wang, L. Yang, _et al._ , _A conceptual model for the outbreak of Coronavirus disease 2019 (COVID-19) in Wuhan, China with individual reaction and governmental action_ , International Journal of Infectious Diseases 93 (2020), pp. 211–216.
  * [16] Z. Liu, P. Magal, O. Seydi, and G. Webb, _Predicting the cumulative number of cases for the COVID-19 epidemic in China from early data_ , Mathematical Biosciences and Engineering 17 (2020), pp. 3040–3051.
  * [17] L. Lu, X. Meng, Z. Mao, and G.E. Karniadakis, _DeepXDE: A deep learning library for solving differential equations_ , arXiv:1907.04502v2 (2020).
  * [18] L. Magri and N.A.K. Doan, _First-principles machine learning modelling of COVID-19_ , arXiv:2004.09478 (2020).
  * [19] M. Rafiq, J. Macías-Díaz, A. Raza, and N. Ahmed, _Design of a nonlinear model for the propagation of COVID-19 and its efficient nonstandard computational implementation_ , Applied Mathematical Modelling 89 (2021), pp. 1835–1846.
  * [20] M. Raissi, P. Perdikaris, and G.E. Karniadakis, _Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations_ , Journal of Computational Physics 378 (2019), pp. 686–707.
  * [21] M. Raissi, N. Ramezani, and P. Seshaiyer, _On parameter estimation approaches for predicting disease transmission through optimization, deep learning and statistical inference methods_ , Letters in Biomathematics 6 (2019), pp. 1–26.
  * [22] B. Ridenhour, J.M. Kowalik, and D.K. Shay, _Unraveling $R_{0}$: Considerations for Public Health Applications_, American Journal of Public Health 108 (2018), pp. S445–S454.
  * [23] F. Shahid, A. Zameer, and M. Muneeb, _Predictions for COVID-19 with deep learning models of LSTM, GRU and Bi-LSTM_ , Chaos, Solitons & Fractals 140 (2020), p. 110212.
  * [24] Y. Shi, _Stochastic dynamic model of SARS spreading_ , Chinese Science Bulletin 48 (2003), pp. 1287–1292.
  * [25] P. van den Driessche, _Reproduction numbers of infectious disease models_ , Infectious Disease Modelling 2 (2017), pp. 288–303.
  * [26] H. Wang, Z. Wang, Y. Dong, R. Chang, C. Xu, X. Yu, S. Zhang, L. Tsamlag, M. Shang, J. Huang, _et al._ , _Phase-adjusted estimation of the number of coronavirus disease 2019 cases in Wuhan, China_ , Cell Discovery 6 (2020), pp. 1–8.
  * [27] J. Wang, A.J. McMichael, B. Meng, N.G. Becker, W. Han, K. Glass, J. Wu, X. Liu, J. Liu, X. Li, _et al._ , _Spatial dynamics of an epidemic of severe acute respiratory syndrome in an urban area_ , Bulletin of the World Health Organization 84 (2006), pp. 965–968.
  * [28] World Health Organizaiton, _Rolling updates on coronavirus disease (COVID-19)_ (2020). Available at https://www.who.int/emergencies/diseases/novel-coronavirus-2019/events-as-they-happen, Updated 31 July 2020.
  * [29] Worldometers, _COVID-19 CORONAVIRUS PANDEMIC_ (2020). Available at https://www.worldometers.info/coronavirus/?#countries.
  * [30] Z. Yang, Z. Zeng, K. Wang, S.S. Wong, W. Liang, M. Zanin, P. Liu, X. Cao, Z. Gao, Z. Mai, _et al._ , _Modified SEIR and AI prediction of the epidemics trend of COVID-19 in China under public health interventions_ , Journal of Thoracic Disease 12 (2020), p. 165.
  * [31] A. Zeroual, F. Harrou, A. Dairi, and Y. Sun, _Deep learning methods for forecasting covid-19 time-series data: A comparative study_ , Chaos, Solitons & Fractals 140 (2020), p. 110121.