# An Inherent Trade-Off in Noisy Neural Communication with Rank-Order Coding

Ibrahim Alsolami and Tomoki Fukai
Okinawa Institute of Science and Technology (OIST)
Onna, Okinawa, Japan 904-0495

###### Abstract

Rank-order coding, a form of temporal coding, has emerged as a promising
scheme to explain the rapid ability of the mammalian brain. Owing to its speed
as well as efficiency, rank-order coding is increasingly gaining interest in
diverse research areas beyond neuroscience. However, much uncertainty still
exists about the performance of rank-order coding under noise. Herein we show
what information rates are fundamentally possible and what trade-offs are at
stake. An unexpected finding in this paper is the emergence of a special class
of errors that, in a regime, increase with less noise.

##  I Introduction

Currently, there is a growing interest in spiking neural networks in physics
and engineering [1, 2]. Mainly because of their potential to improve power
efficiency, which, to a great extent, depends on the coding scheme employed.
For spiking neural networks, a variety of coding schemes are available. Chief
among them is rank-order coding; this coding scheme offers a fundamentally
different approach for neural information transmission. In rank-order coding,
information is encoded in the order of neural spikes; utilizing this degree of
freedom can boost communication speeds and efficiency substantially.

Rank-order coding has been proposed as a faster alternative to the traditional
rate coding scheme [3]. While rate coding is the most widely accepted coding
scheme in neuroscience, it has a subtle problem: some experimental
observations are hard to reconcile with it because it is slow. For example,
primates can respond selectively to the presentation of 3D objects as quickly
as $100$–$150$ ms after the onset of a stimulus [4, 5]. This response is too
fast to be explained with rate coding as it needs, for a reasonable degree of
accuracy, to accumulate spikes over periods much longer than $150$ ms.
Similarly, humans who were asked to determine whether a briefly flashed
picture ($20$ ms) of a natural scene contains a certain category, such as an
animal, could accurately ($\sim 89$–$98\%$) detect whether a category is
present or not within a few to several hundreds of milliseconds [6]. Rate
coding seems to hardly harmonize with such rapid processing speed.

One could, of course, argue that rate coding across a reasonably large number
of neurons could provide speed (in terms of bits/sec). In such an approach, we
would have $n$ neurons firing in parallel, and one counts the number of spikes
generated within a relatively short time window—compared to a long time window
if we had a few neurons. For instance, we would need a longer time window to
accumulate $n$ spikes from a single neuron than $n$ neurons firing in
parallel. This population rate coding scheme can certainly provide speed, but
one must contend with the fact that such an approach is inefficient in terms
of bits/neurons [7]. In fact, the efficiency of this approach is upper-bounded
by $\frac{\log_{2}(n+1)}{n}$ (bits/neuron)—that is, the more neurons, the less
the efficiency. Moreover, the rapid speed of visual processing is likely to be
accomplished with very few spikes [5], and this would require a neural coding
scheme whereby a few neurons can communicate efficiently.

Rank-order coding can offer both speed and efficiency. With $n$ neurons, rank-
order coding can encode $\log_{2}n!$ bits per transmission compared to
$\log_{2}(n+1)$ bits per transmission for rate coding. The encoding ability of
rank-order coding is vast. Take, for instance, a setting with $10$ neurons.
With these neurons, rank-order coding can in principle form $10!=3,628,800$
symbols, i.e., firing orders of neurons (Fig. 1). As $n$ increases, the
encoding ability of rank-order coding rapidly accelerates. This vast amount of
information available in the arrival order of spikes is often forgotten, and
studying neural codes that utilize such arrival order could provide clues on
how neurons can transmit information rapidly and efficiently across brain
regions.

Converging evidence suggests that the relative timing, or rank, of neuronal
firing plays an important part in encoding information. In retinal ganglion
cells of salamanders, the rapid transmission of visual scenes is likely
accomplished by encoding information in the relative timing of spikes [8]. It
was shown later in a population of retinal ganglion cells of mice that the
content of a visual stimulus could be accurately inferred from the wave of the
first stimulus-evoked spikes, indicating the importance of the relative timing
of spikes to encode sensory information [9]. Analysis of odor-evoked responses
of olfactory neurons of Xenopus laevis (African clawed frogs) demonstrated
that the rank of spike latencies is a reliable predictor of odor identity
[10].

In addition to its biological applicability, rank-order coding is gaining
attention in the field of artificial spiking neural networks, which are
becoming popular as they hold great potential in energy-efficient computing
[11]. It was shown that spiking neural networks with rank-order coding can
achieve a high image-classification accuracy with a relatively small number of
spikes in multilayer feedforward networks [12] as well as in recurrent
networks [13]. Rank-order coding is also finding favor in hardware
implementations of spiking neuromorphic processors; it was demonstrated that
rank-order coding can enhance power efficiency [14] and provide a favorable
trade-off between energy consumption and classification accuracy [15].
Recently, it was shown that rank-order coding successfully reduced the on-chip
inference latency in neuromorphic devices [16].

Despite increasing interest in the use of rank-order coding, far too little
attention has been paid to the effect of noise on its performance. Our goal
here is to analytically understand the impact of noise on the performance of
rank-order coding, as noise is unavoidable in any physical system. In rank-
order coding, noise can cause spikes to be swapped with each other, giving
rise to errors. Herein we study how well rank-order coding performs under
noise in terms of information rate (bits/sec) and communication efficiency
(bits/neuron). Contrary to intuition, reducing noise does not necessarily
reduce all types of errors. Moreover, we show that information rate and
communication efficiency cannot be simultaneously maximized due to an
intrinsic trade-off between them.

##  II Methods

We consider a noise model in which spike times of presynaptic neurons exhibit
random delays characterized by an exponential probability density (Fig. 1, Eq.
8). Such noise may arise when neurons do not respond instantly to a stimulus,
at an expected time, but rather with a random delay (Fig. 1). In rank-order
coding, neurons are intensity-to-delay converters: the higher the activation,
the earlier a neuron fires. Traditional integrate-and-fire models have this
intensity-to-delay property: the higher the membrane potential is, the earlier
a neuron will fire. Without loss of generality, we hypothesize that a
postsynaptic neuron responds selectively and reliably to a particular order of
presynaptic spikes. A feed-forward shunting inhibition circuit was suggested
as the underlying mechanism of this precise decoding of temporal patterns [7].
However, exploring the detailed decoding mechanisms is beyond the scope of
this study.

Figure 1: Rank-order coding with temporal noise (random delay). Here $\alpha$
is the spacing between successive spikes before noise is introduced. In this
illustration, the magnitude of a synaptic weight is represented by the size of
the depicted circles. Here postsynaptic neurons integrate-and-fire and are
progressively desensitized by shunting inhibition circuits (red). With
shunting inhibition, the sensitivity of a neuron progressively decreases as
$\beta^{k}$, where $k$ is the arrival order of a spike, and $\beta$ is a
constant that takes values in the range $0<\beta<1$. A postsynaptic neuron is
maximally activated if spikes arrive in the order of its synaptic weights. By
setting the firing threshold to this maximum excitation/activation level, a
postsynaptic neuron becomes selective to a particular temporal pattern. Due to
noise (random delay), an intended spike sequence can be erroneously received.
For instance, in this illustration, noise can cause the sequence ABC ($x=0$)
to be erroneously received as CBA ($y=5$), which impairs both the
communication rate (bits/sec) and efficiency (bits/neuron).

The channel capacity enables us to compute the maximum amount of information
postsynaptic neurons can receive and is defined as [17, 18]

$C=\max_{p(x)}I(X;Y)~{}~{}~{}~{}~{}\text{(bits/symbol)},$ (1)

where $I(X;Y)$ is the mutual information between random variable $X$ (input
symbol) and $Y$ (output symbol), and is given by

$I(X;Y)=H(Y)-H(Y|X)~{}~{}~{}~{}~{}\text{(bits/symbol)}.$ (2)

Here, $H(Y)$ is the entropy of $Y$, and $H(Y|X)$ is the conditional entropy of
$Y$ given $X$.

In Fig. 1, we have $(n!)^{2}$ possible combinations of input and output
symbols, where a symbol is defined as a particular order of neural spikes
(e.g., the sequence ABC). Here $n$ is the number of presynaptic neurons. The
probability of sending symbol $x$ and, because of noise, receiving symbol $y$
is given by the transition probability $p(y|x)$. The following probability
transition matrix describes such communication channel:

$\begin{split}\bm{p(y|x)}=&\left(\begin{smallmatrix}p(0|0)&p(1|0)&\cdots&p(n!-1|0)\\\
p(0|1)&p(1|1)&\cdots&p(n!-1|1)\\\ \vdots&\vdots&\ddots&\vdots\\\
p(0|n!-1)&p(1|n!-1)&\cdots&p(n!-1|n!-1)\end{smallmatrix}\right)_{n!\times
n!}\\\ \\\ =&\left(\begin{smallmatrix}p_{0}&p_{1}&\cdots&p_{n!-1}\\\
p_{1}&p_{0}&\cdots&p_{n!-2}\\\ \vdots&\vdots&\ddots&\vdots\\\
p_{n!-1}&p_{n!-2}&\cdots&p_{0}\end{smallmatrix}\right)_{n!\times n!}.\\\
\end{split}$ (3)

This communication channel is symmetric because rows of the transition matrix
are permutations of each other, and so are the columns. The capacity of this
channel is achieved by a uniform distribution on the input $X$
[18]$\left(p(x)=\dfrac{1}{n!}\right)$, which results in a uniform distribution
on the output $Y$ $\left(p(y)=\dfrac{1}{n!}\right)$, and is given by

$\begin{split}C=&\max_{p(x)}[H(Y)-H(Y|X)]\\\
=&\log_{2}n!-H(\bm{r})~{}~{}~{}~{}~{}\text{(bits/symbol)},\end{split}$ (4)

where $H(\bm{r})=-\sum\limits_{j=0}^{n!-1}p_{j}\log_{2}p_{j}$ is the entropy
of a row of matrix $\bm{p(y|x)}$.

##  III Results

Transition probabilities. Here we determine the transition probabilities to
find the channel capacity in Eq. (4). These probabilities are the likelihood
that a particular neural spike sequence is received under the perturbation of
noise (random delay, Fig. 1). For instance, for three neurons, $p(CBA|ABC)$ is
the probability that the sequence CBA ($y=5$) is erroneously received due to
noise, given that the original noise-free sequence is ABC ($x=0$). It suffices
to compute the transition probability of any row of $\bm{p(y|x)}$ because the
channel is symmetric [18]. Calculations of the transition probabilities are
straightforward but tedious. Therefore, we only evaluate these probabilities
when the number of presynaptic neurons is relatively small (see APPENDIX A for
derivation).

We can obtain:

$\begin{split}p_{0}=p(AB|AB)=&=1-\dfrac{1}{2}e^{-\lambda\alpha}\\\
p_{1}=p(BA|AB)=&=\dfrac{1}{2}e^{-\lambda\alpha}~{}~{},\end{split}$ (5)

for two presynaptic neurons and

$\begin{split}&~{}~{}~{}~{}~{}p_{0}=p(ABC|ABC)=1-e^{-\lambda\alpha}+\dfrac{1}{6}e^{-3\lambda\alpha}\\\
&~{}~{}~{}~{}~{}p_{1}=p(BAC|ABC)=\dfrac{1}{2}e^{-\lambda\alpha}-\dfrac{1}{2}e^{-2\lambda\alpha}+\dfrac{1}{6}e^{-3\lambda\alpha}\\\
&~{}~{}~{}~{}~{}p_{2}=p(ACB|ABC)=\dfrac{1}{2}e^{-\lambda\alpha}-\dfrac{1}{3}e^{-3\lambda\alpha}\\\
&~{}~{}~{}~{}~{}p_{3}=p(CAB|ABC)=\dfrac{1}{6}e^{-3\lambda\alpha}\\\
&~{}~{}~{}~{}~{}p_{4}=p(BCA|ABC)=\dfrac{1}{2}e^{-2\lambda\alpha}-\dfrac{1}{3}e^{-3\lambda\alpha}\\\
&~{}~{}~{}~{}~{}p_{5}=p(CBA|ABC)=\dfrac{1}{6}e^{-3\lambda\alpha}~{}~{},\end{split}$
(6)

for three presynaptic neurons. In the above expressions, $\alpha$ is the
spacing between successive spikes before noise is introduced, and $\lambda$ is
the rate parameter of the exponential distribution of the noise. Results with
four presynaptic neurons are shown in APPENDIX A.

As expected, when $\lambda\alpha$ increases, error probabilities decrease
(Fig. 2(a)). There is a notable exception, however. In the range
$0\leq\lambda\alpha\leq\ln\sqrt{2}$, the error probability $p(ACB|ABC)$
increases. This can be viewed in two different ways: 1) For a fixed value of
$\alpha$, as the noise decreases (that is, $\lambda$ increases), the
probability of this type of error increases (Fig. 2(b)). 2) For a fixed value
of $\lambda$, as the spacing between spikes ($\alpha$) increases the error
probability increases as well. A similar phenomenon is also observed when we
have four ($n=4$) presynaptic neurons (see Fig. 5). Namely, the following
error probabilities:

$\displaystyle p(ABDC|ABCD),~{}p(ACBD|ABCD),$ $\displaystyle
p(ACDB|ABCD),~{}p(ADBC|ABCD),$ $\displaystyle p(ADCB|ABCD),~{}p(BACD|ABCD),$
$\displaystyle p(BADC|ABCD),~{}p(BCAD|ABCD),$ $\displaystyle
p(BCDA|ABCD),~{}p(BDAC|ABCD),$

and $p(BDCA|ABCD)$. These error probabilities momentarily increase with less
noise, which is counter-intuitive: errors typically decrease with less
noise—not the opposite. Throughout this study, we shall refer to this class of
probabilities as atypical probabilities. This type of error is not limited to
exponential noise; it can also be observed, for example, with Gaussian noise
(see Figs. 8 and 9).

Why do errors increase when we have less noise? The emergence of atypical
probabilities can be explained as follows. Let the probability $P(ACB|ABC)$
serve as an example (Fig. 2). Moreover, let random variables $Z_{i}$ ($i=$1,
2, 3) represent a spike’s latency after the perturbation of noise; here $i$
denotes the index of the $i^{\text{th}}$ presynaptic neuron (See APPENDIX for
notation details). For the event $Z_{1}<Z_{3}<Z_{2}$, or equivalently the
sequence $ACB$, to occur, the following two conditions should be
simultaneously satisfied: (i) Random variable $Z_{2}$ needs to be the largest
value and (ii) random variable $Z_{1}$ needs to be the smallest value.

  1. i

