# Mapping and Validating a Point Neuron Model on Intel’s Neuromorphic Hardware
Loihi

Srijanie Dey  Department of Mathematics, Washington State University,
Vancouver, WA, 98686  Alexander Dimitrov  Department of Mathematics,
Washington State University, Vancouver, WA, 98686

###### Abstract

Abstract Neuromorphic hardware is based on emulating the natural biological
structure of the brain. Since its computational model is similar to standard
neural models, it could serve as a computational acceleration for research
projects in the field of neuroscience and artificial intelligence, including
biomedical applications. However, in order to exploit this new generation of
computer chips, rigorous simulation and consequent validation of brain-based
experimental data is imperative. In this work, we investigate the potential of
Intel’s fifth generation neuromorphic chip - ‘Loihi’, which is based on the
novel idea of Spiking Neural Networks (SNNs) emulating the neurons in the
brain. The work is implemented in context of simulating the Leaky Integrate
and Fire (LIF) models based on the mouse primary visual cortex matched to a
rich data set of anatomical, physiological and behavioral constraints.
Simulations on the classical hardware serve as the validation platform for the
neuromorphic implementation. We find that Loihi replicates classical
simulations very efficiently and scales notably well in terms of both time and
energy performance as the networks get larger.
Keywords - Neuromorphic Computing, LIF, Neural Simulations, Validation,
Performance Analysis

##  1 Introduction

The human brain is a rich complex organ made up of numerous neurons and
synapses. Replicating the brain structure and functionality in classical
hardware is an ongoing challenge given the complexity of the brain and
limitations of hardware. The advent of supercomputers now allows for complex
neural model, but at a huge cost of both software complexity and energy
consumption.

A recent intense focus on brain studies, with the BRAIN Insel  _et al._ (2013)
initiative at the US, the Human Brain Project (HBP) Markram  _et al._ (2011)
in Europe, and philanthropic endeavors like Janelia Research Campus Winnubst
_et al._ (2019), and the Allen Institute for Brain Science (AIBS) Lein  _et
al._ (2007), has produced a wealth of new data and knowledge, from records of
neuronal and network dynamics, to fine-grained data on network micro- and
nano-structure, bringing in the era of big neural data. At the same time,
advances in electronics and the search for post-von Neumann computational
paradigms has led to the creation of neuromorphic systems like Intel’s Loihi
Davies  _et al._ (2018), IBM’s TrueNorth DeBole  _et al._ (2019); Löhr  _et
al._ (2020); Akopyan  _et al._ (2015), University of Manchester’s SpiNNaker
Khan  _et al._ (2008) and HBP’s BrainScaleS Grübl  _et al._ (2020).

Neuromorphic chips, as the name suggests - ‘like the brain’ - can mimic the
brain’s function in a truer sense as their design is analogous to the brain
Roy  _et al._ (2019); Thakur  _et al._ (2018). Inspired by its architecture,
we work on developing a principled approach towards obtaining simulations of
biologically realistic neural network models on a novel neuromorphic
commercial hardware platform.

Computers today are limited in this respect because of the way they have been
built historically and the way they process data leading towards more energy
and resource consumption in order to maintain versatility Nawrocki  _et al._
(2016); Ou  _et al._ (2020). Neuromorphic chips on the other hand claim to be
faster and more efficient for a set of specialized tasks Sharp and Furber
(2013); Bhuiyan  _et al._ (2010). In this study, we try to validate this
assertion based on neural models derived from the primary visual cortex (VISp)
of the mouse brain, as seen in recent work done on SpiNNaker Rhodes  _et al._
(2019); Knight and Furber (2016). Intel’s neuromorphic chip Loihi serves as
our neuromorphic platform. Results obtained in Loihi are validated against
classical simulations Rossant  _et al._ (2010); Nandi  _et al._ (2020); Wang
_et al._ (2020) given by AIBS’s software package the Brain Modeling Toolkit
(BMTK) Dai  _et al._ (2020). Overall, our implementation indicates that Loihi
is highly efficient in terms of time and energy in context of large brain
network simulations and thus shapes our central motivation for this work (see
Figure 1 and Table 1). The majority of this manuscript focuses on the trade-
offs necessitated by these improvements, that is, how precise are the Loihi
simulations when validated against BMTK simulations?

Figure 1: As the network size increases, Loihi performs consistently both in terms of time and energy. Simulation of 500$\hskip 2.27626ptms$ of dynamics for 10,000 neurons, results in a Loihi run-time of up to 109 seconds with maximum energy consumption of 5.7 Watts. (See Table 1 for the explicit values.)  Network Size | Loihi Time(s) | BMTK Time(s)
---|---|---
20 | 3.16 | 16.08
100 | 3.94 | 80.59
500 | 8.48 | 459.76
1000 | 12.76 | 570.34
5000 | 51.69 | 830.17
10000 | 109.02 | 1460.40
Table 1: Time Usage in Loihi & BMTK

As a first step in this approach, we focus on a class of neural network
building blocks: point neuronal models as used in large AIBS simulations. The
Generalized Leaky Integrate and Fire Models (GLIFs) Teeter  _et al._ (2018)
have been found to be appropriate for reproducing cellular data under
standardized physiological conditions. The data used for this study is made
available by the AIBS.

The paper is organized as follows. In Section 2, we describe in detail the
features of Loihi and the differences between the neuromorphic and classical
hardware that form the basis for this study. Section 3 explains the
implementation of the continuous LIF equation in BMTK on classical
computational architecture, versus the discrete Loihi setting. Also, we list
the validation methods and the cost function that is used to draw comparisons
between the implementations. In section 4, we list out and explain the various
results leading to a qualitative and quantitative assessment between the two
platforms based in part on methods from Gutzen  _et al._ (2018).

##  2 COMPARISON BETWEEN CLASSICAL AND NEUROMORPHIC PLATFORMS

At present, various simulators are available for implementing spiking neural
networks Brette  _et al._ (2008). In this section, we lay out the details of
the mathematical model and the platforms we use for our work. For the
classical simulation, we use the Brain Modeling Toolkit (BMTK) Dai  _et al._
(2020) developed by the AIBS. Being open source, these resources enable us to
experiment with a varied range of data and thus support our extensive
validation of neuronal models in Loihi. Intel’s fifth-generation chip Loihi
provides us with the tools to implement and test out the various neuromorphic
features. The output provided by Loihi simulations is then compared to the
output of classical simulations implemented in BMTK.

###  2.1 The Brain Modeling Toolkit (BMTK)

The BMTK is a python-based software package for creating and simulating large-
scale neural networks. It supports models of different resolutions, namely,
Biophysical Models, Point Models, Filter Models and Population Models along
with the use of the rich data sets of the Allen Cells Database Lein  _et al._
(2007). It leverages the modeling file format SONATA which includes details on
cell, connectivity and activity properties of a network along with being
compatible with the neurophysiology data format Neurodata Without Borders
(NWB), thus allowing easy access to a vast repertoire of experimental data.

In this study, we work with the Point Neuron Models with simulations supported
by the BMTK module PointNet via NEST 2.11+ Kunkel  _et al._ (2017). For
analysis and visualization we use the HDF5 output format, underlying both
SONATA and NWB’s spike and time series storage.

###  2.2 Loihi

Hardware inspired by the structure and functionality of the brain, envisioned
to provide advantages such as low power consumption, high fault tolerance and
massive parallelism for the next generation of computers, is called
neuromorphic hardware. Towards the end of 2017, Intel Corporation unveiled its
experimental neuromorphic chip called Loihi. We provide a summary of the
platform here, since, unlike NEST and BMTK, Loihi is not so well-known in the
computational neuroscience community.

As of its 2020 rendition, Loihi is a 60-mm2 chip that implements 131,072
leaky-integrate-and-fire neurons. According to Davies  _et al._ (2018), it
uses an asynchronous spiking neural network (SNN), comprising of 128
neuromorphic cores, each with 1024 neural computational units; 3 x86 cores;
along with several off-chip communication interfaces that provide connectivity
to other chips. As Loihi advances the modeling of SNNs in silicon, it
comprises of a large number of features necessary for their implementation
viz., hierarchical connectivity, dendritic compartments, synaptic delays and
synaptic learning rules. Each neuron is represented as a compartment in the
Loihi architecture. The top-level microarchitecture of a Loihi neural
computational unit resembles an actual biological neuron model comprising of
all the functional units. The SYNAPSE unit processes all the incoming spikes
from the previous compartment/neuron and captures the synaptic weight from the
memory. The DENDRITE unit updates the different state variables. The AXON unit
generates the spike message to be carried ahead by the fan out cores. The
LEARNING unit updates the synaptic weights based on a learning rule.

The aim of this study is to establish the groundwork required to execute an
ambitious plan of simulating about $\sim$250,000 neurons with $\sim$500M
synapses, which encapsulates much of the experimentally observed dynamics in
the mouse visual cortex available to the AIBS, thus providing a close
functional replica of the mouse visual cortex. Loihi’s specialized hardware
features holds promise for a real-time, low-powered version of this
implementation.

###  2.3 Leaky Integrate and Fire Model (LIF)

A typical neuron consists of a soma, dendrites and a single axon. Neurons send
signals along an axon to a dendrite through junctions called synapses. The
classical Leaky Integrate and Fire (LIF) equation Gerstner and Kistler (2002)
is a point neuron model which reduces much of the neural geometry and dynamics
in order to achieve computational efficiency. It is one of the simplest and
rather efficient widely used representations of the dynamics of the neuron,
and is given as:

$\displaystyle V^{\prime}(t)$
$\displaystyle=\frac{1}{C}\left[I_{e}(t)-\frac{1}{R}(V(t)-E_{L})\right]$ (1)
$\displaystyle V(t)$ $\displaystyle\leftarrow{V_{r}}\text{, \hskip 2.84544pt
if \hskip 2.84544pt}V(t)>\Theta$ (2)

where,

$\displaystyle V(t)$ $\displaystyle=\text{membrane potential (state)}$
$\displaystyle C$ $\displaystyle=\text{membrane capacitance (parameter)}$
$\displaystyle R$ $\displaystyle=\text{membrane resistance (parameter)}$
$\displaystyle E_{L}$ $\displaystyle=\text{resting potential (parameter)}$
$\displaystyle I_{e}$ $\displaystyle=\text{trans-membrane current (control and
state)}$ $\displaystyle V_{r}$ $\displaystyle=\text{reset membrane potential}$
$\displaystyle\Theta$ $\displaystyle=\text{firing threshold}$

