# Non-spectroscopic sensing enabled by electro-optical reservoir computer

Gleb Anufriev  1 David Furniss  1 Mark Farries  2 and Sendy Phang1  1George
Green Institute for Electromagnetics Research, University of Nottingham, NG7
2RD, UK
2 Gooch & Housego, Broomhill Way, Torquay, TQ2 2QL, UK
<EMAIL_ADDRESS>
###### Abstract

A fully bio-inspired apparatus to perform non-spectroscopic sensing to
discriminate, and determine the constituent concentrations of a chemical
mixture is proposed. Here, fully bio-inspired means that it is comprised of a
biomimetic sensor and a neuromorphic signal processor. The sensor is similar
to the human eye in terms of its trichromat architecture and overlapping
spectral response. The information processor is a neuromorphic system based on
an electro-optical implementation of a reservoir computer. A time-stepping
signal algorithm based on the $\mathcal{Z}$-bilinear transformation has been
developed to realistically simulate the electro-optical reservoir computer;
this demonstrates the discrimination and chemical concentration determination
tasks. We believe such an apparatus offers potential benefits in areas in
which chemical composition needs to be monitored in real time, for example in
chemical processing, and food-beverage (fermentation/brewing) and
environmental monitoring.

††journal: ome††articletype: Research Article

##  1 Introduction

The human nervous system is an efficient and powerful ‘biological computer’.
It consumes 20 Watts of power, $<1\%$ of the power consumption of the average
kettle and can process an enormous stream of information in real-time [1, 2,
3]. Key to the vast capability of the human nervous system lies is how and
where such information is processed. An electronic computer processes
information at a centralised processing unit (von Neumann architecture) with a
separate memory device. The human nervous system is, however, based on a
distributed processing architecture in which each neuron “remembers and
processes” a stream of (electro-chemical) signals as it travels to the brain
for “inferring”. Inspired by the distributed computing architecture of the
human nervous system, several neuromorphic computing architectures have been
developed[4, 5, 6, 7]. The neuromorphic computing system aims to realise the
distributed computing architecture of the human nervous system as hardware,
enabling a hardware-level integration of Artificial Intelligence.

Photonic Reservoir Computer (PhRC) is a relatively new member of the
neuromorphic computing family implemented as a physical photonic/optical
system. Implementing a reservoir computer in photonics leads to a design for
very high speed processing exploiting the availability of ultrawide bandwidth
optical communication commponents. It offers a fundamentally different
approach to other neuromorphic photonic systems, in terms of both architecture
and how the system is trained (or optimised) [4, 5, 6, 7, 8]. PhRC only
trains/optimises the Read-Out whilst retaining the characteristic of the
kernel as a random and semi-chaotic system, whereas other neuromorphic systems
require both kernel and read-out training/optimisation[4, 5, 6, 7]. Such a
training approach makes the PhRC more robust and less sensitive to
perturbation. Consequently, the PhRC offers a more straightforward
implementation and scalability than other neuromorphic photonic
architectures[4, 5, 6, 7, 8]. Several implementations of PhRC have been
reported, for example as a swirl configuration of semiconductor amplifiers [9,
10], or as a photonic crystal cavity [11, 12, 13]. In this paper, we consider
an electro-optical reservoir computer (EORC) implementation with a single
delay line and a Mach-Zehnder modulator (MZM) as the non-linear neuron
activation function [14, 15].

Humans’ powerful computing capacity has mainly been used to interpret sensory
information, $~{}99.99\%$ of the information processing capacity is for
sensory processing while $<0.01\%$ is for attention, decision-making, motion
and language [16]. Motivated by such a fact, in our previous work[12, 13, 17],
we have numerically demonstrated a new methodology to process ultrafast high-
throughput molecular vibrational signals (spectroscopic) using a chaotic
photonic crystal billiard system as its kernel. The present paper reports an
alternative method which is fully-bio-inspired, meaning that the sensor system
is biomimetic of the human eyes and its output will be processed through a
neuromorphic EORC system. Numerical simulation successfully demonstrates that
both chemical discrimination and concentration determination tasks can be
achieved.

The paper is organised as follows: Section 2 describes the system as a whole
and its composition, namely the biomimetic sensor and the EORC and the
modelling approach used. Section 3 discusses the results obtained. The first
part considers the dynamics of the EORC system, which is modelled using a
realistic model considering the filtering effects observed in a practical
setup. The second part demonstrates the application of such a system to
perform the chemical discrimination and concentration determination tasks. The
impact of EORC parameters on the accuracy is investigated and discussed. Some
concluding remarks are then provided in Section 4.

##  2 Biomimetic sensing by neuromorphic signal processing

Figure 1(a) schematically shows the apparatus considered in this paper. It is
comprised of a biomimetic sensor serving as the sensing component, an electro-
optic reservoir computer (EORC) serving as signal processing component and a
read-out component serving as an inferring component of the EORC output.
Figure 1(b) schematically depicts the implementation of Fig. 1(a) using a
single delay-line configuration and a Mach-Zehnder modulator (MZM) as the
nonlinear neuron activation kernel. The biomimetic sensor supplies a three-
channel reading to the arbitrary waveform generator which performs time-
multiplexing on the sensing information for input to the EORC system.

###  2.1 Biomimetic sensor system

This sub-section describes the operating principle and the modelling
procedures of the sensing component used in this work. It is inspired by the
biomimetic design which previously has been demonstrated successfully for
chemical [18, 19, 20] and biological sample [21] detection. The biomimetic
sensor mimics the operation of the human eyes, which are capable of
discriminating over 100,000 shades of (visible) colour over an environmental
background despite being equipped with only three main types of cone cells
[22, 23]. Key to this ability is the spectral sensitivity of these cone cell
types which are sprectrally broad and overlapping [24]. Research on colour
sensitivity of mantis shrimp (Haptosquilla trispinosa), which has 12 different
types of photosensitive cell covering different light wavelengths and
polarisations, however, found that higher number of unique photoreceptors does
not mean an improvement in colour discimination [25]; Mantis’s spectral
discrimination $\Delta\lambda$ is in the region of 15 to 25 nm in comparison
to humans’ which is in the region of 1 to 8 nm [25, 26] over the visible light
region. In [18, 19, 20, 21], a trichromat, three-colour, biomimetic sensor has
been constructed in the mid-infrared range of light and a multi-variate
statistics tool based on a patented comparative-discrimination spectral
detection (CDSD) technique [20, 27] has been used to perform the vectorial
signal processing. In the present work, we have adopted a biomimetic sensor
with a similar trichromat architecture, but the EORC is used for signal
processing. The EORC will be described in sub-section 2.2.

Figure 1: Schematic illustration of the system as a whole; (b) Implementation
of (a) in a single physical node configuration with electro-optic with Mach-
Zehnder modulator (MZM) as nonlinear neuron activation kernel; The multi-
channels output of the biomimetic sensor has been representatively shown in a
Red-Green-Blue (RGB) colouring scheme. (c) Signal flow diagram translating the
flow of signal in (b); (d) Schematic of the architecture of algorithmic
reservoir computer.