The probability of $Z_{2}$ being the largest value (i.e., $Z_{2}>Z_{3}$)
decreases with $\lambda$ (the larger the value of $\lambda$, the less the
noise) because the amount of overlap between the distributions of $Z_{2}$ and
$Z_{3}$ decreases as $\lambda$ increases.

  2. ii

In contrast, the probability of $Z_{1}$ being the smallest value increases
with $\lambda$ in the interval $(0,2\alpha)$ for the following reason. The
event $Z_{2}>Z_{3}$ implies that $Z_{2}>2\alpha$ (Fig. 2(c) and Eq. 8). Thus,
when $Z_{2}>Z_{3}$, more space (from $\alpha$ to $2\alpha$) for $Z_{1}$ has
been made to take the position of the smallest value in the interval
$(0,2\alpha)$, thereby increasing the likelihood of the neural order $ACB$
($Z_{1}<Z_{3}<Z_{2}$).

Factor (i) causes $P(ACB|ABC)$ to decrease, whereas factor (ii) causes
$P(ACB|ABC)$ to increase. The net effect of factors (i) and (ii) is
$P(Z_{1}<2\alpha<Z_{3}<Z_{2})$, which is a concave function. This component
brings about a rare regime in which errors increase with $\lambda$ (or,
equivalently, with $\alpha$). Mathematically, the probability
$P(Z_{1}<2\alpha<Z_{3}<Z_{2})$ can be obtained by splitting the integration
region of $P(ACB|ABC)$ into two parts:

$\displaystyle
P(ACB|ABC)=P(Z_{1}<Z_{3}<Z_{2}~{}|~{}ABC)=\int_{2\alpha}^{\infty}\int_{2\alpha}^{z_{2}}\int_{0}^{z_{3}}f_{1}(z_{1})f_{3}(z_{3})f_{2}(z_{2})~{}\,dz_{1}\,dz_{3}\,dz_{2}$
(7)
$\displaystyle=\underbrace{\left(\int_{2\alpha}^{\infty}\int_{2\alpha}^{z_{2}}\int_{0}^{2\alpha}f_{1}(z_{1})f_{3}(z_{3})f_{2}(z_{2})~{}\,dz_{1}\,dz_{3}\,dz_{2}\right)}_{\text{Concave:
$P(Z_{1}<2\alpha<Z_{3}<Z_{2})$}}+\underbrace{\left(\int_{2\alpha}^{\infty}\int_{2\alpha}^{z_{2}}\int_{2\alpha}^{z_{3}}f_{1}(z_{1})f_{3}(z_{3})f_{2}(z_{2})~{}\,dz_{1}\,dz_{3}\,dz_{2}\right)}_{\text{Convex:
$P(2\alpha<Z_{1}<Z_{3}<Z_{2})$}}$
$\displaystyle=\underbrace{\left(\underbrace{\dfrac{1}{2}e^{-\lambda\alpha}}_{\text{Factor
(i):
$P(2\alpha<Z_{3}<Z_{2})$}}\times~{}~{}~{}\underbrace{\left(1-e^{-2\lambda\alpha}\right)}_{\text{Factor
(ii): $P(Z_{1}<2\alpha)$}}\right)}_{\text{Concave:
$P(Z_{1}<2\alpha<Z_{3}<Z_{2})=P(Z_{1}<2\alpha)~{}\times~{}P(2\alpha<Z_{3}<Z_{2})$}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\underbrace{\left(\dfrac{1}{6}e^{-3\lambda\alpha}\right)}_{\text{Convex:
$P(2\alpha<Z_{1}<Z_{3}<Z_{2})$}}=\dfrac{1}{2}e^{-\lambda\alpha}-\dfrac{1}{3}e^{-3\lambda\alpha}~{}~{}.$

Here $f_{i}(z)$ is the probability density function (pdf) of the
$i^{\text{th}}$ presynaptic neuron, $i\in\\{1,2,\dots,n\\}$. This pdf
describes the likelihood of observing a randomly delayed spike, by noise, at
time $z$ and is given by

$f_{i}(z)=\begin{cases}\lambda~{}e^{-\lambda\bigl{(}z-(i-1)\alpha\bigr{)}}&,~{}\text{if}~{}z\geq(i-1)\alpha\\\
0&,~{}\text{otherwise}~{}~{}~{}.\end{cases}$ (8)

(a)

(b)

(c)

Figure 2: Transition probabilities and spikes of rank-order coding with
temporal noise. (a) Transition probabilities for three presynaptic neurons
($n=3$). (b) Atypical probability $p(ACB|ABC)$; this probability reflects the
likelihood that noise causes the sequence $ABC$ to be erroneously received as
$ACB$. In the range $0<\lambda<\ln\sqrt{2}$, the probability of an error
increases with $\lambda$ (the higher the value of $\lambda$, the less the
noise). Here $\alpha$ is arbitrarily set to 1 (sec), and the number of samples
per point used in the simulation is $10^{9}$. A similar phenomenon, where
errors momentarily increase, can also be observed by fixing $\lambda$ and
varying $\alpha$; that is, errors increase as the spacing between neural
spikes, $\alpha$, increases. (c) Spikes of rank-order coding with temporal
noise (random delay). When $Z_{2}>Z_{3}$, more space for $Z_{1}$ has been made
to take the position of the smallest value in the interval $(0,2\alpha)$; this
causes the probability of the event $ACB$ ($Z_{1}<Z_{3}<Z_{2}$) to increase.

Communication efficiency and information rate. Figure 3(a) shows the
performance of rank-order coding in terms of communication efficiency, which
is defined as

$\gamma=\frac{C}{n}~{}~{}~{}~{}~{}~{}~{}\text{(bits/neuron)},$ (9)

where $C$ can be calculated by using Eq. 4. The communication efficiency
increases monotonically with $\lambda\alpha$. This increase in efficiency
eventually plateaus and is asymptotically bounded by
$\gamma^{*}=\lim_{\lambda\alpha\to\infty}\frac{C}{n}=\frac{\log_{2}(n!)}{n}$
(bits/neuron). Moreover, the higher the value of $n$, the more efficient the
communication is.

We can further evaluate the performance of rank-order coding in terms of
information rate, which is defined as

$R=\frac{C}{\bar{T}}~{}~{}~{}~{}~{}~{}~{}\text{(bits/sec)}.$ (10)

In the absence of noise, the average symbol duration $\bar{T}$ (that is, the
time difference between the first and last spikes of a rank-order coding
symbol) is $(n-1)\alpha$. However, with noise, the average symbol duration
increases and is given by

$\bar{T}=\alpha+\frac{1}{\lambda}e^{-\lambda\alpha}~{}~{}~{}\text{(sec/symbol)}$
(11)

for two presynaptic neurons ($n=2$),

$\bar{T}=2\alpha+\frac{1}{\lambda}e^{-\lambda\alpha}+\frac{1}{2\lambda}e^{-2\lambda\alpha}~{}~{}~{}\text{(sec/symbol)}$
(12)

for three presynaptic neurons ($n=3$), and

$\displaystyle\bar{T}=3\alpha$
$\displaystyle+\frac{1}{\lambda}e^{-\lambda\alpha}+\frac{1}{2\lambda}e^{-2\lambda\alpha}+\frac{1}{2\lambda}e^{-3\lambda\alpha}-\frac{1}{6\lambda}e^{-4\lambda\alpha}$
$\displaystyle-\frac{1}{6\lambda}e^{-5\lambda\alpha}+\frac{1}{6\lambda}e^{-6\lambda\alpha}~{}~{}~{}\text{(sec/symbol)}$
(13)

for four presynaptic neurons ($n=4$) (see APPENDIX C for derivation). In Fig.
3(b), we display the (scaled) information rate as a function of
$\lambda\alpha$. The information rate is a non-monotonic function of
$\lambda\alpha$ and increases with $n$. Moreover, there is an optimal
operating point at which the information rate is maximized; beyond this
critical point, the information rate rapidly diminishes.

In a noisy environment, there exists an inherent trade-off between the
communication efficiency of rank-order coding and its information rate. The
communication efficiency continuously increases with $\lambda\alpha$ (Fig.
3(a)), but this gain in efficiency comes at the cost of a loss in the
information rate once $\lambda\alpha$ is beyond a critical value (Fig. 3(b)).
A range of trade-offs is shown in Fig. 3(c), in which the value of
$\lambda\alpha$ is varied and the pair $(\gamma,\frac{R}{\lambda})$ is
displayed. The resultant curves represent upper-bounds of achievable
information rates and communication efficiencies. Parameter $\alpha$ provides
a means to control the trade-off between information rate and efficiency.

(a)

(b)

(c)

Figure 3: Performance of rank-order coding. (a) Communication efficiency. (b)
Information rates. Here we plot the scaled version ($\frac{R}{\lambda}$) of
$R$ rather than $R$ as it eliminates the need to display $R$ for various
combinations of $\lambda$ and $\alpha$. (c) The trade-off between information
rates and efficiency.

##  IV Discussion

This paper set out to study the impact of noise on the performance of rank-
order coding. Rank-order coding is advantageous as it utilizes the order of
neural spikes, enabling it to boost communication speeds. A disadvantage,
however, is that it is susceptible to temporal noise, which can swap
presynaptic spikes with each other causing errors at postsynaptic neurons. As
such, we considered noise in the form of a random delay to gain insights into
the performance of rank-order coding in terms of information rate and
communication efficiency.

In noisy environments, the information rate and communication efficiency
depend at least on three factors: the spacing between spikes $~{}\alpha$, the
rate parameter $~{}\lambda$, and the number of presynaptic neurons $~{}n$. The
higher the value of $\lambda\alpha$, the more efficient the communication is.
However, increasing $\lambda\alpha$ beyond an optimal operating point has the
adverse effect of reducing the information rate. Additionally, we found a
class of error probabilities that increase with less noise. This result is
counter-intuitive because errors commonly decrease with less noise—not the
opposite. The presence of such error probabilities raises a need for special
care in designing error correction schemes for neuromorphic devices that
employ rank-order coding.

We revealed that rank-order coding has an inherent trade-off between
information rate and communication efficiency. This result could provide
insights to better understand what trade-offs neurons in different brain
regions make (under the rank-order coding hypothesis) between the conflicting
needs to be fast and, at the same time, efficient. For example, it is likely
that neurons at the early stage of the visual-processing pathway (e.g., the
retina) prioritize speed over efficiency. However, efficiency may be favored
over speed at later stages as information is likely at/near its final
(decision-making) destination. The trade-off result also offers a realistic
picture of neuromorphic computing with rank-order coding: information rate and
communication efficiency cannot be simultaneously maximized—a compromising
trade-off between them needs to be made (Fig. 3(c)).

In the present study, we assumed that postsynaptic neurons respond selectively
to a particular order of spikes (temporal pattern). Studies have shown that
cortical neurons exhibit such selectivity to temporal input sequences [19].
Various biological mechanisms of temporal pattern detection have been proposed
(e.g., [20]). A feed-forward shunting inhibition circuit, which progressively
desensitizes a postsynaptic neuron as spikes arrive (see Fig. 1), may
accomplish selectivity to a particular temporal pattern [7]. In such a
setting, a postsynaptic neuron would be maximally activated (and fire only) if
spikes arrive in the order of its synaptic weights. A small portion of
extremely strong synapses observed in log-normally distributed synaptic
weights [21] may enhance this progressive desensitization.

Generalization of our results to an arbitrary number of neurons is a promising
research extension. Another avenue of extension would be to assess the
performance of rank-order coding when not all spikes arrive at postsynaptic
neurons; some presynaptic neurons, for instance, may misfire. The combination
of such noise with temporal noise can significantly affect performance.
Nonetheless, hybrid noise would only affect quantitative results but would not
change the qualitative results of this study.

To conclude, rank-order coding can provide speed and efficiency, but noise
imposes a trade-off between them. The results of this study offer a novel
insight into the performance of rank-order coding.

## Acknowledgment

The authors would like to thank Balashwethan Chockalingam and Thomas Burns for
their suggestions and comments. We are grateful for the help and support
provided by the Scientific Computing and Data Analysis section of the Research
Support Division at OIST. T.F. acknowledges support from KAKENHI grants
JP19H04994 and JP18H05213.

## References

  * [1] D. Marković, A. Mizrahi, D. Querlioz, and J. Grollier, “Physics for neuromorphic computing,” _Nature Reviews Physics_ , vol. 2, no. 9, pp. 499–510, 2020.
  * [2] G. V. Vinci, R. Benzi, and M. Mattia, “Self-consistent stochastic dynamics for finite-size networks of spiking neurons,” _Physical Review Letters_ , vol. 130, no. 9, p. 097402, 2023.
  * [3] R. VanRullen, R. Guyonneau, and S. J. Thorpe, “Spike times make sense,” _Trends in Neurosciences_ , vol. 28, no. 1, pp. 1–4, 2005.
  * [4] S. Thorpe, A. Delorme, and R. Van Rullen, “Spike-based strategies for rapid processing,” _Neural Networks_ , vol. 14, no. 6-7, 2001.
  * [5] S. J. Thorpe and M. Imbert, “Biological constraints on connectionist modelling,” _Connectionism in Perspective_ , 1989.
  * [6] S. Thorpe, D. Fize, and C. Marlot, “Speed of processing in the human visual system,” _Nature_ , vol. 381, no. 6582, pp. 520–522, 1996.
  * [7] L. Bonilla, J. Gautrais, S. Thorpe, and T. Masquelier, “Analyzing time-to-first-spike coding schemes,” _Frontiers in Neuroscience_ , vol. 16, 2022.
  * [8] T. Gollisch and M. Meister, “Rapid neural coding in the retina with relative spike latencies,” _Science_ , vol. 319, no. 5866, 2008.
  * [9] G. Portelli, J. M. Barrett, G. Hilgen, T. Masquelier, A. Maccione, S. Di Marco, L. Berdondini, P. Kornprobst, and E. Sernagor, “Rank order coding: a retinal information decoding strategy revealed by large-scale multielectrode array retinal recordings,” _eNeuro_ , vol. 3, no. 3, 2016.
  * [10] S. Junek, E. Kludt, F. Wolf, and D. Schild, “Olfactory coding with patterns of response latencies,” _Neuron_ , vol. 67, no. 5, 2010.
  * [11] E. Lemaire, L. Cordone, A. Castagnetti, P.-E. Novac, J. Courtois, and B. Miramond, “An analytical estimation of spiking neural networks energy efficiency,” in _International Conference on Neural Information Processing_. Springer, 2022, pp. 574–587.
  * [12] S. R. Kheradpisheh, M. Ganjtabesh, S. J. Thorpe, and T. Masquelier, “STDP-based spiking deep convolutional neural networks for object recognition,” _Neural Networks_ , vol. 99, pp. 56–67, 2018.
  * [13] Y. Yan, H. Chu, X. Chen, Y. Jin, Y. Huan, L. Zheng, and Z. Zou, “Graph-based spatio-temporal backpropagation for training spiking neural networks,” in _2021 IEEE 3rd International Conference on Artificial Intelligence Circuits and Systems (AICAS)_. IEEE, 2021, pp. 1–4.
  * [14] C. Frenkel, M. Lefebvre, J.-D. Legat, and D. Bol, “A 0.086-$\text{mm}^{2}$ 12.7-pJ/SOP 64k-Synapse 256-neuron online-learning digital spiking neuromorphic processor in 28-nm CMOS,” _IEEE Transactions on Biomedical Circuits and Systems_ , vol. 13, no. 1, pp. 145–158, 2018.
  * [15] C. Frenkel, J.-D. Legat, and D. Bol, “MorphIC: A 65-nm 738k-Synapse/$\text{mm}^{2}$ quad-core binary-weight digital neuromorphic processor with stochastic spike-driven online learning,” _IEEE Transactions on Biomedical Circuits and Systems_ , vol. 13, no. 5, pp. 999–1010, 2019.
  * [16] Y. Wu, R. Zhao, J. Zhu, F. Chen, M. Xu, G. Li, S. Song, L. Deng, G. Wang, H. Zheng _et al._ , “Brain-inspired global-local learning incorporated with neuromorphic computing,” _Nature Communications_ , vol. 13, no. 1, p. 65, 2022.
  * [17] C. E. Shannon, “A mathematical theory of communication,” _The Bell System Technical Journal_ , vol. 27, no. 3, pp. 379–423, 1948.
  * [18] T. Cover and J. Thomas, _Elements of Information Theory_. Wiley, 2012.
  * [19] T. Branco, B. A. Clark, and M. Häusser, “Dendritic discrimination of temporal input sequences in cortical neurons,” _Science_ , vol. 329, no. 5999, pp. 1671–1675, 2010.
  * [20] G. R and S. H, “The tempotron: a neuron that learns spike timing-based decisions,” _Nature Neuroscience_ , vol. 9, no. 5866, pp. 420–428, 2006\.
  * [21] B. G and M. K, “The log-dynamic brain: how skewed distributions affect network operations,” _Nature Reviews Neuroscience_ , vol. 15, no. 4, pp. 264–278, 2014.

##  Appendix A Transition Probabilities of Rank-Order Coding with Exponential
Noise (Theory)

Two Neurons. For two presynaptic neurons, the transition probabilities of the
first row of Eq. (3) are

$\begin{split}p_{0}=p(AB|AB)=&\int_{\alpha}^{\infty}\int_{0}^{z_{2}}f_{1}(z_{1})f_{2}(z_{2})\,dz_{1}\,dz_{2}=1-\dfrac{1}{2}e^{-\lambda\alpha}\\\
p_{1}=p(BA|AB)=&\int_{\alpha}^{\infty}\int_{\alpha}^{z_{1}}f_{2}(z_{2})f_{1}(z_{1})\,dz_{2}\,dz_{1}=\dfrac{1}{2}e^{-\lambda\alpha}~{}~{}.\end{split}$

Let $Z_{i}$ be a random variable and $f_{i}(z)$ its probability density
function (pdf); here $f_{i}(z)$ is the pdf of the $i^{\text{th}}$ presynaptic
neuron, $i\in\\{1,2,\dots,n\\}$. This pdf describes the likelihood of
observing a randomly delayed spike, by noise, at time $z$ and is given by

$f_{i}(z)=\begin{cases}\lambda~{}e^{-\lambda\bigl{(}z-(i-1)\alpha\bigr{)}}&,~{}\text{if}~{}z\geq(i-1)\alpha\\\
0&,~{}\text{otherwise}\end{cases}$ (14)

here $\alpha$ is the spacing between successive spikes before noise is
introduced, and $\lambda$ is the rate parameter of the exponential
distribution of the noise. Herein random variables, $Z_{i}$, are independent.

Fig. 4 shows the transition probabilities for two presynaptic neurons ($n=2$).

Figure 4: Transition probabilities for two presynaptic neurons ($n=2$) with
exponential noise. Here $p(y|x)$ is the probability of sending sequence $x$
and receiving sequence $y$.

Three Neurons. For three presynaptic neurons, the transition probabilities of
the first row of Eq. (3) are

$\begin{split}p_{0}=p(ABC|ABC)=&\int_{2\alpha}^{\infty}\int_{\alpha}^{z_{3}}\int_{0}^{z_{2}}f_{1}(z_{1})f_{2}(z_{2})f_{3}(z_{3})~{}\,dz_{1}\,dz_{2}\,dz_{3}=1-e^{-\lambda\alpha}+\dfrac{1}{6}e^{-3\lambda\alpha}\\\
p_{1}=p(BAC|ABC)=&\int_{2\alpha}^{\infty}\int_{\alpha}^{z_{3}}\int_{\alpha}^{z_{1}}f_{2}(z_{2})f_{1}(z_{1})f_{3}(z_{3})~{}\,dz_{2}\,dz_{1}\,d_{3}=\dfrac{1}{2}e^{-\lambda\alpha}-\dfrac{1}{2}e^{-2\lambda\alpha}+\dfrac{1}{6}e^{-3\lambda\alpha}\\\
p_{2}=p(ACB|ABC)=&\int_{2\alpha}^{\infty}\int_{2\alpha}^{z_{2}}\int_{0}^{z_{3}}f_{1}(z_{1})f_{3}(z_{3})f_{2}(z_{2})~{}\,dz_{1}\,dz_{3}\,dz_{2}=\dfrac{1}{2}e^{-\lambda\alpha}-\dfrac{1}{3}e^{-3\lambda\alpha}\\\
p_{3}=p(CAB|ABC)=&\int_{2\alpha}^{\infty}\int_{2\alpha}^{z_{2}}\int_{2\alpha}^{z_{1}}f_{3}(z_{3})f_{1}(z_{1})f_{2}(z_{2})~{}\,dz_{3}\,dz_{1}\,dz_{2}=\dfrac{1}{6}e^{-3\lambda\alpha}\\\
p_{4}=p(BCA|ABC)=&\int_{2\alpha}^{\infty}\int_{2\alpha}^{z_{1}}\int_{\alpha}^{z_{3}}f_{2}(z_{2})f_{3}(z_{3})f_{1}(z_{1})~{}\,dz_{2}\,dz_{3}\,dz_{1}=\dfrac{1}{2}e^{-2\lambda\alpha}-\dfrac{1}{3}e^{-3\lambda\alpha}\\\
p_{5}=p(CBA|ABC)=&\int_{2\alpha}^{\infty}\int_{2\alpha}^{z_{1}}\int_{2\alpha}^{z_{2}}f_{3}(z_{3})f_{2}(z_{2})f_{1}(z_{1})~{}\,dz_{3}\,dz_{2}\,dz_{1}=\dfrac{1}{6}e^{-3\lambda\alpha}~{}~{}~{}.\\\
\end{split}$

Fig. 2(a) shows the transition probabilities for three presynaptic neurons
($n=3$).

Four Neurons. For four presynaptic neurons, the transition probabilities of
the first row of Eq. (3) are

$\begin{split}p_{0}=p(ABCD|ABCD)=&\int_{3\alpha}^{\infty}\int_{2\alpha}^{z_{4}}\int_{\alpha}^{z_{3}}\int_{0}^{z_{2}}f_{1}(z_{1})f_{2}(z_{2})f_{3}(z_{3})f_{4}(z_{4})~{}\,dz_{1}\,dz_{2}\,dz_{3}\,dz_{4}=1-\dfrac{3}{2}e^{-\lambda\alpha}+\dfrac{1}{4}e^{-2\lambda\alpha}+\dfrac{1}{3}e^{-3\lambda\alpha}-\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{1}=p(BACD|ABCD)=&\int_{3\alpha}^{\infty}\int_{2\alpha}^{z_{4}}\int_{\alpha}^{z_{3}}\int_{\alpha}^{z_{1}}f_{2}(z_{2})f_{1}(z_{1})f_{3}(z_{3})f_{4}(z_{4})~{}\,dz_{2}\,dz_{1}\,dz_{3}\,dz_{4}=\dfrac{1}{2}e^{-\lambda\alpha}-\dfrac{3}{4}e^{-2\lambda\alpha}+\dfrac{1}{6}e^{-3\lambda\alpha}+\dfrac{1}{6}e^{-4\lambda\alpha}-\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{2}=p(ACBD|ABCD)=&\int_{3\alpha}^{\infty}\int_{2\alpha}^{z_{4}}\int_{2\alpha}^{z_{2}}\int_{0}^{z_{3}}f_{1}(z_{1})f_{3}(z_{3})f_{2}(z_{2})f_{4}(z_{4})~{}\,dz_{1}\,dz_{3}\,dz_{2}\,dz_{4}=\dfrac{1}{2}e^{-\lambda\alpha}-\dfrac{1}{2}e^{-2\lambda\alpha}-\dfrac{1}{6}e^{-3\lambda\alpha}+\dfrac{1}{4}e^{-4\lambda\alpha}-\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{3}=p(CABD|ABCD)=&\int_{3\alpha}^{\infty}\int_{2\alpha}^{z_{4}}\int_{2\alpha}^{z_{2}}\int_{2\alpha}^{z_{1}}f_{3}(z_{3})f_{1}(z_{1})f_{2}(z_{2})f_{4}(z_{4})~{}\,dz_{3}\,dz_{1}\,dz_{2}\,dz_{4}=\dfrac{1}{6}e^{-3\lambda\alpha}-\dfrac{1}{4}e^{-4\lambda\alpha}+\dfrac{1}{6}e^{-5\lambda\alpha}-\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{4}=p(BCAD|ABCD)=&\int_{3\alpha}^{\infty}\int_{2\alpha}^{z_{4}}\int_{2\alpha}^{z_{1}}\int_{\alpha}^{z_{3}}f_{2}(z_{2})f_{3}(z_{3})f_{1}(z_{1})f_{4}(z_{4})~{}\,dz_{2}\,dz_{3}\,dz_{1}\,dz_{4}=\dfrac{1}{2}e^{-2\lambda\alpha}-\dfrac{5}{6}e^{-3\lambda\alpha}+\dfrac{5}{12}e^{-4\lambda\alpha}-\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{5}=p(CBAD|ABCD)=&\int_{3\alpha}^{\infty}\int_{2\alpha}^{z_{4}}\int_{2\alpha}^{z_{1}}\int_{2\alpha}^{z_{2}}f_{3}(z_{3})f_{2}(z_{2})f_{1}(z_{1})f_{4}(z_{4})~{}\,dz_{3}\,dz_{2}\,dz_{1}\,dz_{4}=\dfrac{1}{6}e^{-3\lambda\alpha}-\dfrac{1}{4}e^{-4\lambda\alpha}+\dfrac{1}{6}e^{-5\lambda\alpha}-\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{6}=p(ABDC|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{3}}\int_{\alpha}^{z_{4}}\int_{0}^{z_{2}}f_{1}(z_{1})f_{2}(z_{2})f_{4}(z_{4})f_{3}(z_{3})~{}\,dz_{1}\,dz_{2}\,dz_{4}\,dz_{3}=\dfrac{1}{2}e^{-\lambda\alpha}-\dfrac{1}{4}e^{-2\lambda\alpha}-\dfrac{1}{3}e^{-3\lambda\alpha}+\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{7}=p(BADC|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{3}}\int_{\alpha}^{z_{4}}\int_{\alpha}^{z_{1}}f_{2}(z_{2})f_{1}(z_{1})f_{4}(z_{4})f_{3}(z_{3})~{}\,dz_{2}\,dz_{1}\,dz_{4}\,dz_{3}=\dfrac{1}{4}e^{-2\lambda\alpha}-\dfrac{1}{3}e^{-4\lambda\alpha}+\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{8}=p(ADBC|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{3}}\int_{3\alpha}^{z_{2}}\int_{0}^{z_{4}}f_{1}(z_{1})f_{4}(z_{4})f_{2}(z_{2})f_{3}(z_{3})~{}\,dz_{1}\,dz_{4}\,dz_{2}\,dz_{3}=\dfrac{1}{6}e^{-3\lambda\alpha}-\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{9}=p(DABC|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{3}}\int_{3\alpha}^{z_{2}}\int_{3\alpha}^{z_{1}}f_{4}(z_{4})f_{1}(z_{1})f_{2}(z_{2})f_{3}(z_{3})~{}\,dz_{4}\,dz_{1}\,dz_{2}\,dz_{3}=\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{10}=p(BDAC|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{3}}\int_{3\alpha}^{z_{1}}\int_{\alpha}^{z_{4}}f_{2}(z_{2})f_{4}(z_{4})f_{1}(z_{1})f_{3}(z_{3})~{}\,dz_{2}\,dz_{4}\,dz_{1}\,dz_{3}=\dfrac{1}{6}e^{-4\lambda\alpha}-\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{11}=p(DBAC|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{3}}\int_{3\alpha}^{z_{1}}\int_{3\alpha}^{z_{2}}f_{4}(z_{4})f_{2}(z_{2})f_{1}(z_{1})f_{3}(z_{3})~{}\,dz_{4}\,dz_{2}\,dz_{1}\,dz_{3}=\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{12}=p(ACDB|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{2}}\int_{2\alpha}^{z_{4}}\int_{0}^{z_{3}}f_{1}(z_{1})f_{3}(z_{3})f_{4}(z_{4})f_{2}(z_{2})~{}\,dz_{1}\,dz_{3}\,dz_{4}\,dz_{2}=\dfrac{1}{2}e^{-2\lambda\alpha}-\dfrac{1}{3}e^{-3\lambda\alpha}-\dfrac{1}{4}e^{-4\lambda\alpha}+\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{13}=p(CADB|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{2}}\int_{2\alpha}^{z_{4}}\int_{2\alpha}^{z_{1}}f_{3}(z_{3})f_{1}(z_{1})f_{4}(z_{4})f_{2}(z_{2})~{}\,dz_{3}\,dz_{1}\,dz_{4}\,dz_{2}=\dfrac{1}{4}e^{-4\lambda\alpha}-\dfrac{1}{3}e^{-5\lambda\alpha}+\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{14}=p(ADCB|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{2}}\int_{3\alpha}^{z_{3}}\int_{0}^{z_{4}}f_{1}(z_{1})f_{4}(z_{4})f_{3}(z_{3})f_{2}(z_{2})~{}\,dz_{1}\,dz_{4}\,dz_{3}\,dz_{2}=\dfrac{1}{6}e^{-3\lambda\alpha}-\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{15}=p(DACB|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{2}}\int_{3\alpha}^{z_{3}}\int_{3\alpha}^{z_{1}}f_{4}(z_{4})f_{1}(z_{1})f_{3}(z_{3})f_{2}(z_{2})~{}\,dz_{4}\,dz_{1}\,dz_{3}\,dz_{2}=\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{16}=p(CDAB|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{2}}\int_{3\alpha}^{z_{1}}\int_{2\alpha}^{z_{4}}f_{3}(z_{3})f_{4}(z_{4})f_{1}(z_{1})f_{2}(z_{2})~{}\,dz_{1}\,dz_{3}\,dz_{1}\,dz_{2}=\dfrac{1}{6}e^{-5\lambda\alpha}-\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{17}=p(DCAB|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{2}}\int_{3\alpha}^{z_{1}}\int_{3\alpha}^{z_{3}}f_{4}(z_{4})f_{3}(z_{3})f_{1}(z_{1})f_{2}(z_{2})~{}\,dz_{3}\,dz_{4}\,dz_{1}\,dz_{2}=\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{18}=p(BCDA|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{1}}\int_{2\alpha}^{z_{4}}\int_{\alpha}^{z_{3}}f_{2}(z_{2})f_{3}(z_{3})f_{4}(z_{4})f_{1}(z_{1})~{}\,dz_{2}\,dz_{3}\,dz_{4}\,dz_{1}=\dfrac{1}{2}e^{-3\lambda\alpha}-\dfrac{7}{12}e^{-4\lambda\alpha}+\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{19}=p(CBDA|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{1}}\int_{2\alpha}^{z_{4}}\int_{2\alpha}^{z_{2}}f_{3}(z_{3})f_{2}(z_{2})f_{4}(z_{4})f_{1}(z_{1})~{}\,dz_{3}\,dz_{2}\,dz_{4}\,dz_{1}=\dfrac{1}{4}e^{-4\lambda\alpha}-\dfrac{1}{3}e^{-5\lambda\alpha}+\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{20}=p(BDCA|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{1}}\int_{3\alpha}^{z_{3}}\int_{\alpha}^{z_{4}}f_{2}(z_{2})f_{4}(z_{4})f_{3}(z_{3})f_{1}(z_{1})~{}\,dz_{2}\,dz_{4}\,dz_{3}\,dz_{1}=\dfrac{1}{6}e^{-4\lambda\alpha}-\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{21}=p(DBCA|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{1}}\int_{3\alpha}^{z_{3}}\int_{3\alpha}^{z_{2}}f_{4}(z_{4})f_{2}(z_{2})f_{3}(z_{3})f_{1}(z_{1})~{}\,dz_{4}\,dz_{2}\,dz_{3}\,dz_{1}=\dfrac{1}{24}e^{-6\lambda\alpha}\\\
p_{22}=p(CDBA|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{1}}\int_{3\alpha}^{z_{2}}\int_{2\alpha}^{z_{4}}f_{3}(z_{3})f_{4}(z_{4})f_{2}(z_{2})f_{1}(z_{1})~{}\,dz_{3}\,dz_{4}\,dz_{2}\,dz_{1}=\dfrac{1}{6}e^{-5\lambda\alpha}-\dfrac{1}{8}e^{-6\lambda\alpha}\\\
p_{23}=p(DCBA|ABCD)=&\int_{3\alpha}^{\infty}\int_{3\alpha}^{z_{1}}\int_{3\alpha}^{z_{2}}\int_{3\alpha}^{z_{3}}f_{4}(z_{4})f_{3}(z_{3})f_{2}(z_{2})f_{1}(z_{1})~{}\,dz_{4}\,dz_{3}\,dz_{2}\,dz_{1}=\dfrac{1}{24}e^{-6\lambda\alpha}~{}~{}~{}.\end{split}$

Fig. 5 shows the transition probabilities for four presynaptic neurons
($n=4$).

Figure 5: Transition probabilities for four presynaptic neurons ($n=4$) with
exponential noise. Here $p(y|x)$ is the probability of sending sequence $x$
and receiving sequence $y$.

##  Appendix B Transition Probabilities of Rank-Order Coding with Gaussian
Noise (Simulation)

Consider Fig. 6. In this section, we provide simulation results of transition
probability when the temporal noise is Gaussian,
$Z_{i}\sim\mathcal{N}\left((i-1)\alpha,\,\sigma^{2}\right)\,$. Here random
variable $Z_{i}$ describes the random jitter caused by noise. The variance of
the Gaussian noise is $\sigma^{2}$ and its mean is $(i-1)\alpha$, where index
$i$ is for the $i^{\text{th}}$ number of presynaptic neuron, and $\alpha$ is
the original spacing between successive spikes before noise is introduced.

Figure 6: Spikes of rank-order coding with temporal noise (Guassian). Figure
7: Simulation results of transition probabilities for two presynaptic neurons
($n=2$) with Gaussian noise. Here $p(y|x)$ is the probability of sending
sequence $x$ and receiving sequence $y$. The variance of the Gaussian
distribution, $\sigma^{2}$, is arbitrarily set to 1 (sec2). Additionally,
$\alpha$ is the original spacing between successive spikes before noise is
introduced. The number of samples per point used in the simulation is
$10^{9}$. Figure 8: Simulation results of transition probabilities for three
presynaptic neurons ($n=3$) with Gaussian noise. Here $p(y|x)$ is the
probability of sending sequence $x$ and receiving sequence $y$. The variance
of the Gaussian distribution, $\sigma^{2}$, is arbitrarily set to 1 (sec2).
Additionally, $\alpha$ is the original spacing between successive spikes
before noise is introduced. The number of samples per point used in the
simulation is $10^{9}$. Figure 9: Simulation results of transition
probabilities for four presynaptic neurons ($n=4$) with Gaussian noise. Here
$p(y|x)$ is the probability of sending sequence $x$ and receiving sequence
$y$. The variance of the Gaussian distribution, $\sigma^{2}$, is arbitrarily
set to 1 (sec2). Additionally, $\alpha$ is the original spacing between
successive spikes before noise is introduced. The number of samples per point
used in the simulation is $10^{9}$.

##  Appendix C Average Symbol Duration, $\bar{T}$

Due to noise, the duration between the first ($Z_{(n)}$) and last spike
($Z_{(1)}$) of a rank-order coding symbol is a random variable
$T=Z_{(n)}-Z_{(1)}$. Our goal in this section is to find the average value of
$T$ for the set $\\{2,3,4\\}$ of presynaptic neurons. We first start by
finding the pdfs of $Z_{(1)}$ and $Z_{(n)}$, and then their mean values.

Let $Z_{i}$ be a random variable and $f_{Z_{i}}(z)$ its pdf; here index $i$
indicates the $i^{\text{th}}$ presynaptic neuron. This pdf describes the
likelihood of observing a randomly jittered spike, by noise, at time $z$ and
is given by

$f_{Z_{i}}(z)=\begin{cases}\lambda~{}e^{-\lambda\bigl{(}z-(i-1)\alpha\bigr{)}}&,~{}\text{if}~{}z\geq(i-1)\alpha\\\
0&,~{}\text{otherwise}\end{cases}$ (15)

here $\alpha$ is the spacing between successive spikes before noise is
introduced, and $\lambda$ is the rate parameter of the exponential
distribution of the noise. Herein random variables, $Z_{i}$, are independent.

###  C-A Two Presynaptic Neurons

Minimum. Let random variable $Z_{(1)}$ be the first-order statistic defined as

$Z_{(1)}=\min\\{Z_{1},Z_{2}\\}~{}~{}.$ (16)

The pdf of $Z_{(1)}$ has three cases:

Case 1: $z<0$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2}\\}=\varnothing~{}\\\\[7.11317pt]
\end{split}$ (17)