Here, ${}^{\prime}=d/dt$, $t$ is time in $ms$, the membrane potential $V(t)$
of the neuron is in mV. A LIF neuron fires when $V(t)>\Theta$, i.e. the
membrane potential exceeds the firing threshold $\Theta$ and subsequently the
membrane potential is set to a reset value $V_{r}$.

The classical LIF model (point generalized LIF) has been shown to match
closely the dynamics of real neurons under a variety of conditions Teeter  _et
al._ (2018), as listed in the Allen Cell Types Database Lein  _et al._ (2007).
In addition, this model matches the LIF abstraction in Loihi to some extent.
Thus, we work with this model throughout this study to establish the basis for
comparison for the two platforms, validate the neuromorphic implementation
against the ground truth of a standard implementation, and provide evidence
that our neuromorphic platform performs more efficiently.

###  2.4 Loihi LIF Model

In an SNN, spiking neurons form the primary processing elements. The
individual neurons are connected through junctions called synapses and
interact with each other through single-bit events called spikes. Each spike
train can be represented as a list of event times, e.g. as a sum of Dirac
delta functions $\sigma(t)=\sum_{i}\delta(t-t_{i})$ where $t_{i}$ is the time
of the $i$-th spike.

Since Loihi encapsulates the working of an SNN, one of the computational
models it implements is a variation of the LIF model based on two internal
state variables : the synaptic current and the membrane potential Davies  _et
al._ (2018).

$\displaystyle u(t)$ $\displaystyle=\sum_{j}w_{j}(\alpha_{j}*\sigma_{j})(t)+b$
(3) $\displaystyle\dot{v}(t)$ $\displaystyle=-\frac{1}{\tau_{v}}v(t)+u(t)$ (4)
$\displaystyle v(t)$ $\displaystyle\leftarrow{0}\text{,\hskip 2.84544pt if
\hskip 2.84544pt }v(t)>\theta$ (5)

where,

$\displaystyle v(t)$ $\displaystyle=\text{membrane potential}$ $\displaystyle
u(t)$ $\displaystyle=\text{synaptic current}$ $\displaystyle w$
$\displaystyle=\text{synaptic weight}$ $\displaystyle\alpha$
$\displaystyle=\text{synaptic response function}$ $\displaystyle b$
$\displaystyle=\text{constant bias current}$ $\displaystyle\theta$
$\displaystyle=\text{firing threshold}$

A neuron sends out a spike when its membrane potential exceeds its firing
threshold $\theta$, i.e., $v(t)>\theta$. After a spike occurs, $v(t)$ is reset
to $0$. As in the classical LIF model, here $\dot{}=d/dt$. However, time and
membrane potential values here are in arbitrary units.

Loihi follows a fixed-size discrete time-step model, similar to an explicit
Euler integration scheme, where the time steps relate to the algorithmic time
of the computation. This algorithmic time may differ from the hardware
execution time. Moreover, to increase the efficiency of the chip, specific
bit-size constraints are imposed on the state variables. We discuss the ones
relevant for the LIF model implementation in the following section.

##  3 METHODS

###  3.1 Model setup and integration

The classical LIF model as represented in equation (1) and (2) can be
rewritten as :

$\displaystyle V^{\prime}(t)$
$\displaystyle=-\frac{1}{\tau}V(t)+\frac{1}{C}\left[I_{e}(t)+\frac{1}{R}E_{L}\right]$
(6)

where $\tau=RC$ is membrane time constant of the neuron.

For a non-homogeneous linear differential equation,

$\displaystyle\frac{df}{dt}=af+g$ (7)

the solution is given by the ‘variation of constants’ method as :

$\displaystyle f(t)=e^{ct}\int_{0}^{t}g(s)e^{-cs}ds$

Comparing equation (6) to equation (7), we have,

$\displaystyle f$ $\displaystyle=V(t)$ $\displaystyle g$
$\displaystyle=\frac{1}{C}(I_{e})+\frac{1}{C}(E_{L})$

Here, the postsynaptic current $I_{e}$ is in the form of an exponent function.
However, calculating the above integral at every step i.e. at all grid points
$t_{i}\leq t$ proves to be quite expensive.

BMTK uses NEST as backend to implement the above membrane potential dynamics.
To avoid the expensive computations, NEST chooses to use the linear exact
integration method Rotter and Diesmann (1999), given below as follows :

Equation (6) is rewritten as a multidimensional homogeneous differential
equation:

$\displaystyle\frac{d}{dt}y=Ay$ (8)

where,

$A=\begin{pmatrix}a_{n}&a_{n-1}&\cdots&\cdots&a_{1}&0\\\ 1&0&\cdots&0&0&0\\\
0&\ddots&\ddots&\vdots&\vdots&\vdots\\\ \vdots&\ddots&\ddots&0&0&0\\\
0&0&\ddots&1&0&0\\\ 0&0&\cdots&0&\frac{1}{C}&\frac{1}{-\tau}\end{pmatrix}$

The solution is given by :

$\displaystyle y(t)$ $\displaystyle=e^{At}y_{0}$ (9) $\displaystyle y_{t+h}$
$\displaystyle=y(t+h)=e^{A(t+h)}y(0)=e^{Ah}\cdot y_{t}$ (10)

for a fixed time-step $h$. It saves exorbitant computations since each step
involves multiplication only.

####  3.1.1 Mapping between BMTK and Loihi models

In this section, we illustrate the primary step of implementing the BMTK-NEST
LIF integration into the fixed-size discrete time-step Loihi dynamics.
Following the linear exact numerics in the NEST implementation, we implement
our model in the Loihi discrete setting using the forward Euler method as
discussed below.

First, we rewrite equation (1), the standard LIF model to resemble the Loihi
form as given in equation (4). Since Loihi parameters are unit-less, we
introduce a re-scaling parameter $V_{s}$, which converts standard physical
units used in BMTK to Loihi units.

As we compare equation (2) and equation (5), it can be seen that for BMTK the
membrane potential reset value is set to $V_{r}$ whereas it is set to zero for
Loihi.

Thus, the forward transformation from BMTK to Loihi looks as follows :

$\displaystyle v=(V-V_{r})/V_{s}$ (11)

which produces an inverse transformation, to arrive back at the BMTK values,
given by :

$\displaystyle V=v\cdot V_{s}+V_{r}$ (12)

Substituting the expression in (12) in (1) and isolating $v$, we get :

$\displaystyle V^{\prime}(t)$
$\displaystyle=\frac{1}{C}\left[I_{e}(t)-\frac{1}{R}(V(t)-E_{L})\right]$
$\displaystyle|V=vV_{s}+V_{r}\implies$ (13) $\displaystyle v^{\prime}(t)V_{s}$
$\displaystyle=-\frac{1}{RC}(vV_{s}+V_{r}-E_{L})+\frac{1}{C}I_{e}(t)$
$\displaystyle|/V_{s}\implies$ (14) $\displaystyle v^{\prime}(t)$
$\displaystyle=-\frac{1}{\tau_{v}}v(t)+\frac{1}{\tau_{v}}\frac{E_{L}-V_{r}}{V_{s}}+\frac{1}{C}\frac{I_{e}}{V_{s}}$
(15) $\displaystyle=-\frac{1}{\tau_{v}}v(t)+u(t)$ (16)

with

$\displaystyle v(t)$ $\displaystyle=\frac{V(t)-V_{r}}{V_{s}},$ (17)
$\displaystyle u(t)$
$\displaystyle=\frac{1}{CV_{s}}I_{e}(t)+\frac{1}{\tau_{v}}\frac{E_{L}-V_{r}}{V_{s}},$
(18) $\displaystyle\tau_{v}$ $\displaystyle=RC,$ (19) $\displaystyle\theta$
$\displaystyle=\frac{\Theta-V_{r}}{V_{s}}$ (20)

Here, we reintroduce the LIF threshold $\Theta$ and the corresponding Loihi
threshold $\theta$ in equation (20).

Loihi implements the continuous LIF as a discrete finite state machine model
Xin Jin  _et al._ (2008); Mikaitis  _et al._ (2018) implemented in silicon.
The actual computation is similar to a forward Euler scheme with some
peculiarities reflecting engineering design trade-offs. Specifically, the
$v(t)$ state evolves on-chip according to the update rule,

$\displaystyle v(t+1)=v(t)\left[1-\frac{\delta_{v}}{2^{12}}\right]+b+u(t)$
(21)

where $\delta_{v}$ is the membrane decay constant.

Using the forward Euler method :

$\displaystyle y_{n+1}=y_{n}+f(t_{n},y_{n}).dt$

where $y_{n+1}=y(t_{n+1})$ and $t_{n+1}=t_{n}+dt$ for a fixed time-step $dt$,
we transform the classical LIF model into a form followed in equation (21).
Thus, transforming the LIF model into the discrete form and grouping terms to
match the Loihi integration (9) yields the following :

$\displaystyle\frac{v(t+dt)-v(t)}{dt}$
$\displaystyle=-\frac{1}{\tau_{v}}v(t)+u(t)$ (22) $\displaystyle\implies
v(t+dt)$ $\displaystyle=v(t)(1-\frac{dt}{\tau_{v}})+u(t)dt$ (23)

where $dt$ is the fixed time-step.

In order to equate the Loihi computation (21) with the Euler scheme (23), we
use $dt$ with units _ms/Loihi timestep_ i.e. 1 BMTK millisecond per Loihi
timestep that leads to a Loihi voltage decay parameter $\delta_{v}$ such that,

$\displaystyle(2^{12}-\delta_{v})~{}2^{-12}$
$\displaystyle=(1-\frac{dt}{\tau_{v}})$ (24) $\displaystyle\implies\delta_{v}$
$\displaystyle=\frac{dt}{\tau_{v}}~{}2^{12}=\frac{dt}{RC}~{}2^{12}$ (25)

###  3.2 Validation Methods

  * Data
The data used here is provided by the Allen Mouse Brain Atlas Lein  _et al._
(2007), which is a survey of single cells from the mouse brain, obtained via
intracellular electrophysiological recordings done through a highly
standardized process. We focus on neurons with different morphologies with
available GLIF parameters. The data used can be accessed in the Allen Cell
Types Database. Our LIF model is implemented and simulated on BMTK based on
this data, and these simulations form the ground truth for validating the
Loihi implementations.

  * Methods and Cost Function
