# Photonic Spiking Neural Networks with Highly Efficient Training Protocols
for Ultrafast Neuromorphic Computing Systems

Dafydd Owen-Newns  Joshua Robertson  Matěj Hejda  Antonio Hurtado

###### Abstract

Photonic technologies offer great prospects for novel ultrafast, energy-
efficient and hardware-friendly neuromorphic (brain-like) computing platforms.
Moreover, neuromorphic photonic approaches based upon ubiquitous, technology-
mature and low-cost Vertical-Cavity Surface Emitting Lasers (VCSELs) (devices
found in fibre-optic transmitters, mobile phones, automotive sensors, etc.)
are of particular interest. Given VCSELs have shown the ability to realise
neuronal optical spiking responses (at ultrafast GHz rates), their use for
spike-based information processing systems has been proposed. In this work,
Spiking Neural Network (SNN) operation, based on a hardware-friendly photonic
system of just one Vertical Cavity Surface Emitting Laser (VCSEL), is reported
alongside a novel binary weight ’significance’ training scheme that fully
capitalises on the discrete nature of the optical spikes used by the SNN to
process input information. The VCSEL-based photonic SNN is tested with a
highly complex, multivariate, classification task (MADELON) before performance
is compared using a traditional least-squares training method and the
alternative novel binary weighting scheme. Excellent classification accuracies
of $>$ 94$\%$ are reached by both training methods, exceeding the benchmark
performance of the dataset in a fraction of processing time. The newly
reported training scheme also dramatically reduces training set size
requirements as well as the number of trained nodes ($\leq$ 1$\%$ of the total
network node count). This VCSEL-based photonic SNN, in combination with the
reported ’significance’ weighting scheme, therefore grants ultrafast spike-
based optical processing with highly reduced training requirements and
hardware complexity for potential application in future neuromorphic systems
and artificial intelligence applications.

##  1 Introduction

###  1.1 Photonic ANNs and Spiking VCSELs

Artificial Neural Networks (NNs) are more frequently appearing in the fast-
developing field of Artificial Intelligence (AI) thanks to their proven
ability to provide high level performance at numerous complex information
processing tasks [1, 2] such as computer vision and data classification. NNs,
that draw inspiration from the networks of biological neurons in the brain,
are formed of highly parallel structures of nodes (neurons) that realize
numerous non-linear transformations to achieve efficient processing and
decision making. The desire to move towards a more-than-Moore era of computing
and diverge from traditional computing paradigms has also spurred on reports
of brain-inspired computing, where now multiple systems of large-scale
electronic ANN systems, such as the Loihi [3], TrueNorth [4] and BrainScaleS
[5], amongst others, exist. However, recently the number of works
demonstrating optical approaches to beyond-Von Neumann computing is also
growing, given the appeal of the inherent properties of the optical medium.
Optical-based systems can operate with very low component crosstalk, perform
across large data bandwidth ranges, and utilise non-interacting frequencies of
lights to achieve highly parallel division-multiplexing applications [6]. More
importantly however, the photonic platform can allow NNs to operate with low
power consumption and ultrafast speed, overcoming the squeeze of power
efficiency and clock speed faced by state of the art chip-scale electronic
components [7, 8]. The realisation of light-based NNs, or so called Photonic
Neural Network (PNN) systems, is therefore of key interest to all future AI
applications where the benefits of increased operation speed and power
efficiency is directly felt in the training and processing of large data
volumes.

Despite the relative infancy of brain-inspired (neuromorphic) photonic
systems, the field is advancing quickly with the attractive properties of
optics already inspiring the production of multiple PNN accelerators [9, 10,
11, 12, 13, 14, 15, 16]. These system, based on technologies such as phase
change materials [9, 10], optical modulators [11, 16], and micro-ring weight
banks [12, 13, 14, 15], amongst others, have risen to demonstrate different
information processing tasks with photonic components. Neuromorphic photonic
systems, built with semiconductor laser (SL) technologies, have also
demonstrated intriguing neuron-like non-linear dynamics key to the operation
of PNN systems [17, 18]. Vertical-Cavity Surface-Emitting lasers (VCSELs),
well-established, commercially available and increasingly ubiquitous devices
(found in mobile phones, data centres, supermarket barcode scanners, etc.),
are one such neuromorphic SL technology. The ability of VCSELs to imitate the
behaviours of biological neurons was first proposed in [19], where it was
shown that controllable spike activation could be achieved with incoming
optical signals by exploiting non-linear dynamics surrounding the injection
locking condition. The neuron-like spiking responses were achieved at
ultrafast (GHz) rates with spike widths of approximately 100 ps, making the
spiking responses $>6$ orders of magnitude faster than biological neurons, and
multiple orders of magnitude faster than some electronic spiking systems.
VCSELs have since showcased several common neuronal behaviours (such as
integrate and fire spiking [20, 21], refractoriness and rate coding [22]) in
experimental realisations. Further, early demonstrations of neuromorphic
information processing functionalities with VCSELs has been achieved using the
spiking neuronal behaviours of VCSELs both experimentally, including image
processing (edge feature detection) [23, 24, 25], pattern recognition [21] and
exclusive OR (XOR) operation [26], as well as theoretically through
simulations of the Yamada and Spin-flip models (see [27] for a review).
Moreover, single VCSEL devices have shown processing functionality when
operated in combination with software-based Spiking Neural Networks (SNNs),
achieving high classification accuracy at the MNIST hand written digit
recognition task [28]. VCSELs therefore represent an increasingly promising
technological approach to spike-based PNNs that operate with ultrafast,
telecommunication wavelength, optical signals, that is hardware friendly, low
power and low cost.

Beyond the photonic imitation of spiking neuronal behaviours, the technique
reservoir computing (RC) has been demonstrated as a powerful method for
creating PNNs, yielding excellent performance in complex tasks while also
benefiting from highly hardware-friendly architectures. First developed in the
early 2000’s [29, 30], it was shown that in RC architectures (NN architectures
with fixed-weight hidden layer connections) only the output layer required
training to achieve successful performance. In these RC architectures, the
nodes within the hidden layer are referred to as the ”reservoir”. The
reservoir is therefore formed of unknown, fixed-strength, interconnected non-
linear elements that are coupled to an output layer. Given the fixed nature of
the connections, RC limits the training requirements of NNs, helping
significantly reduce the resources (computational power and time) needed to
train large networks of nodes for successful operation [29]. There have been
numerous theoretical [31, 32] and experimental [33, 34, 35] reports of
reservoir computers based on the photonics platform with VCSELs again
featuring frequently [36, 37, 38, 39, 40, 41]. VCSELs are one of the photonics
devices where both time-delay reservoirs (TDRs) [38, 39, 40, 41], which
multiplex nodes over time to construct virtual networks and create memory
through feedback connections, and spatial temporal reservoirs [42, 43], which
spatially multiplex nodes, have demonstrated successful operation. These
systems have revealed impressive performance on numerous complex processing
tasks, utilising both off-the-shelf [41], and bespoke large area (LA)-VCSEL
designs [42], to realise photonic computing with continuous (non-spiking)
VCSEL signals (see [44] for a review). More recently, we demonstrated for the
first time the combination of the reservoir computing technique with a
neuromorphic spiking VCSEL-neuron [45]. In that report we revealed that by
using the time-multiplexing technique and the masking of inputs, we could
interpret a single VCSEL neuron as an fixed-weight interconnected virtual SNN.
This photonic system realised both the reduced training requirements of an RC
system as well as the sparse, all-or-nothing, binary spike-based
representation of a full SNN. The VCSEL-based SNN operated with short temporal
nodes (250 ps), fast optical spikes (100 ps-long), and a configurable network
node number, that allowed it to demonstrate successful operation on the
benchmark Iris Flower classification task [46] with excellent performance.

In this work, we show that we can not only use time-multiplexing to create a
VCSEL-based photonic SNN capable of a highly complex processing task, but more
importantly we can, through alternative training schemes, achieve high
performance with significantly reduced training requirements (namely reduced
training set sizes and number of trained nodes). This paper adopts the
following structure. Firstly, in Section 1.2 we will introduce the
experimental setup used in this work to create the photonic SNN and discuss
its operating principles. Next, in Section 1.3 we will introduce the
experimental classification task and explain how data is prepared for
injection into the VCSEL-based SNN. Then, in Section 2.1 we will discuss the
results of the classification task and explore training through the standard
linear least-square regression method. In Section 2.2 we will then introduce
the alternative binary weight training scheme and discuss its performance
relative for various training set sizes and node requirements. Finally, we
will provide conclusions in Section 3.

Figure 1: a) The experimental setup used to create the photonic SNN with a
single VCSEL device. b) The lasing threshold (Light-Current) characteristics
of the experimental VCSEL-neuron reveals a threshold current of 1.423 mA. c)
The spectra of the experimental VCSEL-neuron at different bias currents. At
the operating current (3.5 mA - red) the device exhibited two orthogonal
polarisation peaks at 1287.44 nm and 1287.59 nm. Measurements were taken with
a room temperature (293 K) stabilised device. Figure 2: a) The photonic SNN
operating principle. Masked input data is injected into the VCSEL forming the
photonic SNN. The input information is time-multiplexed with each $\theta$
time-slot representing a virtual node (neuron) in the SNN. The time-
multiplexed output of the VCSEL is interpreted as a binary node output, either
spiking (1) or non-spiking (0). The output layer nodes are weighted offline to
provide the photonic SNN readout. b) The architecture of the photonic SNN. 500
features per datapoint are input into the SNN system. In this experimental
demonstration the hidden layer was selected to be made up of either 2048 or
4096 virtual nodes. The output layer is trained offline to provide a datapoint
classification (either +1 or -1). c) The MADELON classification task.
Datapoints are placed on the corners of a 5-dimensional hypercube and randomly
labelled -1 or +1. Each artificial data point has 5 informative features (each
corresponding to 1 of 5 dimensions). 15 linear combinations of the original
features are made creating 20 total features. 480 non-informative features
(that carry no information of the label) are then added to the feature set.
The task requires all 500 data features to be used to classify the initial
data point label -1 or +1.

###  1.2 Photonic SNN Setup and Operation