As mentioned, the biomimetic sensor used in this paper comprised of three
channels whose spectral responses are different, see Fig. 2. Figure 2(a)
schematically shows the implementation of the biomimetic sensor constructed
in[18, 19, 20]. The photodetectors are equipped with a filter at different
central wavelength, in Fig. 2(a) the filters have been representatively shown
in Red-Green-Blue (RGB) layer in front of the detectors. In the present work,
we have specifically considered that each channel has a Gaussian spectral
response peaking at different wavelengths, i.e., 3.37, 3.45 and 3.53 $\mu$m
and with a half-power bandwidth (HPBW) of 60 nm. Such sensors are readily
available, for example as dielectric filtered photodetectors[28]. The specific
parameter values of the sensor have been chosen as such they are broad and
spectrally overlapped[18, 19, 20, 21], see Fig. 2(b). The choice of peak
spectral response around the 3 µm of the mid-infrared range is motivated by
the richness of carbon-hydrogen stretching fingerprints in this region[29].
Note that in previous works[18, 19, 20, 21], the sensor has also operated at
similar wavelengths of interest.

For demonstration, in this work methanol, ethanol and isopropanol and their
mixture will be considered as chemicals to be detected. The complex-valued
frequency-dependent refractive indices of the pure substances have been taken
from [30] whereas the refractive indices of mixture are calculated using the
Lorentz-Lorenz rule [31, 32, 33] by,

$\displaystyle\frac{n^{2}-1}{n^{2}+2}=\sum_{k}\varrho_{k}\frac{n_{k}^{2}-1}{n_{k}^{2}+2}$
(1)

in (1), $n$ and $n_{k}$ are the refractive indices of the mixture and the pure
components, respectively. The $\varrho_{k}$ denotes the relative volume
mixture of component $k$, with $\sum_{k}\varrho_{k}=1$. The transmitted power
ratio, $T$, upon absorption by the sample, is calculated in the absence of any
scattering and loss due to the sample holder as

$\displaystyle T(f,n)=\exp\left(-\frac{4\pi f|n_{I}(f)|\ell}{c_{0}}\right)$
(2)

In (2), $f$ and $c_{0}$ are the frequency and speed of light respectively.
Note that the refractive index $n=n_{R}+i|n_{I}|$ convention is used. The
thickness of the sample $\ell=5$ $\mu$m is fixed throughout this work. It is
important to note that the treatment used here to calculate the transmitted
light in the mixture sample is consistent with the Maxwell’s equations of
electromagnetics; for a detailed discussion readers are referred to Review
[33]. For reference, the absorption, $\alpha=1-T$, of pure methanol, ethanol
and isopropanol is overlaid over the spectral response of the biomimetic
sensors in Fig. 2(b). The power observed by the photodetector is calculated
by,

$\displaystyle
P_{m}(n)=\int_{0}^{\infty}\mathcal{S}_{\mathrm{in}}(f)\mathcal{D}_{m}(f)T(f,n)\mathrm{d}f,\qquad
m=1,2,3$ (3)

where, $\mathcal{D}_{m}(f)$ defines the spectral response of the photodetector
at channel $m$ of the sensor. For demonstration purposes, the sensor is
assumed to be equipped with a broadband light source with a top-hat spectral
power density $\mathcal{S}_{\mathrm{in}}$ profile from 3.1 to 3.7 $\mu$m in
wavelength that is normalised, i.e.,
$\int\mathcal{S}_{\mathrm{in}}(f)\mathrm{d}f=1$ W. For practical realisations,
photo-luminescent emission from rare earth doped optical fibre [34] or a
supercontinuum light-source [35, 36] can be used. Figure 8(a) shows the
pseudo-colouring of the biomimetic sensor output
$\mathbf{P}(n)=[P_{1};P_{2};P_{3}]^{\mathrm{T}}$ following the Red-Green-Blue
(RGB) colour scheme. To allow better contrast for visualisation, the sensor
output has been normalised during the pseudo-colouring, that is
$\mathbf{P}_{\mathrm{visual}}(n)=\mathbf{P}(n)/\max[\\{\mathbf{P}\\}]$, where
$\\{\mathbf{P}\\}$ denotes the whole set of $\mathbf{P}(n)$. Note that only on
Figure 8(a), normalisation has been performed. Figure 8 will be discussed in
detail in Section 3.3.

Figure 2: Schematic design of the three-channel biomimetic sensor; (b) the
absorption spectra $\alpha$ and the filter response $\mathcal{D}_{m}(f)$ of
the photodetector (PD).

###  2.2 Electro-optical reservoir computer (EORC)

In the previous sub-section, the operating principle and model description of
the biomimetic sensor was described. In this sub-Section, its physical
implementation using electro-optical (EO) modulator system and its numerical
model are described. The EORC is a physical implementation of the Reservoir
Computer (RC) computational framework which is a generalisation of a
traditional artificial neural-network (ANN) based on the recurrent neural-
network (RNN) [7, 37, 38, 39]. For detail on the framework of the RC, the
reader is referred to [39]. Figure 1(d) shows a schematic of the architecture
of an RC scheme. An RC system consists of three main parts, namely an input, a
kernel and a Read-Out layer. The origin of the name ‘Reservoir Computer’ comes
from the architecture of the RC which features a ‘kernel’ that is a reservoir
containing a large number of neurons.

Figure 1(b) schematically shows the EORC system considered which is a
combination of optical (in blue) and (electronic) radio frequency (RF) (in
black) systems. The EORC is a feedback system in which the signal is delayed
by a time $\tau$ due to the presence of the long delay line. The delayed
feedback signal is combined with the present-time input signal to modulate the
MZM completing the feedback loop. The following sub-sections describe the
implementation of each part of the RC as an EORC system.

####  2.2.1 Input node by signal time multiplexing

The role of the input node is to distribute the input signal $u(t)$ to the
reservoir kernel. In the algorithmic reservoir computer case, one introduces a
numerical input weight $\mathbf{W}^{\mathrm{in}}\in\mathbb{R}^{N_{x}\times
N_{u}}$ which is a matrix of size, $N_{x}\times N_{u}$, where $N_{x}$ and
$N_{u}$ are the number of neuron nodes and the number of input signals,
respectively. For the translation system as an EORC, that is considered here,
neuron nodes are created artificially by means of time multiplexing and
masking during the information encoding process; such artificial temporal
nodes have also been referred to as virtual neuron nodes [14]. Practically
this process done by the Arbitrary Waveform Generator (AWG), see Fig. 1(b).

Figure 3: (a) Serialised information signal $j(t)$ which is based on the
outputs of the biomimetic sensor $P_{1}$, $P_{2}$ and $P_{3}$; (b)
Illustration of the mask signal $m(t)=m(t+\tau)$; (c) Input signal generated
upon application of mask on serialised information signal; (d) Illustration of
the time-demultiplexing process to generate virtual nodes.

The signal time-multiplexing is achieved by performing a sample and hold
operation, that is the measurement from each channel is serialised for a
duration of $\tau$. For the three-channel output biomimetic sensor considered
in this paper, i.e., $\mathbf{P}(n)=[P_{1};P_{2};P_{3}]^{\mathrm{T}}$, the
serialised signal $j(t)$ is illustrated in Fig. 3(a). This time-multiplexing
process allows us to encode the sensor information as a time-series signal.

Subsequently, a periodic mask signal $m(t)=m(t+\tau)$ is generated whose
values are randomly chosen from some random generator with mean of 1. This
mask is piece-wise constant over the period of $\theta=\tau/N_{x}$. Here, the
number of mask points is equal to the number of neuron nodes $N_{x}$, because,
as mentioned, the mask serves as the serialised input weight of the
information signal with the virtual node. Figure 3(b) illustrates the mask
signal for $N_{x}=10$. The input signal is obtained by application of mask
signal to the time-multiplexed signal, i.e., $u(t)=m(t)j(t)$, see Fig. 3(c).
The subsequent process applied upon input signal $u(t)$ by the ‘neuron’ in the
reservoir is described in the following.