A spike-to-spike equivalence between the simulated spiking activity and the
experimental data only allows for a qualitative validation. Since BMTK and
Loihi run on two different computing environments, a statistical comparison
becomes imperative to quantify the level of similarity between the two
implementations.

We compare the simulation dynamics for both implementations based on the
following -

    1. Raster Plot: We evaluate the membrane potential response at each time-step. The X-axis represents the membrane potential and the Y-axis represents the time-step. Raster plot helps to highlight the difference in the spike dynamics at each step and thus tune the parameters as needed.

    2. Scatter Plot: For examining association between the two implementations, we use color-coded scatter plots identifying the correlation relationships. We add a trend line to illustrate the strength of the relationship and pin down the outliers to improve the simulation results. Since we anticipate an almost perfect linear relationship, we quantify the match with its Pearson correlation coefficient,

$\displaystyle
r=\frac{\sum_{i=1}^{n}\left({y}_{L}^{i}-\bar{y_{L}}\right)\left({y}_{B}^{i}-\bar{y_{B}}\right)}{\sqrt{\sum_{i=1}^{n}\left({y}_{L}^{i}-\bar{y_{L}}\right)^{2}\sum_{i=1}^{n}\left({y}_{B}^{i}-\bar{y_{B}}\right)^{2}}}$
(26)

where,

$\displaystyle i$ $\displaystyle=\text{index of data point}$
$\displaystyle{y}_{L}$ $\displaystyle=\text{transformed Loihi values}$
$\displaystyle y_{B}$ $\displaystyle=\text{original BMTK values}$
$\displaystyle n$ $\displaystyle=\text{number of data points}$

    3. Distribution Function: Another approach we take is to compare the distributions of attained state values in the two cases. We use density plot as a representation of those distributions, thus allowing us to compare the two implementations in terms of concentration and spread of the values and provide a basis for comparing the collective dynamics of the implementations.

    4. Cost Function: To quantify the error between the BMTK and Loihi membrane potential values, we base our cost function on Root Mean Square Error (RMSE) values as follows :

$\displaystyle
RMSE=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left({y}_{L}^{i}-y_{B}^{i}\right)^{2}}$
(27)

##  4 RESULTS

In order to be able to simulate a network of over 250,000 neurons with a
connectivity of over 500M synapses in the neuromorphic hardware, we begin by
ensuring an high quality replication of individual neural and smaller network
models. The replication performance here is evaluated based on membrane
potential and current responses, the two state variables. We conjecture that
securing a good replica for smaller models will ensure that parameters can be
calibrated correctly and thus can be carried forward for the bigger networks
needed in biological context Gutzen  _et al._ (2018); Trensch  _et al._
(2018); Herz  _et al._ (2006). We begin our work on a single-neuron network
sub threshold dynamics driven by both bias current and external spikes to
ensure Loihi is able to handle both stimuli efficiently. Our test suite
consists of LIF models based on 50 different parameter sets. We perform
rigorous quantitative analysis of our results based on various statistical
measures to demonstrate that we have replication of high quality. It is
important to restate here that we test our results based on neurons with
different morphologies and biophysics, which attribute to the different
parameter sets.

###  4.1 Simulations of a Single Neuron

We begin by simulating a single-neuron network in BMTK. The simulation is run
for 500$\hskip 1.42271ptms$. The classical parameters are translated to Loihi
values and the corresponding LIF model is implemented as an one-neuron SNN
executed for 500 time steps in Loihi. The simulations are driven either by
bias current or external spikes. In the Loihi network, neurons are denoted by
compartments. The compartment dynamics are guided by the parameters - bias
current mantissa, membrane potential threshold, membrane potential decay and
current decay. It is worth iterating here that the membrane potential values
in Loihi are unit-less as opposed to the BMTK values which are assigned units
of millivolts (mV). We test the efficacy of our replication, both
qualitatively and quantitatively, for all 50 parameters sets and find that the
results are consistent with the ones described below.

In Figure 2, we illustrate the implementations achieved through bias current
and external spikes on two different parameter sets.

Figure 2: Membrane potential response for single-neuron network based on two
different neuron parameters. The first simulation is driven by bias current
and the second one by external spikes.

For a qualitative comparison, it can be seen from Figure 2 that Loihi
implementations simulate BMTK results very closely. We have exactness in terms
of spike frequency, spike amplitude, and response values. Since Loihi membrane
potential values are unit-less, we map them back to the BMTK values (mV)
before performing the comparison. The inverse mapping from Loihi to BMTK is
performed based on equation (12), i.e.,

$V=v\cdot V_{s}+V_{r}$

We perform a quantitative assessment of the replication using several
statistical measures. As seen from Table 2, the values on the two platforms
are highly correlated with a relatively small RMSE.

Stimulus | Correlation | RMSE
---|---|---
Bias Current | 0.999992 |  $1.1374\times 10^{-4}$ mV/$ms$
External Spikes | 0.999942 |  $4.208\times 10^{-5}$ mV/$ms$
Table 2: Correlation and RMSE between BMTK and Loihi membrane potential values

Figure 3 illustrates the comparison of Loihi implementations against the BMTK
implementations for the two different stimuli using various statistical
measures - (a) Distribution function approximating the membrane potential
dynamics (b) Raster plot of the spiking network activity (c) Scatter Plot
highlighting the positive coefficient between the two implementations.

(a) Density Plot

(b) Raster Plot

(c) Scatter Plot

Figure 3: Validation Plots for simulations based on two different stimuli. The
first and the second row show the results for Bias Current and External Spike
stimulus respectively.

###  4.2 Simulation Using Varied Precision

As already stated, Loihi follows a fixed-size discrete time-step model along
with bit-size constraints for the different parameters. Thus, we investigate
how changing the precision of the time scale and the neuron state values
affects the accuracy of the simulations. We investigate this property for the
two state variables - membrane potential and current.

Figure 4 below illustrates the membrane potential and current responses of a
single neuron model in BMTK which form the basis of our comparison for the
results below.

Figure 4: Membrane potential response and current response for a single neuron
model in BMTK

####  4.2.1 Simulation Using Varied Temporal Precision

For Loihi’s fixed simulation time-step, we assign different time units to each
step and test the corresponding simulation precision. This is achieved through
the ‘$dt$’ available in our equations while transforming the classical LIF
model to Loihi neural model. It enables us to experiment with several time
units Hopkins and Furber (2015). Following equation (24), the change of a
time-step while working with the Loihi neural model necessitates a
corresponding variation of the time constant ‘$\tau$’ to yield the desired
results.

We check the results for $dt=10.0,1.0\text{ and }0.1(ms/timestep).$ As
mentioned earlier, we run the simulation for 500$\hskip 2.27626ptms$, thus the
corresponding number of time steps in Loihi for $dt=10$ and $dt=0.1$ becomes
50 and 5000 respectively, and for $dt=1.0$ it remains the same. Figure 5
illustrates the related Loihi simulation for membrane potential and current.

Figure 5: Comparison of membrane potential and current plots with different
temporal precisions in Loihi - (a) $dt=0.1$ (b) $dt=1.0$ (c) $dt=10.0$. For
$dt=10.0$, number of time-steps are 50 and for $dt=0.1$, number of time-steps
are 5000.

Error comparison for temporal precision

We compare the simulations in Loihi with different temporal precisions against
the simulations in BMTK. We calculate the RMSE to be able to deduce the
result. As can be seen from Figure 6, the error is lowest when 1$\hskip
2.27626ptms$ of simulation time in BMTK equates to 1 time-step in Loihi for
membrane potential and current. Thus, for the LIF model simulations, Loihi
hardware time is close to the actual timescale of the simulations.

Figure 6: Error comparison for different temporal precisions. In both panels,
the RMSE for the corresponding state is plotted against the $log$ of the
temporal precision $dt$.

####  4.2.2 Simulation Using Varied Voltage Precision

We repeat the precision study by changing the voltage precision values using
the re-scaling parameter $V_{s}$. To check different precision results, we try
$1$K$/$mV, $10$K$/$mV  and  $100$K$/$mV (state level/mV) by using
$V_{s}=1.0\times 10^{-3},1.0\times 10^{-4}$ and $1.0\times 10^{-5}$
respectively. Figure 7 illustrates the neuron state simulations based on
different voltage precisions.

Figure 7: Comparison of membrane potential and current plots with different
voltage precisions -
(a) $V_{s}=1.0\times 10^{-3}$ (b) $V_{s}=1.0\times 10^{-4}$ (c)
$V_{s}=1.0\times 10^{-5}.$

Error comparison for voltage precision

As can be seen from Figure 8, for membrane potential - error decreases
significantly as the precision increases from $V_{s}=1.0\times 10^{-3}$ to
$V_{s}=1.0\times 10^{-4}$. However, the error is extremely small for the
current simulation and remains the same for $V_{s}=1.0\times 10^{-4}$ and
$V_{s}=1.0\times{10}^{-5}$.

Figure 8: Error comparison for different membrane potential precisions. In
both panels, the RMSE for the corresponding state is plotted against the
$-log$ of the voltage scale $V_{s}$.

In conclusion, the effect of precision on both scales depends on the model
parameters and the information needing to be preserved. However, there are
important performance differences. Increased voltage precision is essentially
free, as it does not tax the hardware resources any further, and the sole risk
is from computation overflow in cases of the Loihi voltage state nearing the
precision of the voltage register. Increased time precision on the other hand
has two important drawbacks: it increases simulation time (proportionately to
increased precision), and it decreases the range of voltage decay timescale
that can be represented (again, proportionately to increased precision). Thus,
the choice of simulation time step and corresponding precision should be
weighed against these trade offs.

###  4.3 Simulation of Larger Networks

After establishing and verifying the calibrated Loihi parameters, we extend
our network to larger number of neurons. Currently, we have successfully
implemented a network of $\sim$10,000 neurons. It is done by increasing the
number of cores that are involved in the processing of the network, as each
core supports only 1024 neurons.

(a) BMTK network of 20 neurons

(b) Loihi network of 20 neurons

Figure 9: Loihi replicates large network response of BMTK