In this work, we report a hardware-friendly photonic SNN based on a single
time-multiplexed VCSEL-neuron. The experimental setup used to build and
investigate the photonic SNN architecture of this work is shown graphically in
Fig. 1 (a). The photonic SNN is built with a single telecom-wavelength VCSEL
and commercially-sourced fibre-optic telecom components. The selected VCSEL
had wavelength emission within the key optical telecommunication O-band window
(centred around 1300 nm) to showcase full compatibility of the system with
optical networking and data centre technology. In the experimental setup,
external light injection is provided by a tuneable laser (TL) source which is
intensity modulated by a 10 GHz Mach Zehnder modulator (MZ). The MZ is
responsible for encoding the input information in the intensity of the optical
injection before it enters the VCSEL-neuron. The MZ is driven by a DC power
supply (PS) and a fast 12 GSa/s 5 GHz arbitrary waveform generator (AWG) that
generates the required input information. A 10 dB electrical amplifier (AMP)
is used to increase the signal from the AWG, and a polarisation controller
(PC) is used to maximise TL light coupling into the MZ. The modulated
injection is split by a fibre coupler to make a measurement of injection power
before the TL light enters the circulator and the VCSEL. The output of the
VCSEL is captured by the circulator and fed to analysis, where it was measured
by an optical spectrum analyser (OSA) and a fast 40 GSa/s 16 GHz real-time
oscilloscope (OSC). The threshold and spectral characteristics of the O-band
VCSEL used in the photonic SNN system are shown in Figs 1 (b) & (c). When
stabilised at room temperature (293 K), the VCSEL exhibited a threshold
current of 1.423 mA and produced two orthogonal polarisation peaks in the
fundamental transverse mode. In this work, with the selected operating current
of 3.5 mA, these orthogonal polarisation peaks occurred at 1287.44 nm and
1287.59 nm, producing a total output power of $\sim$220 $\mu$W. Using a PC,
optical injection was polarisation-matched to the subsidiary (1287.59 nm)
polarisation mode of the VCSEL, with injection made using a negative frequency
detuning to induce polarisation switching and injection locking. The
modulating intensity is responsible for the triggering of high speed (approx.
100 ps-long) neuron-like spiking responses near the injection locking
boundary.

In this work, as in [45], we combine the concepts of reservoir computing with
the photonic spiking VCSEL-neuron to create a fully photonic SNN using
significantly reduced hardware (a single laser device). As discussed in
Section 1.1, reservoir computers are a type of artificial NN that host a
hidden layer of interconnected nodes whose connections and weights are fixed.
This means that in reservoir computers only the output layer of the network
needs to be trained to achieve high performance [29]. Further, reservoir
computers can make use of time-multiplexing to interpret the output of one
node as the output of many nodes, by sampling at set times ($\theta$). Here,
using the VCSEL-neuron, we apply the same concepts to interpret the output of
one spiking VCSEL as an entire photonic SNN, see Fig 2 (a). The output of the
VCSEL is time-multiplexed at intervals of $\theta$ to produce a controllable
number of network nodes ($N_{v}$). The coupling of network nodes is achieved
given the duration of the nodes ($\theta$) is less than the timescale of the
neuron-like spiking non-linearity in the VCSEL (typically $\sim$ 1 ns [22]).
This condition allows the neuronal leaky integrate-and-fire and refractory
behaviours of the VCSEL-neuron to link multiple nodes coupled in time, forming
a complex interconnected virtual network structure. As in [45] (and those of
other VCSEL-based reservoir computers [41]), the input data is prepared by
randomly masking data point features and injecting them into the VCSEL as
continuously modulated light. Following injection, each $\theta$-long node
will deliver an all-or-nothing binary output, determined by the presence (or
absence) of a fast 150 ps optical spike, effectively realising the output of a
fast photonic SNN. In this system the output weights (the weights applied to
each node) are trained and calculated offline to provide the readout of the
system for the prescribed processing task (the task and training schemes are
discussed in Sections 1.3 & 2).

In [45], we recently demonstrated experimentally that a VCSEL-neuron system
can successfully realise SNN operation with fast optical spiking signals.
There we showcased the benchmark nonlinear Iris flower classification task,
comprising of 150 non-linearly separable data points (flowers), 4 features per
data point (sepal & petal length/height) and 3 flower classes. Yet, using our
VCSEL-based photonic SNN we demonstrated that a very high overall average
classification accuracy of $>$97 % could be achieved, showcasing the photonic
SNNs powerful flexibility and performance. Importantly, in that work the
training of the output weights was achieved using an ordinary least squares
regression method, and hence the node outputs were treated as if it were
continuous (spiking nodes were assigned the value 1.0, and others given 0.0).
This result required the offline calculation of all node weights as real
number (floating point) values. In this report we focus now on completing a
much more complex problem and apply an alternative training scheme, whereby
the discrete nature of the binary all-or-nothing spiking responses can be used
to improve the speed and efficiency of training. In the following sections we
will introduce the new, high-complexity, non-linear classification task
(MADELON) [47] and demonstrate successful performance of the photonic SNN with
two training methods; the aforementioned least-square regression training
method, and an alternative node-significance training method, to highlight the
attainment of notable improvements in both training speed and resource
requirements.

###  1.3 Data Set Preparation

We further demonstrate the capability of the VCSEL-based photonic SNN by
applying the system to a highly complex classification task. The task used in
this work is an artificial dataset named ’MADELON’, which was created to test
feature selection methods as part of a NIPS 2003 challenge [47]. A simplified
schematic of the MADELON data creation process is shown in Fig. 2 (c). The
dataset is made up of 32 clusters of datapoints placed on the vertices of a
5-dimensional (5D) hypercube. Two classes of data are created by randomly
labelling each datapoint with the value -1 or +1. Each data point has 5
informative features (one for each dimension) that can be used to identify the
class of the data. Using the 5 informative features, 15 linear combinations
are created, increasing the feature set size to 20. Finally, the 20 features
are randomly scattered between 480 non-informative (no predictive power)
features creating a total of 500 features per data point. The task is
therefore a multivariate, nonlinear, two-class problem with continuous input
variables. The task is of high complexity given a classification must be made
using a feature set where no feature by itself is informative. The problem is
further complicated by the addition of non-informative features, which require
the system to eliminate and ignore distractions.

The entire MADELON dataset contains 2000 data points intended for training.
However, here to demonstrate the operation of the VCSEL-based photonic SNN, we
used 300 data points in total (150 of each class). Limiting the number of
input points reduced the length of the time-multiplexed output generated
during the experimental runs, hence allowing us to test different network
architectures with diverse node numbers while remaining within the memory of
our experimental equipment. As the MADELON task is of much higher complexity
than the preliminary Iris Flower classification task [45], a higher number of
virtual nodes ($N_{v}$) were selected, namely 2048 and 4096, to test the
classification performance of the photonic SNN.

The photonic SNN architecture used to complete the task is shown in Fig. 2
(b). Each of the 500 features are fed into the $N_{v}$ input layer nodes via a
random mask. The masking of the input data is achieved by multiplying a column
vector of all 500 features with a 500 x $N_{v}$ random matrix, resulting in a
vector of $N_{v}$ components. The masked data is fed into the hidden layer of
the photonic SNN, where fixed connections and weights are created by the non-
linear temporal dynamical behaviour in the VCSEL. Experimentally, the randomly
masked data is generated by the AWG and modulated by the MZ into the optical
injection of the VCSEL. The output layer of $N_{v}$ nodes is then provided by
the segmentation of the VCSEL’s spiking output timeseries at $\theta$
intervals. In this work a node separation ($\theta$) of 250 ps was used to
time-multiplex network nodes, creating multiple interacting nodes within the 1
ns timescale of non-linearity in the VCSEL [22, 21]. The final classification
layer of the SNN then provides a -1 or 1 label to the data following the
offline training and weighting of output nodes. The following section will
discuss the two training methods implemented in this work.

##  2 Results

First, the MADELON task was run in the photonic VCSEL-based SNN using a
network node count of $N_{v}=2048$. This results in a total processing time of
512 ns per datapoint, in the photonic SNN. The VCSEL was driven with a current
of 3.5 mA ($\sim$ 2.5 times the lasing threshold current), at a temperature of
293 K; the injected light had mean power of 142.6 $\mu$W, and was injected
into the VCSEL with an initial frequency detuning of -5.4 GHz with respect to
the resonant wavelength of the VCSEL’s subsidiary polarisation mode peak.
These conditions put the VCSEL into a state in which it was injection-locked
to the externally-injected signal. Here, near the injection locking boundary,
optical spikes could be triggered by the varying intensity of the modulated
optical injection light.

On a second experimental run, the number of virtual (spiking) nodes in the
network was increased to $N_{v}=4096$ (1024 ns processing time per datapoint).
In this second run the VCSEL was biased with a current equal to 3.45 mA, at a
set temperature of 293 K. The injected light signal, encoded with the input
data, had a mean optical power of 161 $\mu$W and was injected with an initial
frequency detuning of -4.6 GHz (with respect to the resonance of the VCSEL’s
subsidiary polarisation mode).

The output of the photonic SNN, for a given datapoint, is a vector $S$ of
length $N_{v}+8$. The 8 additional virtual nodes are the result of zero
padding at the end of the input data sequence, equivalent to a 2 ns reset of
the system between inputs. Resetting the system between consecutive datapoint
inputs prevented the previous datapoint from influencing the new datapoint and
generating undesired spike activations. In the SNN output vector $S$, the ith
element was readout as 1 if the ith node spiked within the 250 ps-long time-
multiplexed segment, and was assigned value 0 otherwise.

Figure 3 shows experimentally-measured time-series depicting the optical input
(Fig. 3 (a-b)) and output signals (Fig. 3 (c)) of the photonic SNN. The
example shown in Fig. 3 is for a network architecture with 2048 nodes. Time-
multiplexing is used to inject each masked datapoint into all the (virtual)
nodes of the photonic SNN (at a rate of 250 ps per virtual node). This creates
a time varying signal for each datapoint, which is applied as an intensity
modulation of the light injected into the SNN. During the experiment, all 300
masked datapoints (time varying signals) are injected into the SNN
sequentially, creating a large optical input timeseries. Figure 3 (a) depicts
part of the large optical input time-series, showing 50 sequentially encoded
masked datapoints of the MADELON dataset. Figure 3 (b) plots the first 512 ns
of the optical input time-series in Fig. 3 (a), which corresponds to the input
signal encoding the first MADELON datapoint. As the amplitude of the optical
input signal varies, it will trigger the firing of optical spike events from
the VCSEL at specific times (virtual nodes). Therefore, different input
signals (masked datapoints) entering the photonic SNN will elicit different
optical spike trains from the network. The optical spiking time-series,
measured at the output of the SNN in response to the first injected MADELON
dataset point (input signal in Fig. 3b), is shown in Fig. 3 (c).

Figure 4 plots a 2D temporal map that collects in one plot the spiking output
from all virtual nodes of the photonic SNN for all 300 consecutively injected
data points in the MADELON task. In Fig. 4 a yellow (blue) dot indicates a
spiking (non-spiking) response from the VCSEL-based SNN, and the resulting 1
(0) binary output of the virtual node in vector $S$. The output vector $S$ is
plotted horizontally for each input datapoint, revealing trends in node
activation across inputs and classes. The results for photonic SNNs configured
with 2048, and 4096, virtual (spiking) nodes are plotted in Figures 4 (a) and
(b), respectively. The red line through the centre of the map separates the
results for each distinct class (-1 & +1) of input datapoint. Figures 4 (a)
and (b) show that each input datapoint of the same class elicits a similar
spike train at the output of the photonic SNN. These spike trains are
sufficiently distinct from those obtained by the other class of datapoint, to
allow for a classification operation.