because if $z<0$, then $f_{Z_{1}}(z)=f_{Z_{2}}(z)=0$, (Eq. 15). Here,
$\varnothing$ denotes the empty set.
Therefore,

$f_{Z_{(1)}}(z)=0~{}.$ (18)

Case 2: $0\leq z<\alpha$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2}\\}=Z_{1}~{}~{}\\\\[7.11317pt]
\end{split}$ (19)

because in the range $0\leq z<\alpha$, $f_{Z_{2}}(z)=0$, (Eq. 15). Therefore,

$f_{Z_{(1)}}(z)=\lambda e^{-\lambda{z}}~{}.$ (20)

Case 3: $\alpha\leq z$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2}\\}.\end{split}$ (21)

The cdf of the minimum is

$\begin{split}F_{Z_{(1)}}(z)=&P\left(Z_{(1)}\leq z\right)\\\
=&1-P\left(Z_{(1)}>z\right)\\\ =&1-P(Z_{1}>z,Z_{2}>z)\\\
=&1-P(Z_{1}>z)(Z_{2}>z)~{}~{}~{}~{}\text{because of independence}\\\
=&1-\Bigl{[}\Bigl{(}1-P(Z_{1}\leq z)\Bigr{)}\Bigl{(}1-P(Z_{2}\leq
z)\Bigr{)}\Bigr{]}\\\
=&1-e^{-2\lambda(z-\frac{1}{2}\alpha)}~{}~{}.\end{split}$ (22)

Accordingly, the pdf of the minimum is

$f_{Z_{(1)}}(z)=\frac{d}{dz}F_{Z_{(1)}}(z)=2\lambda
e^{-2\lambda(z-\frac{1}{2}\alpha)}~{}~{}.$ (23)

Combining Case 1 to 3:

$f_{Z_{(1)}}(z)=\begin{cases}0&,~{}\text{if}~{}z<0\\\\[7.11317pt] \lambda
e^{-\lambda{z}}&,~{}\text{if}~{}0\leq z<\alpha\\\\[7.11317pt] 2\lambda
e^{-2\lambda(z-\frac{1}{2}\alpha)}&,~{}\text{if}~{}\alpha\leq
z~{}~{}~{}.\end{cases}$ (24)

Using Eq. 24, the mean value of $Z_{(1)}$ is

$\mathbb{E}[Z_{(1)}]=\int_{-\infty}^{\infty}zf_{Z_{(1)}}(z)\,dz=\frac{1}{\lambda}-\frac{1}{2\lambda}e^{-\lambda\alpha}~{}~{}~{}\text{(sec/symbol)}.$
(25)

Maximum. Let random variable $Z_{(2)}$ be the second-order statistic defined
as

$Z_{(2)}=\max\\{Z_{1},Z_{2}\\}~{}.$ (26)

The pdf of $Z_{(2)}$ has two cases:

Case 1: $z<\alpha$

$Z_{(2)}=\max\\{Z_{1},Z_{2}\\}=\varnothing~{}$ (27)

because for two neurons the maximum, that is the second-order statistic, needs
to be greater or equal to $\alpha$ (Eq. 15). Therefore,

$\begin{split}f_{Z_{(2)}}(z)=&~{}0~{}.\end{split}$ (28)

Case 2: $z\geq\alpha$

$Z_{(2)}=\max\\{Z_{1},Z_{2}\\}.$ (29)

The cdf of the maximum is

$\begin{split}F_{Z_{(2)}}(z)=&P\left(Z_{(2)}\leq z\right)\\\ =&P(Z_{1}\leq
z,Z_{2}\leq z)\\\ =&P(Z_{1}\leq z)(Z_{2}\leq z)~{}~{}~{}~{}\text{because of
independence}\\\ =&1-e^{-\lambda
z}-e^{-\lambda(z-\alpha)}+e^{-2\lambda(z-\frac{1}{2}\alpha)}~{}~{}.\end{split}$
(30)

Accordingly, the pdf of the maximum is

$f_{Z_{(2)}}(z)=\frac{d}{dz}F_{Z_{(2)}}(z)=\lambda e^{-\lambda z}+\lambda
e^{-\lambda(z-\alpha)}-2\lambda e^{-2\lambda(z-\frac{1}{2}\alpha)}~{}~{}.$
(31)

Combining Case 1 and 2:

$f_{Z_{(2)}}(z)=\begin{cases}0&,~{}\text{if}~{}z<\alpha\\\\[7.11317pt] \lambda
e^{-\lambda z}+\lambda e^{-\lambda(z-\alpha)}-2\lambda
e^{-2\lambda(z-\frac{1}{2}\alpha)}&,~{}\text{if}~{}\alpha\leq
z~{}~{}~{}.\end{cases}$ (32)

Using Eq. 32, the mean value of $Z_{(2)}$ is

$\mathbb{E}[Z_{(2)}]=\int_{-\infty}^{\infty}zf_{Z_{(2)}}(z)\,dz=\frac{\lambda\alpha+1}{\lambda}+\frac{1}{2\lambda}e^{-\lambda\alpha}~{}~{}~{}\text{(sec/symbol)}.$
(33)

Using Eq. 25 and 33, the average symbol duration for two neurons is

$\bar{T}=\mathbb{E}~{}[T]=\mathbb{E}~{}[Z_{(2)}-Z_{(1)}]=\mathbb{E}~{}[Z_{(2)}]-\mathbb{E}~{}[Z_{(1)}]=\alpha+\frac{1}{\lambda}e^{-\lambda\alpha}~{}~{}\text{(sec/symbol)}.$
(34)

Fig. 10 compares theory (Eq. 34) with simulation.