Figure 9 illustrates an equivalent simulation for a network of 20 neurons
between BMTK and Loihi indicating that Loihi is capable of emulating larger
BMTK networks. (Figure 9 (b) shows the raw membrane potential values for
Loihi, before the inverse mapping in equation (12) is applied). Here, we have
a correlation of 0.99985 with an RMSE of 0.24 $\times 10^{-3}$ mV/ $ms$. This
also validates the fact that the calibration of parameters for smaller
networks done earlier is valid. This lays the foundation for building bigger
and more complex networks on Loihi.

###  4.4 Comparison of Performance

The performance of Loihi far exceeds that of BMTK. For smaller networks, the
time taken to build, compile and execute the networks is approximately similar
in both the platforms. However, as the networks get larger, BMTK takes much
longer than Loihi to complete the simulations.

(a) Simulation time vs. logarithm of network size

(b) Logarithm of simulation time vs. logarithm of network size

Figure 10: Performance comparison between BMTK and Loihi for network sizes
ranging from 1 to 10000, expressed in $log$ scale.

Figure 10 compares the time taken by BMTK and Loihi to complete the
simulations. We measure the system and user CPU time of each run, excluding
the time elapsed in queue or sleep.

The irregularities in the BMTK curve can be attributed to the initial loading
of the data through ‘h5py’ files which have an initial large load time but
scale well as the network gets larger. Moreover, NEST initializes neurons in
batches and memory pre-allocation follows the $log_{2}$ scale. Care was taken
to run the simulations on a single core on the same node to ensure no
variability in terms of system performance.

Evidently, Loihi could prove to be a huge advantage as we aim to undertake the
mammoth task of emulating the mouse primary visual cortex.

##  5 Conclusion

Inspired by the brain, neuromorphic computing holds great potential in
tackling tasks with extremely low power and high efficiency. Many large-scale
efforts including the TrueNorth, SpiNNaker and BrainScaleS have been
demonstrated as a tool for neural simulations, each replete with its own
strengths and constraints. Fabricated on Intel’s 14nm technology node, Loihi
is a forward-looking and continuously evolving state-of-the-art architecture
for modeling spiking neural networks in silicon. As opposed to its
predecessors, Loihi encompasses a wide range of novel features such as
hierarchical connectivity, dendritic compartments, synaptic delays and
programming synaptic learning rule. These features, together with solid SDK
support by Intel, and a growing research community, make Loihi an effective
platform to explore a wealth of neuromorphic features in more detail than
before.

In this work, we have demonstrated that Loihi is capable of replicating the
continuous dynamics of point neuronal models with high degree of precision and
does so with much greater efficiency in terms of time and energy. The work
comes with its challenges as simulations built on the conventional chips
cannot be trivially mapped to the neuromorphic platform as its architecture
differs remarkably from the conventional hardware. Classical simulations from
the Brain Modeling Toolkit (BMTK) developed by the Allen Institute of Brain
Science (AIBS) serves as the foundation of our neuromorphic validation.

For comparison between the conventional and the neuromorphic platforms, we use
both qualitative and quantitative measures. It can be seen that Loihi
replicates BMTK very closely in terms of both membrane potential and current,
the two state variables on which the Loihi LIF model evolves. We use different
validation methods and quantitative measures to assess the equivalence and
identify sets of parameters which maximize precision while retaining high
performance levels. Furthermore, simulation results indicate Loihi is highly
efficient in terms of speed and scalability as compared to BMTK. As indicated
by the Figure 10, Loihi is almost a thousand times faster than BMTK for
smaller networks and scales markedly well with network size.

In conclusion, this work demonstrates that classical simulations based on
Generalized Leaky Integrate-and-Fire (GLIF) point neuronal models can be
successfully replicated on Loihi with an exceptional degree of precision.
Additionally, the high efficiency and low-power consumption of the
neuromorphic platform with increasing network size paves the way for a
complete replication of the mouse visual cortex dynamics, comprising hundreds
of thousands of neurons and millions of synapses.

## Acknowledgements

We thank Dr. Kael Dai for his helpful suggestions and advice regarding the use
of BMTK. We thank Intel for providing us with access to Loihi and the Intel
technical support team for helpful feedback on technical issues.

## References

  * Insel  _et al._ (2013) T. R. Insel, S. C. Landis, and F. S. Collins, Science 340, 687 (2013), https://science.sciencemag.org/content/340/6133/687.full.pdf .
  * Markram  _et al._ (2011) H. Markram, K. Meier, T. Lippert, S. Grillner, R. Frackowiak, S. Dehaene, A. Knoll, H. Sompolinsky, K. Verstreken, J. DeFelipe, S. Grant, J.-P. Changeux, and A. Saria, Procedia Computer Science 7, 39 (2011), proceedings of the 2nd European Future Technologies Conference and Exhibition 2011 (FET 11).
  * Winnubst  _et al._ (2019) J. Winnubst, E. Bas, T. A. Ferreira, Z. Wu, M. N. Economo, P. Edson, B. J. Arthur, C. Bruns, K. Rokicki, D. Schauder, D. J. Olbris, S. D. Murphy, D. G. Ackerman, C. Arshadi, P. Baldwin, R. Blake, A. Elsayed, M. Hasan, D. Ramirez, B. Dos Santos, M. Weldon, A. Zafar, J. T. Dudman, C. R. Gerfen, A. W. Hantman, W. Korff, S. M. Sternson, N. Spruston, K. Svoboda, and J. Chandrashekar, Cell 179, 268 (2019).
  * Lein  _et al._ (2007) E. Lein, M. Hawrylycz, N. Ao, M. Ayres, A. Bensinger, A. Bernard, A. Boe, M. Boguski, K. Brockway, E. Byrnes, L. Chen, L. Chen, T.-M. Chen, M. Chin, J. Chong, B. Crook, A. Czaplinska, C. Dang, S. Datta, and A. Jones, Nature 445, 168 (2007).
  * Davies  _et al._ (2018) M. Davies, N. Srinivasa, T. H. Lin, G. Chinya, Y. Cao, S. H. Choday, G. Dimou, P. Joshi, N. Imam, S. Jain, Y. Liao, C. K. Lin, A. Lines, R. Liu, D. Mathaikutty, S. McCoy, A. Paul, J. Tse, G. Venkataramanan, Y. H. Weng, A. Wild, Y. Yang, and H. Wang, IEEE Micro 38, 82 (2018).
  * DeBole  _et al._ (2019) M. DeBole, B. Taba, A. Amir, F. Akopyan, A. Andreopoulos, W. Risk, J. Kusnitz, C. O. Otero, T. Nayak, R. Appuswamy, P. J. Carlson, A. Cassidy, P. Datta, S. K. Esser, G. Garreau, K. L. Holland, S. Lekuch, M. Mastro, J. McKinstry, C. di Nolfo, B. Paulovicks, J. Sawada, K. Schleupen, B. Shaw, J. L. Klamo, M. Flickner, J. Arthur, and D. Modha, Computer 52, 20 (2019).
  * Löhr  _et al._ (2020) M. P. R. Löhr, C. Jarvers, and H. Neumann, 2020 2nd IEEE International Conference on Artificial Intelligence Circuits and Systems (AICAS) , 113 (2020).
  * Akopyan  _et al._ (2015) F. Akopyan, J. Sawada, A. Cassidy, R. Alvarez-Icaza, J. Arthur, P. Merolla, N. Imam, Y. Nakamura, P. Datta, G.-J. Nam, B. Taba, M. Beakes, B. Brezzo, J. B. Kuang, R. Manohar, W. Risk, B. Jackson, and D. Modha, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 34, 1537 (2015).
  * Khan  _et al._ (2008) M. M. Khan, D. R. Lester, L. A. Plana, A. Rast, X. Jin, E. Painkras, and S. B. Furber, in  _2008 IEEE International Joint Conference on Neural Networks (IEEE World Congress on Computational Intelligence)_ (2008) pp. 2849–2856.
  * Grübl  _et al._ (2020) A. Grübl, S. Billaudelle, B. Cramer, V. Karasenko, and J. Schemmel, “Verification and design methods for the brainscales neuromorphic hardware system,” (2020), arXiv:2003.11455 [cs.AR] .
  * Roy  _et al._ (2019) K. Roy, A. Jaiswal, and P. Panda, Nature 575, 607 (2019).
  * Thakur  _et al._ (2018) C. S. Thakur, J. L. Molin, G. Cauwenberghs, G. Indiveri, K. Kumar, N. Qiao, J. Schemmel, R. Wang, E. Chicca, J. Olson Hasler, J.-s. Seo, S. Yu, Y. Cao, A. van Schaik, and R. Etienne-Cummings, Frontiers in Neuroscience 12, 891 (2018).
  * Nawrocki  _et al._ (2016) R. A. Nawrocki, R. M. Voyles, and S. E. Shaheen, IEEE Transactions on Electron Devices 63, 3819 (2016).
  * Ou  _et al._ (2020) Q.-F. Ou, B.-S. Xiong, L. Yu, J. Wen, L. Wang, and Y. Tong, Materials 13, 3532 (2020).
  * Sharp and Furber (2013) T. Sharp and S. Furber, The 2013 International Joint Conference on Neural Networks (IJCNN) , 1 (2013).
  * Bhuiyan  _et al._ (2010) M. A. Bhuiyan, A. Nallamuthu, M. C. Smith, and V. K. Pallipuram, in  _2010 FOURTH INTERNATIONAL WORKSHOP ON HIGH-PERFORMANCE RECONFIGURABLE COMPUTING TECHNOLOGY AND APPLICATIONS (HPRCTA)_ (2010) pp. 1–9.
  * Rhodes  _et al._ (2019) O. Rhodes, L. Peres, A. Rowley, A. Gait, L. Plana, C. Y. A. Brenninkmeijer, and S. Furber, Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 378 (2019).
  * Knight and Furber (2016) J. C. Knight and S. Furber, Frontiers in Neuroscience 10 (2016).
  * Rossant  _et al._ (2010) C. Rossant, D. Goodman, J. Platkiewicz, and R. Brette, Frontiers in neuroinformatics 4, 2 (2010).
  * Nandi  _et al._ (2020) A. Nandi, T. Chartrand, W. V. Geit, A. Buchin, Z. Yao, S. Y. Lee, Y. Wei, B. Kalmbach, B. Lee, E. Lein, J. Berg, U. Sümbül, C. Koch, B. Tasic, and C. A. Anastassiou, bioRxiv (2020), 10.1101/2020.04.09.030239, https://www.biorxiv.org/content/early/2020/04/10/2020.04.09.030239.full.pdf .
  * Wang  _et al._ (2020) Q. Wang, S.-L. Ding, Y. Li, J. Royall, D. Feng, P. Lesnar, N. Graddis, M. Naeemi, B. Facer, A. Ho, T. Dolbeare, B. Blanchard, N. Dee, W. Wakeman, K. E. Hirokawa, A. Szafer, S. M. Sunkin, S. W. Oh, A. Bernard, J. W. Phillips, M. Hawrylycz, C. Koch, H. Zeng, J. A. Harris, and L. Ng, Cell 181, 936 (2020).
  * Dai  _et al._ (2020) K. Dai, S. L. Gratiy, Y. N. Billeh, R. Xu, B. Cai, N. Cain, A. E. Rimehaug, A. J. Stasik, G. T. Einevoll, S. Mihalas, C. Koch, and A. Arkhipov, PLOS Computational Biology 16, 1 (2020).
  * Teeter  _et al._ (2018) C. Teeter, R. Iyer, V. Menon, N. Gouwens, D. Feng, J. Berg, A. Szafer, N. Cain, H. Zeng, M. Hawrylycz, _et al._ , Nature communications 9, 709 (2018).
  * Gutzen  _et al._ (2018) R. Gutzen, M. von Papen, G. Trensch, P. Quaglio, S. Grün, and M. Denker, Frontiers in neuroinformatics 12 (2018).
  * Brette  _et al._ (2008) R. Brette, M. Lilith, T. Carnevale, M. Hines, D. Beeman, J. Bower, M. Diesmann, A. Morrison, P. Goodman, F. Harris, M. Zirpe, T. Natschläger, D. Pecevski, B. Ermentrout, M. Djurfeldt, A. Lansner, O. Rochel, T. Viéville, E. Muller, and A. Destexhe, Journal of computational neuroscience 23, 349 (2008).
  * Kunkel  _et al._ (2017) S. Kunkel, A. Morrison, P. Weidel, J. M. Eppler, A. Sinha, W. Schenck, M. Schmidt, S. B. Vennemo, J. Jordan, A. Peyser, D. Plotnikov, S. Graber, T. Fardet, D. Terhorst, H. Mørk, G. Trensch, A. Seeholzer, R. Deepu, J. Hahne, I. Blundell, T. Ippen, J. Schuecker, H. Bos, S. Diaz, E. Hagen, C. Mahmoudian, Sepehr andBachmann, M. E. Lepperød, O. Breitwieser, B. Golosio, H. Rothe, H. Setareh, M. Djurfeldt, T. Schumann, A. Shusharin, J. Garrido, E. B. Muller, A. Rao, J. H. Vieites, and H. E. Plesser, “Nest 2.12.0,”  (2017).
  * Gerstner and Kistler (2002) W. Gerstner and W. M. Kistler,  _Spiking Neuron Models: Single Neurons, Populations, Plasticity_ (Cambridge University Press, 2002).
  * Rotter and Diesmann (1999) S. Rotter and M. Diesmann, Biological Cybernetics 81, 381 (1999).
  * Xin Jin  _et al._ (2008) Xin Jin, S. B. Furber, and J. V. Woods, in  _2008 IEEE International Joint Conference on Neural Networks (IEEE World Congress on Computational Intelligence)_ (2008) pp. 2812–2819.
  * Mikaitis  _et al._ (2018) M. Mikaitis, D. Lester, D. Shang, S. Furber, G. Liu, J. Garside, S. Scholze, S. Höppner, and A. Dixius (2018).
  * Trensch  _et al._ (2018) G. Trensch, R. Gutzen, I. Blundell, M. Denker, and A. Morrison, Frontiers in Neuroinformatics 12, 81 (2018).
  * Herz  _et al._ (2006) A. V. M. Herz, T. Gollisch, C. K. Machens, and D. Jaeger, Science 314, 80 (2006), https://science.sciencemag.org/content/314/5796/80.full.pdf .
  * Hopkins and Furber (2015) M. Hopkins and S. Furber, Neural Computation 27, 2148 (2015).