####  2.2.2 Kernel by a single delay-line non-linear node

The input signal $u(t)$, here simulated by the Arbitrary Waveform Generator,
will be combined with the delayed feedback signal $x(t-\tau)$ and transformed
via a non-linear modulation process occurring in the neuron kernel. Throughout
this paper, we consider a neuron kernel implementation based on the single
physical node delay line with a Mach-Zehnder modulator (MZM) as the non-linear
neuron activation function, see Fig. 1(b). Implementation of a RC based on
this system and its variant have been reported previously in the context of
all-optical computing and communication applications [14, 40]. Initially, a
similar setup to Fig. 1(a) was proposed by Ikeda [41, 42, 43] to demonstrate
the existence of multiple stability state and chaos in a non-linear optical
ring resonator system, mathematically, modelled by the Ikeda map given by,

$\displaystyle v(t)=f_{\mathrm{NL}}[v(t-\tau)]$ (4)

where, the $f_{\mathrm{NL}}[\cdot]$ denotes a non-linear function and $\tau$
denotes the time delay parameter. Considering the specific implementation
based on MZM as the non-linear function kernel, here, we have

$\displaystyle
f_{\mathrm{NL}}[v(t)]=I_{0}\cos^{2}\left(\frac{\pi}{2V_{\pi}^{\mathrm{RF}}}v(t)+\Phi_{0}\right)$
(5)

Equation (5) is the well-known non-linear transmittance of a MZM due to two-
travelling wave interference process [44]. Parameter $I_{0}$ is the intensity
of continuous-wave laser which is being modulated by the driving RF signal
$v(t)$. Parameter $V_{\pi}^{\mathrm{RF}}$ is the half-wave voltage of the MZM
i.e. the voltage which would phase shift the light by $\pi$. Parameter
$\Phi_{0}$ is the initial phase difference between the two-travelling wave in
the MZM. Commercial MZMs are usually fabricated with $\Phi_{0}\approx-\pi/4$,
and so operates near the (almost) linear region of (5); this operation is
commonly referred to as the ‘Quad’ point. For numerical simplicity, here, the
half-wave voltage is taken as $V_{\pi}^{\mathrm{RF}}=\pi/2$.

The original Ikeda map (4) has been adapted to model a realistic experimental
setup [15, 45]. Specifically, low-pass and bandpass filter effects are
included to describe the inclusion of the photodetector and RF signal
amplifier, respectively [15, 45]. This leads to the following delay
differential equations (DDEs),

$\displaystyle\delta\frac{\mathrm{d}v}{\mathrm{d}t}+v=r,$ (6) $\displaystyle
v+\frac{2m}{\omega_{0}}\frac{\mathrm{d}v}{\mathrm{d}t}+\frac{1}{\omega_{0}^{2}}\frac{\mathrm{d}^{2}v}{\mathrm{d}t^{2}}=\frac{2m}{\omega_{0}}\frac{\mathrm{d}r}{\mathrm{d}t},$
(7)

where $r(t):=x(t-\tau)+u(t)=f_{\mathrm{NL}}[\gamma v(t-\tau)]+u(t)$. The Ikeda
map with the inclusion of the low-pass filter is given by (6) and the bandpass
filter effect is described by (7). Parameter $\delta$ in (6) defines the high-
frequency cut-off of the low-pass filter $f_{H}=1/(2\pi\delta)$. In the
bandpass filter case of (7), the filter parameters
$m=\frac{1}{2}\sqrt{f_{H}/f_{L}}$ denotes the damping factor and
$\omega_{0}=2\pi\sqrt{f_{L}f_{H}}$ denotes the internal pseudo-oscillator
resonant frequency of the bandpass filter whose values depend on the low-
frequency $f_{L}$ and high-frequency $f_{H}$ cut-off of the filter, see
Appendix 5.1. We note that while current photodetectors can have a very high
cut-off frequency, e.g., up to 20 GHz[46], in practice the bandwidth is
further limited by the electronics setup, e.g., RF amplifier which is a
bandpass filter[47]. It is noted here that (7) has also been formulated as an
integro-differential delay equation (iDDE) as shown in [45].The impact of
filtering will be discussed in detail in Section 3.1.

The DDEs (6,7) have been represented as a signal flow diagram in Fig. 1(c).
Figure 1(c) translates the mathematical operation performed at each stage of
the signal’s flow and uses the notation that will be used throughout this
paper. Parameter $\gamma$ has been added to represent the effective gain of
the system. Operator $h(t)$ represents the time-domain operation of the filter
effect, implementing DDE (6) or (7). In this paper, we have implemented a
time-stepping algorithm, performing the integration of (6) or (7) up to
second-order accuracy, based on the $\mathcal{Z}$-bilinear transformation
technique. This technique has previously used in conjunction with the
transmission-line method to model non-linear photonic materials described by a
Duffing model [48, 49, 50, 51]. Specific implementation of the filters (6) and
(7) is provided in Appendix 5.1. Other methods [43, 52, 53] could also be
used.

The activation state signal $x(t)$, as the observable, is monitored and
demultiplexed for further weighting by the Read-out layer. The subsequent sub-
section describes the time-demultiplexing process, and the weight applied to
the activation state signal.

####  2.2.3 Read-Out layer

To allow weighting by the Read-Out layer, the activation state signal $x(t)$
is sampled and demultiplexed. Figure 3(d) schematically illustrates the time-
demultiplexing process. The neuron activation state signal $x(t)$ is sampled
with a sampling-time of $\theta$ and spun with a winding period of $\tau$. The
demultiplexed activation signal is partitioned as
$\mathbf{x}(\tilde{t})=[x_{1}(\tilde{t}),x_{2}(\tilde{t}),\cdots,x_{k}(\tilde{t})]^{\mathrm{T}}$
where $k=1,2,3,\cdots,N_{x}$ is the index of the virtual nodes and each
entries of matrix $\mathbf{x}$ are given by
$x_{k}(\tilde{t}\tau+k\theta)=x(t)$ with $\tilde{t}=0,1,\cdots$ being the
winding index.

In the following step, the demultiplexed signal is weighted by the read-out
layer,

$\displaystyle\mathbf{y}(\tilde{t})=\mathbf{W}^{\mathrm{out}}[1;\mathbf{x}(\tilde{t});\mathbf{y}(\tilde{t}-1)]^{\mathrm{T}},$
(8)

where $\mathbf{W}^{\mathrm{out}}\in\mathbb{R}^{N_{y}\times(1+N_{x}+N_{y})}$ is
the read-out weighting. Note, that in (8) there is constant bias of 1 which
has a similar role as in a feed-forward neural network to allow shifting of
discrimination plane [54]. The read-out weight $\mathbf{W}^{\mathrm{out}}$ is
obtained following a training session. Figure 1(c) shows the flow of signals
in the normal operation in black whilst the flow of signals in the training
operation is shown in red. During the training operation, given the pair of
activation state signal $\mathbf{x}(\tilde{t})$ and the teacher signal
$\mathbf{y}_{\mathrm{teacher}}(\tilde{t})$, the weight of the read-out layer
is optimised. For details of the optimisation algorithm used to calculate the
read-out weight via Tikhonov regularisation with cross-validation, the reader
is referred to [13].

