SIMPEL: Circuit model for photonic spike processing laser neurons

We propose an equivalent circuit model for photonic spike processing laser neurons with an embedded saturable absorber---a simulation model for photonic excitable lasers (SIMPEL). We show that by mapping the laser neuron rate equations into a circuit model, SPICE analysis can be used as an efficient and accurate engine for numerical calculations, capable of generalization to a variety of different laser neuron types found in literature. The development of this model parallels the Hodgkin--Huxley model of neuron biophysics, a circuit framework which brought efficiency, modularity, and generalizability to the study of neural dynamics. We employ the model to study various signal-processing effects such as excitability with excitatory and inhibitory pulses, binary all-or-nothing response, and bistable dynamics.


Introduction
Photonics has recently witnessed [1][2][3][4][5][6][7][8][9][10][11][12] a deeply committed exploration of neuro-inspired unconventional computing paradigms that promise to outperform conventional technology in certain problem domains. This field of neuromorphic engineering [13][14][15] aims to build machines that better interact with natural environments by applying the circuit and system principles of neuronal computation, including robust analog signaling, physics-based dynamics, distributed complexity, and learning. Cognitive computing platforms [16,17] inspired by the architecture of the brain promise potent advantages in efficiency, fault tolerance and adaptability over von Neumann architectures for tasks involving pattern analysis, decision making, optimization, learning, and real-time control of many-sensor, many-actuator systems. These neural-inspired systems are typified by a set of computational principles, including hybrid analog-digital signal representations [18], co-location of memory and processing [19], unsupervised statistical learning [20], and distributed representations of information [21].
A sparse coding scheme called spiking [22,23] has been recognized as a cortical encoding strategy [24,25] with firm code-theoretic justifications [26,27] and promises extreme improvements to computational power efficiency [14]. Since spikes are discrete events that occur at analog times, this encoding scheme represents a hybrid between traditional analog and digital approaches, capable of both expressiveness and robustness to noise [28]. This distributed, asynchronous model processes information using both space and time [21,23], is naturally robust, and is amenable to algorithms for unsupervised adaptation [19,20]. The marriage of photonics with spike processing is fundamentally enabled by the strong analogies [7] of the underlying physics between the dynamics of biological neurons and lasers; they both can be understood within the framework of dynamical systems theory, and can display the crucial property of excitability-systems that can be excited from their in a stable steady rest state to emit a spike by a super-threshold followed by a refractory period. The rate equations of photonics devices, however operate approximately eight orders of magnitude faster than biological time scales. In addition to the high switching speeds and high communication bandwidth, the low cross-talk achievable in photonics are very well suited for an ultrafast spike-based information scheme with high interconnection densities. Furthermore, the high wall-plug efficiencies of photonic devices may allow such implementations to match or eclipse equivalent electronic systems in low energy usage. Consequently, a network of photonic neurons-photonic spike processorscould access a picosecond and low-power computationally rich domain that is inaccessible by other technologies. This novel processing domain-ultrafast cognitive computing-represents a broad domain of applications where quick, temporally precise and robust systems are necessary, including: adaptive control, learning, perception, motion control, sensory processing, autonomous robotics, and cognitive processing of the radio frequency (RF) spectrum.
Despite their potential in this regard, all the excitable laser systems studied in the context of spike processing, either with the tools of bifurcation theory [8,29,30] or experimentally [5,31], have been limited to a couple of devices. This has been largely due to the lack of a dynamical simulation platform for photonic neurons. Due to the complexity associated with these systems, it is impossible to apply simplifying assumptions for purely analytic approaches which would risk not capturing all the rich dynamics associated with such systems. A platform for simulating photonics neurons must simultaneously capture their underlying physics (photoncarrier dynamics) and be amenable for studying large scale on-chip optical network architecture to support massively parallel communication between high-performance spiking laser neurons [32][33][34] (also highlighted in [35]). Furthermore, such a simulation platform would be crucial to gain critical insight into the behavior of a photonic spike processor under a variety of conditions.
The photon-carrier dynamics in a excitable laser are strongly coupled, making the simulation of transient spiking phenomena in laser neurons much more involved than steady-state simulation of typical CW lasers. Multidimensional device-level programs (for example, based on finite element method (FEM)) are too computationally intensive to accurately account for these complex photon-carrier interactions, and they are also impossible to utilize for simulating a scalable integrated optoelectronic system with hundreds or thousands of such devices. Numerical methods such as the Runge-Kutta methods can be employed for calculating solutions of the rate equations. However, they have a number of limitations for modeling practical situations; for example, they cannot easily incorporate parasitic networks for studying device interactions and determining the frequency response. Other simulation techniques take advantage of the sparse nature of spiking signals and adopt event-based models for extremely large-scale networks. While computationally efficient, event-based models must overlook many subtleties in the dynamics of individual units, which are fundamentally driven by differential equation mod-els.