$\langle$reception date$\rangle$ $\langle$acception date$\rangle$
$\langle$publication date$\rangle$

††affiliationtext: Research Institute of Stellar Explosive Phenomena, Fukuoka
University, 8-19-1 Nanakuma, Jonan-ku, Fukuoka-shi, Fukuoka 814-0180, Japan

astroparticle physics — elementary particles — supernovae: general — ISM:
supernova remnants

# Heavy axion-like particles and MeV decay photons from nearby type Ia
supernovae

Kanji Mori<EMAIL_ADDRESS>
###### Abstract

Axion-like particles (ALPs) are hypothetical bosons which may couple with
photons. Since many ALPs can be emitted from hot and dense astrophysical
plasma, nearby supernovae (SNe) are a possible probe into their properties
including the ALP mass $m_{a}$ and the coupling constant $g_{a\gamma}$ between
ALPs and photons. I calculated ALP emission from a type Ia SN (SN Ia) model
with the near-Chandrasekhar mass. It is found that the ALP luminosity from SNe
Ia reaches $\sim 10^{43}(g_{a\gamma}/10^{-10}\;\mathrm{GeV}^{-1})^{2}$ erg s-1
if $m_{a}\lesssim 1$ MeV. Heavy ALPs emitted from SNe are unstable and decay
into photons. I predict the time delay and the flux of decay photons that
reach Earth from a nearby SN Ia. It is found that the decay photons may
provide a constraint on $g_{a\gamma}$ which is as stringent as an SN 1987A
limit if an SN Ia is located 1 kpc away or closer and next-generation MeV
$\gamma$-ray satellites observe it $\sim 1-10$ years after the explosion.

##  1 Introduction

Axion-like particles (ALPs) are a class of hypothetical bosons beyond the
standard model of particle physics (e.g. [Choi, Im, & Shin (2020), Di Luzio et
al. (2020)]). Axions are pseudo Nambu-Goldstone bosons originally introduced
to solve the strong CP problem of quantum chromodynamics (Weinberg, 1978;
Wilczek, 1978). ALPs, which share similar properties with axions, emerge as a
result of string theory compactifications (Svrcek & Witten, 2006; Arvanitaki
et al., 2010). It is notable that particle models with ALPs may solve the
electroweak hierarchy problem (Graham, Kaplan, & Rajendran, 2015).

The interaction between an ALP and photons is described by the Lagrangian
(Raffelt & Stodolsky, 1988)

$\mathcal{L}_{a\gamma\gamma}=-\frac{1}{4}g_{a\gamma}a\tilde{F}^{\mu\nu}F_{\mu\nu},$
(1)

where $g_{a\gamma}$ is the coupling constant, $a$ is the ALP field,
$F_{\mu\nu}$ is the electromagnetic tensor, and $\tilde{F}_{\mu\nu}$ is its
dual. The mass $m_{a}$ and the coupling $g_{a\gamma}$ of heavy ALPs with
$m_{a}\gtrsim 10$ keV have been explored by beam dump experiments (Jaeckel &
Spannowsky, 2016; Dolan et al., 2017; Döbrich, Jaeckel, & Spadaro, 2019;
Banerjee et al., 2020) and astronomical objects including globular clusters
(Carenza et al., 2020), white dwarfs (Dolan, Hiskens, & Volkas, 2021), and the
neutrino signal (Lee, 2018; Lucente et al., 2020) and $\gamma$-rays
(Giannotti, Duffy, & Nita, 2011; Jaeckel, Malta, & Redondo, 2018) from SN
1987A. Also, implications of heavy ALPs on cosmological phenomena including
Big Bang nucleosynthesis and the cosmic microwave background have been
discussed (Cadamuro et al., 2011; Depta, Hufnagel, & Schmidt-Hoberg, 2020).

Since many ALPs may be produced in hot plasma in supernovae (SNe) depending on
$m_{a}$ and $g_{a\gamma}$, a nearby SN can be used as a probe of ALPs.
Although core-collapse SNe produce more ALPs than type Ia SNe (SNe Ia) do
because of the higher temperature, the next nearby SN after SN 1987A can be
type Ia (Li et al., 2011). It is therefore desirable to predict the ALP
emission from SNe Ia to discuss a possible constraint that can be obtained
from the event.

A heavy ALP emitted from an SN decays into two photons during propagation to
the Solar System. It is hence possible to constrain the ALP parameters by
$\gamma$-ray observations of a nearby astronomical object which emits a lot of
ALPs. In this paper, I calculate the ALP production in a nearby SN Ia and
predict the $\gamma$-ray flux which originates from the ALP decay.

This paper is organized as follows. In Section 2, I describe the prescription
of the ALP production and decay and the SN Ia model. In Section 3, the
intrinsic ALP emission from the SN Ia is shown. In Section 4, the $\gamma$-ray
flux on Earth is estimated and its observability is discussed. In Section 5,
the results are summarized and future prospects are explained.

##  2 Method

In this Section, I describe the ALP production and decay rates and the SN Ia
model adopted in this work.

###  2.1 ALP production

In this study, we consider the Primakoff process (Primakoff, 1951)

$\displaystyle\gamma+e^{-}$ $\displaystyle\rightarrow$ $\displaystyle a+e^{-}$
(2) $\displaystyle\gamma+Ze$ $\displaystyle\rightarrow$ $\displaystyle a+Ze$
(3)

induced by electrons and ions and photon coalescence

$\displaystyle\gamma+\gamma\rightarrow a$ (4)

as the ALP production processes.

The Primakoff rate is given by (di Lella et al., 2000; Carenza et al., 2020;
Lucente et al., 2020)

$\displaystyle\Gamma_{\gamma\rightarrow
a}=g_{a\gamma}^{2}\frac{T\kappa^{2}}{32\pi}\frac{p}{E}\left(\frac{((k+p)^{2}+\kappa^{2})((k-p)^{2}+\kappa^{2})}{4kp\kappa^{2}}\right.$
$\displaystyle\left.\times\ln\left(\frac{(k+p)^{2}+\kappa^{2}}{(k-p)^{2}+\kappa^{2}}\right)-\frac{(k^{2}-p^{2})^{2}}{4kp\kappa^{2}}\ln\left(\frac{(k+p)^{2}}{(k-p)^{2}}\right)-1\right),$
(5)