Figure 3: Time series showing the measured input (blue) and output (red) of
the photonic SNN for a network architecture of 2048 virtual nodes. a) A
section of injection waveform showing the time-multiplexed intensity
variations of 50 consecutive masked datapoints of the MADELON dataset. b) A
512 ns segment of the input timeseries corresponding to the first MADELON
datapoint. Each of the 2048 virtual nodes have a duration of 250 ps, creating
the 512 ns input timeseries. c) The corresponding measurement of the SNN
output when injected with the first MADELON datapoint. Figure 4: Temporal maps
showing the optical spiking patterns produced by the SNN in response to the
input dataset using a) 2048 virtual nodes and b) 4096 virtual nodes.

For the case of the SNN with a 2048 node architecture (Fig. 4 (a)), a total of
298 datapoints (149 per class) were used in the experimental run. The reason
for this is that the first datapoint in class $-1$ was affected by an
experimental artifact, and was hence remove from analysis. One point from
class $+1$ was therefore also removed to keep the size of each set equal,
resulting in a total of 298 points. For the network architecture with
$Nv=4096$ (Fig. 4 (b)), all 300 data points (150 per class) of the MADELON
task were used to test the network operation.

Training is completed using a randomly selected subset of the input
datapoints, and their corresponding spiking responses. The randomly selected
subset contains an equal number of points ($N_{t}$) of each class, for a total
training set size of $2N_{t}$. In our photonic SNN, the training is used to
calculate the weights applied to the output layer nodes (i.e. the temporal
spiking patterns produced by the network and the corresponding vector $S$).
The output weights are determined in such a way to achieve the best
performance on this set, with the assumption that the training set is a good
representation of the entire task dataset. The input datapoints unused in the
training process are subsequently used to test the performance of the MADELON
task. The output layer weights calculated during training are applied to the
output layer to infer the class of a given test datapoint. The inferred
datapoint class is then compared to the true class label and performance is
measured using the classification accuracy, the fraction of correctly
classified datapoints across all tested datapoints. The training/testing
process was repeated several times, randomly selecting the training set each
time in order to find an average value for the performance (a process known as
random cross-validation).

###  2.1 Ordinary Least Squares

In this work, we first used the Ordinary-Least-Squares (OLS) training method
to test the performance of the photonic SNN in the MADELON task, a standard
training method applied in other photonic RC systems [48, 45]. Here, the
output spiking patterns from the VCSEL (and corresponding vector $S$) are
treated as a binary sequence of floating-point 1.0s and 0.0s. Then, the output
layer weights are found via linear regression where the values of each node
(1.0 or 0.0) are arranged into matrix $X$ (where the rows of $X$ are the
output vectors $S$). The labels (expected values) of each datapoint form a
second column vector $Y$. The weight calculation then consists of finding the
weight matrix $W$ that minimises Equation (1):

$E_{rror}(W)=\left|Y-XW\right|^{2}.$ (1)

Equation (1) is solved using $W=\left(X^{T}X\right)^{-1}X^{T}Y$ (the Penrose
inverse of $X$ multiplied by $Y$). For this task there are two output classes
(-1 & +1), so the labels are two-element vectors, in which the nth element is
set to $1$ to denote the nth class, with the other elements being zero. This
also means the resulting weight matrix $W$ has dimension of $\left(N_{v}\times
2\right)$. Using $W$, the prediction for the class of a datapoint (produces a
spike pattern $S$) can be found by taking $argmax(S*W)$. Applying the weights
to the spike pattern vector gives a score towards each class label, where the
class with the highest score is the inferred label (prediction) of the
photonic SNN system.

Using only 15 training datapoints per class (10 $\%$ of the total 300 input
datapoints), the OLS training method achieved very high performance with a
peak accuracy of $91\%$ (for a 2048 node architecture) and $94.5\%$ (for a
4096 node architecture). Example confusion matrices for $N_{v}=2048$ and
$N_{v}=4096$, trained using this method and $N_{t}=15$ are shown in Fig. 5. As
expected, the network architecture with the higher node number ($N_{v}$)
provide the better classification performance overall. The larger number of
nodes creates more spiking responses from the VCSEL, allowing for more class-
specific responses and a better classification. The capability to adjust node
number at will is an inherent advantage of the photonic SNN of this work, as
no physical/hardware changes are required to increase the network
architecture. Instead, all that is require is the alteration of the input mask
and a longer timeseries measurement. This makes the reported VCSEL-based
photonic SNN an attractive and flexible system for the fast, optical spike-
based, implementation of numerous processing tasks.

Figure 5: Confusion matrices depicting the performance of the photonic SNN
when trained using the ordinary least squares method. A training set size of
15 ($N_{t}=15$) was used to achieved accuracies of 91 $\%$ and 94.4 $\%$ for
architectures of 2048 and 4096 nodes, respectively.

###  2.2 Binary Weight Training

As discussed in the previous section (Section 2.1) the standard OLS training
method can deliver high classification performance during the MADELON task.
However, the OLS method does not take full advantage of the discrete binary
nature of the fast optical spiking responses in the photonic SNN. For example,
it can be observed that in the ideal case of an SNN response, spiking patterns
would be perfectly consistent between data of the same class, as well as
completely distinct from data of alternative classes. In this ideal case,
$2^{k}$ classes would only require the use of $k$ nodes to be complete
separate and identify datapoints. In the ideal case of the MADELON task this
would mean only 1 node per class is required for identification. In the case
of the experimental data, initally identifying nodes that have more predictive
power (i.e. spike only for a particular class) prior to training, would
therefore allow for creation of weight matrices that retain high
classification accuracy, despite the use and calculation of only a few nodes
(a fraction of the total node count).

We introduce the following process for finding weights for this SNN:

  1. 1.

Count the number of spikes that occurred for each class in each node, denoted
$s_{n,i}$ ($\\#$ spikes from data of class $i$ occurring in node $n$).

  2. 2.

Calculate the significance score of each node towards each class according to
Equation (2):

$z_{n,i}=s_{n,i}^{2}/\sum_{i}s_{n,i}$ (2)

  3. 3.

For each class, choose the top $N_{n}$ scoring nodes and set their weight
towards that class ($W_{n,i}$) to 1, and all remaining node weights to zero.

Applying the aforementioned algorithm provides the ’significance scores’ for
all network nodes, with higher scores obtained for nodes that contain not only
more spikes, but nodes that have their spikes more concentrated in one
specific class. The resulting output weight matrix $W$ will only have nonzero
entries on nodes that are strongly indicative of one specific class. The 2D
temporal maps shown in Fig. 6 reveal the spike trains of a subset of network
nodes, for the case when the photonic SNN was operated with a total node count
of 2048. This plot demonstrates the significance scoring system applied to
different selected nodes. For example, Fig. 6 shows that node 844 (marked with
an amber arrow) elicited an optical spike for 57 of the total 298 input
datapoints used in this experimental run. In this case, the spiking responses
from $n=844$ are evenly distributed across each class (31 spikes for class -1,
and 26 for class +1). For this specific case, the significance scores for this
node are $z_{844,-1}=16.9$ (class -1) and $z_{844,+1}=11.9$ (class +1). In
parallel, Fig. 6 shows that Node 885 (marked with a red arrow) fires fewer
optical spikes overall. In total 30 responses are elicited, however because
they are less evenly distributed (and hence more indicative of a particular
class), the node scores $z_{885,-1}=2.1$ (class -1) and $z_{885,+1}=16.1$
(class +1), a class +1 significance score similar to node 844. Finally, when
the appearance of fast optical spikes is much more biased towards one specific
class, and occur more likely overall, the resulting significance score $z$ is
much higher. For example, in node 853 (marked with a green arrow) a total of
103 spikes occurred, with 88 of those occurring for class -1, yielding scores
of $z_{853,-1}=75.1$ and $z_{853,+1}=2.1$. With this method, the weight
$W_{853,-1}$ of node 853 would be set to 1, as if a spike is detected in node
853, it is likely that the datapoint is of class -1. In this example the other
highlighted nodes would have their weights ($W_{844,-1}$ & $W_{885,-1}$) set
to zero.

For a data point to be classified the spike output vector $S$ must be
multiplied by the weight matrix $W_{i}$ of each class. The class with the
highest sum value is the prediction for the class of the data point. Since the
output layer weights are binary, this is equivalent to applying a filter to
the optical spike trains and counting the number of spikes remaining. Using
this method, the performance can be tuned by choosing the number of nodes
($N_{n}$) that are selected to be used for training (training node number).
Using a fraction of the total number of available nodes can result in equal
(if not better) performance, while keeping the weight matrix and the number
weight calculations, sparse.

Figure 6: Section of a 2D temporal map for a network architecture of 2048
nodes (between nodes 842 and 890). The arrows mark nodes 844 (orange), 853
(green) and 885 (red). In node 853, 88 spikes occur for class -1, and 15 for
class +1, giving it a high significance score of 75.1 (for class -1) and 2.1
(for class +1). Node 844 has significance scores of 16.9 (for class -1) and
11.9 (for class +1). Node 885 has significance scores of 2.1 (for class -1)
and 16.1 (for class 1). Node 853 has the highest significance score meaning it
would likely be used to perform a classification. Figure 7: Confusion matrices
depicting the performance of the photonic SNN when trained using the
significance score method. A training set size of 15 ($N_{t}=15$) was used to
achieved accuracies of 92.5 $\%$ and 95.2 $\%$ for architectures of 2048 and
4096 nodes, respectively. Only the 20 most significant nodes ($N_{n}=20$) were
used to train the SNN.

Figure 7 shows the resulting confusion matrices of the MADELON task when the
significance training algorithm is run with 15 training datapoints
($N_{t}=15$). For the photonic SNN architecture with 2048 nodes a peak
accuracy of 92.5 $\%$ was obtained and for the 4096 node architecture a peak
accuracy of 95.2 $\%$ was achieved. Again, as expected the larger node number
architecture provided the higher classification performance. Comparing this
performance to the OLS training method, also run at $N_{t}=15$, we find the
significance training shows small improvement on peak accuracy in both cases
(2048 and 4096 node architectures). This improvement to performance shows that
despite the same training set size, by considering only nodes significant to
specific classes, we can achieve similar or higher classification performance.
Furthermore, it is important to highlight that the training node number used
in the significance scoring method was substantially lower than that used in
the OLS method. In the novel training method only the top 20 most significant
nodes ($N_{n}=20$) were used in both the reported SNN architecture cases (2048
and 4096). The novel training method therefore only uses 1$\%$ and 0.5$\%$
(for 2048 and 4096 nodes respectively) of the total nodes used for OLS
training where every node in the network is trained ($N_{n}=N_{v}$). This
result indicates that training less nodes overall is not only more efficient
computationally, but does not directly hinder the performance of the system
(in this case increase classification performance). These newly presented
classification performance values are not however the limit of the photonic
SNN with this training scheme. In Figs. 8 and 9, we consider a larger range of
training set sizes ($N_{t}$), and compare their optimal performance as well as
their optimal training node number ($N_{n}$).

Figures 8 (a) and 9 (a) show the peak accuracy achieved in the MADELON
classification task as a function of the training set size used, for both 2048
and 4096 network architectures, respectively. Figures 8 (b) and 9 (b) show in
turn the optimal training node number ($N_{n}$) used to attain the shown
optimal accuracy.

Both Figs. 8 and 9 reveal two key results. First, the photonic SNN achieves
very high accuracy levels despite the complexity of the MADELON task and its
very large number of features. Specifically, a maximum accuracy of 94.4$\%$ is
achieved for the case of a 2048 node architecture, and a maximum accuracy of
95.7$\%$ is achieved for the case of a 4096 node architecture. The performance
of this photonic SNN system is in fact higher than, the results of previous
software-based machine learning algorithms reported following the NIPS 2003
challenge [47], which reached accuracy levels of up to 93.78 $\%$.
Importantly, one key difference here is that a photonic SNN is applied to the
MADELON task with a very hardware friendly implementation (using just one
VCSEL), low-power operation (sub-pJ energy per spike, $\sim$ 150 $\mu$W
average optical powers, and 3.5 mA of applied bias current), and ultrafast
performance (250 ps/node yielding a total processing time of 512 ns and 1024
ns per datapoint, for 2048 and 4096 nodes respectively). Secondly, the
performance of the photonic SNN increases with the training set size
(associated with best-fit training), and the optimal value of $N_{n}$ may
slightly decrease with training set size. This result is somewhat expected as
training with more datapoints typically improves system performance, resulting
in a lower dependence on high node numbers. Thirdly, Figs. 8 and 9 show that
the optimal training node number ($N_{n}$) remains consistently low while the
system retains very high performance, across all training set sizes ($N_{t}$)
(from 1 to 100). Only small node numbers, as low as $<$ 10 nodes (out of the
total node count $N_{v}$) need to be considered to successfully achieve high
performance in the MADELON task. This means that independently of the training
set size the high reduced training requirements of the significance scoring
method remain a key benefit to the alternative training scheme. Finally,
according to Figs. 8 and 9, very small training set sizes, $N_{t}\leq\,10$
($\leq$ 3 $\%$ of the total dataset inputs), are capable of effectively
training the photonic SNN system and achieving high accuracy ($>90\%$) in the
MADELON task. Comparing the previous OLS results for a network architecture of
4096 (from 2.1) to those of Fig. 9, we find that at training set sizes $<$ 15,
the significance scoring method achieves similar levels of classification
accuracy. Making use of the novel training approach can therefore enable
similar performance to the OLS method but with the remarkable benefit of being
more flexible in the number of training datapoints.

Figure 8: Photonic SNN performance versus the training set size, using a 2048
node architecture. a) The optimal classification accuracy (averaged over 10
random selections of the training set) is plotted against training set size
($N_{t}$, blue). The inset highlights the consistently high performance over
various training set sizes ($N_{t}$ ranges from 1 to 100). b) The number of
training nodes ($N_{n}$) used to achieve the corresponding optimal accuracy,
plotted against increasing training set size (red) . The peak accuracy,
94.4$\%$, occurs when $N_{n}=6$ and $N_{t}=76$ (circled). Figure 9: Photonic
SNN performance versus the training set size, using a 4096 node architecture.
a) The optimal classification accuracy (averaged over 10 random selections of
the training set) is plotted against training set size ($N_{t}$, blue). The
inset highlights the consistently high performance over various training set
sizes ($N_{t}$ ranges from 1 to 100). b) The number of training nodes
($N_{n}$) used to achieve the corresponding optimal accuracy, plotted against
increasing training set size (red). The peak accuracy, 95.7$\%$, occurs when
$N_{n}=18$ and $N_{t}=91$ (circled).