Figure 10: Theory (Eq. 34) versus simulation for the average symbol duration,
$\bar{T}$. The number of presynaptic neurons is $n=2$. Here $\alpha$ is
arbitrarily set to 1 (sec), and the number of samples per point used in the
simulation is $10^{9}$.

###  C-B Three Presynaptic Neurons

Minimum. Let random variable $Z_{(1)}$ be the first-order statistic defined as

$Z_{(1)}=\min\\{Z_{1},Z_{2},Z_{3}\\}~{}.$ (35)

The pdf of $Z_{(1)}$ has four cases:

Case 1: $z<0$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2},Z_{3}\\}=\varnothing~{}~{}\\\\[7.11317pt]
\end{split}$ (36)

because if $z<0$, then $f_{Z_{1}}(z)=f_{Z_{2}}(z)=f_{Z_{3}}(z)=0$, (Eq. 15).
Therefore,

$f_{Z_{(1)}}(z)=0~{}.$ (37)

Case 2: $0\leq z<\alpha$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2},Z_{3}\\}=Z_{1}~{}~{}\\\\[7.11317pt]
\end{split}$ (38)

because in the range $0\leq z<\alpha$, $f_{Z_{2}}(z)=f_{Z_{3}}(z)=0$, (Eq.
15). Therefore,

$f_{Z_{(1)}}(z)=\lambda e^{-\lambda{z}}~{}.$ (39)

Case 3: $\alpha\leq z<2\alpha$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2},Z_{3}\\}=\min\\{Z_{1},Z_{2}\\}~{}\end{split}$
(40)

because in the range $\alpha\leq z<2\alpha$, $f_{Z_{3}}(z)=0$, (Eq. 15).

The cdf of the minimum is

$\begin{split}F_{Z_{(1)}}(z)=&P\left(Z_{(1)}\leq z\right)\\\
=&1-P\left(Z_{(1)}>z\right)\\\ =&1-P(Z_{1}>z,Z_{2}>z)\\\
=&1-P(Z_{1}>z)(Z_{2}>z)~{}~{}~{}~{}\text{because of independence}\\\
=&1-\Bigl{[}\Bigl{(}1-P(Z_{1}\leq z)\Bigr{)}\Bigl{(}1-P(Z_{2}\leq
z)\Bigr{)}\Bigr{]}\\\
=&1-e^{-2\lambda(z-\frac{1}{2}\alpha)}~{}~{}.\end{split}$ (41)

Accordingly, the pdf of the minimum is

$f_{Z_{(1)}}(z)=\frac{d}{dz}F_{Z_{(1)}}(z)=2\lambda
e^{-2\lambda(z-\frac{1}{2}\alpha)}~{}~{}.$ (42)

Case 4: $2\alpha\leq z$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2},Z_{3}\\}.\end{split}$ (43)

The cdf of the minimum is

$\begin{split}F_{Z_{(1)}}(z)=&P\left(Z_{(1)}\leq z\right)\\\
=&1-P\left(Z_{(1)}>z\right)\\\ =&1-P(Z_{1}>z,Z_{2}>z,Z_{3}>z)\\\
=&1-P(Z_{1}>z)(Z_{2}>z)(Z_{3}>z)~{}~{}~{}~{}\text{because of independence}\\\
=&1-\Bigl{[}\Bigl{(}1-P(Z_{1}\leq z)\Bigr{)}\Bigl{(}1-P(Z_{2}\leq
z)\Bigr{)}\Bigl{(}1-P(Z_{3}\leq z)\Bigr{)}\Bigr{]}\\\
=&1-e^{-3\lambda(z-\alpha)}~{}~{}.\end{split}$ (44)

Accordingly, the pdf of the minimum is

$f_{Z_{(1)}}(z)=\frac{d}{dz}F_{Z_{(1)}}(z)=3\lambda
e^{-3\lambda(z-\alpha)}~{}~{}.$ (45)

Combining Case 1 to 4:

$f_{Z_{(1)}}(z)=\begin{cases}0&,~{}\text{if}~{}z<0\\\\[7.11317pt] \lambda
e^{-\lambda{z}}&,~{}\text{if}~{}0\leq z<\alpha\\\\[7.11317pt] 2\lambda
e^{-2\lambda(z-\frac{1}{2}\alpha)}&,~{}\text{if}~{}\alpha\leq
z<2\alpha\\\\[7.11317pt] 3\lambda
e^{-3\lambda(z-\alpha)}&,~{}\text{if}~{}2\alpha\leq z~{}~{}.\par\end{cases}$
(46)

Using Eq. 46, the mean value of $Z_{(1)}$ is

$\mathbb{E}[Z_{(1)}]=\int_{-\infty}^{\infty}zf_{Z_{(1)}}(z)\,dz=\frac{1}{\lambda}-\frac{1}{2\lambda}e^{-\lambda\alpha}-\frac{1}{6\lambda}e^{-3\lambda\alpha}~{}~{}~{}\text{(sec/symbol)}.$
(47)

Maximum. Let random variable $Z_{(3)}$ be the third-order statistic defined as

$Z_{(3)}=\max\\{Z_{1},Z_{2},Z_{3}\\}~{}.$ (48)

The pdf of $Z_{(3)}$ has two cases:

Case 1: $z<2\alpha$

$Z_{(3)}=\max\\{Z_{1},Z_{2},Z_{3}\\}=\varnothing~{}$ (49)

because for three neurons the maximum, that is the third-order statistic,
needs to be greateror equal to $2\alpha$ (Eq. 15). Therefore,

$\begin{split}f_{Z_{(3)}}(z)=&~{}0~{}.\end{split}$ (50)

Case 2: $z\geq 2\alpha$

$Z_{(3)}=\max\\{Z_{1},Z_{2},Z_{3}\\}.$ (51)

The cdf of the maximum is

$\begin{split}F_{Z_{(3)}}(z)=&P\left(Z_{(3)}\leq z\right)\\\ =&P(Z_{1}\leq
z,Z_{2}\leq z,Z_{3}\leq z)\\\ =&P(Z_{1}\leq z)(Z_{2}\leq z)(Z_{3}\leq
z)~{}~{}~{}~{}\text{because of independence}\\\
=&1-e^{-\lambda(z-\alpha)}-e^{-\lambda
z}+e^{-\lambda(2z-\alpha)}-e^{-\lambda(z-2\alpha)}+e^{-\lambda(2z-3\alpha)}+e^{-2\lambda(z-\alpha)}-e^{-3\lambda(z-\alpha)}~{}~{}.\end{split}$
(52)

Accordingly, the pdf of the maximum is

$f_{Z_{(3)}}(z)=\frac{d}{dz}F_{Z_{(3)}}(z)=\lambda
e^{-\lambda(z-\alpha)}+\lambda e^{-\lambda z}-2\lambda
e^{-\lambda(2z-\alpha)}+\lambda e^{-\lambda(z-2\alpha)}-2\lambda
e^{-\lambda(2z-3\alpha)}-2\lambda e^{-2\lambda(z-\alpha)}+3\lambda
e^{-3\lambda(z-\alpha)}.$ (53)

Combining Case 1 and 2:

$f_{Z_{(3)}}(z)=\begin{cases}0&,~{}\text{if}~{}z<2\alpha\\\\[7.11317pt]
\lambda e^{-\lambda(z-\alpha)}+\lambda e^{-\lambda z}-2\lambda
e^{-\lambda(2z-\alpha)}+\lambda e^{-\lambda(z-2\alpha)}-2\lambda
e^{-\lambda(2z-3\alpha)}-2\lambda e^{-2\lambda(z-\alpha)}+3\lambda
e^{-3\lambda(z-\alpha)}&,~{}\text{if}~{}2\alpha\leq z\end{cases}$ (54)

Using Eq. 54, the mean value of $Z_{(3)}$ is

$\mathbb{E}[Z_{(3)}]=\int_{-\infty}^{\infty}zf_{Z_{(3)}}(z)\,dz=\frac{2\lambda\alpha+1}{\lambda}+\frac{1}{2\lambda}e^{-\lambda\alpha}+\frac{1}{2\lambda}e^{-2\lambda\alpha}-\frac{1}{6\lambda}e^{-3\lambda\alpha}~{}~{}~{}\text{(sec/symbol)}.$
(55)

Using Eq. 47 and 55, the average symbol duration for three neurons is

$\bar{T}=\mathbb{E}~{}[T]=\mathbb{E}~{}[Z_{(3)}-Z_{(1)}]=\mathbb{E}~{}[Z_{(3)}]-\mathbb{E}~{}[Z_{(1)}]=2\alpha+\frac{1}{\lambda}e^{-\lambda\alpha}+\frac{1}{2\lambda}e^{-2\lambda\alpha}~{}~{}\text{(sec/symbol)}.$
(56)

Fig. 11 compares theory (Eq. 56) with simulation.

Figure 11: Theory (Eq. 56) versus simulation for the average symbol duration,
$\bar{T}$. The number of presynaptic neurons is $n=3$. Here $\alpha$ is
arbitrarily set to 1 (sec), and the number of samples per point used in the
simulation is $10^{9}$.

###  C-C Four Presynaptic Neurons

Minimum. Let random variable $Z_{(1)}$ be the first-order statistic defined as

$Z_{(1)}=\min\\{Z_{1},Z_{2},Z_{3},Z_{4}\\}~{}.$ (57)

The pdf of $Z_{(1)}$ has five cases:

Case 1: $z<0$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2},Z_{3},Z_{4}\\}=\varnothing~{}~{}\\\\[7.11317pt]
\end{split}$ (58)

because if $z<0$, then
$f_{Z_{1}}(z)=f_{Z_{2}}(z)=f_{Z_{3}}(z)=f_{Z_{4}}(z)=0$, (Eq. 15). Therefore,

$f_{Z_{(1)}}(z)=0~{}.$ (59)

Case 2: $0\leq z<\alpha$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2},Z_{3},Z_{4}\\}=Z_{1}~{}~{}\\\\[7.11317pt]
\end{split}$ (60)

because in the range $0\leq z<\alpha$,
$f_{Z_{2}}(z)=f_{Z_{3}}(z)=f_{Z_{4}}(z)=0$, (Eq. 15). Therefore,

$f_{Z_{(1)}}(z)=\lambda e^{-\lambda{z}}~{}.$ (61)

Case 3: $\alpha\leq z<2\alpha$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2},Z_{3},Z_{4}\\}=\min\\{Z_{1},Z_{2}\\}~{}\end{split}$
(62)

because in the range $\alpha\leq z<2\alpha$, $f_{Z_{3}}(z)=f_{Z_{4}}(z)=0$ .

The cdf of the minimum is

$\begin{split}F_{Z_{(1)}}(z)=&P\left(Z_{(1)}\leq z\right)\\\
=&1-P\left(Z_{(1)}>z\right)\\\ =&1-P(Z_{1}>z,Z_{2}>z)\\\
=&1-P(Z_{1}>z)(Z_{2}>z)~{}~{}~{}~{}\text{because of independence}\\\
=&1-\Bigl{[}\Bigl{(}1-P(Z_{1}\leq z)\Bigr{)}\Bigl{(}1-P(Z_{2}\leq
z)\Bigr{)}\Bigr{]}\\\
=&1-e^{-2\lambda(z-\frac{1}{2}\alpha)}~{}~{}.\end{split}$ (63)

Accordingly, the pdf of the minimum is

$f_{Z_{(1)}}(z)=\frac{d}{dz}F_{Z_{(1)}}(z)=2\lambda
e^{-2\lambda(z-\frac{1}{2}\alpha)}~{}~{}~{}.$ (64)

Case 4: $2\alpha\leq z<3\alpha$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2},Z_{3},Z_{4}\\}=\min\\{Z_{1},Z_{2},Z_{3}\\}\end{split}$
(65)

because in the range $2\alpha\leq z<3\alpha$, $f_{Z_{4}}(z)=0$.
The cdf of the minimum is

$\begin{split}F_{Z_{(1)}}(z)=&P\left(Z_{(1)}\leq z\right)\\\
=&1-P\left(Z_{(1)}>z\right)\\\ =&1-P(Z_{1}>z,Z_{2}>z,Z_{3}>z)\\\
=&1-P(Z_{1}>z)(Z_{2}>z)(Z_{3}>z)~{}~{}~{}~{}\text{because of independence}\\\
=&1-\Bigl{[}\Bigl{(}1-P(Z_{1}\leq z)\Bigr{)}\Bigl{(}1-P(Z_{2}\leq
z)\Bigr{)}\Bigl{(}1-P(Z_{3}\leq z)\Bigr{)}\Bigr{]}\\\
=&1-e^{-3\lambda(z-\alpha)}~{}~{}.\end{split}$ (66)

Accordingly, the pdf of the minimum is

$f_{Z_{(1)}}(z)=\frac{d}{dz}F_{Z_{(1)}}(z)=3\lambda
e^{-3\lambda(z-\alpha)}~{}~{}.$ (67)

Case 5: $3\alpha\leq z$

$\begin{split}Z_{(1)}=&\min\\{Z_{1},Z_{2},Z_{3},Z_{4}\\}.\end{split}$ (68)

The cdf of the minimum is

$\begin{split}F_{Z_{(1)}}(z)=&P\left(Z_{(1)}\leq z\right)\\\
=&1-P\left(Z_{(1)}>z\right)\\\ =&1-P(Z_{1}>z,Z_{2}>z,Z_{3}>z,Z_{4}>z)\\\
=&1-P(Z_{1}>z)(Z_{2}>z)(Z_{3}>z)(Z_{4}>z)~{}~{}~{}~{}\text{because of
independence}\\\ =&1-\Bigl{[}\Bigl{(}1-P(Z_{1}\leq
z)\Bigr{)}\Bigl{(}1-P(Z_{2}\leq z)\Bigr{)}\Bigl{(}1-P(Z_{3}\leq
z)\Bigr{)}\Bigl{(}1-P(Z_{4}\leq z)\Bigr{)}\Bigr{]}\\\
=&1-e^{-4\lambda(z-\frac{3}{2}\alpha)}~{}~{}.\end{split}$ (69)

Accordingly, the pdf of the minimum is

$f_{Z_{(1)}}(z)=\frac{d}{dz}F_{Z_{(1)}}(z)=4\lambda
e^{-4\lambda(z-\frac{3}{2}\alpha)}~{}~{}.$ (70)

Combining Case 1 to 5:

$f_{Z_{(1)}}(z)=\begin{cases}0&,~{}\text{if}~{}z<0\\\\[7.11317pt] \lambda
e^{-\lambda{z}}&,~{}\text{if}~{}0\leq z<\alpha\\\\[7.11317pt] 2\lambda
e^{-2\lambda(z-\frac{1}{2}\alpha)}&,~{}\text{if}~{}\alpha\leq
z<2\alpha\\\\[7.11317pt] 3\lambda
e^{-3\lambda(z-\alpha)}&,~{}\text{if}~{}2\alpha\leq z<3\alpha\\\\[7.11317pt]
4\lambda e^{-4\lambda(z-\frac{3}{2}\alpha)}&,~{}\text{if}~{}3\alpha\leq
z~{}~{}.\end{cases}$ (71)

Using Eq. 71, the mean value of $Z_{(1)}$ is

$\mathbb{E}[Z_{(1)}]=\int_{-\infty}^{\infty}zf_{Z_{(1)}}(z)\,dz=\frac{1}{\lambda}-\frac{1}{2\lambda}e^{-\lambda\alpha}-\frac{1}{6\lambda}e^{-3\lambda\alpha}-\frac{1}{12\lambda}e^{-6\lambda\alpha}~{}~{}~{}\text{(sec/symbol)}.$
(72)

Maximum. Let random variable $Z_{(4)}$ be the fourth-order statistic defined
as

$Z_{(4)}=\max\\{Z_{1},Z_{2},Z_{3},Z_{4}\\}~{}.$ (73)

The pdf of $Z_{(4)}$ has two cases:

Case 1: $z<3\alpha$

$Z_{(4)}=\max\\{Z_{1},Z_{2},Z_{3},Z_{4}\\}=\varnothing~{}$ (74)

because for four neurons the maximum, that is the fourth-order statistic,
needs to be greater or equal to $3\alpha$ (Eq. 15). Therefore,

$\begin{split}f_{Z_{(4)}}(z)=&~{}0~{}~{}.\end{split}$ (75)

Case 2: $z\geq 3\alpha$

$Z_{(4)}=\max\\{Z_{1},Z_{2},Z_{3},Z_{4}\\}.$ (76)

The cdf of the maximum is

$\displaystyle F_{Z_{(4)}}(z)=$ $\displaystyle P\left(Z_{(4)}\leq z\right)$
$\displaystyle=$ $\displaystyle P(Z_{1}\leq z,Z_{2}\leq z,Z_{3}\leq
z,Z_{4}\leq z)$ $\displaystyle=$ $\displaystyle P(Z_{1}\leq z)(Z_{2}\leq
z)(Z_{3}\leq z)(Z_{4}\leq z)~{}~{}~{}~{}\text{because of independence}$
$\displaystyle=$ $\displaystyle 1-e^{-\lambda(z-\alpha)}-e^{-\lambda
z}+e^{-\lambda(2z-\alpha)}-e^{-\lambda(z-2\alpha)}+2e^{-\lambda(2z-3\alpha)}+e^{-2\lambda(z-\alpha)}-e^{-3\lambda(z-\alpha)}$
$\displaystyle\qquad+e^{-2\lambda(z-2\alpha)}-e^{-\lambda(3z-4\alpha)}+e^{-\lambda(2z-5\alpha)}-e^{-3\lambda(z-2\alpha)}-e^{-\lambda(3z-5\alpha)}+e^{-2\lambda(2z-3\alpha)}-e^{-\lambda(z-3\alpha)}.$
(77)

Accordingly, the pdf of the maximum is

$\displaystyle f_{Z_{(4)}}(z)=\frac{d}{dz}F_{Z_{(4)}}(z)=$
$\displaystyle\lambda e^{-\lambda(z-\alpha)}+\lambda e^{-\lambda z}-2\lambda
e^{-\lambda(2z-\alpha)}+\lambda e^{-\lambda(z-2\alpha)}-4\lambda
e^{-\lambda(2z-3\alpha)}$ $\displaystyle\qquad-2\lambda
e^{-2\lambda(z-\alpha)}+3\lambda e^{-3\lambda(z-\alpha)}-2\lambda
e^{-2\lambda(z-2\alpha)}+3\lambda e^{-\lambda(3z-4\alpha)}$
$\displaystyle\qquad-2\lambda e^{-\lambda(2z-5\alpha)}+3\lambda
e^{-3\lambda(z-2\alpha)}+3\lambda e^{-\lambda(3z-5\alpha)}-4\lambda
e^{-2\lambda(2z-3\alpha)}+\lambda e^{-\lambda(z-3\alpha)}~{}~{}~{}.$

Combining Case 1 and 2:

$f_{Z_{(4)}}(z)=\begin{cases}0&,~{}\text{if}~{}z<3\alpha\\\\[42.67912pt]
\lambda e^{-\lambda(z-\alpha)}+\lambda e^{-\lambda z}-2\lambda
e^{-\lambda(2z-\alpha)}+\lambda e^{-\lambda(z-2\alpha)}-4\lambda
e^{-\lambda(2z-3\alpha)}&,~{}\text{if}~{}3\alpha\leq z\\\\[7.11317pt]
\qquad-2\lambda e^{-2\lambda(z-\alpha)}+3\lambda
e^{-3\lambda(z-\alpha)}-2\lambda e^{-2\lambda(z-2\alpha)}+3\lambda
e^{-\lambda(3z-4\alpha)}&\\\\[7.11317pt] \qquad-2\lambda
e^{-\lambda(2z-5\alpha)}+3\lambda e^{-3\lambda(z-2\alpha)}+3\lambda
e^{-\lambda(3z-5\alpha)}-4\lambda e^{-2\lambda(2z-3\alpha)}+\lambda
e^{-\lambda(z-3\alpha)}~{}~{}.&\end{cases}$ (79)

Using Eq. 79, the mean value of $Z_{(4)}$ is

$\mathbb{E}[Z_{(4)}]=\int_{-\infty}^{\infty}zf_{Z_{(4)}}(z)\,dz=\frac{3\lambda\alpha+1}{\lambda}+\frac{1}{2\lambda}e^{-\lambda\alpha}+\frac{1}{2\lambda}e^{-2\lambda\alpha}+\frac{1}{3\lambda}e^{-3\lambda\alpha}-\frac{1}{6\lambda}e^{-4\lambda\alpha}-\frac{1}{6\lambda}e^{-5\lambda\alpha}+\frac{1}{12\lambda}e^{-6\lambda\alpha}~{}~{}\text{(sec/symbol)}.$
(80)

Using Eq. 72 and 80, the average symbol duration for four neurons is

$\displaystyle\bar{T}=\mathbb{E}~{}[T]=\mathbb{E}~{}[Z_{(4)}-Z_{(1)}]=\mathbb{E}~{}[Z_{(4)}]-\mathbb{E}~{}[Z_{(1)}]=3\alpha+\frac{1}{\lambda}e^{-\lambda\alpha}+\frac{1}{2\lambda}e^{-2\lambda\alpha}+\frac{1}{2\lambda}e^{-3\lambda\alpha}$
$\displaystyle-\frac{1}{6\lambda}e^{-4\lambda\alpha}-\frac{1}{6\lambda}e^{-5\lambda\alpha}+\frac{1}{6\lambda}e^{-6\lambda\alpha}~{}~{}~{}\text{(sec/symbol)}.$
(81)

Fig. 12 compares theory (Eq. C-C) with simulation.

Figure 12: Theory (Eq. C-C) versus simulation for the average symbol duration,
$\bar{T}$. The number of presynaptic neurons is $n=4$. Here $\alpha$ is
arbitrarily set to 1 (sec), and the number of samples per point used in the
simulation is $10^{9}$.


# Comparisons of Type Ia Supernova Light Curves in the UV and Optical
with the Swift Ultra-Violet/Optical Telescope

Yaswant Devarakonda  George P. and Cynthia Woods Mitchell Institute for
Fundamental Physics $\&$ Astronomy,
Texas A$\&$M University, Department of Physics and Astronomy,
4242 TAMU, College Station, TX 77843, USA  Peter J. Brown  George P. and
Cynthia Woods Mitchell Institute for Fundamental Physics $\&$ Astronomy,
Texas A$\&$M University, Department of Physics and Astronomy,
4242 TAMU, College Station, TX 77843, USA

(Accepted April 10, 2022)

###### Abstract

We examine the light curve parameters of 97 nearby Type Ia supernovae in the
ultraviolet and optical using observations from the Swift Ultra-Violet/Optical
Telescope. Our light curve models used a linear combinations of templates,
which were based on Functional Principal Component Analysis of the optical
photometry of SNe Ia. The weights for each principal component used in the fit
capture certain aspects of the light curve properties. We find that there is
little difference in the overall variability of principal component scores
across filters. We also find correlations between certain filters and PC
components, indicating that the UV and optical properties seem to be tied
together. This is a promising step in UV SNe Ia analysis, and suggests that UV
light curves may be used to infer optical properties, paving the way for
future cosmological studies.

Type Ia Supernovae

††journal: AJ

##  1 Introduction

Type Ia Supernovae (SNe Ia) are useful as cosmological tools due to their
highly uniform optical luminosity, which correlates strongly with their
optical light curve profile and color (Pskovskii, 1977; Phillips, 1993; Riess
et al., 1996; Phillips et al., 1999). This has led to their use in studies for
distance measurements and constraints on cosmological parameters (Riess et
al., 1998; Perlmutter et al., 1999).

In optical wavelengths, SNe Ia have a strong relationship between their peak
luminosity and luminosity-decline rate, leading to their widespread usage as
standardizable candles (Branch, 1998; Leibundgut, 2001). This process is
generally done by calibrating the peak luminosity with distance-independent
metrics such as the light curve shape. This can be done by measuring the
decline rate for the first 15 days after maximum brightness ($\Delta m_{15}$),
often in the B band (Phillips, 1993; Phillips et al., 1999; Garnavich et al.,
2004). Another common method is measuring the stretch of a light curve when
fitting to a template (Goldhaber et al., 2001).

The luminosity-width relation shows that more luminous SNe Ia have broader
light curves. The cause of this relation is believed to be the production of
${}^{56}Ni$ during the SNIa explosion, and the change in the temperature and
ionization evolution as a result (Nugent et al., 1995; Höflich et al., 1996;
Mazzali et al., 2001; Kasen & Woosley, 2007). ${}^{56}Ni$ drives much of the
SNIa luminosity, and the light curve shape is greatly affected by opacity
effects from${}^{56}Ni$ and Fe-group elements (Mazzali et al., 2007).

Subsequent studies have looked beyond the optical to see if these strong
correlations still hold true at other wavelengths. For example, the near-
infrared (NIR) luminosity has been shown to be less susceptible to the effects
of interstellar extinction while also having a diminished dependence on the
shape of the light curve (Krisciunas et al., 2000; Wood-Vasey et al., 2008;
Mandel et al., 2011; Kattner et al., 2012). Meanwhile, studies in the near-
ultraviolet have shown that the correlations have much larger scatter (Jha et
al., 2006; Kessler et al., 2009; Brown et al., 2010; Milne et al., 2010; Brown
et al., 2017). Understanding the scatter at shorter wavelengths will be an
important task, as it enables us to better predict SNe Ia properties at larger
distances where their rest-frame UV emission will be redshifted into the
optical.

One issue with analyzing the UV light curves of SNe Ia is the lack of
empirical templates at those wavelengths. In the optical, the most common
optical fitting techniques are MLCS (Riess et al., 1996; Jha et al., 2007) and
SALT3 (Kenworthy et al., 2021). MLCS directly models the light curve data
using vectorized templates over a hypothesized grid, whereas SALT3 uses
spectral energy distributions to model light curves. Here, we use templates
based on functional principal component analysis (FPCA) of optical SNe Ia
light curves (He et al., 2018). We apply these templates to UV and optical
observations of 97 SNe Ia to examine how well this technique can characterize
light curves in a regime that it wasn’t trained in, and if the fitted
parameters can shed new light into the diversity of SNe Ia in the UV.

In Section 2 we describe our observations and methodology for template
fitting. Then in Section 3 we apply statistical tests to compare the diversity
and correlation of our fitted parameters in the UV and optical regimes.
Finally, in Section 4 we examine our results in detail and summarize our work
in Section 5.

##  2 Methods

###  2.1 Observations

Observations were done with the Ultra-Violet Optical Telescope (UVOT; Roming
et al., 2005)on the Neil Gehrels Swift Observatory (Gehrels et al., 2004).
Observations were usually performed in all six medium-band UVOT filters:
$UVW2$, $UVM2$, $UVW1$, $U$, $B$, and $V$. See Poole et al. (2008) for
calibration and filter curve details, and Brown et al. (2010) for the
effective throughput for a SN Ia. In total, 219 SNe Ia were observed by Swift
UVOT between March 2005 and October 2016 (SN 2005am through 2016gxp,
inclusive) with photometry available through the Swift Optical Ultra-Violet
Supernova Archive (SOUSA; Brown et al., 2014). A 3” or 5” aperture was used to
maximise the signal-to-noise and zeropoints were based on the updated
calibration of Breeveld et al. (2011) and the time-sensitivity correction from
July 2015.

Upon manual inspection of the galaxy-subtracted UVOT light curves, we selected
110 SNe Ia which had at least four data points in the $UVW1$ filter covering
the peak and early decline. Due to the faintness of SNe Ia in the mid-UV, some
did not have adequate light curves in $UVW2$ and $UVM2$, and some were
saturated in the $U$, $B$, and $V$ filters at peak. The light curves were
initially fit following the procedure in Brown et al. (2017). Template light
curves were fit to our data, using the peak magnitude, time of maximum, and
stretch (linear scaling in the time axis) as free parameters. We used the $B$
and $V$ template curves from MLCS2k2 (Jha et al., 2007) and the $UVW2$,
$UVM2$, and $UVW1$ light curves of SN 2011fe (Brown et al., 2012). The $UVW1$
light curve of SN 2011fe was also used to fit the $U$ band light curves due to
the similarity of the mean light curves of $U$ and $UVW1$ bands found by Milne
et al. (2010) and the saturation of the Swift $U$ band light curve of SN
2011fe near peak.

###  2.2 FPCA Model