where $p=\sqrt{E_{a}^{2}-m_{a}^{2}}$ and
$k=\sqrt{\omega^{2}-\omega_{\mathrm{pl}}^{2}}$ are the momenta and $E_{a}$ and
$\omega$ are the energies of ALPs and photons, respectively. The energy
conservation leads to a relation $E_{a}=\omega$. The plasma frequency is
denoted by $\omega_{\mathrm{pl}}$ and

$\displaystyle\kappa=\sqrt{\frac{4\pi\alpha n^{\mathrm{eff}}}{T}}$ (6)

is the inverse Debye-Hückel length, where $\alpha$ is the fine structure
constant and
$n^{\mathrm{eff}}=n_{e}^{\mathrm{eff}}+n^{\mathrm{eff}}_{\mathrm{ion}}$ is the
sum of the effective electron and ion number densities. The plasma frequency
is evaluated as (Chapter 6 and Appendix D in Raffelt (1996))

$\displaystyle\omega_{\mathrm{pl}}=\sqrt{\frac{4\pi\alpha
n_{e}}{E_{\mathrm{F}}}}\approx
28.7\;\mathrm{eV}\frac{(Y_{e}\rho)^{\frac{1}{2}}}{(1+(1.019\times
10^{-6}Y_{e}\rho)^{\frac{2}{3}})^{\frac{1}{4}}},$ (7)

where $n_{e}$ is the electron number density, $E_{\mathrm{F}}$ is the electron
Fermi energy, $Y_{e}$ is the electron mole fraction, and $\rho$ is the density
in units of g cm-3. Equation (7) is an asymptotic formula for degenerate
plasma and holds for both non-relativistic ($p_{\mathrm{F}}\ll m_{e}$) and
relativistic ($p_{\mathrm{F}}\gg m_{e}$) cases, where $m_{e}$ is the electron
mass and $p_{\mathrm{F}}$ is the electron Fermi momentum. Since the gas is
always degenerate in the region of interest, $\omega_{\mathrm{pl}}$ is
reasonably approximated by the equation.

The effective ion density is defined as

$\displaystyle
n^{\mathrm{eff}}_{\mathrm{ion}}=\frac{\rho}{m_{\mathrm{u}}}\sum_{j}\frac{Z_{j}^{2}X_{j}}{A_{j}},$
(8)

where $\rho$ is the density, $m_{\mathrm{u}}$ is the atomic mass unit, and
$Z_{j}$, $X_{j}$, and $A_{j}$ are the atomic number, the mass fraction, and
the mass number of the $j$-th nucleus, respectively.

The effective electron density $n_{e}^{\mathrm{eff}}$ takes degeneracy into
account. It is calculated as $n_{e}^{\mathrm{eff}}=R_{\mathrm{deg}}n_{e}$,
where $n_{e}=\rho Y_{e}/m_{\mathrm{u}}$, and $R_{\mathrm{deg}}$ is a
coefficient defined as follows (Payez et al., 2015; Straniero et al., 2019).
We introduce a parameter $z=\rho Y_{e}/T^{\frac{3}{2}}$ where $T$ is the
temperature. When $z$ with $T$ in K and $\rho$ in g cm-3 is less than or equal
to $5.45\times 10^{-11}$, $R_{\mathrm{deg}}=1$. When $z\in(5.45\times
10^{-11},\;7.2\times 10^{-8}]$,

$\displaystyle
R_{\mathrm{deg}}=0.63+0.3\arctan\left(0.65-9316z^{0.48}+\frac{0.019}{z^{0.212}}\right).$
(9)

Finally, when $z>7.2\times 10^{-8}$,

$\displaystyle R_{\mathrm{deg}}=4.78\times 10^{-6}z^{-0.667}.$ (10)

The ALP production rate per a unit volume and time is then written as

$\displaystyle Q_{a}$ $\displaystyle=$ $\displaystyle
2\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}}\Gamma_{\gamma\rightarrow a}\omega
f(\omega)$ (11) $\displaystyle=$
$\displaystyle\frac{1}{\pi^{2}}\int^{\infty}_{m_{a}}d\omega\omega^{2}\sqrt{\omega^{2}-\omega_{\mathrm{pl}}^{2}}\Gamma_{\gamma\rightarrow
a}f(\omega),$

where $f(\omega)$ is the photon Bose-Einstein distribution. The factor 2 in
the second term stems from the degrees of freedom of photons.

The ALP production rate due to the photon coalescence is given by (di Lella et
al., 2000; Carenza et al., 2020; Lucente et al., 2020)

$\frac{d^{2}n_{a}}{dtdE}=g_{a\gamma}^{2}\frac{m_{a}^{4}}{128\pi^{3}}p\left(1-\frac{4\omega_{\mathrm{pl}}^{2}}{m_{a}^{2}}\right)^{\frac{3}{2}}e^{-\frac{E_{a}}{T}},$
(12)

and the ALP production rate for this reaction is then

$Q_{a}=\int^{\infty}_{m_{a}}dE_{a}E_{a}\frac{d^{2}n_{a}}{dtdE_{a}}.$ (13)

It should be noted that this reaction is possible only when the kinematic
condition $m_{a}>2\omega_{\mathrm{pl}}$ is satisfied.

When the ALP kinetic energy is smaller than the gravitational energy, the
particle is trapped (Dreiner et al., 2003; Lucente et al., 2020). In this
case, the decay photons cannot be observed. However, the gravitational energy
$E_{\mathrm{G}}=GMm_{a}/r$ between the star and an ALP is only $\sim 7$ keV,
where $G$ is the gravitational constant, $M\sim 1.4M_{\odot}$ is the total
mass of the star, and $r$ is the stellar radius. I here assumed $r\sim 3\times
10^{3}$ km. Since the gravitational energy is much less than ALP energies
($\sim$ MeV), the trapping effect is negligible.

###  2.2 ALP decay

Heavy ALPs are unstable particles. The mean free path (MFP) due to the ALP
decay ($a\rightarrow\gamma\gamma$) is (Jaeckel, Malta, & Redondo, 2018;
Lucente et al., 2020)

$\displaystyle\lambda_{a\rightarrow\gamma\gamma}(E)=\frac{\beta\gamma}{\Gamma_{a\rightarrow\gamma\gamma}}$
(14)

where $\gamma$ is the Lorentz factor of ALPs, $\beta=\sqrt{1-\gamma^{-2}}$,
and

$\Gamma_{a\rightarrow\gamma\gamma}=g_{a\gamma}^{2}\frac{m_{a}^{3}}{64\pi}.$
(15)

The MFP is then evaluated as

$\lambda_{a\rightarrow\gamma\gamma}\sim
0.4\;\mathrm{ly}\;g_{10}^{-2}\left(\frac{E_{a}}{1\;\mathrm{MeV}}\right)\left(\frac{m_{a}}{1\;\mathrm{MeV}}\right)^{-4},$
(16)

where $g_{10}=g_{a\gamma}/(10^{-10}\;\mathrm{GeV}^{-1})$.

###  2.3 SN Ia model

Figure 1: The profiles of (a) temperature, (b) density, and (c) plasma
frequency in the SN Ia model as a function of the mass coordinate $M_{r}$.

SNe Ia are a thermonuclear explosion of a carbon-oxygen white dwarf (WD). In
near-Chandrasekhar mass models of SNe Ia, accretion occurs on the WD and its
mass becomes close to the Chandrasekhar mass. As a result, carbon fusion is
ignited at the center and the whole star becomes unbound. In this study, we
adopt W7 (Nomoto, Thielemann, & Yokoi, 1984) as the SN Ia model, which is a
one-dimensional model with pure deflagration.

Figure 1 shows the profiles of the temperature, the density, and the plasma
frequency in the W7 model. Here $M_{r}$ is the mass coordinate defined as

$\displaystyle M_{r}=4\pi\int^{r}_{0}\rho(r^{\prime})r^{\prime 2}dr^{\prime}$
(17)

and $T_{9}$ is the temperature in units of $10^{9}$ K. The deflagration wave
starts at the center and propagates toward the surface. When the flame passes,
the material is heated and the temperature rises. The initial density is as
high as $\sim 10^{9}$ g cm-3 because the star is close to the Chandrasekhar
mass. Since the WD expands as the star explodes, the density decreases as a
function of time.

The Primakoff rate shown in equation (5) is dependent on the chemical
composition. In this calculation, I adopt the nuclear statistical equilibrium
(NSE) prescription (Seitenzahl et al., 2009) when $T_{9}\geq 3$ and carbon-
oxygen-neon plasma with $X(^{12}\mathrm{C})=0.475$,
$X(^{16}\mathrm{O})=0.500$, and $X(^{22}\mathrm{Ne})=0.015$ when $T_{9}<3$.
The NSE is a good approximation for hot plasma with $T_{9}\geq 3$ where the
majority of ALPs is produced. Because the NSE approximation is not valid when
$T_{9}<3$, the composition is switched to the carbon-oxygen-neon plasma for
cooler gas. The NSE table shown by Seitenzahl et al. (2009) considers 443
nuclei from neutron to krypton. All of these nuclei are considered in the
calculation of the Primakoff rate. In the NSE, the nuclear abundances are
determined only by $T$, $\rho$, and $Y_{e}$. I adopt a constant electron
fraction $Y_{e}=0.4575$ on the basis of a network calculation (Mori et al.,
2020).

##  3 Intrinsic ALP production in the SN

Figure 2: The ALP production rate per a unit volume and time as a function of
the mass coordinate. The coupling constant and the ALP mass are fixed to
$g_{10}=1$ and $m_{a}=0.1$ MeV, respectively. Figure 3: The ALP luminosity as
a function of time. The coupling constant is fixed to $g_{10}=1$.

Using the SN model shown in figure 1 and equation (11), we can calculate the
ALP emissivity from an SN Ia. Figure 2 shows the profile of the ALP emissivity
$Q_{a}$ defined in equations (11) and (13). In this figure, the ALP mass is
fixed to $m_{a}=0.1$ MeV and the coupling constant is fixed to $g_{10}=1$.
Since $Q_{a}$ rapidly increases as a function of $T$, ALPs are emitted behind
the flame. Although $Q_{a}$ becomes smaller as the SN cools, the region where
ALPs are emitted becomes wider because of the propagation of the flame.

In figure 3, I show the ALP luminosity