Our contribution
Here, we propose SIMPEL-SImulation Model for Photonic Excitable Lasers-which bridges the gap between the underlying physics and relevant dynamics of excitable lasers by transforming its rate equations to an equivalent circuit representation. We show that this circuit model leads to a highly efficient simulation framework that expedites the computation of the rate equations, and also facilitates computer-aided design (CAD) and analysis of an integrated optoelectronic system. The circuit modeling approach captures the important dynamical behaviors of a photonic neuron in a way that is extendable to more refined studies of parasitics and networks while also being agnostic to specifics of device implementation. In our approach, the carrier density and the photon density are converted to current and voltage representations, respectively. The input current (carrier) and output voltage (light) form a two-port equivalent model. This two-port circuit can be coupled with other such circuits for an efficient and scalable simulation platform to study for example an architecture of an interconnected network of laser neurons at the system level for wide range of applications in high-performance computing, adaptive control, and learning. Our model is compatible with SPICE, the de facto industrial standard for computer-aided circuit analysis, which is hyper optimized for circuit analysis.
By recasting excitable laser dynamics in a common language of circuit analysis, SIMPEL manifests several advantageous features that will aid in the design and characterization of laser neuron devices for signal processing. This approach parallels the work of Hodgkin and Huxley [36], who developed an equivalent circuit model for biological neuron dynamics in order to express neuron behavior in a common language of engineering that can abstract away underlying biochemical mechanisms. A circuit analysis framework is a modular abstraction, which is important for incorporating second-order effects, such as parasitics, and for studying small networks of interconnected laser devices. Compatibility with the SPICE engine means that the simulation is highly optimized. The resulting increase in speed over physical and Runge-Kutta methods is particularly relevant for Monte-Carlo studies of noise in which stochastic effects are studied with many calls of the same simulation. Finally, SIMPEL's phenomenological (a.k.a. behavioral) nature makes it applicable to different physical implementations of an excitable laser with embedded saturable absorber (SA) without resorting to full-scale device simulation. Although the model and parameters must be derived from the underlying optoelectronic mechanisms, SIMPEL can more readily grant insight about the relationship between these parameters and behavioral dynamics relevant for signal processing, such as energy usage and noise. One example of an application is parameter fitting data measured from a laser neuron circuit. In this case, the effects of noise and parasitics must be accounted for, and the internal physical mechanism is not observable. SIMPEL could allow a user to learn the most about circuit operation by improving control over the fit degrees-of-freedom and, more generally, over the synthesis of a priori knowledge with observed behavior.  The LIF model is one of the most ubiquitous models in computational neuroscience and is the simplest known model for spike processing [37]. The dendrite tree collects N inputs which represent induced currents in input synapses x j that are continuous time series consisting either of spikes or continuous analog values. Each input is independently weighted by ω j , which can be positive or negative, and delayed by τ j resulting in a time series that is spatially summed. This aggregate input induces an electrical current, Fig. 1. Schematic of a biological neuron and a two-section excitable laser that share key dynamical properties. In the LIF neuron model, weighted and delayed input signals are spatially summed at the dendritic tree into an input current, which travel to the soma and perturb the internal state variable, the voltage. The soma performs integration and then applies a threshold to make a spike or no-spike decision. After a spike is released, the voltage is reset. The resulting spike is sent to other neurons in the network. The excitable laser is composed of a gain section, SA, and mirrors for cavity feedback. The inputs selectively perturb the gain optically or electrically. The gain medium acts as a temporal integrator while the SA acts as a threshold detector; it extracts most of the stored energy from the gain medium into the optical mode. These dynamics emulate excitability, one of the most critical properties of a spiking neuron.
where the membrane potential V m (t), the voltage difference across their membrane, acts as the primary internal (activation) state variable. The weights and delays determine the dynamics of network, providing a way of programming a neuromorphic system. The soma acts as a firstorder low-pass filter or a leaky integrator, with the integration time constant τ m = R m C m that determines the exponential decay rate of the impulse response function and where R m and C m model the resistance and capacitance associated with the membrane, respectively. The leakage current through R m drives the membrane voltage V m (t) to 0, but an active membrane pumping current counteracts it and maintains a resting membrane voltage at a value of where t 0 is the last time the neuron spiked. Finally, the axon carries an action potential, or spike, to other neurons in the network when the integrated signal exceeds a threshold; that is, if V m (t) ≥ V thresh , then the neuron outputs a spike and V m (t) is set to V reset . This is followed by a relative refractory period, during which V m (t) recovers from V reset to the resting potential V L in which is difficult to induce the firing of a spike. Consequently, the output of the neuron consists of a series of spikes that occur at continuously valued times. There are three influences on V m (t): passive leakage of current, an active pumping current, and external inputs generating time-varying membrane conductance changes. Including a set of digital conditions, we arrive at a typical LIF model for an individual neuron: release a pulse and set V m (t) → V reset .