##  3 Conclusion

This work demonstrates the high classification performance of a novel laser-
based photonic spiking neural network in tackling a highly complex,
multivariate, non-linear classification task (MADELON) with 500 datapoint
features. Importantly, this work also introduces a novel ’significance’
training approach that makes use of binary weights (0 or 1), and that
leverages on the advantages of the discrete optical spiking signals found in
the photonic SNN. The experimental approach to an SNN combines the spiking
dynamics of a VCSEL with a novel network architecture inspired by the
reservoir computing paradigm to process data entirely optically at very high
speed (GHz rates). The SNN uses all-optical neuron-like spikes to create a
time multiplexed feed-forward spiking neural network, in which the values of
each time-multiplexed (virtual) node are linked through the VCSEL’s non-linear
temporal dynamics. The computational power of the photonic SNN is demonstrated
first using an OLS method of weight training. We show that by training the
output layer weights with this OLS approach we can achieve a high accuracies
up to 91$\%$ and 94.4$\%$ for the MADELON classification task using SNN
architectures of 2048 and 4096 nodes.

Additionally, in this work, we introduce a new ’significance’ training
approach, which assigns binary weights to (optical spiking) nodes according to
their overall usefulness/significance score. In this approach, only high
significance scoring nodes, nodes that spike frequently for one class but not
others, are considered and used for the network training and performance
evaluation. We show that only a very small fraction ($\leq$ 1$\%$ in the
presented case) of the total number of nodes (in the output layer) are
required to successfully classify data. We show that classification accuracies
of 94.4$\%$ and 95.7$\%$ can be achieved by this new training method. The
accuracies provided by the significance training approach show an improvement
on those achieved by the OLS method, while also significantly reducing the
number of training nodes. Moreover, the photonic SNN demonstrated
classification performance that improved upon the benchmark accuracy
(93.78$\%$) achieved by software-implemented NNs reported by the dataset
authors [47]. Additionally, we demonstrated that the photonic SNN, trained
with the new significance method, can also realise high level performance with
small training set sizes ($<$10 datapoints), further reducing the overall
resources necessary for training the optical system.

Finally, the presented photonic SNN also offers several inherent physical and
computational benefits over traditional digital semiconductor processing
systems, notably ultrafast performance (250 ps/node), low-power usage ($\sim$
150 $\mu$W average optical powers, and 3.5 mA of applied bias current)
importantly, a hardware friendly implementation (using just one VCSEL to
process all virtual nodes). Furthermore, the VCSEL-based photonic SNN can
adjust performance and processing rate by changing the number of virtual nodes
used in the system, which can be done by arbitrarily and on the fly during
pre-processing. In conclusion we believe these results open possibilities for
further photonics-based processing systems that run and operate entirely on
optical hardware, and that are capable of solving highly complex tasks with
very high accuracy and ultrafast, energy-efficient operation.

## Acknowledgments

### Author Contributions

D. Owen-Newns performed the pre- and post-processing of all data. Both D.
Owen-Newns & J. Robertson performed the experimental runs of the photonic
system. A. Hurtado supervised all research efforts. All authors contributed
equally to the writing of the manuscript.

### Funding

The authors acknowledge this work was supported by the UKRI Turing AI
Acceleration Fellowships Programme (EP$/$V025198$/$1), by the European
Commission (Grant 828841-ChipAI-H2020-FETOPEN-2018-2020), and by the UK EPSRC
(EP$/$N509760$/$1, EP$/$P006973$/$1).

### Conflicts of Interest

The author(s) declare(s) that there is no conflict of interest regarding the
publication of this article.

### Data Availability

All data underpinning this publication are openly available from the
University of Strathclyde KnowledgeBase at https://doi.org/x.xxxxxxx. For the
purpose of Open Access, the author has applied a CC BY public copyright
licence to any Author Accepted Manuscript (AAM) version arising from this
submission.