$L_{a}=4\pi\int Q_{a}(r)r^{2}dr,$ (18)

where the integral is performed from the center to the surface of the star. I
adopt $g_{10}=1$ in this plot. It is seen that $L_{a}$ reaches its peak at
$0.5-1$ s after the ignition. When $m_{a}\lesssim 1$ MeV, $L_{a}$ reaches
$\sim 10^{43}$ erg s-1 at its peak. Coincidentally, this value is close to
peak photon luminosities of a typical SN Ia. The ALP luminosity is much
smaller when $m_{a}\gtrsim 5$ MeV because the thermal photons are not
energetic enough to produce such heavy ALPs.

In the case of $m_{a}=0.1$ MeV, we can see a bend in $L_{a}$ at $t\sim 1.4$ s.
This is an artifact that originates from a sudden switch from the NSE to the
carbon-oxygen-neon plasma and will disappear if dynamical nucleosynthetic
processes are considered.

Figure 4: The contribution of the Primakoff process and photon coalescence on
the ALP luminosity as a function of the ALP mass. The result at $t=0.90$ s
after the ignition is shown. The coupling constant is fixed to $g_{10}=1$.

Figure 3 is the total ALP luminosity induced by the Primakoff process and
photon coalescence. Figure 4 shows the contribution of each process at
$t=0.90$ s as a function of the ALP mass. It is seen that the Primakoff
process is dominant if $m_{a}\lesssim 2$ MeV while photon coalescence is
dominant if $m_{a}\gtrsim 2$ MeV.

In figures 2, 3 and 4, $g_{10}=1$ is assumed. However, it is notable that
$Q_{a}$ and $L_{a}$ are simply proportional to $g_{10}^{2}$ as we can see from
equations (5) and (12).

##  4 Decay photons observed on Earth

As described in Section 2.2, a heavy ALP can decay into two photons, which can
be detected on Earth by $\gamma$-ray telescopes. In this Section, I predict
the expected number and the flux of the decay photons.

Figure 5: A schematic picture of geometry between a nearby SN and Earth.

The configuration discussed here is shown in figure 5. In this Section, we
assume that the distance to the SN Ia is $d_{\mathrm{SN}}=1$ kpc if not stated
otherwise. The distance which an ALP travels before its decay is denoted by
$L_{1}$ and the distance which the decay photon travels is denoted by $L_{2}$.
The angle between the paths of the ALP and the photon is denoted by $\alpha$.
It is easily shown that

$L_{2}=-L_{1}\cos\alpha+\sqrt{d_{\mathrm{SN}}^{2}-L_{1}^{2}\sin^{2}\alpha}.$
(19)

The angular distribution of decay photons is isotropic if an ALP is in its
rest frame. However, since the SN ALP travels at a relativistic velocity, the
decay photons are beamed toward the forward direction in our frame. The most
likely angle is $\alpha\sim\gamma^{-1}$ (Jaeckel, Malta, & Redondo, 2018;
Caputo et al., 2021).

Since ALPs are massive, they arrive on Earth later than photons. The time
delay is written as

$t(E_{\gamma})=\frac{L_{1}}{\beta}+L_{2}-d_{\mathrm{SN}},$ (20)

where $E_{\gamma}=E_{a}/2$ is the photon energy. Because $L_{1}$, $L_{2}$, and
$\beta$ are dependent on energy, the decay photons with different energies
arrive on Earth at different times. In the region of the parameter space
discussed in this study, $L_{1}\approx\lambda_{a\rightarrow\gamma\gamma}$ is
always shorter than $d_{\mathrm{SN}}$.

Figure 6 shows the time delay $t$ in the case of $m_{a}=1$ MeV as a function
of $E_{\gamma}=E_{a}/2$. It is seen that more energetic photons arrive earlier
regardless of $m_{a}$ because energetic ALPs run faster. The delay time
becomes shorter when $g_{10}$ is larger since the MFP is shorter when $g_{10}$
is larger and hence the signal runs at the speed of light longer. The low-
energy cutoff in the plot is the lowest possible photon energy which is
equivalent to $m_{a}/2$.

Figure 6: The delay time $t(E_{\gamma})$ of decay photons observed on Earth.
The ALP mass is fixed to $m_{a}=1$ MeV.

Given the delay time of decay photons, it is possible to estimate the energy
flux observed on Earth which is defined as

$\displaystyle F(t)=E_{\gamma}(t)\frac{dn_{\gamma}(t)}{dt},$ (21)

where $dn_{\gamma}/dt$ is the number of photons with the energy $E_{\gamma}$
per a unit area and time. The photon number is calculated as (Caputo et al.,
2021)

$\displaystyle\frac{dn_{\gamma}}{dt}=\frac{2}{4\pi
d_{\mathrm{SN}}^{2}}\frac{dn_{a}}{dE_{a}}\left(\frac{dt}{dE_{a}}\right)^{-1}e^{-\frac{R_{\mathrm{SN}}}{\lambda_{a\rightarrow\gamma\gamma}(E_{a})}}$
$\displaystyle\times\left(1-e^{-\frac{d_{\mathrm{SN}}}{\lambda_{a\rightarrow\gamma\gamma}(E_{a})}}\right).$
(22)

The overall factor 2 corresponds to two decay photons that are produced by an
ALP. The factor $dn_{a}/dE_{a}$ is the spectrum of ALPs emitted from the SN
and $dt/dE_{a}$ is the energy derivative of the delay time. The factor
$\exp(-R_{\mathrm{SN}}/\lambda_{a\rightarrow\gamma\gamma}(E_{a}))$ is the
survival probability of ALPs inside the SN, where $R_{\mathrm{SN}}=3\times
10^{3}$ km is assumed on the basis of the mass-radius relation of massive
white dwarfs (Hamada & Salpeter, 1961) If an ALP decays in the SN, the decay
photon cannot escape to the interstellar space. Another factor
$1-\exp(-{d_{\mathrm{SN}}}/{\lambda_{a\rightarrow\gamma\gamma}(E_{a})})$ is
the ALP decay probability between the SN and Earth. Since we are focusing on
relativistic ALPs, the decay photons are strongly beamed toward the forward
direction. We therefore ignore photons that propagate backward with this
factor.

Figure 7 shows the flux $F$ of decay photons as a function of time after the
SN explosion. It is seen that the photon flux is diluted in time and the peak
is expected 1-10 years after the explosion. The flux is larger if the ALP-
photon coupling is larger because the intrinsic ALP emission from an SN
increases as a function of $g_{10}$. The arrival time of the peak becomes
earlier if $g_{10}$ is larger as we saw in figure 6. The photon flux and the
delay time at the peak are summarized in table 1.

Figure 7: The flux $F(t)$ of decay photons observed on Earth. The ALP mass is fixed to $m_{a}=1$ MeV. $g_{10}$ | Delay Time [yr] | Peak Flux [erg cm-2 s-1] |  $E_{\gamma}$ [MeV] |  $\theta$ [deg]
---|---|---|---|---
0.1 | 12 | $1.2\times 10^{-12}$ | 1.6 | 0.55
0.2 | 3.1 | $1.9\times 10^{-11}$ | 1.6 | 0.14
0.3 | 1.4 | 9.7$\times 10^{-11}$ | 1.6 | 0.061
Table 1: The summary on decay photons and the angular radius of the halo when
the flux reaches its peak. The ALP mass is fixed to 1 MeV. Figure 8: The
spectrum $dn_{a}/dE_{a}$ of ALPs emitted from the SN. The coupling constant is
fixed to $g_{10}=1$.

Figure 8 shows the spectrum of ALPs produced in the SN as a function of the
ALP energy $E_{a}$. In the case of $m_{a}=1$ MeV, the peak in the ALP spectrum
is at $E_{a}\sim 1.5$ MeV, while table 1 shows that ALPs that produce the
photon flux peak are as energetic as $E_{a}\sim 3.2$ MeV. The difference stems
from the time delay. As seen from figure 6, low-energy ALPs arrive on Earth
later and thus the flux of their decay photons is diluted in time. As a
result, the peak in the photon flux is more energetic than the peak in the ALP
spectrum.

Since the incidence direction of ALPs does not coincide with the line of sight
between the SN and Earth, the decay photons would form a halo around the SN
when they are observed on Earth. When $L_{1}$ is close to $d_{\mathrm{SN}}$,
the size of the halo can be so large that it exceeds a finite angular
acceptance of detectors (Jaeckel, Malta, & Redondo, 2018). When $m_{a}$ is
fixed, the halo is larger if $g_{10}$ is smaller because $L_{1}$ is longer. It
is easily shown from figure 5 that the angular radius $\theta$ of the halo is
written as

$\sin\theta=\frac{L_{1}}{d_{\mathrm{SN}}}\sin\alpha.$ (23)

The values of $\theta$ at the peak of the photon flux are shown in table 1. It
is seen that the halo is larger than Moon if $g_{10}\lesssim 0.1$. Such a
large halo might make observations of decay photons challenging.

The photon flux $\sim 1$ year after the SN event is as high as $\sim 10^{-11}$
erg cm-2 s-1 if $g_{10}\gtrsim 0.3$ and $m_{a}=1$ MeV. We can see from figure
6 that $E_{\gamma}=O(\mathrm{MeV})$ at the peak of the flux. The sensitivity
of instruments for MeV photons on current $\gamma$-ray satellites is not
enough to detect the signal because the instrumental background prevents one
from improving it. (Schönfelder, 2004). However, next-generation instruments
including the Electron-Tracking Compton Camera (ETCC: Tanimori et al. (2015))
and All-sky Medium Energy Gamma-ray Observatory (AMEGO: McEnery et al. (2019))
plan to achieve sensitivities as high as $\sim 10^{-12}-10^{-11}$ erg cm-2
s-1. Observing MeV $\gamma$-rays from a nearby SN Ia with the next-generation
instruments may give an independent constraint on the ALP parameters.