##  3 Results and discussion

In Section 2, the sensing component of the system in the form of a biomimetic
sensor with 3 output channels, and the architecture of the EORC system and its
model which is used to process the sensor signal have been described. This
Section starts with results and discussion of the dynamics of the single delay
line feedback system (sub-sections 3.1) before describing in sub-sections 3.2
and 3.3 the application of such a system as an EORC for chemical sensing
detection.

###  3.1 Dynamics of electro-optic reservoir computer under limited bandwidth
operation

A delay feedback system, including the Ikeda map [43], is known to exhibits
chaotic behaviour. This sub-section will discuss the influence of the
filtering effects that would be present in a practical EORC system, due to the
inclusion of photodetector and RF signal amplifier. For this purpose, we
consider the case that there is no information signal to be processed, i.e.,
$u(t)=0$.

Figure 4: (a,b) Bifurcation diagram of Ikeda map (4) and (6), respectively.
(c,d,e,f,g,h,i) The temporal signal $x(t)$ at the steady-state for different
laser intensity parameters of $I_{1}=1.5$, $I_{2}=1.9$, $I_{3}=2.1$,
$I_{4}=2.14$, $I_{5}=2.3$, $I_{6}=2.38$ and $I_{7}=2.95$ for Ikeda map (4).
Subfigure (j) is an enlarged view of (n). (k,l,m) Temporal signal $x(t)$ at
the steady-state for $I_{4}$, $I_{6}$ and $I_{7}$ for low-pass feedback Ikeda
system (6). (n) Transient state of $x(t)$ for Ikeda map (4) in blue and low-
pass feedback Ikeda system (6) in red. For convenience, the intensity
parameters considered are marked in (a,b).

Figure 4(a) shows the bifurcation diagram of the original Ikeda map (4). The
bifurcation diagram shows all possible values for neuron activation signal
$x(t)$ at the steady state as a function of laser intensity parameter $I_{0}$,
also see Fig. 1(c) for notations used. Specific parameters used for Fig. 4(a)
are $\tau=1$ $\mu$s, $\gamma=1$ and $\Phi_{0}=-0.85\pi$ similar to that use in
[45, Fig. 3]. Such a delay $\tau=$ 1 $\mu$s can, in practice, be achieved by
using a $\approx$200 m spool of a standard SMF-28 fibre. Signals $x(t)$ for
intensities of $I_{1}=1.5$, $I_{2}=1.9$ and $I_{3}=2.1$ are shown in Fig.
4(c,d,e); these specific intensity points are marked in Fig. 4(a) for
reference. When operating at $I_{1}$, the neuron signal is single valued
around 1.5, while for the case of $I_{2}$ and $I_{3}$ the signal $x(t)$ are
periodic and multivalued, i.e., two and four stationary states, respectively.
Figure 4(a) further shows that while at low intensity the steady state of
neuron activation signal $x(t)$ is single valued, it splits as $I_{0}$
increases; the splitting point is called the bifurcation point. Furthermore,
Fig. 4(a) shows that by increasing the parameter $I_{0}$ the steady state
solution splits further and begins to behave chaotically, also see Fig. 4(f)
for the case of $I_{5}=2.3$ and also, noting the bistable response at
$I_{0}=1.16$. For detail on the bifurcation and bistablilty phenomena, and the
route to chaos of the Ikeda map system, see [45].

Figure 5: Bifurcation diagram of the bandpass filter feedback of eq. (7). Blue
bullets points are the ‘fixed points’ solution of (7), see Appendix 5.2. The
line $x=I_{0}$ is included for reference; (b,c,d,e,f,g) Temporal signal
$x(t)-I_{0}/2$ at different intensity parameter of $I_{1}=0.5$, $I_{2}=0.7$,
$I_{3}=0.9$, $I_{4}=1.2$, $I_{5}=2$ and $I_{6}=3$. These intensity parameters
also are marked in (a) for reference.

The bifurcation diagram for the low-pass feedback (6) case is shown in Fig.
4(b). Figure 4(b) used the same parameters as for Fig. 4(a) in addition to the
high-frequency cut-off (-3dB) parameter of $f_{H}=10$ GHz; such a high cut-off
is currently possible due to the availability of ultrafast photodetectors,
e.g., InGaAs detectors [46]. It is noted that in the absence of any filtering
effect, the original Ikeda map (4) operates at infinite bandwidth and as such
the response of the system is instantaneous. In contrast, the DDE (6), i.e.,
low-pass feedback system, operates with a limited bandwidth. Note that DDE (6)
is a causal system whose response is dependent on a past event. Visual
comparison of Fig. 4(a) and Fig. 4(b) shows that the low-pass filter effect
has small impact for low intensity $I_{0}<2.1$. The temporal signals $x(t)$
for intensity parameters of $I_{1}$, $I_{2}$ and $I_{3}$ are also overlapped
with that of the original Ikeda map in Fig. 4(c,d,e). Deviation from the
original Ikeda map in the low-pass feedback case is noted at a high-intensity
whose value depends on the cut-off frequency $f_{H}$; for $f_{H}=10$ GHz in
Fig. 4(b), this occurs at $I_{4}\approx 2.14$. Comparing Fig. 4(g) and 4(k),
there is a noticeable presence of ripples at the stationary point in the low-
pass feedback case but not in the infinite bandwidth case. This ripple
originates from the causal response of the low-pass filter which is not
instantaneous resulting a dispersion effect of the signal. This rationale is
confirmed by analysing Fig. 4(n,j). Figure 4(j) is an enlargement of Fig. 4(n)
which shows the transient of $x(t)$ when operated at $I_{4}$. It shows that,
in the low-pass feedback system, signal $x(t)$ exhibits a typical dispersive
ripple at signal’s transition whose impact accumulates as $x(t)$ evolves over
time.

The bifurcation diagram for the case of the bandpass feedback system is shown
in Fig. 5(a); the temporal signal $x(t)$ for different intensity parameters
$I_{i}$, $i=1,2,3,4,5,6$ are shown in Fig. 5(b-g). The bandpass filter
parameters are $f_{L}=10$ kHz and $f_{H}=10$ GHz, for the low and high-
frequency cut-off, see Appendix 5.1. In addition, the MZM is operated at the
Quad operational condition, i.e., $\Phi_{0}=-\pi/4$. At low intensity $I_{0}$,
signal $x(t)$ reaches a steady state value of $x(t\rightarrow\infty)=I_{0}/2$.
This is because bandpass filter removes the zero-frequency (DC) component of
signal $v(t)$, as such $v(t\rightarrow\infty)=0$ leads to
$x(t\rightarrow\infty)=I_{0}\cos^{2}(\pi/4)$. Visual inspection of Fig. 5(b –
d) shows that the time required for the system to reach steady state becomes
longer as the intensity increases. Then, at $I_{0}=1$, the neuron activation
signal begins to oscillate. Figure 5(e) which depicts the signal $x(t)$, at
$I_{0}=1.2$, oscillates at the steady state and is exhibits the skewed
waveform typically found in strong non-linear media due to the presence of
high-order harmonics [55].