## References

  * [1] Laith Alzubaidi et al.  “Review of deep learning: concepts, CNN architectures, challenges, applications, future directions”  In _Journal of Big Data_ 8.1, 2021, pp. 53  DOI: 10.1186/s40537-021-00444-8
  * [2] Nikolay L. Kazanskiy, Muhammad A. Butt and Svetlana N. Khonina  “Optical Computing: Status and Perspectives”  In _Nanomaterials_ 12.13, 2022, pp. 2171  DOI: 10.3390/nano12132171
  * [3] Mike Davies et al.  “Loihi: A Neuromorphic Manycore Processor with On-Chip Learning”  In _IEEE Micro_ 38.1, 2018, pp. 82–99  DOI: 10.1109/MM.2018.112130359
  * [4] Michael V. DeBole et al.  “TrueNorth: Accelerating From Zero to 64 Million Neurons in 10 Years”  In _Computer_ 52.5, 2019, pp. 20–29  DOI: 10.1109/MC.2019.2903009
  * [5] Christian Pehle et al.  “The BrainScaleS-2 accelerated neuromorphic system with hybrid plasticity”, 2022  arXiv: http://arxiv.org/abs/2201.11063
  * [6] C.A. Brackett  “Dense wavelength division multiplexing networks: principles and applications”  In _IEEE Journal on Selected Areas in Communications_ 8.6, 1990, pp. 948–964  DOI: 10.1109/49.57798
  * [7] Paul R. Prucnal, Bhavin J. Shastri and Malvin Carl Teich  “Neuromorphic Photonics”  CRC Press, 2017  DOI: 10.1201/9781315370590
  * [8] David A B Miller  “Attojoule Optoelectronics for Low-Energy Information Processing and Communications”  In _J. Lightwave Technol._ 35.3 OSA, 2017, pp. 346–396  URL: http://www.osapublishing.org/jlt/abstract.cfm?URI=jlt-35-3-346
  * [9] J. Feldmann et al.  “All-optical spiking neurosynaptic networks with self-learning capabilities”  In _Nature_ 569.7755, 2019, pp. 208–214  DOI: 10.1038/s41586-019-1157-8
  * [10] J. Feldmann et al.  “Parallel convolutional processing using an integrated photonic tensor core”  In _Nature_ 589.7840, 2021, pp. 52–58  DOI: 10.1038/s41586-020-03070-1
  * [11] Shaofu Xu et al.  “High-accuracy optical convolution unit architecture for convolutional neural networks by cascaded acousto-optical modulator arrays”  In _Optics Express_ 27.14, 2019, pp. 19778  DOI: 10.1364/OE.27.019778
  * [12] Alexander N. Tait et al.  “Neuromorphic photonic networks using silicon photonic weight banks”  In _Scientific Reports_ 7.1, 2017, pp. 7430  DOI: 10.1038/s41598-017-07754-z
  * [13] Thomas Ferreira Lima et al.  “Machine Learning With Neuromorphic Photonics”  In _Journal of Lightwave Technology_ 37.5, 2019, pp. 1515–1534  DOI: 10.1109/JLT.2019.2903474
  * [14] Armin Mehrabian, Yousra Al-Kabani, Volker J Sorger and Tarek El-Ghazawi  “PCNNA: A Photonic Convolutional Neural Network Accelerator”  In _2018 31st IEEE International System-on-Chip Conference (SOCC)_ IEEE, 2018, pp. 169–173  DOI: 10.1109/SOCC.2018.8618542
  * [15] Farshid Ashtiani, Alexander J. Geers and Firooz Aflatouni  “An on-chip photonic deep neural network for image classification”  In _Nature_ 606.7914, 2022, pp. 501–506  DOI: 10.1038/s41586-022-04714-0
  * [16] H. Zhang et al.  “An optical neural chip for implementing complex-valued neural network”  In _Nature Communications_ 12.1, 2021, pp. 457  DOI: 10.1038/s41467-020-20719-7
  * [17] Zaijun Chen et al.  “Deep Learning with Coherent VCSEL Neural Networks”, 2022  arXiv: http://arxiv.org/abs/2207.05329
  * [18] Paul R. Prucnal et al.  “Recent progress in semiconductor excitable lasers for photonic spike processing”  In _Advances in Optics and Photonics_ 8.2, 2016, pp. 228  DOI: 10.1364/AOP.8.000228
  * [19] Antonio Hurtado and Julien Javaloyes  “Controllable spiking patterns in long-wavelength vertical cavity surface emitting lasers for neuromorphic photonics systems”  In _Applied Physics Letters_ 107.24, 2015, pp. 241103  DOI: 10.1063/1.4937730
  * [20] J Robertson et al.  “Towards Neuromorphic Photonic Networks of Ultrafast Spiking Laser Neurons”  In _IEEE Journal of Selected Topics in Quantum Electronics_ , 2019  DOI: 10.1109/JSTQE.2019.2931215
  * [21] J Robertson, M Hejda, J Bueno and A Hurtado  “Ultrafast optical integration and pattern classification for neuromorphic photonics based on spiking VCSEL neurons”  In _Scientific Reports_ 10.1, 2020  DOI: 10.1038/s41598-020-62945-5
  * [22] Matěj Hejda, Joshua Robertson, Julián Bueno and Antonio Hurtado  “Spike-based information encoding in vertical cavity surface emitting lasers for neuromorphic photonic systems”  In _JPhys Photonics_ 2.4 IOP Publishing, 2020, pp. 44001  DOI: 10.1088/2515-7647/aba670
  * [23] J Robertson et al.  “Image edge detection with a photonic spiking VCSEL-neuron”  In _Optics Express_ 28.25, 2020, pp. 37526–37537  DOI: 10.1364/OE.408747
  * [24] Yahui Zhang et al.  “All-optical neuromorphic binary convolution with a spiking VCSEL neuron for image gradient magnitudes”  In _Photon. Res._ 9.5 OSA, 2021, pp. B201–B209  DOI: 10.1364/PRJ.412141
  * [25] Joshua Robertson et al.  “Image edge detection with a photonic spiking VCSEL-neuron”  In _Optics Express_ 28.25, 2020, pp. 37526–37537  DOI: 10.1364/OE.408747
  * [26] Yahui Zhang et al.  “Experimental demonstration of pyramidal neuron-like dynamics dominated by dendritic action potentials based on a VCSEL for all-optical XOR classification task”  In _Photonics Research_ 9.6, 2021, pp. 1055  DOI: 10.1364/PRJ.422628
  * [27] Shuiying Xiang et al.  “A review: Photonics devices, architectures, and algorithms for optical neural computing”  In _Journal of Semiconductors_ 42.2, 2021, pp. 023105  DOI: 10.1088/1674-4926/42/2/023105
  * [28] Joshua Robertson et al.  “Ultrafast neuromorphic photonic image processing with a VCSEL neuron”  In _Scientific Reports_ 12.1, 2022, pp. 4874  DOI: 10.1038/s41598-022-08703-1
  * [29] Herbert Jaeger  “The” echo state” approach to analysing and training recurrent neural networks-with an erratum note”’  In _Bonn, Germany: German National Research Center for Information Technology GMD Technical Report_ 148, 2001
  * [30] Wolfgang Maass, Thomas Natschläger and Henry Markram  “Real-Time Computing Without Stable States: A New Framework for Neural Computation Based on Perturbations”  In _Neural Computation_ 14.11, 2002, pp. 2531–2560  DOI: 10.1162/089976602760407955
  * [31] Tobias Hülser, Felix Köster, Lina Jaurigue and Kathy Lüdge  “Role of delay-times in delay-based photonic reservoir computing [Invited]”  In _Optical Materials Express_ 12.3, 2022, pp. 1214  DOI: 10.1364/OME.451016
  * [32] Andre Rohm, Lina Jaurigue and Kathy Ludge  “Reservoir Computing Using Laser Networks”  In _IEEE Journal of Selected Topics in Quantum Electronics_ 26.1, 2020, pp. 1–8  DOI: 10.1109/JSTQE.2019.2927578
  * [33] Kristof Vandoorne et al.  “Experimental demonstration of reservoir computing on a silicon photonics chip”  In _Nature Communications_ 5.1, 2014, pp. 3541  DOI: 10.1038/ncomms4541
  * [34] Daniel Brunner, Miguel C. Soriano, Claudio R. Mirasso and Ingo Fischer  “Parallel photonic information processing at gigabyte per second data rates using transient states”  In _Nature Communications_ 4.1, 2013, pp. 1364  DOI: 10.1038/ncomms2368
  * [35] Quentin Vinckier et al.  “High-performance photonic reservoir computer based on a coherently driven passive cavity”  In _Optica_ 2.5, 2015, pp. 438  DOI: 10.1364/OPTICA.2.000438
  * [36] Julián Bueno, Daniel Brunner, Miguel C. Soriano and Ingo Fischer  “Conditions for reservoir computing performance using semiconductor lasers with delayed optical feedback”  In _Optics Express_ 25.3, 2017, pp. 2401  DOI: 10.1364/OE.25.002401
  * [37] Apostolos Argyris, Julian Bueno and Ingo Fischer  “PAM-4 Transmission at 1550 nm Using Photonic Reservoir Computing Post-Processing”  In _IEEE Access_ 7, 2019, pp. 37017–37025  DOI: 10.1109/ACCESS.2019.2905422
  * [38] Jeremy Vatin, Damien Rontani and Marc Sciamanna  “Enhanced performance of a reservoir computer using polarization dynamics in VCSELs”  In _Optics Letters_ 43.18, 2018, pp. 4497  DOI: 10.1364/OL.43.004497
  * [39] Jeremy Vatin, Damien Rontani and Marc Sciamanna  “Experimental reservoir computing using VCSEL polarization dynamics”  In _Optics Express_ 27.13, 2019, pp. 18579  DOI: 10.1364/OE.27.018579
  * [40] Jeremy Vatin, Damien Rontani and Marc Sciamanna  “Experimental realization of dual task processing with a photonic reservoir computer”  In _APL Photonics_ 5.8, 2020, pp. 086105  DOI: 10.1063/5.0017574
  * [41] Julian Bueno, Joshua Robertson, Matej Hejda and Antonio Hurtado  “Comprehensive Performance Analysis of a VCSEL-Based Photonic Reservoir Computer”  In _IEEE Photonics Technology Letters_ 33.16, 2021, pp. 920–923  DOI: 10.1109/LPT.2021.3075095
  * [42] Xavier Porte et al.  “A complete, parallel and autonomous photonic neural network in a semiconductor multimode laser”  In _Journal of Physics: Photonics_ 3.2, 2021, pp. 024017  DOI: 10.1088/2515-7647/abf6bd
  * [43] Anas Skalli et al.  “Computational metrics and parameters of an injection-locked large area semiconductor laser for neural network computing”  In _arXiv e-prints_ , 2021, pp. arXiv:2112.08947  arXiv:2112.08947 [cs.ET]
  * [44] Anas Skalli et al.  “Photonic neuromorphic computing using vertical cavity semiconductor lasers”  In _Optical Materials Express_ 12.6, 2022, pp. 2395  DOI: 10.1364/OME.450926
  * [45] Dafydd Owen-Newns, Joshua Robertson, Matej Hejda and Antonio Hurtado  “GHz Rate Neuromorphic Photonic Spiking Neural Network With a Single Vertical-Cavity Surface-Emitting Laser (VCSEL)”  In _IEEE Journal of Selected Topics in Quantum Electronics_ 29.2: Optical Computing, 2023, pp. 1–10  DOI: 10.1109/JSTQE.2022.3205716
  * [46] R.. Fisher  “THE USE OF MULTIPLE MEASUREMENTS IN TAXONOMIC PROBLEMS”  In _Annals of Eugenics_ 7.2, 1936, pp. 179–188  DOI: https://doi.org/10.1111/j.1469-1809.1936.tb02137.x
  * [47] Isabelle Guyon et al.  “Competitive baseline methods set new standards for the NIPS 2003 feature selection benchmark”  In _Pattern Recognition Letters_ 28.12, 2007, pp. 1438–1444  DOI: https://doi.org/10.1016/j.patrec.2007.02.014
  * [48] J Bueno, J Robertson, M Hejda and A Hurtado  “Experimental Implementation of a Photonic Neural Network with a 1550nm-VCSEL subject to Optical Injection and Delayed Optical Feedback”  In _2020 IEEE Photonics Conference, IPC 2020 - Proceedings_ , 2020  DOI: 10.1109/IPC47351.2020.9252399


#  Quark cores in extensions of the MIT Bag model

Salil Joshi  School of Physics, University of Hyderabad, Gachibowli,
Hyderabad, India 500046  Sovan Sau  School of Physics, University of
Hyderabad, Gachibowli, Hyderabad, India 500046  Soma Sanyal  School of
Physics, University of Hyderabad, Gachibowli, Hyderabad, India 500046

###### Abstract

Recent observations of massive pulsars having masses of the order of two solar
mass pose a new challenge for compact objects such as hybrid stars and neutron
stars. Extensions of the bag model and the Nambu-Jona-Lasino model have been
used to model these stars to get higher mass stars. Quark matter has been
predicted in the cores of these massive stars. In this work we show that an
extension of the bag model, with a chemical potential dependent bag parameter
can lead to an isentropic phase transition in the core of the neutron star.
Our model shows that an EoS having all three quarks $u$, $d$ and $s$ would
lead to massive stars with stable quark matter. We find that the mass of the
stars not only depends on the bag constant but also on the mass of the strange
quark. The mass - radius ratio which determines the redshift values on the
surface, indicates that it is possible to obtain stable self bound strange
quark matter stars with reasonable values of the bag pressure which correspond
to the recent observations of large mass stars.

strange matter, stability, Equation of state.

##  I Introduction

Among the various hypothetical compact objects discussed in the literature,
the quark star and the hybrid star hold a very important place, as they are
objects capable of generating gravitational waves. The detection of
gravitational waves by recent experiments raise the possibility of
experimental determination of these special stars. Stable quark stars had been
first predicted in ref.witten ; this was followed by various other studies on
quark stars quarkstars . The possibility of quark matter inside the cores of
neutron stars has also been discussed Collins ; Baym ; alvarez . Recent
observations of relatively massive pulsars like PSR J1614-2230 and
PSRJ0348+0432, have further contributed to the challenges presented by these
stars.