Excitable laser model
Next, we briefly summarize the recently discovered mathematical analogy between the LIF neuron model and an excitable laser composed of a gain section with an embedded SA [7,38] as illustrated in Fig. 1. The gain medium acts as a temporal integrator with a time constant that is equal to the carrier recombination lifetime. The SA extracts most of the stored energy from the gain medium into the optical mode and performs the function of a threshold detector. This gain-absorber interplay emulates one of the most critical dynamical properties of a spiking neuron-excitability. The Yamada model [39], describes the behavior of lasers with independent gain and SA sections with an approximately constant intensity profile across the cavity. We assume that the dynamics operate such that the gain is a slow variable, while the intensity and loss are both fast. This three-dimensional dynamical system can be described with the following equations: is the laser intensity, A is the bias current of the gain, B is the level of absorption, a describes the differential absorption relative to the differential gain, γ G is the relaxation rate of the gain, γ Q is the relaxation rate of the absorber, γ I is the reverse photon lifetime, and ε f (G) represents the small contributions to the intensity made by spontaneous emission, (noise term) where ε is very small.
It has been shown in [7] that, in certain parameter regimes, the behavior of the system closely approximates the spiking LIF model. Assuming that the inputs to the system cause perturbations to the gain G(t) only, and that the fast dynamics are nearly instantaneous, we can compress the behavior of this system into the following set of equations and conditions: where θ (t) represent input perturbations. The conditional statements account for the fast dynamics of the system that occur on times scales of order 1/γ I , and other various assumptionsincluding the fast Q(t) variable and operation close to threshold-assure that G thresh , G reset and the pulse amplitude remain constant. Comparing this to the LIF model, or equation (1), the analogy between the equations becomes clear. Setting the variables γ G = 1/R m C m , A = V L , θ (t) = I app (t)/R m C m , and G(t) = V m (t) shows their algebraic equivalence. Thus, the gain of the laser G(t) can be thought of as a virtual membrane voltage, the input current A as a virtual leakage voltage, etc. There is a key difference, however-both dynamical systems operate on vastly different time scales. Whereas biological neurons have time constants τ m = C m R m on order of milliseconds, carrier lifetimes of laser gain sections are typically in the nanosecond range and can go down to picosecond.