Figure 6: Impact of low-frequency cut-off $f_{L}$ on the dynamic of bandpass
feedback system beyond the bifurcation point. (a) Spectra of $x(t)$ for
different cases of $f_{L}=10$, 20 and 30 kHz. The temporal evolution of $x(t)$
is shown in (b,c,d), respectively. Note that in (b,c,d), $x(t)$ is shown from
$t=0$ to steady state. In all cases the intensity parameter is $I_{0}=1.2$.

Furthermore, visual inspection of Fig. 5(e) shows that the signal $x(t)$
oscillates at longer period $\approx$ 25 $\mu$s than the round-trip period
$\tau=1$ $\mu$s. To investigate the origin of this low-frequency oscillation,
the spectrum of $x(t)$ is displayed in Fig. 6(a) for different low-frequency
cut-offs of the bandpass filter, namely $f_{L}=10$, 20 and 30 kHz. The
resulting temporal signals $x(t)$ are shown in Fig. 6(b-d), respectively.
Figure 6(a) displays different characteristic behaviour in low- and high-
frequency regions. The high-frequency region is characterised by the resonance
of the feedback system whose fundamental frequency occurs at $1/\tau=1$ MHz
and harmonics of which are produced from the nonlinear process of the MZM.
Note that these harmonics are similar regardless the different $f_{L}$
parameters.

The characteristic behaviour in the low-frequency region is, however,
significantly different for the different cases of $f_{L}$. Figure 6(a) shows
that the fundamental frequencies of the low-frequency region are increasing,
i.e., 38, 55, 67 kHz for increasing $f_{L}=10$, 20 and 30 kHz, confirming the
origin of the low-frequency oscillation as the interplay of the $f_{L}$
parameter and the nonlinear feedback effect. Analytical investigations on
similar bandpass feedback systems [56, 57] have also noted the existence of
this novel interplay effect, which is attributed to the stability of the zero
‘fixed point’. As such, the fixed points of (7) are shown in Fig. 5(a) as
bullets points; Appendix 5.2 describes the fixed points and the bifurcation
point of (7) for Quad operation. For $\gamma=1$, the bifurcation point, which
occurs at $I_{0}=1$, marks the beginning of such oscillation in $x(t)$.
Further increase of intensity leads to the generation of more harmonics which
manifest themselves as ripples, as shown in Fig. 5(f). Similar to the case
shown in Fig. 4(a and b), further increase of intensity leads to chaotic
dynamics of the system. Although not shown in Fig. 5(a), we note the existence
of a metastable regime at $I_{0}\approx 5$ after the chaotic regime.

###  3.2 Sample discrimination task

In sub-section 3.1, the impact of filter (low-pass and bandpass) on the
dynamics of the EORC was discussed. In this sub-section, the EORC is trained
and tasked to discriminate three chemicals, namely methanol, ethanol and
isopropanol.

To achieve such a task, the input signal has been generated following the
scheme described in Section 2.2. Moreover, the teacher signal
$\mathbf{y}_{\mathrm{teacher}}$ has been defined using one-hot-encoding, as

$\displaystyle\mathbf{y}_{\mathrm{teacher}}:=[y_{1};y_{2};y_{3}]^{\mathrm{T}}=\begin{cases}[1;0;0]^{\mathrm{T}},\qquad\text{for
methanol}\\\ [0;1;0]^{\mathrm{T}},\qquad\text{for ethanol}\\\
[0;0;1]^{\mathrm{T}},\qquad\text{for isopropanol}\end{cases}.$ (9)

Moreover, we only consider an EORC based on bandpass feedback system modelled
via (7); this is a realistic scenario as practical loss in the experimental
system [45, 15, 14] demands the use of an RF amplifier. The bandpass filter
parameters used throughout are $f_{L}=10$ kHz and $f_{H}=10$ GHz; such a
broadband filter is achievable in practice noting the availability of
ultrawide bandwidth photodetectors, e.g., [46] and RF amplifiers, e.g., [47].
The delay time of the EORC used is the same as in the previous Sections
$\tau=1$ $\mu$s. The impact of different numbers of virtual nodes (or number
of mask points) $N_{x}$ is now investigated.

Figure 7: Impact of EORC parameters $I_{0}$ and $\gamma$ to the accuracy of
the trichotomy task. (a) The NMSE $\varepsilon$, in logarithmic scale, is
shown as the colour as function of $I_{0}$ and $\gamma$. The Bifurcation
points are included for reference. (b) Error as function of intensity
parameter $I_{0}$ for gain parameter $\gamma=5$; (c) Impact of number of nodes
$N_{x}$ to the NMSE $\varepsilon$.

Optimum operation of the EORC, which leads to the minimum Normalised Mean
Square Error (NMSE) $\varepsilon$, occurs around the bifurcation point
$I_{0}=1/\gamma$. It is because that around the bifurcation point the EORC is
sensitive to perturbation as such is suitable for sensing/detection
application. To show the relation between the intensity and gain parameter for
the optimum operation of the EORC, Fig. 7(a) depicts the error $\varepsilon$
on a logarithmic scale as a function of intensity $I_{0}$ and effective gain
$\gamma$ parameters for $N_{x}=30$, which is calculated by

$\displaystyle\varepsilon=\frac{\langle|\mathbf{y}-\mathbf{y}_{\mathrm{teacher}}|^{2}\rangle}{\langle|\mathbf{y}_{\mathrm{teacher}}|^{2}\rangle},$
(10)

where $|\cdot|$ denotes the Euclidean norm operator and $\langle\cdot\rangle$
is the assembled averaging operation over the datapoints $\tilde{t}$ and
chemical samples. Figure 7(b) shows the error $\varepsilon$ as a function of
$I_{0}$ for $N_{x}=20$, 40 and 60; in all cases $\gamma=5$. First, it confirms
that for any $N_{x}$ the optimum operation of the system which leads to
minimum error occurs around the bifurcation point $I_{0}=1/\gamma=0.2$. To
investigate the impact of number of virtual nodes $N_{x}$, Fig. 7(c) shows
error $\varepsilon$ for $I_{0}=0.2$ and $\gamma=5$ for different $N_{x}$
parameters. It confirms that one can further improved accuracy of the system
by increasing the number of nodes $N_{x}$. Hereafter, we will use $N_{x}=50$
as this is shown to be sufficient to provide good accuracy operation.

###  3.3 Chemical concentration determination task

In sub-section 3.2, the EORC has been trained and tasked to detecting pure
methanol, ethanol and isopropanol. In this sub-section, the EORC is trained
and tasked to determine the concentration of methanol, ethanol and isopropanol
mixture. To achieve such a task, the teacher signal (9) is modified to,

$\displaystyle\mathbf{y}_{\mathrm{teacher}}:=[y_{1};y_{2};y_{3}]^{\mathrm{T}}=[\varrho_{1};\varrho_{2};\varrho_{3}]^{\mathrm{T}},$
(11)

where $\varrho_{k}$ is the vol/vol fractions as defined in (1) in Section 2.1;
as such $\varrho_{1}+\varrho_{2}+\varrho_{2}=1$.

Figure 8: The EORC application for concentration determination task. (a) RGB
pseudo-colouring of the output of the biomimetic sensor
$\mathbf{P}(n)=[P_{1};P_{2}:P_{3}]^{\mathrm{T}}$ for different concentration
of chemical mixture shown as a ternary diagram. Note that for visual contrast
$\mathbf{P}(n)=[P_{1};P_{2};P_{3}]^{\mathrm{T}}$ has been normalised, see sub-
section 2.1. (b) Training samples of the EORC at discrete concentration
combinations. (c) RGB pseudo-colouring of the steady state output signal
$\mathbf{y}(t)$ for different concentration of chemical mixtures.