To model our sample, we used templates constructed from functional principal
component analysis (FPCA, He et al., 2018), a statistical method that uses a
linear combination of a mean function and a few principal component (PC)
functions to represent the light curve. The coefficients associated with these
functions are called PC scores. The mean light curve function and PC functions
were trained using a data set of light curve photometry for 111 SNe Ia in the
$BVRI$ bands. A four PC template was created for each band, as well as a
generic (band-vague) template that combined the training data across bands
that can be used to model light curves in any filter.

To retrain the FPCA model to include UV photometry would require more quality
observations and an amount of work beyond the scope of this paper. Rather, we
take the simpler step of using the band-vague template to fit the UV and
optical light curves simultaneously. While the inclusion of the R and I band
in the training set leads to an over-prediction of flux at later phases, we
find (2.4) that the quality of the fit is comparable between the UV and
optical bands in the regions near peak flux. Using the same template to fit
across all bands has the added benefit of each PC controlling the same aspects
of each light curve, allowing for direct comparisons. In Section 4 we examine
possible future work with UV trained modeling.

Per He et al. (2018), the first PC reflects the decline rate after the peak,
and the second PC reflects the light curve width around the peak. The third PC
adjusts the pre-peak rise as well as a bump 20 days after peak that is
commonly seen in the near-IR bands, while the fourth PC accounts for more
complex fluctuations. While the first three components are tied to certain
light curve properties, they also describe a wider range of parameters. This
means that while the first PC is related to the decline rate, it is not a
perfect one-to-one correlation. In the original data set that this template
was tested on the first PC score explained 92$\%$ of variability, the first
two explained 96$\%$, and all four explained 99$\%$ of variability. In Figure
1 we replicate the effects of a typical range of PC weights for visualization.

Figure 1: Examples of the effect that each PC template has on the light curve
for a range of common B band PC score weights. The PC1 template mostly
controls the late-phase decline rate, the PC2 template controls the overall
stretch of the light curve and the post-peak decline rate, and the PC3
template has a plays a role in the pre-peak rise rate and the late phase flux.

###  2.3 Template Fitting

We calculated the phase of each light curve in all six filters, corrected for
redshift time dilation, by using the estimated peak from the preliminary
modeling as our zero point.

$Phase_{i}=\frac{MJD_{i}-MJD_{peak}}{1+z}$ (1)

We then removed data points that occur 10 days before the peak or 40 days
after the peak, as those are outside of the domain of the template. We ensured
that there was at least one data point before the peak and one at least 15
days after the peak, and that there was a minimum of five data points. The
template models were then linearly interpolated to match the phases of the
observations, and linearly combined to match the observed light curve with
weights for the templates determined by the curve_fit package from the SciPy
library in python (Virtanen et al., 2020). The program attempts to fit the
function to the data by minimizing the non-linear least squared errors. Below
is an example of the function,

$LC(p)=M_{peak}+Mean_{template}(p)+\sum_{i=1}^{N}W_{i}*PC_{i}(p)$ (2)

where $p$ is the phase, $M_{peak}$ is the peak magnitude (a free parameter,
but given an initial guess based on the preliminary modeling estimates),
$Mean_{template}$ is the mean light curve from the template, and $W_{i}$ is
the weight of each principal component as determined by curve_fit.

Figure 2: Swift observations of SN 2011by (blue) and the 3 template model
(orange) with 1 $\sigma$ errors (red). A phase of 0 in each subplot
corresponds to the date of the modeled peak brightness for each band.

In total, there were 97 SNe Ia that had a fitted light curve in at least one
filter, and 14 SNe Ia that had a fitted light curve in every filter. The
$UVM2$ filter had the fewest well-fit light curves, due to the lower S/N of
observations compared to other filters. Uncertainties in our model fits were
determined by examining the full covariance matrix that results output by the
curve_fit package. We measured the $R^{2}$ values for the modeled SNe Ia
(Figure 3) as a test of goodness-of-fit, and find that the majority of SNe Ia
behavior is well explained by our model. The $U$ and $UVW1$ curves (which were
not included in the training of the band-vague fitting) are fit with similar
accuracy to the $B$ and $V$ band, respectively, and the PC fits explain 95% of
the points for 75% of the mid-UV ($UVM2$ and $UVW2$) light curves.

Figure 3: The $R^{2}$ values of modeled SNe Ia in each band. This value
measures the amount of variance in the data that is explained by the model on
a scale of 0-1. This measurement is done without accounting for the known
observational error in the data, and is thus a lower limit on the true $R^{2}$
value. The notches on each box plot indicate the uncertainty in the median,
and the whiskers indicate 1.5$\times$IQR distance.

###  2.4 Model Selection

To determine how many principal components were needed to ensure a good
quality fit for our data set, we employed several tests of model selection.
From the 14 SNe Ia with fits in all six bands, we selected four that had a
relatively large number of observations (SN 2005ke, SN 2007af, SN 2007on, SN
2011by) to serve as test cases. Using the train_test_split package in the
scikit-learn library in Python (Pedregosa et al., 2011), we performed 100
bootstrap draws from each SNIa with a 70/30 split between the training and
testing sample sizes. This allowed us to test the robustness of our fit with
light curves that have sparse or inconsistent sampling, which is a more
realistic case for the majority of our observations.

Figure 4: An example of one of the bootstrap samples used for estimating the
robustness of a three-PC model vs a four-PC model, using SN 2011by. After the
model parameters are estimated using a training set, the model is then fit to
a testing set and checked for goodness of fit. Increasing the number of
templates used in the model from three to four does little to improve the fit
of the model to the data.

The Akaike Information Criterion (AIC) (Akaike, 1974) was calculated for each
bootstrap sample and the mean AIC score for the two-, three-, and four-PC
model were compared. The AIC score is a comparison test that measures how well
a model fits a data set given a certain level of model complexity. In general,
a lower AIC score indicates a better model. We employed a modified version of
the AIC to account for the relatively small sample size, as shown in Equation
3.

$\begin{split}AIC&=2k-2ln(\hat{L})\\\ AIC_{small\
sample}&=AIC+\frac{2k^{2}+2k}{n-k-1}\end{split}$ (3)

$\hat{L}$ is the maximum value of the likelihood function of the model, $k$ is
the number of parameters, and $n$ is the sample size. We also explored using
the Bayesian Information Criterion, but determined that our sample sizes were
too small for an accurate assessment.

Overall, we found that the two-PC and three-PC models had the lowest mean AIC
scores in all six filters, signifying them as more appropriate models than the
four PC model to use for our analysis. Between the two-PC and three-PC models,
neither one had a significant difference when examining the larger data set.
Future work with the FPCA template fitting could include individualized
fitting for each SNIa and filter with more rigorous testing, but for uniform
fitting we chose the three PC model. It added robustness to our modeling by
improving our estimates of the rise of the light curves, with a minor risk of
overfitting certain SNe Ia.

Table 1: A comparison of the mean and standard deviation of the AIC scores taken from 100 random samples of SN 2011by data. In general, a lower AIC score indicates a more favorable model.  Filter | Training Set | Testing Set | 3 Template | 3 Template | 4 Template | 4 Template
---|---|---|---|---|---|---
Type | Sample Size | Sample Size | Mean | $\sigma$ | Mean | $\sigma$
W2 | 22 | 10 | 5.4 | 0.7 | 12.3 | 0.7
M2 | 18 | 8 | 6.3 | 1.4 | 17.0 | 1.8
W1 | 18 | 9 | 4.9 | 1.3 | 12.8 | 1.4
U | 23 | 10 | 3.5 | 1.4 | 10.4 | 1.5
B | 23 | 10 | 3.0 | 1.7 | 9.6 | 1.9
V | 16 | 8 | 4.0 | 1.5 | 14.0 | 1.6

##  3 Analysis and Results

We employed a series of statistical tests to examine the similarity in the
distribution of PC scores across the different bands. The Brown-Forsythe test
(Brown & Forsythe, 1974) examines the differences in the variance of two data
sets by examining the distribution around the median. The null hypothesis for
this test is that there is no statistical difference in the variance, while a
p-value below 0.05 would indicate there is a difference in variance (Section
7).

Additionally, we visually inspected the data variability with box plots (Fig.
5)(Tukey, 1977). Boxplots allow us to identify quickly and assess outliers,
while also being less susceptible to issues such as sample size and skewness
that may bias p-values. In particular, we examined the interquantile range
(IQR) of the different PC scores in each band. Similar IQR scores would
suggest that the populations have similar variances. Visual analysis of the
outliers revealed that many had models with acceptable fit quality, indicating
that our model is capable of fitting the full diversity of our sample.

Figure 5: Box plots depicting the spread of parameters in each band. The
notches indicate the uncertainty in the median, and the whiskers indicate
1.5$\times$IQR distance. Light curves with parameter estimates beyond the
whiskers are labeled as circles.

We find little statistical difference between the spread of PC1 and PC2 scores
for both the UV and optical bands when compared to the $B$ band PC scores. For
PC3 weights, only the $V$ band has a similar spread in values to the $B$ band,
all other bands show a much larger spread in potential values. As both the PC1
and PC2 scores relate to the near-peak decline rate, nether can be solely used
as a proxy for $\Delta m_{15}$. However the full model can be used to directly
estimate the $\Delta m_{15}$, and the derived values show similar diversity
with the $B$ band values (except for the $V$ band, which has a much tighter
spread of values).

To examine the one-to-one correlations between filters, we plotted the
parameter values in each filter versus the $B$ band value. We then estimated
the Spearman rank-order coefficient (Kokoska & Zwilinger, 2000), a non-
parametric measure of two sets for a monotonic relationship. For parameters
that have a p-value below 0.05, we test out whether a linear model can
describe the relationship. While the true relationship may not be linear,
having the linear fit can be a useful visual representation of the correlation
and distribution of parameters. To perform this linear regression we employed
the Python routine LINMIX, which uses Bayesian inference to return random
draws from the posterior distribution (Kelly, 2007). The routine uses a 1,000
iteration Markov Chain Monte Carlo procedure to converge on the posterior. The
benefit of using this routine is its ability to handle uncertainties in both
the independent and dependent variables as well as intrinsic scatter. From the
iterations, we take the average value of the slope and intercept as our best
fit parameters, and estimate the $r^{2}$ value to see how well the linear fit
describes the correlation. The $r^{2}$ value is based on the total sum of
squares ($SS_{tot}$) and the sum of squared residuals ($SS_{res}$) of the
observed ($y$) and modeled ($f$) data, where

$\begin{split}SS_{tot}&=\sum_{i}(y_{i}-\hat{y})^{2}\\\
SS_{res}&=\sum_{i}(y_{i}-f_{i})^{2}\\\
r^{2}&=1-\frac{SS_{res}}{SS_{tot}}\end{split}$ (4)

We find that while the optical properties ($U$ and $V$ bands) are highly
related to the $B$ band properties, only the PC2 weights (Fig. 9) show a
consistent correlation between the mid-UV and near-UV bands and the $B$ band.
For PC1 (Fig. 8), there is a strong correlation between the $UVW2$ band and
the $B$ band, which may be due to the relatively higher quality and quantity
of data compared to the $UVM2$ and $UVW1$ bands.

There is little to no statistical correlation for the PC3 weights (Fig. 10)
between the UV bands and the $B$ band, although the $UVW1$ could be a
statistical edge case that lacks enough quality data to make a solid
statement. The $\Delta m_{15}$ (Fig. 11) shows a strong correlation and a
relatively tight linear relation with the $B$ band in both the UV and the
optical.

Figure 6: The comparison between the B band and W2 band PC2 weights for our
sample of selected SNe Ia. The red lines indicate the range of possible
regression fits from the LINMIX routine, and the black line indicates the
linear fit based on the mean of the LINMIX coefficients. See Figures 8, 9, 10,
and 11 in the appendix for the full set of correlated properties.

Finally, we examined the differences in the peak magnitude phase between the
$B$ band and the other filters, using a similar analysis to the PC score
distributions. We find (Fig 7) that the $U$ band and $UVW2$ band light curves
have a tight spread of peaks phases at around 2 days before the $B$ band peak,
and the $V$ band has a tight spread at around 2 days after the $B$ band peak
(Figure 7). The $UVM2$ and $UVW1$ band light curves also seem to peak earlier
than the $B$ band light curves, although there is a much larger spread in the
distribution of peak phases.

Examining the outliers, we find that SNe Ia that peak in the UV 4 or more days
before the B band peak were more likely to be classified as over-luminous in
the literature. The most extreme outliers, SN 2011aa (Kamiya, 2012), SN 2012dn
(Copin et al., 2012), and LSQ12gdj (Scalzo et al., 2014), are all super-
Chandrasekhar mass candidates. Thus, the difference in the time of peak in the
UV and optical may be a useful metric for sub-type classification.

##  4 Discussion

We successfully fit UV observations of SNe Ia using templates constructed from
functional principal component analysis of optical photometry. While the
quality of fit is better for the optical observations, quality of the UV fits
were sufficient to perform a statistical analysis of the SNe Ia light curve
properties. With improved UV observations and a larger data set, it should be
possible to create filter-specific UV templates that will greatly improve
future analysis of SNe Ia light curves.

We find the the overall spread of the SNe Ia parameters is similar between the
UV and optical bands in our sample, particularly for the PC2 template weights
which are most directly tied to the stretch and decline rates of the SNe Ia
light curve. In addition, we offer strong statistical evidence for a
correlation between the UV and optical PC2 weights, suggesting that we can
directly interpolate the spread of the SNe Ia light curve in one band based on
the other. This will be highly useful with upcoming large survey telescopes
such as the Rubin observatory, which is expected to find millions of SNe Ia at
a wide range of redshifts. Interpreting the rest frame UV light from distant
SNe Ia will require further understanding of the UV emission from nearby
objects.

The lack of a strong correlation in the PC1 weights could suggest that physics
behind the UV and optical light curves may be different at later phases
($\geq$20), where the PC1 weights have the largest impact. For phases
immediately after the peak the PC2 weights have the strongest effect on the
decline rate, resulting in the strong correlation with the $\Delta m_{15}$
rates. The correlations between PC2 weights in different bands suggests a
common origin for the near-peak variations in the UV and optical light curves.
Further studies with multi-epoch UV spectra of SNe Ia could provide more
details into the emission mechanisms and the variations of SNe Ia light
curves.

Previous work in the field has focused on rest-frame optical emission, leading
to a wealth of robust models and large data sets. In comparison, there is a
lack of quality observations and modeling in the UV. Our findings suggest that
the UV and optical light curve properties may be more similar than previously
suggested, which suggests the worth of further study into this regime. Future
work will greatly constrain the true behavior of UV SNe Ia light curves
relative to the better understood optical regime.