Rate equations
We start with the coupled rate equations for the photon and carrier densities for a two-section excitable laser with gain and SA regions, assuming only one longitudinal and one transverse optical mode is lasing [40]: Equations (3) and (4) relate the rate of change in the gain and SA regions' carrier concentration n χ to the injection current i χ , the carrier recombination rate, and the stimulated-emission rate. Note that the gain and SA cavities are represented by the subscript χ = a or s, respectively. The gain current term i a = I a + i ea accounts for the pump current I a and the electrical modulation i ea of the gain, whereas the SA current term i s = I s allows only for an adjustable threshold. Equation (5) relates the rate of change in photon number N ph that is common to the gain and SA regions, to photon decay, the rate of stimulated-emission, and the rate of recombination into the lasing mode. In addition, η i,χ is the current-injection efficiency, V χ is the cavity volume, q is the electron charge, τ χ is the rate of carrier recombination, Γ χ is the optical confinement factor, β is the spontaneous-emission coupling factor, B r is the bimolecular recombination term, and τ ph is the photon lifetime given as τ −1 ph = (c/n r )[α + ln(1/R 1 R 2 )/2L], where R 1 and R 2 are the cavity mirror reflectivities, L is the cavity length, and α is the internal loss of the cavity. Note that, for the active and passive regions to stay as gain and absorber regions, I a > qV a n a /τ a and I s < qV s n s /τ s , respectively. The stimulated-emission rate includes a carrier-dependent gain term g(n χ ) which can take on either a linear [41] or logarithmic form [42]. We chose the former for simplicity; that is, g(n χ ) = g χ (n χ − n 0,χ ), where g χ is the gain and SA region differential gain and loss coefficient, respectively, and n 0,χ is the optical transparency carrier density. Furthermore, these rate equations can be further generalized by adding a gain-saturation term [43], where ε χ is the phenomenolgical gain-compression factor. Finally, the laser output power P out is related to the photon number inside the cavity via: where λ is the lasing wavelength, η c is the output power coupling coefficient, h is Planck's constant, and c is the speed of light in a vacuum.

Equivalent circuit
For a given set of injection currents {I a , I s } in the gain and SA regions, operating point analysis of the excitable laser described by the rate equations (3)-(5) and output power (6) leads to four solutions. In addition to the correct nonnegative solution regime, in which the solutions for the photon density N ph , and carrier densities n a and n s , are all nonnegative when I a ≥ 0 and I s ≥ 0, there are also negative-power and a high-power regimes. It is therefore necessary to transform the carrier population density in the respective cavities n χ and the laser output power P out via the following pair of transformations, respectively [41]: n a = n eq,a exp qv a nkT (7) n s = n eq,s exp qv s nkT (8) where n eq,χ is the equilibrium carrier density, v χ is the voltage across the gain and SA region of the laser, n is a diode ideality factor (typically set to two for III-V devices [44,45]), v m is a new variable for parameterizing P out , δ is a small arbitrary constant set to 10 −60 , k is Boltzmann's constant, and T is the temperature of the excitable laser. It has been shown in [46] that these transformations eliminate the nonphysical solutions (negative-power and a high-power regime) and improve the convergence properties of the model during simulation. We model the carrier's dynamics dn χ /dt, by substituting the set of variable transformations (7)-(9), and the output power (6), into the rate equations (3) and (4). After applying the appropriate manipulations, we obtain: (10) With some additional rearrangement (10) can be written in terms of the respective cavity currents as where i T 1 with  Similarly, to model the photon dynamics dN ph /dt, we substitute the transformations (7)-(9), and the output power (6), into the rate equation (5). After applying the appropriate manipulations, we obtain With some additional rearrangements and the definition of suitable circuit elements, (18) can be written as where G r,a = τ ph Γ a g Θ a I T 1 and C ph = 2τ ph and R ph = 1Ω.
Finally, E pf transforms the node voltage v m , into the output power P out , by These equations can also be mapped directly into an equivalent laser neuron model as shown in Fig. 2. C ph and R ph help model the time-variation of the photon density under the effect of spontaneous and stimulated emission, which are accounted for by the nonlinear dependent current sources {G r,a , G r,s } and B, respectively. Finally, the nonlinear voltage source E pf , produces the laser neuron optical output power in the form of a node voltage at o.