A recent nearby core-collapse SN, namely SN 1987A, has provided a stringent
constraint on ALPs with $m_{a}\sim 10$ keV-100 MeV (Giannotti, Duffy, & Nita,
2011; Jaeckel, Malta, & Redondo, 2018). Jaeckel, Malta, & Redondo (2018) used
an upper limit on the $\gamma$-ray fluence obtained by the Solar Maximum
Mission (Chupp, Vestrand, & Reppin, 1989) to constrain $g_{10}$. They
concluded that $g_{10}\lesssim 0.2$ for ALPs with $m_{a}=1$ MeV. If the ALP
parameters are on this limit, the $\gamma$-ray flux from a nearby SN Ia can be
as high as $10^{-12}-10^{-11}$ erg cm-2 s-1, as can be seen from figure 7.
Since this flux may be observable with next-generation instruments, SNe Ia can
potentially provide an independent constraint which is comparable to the
constraints obtained with core-collapse SNe.

##  5 Summary and Discussion

In this study, I calculated the emission of heavy ALPs from a nearby SN Ia for
the first time. It is found that the ALP luminosity reaches $L_{a}\sim
10^{43}g_{10}^{2}$ erg s-1 for $m_{a}\lesssim 1$ MeV. The emission is
suppressed for heavier ALPs because thermal photons are not energetic enough.
I then estimated the flux of MeV decay photons from an SN Ia that is located 1
kpc away. The photons observed on Earth are diluted in time: energetic photons
arrive earlier. The flux can be as high as $10^{-12}-10^{-11}$ erg cm-2 s-1 if
the ALP parameters are on the SN 1987A limit (Jaeckel, Malta, & Redondo,
2018). Next-generation MeV $\gamma$-ray telescopes including ETCC and AMEGO
may provide an independent constraint on $m_{a}$ and $g_{10}$ which can be as
stringent as the SN 1987A limit if an SN Ia happens $\sim 1$ kpc away or
closer.

In my calculation, a near-Chandrasekhar mass model is adopted (Nomoto,
Thielemann, & Yokoi, 1984). However, the nature of SN Ia progenitors is still
under debate. Recently, sub-Chandrasekhar mass progenitor models have been
developed as well (Fink et al., 2010; Woosley & Kasen, 2011; Shen et al.,
2018; Leung & Nomoto, 2020). Since the central density in a sub-Chandrasekhar
mass WD can be as low as $\sim 10^{7}$ g cm-3, the ALP emissivity can be
significantly affected. It is hence desirable to perform a systematic study
that focuses on the model dependence.

In this study, SNe are not assumed to be magnetized. However, magnetic fields
in a white dwarf can be as strong as $\sim 100$ MG (Schmidt et al., 2003). In
such a case, the flux of decay photons can be significantly higher than the
unmagnetized case because of two reasons. One is that plasmon-ALP conversion
in plasma by external magnetic fields can enhance the intrinsic ALP production
in a white dwarf (Mikheev, Raffelt, & Vassilevskaya, 1998; Caputo, Millar, &
Vitagliano, 2020; Caputo et al., 2021). Another reason is that magnetic fields
can convert ALPs into photons to enhance the photon flux (Dessert, Long, &
Safdi, 2019). Systematic studies on the effects of magnetic fields are
desirable.

It has been reported that neutrinos from a nearby SN Ia may be observable and
provide information on the explosion mechanism (Kunugise & Iwamoto, 2007;
Odrzywolek & Plewa, 2011; Wright et al., 2016). The predicted peak luminosity
of neutrinos reaches $\sim 10^{50}$ erg s-1, which is much higher than the ALP
luminosity. The ALP emission is thus not likely to affect the neutrino
emission. Multi-messenger astronomy of SNe Ia would shed light on both of the
nature of SNe and fundamental physics.

The author thanks Kei Kotake and Tomoya Takiwaki for fruitful discussions.
This work is supported by Research Institute of Stellar Explosive Phenomena at
Fukuoka University.

## References

  * Arvanitaki et al. (2010) Arvanitaki, A., Dimopoulos, S., Dubovsky, S., et al. 2010, Phys. Rev. D, 81, 123530
  * Banerjee et al. (2020) Banerjee, D., Bernhard, J., Burtsev, V. E., et al. 2020, Phys. Rev. Lett., 125, 081801
  * Cadamuro et al. (2011) Cadamuro, D., Hannestad, S., Raffelt, G., et al. 2011, J. Cosmology Astropart. Phys, 2011, 003
  * Caputo, Millar, & Vitagliano (2020) Caputo, A., Millar, A. J., & Vitagliano, E. 2020, Phys. Rev. D, 101, 123004
  * Caputo et al. (2021) Caputo, A., Carenza, P., Lucente, G., et al. 2021, arXiv:2104.05727
  * Carenza et al. (2020) Carenza, P., Straniero, O., Döbrich, B., et al. 2020, Physics Letters B, 809, 135709
  * Choi, Im, & Shin (2020) Choi, K., Im, S. H., & Shin, C. S. 2020, arXiv:2012.05029
  * Chupp, Vestrand, & Reppin (1989) Chupp, E. L., Vestrand, W. T., & Reppin, C. 1989, Phys. Rev. Lett., 62, 505
  * Depta, Hufnagel, & Schmidt-Hoberg (2020) Depta, P. F., Hufnagel, M., & Schmidt-Hoberg, K. 2020, J. Cosmology Astropart. Phys, 2020, 009
  * Dessert, Long, & Safdi (2019) Dessert, C., Long, A. J., & Safdi, B. R. 2019, Phys. Rev. Lett., 123, 061104
  * di Lella et al. (2000) di Lella, L., Pilaftsis, A., Raffelt, G., et al. 2000, Phys. Rev. D, 62, 125011
  * Di Luzio et al. (2020) Di Luzio, L., Giannotti, M., Nardi, E., et al. 2020, Phys. Rep., 870, 1
  * Döbrich, Jaeckel, & Spadaro (2019) Döbrich, B., Jaeckel, J., & Spadaro, T. 2019, Journal of High Energy Physics, 2019, 213
  * Dolan et al. (2017) Dolan, M. J., Ferber, T., Hearty, C., et al. 2017, Journal of High Energy Physics, 2017, 94
  * Dolan, Hiskens, & Volkas (2021) Dolan, M. J., Hiskens, F. J., & Volkas, R. R. 2021, arXiv:2102.00379
  * Dreiner et al. (2003) Dreiner, H. K., Hanhart, C., Langenfeld, U., et al. 2003, Phys. Rev. D, 68, 055004
  * Fink et al. (2010) Fink, M., Röpke, F. K., Hillebrandt, W., et al. 2010, A&A, 514, A53
  * Giannotti, Duffy, & Nita (2011) Giannotti, M., Duffy, L. D., & Nita, R. 2011, J. Cosmology Astropart. Phys, 2011, 015
  * Graham, Kaplan, & Rajendran (2015) Graham, P. W., Kaplan, D. E., & Rajendran, S. 2015, Phys. Rev. Lett., 115, 221801
  * Hamada & Salpeter (1961) Hamada, T. & Salpeter, E. E. 1961, ApJ, 134, 683
  * Jaeckel & Spannowsky (2016) Jaeckel, J. & Spannowsky, M. 2016, Physics Letters B, 753, 482
  * Jaeckel, Malta, & Redondo (2018) Jaeckel, J., Malta, P. C., & Redondo, J. 2018, Phys. Rev. D, 98, 055032
  * Kunugise & Iwamoto (2007) Kunugise, T. & Iwamoto, K. 2007, PASJ, 59, L57
  * Lee (2018) Lee, J. S. 2018, arXiv:1808.10136
  * Leung & Nomoto (2020) Leung, S.-C. & Nomoto, K. 2020, ApJ, 888, 80
  * Li et al. (2011) Li, W., Chornock, R., Leaman, J., et al. 2011, MNRAS, 412, 1473
  * Lucente et al. (2020) Lucente, G., Carenza, P., Fischer, T., et al. 2020, J. Cosmology Astropart. Phys, 2020, 008
  * McEnery et al. (2019) McEnery, J., van der Horst, A., Dominguez, A., et al. 2019, BAAS
  * Mikheev, Raffelt, & Vassilevskaya (1998) Mikheev, N. V., Raffelt, G., & Vassilevskaya, L. A. 1998, Phys. Rev. D, 58, 055008
  * Mori et al. (2020) Mori, K., Suzuki, T., Honma, M., et al. 2020, ApJ, 904, 29
  * Nomoto, Thielemann, & Yokoi (1984) Nomoto, K., Thielemann, F.-K., & Yokoi, K. 1984, ApJ, 286, 644
  * Odrzywolek & Plewa (2011) Odrzywolek, A. & Plewa, T. 2011, A&A, 529, A156
  * Payez et al. (2015) Payez, A., Evoli, C., Fischer, T., et al. 2015, J. Cosmology Astropart. Phys, 2015, 006
  * Primakoff (1951) Primakoff, H. 1951, Physical Review, 81, 899
  * Raffelt & Stodolsky (1988) Raffelt, G. & Stodolsky, L. 1988, Phys. Rev. D, 37, 1237
  * Raffelt (1996) Raffelt, G. G. 1996, Stars as laboratories for fundamental physics : the astrophysics of neutrinos, axions, and other weakly interacting particles / Georg G. Raffelt. Chicago : University of Chicago Press, 1996.
  * Schmidt et al. (2003) Schmidt, G. D., Harris, H. C., Liebert, J., et al. 2003, ApJ, 595, 1101
  * Schönfelder (2004) Schönfelder, V. 2004, New Astronomy Reviews, 48, 193
  * Seitenzahl et al. (2009) Seitenzahl, I. R., Townsley, D. M., Peng, F., et al. 2009, Atomic Data and Nuclear Data Tables, 95, 96
  * Shen et al. (2018) Shen, K. J., Kasen, D., Miles, B. J., et al. 2018, ApJ, 854, 52
  * Straniero et al. (2019) Straniero, O., Dominguez, I., Piersanti, L., et al. 2019, ApJ, 881, 158
  * Svrcek & Witten (2006) Svrcek, P. & Witten, E. 2006, Journal of High Energy Physics, 2006, 051
  * Tanimori et al. (2015) Tanimori, T., Kubo, H., Takada, A., et al. 2015, ApJ, 810, 28
  * Weinberg (1978) Weinberg, S. 1978, Phys. Rev. Lett., 40, 223
  * Wilczek (1978) Wilczek, F. 1978, Phys. Rev. Lett., 40, 279
  * Woosley & Kasen (2011) Woosley, S. E. & Kasen, D. 2011, ApJ, 734, 38
  * Wright et al. (2016) Wright, W. P., Nagaraj, G., Kneller, J. P., et al. 2016, Phys. Rev. D, 94, 025026