The Equation Of State (EoS) of these stars present a challenge in the current
scientific scenario. The most commonly used EoS is the bag model bagmodel ;
bagmodel1 ; bagmodel2 ; bagmodel3 . Apart from the standard bag model there
have been discussions in the literature of several extensions of the bag model
extensions . Recent studies of hybrid stars have shown that it is possible to
have stars with masses of the order of two solar masses and higher. These
higher mass limits are usually set by observations of pulsar glitches. The
masses of these stars (PSR J1614-2230 and PSRJ0348+0432) star1 ; star2  can be
successfully modelled in extensions of the bag modelritammalik . In these
studies, the author had used a density dependent bag constant. Such density
dependent bag parameter has also been used by other groups to study strange
stars nprasad . Various modified bag constants have been used to understand
the quark hadron phase transition in baryon dense environments. It has been
established that bulk strange matter can be stable farhi . They can be formed
by the trapping of strange matter due to phase transitions in the early
universe or by the slow burning of the core of the neutron star which
gradually changes it to a quark star harko . Thus strange bulk matter can
exist as cores in neutron stars, as quark stars or as hybrid stars. Apart from
these there have been discussions of a color flavor locked phase in the core
of the neutron stars and strange neutron stars cfl . The Nambu Jona Lasino
model has also been used to show that the core of massive neutron stars may
contain a region of quark hybrid mass njl . All these models have been
motivated by various observational data on the mass and radius of these stars.

There is also the question of the phase boundary between the confined and the
deconfined phases. The nature of the quark hadron phase transition usually
gives the kind of boundary between the two phases. In recent times, there have
been quite a few models where a fixed entropy per baryon has been suggested at
the phase boundary mariani . The different stages of cooling of proto-neutron
stars are modelled as isentropic phase transitions. This has been done using
the Maxwell construction. Other than this, the isentropic phase transition is
also considered for neutron stars at birth using detailed simulations masuda .
The isentropic quark-hadron phase transition had been considered for the
extended bag model previously in the case of relativistic heavy ion collision
as well as the early universe leonidov ; patra . It has been conjectured that
an extended bag model with a temperature and baryon chemical potential
dependent bag pressure can give an isentropic phase boundary between the two
phases. Such modified bag models have been discussed in some detail in the
literature . In this work we use the method introduced in ref. leonidov  by
Leonidov et. al. to see whether such extension of the bag models can lead to
stable strange bulk matter with mass and radius consistent with experimental
observations. Our aim is to demonstrate that an extension of the bag model can
explain the current observables related to large mass stars with a deconfined
quark gluon plasma at its core.

Here we have characterized the properties of these massive objects using the
bag constant and the mass of the strange quark. Though we have studied other
parameters too, these two properties are the ones that we focus on. Since the
strange quark mass is already constrained from other experiments, this helps
in obtaining very specific and limited ranges for our compact objects. The
quark matter that we get in this model is energetically favourable over the 56
Fe crystal. The other requirement of any such model is the gravitational
stability. The standard way of checking the stability is through the Tolman -
Oppenheimer - Volkoff (TOV) equations. The partition function usually gives
the energy density, pressure and number density of the quark gluon plasma.
These are then substituted in the TOV equations. The TOV equations are
routinely used to study the stability of strong gravitational objects like the
neutron stars, the hybrid stars and the quark stars. The TOV equations
corresponding to the EoS of the modified bag model are then solved numerically
to obtain the mass radius ratio corresponding to various different chemical
potentials at different temperatures. We find that the EoS leads to structures
that are quite stable for a large range of chemical potentials at low
temperatures.

We also study an experimental observable that is directly related to the
compactness of the star. The mass and radius of the star allow us to calculate
the surface red shift of the star. This can be measured through direct
observation and is well documented for most stars. There are several massive
stars which are considered as candidates for strange stars. We calculate the
red shift from our model and find that it is within the range of several
candidates for strange stars. We will discuss these candidates in detail in
section IV of the paper.

In section II we present the EoS of the extended bag model used for modelling
the quark cores of the strange matter in more detail. We have divided this
section into three subsections. In the first subsection, we discuss the bag
model EoS with two massless flavors of quarks, the $u$ and the $d$. In the
second subsections, we discuss the EoS with three massless flavors of quarks,
the $u$, $d$ and the $s$. In the third subsection we discuss the EoS where the
strange quark is taken to be massive compared to the other two quarks. In
section III, we discuss the TOV equations and their numerical solutions, In
section IV, we present the results and discuss the possible candidates which
fit our model. Finally we present our conclusions in section V.

##  II The extended bag model

We first briefly discuss the extended bag model that we use in our
calculations. The original MIT bag model was the first phenomenological model
which successfully modelled the phase transition of the quark gluon plasma
(QGP) state to the hadronic phase. The phase transition was a first order
phase transition which led to inhomogeneous baryon densities in the plasma
witten . The equation of state consisted of a bag constant $B$ which was
actually the pressure difference between the two phases. Leonidov et. al
leonidov  modified the bag model by changing the bag constant to a variable
dependent on chemical potential $\mu$ and temperature $T$. Their model was
consistent with the Gibbs equilibrium criteria for a phase transition. The
model also conserved baryon number and entropy at the phase boundary.
Following Leonidov’s approach others also modified the bag constant by making
it dependent on the chemical potential and temperature patra . The discovery
of the massive neutron stars indicates that the EoS leading to the quark cores
in the neutron stars has to be reasonably stiff. We find that the isentropic
phase transition in the core with a chemical potential and temperature
dependent bag constant does provide this stiffness to the EoS. The entropy per
baryon number is also maintained at the boundary. We find that it is possible
to have self bound strange quark matter in the core within reasonable values
of bag pressure. Though an isentropic phase transition is not a strong
requirement for phase transitions in these conditions, the grand canonical
partition function used by Leonidov leonidov  and Patra patra  do generate the
stiffness that is required to have stars with a mass greater than $1.5M_{0}$
and radii of the order of $12$ kms.

We use their grand canonical partition function and calculate the number
density and the pressure of these stars. For massless quarks, both two flavors
($u$ and $d$) as well as three flavors ($u$, $d$ and $s$), the energy density
and pressure is given by,

$\begin{split}\epsilon=\frac{N_{c}N_{f}}{2}(\frac{7}{30}\pi^{2}T^{4}+\mu^{2}T^{2}+\frac{\mu^{4}}{2\pi^{2}})+\frac{\pi^{2}}{15}N_{g}T^{4}+\\\
B(\mu,T)-\frac{\partial B(\mu,T)}{\partial\mu}T-\frac{\partial
B(\mu,T)}{\partial T}\mu\end{split}$ (1)

Here $N_{c}$ and $N_{f}$ are the number of color and flavor degrees of
freedom.

$\begin{split}P=\frac{N_{c}N_{f}}{6}(\frac{7}{30}\pi^{2}T^{4}+\mu^{2}T^{2}+\frac{1}{2\pi}\mu^{4})+\\\
\frac{\pi^{2}}{45}N_{g}VT^{4}-B(\mu,T)\end{split}$ (2)

In our model, we consider the EoS in the limit of high baryon density. As was
shown by Leonidov leonidov , in this limit, the thermodynamical quantities can
be written as the sum of two parts. One part gives the zero temperature
contribution and the second part gives the finite temperature contribution.
When we consider the $s$ quark to be massive, these contributions change
significantly. Hence in the next subsections, we will first consider the
quarks to be massless and obtain the mass and radius relations for the three
flavor case. Then we will consider the strange quark to be massive and obtain
the mass and radius relations for that case separately. For the two flavor
case, the plasma is found to be meta stable.

###  II.1 Case 1: Two massless flavors ($u$ and $d$)

We first look at the two flavor $u$ and $d$ case. Here we do not consider the
leptons but maintain charge neutrality by assuming that the number of $d$
quarks is twice that of the $u$ quarks. As has been shown in ref.madsen , the
plasma stability parameters are only weakly dependent on lepton contributions
in the case of the two flavor plasma hence the lepton contribution to the
charge neutrality can be ignored. The general way to construct the EoS is to
describe the EoS in each of the two phases and then perform the Maxwell
construction to join the two phases along their common boundary. We follow
Leonidov’s approach in modifying the bag constant. This means that at the
phase boundary the following equation must hold,

$\frac{S_{q}-\frac{\partial B(\mu,T)}{\partial T}}{n_{q}-\frac{\partial
B(\mu,T)}{\partial\mu}}=\frac{S_{H}}{n_{H}}$ (3)

Here the suffix ”q” denotes the thermodynamical variables in the quark phase,
while the suffix ”H” denotes the same variables in the hadronic phase. The
modified bag’s constant in our case is the following,

$B(\mu,T)\simeq B_{0}+\mu^{2}T^{2}-\frac{\mu^{4}T^{2}}{\theta_{H}^{2}}$ (4)

where $\theta_{H}=(\mu^{2}-m_{H}^{2})^{1/2}$ and $m_{H}$ is hadron mass.

While the core consists of the quarks, the other side of the boundary is the
hadronic phase. In this case, the hadronic phase consists of a non-interacting
neutron - pion gas. The hadronic partition function, is based on a hard core
neutron neutron repulsion. Due to this, it is necessary to divide any
thermodynamical quantity for point like particles by a volume factor leonidov
. The volume factor is given by, $(1+\frac{4}{3}\pi r_{n}^{3})$, where $r_{n}$
is the radius of the neutrons. For the case of the two flavor massless quarks,
a change in the volume factor does not give any significant change in the
results. However,in the case of three flavor massless quarks, a change in the
volume factor affects the values of the chemical potential of the quarks and
the value of the corresponding bag constant for which the star is stable. The
important quantity that can bring in a change in the volume factor is the
radius of the neutron. We will discuss this later when we discuss the three
flavor massless case. As mentioned in the introduction, we need to check the
gravitational stability of the massive objects. For this we need to solve the
TOV equations using the pressure obtained from our EoS. So the pressure, in
our case turns out to be,

$P=\frac{1}{3}(\epsilon-4B_{0})-\frac{1}{3}\left(\frac{2\mu^{4}T^{2}}{\theta_{H}^{2}}-\frac{18\mu^{6}T^{2}}{\theta_{H}^{4}}\right).$
(5)

Apart from the gravitational stability, to form a stable plasma state the
energy per baryon number of the quark gluon plasma has to be calculated. For
the two flavor quark matter at low pressures, the energy per baryon number
should be larger than that of the nuclear plasma (which is around $940$ MeV)
and also a ${}^{56}Fe$ crystal which has the minimum binding energy per
nucleon madsen ; olinto . Substituting all these constraints, we find the
chemical potential is constrained between $313.5$ MeV to $331$ MeV. The bag
constant is also constrained within a narrow range between $149$ MeV to $154$
MeV. The temperature does not have much effect on these constraints as long as
it is below $1$ MeV. The maximum size, can be obtained for the minimum bag
constant and highest value of chemical potential. The baryon density turns out
to be $0.32/fm^{3}$. However, the plasma in this case is meta stable and one
cannot give a definite mass and radius to this plasma, hence we introduce the
strange quark and look at the plasma composed of three massless quark flavors.