Laser neuron device structures
As a demonstration of the utility of our circuit approach, we simulate our model using realistic parameters for the recently proposed vertical-cavity surface-emitting laser (VCSEL) photonic neuron [7] and a distributed feedback (DFB) laser photonic neuron [38]. Fig. 3 illustrates the cross sections of the VCSEL and DFB excitable lasers with an embedded SA. Despite their differences, because of the simplicity and generality of our approach, the dynamics of both models can be represented under the same theoretical framework and are two viable candidates for large, integrated photonic neural networks. Each model possesses complementary advantages and disadvantages to the other. VCSEL photonic neurons occupy small footprints, can be fabricated in large arrays allowing for massive scalability and low power use [47]. In addition, the model is amenable to a variety of different interconnection schemes: VCSELs can send signals upward to form 3D interconnects [48], can emit downward into an interconnection layer via grating couplers [49], or connect monolithically through intra-cavity holographic gratings [50]. An excitable VCSEL with an intra-cavity SA that operates using the same rate equation model described above has already been experimentally realized [51].
DFB laser photonic neurons, in contrast, emit light in the planar direction. Although they have lower wall plug efficiencies, use more power and occupy larger spatial footprints than their VCSEL counterparts, their natural affinity for waveguide coupling and lithographically defined operating wavelengths post growth makes them a strong candidate for integrating photonic neural networks on a single chip. DFB lasers can be coupled into passive waveguides in III-V materials [52], can butt couple into waveguides on other materials such as silicon [53], or can be defined within a silicon on insulator (SOI) substrate using techniques such as wafer bonding [54].

Results and discussion
We simulate our laser neuron circuit model using the HSPICE circuit simulator from Synopsys 1 . For our analysis, we consider typical VCSEL-SA and DFB-SA excitable lasers with material and geometrical parameters as given in Table 1. Fig. 4 depicts the simulation setup used to test the laser neuron equivalent circuit model in Fig. 2. In the simulation setup, the dc current sources I a and I s provide the bias conditions for gain and SA regions, respectively. The pulsed current sources summed as i ea , model the excitatory and inhibitory pulses (from other laser neurons), and modulate the gain region. The dummy load R L , enables the measurement of output power from the laser neuron circuit model. As stated earlier, in certain parameter regimes, a simple model of a single-mode laser with an SA section has been proven to be analogous to the equations governing an LIF neuron [7,38]. For such a laser, the active and passive cavities must stay as gain and absorber regions. This means that I a > qV a n a /τ a and I s < qV s n s /τ s , respectively. Fig. 5 depicts the simulated excitability of the VCSEL-SA system. This behavior resembles neural spiking behavior. The VCSEL-SA model is biased just below the threshold: gain current, I a = 2.7 mA, and SA current, I s = 0 mA. In Fig. 5(a)  proportional to its energy-gain enhancement. The first pair of pulses causes the laser to fire and emit a pulse, after which a refractory period occurs. During this period, a second pair of pulses is unable to cause the laser to fire. After several nanoseconds, the laser has recovered to its equilibrium state and a third pair of excitatory spikes causes it to fire again. In Fig. 5(b) an excitatory and inhibitory pair replaces the third pair of excitatory spikes. Inhibition, which decreases the carrier concentration within the gain region-gain depletion-cancels the excitatory activity and prevents the laser from firing. A crucial property of dynamical systems that arises out of the formation of hysteric attractors is bistability, which plays an important role in the formation of memory in processing systems [7]. Here, the system is recursive rather than feedforward, possessing a network path that contains a closed loop allowing the system to exhibit hysteresis, and it is essentially an extension of an autapse. Fig. 6(a) illustrates an excitable laser whose output is fed back to the input with a delay element. This circuit represents a test of the network's ability to handle recursive feedback. The VCSEL-SA excitable laser model is employed for this simulation with the biasing conditions for the gain and SA regions: I a = 2.7 mA (just below threshold) and I s = 0 mA. Fig. 6(b) shows the result for the simulation. An excitatory pulse input to the laser neuron at t = 10 ns, initiates the system to settle to a new attractor. The first output pulse is fed back to the input after being delayed by t delay = 4 ns, which initiates another excitatory pulse at the output. This recursive process results in a train of output pulses at fixed intervals before being deactivated by a precisely timed inhibitory pulse at t = 40 ns. This system successfully maintains the stability of self-pulsations. The system is also capable of stabilizing to other states, including those with multiple pulses or different pulse intervals and thus acts as an optical pattern buffer over longer time scales [7].
Note, since the output of the laser neuron is light represented in the circuit model as a voltage and its input is a current, a voltage-controlled current source is employed for the feedback connection. A practical implementation for interconnecting laser neurons requires a photodetector at its input which can couple the optical outputs of other laser neurons and electrically modulate the gain section of the laser neuron. We recently proposed such an excitable laser and photodetector system that can emulate both a LIF neuron and a synaptic variable, completing a computational paradigm for scalable optical computing [38]. Here, two photodetectors receive optical pulses (excitatory and inhibitory) from a network and are subtracted passively by a push-pull wire junction. The resulting photocurrent signal conducts over a short wire to modulate the laser gain section. Dynamics introduced by the photodetectors are analogous to synaptic dynamics governing the concentration of neurotransmitters in between signaling biological neurons. Future work will include a first-order low-pass filter to model the photocurrent flow in photodiodes and neural synaptic dynamics. The conversion between optical and electronic domains also restricts the propagation of optical phase noise and the need for direct wavelength conversion, thus eliminating two major barriers facing scalable optical computing [38].
We now consider the response of the excitable laser neurons to a single impulse over a range of amplitudes. The underlying feedback dynamics of the laser cavity yield an ideal allor-nothing energy transfer function that resembles a step as shown in Fig. 7(a). More than a logic level buffer that makes input assumptions, the output of the laser neuron is binary for any analog input. For neural systems, a binary all-or-nothing pulse output is critical to ensure amplitude robustness to channel noise. Other systems, including analog-to-digital converters and comparators, also rely on an all-or-nothing response. When used as a thresholder, the laser neuron has several uncommon advantages over prior photonic thresholding technologies that do not use laser cavity feedback [61,62]. Favorable properties include an all-or-nothing response, clean pulse generation regardless of input pulse shape, and low threshold amplitude less than a few mAs.
In phase space, the decision of output pulse or no pulse is determined by an attractive manifold (a.k.a separatrix), which separates the basin of attraction of the at rest steady-state from a large transient trajectory. There are no intermediate output possibilities, but the decision latency is greater when the input causes the dynamical state to fall close to the separatrix as depicted in Fig. 7(b). This near-threshold spike latency, characteristic of separatrix dynamics, is also observed in the Hodgkin-Huxley model [36]. While variable latency can lead to jitter in synchronous systems, it has also been proposed as a mechanism for converting continuous amplitude to spike-timing codes in certain neural systems such as the retina [63]. Further study will focus on the effects of noise on the energy transfer function and decision latency.