These correlations will be important for the study of high redshift SNe Ia,
whose rest-frame UV and optical emission will be simultaneously observable
with upcoming ground-based telescopes such as the Vera Rubin Observatory. In
the case that the optical light curve is poorly sampled (regardless of
distance), the UV light curve may be able to provide constraints on optical
properties such as the decline rate or stretch. The UV emission’s potential
for sub-type classification may, with further study, also assist in
standardizing the optical relations. The UV has also been shown to
significantly improve estimates on the amount of reddening due to dust
extinction when used in conjunction with optical and NIR data (Amanullah et
al., 2015), further showing the usefulness of multi-band studies.The diversity
of the UV absolute magnitudes remains poorly constrained, however, hindering
cosmological studies based only on UV observations.

Figure 7: A boxplot of the peak phases for the light curves in each band
relative to the $B$ band peak phase. Both the $UVW2$ and $U$ band light curves
have a narrow distribution of peak phases at around 2 days before the B band
peak phase, while the $V$ band light curves have a narrow distribution at
around 2 days after the $B$ band peak phase. The $UVW1$ band light curves seem
to peak earlier than the $B$ band light curves, but with a larger spread in
the distribution.

##  5 Conclusion

We modeled the light curves of 97 SNe Ia in the UV and optical from the Ultra-
Violet Optical Telescope on the Neil Gehrels Swift Observatory in order to
better understand the correlation of their properties in different wavelength
regimes. Modeling was done using a linear combination of band-vague templates,
based on functional principal component analysis (FPCA) of a sample of SNe Ia
in the optical bands. Each principal component is associated with a different
light curve property, making it a useful tool for finding correlations between
the filters.

We find that the overall dispersion in the PC1 and PC2 weights in all filters
were similar, suggesting that the UV light curves are not as variable as the
optical light curves. There is not a strong correlation in the PC1 weights for
the UV and optical filters, which would indicate the lack of a strong
correlation in the late-phase decline rate. There is a statistically strong
correlation between the UV and optical filters for the PC2 weights, which
would suggest that the stretch of the light curves and the post-peak decline
rates are correlated across the wavelength regimes. This could point to
underlying physics that ties the UV and optical flux together, with the bulk
of the light curve shape likely being driven by ${}^{56}Ni$ decay. The
differences that do exist are likely due to wavelength-dependent metallicity
effects and the opacity of the SNe Ia atmosphere. This is the first study to
statistically show that the UV and optical light curve properties have similar
variability, and that the stretch and decline rate of the UV and optical light
curves are significantly correlated.

##  6 Acknowledgements

The authors are grateful for the continuing support of the Mitchell Institute
for Fundamental Physics and Astrophysics. This work was supported by NASA
grant 80NSSC20K0456, “SOUSA’s Sequel: Improving Standard Candles by Improving
UV Calibration”. The authors would like to thank Minjee Kim, Jianhua, Huang,
and Xiaomeng Yan for lessons learned from related work. We also thank Nicholas
Suntzeff, Lauren Aldoroty, and Jiawen Yang for their assistance and expertise,
as well as the anonymous referee for their helpful comments. This work made
extensive use of the Python Programming language and the following community-
developed/maintained software packages: Jupyter (Kluyver et al., 2016), NumPy
(Harris et al., 2020), Matplotlib (Hunter, 2007), Pandas (Wes McKinney, 2010),
and SciPy (Virtanen et al., 2020). This work also used arXiv.org and NASA’s
Astrophysics Data System for bibliographic information.

## References

  * Akaike (1974) Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716, doi: 10.1109/TAC.1974.1100705
  * Amanullah et al. (2015) Amanullah, R., Johansson, J., Goobar, A., et al. 2015, MNRAS, 453, 3300, doi: 10.1093/mnras/stv1505
  * Branch (1998) Branch, D. 1998, ARAA, 36, 17, doi: 10.1146/annurev.astro.36.1.17
  * Breeveld et al. (2011) Breeveld, A. A., Landsman, W., Holland, S. T., et al. 2011, in American Institute of Physics Conference Series, Vol. 1358, GAMMA RAY BURSTS 2010. AIP Conference Proceedings, ed. J. E. McEnery, J. L. Racusin, & N. Gehrels, 373–376, doi: 10.1063/1.3621807
  * Brown & Forsythe (1974) Brown, M. B., & Forsythe, A. B. 1974, Journal of the American Statistical Association, 69, 364, doi: 10.1080/01621459.1974.10482955
  * Brown et al. (2017) Brown, P. J., Landez, N. J., Milne, P. A., & Stritzinger, M. D. 2017, ApJ, 836, 232, doi: 10.3847/1538-4357/aa5f5a
  * Brown et al. (2010) Brown, P. J., Roming, P. W. A., Milne, P., et al. 2010, ApJ, 721, 1608, doi: 10.1088/0004-637X/721/2/1608
  * Brown et al. (2012) Brown, P. J., Dawson, K. S., de Pasquale, M., et al. 2012, ApJ, 753, 22, doi: 10.1088/0004-637X/753/1/22
  * Brown et al. (2014) Brown, P. J., Kuin, P., Scalzo, R., et al. 2014, ApJ, 787, 29, doi: 10.1088/0004-637X/787/1/29
  * Copin et al. (2012) Copin, Y., Gangler, E., Pereira, R., et al. 2012, The Astronomer’s Telegram, 4253, 1
  * Garnavich et al. (2004) Garnavich, P. M., Bonanos, A. Z., Krisciunas, K., et al. 2004, ApJ, 613, 1120, doi: 10.1086/422986
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005, doi: 10.1086/422091
  * Goldhaber et al. (2001) Goldhaber, G., Groom, D. E., Kim, A., et al. 2001, ApJ, 558, 359, doi: 10.1086/322460
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * He et al. (2018) He, S., Wang, L., & Huang, J. Z. 2018, ApJ, 857, 110, doi: 10.3847/1538-4357/aab0a8
  * Höflich et al. (1996) Höflich, P., Khokhlov, A., Wheeler, J. C., et al. 1996, ApJL, 472, L81+, doi: 10.1086/310363
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jha et al. (2007) Jha, S., Riess, A. G., & Kirshner, R. P. 2007, ApJ, 659, 122, doi: 10.1086/512054
  * Jha et al. (2006) Jha, S., Kirshner, R. P., Challis, P., et al. 2006, AJ, 131, 527, doi: 10.1086/497989
  * Kamiya (2012) Kamiya, Y. 2012, Light Curve Model for SN 2011aa, Poster presented at Supernovae Illuminating the Universe: from Individuals to Populations, September 10–14, Garching bei Muenchen, Germany
  * Kasen & Woosley (2007) Kasen, D., & Woosley, S. E. 2007, ApJ, 656, 661, doi: 10.1086/510375
  * Kattner et al. (2012) Kattner, S., Leonard, D. C., Burns, C. R., et al. 2012, PASP, 124, 114, doi: 10.1086/664734
  * Kelly (2007) Kelly, B. C. 2007, ApJ, 665, 1489, doi: 10.1086/519947
  * Kenworthy et al. (2021) Kenworthy, W. D., Jones, D. O., Dai, M., et al. 2021, arXiv e-prints, arXiv:2104.07795.  https://arxiv.org/abs/2104.07795
  * Kessler et al. (2009) Kessler, R., Becker, A. C., Cinabro, D., et al. 2009, ApJs, 185, 32, doi: 10.1088/0067-0049/185/1/32
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Scmidt (IOS Press), 87–90.  https://eprints.soton.ac.uk/403913/
  * Kokoska & Zwilinger (2000) Kokoska, S., & Zwilinger, D. 2000, Chapman & Hall: New York, doi: 10.1201/b16923
  * Krisciunas et al. (2000) Krisciunas, K., Hastings, N. C., Loomis, K., et al. 2000, ApJ, 539, 658, doi: 10.1086/309263
  * Leibundgut (2001) Leibundgut, B. 2001, ARAA, 39, 67, doi: 10.1146/annurev.astro.39.1.67
  * Mandel et al. (2011) Mandel, K. S., Narayan, G., & Kirshner, R. P. 2011, ApJ, 731, 120, doi: 10.1088/0004-637X/731/2/120
  * Mazzali et al. (2001) Mazzali, P. A., Nomoto, K., Cappellaro, E., et al. 2001, ApJ, 547, 988, doi: 10.1086/318428
  * Mazzali et al. (2007) Mazzali, P. A., Röpke, F. K., Benetti, S., & Hillebrandt, W. 2007, Science, 315, 825, doi: 10.1126/science.1136259
  * Milne et al. (2010) Milne, P. A., Brown, P. J., Roming, P. W. A., et al. 2010, ApJ, 721, 1627, doi: 10.1088/0004-637X/721/2/1627
  * Nugent et al. (1995) Nugent, P., Phillips, M., Baron, E., Branch, D., & Hauschildt, P. 1995, ApJL, 455, L147+, doi: 10.1086/309846
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825.  http://jmlr.org/papers/v12/pedregosa11a.html
  * Perlmutter et al. (1999) Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565, doi: 10.1086/307221
  * Phillips (1993) Phillips, M. M. 1993, ApJL, 413, L105, doi: 10.1086/186970
  * Phillips et al. (1999) Phillips, M. M., Lira, P., Suntzeff, N. B., et al. 1999, AJ, 118, 1766, doi: 10.1086/301032
  * Poole et al. (2008) Poole, T. S., Breeveld, A. A., Page, M. J., et al. 2008, MNRAS, 383, 627, doi: 10.1111/j.1365-2966.2007.12563.x
  * Pskovskii (1977) Pskovskii, I. P. 1977, Soviet Astronomy, 21, 675
  * Riess et al. (1996) Riess, A. G., Press, W. H., & Kirshner, R. P. 1996, ApJ, 473, 88, doi: 10.1086/178129
  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009, doi: 10.1086/300499
  * Roming et al. (2005) Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, Space Science Reviews, 120, 95, doi: 10.1007/s11214-005-5095-4
  * Scalzo et al. (2014) Scalzo, R. A., Childress, M., Tucker, B., et al. 2014, MNRAS, 445, 30, doi: 10.1093/mnras/stu1723
  * Tukey (1977) Tukey, J. W. 1977, Exploratory Data Analysis (Addison-Wesley)
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * Wood-Vasey et al. (2008) Wood-Vasey, W. M., Friedman, A. S., Bloom, J. S., et al. 2008, ApJ, 689, 377, doi: 10.1086/592374

##  7 Appendix

Band | Interquartile Range | Brown-Forsythe Test Statistic | Brown-Forsythe P-Value
---|---|---|---
W2 | 2.062 | 0.169 | 0.682
M2 | 2.235 | 0.012 | 0.913
W1 | 1.913 | 0.112 | 0.738
U | 2.679 | 0.118 | 0.732
B | 2.320 | - | -
V | 2.597 | 1.516 | 0.221
Table 2: A comparison of the variability in the PC1 values as compared to the $B$ band. Similar IQR scores and P values $>$ 0.05 would suggest that there is not statistical difference in the variability between the bands. Band | Interquartile Range | Brown-Forsythe Test Statistic | Brown-Forsythe P-Value
---|---|---|---
W2 | 1.304 | 0.426 | 0.516
M2 | 1.731 | 3.118 | 0.082
W1 | 1.110 | 0.359 | 0.550
U | 1.095 | 0 | 0.983
B | 0.900 | - | -
V | 0.581 | 10.731 | 0.001
Table 3: A comparison of the variability in the PC2 values as compared to the $B$ band. Similar IQR scores and P values $>$ 0.05 would suggest that there is not statistical difference in the variability between the bands. Band | Interquartile Range | Brown-Forsythe Test Statistic | Brown-Forsythe P-Value
---|---|---|---
W2 | 0.861 | 7.804 | 0.007
M2 | 0.925 | 13.237 | 0.001
W1 | 1.444 | 28.489 | 0
U | 1.075 | 21.533 | 0
B | 0.448 | - | -
V | 0.481 | 0.062 | 0.804
Table 4: A comparison of the variability in the PC3 values as compared to the $B$ band. Similar IQR scores and P values $>$ 0.05 would suggest that there is not statistical difference in the variability between the bands. Band | Interquartile Range | Brown-Forsythe Test Statistic | Brown-Forsythe P-Value
---|---|---|---
W2 | 0.350 | 0.008 | 0.930
M2 | 0.473 | 0.107 | 0.745
W1 | 0.367 | 1.533 | 0.219
U | 0.337 | 0.031 | 0.861
B | 0.363 | - | -
V | 0.188 | 9.904 | 0.002
Table 5: A comparison of the variability in the $\Delta m_{15}$ values as
compared to the $B$ band. Similar IQR scores and P values $>$ 0.05 would
suggest that there is not statistical difference in the variability between
the bands. Figure 8: The correlation between the $B$ band and other Swift UVOT
bands for the PC1 weights. Only the $UVW2$, $U$, and $V$ bands have a
statistically significant correlation based on the Spearman rank-order
coefficient test. The red lines indicate the spread of the MCMC draws from
Linmix, and the black line is the mean linear fit from which we calulate the
$r^{2}$ value. The PC1 weights are most directly related to the decline rate
of the light curves. Figure 9: The correlation between the B band and the
other Swift UVOT bands for the PC2 weights. The statistical significance of
the correlation is based on the Spearman rank-order coefficient test. The red
lines indicate the spread of the MCMC draws from Linmix, and the black line is
the mean linear fit from which we calulate the $r^{2}$ value. The PC2 weights
are most directly related to the spread of the light curves. Figure 10: The
correlation between the B band and the other Swift UVOT bands for the PC3
weights. Only the $U$ and $V$ bands have a statistically significant
correlation based on the Spearman rank-order coefficient test. The red lines
indicate the spread of the MCMC draws from Linmix, and the black line is the
mean linear fit from which we calulate the $r^{2}$ value. The PC3 weights are
most directly related to the pre-peak rise in the light curve and late phase
adjustments. Figure 11: The correlation between the B band and the other Swift
UVOT bands for the 15 day decline rate post-peak. The red lines indicate the
spread of the MCMC draws from Linmix, and the black line is the mean linear
fit from which we calulate the $r^{2}$ value. The decline rate is calculated
from the best fit FPCA model for each individual light curve.