###  II.2 Case 2: Three massless flavors ($u$, $d$ and $s$)

We now consider the plasma to have three massless flavors of quarks. We treat
the plasma as a degenerate Fermi gas with equal numbers of quarks of different
flavors. This and the fact that the plasma is beta equilibrated gives the
necessary charge neutrality to the plasma. The only change in the expressions
for the bag constant and the pressure come from the fact that the number of
flavors have been increased to three. As of now the assumption of
$\mu_{u}=\mu_{d}=\mu_{s}$ holds as all the flavors are considered to be
massless. The bag parameter becomes,

$B(\mu,T)\simeq
B_{0}+\frac{3}{2}\left(\frac{1}{9}\mu^{2}T^{2}-\frac{1}{81}\frac{\mu^{4}T^{2}}{\theta_{H}^{2}}\right)$
(6)

and the corresponding pressure equation is given by,

$P=\frac{1}{3}(\epsilon-4B_{0})-\frac{1}{2}\left(\frac{2\mu^{4}T^{2}}{\theta_{H}^{2}}+\frac{18\mu^{6}T^{2}}{\theta_{H}^{4}}\right)$
(7)

Here $B_{0}$ is the bag parameter at $T=0$. For this case too, we have to
check the stability of the plasma. Since the number of flavors has increased,
the energy per baryon will also change accordingly. The preferred state of a
stable plasma comprising of all the quark flavors should have energy per
baryon less than that the nuclear energy of $940$ MeV so that the plasma does
not hadronize. Similarly, the energy per baryon should be less than the
binding energy of the ${}^{56}Fe$ crystal. The lower energy per baryon will
make the quark gluon plasma the preferred state for the bulk matter madsen ;
olinto . As mentioned before, the stability of the plasma puts a constraint on
the values of the chemical potential as well as the bag constant. We see that
the stability of the plasma increases with smaller bag constant values.

Again similar to the previous case, since the hadronic phase is also present,
there are some parameters in the hadronic phase which also contribute to the
stability of the plasma. We find that in the three flavor case, the radius of
the neutron which contributes to the volume correction in the neutron - pion
gas does change the stability values of the core pressure. This is because we
have equated the pressure of the hadronic phase and the quark phase to obtain
the bag constant. We have plotted the pressure in the core for different
values of bag constant for different values of the neutron radius for which
the quark plasma is stable in fig 1. As can be seen from the figure, the
increasing bag constant decreases the pressure in the core.

Figure 1: Core pressure for stable masses at different bag constants for
different values of the neutron radius.

We thus find that the change in the volume corrections to the non-interacting
neutron-pion gas affect the limiting values of the pressure and the density at
the core. In all the cases, we have considered we have kept the hadronic phase
to be the same as we are interested in the quark phase. We find that as long
as the hadronic part consists of a neutron - pion gas, the only parameter that
affects the stability of the plasma is the radius of the neutron through the
volume corrections. Though this effect also occurs when the strange quark mass
is taken into account, the dominance of the strange quark mass leads to very
small variations in the core pressure due to the change in the volume
corrections.

###  II.3 Case 3: Massless $u$,$d$ and massive $s$ quarks

If we consider the $s$ quark to be massive then we cannot have
$\mu_{u}=\mu_{d}=\mu_{s}$. So it is important to consider the first order
correction terms to the temperature and the entropy. As in the previous case,
strange matter is modelled as a Fermi gas of up, down, and strange quarks.
This time the charge neutrality of the system is maintained by the electrons.
Weak interactions ($\beta$ equilibrium) will then lead to
$\mu_{s}=\mu_{d}=\mu$ and $\mu_{u}+\mu_{e}=\mu$ alcock . The number densities
are related by, $\frac{2}{3}n_{u}-\frac{1}{3}n_{u}-\frac{1}{3}n_{s}-n_{e}=0$.
This means that effectively we are dealing with only two chemical potentials
$\mu_{u}$ and $\mu_{s}$.

To deal with high baryon densities and low temperatures, finite temperature
corrections are added to the zero temperature terms. This means that,

$P_{s}\simeq P_{s}^{0}+P_{s}^{1}T^{2}+P_{s}^{2}T^{4}$ (8) $n_{s}\simeq
n_{s}^{0}+n_{s}^{1}T^{2}+n_{s}^{2}T^{4}$ (9) $S_{s}\simeq
S_{s}^{1}T+S_{s}^{2}T^{3}$ (10)

The zero temperature terms are given by,

$P_{s}^{0}=\frac{1}{6\pi^{2}}(\mu_{s}\theta_{s}(\theta_{s}^{2}-\frac{3m_{s}^{2}}{2})+\frac{3m_{s}^{4}}{2}ln(\frac{\mu_{s}+\theta_{s}}{m_{s}}))$
(11) $n_{s}^{0}=\frac{2\theta^{2}}{9\pi^{2}}$ (12)

Here $m_{s}$ is the mass of the $s$ quark and
$\theta_{s}=\sqrt{(\mu_{s}^{2}-m_{s}^{2})}$. At high baryon densities and low
temperatures $\frac{\partial B}{\partial T}$ dominates over $\frac{\partial
B}{\partial\mu}$, based on that at the phase boundary one can obtain the bag
constant as,

$B(\mu,T)\simeq
B_{0}+(\mu_{u}^{2}+\mu_{d}^{2}-\frac{2}{3}\mu_{s}\theta_{s})\frac{T^{2}}{2}-\frac{\mu_{s}T^{2}}{2\theta_{H}^{2}}(\mu_{u}^{3}+\mu_{d}^{3}+\frac{2}{3}\theta_{s}^{3})$
(13)

Once the bag constant is obtained the pressure can be expressed as,

$\begin{split}P=\frac{1}{3}(\epsilon-4B_{0}\frac{4\mu_{e}^{4}}{3\pi^{2}}+\frac{m_{s}^{4}}{\pi^{2}}ln(\frac{\mu_{s}+\theta_{s}}{m_{s}})-(\mu_{u}^{3}+\mu_{d}^{3}-\frac{1}{3}\theta_{s}^{3})(\frac{\mu_{s}T^{2}}{\theta_{H}^{2}})+\frac{\mu_{s}\theta_{s}}{2\pi^{2}}(\theta_{s}^{2}-\frac{3m_{s}^{2}}{2})\\\
-\frac{\mu_{s}^{3}\theta_{s}T^{2}}{\theta_{H}^{2}}-\frac{\mu_{s}}{2\pi^{2}}(\mu_{s}^{2}\theta_{s}+\frac{m_{s}^{2}}{2\theta_{s}(m_{s}^{2}-\mu_{s}^{2})})+\frac{\mu_{s}T^{2}}{2}(\frac{18\mu_{s}^{2}}{\theta_{H}^{4}}(\mu_{u}^{3}+\mu_{s}^{3}+\frac{2\theta_{s}^{2}}{\theta_{H}^{4}})))\end{split}$
(14)

As we can see apart from the chemical potentials of the $u$ and $s$ quark, the
main parameters are the $s$ quark mass, the temperature and the bag constant.
The $s$ quark mass is already constrained by other experiments. So though we
do get large stable mass stars with a smaller value of the $s$ quark mass, we
do not vary the $s$ quark mass below $90$ MeV . Generally this is the lower
limit of the $s$ quark mass from various other sources olive . The temperature
is also below $1$ MeV. The constraint on the temperature comes as the entropy
per baryon number has to be continuous even in the bulk. Increasing the
temperature violates this continuity and hence we are constrained to
temperatures below $1$ MeV. The electron fraction in the core increases in
this temperature range reaching a maximum at $1$ MeV, consequently the
strangeness fraction goes down. So we find that lower the temperature, higher
is the strangeness fraction in the core. However, the strangeness fraction
saturates to a value of $0.33$ and becomes independent of temperature at lower
values.

In the previous section, we had seen that some parameters in the hadronic
phase also affect the stability of the quark cores. We had found that changes
in the volume corrections to the neutron - pion gas cause variation to the
core pressure. In this case, the chemical potential of the quarks $u$ and $s$
along with the mass of the strange quark dominate the core pressure. In fact
as the chemical potential of the $u$ quark is increased, the electron fraction
and the strangeness fraction required for a stable core goes down while the
baryon number density increases to maintain the stability. The rate of
increase(decrease) of the electron fraction, the strangeness fraction and the
baryon number density has been plotted against the chemical potential of the
$u$ quark in fig 2. Generally a larger value of the chemical potential of the
$u$ quark ($\mu_{u}$) results in more stable quark cores provided the strange
quark mass is on the lower side. This is seen by calculating the energy
density of the quark core, which has to be more stable than the ${}^{56}Fe$
nucleus. For all the values in fig 2, we get stable quark matter. The
stability is checked by calculating the energy per baryon in the plasma.

Figure 2: The rate of change of the electron fraction, the strangeness
fraction and the baryon number density with respect to the chemical potential
of the $u$ quark.

Thus we have found that it is possible to have energetically stable quark
matter using three flavors of the quark gluon plasma. As these are massive
objects, we know that gravitation plays an important role in establishing the
stability of these stars. The macroscopic quantities of mass and radius are
thus calculated using the Tolman - Oppenheimer - Volkoff equations TOV . In
the next sections, we briefly describe these equations and their solutions.

##  III The Tolman - Oppenheimer - Volkoff equations

We consider the TOV equations for a spherical and isotropic metric. We assume
that the core of the quark star has a uniform density in all the directions.
The TOV equations are then given by,