Conclusion
There has been a recent surge in unconventional computing paradigms that leverage the underlying physics of devices to breach the limitations inherent in traditional von Neumann architectures and CMOS device technologies. Inspired by biological neural networks, cognitive computers could outperform current technology in both complexity and power efficiency. The computational primitive for such a platform is the well-studied and paradigmatic spiking neuron. A photonic realization of spiking neuron dynamics-laser neurons-harnesses the highswitching speeds, wide communication bandwidth of optics, and low cross-talk achievable in photonics, making it well suited for an ultrafast spike-based information scheme.
In this paper, we have proposed an equivalent circuit model for laser neurons based on excitable lasers with an SA and direct gain injection. We show that by mapping the laser neuron rate equations into a circuit model, SPICE analysis can be used as an efficient and accurate engine for numerical calculations. Due to the strongly-coupled photon and carrier interactions, the laser neuron rate equations can be only solved using numerical methods such as the Runge-Kutta method. In the same way, the non-perturbative dynamics of neuron biophysics prevent any analytical solution to their behavior. Hodgkin and Huxley were the first to take an approach of mapping the underlying physics present in the neuron to equivalent circuit representations, in order to describe and simulate them within a universal and powerful engineering framework. We take a parallel approach in mapping the key physics of excitable lasers to an equivalent circuit, in hopes of establishing a foundation of comparable utility for laser neuron research.
By leveraging the modularity advantages of a circuit abstraction, SIMPEL will serve as a powerful tool in future studies. We have applied the reported model to different excitable laser types (i.e. VCSEL-SA and DFB-SA) and employed it to study signal-processing behaviors, including excitability with excitatory and inhibitory pulses, binary all-or-nothing response, and bistability. As a next step, we will investigate the implementation of spike-timing-dependent plasticity (STDP)-one of the most important algorithms for spike-based learning [20,64]with our model. Optical STDP operating on unprecedented time scales is potentially useful for applications including coincidence detection, sequence learning, path learning, and directional selectivity in visual response. Further work could also investigate power consumption, temperature, and noise, for studying the architecture of an interconnected network of laser neurons at the system level that could have a wide range of applications in high-performance computing, adaptive control, and RF spectrum processing.