Here, the EORC is only trained to recognise the concentration of few possible
combinations. Figure 8(b) shows the training points used which are spaced by
20$\%$ within the ternary diagram. This specific task highlights the key
capability of the artificial neural network system to perform data-
interpolation and regression tasks.

To find the optimum operation of EORC, the error $\varepsilon$ is shown as a
function of intensity $I_{0}$ and gain $\gamma$ parameters in Fig. 8. It is
similar to Fig. 8(a), with the difference being that it is based on the
training at points shown in Fig. 8(b). From Fig. 8(b), the optimum operational
parameters of the EORC interestingly are observed to occur at the low-
intensity and high-gain $\gamma$ parameter along the bifurcation line
$\gamma=1/I_{0}$. As such, for this task the following parameters of the EORC
will be used, $I_{0}=0.197$, $\gamma=5$ and $N_{x}=50$.

Ternary diagram, Fig. 8(d), shows the steady state of output signal
$\mathbf{y}(\tilde{t})$ displayed using a RGB pseudo-colouring scheme. It is
important to note that no post-processing was applied to generate Fig. 8(d).
By visual inspection, one can observe three distinct regions in Fig. 8(d)
clearly indicating the presence of three underlying chemical specimens.
Furthermore, one can see the change of shades of the colour-encoding as the
concentration of the underlying chemicals are varied. In contrast with Fig.
8(a), there is almost no colour contrast for the case of small methanol
concentration, indicating that discriminating ethanol and isopropanol mixture
is a challenging task. This is perhaps understandable since, as, as seen in
Fig. 2(b), the absorption spectra $\alpha$ of pure ethanol and isopropanol are
similar.

##  4 Conclusions

This paper proposed a novel fully bio-inspired apparatus to monitor and
determine the concentration of chemical mixture. It implements a human eyes-
like sensor and neuromorphic optical system based on electro-optical
implementation of reservoir computer (EORC) with the single physical node
delay line and Mach-Zehnder nonlinear element for information processing.
Here, a delay differential equation including a realistic bandpass filtering
effect was used to demonstrate this sensing application. Two kinds of sensing
application were simulated successfully, namely, to discriminate known
(trained) samples of pure chemicals and to determine the concentration of
mixture chemical. The latter application highlights the true capabilities of
EORC as member of artificial neural network family to perform interpolation
and regression task. Experimental demonstration is currently in progress and
soon will be reported. Time-stepping algorithms based on
$\mathcal{Z}$-bilinear transform technique have been developed to numerically
solved the governing equations. Further optimisation that was not explored
here is the optimisation of the biomimetic sensor, for example the number of
channels and the wavelength of the filter, we believe such an optimisation
will further increase the accuracy of the whole system.

##  5 Appendices

###  5.1 Digital filter model of (6 and 7) by $\mathcal{Z}$-bilinear transform
technique

Low-pass digital filter: Performing Laplace transform on the DDE (6),

$\displaystyle\mathcal{L}[h(t)]=\frac{v(s)}{r(s)}=H(s)=\frac{1}{1+\delta s}.$
(12)

Performing $\mathcal{Z}$-bilinear transformation $s\leftarrow\frac{2}{\Delta
T}\left(\frac{1-z^{-1}}{1+z^{-1}}\right)$ on (12), it can be shown

$\displaystyle v(z)=K_{1}r(z)+z^{-1}\\{K_{1}K_{2}v(z)+K_{1}r(z)\\},$ (13)

where

$\displaystyle K_{1}=\left(\frac{2\delta}{\Delta T}+1\right)^{-1};\qquad
K_{2}=\frac{2\delta}{\Delta T}-1,$ (14)

and $\Delta T$ denote the numerical accuracy of the model, throughout this
paper $\Delta T=\tau/N_{x}$, with $N_{x}$ denote the number of virtual nodes.

Band-pass digital filter: Performing Laplace transform on the DDE (7),

$\displaystyle\mathcal{L}[h(t)]=\frac{v(s)}{r(s)}=H(s)=\frac{\frac{2m}{\omega_{0}}s}{1+\frac{2m}{\omega_{0}}s+\frac{1}{\omega_{0}^{2}}s^{2}}.$
(15)

Performing $\mathcal{Z}$-bilinear transformation $s\leftarrow\frac{2}{\Delta
T}\left(\frac{1-z^{-1}}{1+z^{-1}}\right)$ on (15), it can be shown

$\displaystyle
v(z)=K_{1}K_{2}r(z)+z^{-1}\left\\{K_{2}K_{3}v(z)+z^{-1}\\{K_{2}K_{4}v(z)-K_{1}K_{2}r(z)\\}\right\\},$
(16)

where

$\displaystyle\begin{split}&K_{1}=\frac{4m}{\omega_{0}\Delta T};\quad
K_{2}=\left(1+\frac{4m}{\omega_{0}\Delta T}+\frac{4}{\omega_{0}^{2}\Delta
T^{2}}\right)^{-1};\\\ &K_{3}=-2+\frac{8}{\omega_{0}^{2}\Delta T^{2}};\quad
K_{4}=-1+\frac{4m}{\omega_{0}\Delta T}-\frac{4}{\omega_{0}^{2}\Delta
T^{2}},\end{split}$ (17)

and $\Delta T$ denote the numerical accuracy of the model, throughout this
paper $\Delta T=\tau/N_{x}$, with $N_{x}$ denote the number of virtual nodes.

Depending on the value of parameter $m$, (15) can model either a narrow or
broad band-pass filter [45]. With $m<1$, (15) has complex eigenvalues
exhibiting a high-Q factor oscillator, as such operates as a narrow band-pass
filter centred at $\omega_{0}$. For $m>1$, the eigenvalues of (15) are real
and negative valued, such system is a strongly damped oscillator. The transfer
function (15) with $m>1$ exhibits a broad band-pass filter effect with low-
frequency cut-off (-3 dB), $f_{L}=\omega_{0}/(4\pi m)$ and the high-frequency
cut-off (-3 dB) $f_{H}=(m\omega_{0})/\pi$, from which filter parameters
$m=\frac{1}{2}\sqrt{f_{H}/f_{L}}$ and $\omega_{0}=2\pi\sqrt{f_{L}f_{H}}$. The
digital filter circuit of (13) for the low-pass and (16) for the bandpass
filter is depicted in Fig. 9(a and b), respectively.

Figure 9: Digital filter design for the low-pass filter (a) and the bandpass
(b) filter used to model DDEs (6) and (7), respectively.

###  5.2 Stationary and bifurcation point of (7) at the Quad operational
condition

Taking an adiabatic limit at the steady state reduces (7) to

$\displaystyle
x(t+\tau)=v(t+\tau)+I_{0}\cos^{2}\Phi_{0}=I_{0}\cos^{2}\left(\gamma
v(t)+\Phi_{0}\right).$ (18)

For specific ‘Quad’ operational point at positive slope of (5), one can
calculate the ‘fixed point’ $\bar{v}(I_{0},\Phi_{0}=-\pi/4)$ by solving,

$\displaystyle\bar{v}=I_{0}\cos^{2}\left(\gamma\bar{v}-\frac{\pi}{4}\right)-\frac{I_{0}}{2}.$
(19)