$\frac{dP}{dr}=\frac{G}{r}\frac{(\epsilon(r)+P(r))(M(r)c^{2}+4\pi
r^{3}P(r)}{(rc^{2}-2GM(r))}$ (15)

and

$\frac{dM}{dr}=\frac{4\pi r^{2}\epsilon(r)}{c^{2}}$ (16)

where the mass density $\rho(r)=\frac{\epsilon(r)}{c^{2}}$ and $c$ (cm/sec) is
the speed of light and $\epsilon(r)$ is the energy density of the plasma in
MeV. The TOV equations are solved numerically with appropriate boundary
conditions to obtain the mass distribution of the stars. The solution also
gives us the gravitationally allowed values of the mass and radius of the
stars. We solve the TOV equations numerically for the high baryon density
regime. We find that stable stars are possible for various values of the bag
constant. Detailed results and graphs are shown in the next section.

Apart from the mass radius ratio, we would also like to look at a direct
observable that is affected by the compactness of the star. One such
observable is the surface redshift values of the stars. It is a parameter that
corresponds to the redshift experienced by a radially propagating photon
travelling from the star’s surface to infinity. It can be directly calculated
from the mass - radius ratio of a star.

$z_{s}=\left(1-\frac{2GM}{Rc^{2}}\right)^{-1/2}-1$ (17)

For all the cases, where we have a stable configuration from plasma stability
as well as gravitational stability, we have calculated the red shift values of
the stars. In the next section, we will show that the surface redshift values
obtained for the strange stars in our case are consistent with experimentally
observed red shift values for different strange star candidates.

##  IV Results

###  IV.1 Mass-radius ratios with three massless flavors ($u$,$d$,$s$)

Since the two flavor plasma only gives us a metastable plasma we do not solve
the TOV equations for the two flavor plasma. In the case of the three massless
flavors of quarks we vary the various parameters and do a systematic search to
put constraints on the various parameters of the model.

Figure 3: $M$ (in solar masses) Vs $R$ at different bag constants for massless
quarks.

For massless quarks, we do get large mass values close to $2.2M_{0}$ with a
radius of $14$ km. This occurs for lower bag constants. The strange quark is
considered massless here and therefore we get the mass values greater than two
solar masses. We then obtain the surface red shift values for these large mass
stars and the resultant plot is shown in fig 4.

Figure 4: Red shift values for cores with massless quarks.

As seen from fig 4., we obtain redshift values of $0.4$ for massive stars of
$2M_{0}$ or higher. It has been predicted that PSR J0348+0432 has a radius of
about $12$ kms and mass of $2.1M_{0}$ star1  with a red shift value of $0.4$.
So our model with massless quarks is able to reproduce the mass - ratio
required to have stable quark stars as seen from the surface red shift data.
To obtain the large mass stars, our bag constant should be around $130-140$
MeV.

###  IV.2 Mass-radius ratios for the massive strange quark

We now present the results of the TOV equations with the pressure obtained
from the EoS corresponding to the massive strange quarks. The mass of the
quark puts a constraint on the quark chemical potentials. We plot the mass
radius ratio for stable stars for two different masses of the strange quark.
Since we know that the larger mass stars would occur in the range of bag
constant values between $130-140$ MeV, we give the results for two values
within this range.

Figure 5:  $M$ (in solar masses) vs $R$ at different bag constants for
$m_{s}=100$ MeV, $\mu_{u}=100$ MeV for a u-d-s massive plasma  Figure 6: $M$
(in solar masses) vs $R$ at different bag constants for $m_{s}=150$ MeV,
$\mu_{u}=100$ MeV for a u-d-s massive plasma

As we see from fig 5 and fig 6, we get larger mass with smaller values of
strange quark mass. So if the star has massive strange quarks, it’s size will
be constrained by the strange quark mass. We have plotted the redshift values
for different masses of the stars.

Figure 7: Redshift values with the plasma having a massive strange quark of
mass $m_{s}=100$ MeV  Figure 8:  Redshift values with the plasma having a
massive strange quark of mass $m_{s}=150$ MeV

Here since we have constrained the strange quark value to a lower limit of
$90$ MeV, we do not get massive stars with masses greater than $2M_{0}$. We
can get such values if we lower the strange quark values to around $75$ MeV.
We find that stars with masses of the order of $1.5M_{0}$ would typically have
surface red shift values around 0.22. This has also been seen in other models
star1 ; star2 . Thus the observational constraints given by red - shift values
are all satisfied with different values of strange quark mass in our model. We
now discuss some observational data in support of our model.

There are many observational instances of stars with these kind of parameters.
The XTE J1739-285 with a mass of $1.51$ solar mass and 10.9 km radius star3 ,
is often referred to as a quark star. It is rapidly spinning and has a red
shift ranging from $1.8$ to values of $2.4$ depending on its radius. Our model
for the quark star with a massive strange quark of $100$ MeV and bag constants
between $130-135$ MeV gives red shifts in this range. The EXO 0748-676 with a
mass of $2.1$ solar masses and $13.8$ km in radius star4  has a red shift
value of $0.35$. The EXO 0748-676 can be a quark star with a bag constant of
$135$ MeV as per this model. For stars like PSR J1614-2230, the gravitational
redshift values are in the range $0.41-0.5$, such high redshift values occur
for the quark plasma with massless quarks. As we have mentioned before the
mass - radius ratio and the red shift values depend crucially on the masses of
the strange quark.

##  V Conclusions

In conclusion, we have looked at bulk strange matter using a bag constant
which is dependent on chemical potential and temperature. The bag constant in
this model had originally been derived by Leonidov et. al. for two massless
flavors of quarks. In this work we have extended their model to include the
strange quark. Though an isentropic phase transition is not essential in the
core of the star, we show that such a phase transition can lead to stable
quark matter inside the stars. We have done a systematic study of the
parameters of the model. We have shown that the presence of the massless
strange quark increases the stability of the plasma and the mass limits of
these stars. The stability of the plasma limits our bag constant to the lowest
value of $130$ MeV. As is already known, the stability of the quark matter
increases with decrease of the bag constant. However, in this work, we have
also shown that apart from the bag constant, it is the mass of the strange
quark that determines the radius and mass of the massive object. Since the
mass-radius ratio determines the surface redshift of the star, we have also
obtained the surface red shift value for both massless and massive strange
quark stars. We have shown that our model can reproduce the mass- radius ratio
and the surface red shift of stars like PSR J1614-2230 and PSR J0348-0432 for
bag constant values in the range $130-140$ MeV. The mass of the strange quark
is important and we get larger mass stars for lower values of strange quark
mass. Though we have maintained the lower value of strange quark mass at $90$
MeV due to constraints from other sources, there are previous studies where
the strange quark mass is considered from $80$ MeV nanapan  onwards. Lowering
the strange quark mass will result in higher mass stars. Thus we have
established that it is quite possible to have a stable quark gluon plasma bulk
phase in the core of a massive star with no other exotic phases.

Apart from PSR J1614-2230 and PSR J0348-0432 already mentioned before, both
mass-radius ratios and the red shifts of XTE J1739-285 and EXO 0748-676 can be
explained by our model. Thus we do find a large number of candidates with
larger masses and surface red shifts which fit our model. We have not
considered the rotational effects of this model in this work, we have studied
the stability of the quark mass only in the bulk. We have established that it
is possible to model masses of the order of two solar masses using extensions
of the three flavour bag model with a $s$ quark whose mass essentially
determines the mass-radius ratio of the star. Our future plan is to study the
rotational effects of these stars.

Acknowledgments

S.J would like to acknowledge partial financial support (through long term
project scheme) from the UGC Networking Resource Center, School of Physics,
University of Hyderabad.

## References

  * (1) E. Witten, Phys. Rev. D 30, 272 (1984);
  * (2) S. Banerjee, S. K. Ghosh and S. Raha, J. Phys. G26, L1 (2000);
  * (3) J. C. Collins and M. J. Perry, Phys. Rev. Lett. 34, 1353 (1975);
  * (4) G. Baym , Physica A, 96, 131-135 (1979);
  * (5) E. Alvarez and J. M. Ibáñez, Astron. Astrophys., 98, 390-396 (1981);
  * (6) A. Rosenhauer, E. F. Staubo, L. P. Csernai, T. Øvergård and E. Østgaard, Nucl. Phys. A, 540 (1992) 630-645;
  * (7) M. Kutschera and A. Kotlorz; Astrophysical Journal, v 419, 752 (1993);
  * (8) G. H. Bordbar, M. Bigdeli and T. Yazdizadeh, IJMPA vol 21, No 28, 5991 (2006);
  * (9) Hua Li, Xin-Lian Luo, and Hong-Shi Zong, Phys. Rev. D 82, 065017 (2010).
  * (10) A. Aziz, S. Ray, F. Rahaman, M. Khlopov and B. K. Guha, Int.J.Mod.Phys. D, 28, 1941006 (2019).
  * (11) F. Özel, D. Psaltis, S. Ransom, P. Demorest, M. Alford, Astrophys. J. 724, L199 (2010); X.F. Zhao, Acta Phys. Pol. B 44, 211 (2013); X.F. Zhao, A.J. Dong, H.Y. Jia, Acta Phys. Pol. B 43, 1783 (2012)
  * (12) J. L. Huo, X. F. Zhao, Chinese Journal of Physics, 56, 292 (2018); X. F. Zhao, Int. J. Mod. Phys. D 24, 1550058 (2015).
  * (13) R. Mallick, Phys. Rev. C 87, 025804, (2013).
  * (14) N. Prasad and R. S. Bhalerao, Phys. Rev. D 69, 103001 (2004).
  * (15) E. Farhi and R. L. Jaffe, Phys. Rev. D. 30, 2379, (1984).
  * (16) T.Harko, K.S. Cheng and P.S. Tang, The Astrophysical Journal 608, No 2, 945 (2004).
  * (17) J. L. Zdunik and P. Haensel, Astronomy and Astrophysics, A61, 555 (2013); V. Soni and D. Bhattacharya, Physics Letters B 643, 158 (2006).
  * (18) M. Orsaria, H. Rodrigues, F. Weber and G. A. Contrera, Phys.Rev. D87, 023001, (2013).
  * (19) M. Mariani, M. Orsaria, and H. Vucetich, Astronomy and Astrophysics, Volume 601, A21, Pgs 12, (2017).
  * (20) K. Masuda, T. Hatsuda and T. Takatsuka, Prog. Theor. Exp. Phys. 021, (2016); K. Masuda, T. Hatsuda and T. Takatsuka, Eur. Phys. J. A 52: 65, (2016).
  * (21) A. Leonidov, K. Redlich, H. Satz, E. Suhonen, and G. Weber, Phys. Rev. D 50, 4657 (1994);
  * (22) B. K. Patra and C. P. Singh, Z. Phys. C 74, 699—703 (1997); B. K. Patra and C. P. Singh, Phys. Rev. D 54, 3551 (1996).
  * (23) J. Madsen, Lect.Notes Phys.516:162-203, (1999).
  * (24) A. V. Olinto, The Physics of Strange Matter, Proceedings of the 2nd International Workshop on Relativistic Aspects of Nuclear Physics, (1991).
  * (25) C. Alcock, E. Farhi, and A. Olinto, Astrophys. J. 310, 261 (1986).
  * (26) K. A. Olive and Particle Data Group, Chinese Physics C 38, 090001 (2014).
  * (27) R. C. Tolman, Phys. Rev. 55, 364; J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374.
  * (28) C. M. Zhang, H. X. Yin, Y. H. Zhao, Y. C. Wei, X. D. Li, Publ.Astron.Soc.Pac., vol. 119 pg. 1108 (2007).
  * (29) G. Zhang, M. Méndez, P. Jonker, B. Hiemstra, MNRAS, 414, 1077 (2011);
  * (30) Na-Na Pan and Xiao - Ping Zheng, Chinese Journal of Astronomy and Astrophysics, Vol:7 Pgs 675-684, (2007).