In the presence of band-pass filter, the steady state of signal $v(t)$ will
has a mean value of zeros. Thus, one can expand the non-linear term in its
series form, from which the stationary point of (18) can be approximated by,

$\displaystyle\bar{v}=I_{0}\left(\frac{1}{2}+\gamma\bar{v}-\frac{2}{3}\gamma^{3}\bar{v}^{3}\right)-\frac{I_{0}}{2}\rightarrow\bar{v}=\left\\{0;\pm\frac{3}{2\gamma}\sqrt{\gamma-\frac{1}{I_{0}}}\right\\}.$
(20)

Substituting (20) back to (18),

$\displaystyle\bar{x}=\left\\{\frac{I_{0}}{2};\frac{I_{0}}{2}\pm\frac{3}{2\gamma}\sqrt{\gamma-\frac{1}{I_{0}}}\right\\}.$
(21)

Specifically, the bifurcation point occurs at,

$\displaystyle\gamma=\frac{1}{I_{0}}.$ (22)

The analytical solution for the general case of (7) can be found in [56].

Funding D.F. acknowledges the support of the Engineering and Physical Sciences
Research Council (EP/P013708/1, EP/T010762/1).

Disclosures The authors declare no conflicts of interest.

Data availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

## References

  * [1] W. B. Levy and V. G. Calvert, “Computation in the human cerebral cortex uses less than 0.2 watts yet this great expense is optimal when considering communication costs,” BioRxiv (2020).
  * [2] C. Wilson, “Your brain is like 100 billion mini-computers all working together,” New Scientist 18 (2018).
  * [3] L. Beaulieu-Laroche, E. H. Toloza, M.-S. van der Goes, M. Lafourcade, D. Barnagian, Z. M. Williams, E. N. Eskandar, M. P. Frosch, S. S. Cash, and M. T. Harnett, “Enhanced dendritic compartmentalization in human cortical neurons,” Cell 175, 643–651.e14 (2018).
  * [4] P. R. Prucnal, B. J. Shastri, and M. C. Teich, _Neuromorphic Photonics_ (CRC Press, 2017).
  * [5] B. J. Shastri, A. N. Tait, T. F. de Lima, M. A. Nahmias, H.-T. Peng, and P. R. Prucnal, “Principles of neuromorphic photonics,” arXiv preprint arXiv:1801.00016 (2017).
  * [6] T. F. De Lima, B. J. Shastri, A. N. Tait, M. A. Nahmias, and P. R. Prucnal, “Progress in neuromorphic photonics,” Nanophotonics 6, 577–599 (2017).
  * [7] G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, “Recent advances in physical reservoir computing: A review,” Neural Networks 115, 100–123 (2019).
  * [8] D. V. Christensen, R. Dittmann, B. Linares-Barranco, A. Sebastian, M. L. Gallo, A. Redaelli, S. Slesazeck, T. Mikolajick, S. Spiga, S. Menzel _et al._ , “2021 roadmap on neuromorphic computing and engineering,” arXiv preprint arXiv:2105.05956 (2021).
  * [9] K. Vandoorne, J. Dambre, D. Verstraeten, B. Schrauwen, and P. Bienstman, “Parallel reservoir computing using optical amplifiers,” IEEE transactions on neural networks 22, 1469–1481 (2011).
  * [10] K. Vandoorne, W. Dierckx, B. Schrauwen, D. Verstraeten, R. Baets, P. Bienstman, and J. V. Campenhout, “Toward optical signal processing using photonic reservoir computing,” Optics Express 16, 11182 (2008).
  * [11] F. Laporte, A. Katumba, J. Dambre, and P. Bienstman, “Numerical demonstration of neuromorphic computing with photonic crystal cavities,” Optics Express 26, 7955 (2018).
  * [12] S. Phang, P. Bienstman, A. B. Seddon, C. J. Mellor, and T. M. Benson, “A chaotic microresonator structure for an optical implementation of an artificial neural network,” in _Integrated Optics: Devices, Materials, and Technologies XXIV,_ S. M. García-Blanco and P. Cheben, eds. (SPIE, 2020).
  * [13] S. Phang, P. D. Sewell, A. Vukovic, and T. M. Benson, “The optical reservoir computer: a new approach to a programmable integrated optics system based on an artificial neural network,” in _Integrated Optics Volume 2: Characterization, devices and applications,_ (Institution of Engineering and Technology, 2020), pp. 361–380.
  * [14] L. Appeltant, M. C. Soriano, G. Van der Sande, J. Danckaert, S. Massar, J. Dambre, B. Schrauwen, C. R. Mirasso, and I. Fischer, “Information processing using a single dynamical node as complex system,” Nature communications 2, 1–6 (2011).
  * [15] Y. Paquot, F. Duport, A. Smerieri, J. Dambre, B. Schrauwen, M. Haelterman, and S. Massar, “Optoelectronic reservoir computing,” Scientific Reports 2 (2012).
  * [16] T. Wu, A. J. Dufford, M.-A. Mackie, L. J. Egan, and J. Fan, “The capacity of cognitive control estimated from a perceptual decision making task,” Scientific Reports 6 (2016).
  * [17] S. Phang, D. Furniss, C. Mellor, G. Roelkens, A. B. Seddon, P. Bientsman, and T. M. Benson, “Neuromorphic sensing via temporal signal signature processed by photonic reservoir computer,” in _Optical Biopsy XIX: Toward Real-Time Spectroscopic Imaging and Diagnosis,_ R. R. Alfano, S. G. Demos, and A. B. Seddon, eds. (SPIE, 2021).
  * [18] K. J. Major, K. J. Ewing, M. K. Poutous, J. S. Sanghera, and I. D. Aggarwal, “Filter-based chemical sensors for hazardous materials,” in _Detection and Sensing of Mines, Explosive Objects, and Obscured Targets XIX,_ S. S. Bishop and J. C. Isaacs, eds. (SPIE, 2014).
  * [19] K. J. Major, M. K. Poutous, K. F. Dunnill, P. C. Deguzman, J. S. Sanghera, I. D. Aggarwal, and K. J. Ewing, “Biomimetic optical-filter detection system for discrimination of infrared chemical signatures,” Analytical Chemistry 88, 11491–11497 (2016).
  * [20] K. J. Major, M. K. Poutous, I. D. Aggarwal, J. S. Sanghera, and K. J. Ewing, “Analytical procedure to assess the performance characteristics of a non-spectroscopic infrared optical sensor for discrimination of chemical vapors,” Applied Optics 57, 8903 (2018).
  * [21] K. J. Ewing, K. Major, and J. S. Sanghera, “A biomimetic optical approach to skin cancer detection,” in _Optical Fibers and Sensors for Medical Diagnostics and Treatment Applications XX,_ I. Gannot, ed. (SPIE, 2020).
  * [22] M. Flinkman, H. Laamanen, P. Vahimaa, and M. Hauta-Kasari, “Number of colors generated by smooth nonfluorescent reflectance spectra,” Journal of the Optical Society of America A 29, 2566 (2012).
  * [23] J. M. M. Linhares, P. D. Pinto, and S. M. C. Nascimento, “The number of discernible colors in natural scenes,” Journal of the Optical Society of America A 25, 2918 (2008).
  * [24] J. L. Schnapf, T. W. Kraft, and D. A. Baylor, “Spectral sensitivity of human cone photoreceptors,” Nature 325, 439–441 (1987).
  * [25] H. H. Thoen, M. J. How, T.-H. Chiou, and J. Marshall, “A different form of color vision in mantis shrimp,” Science 343, 411–413 (2014).
  * [26] R. L. D. Valois and G. H. Jacobs, “Primate color vision,” Science 162, 533–540 (1968).
  * [27] M. K. Poutous, I. D. Aggarwal, K. J. Major, J. S. Sanghera, and K. Ewing, “Comparative discrimination spectral detection system and method for the identification of chemicals with overlapping spectral signatures,” (2018). US Patent 9,857,295.
  * [28] Infrarec, “FPI Detectors,” https://www.infratec.co.uk/downloads/en/sensor-division/flyer/infratec-fpi-b-en-mail.pdf. Accessed: 2021-11-25.
  * [29] F. Borondics, M. Jossent, C. Sandt, L. Lavoute, D. Gaponov, A. Hideur, P. Dumas, and S. Février, “Supercontinuum-based fourier transform infrared spectromicroscopy,” Optica 5, 378 (2018).
  * [30] T. L. Myers, R. G. Tonkyn, T. O. Danby, M. S. Taubman, B. E. Bernacki, J. C. Birnbaum, S. W. Sharpe, and T. J. Johnson, “Accurate measurement of the optical constants n and k for a series of 57 inorganic and organic liquids for optical modeling and detection,” Applied spectroscopy 72, 535–550 (2018).
  * [31] H. A. Lorentz, _The theory of electrons and its applications to the phenomena of light and radiant heat_ , vol. 29 (GE Stechert & Company, 1916).
  * [32] N. An, B. Zhuang, M. Li, Y. Lu, and Z.-G. Wang, “Combined theoretical and experimental study of refractive indices of water–acetonitrile–salt systems,” The Journal of Physical Chemistry B 119, 10701–10709 (2015).
  * [33] T. G. Mayerhöfer, S. Pahlow, and J. Popp, “The bouguer-beer-lambert law: Shining light on the obscure,” ChemPhysChem 21, 2029–2046 (2020).
  * [34] L. Sojka, Z. Tang, D. Jayasuriya, M. Shen, D. Furniss, E. Barney, T. M. Benson, A. B. Seddon, and S. Sujecki, “Ultra-broadband mid-infrared emission from a pr3$+$/dy3$+$ co-doped selenide-chalcogenide glass fiber spectrally shaped by varying the pumping arrangement [invited],” Optical Materials Express 9, 2291 (2019).
  * [35] C. R. Petersen, U. Møller, I. Kubat, B. Zhou, S. Dupont, J. Ramsay, T. Benson, S. Sujecki, N. Abdel-Moneim, Z. Tang _et al._ , “Mid-infrared supercontinuum covering the 1.4–13.3 $\mu$m molecular fingerprint region using ultra-high na chalcogenide step-index fibre,” Nature Photonics 8, 830–834 (2014).
  * [36] T. Sylvestre, E. Genier, A. N. Ghosh, P. Bowen, G. Genty, J. Troles, A. Mussot, A. C. Peacock, M. Klimczak, A. M. Heidt, J. C. Travers, O. Bang, and J. M. Dudley, “Recent advances in supercontinuum generation in specialty optical fibers [invited],” Journal of the Optical Society of America B 38, F90 (2021).
  * [37] N. Ganesh, “Rebooting neuromorphic design-a complexity engineering approach,” in _2020 International Conference on Rebooting Computing (ICRC),_ (IEEE, 2020), pp. 80–89.
  * [38] G. V. der Sande, D. Brunner, and M. C. Soriano, “Advances in photonic reservoir computing,” Nanophotonics 6, 561–576 (2017).
  * [39] M. Lukoševičius and H. Jaeger, “Reservoir computing approaches to recurrent neural network training,” Computer Science Review 3, 127–149 (2009).
  * [40] A. Argyris, J. Bueno, and I. Fischer, “Photonic machine learning implementation for signal recovery in optical communications,” Scientific reports 8, 1–13 (2018).
  * [41] K. Ikeda, “Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system,” Optics Communications 30, 257–261 (1979).
  * [42] K. Ikeda, H. Daido, and O. Akimoto, “Optical turbulence: Chaotic behavior of transmitted light from a ring cavity,” Physical Review Letters 45, 709–712 (1980).
  * [43] J. C. Sprott and J. C. Sprott, _Chaos and time-series analysis_ (Oxford University Press, 2001).
  * [44] B. E. Saleh and M. C. Teich, _Fundamentals of photonics_ (John Wiley & sons, 2019).
  * [45] L. Larger, “Complexity in electro-optic delay dynamics: modelling, design and applications,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371, 20120464 (2013).
  * [46] G&H, “High speed detectors,” https://www.gandh.com/wp-content/uploads/2021/10/GH-DS-FO-High-Speed-Detector-EM169-7010-Rev5-1.pdf. Accessed: 2021-11-25.
  * [47] iXblue, “DR-AN-20-MO 20 GHz Analog Driver,” https://photonics.ixblue.com/sites/default/files/2021-02/DR-AN-20-MO.pdf. Accessed: 2021-11-25.
  * [48] S. Phang, T. M. Benson, H. Susanto, S. C. Creagh, G. Gradoni, P. D. Sewell, and A. Vukovic, “Theory and numerical modelling of parity-time symmetric structures in photonics: Introduction and grating structures in one dimension,” in _Recent Trends in Computational Photonics,_ (SpringerNature, 2017), pp. 161–205.
  * [49] S. Phang, A. Vukovic, G. Gradoni, P. D. Sewell, T. M. Benson, and S. C. Creagh, “Theory and numerical modelling of parity-time symmetric structures in photonics: Boundary integral equation for coupled microresonator structures,” in _Recent Trends in Computational Photonics,_ (SpringerNature, 2017), pp. 207–233.
  * [50] S. Phang, A. Vukovic, H. Susanto, T. M. Benson, and P. Sewell, “Impact of dispersive and saturable gain/loss on bistability of nonlinear parity–time bragg gratings,” Optics Letters 39, 2603 (2014).
  * [51] S. Phang, “Theory and numerical modelling of parity-time symmetric structures for photonics,” Ph.D. thesis, University of Nottingham (2016).
  * [52] C. W. Cryer, “Numerical methods for functional differential equations,” in _Delay and Functional Differential Equations and their Applications,_ (Elsevier, 1972), pp. 17–101.
  * [53] A. Bellen and M. Zennaro, _Numerical methods for delay differential equations_ (Oxford university press, 2013).
  * [54] R. Reed and R. J. MarksII, _Neural smithing: supervised learning in feedforward artificial neural networks_ (Mit Press, 1999).
  * [55] V. Janyani, A. Vukovic, J. D. Paul, P. Sewell, and T. M. Benson, “Time domain simulation in photonics: A comparison of nonlinear dispersive polarisation models,” Optical and Quantum Electronics 37, 3–24 (2005).
  * [56] L. Weicker, T. Erneux, O. D'Huys, J. Danckaert, M. Jacquot, Y. Chembo, and L. Larger, “Slow–fast dynamics of a time-delayed electro-optic oscillator,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371, 20120459 (2013).
  * [57] L. Weicker, T. Erneux, O. d’Huys, J. Danckaert, M. Jacquot, Y. Chembo, and L. Larger, “Strongly asymmetric square waves in a time-delayed system,” Physical Review E 86, 055201 (2012).