Neuromorphic photonic circuit modeling in Verilog-A

The architecture of conventional computers is based on a separate memory and processor, which is not well suited for massively parallel and distributed computational models such as neural networks for machine learning.¹ The rise in machine learning applications and the availability of larger datasets have motivated the investigation of new computing architectures that can solely be dedicated to performing complex calculations at a much higher speed than conventional computers or graphics processing units (GPUs). Neuromorphic engineering attempts to implement neural networks directly on hardware that reflects their massively distributed nature for improved computational efficiency. Recent advances in integrated photonic neuromorphic architectures promise to deliver computations with higher speed (higher bandwidth and lower latency) than electronics by exploiting the parallel nature of light through wavelength-division multiplexing (WDM) and the passive nature of optical waveguides.^1,2

For a successful experimental demonstration of a large-scale neuromorphic photonic integrated circuit (PIC), it is imperative to first simulate the PIC by capturing the underlying physics of the individual devices and their interactions and predict the system behavior in the presence of an external stimulus, namely, electrical and optical signals. Before fabricating the PIC, the circuit-level simulations need to predict the system response.

PICs have both photonic and electronic devices, which are heterogeneously or monolithically integrated, including packaging components (e.g., wire bonds, transmission lines, and biasing circuits). Therefore, it is critical to perform an electro-optical co-simulation that could also simulate non-ideal effects, such as parasitics and noise. Different kinds of commercial software are available to simulate photonic devices at the physical level such as Ansys Lumerical FDTD,³ MEEP,⁴ MODE,⁵ etc. Most of these photonic software are based on finite-difference time-domain (FDTD) methods that allow the maximum of detail to be captured, but are time-consuming and computationally expensive when dealing with larger scale networks. Circuit level simulators, such as Ansys Lumerical INTERCONNECT,³ primarily offer photonic device modeling with minimal electronics. INTERCONNECT offers interoperability with SPICE-level simulators using the Cadence-INTERCONNECT co-simulation environment, where the electrical circuit is simulated in Cadence Spectre and the optical circuit is simulated in INTERCONNECT. In every time step, Spectre and INTERCONNECT communicate data bidirectionally between the electrical and optical domains. This makes it significantly slower than an electronic design automation (EDA)-only simulation with Verilog-A models for the photonic components. When working with feedback systems, such as continuous-time recurrent neural networks (CTRNNs), and reservoir computing,⁶ where the outputs depend on the history of inputs and time dynamics plays an essential role in achieving a specific task such as quadratic programming (QP),⁷ chaos generation,⁸ long short-term memory (LSTM),⁹ etc., it becomes crucial to perform any co-simulation on the same platform. Instead of dealing with complex FDTD calculations, an efficient approach is to work with compact models (accurate, computationally efficient, and easy for parametric extraction) of photonic devices, which can provide maximum detail and incorporate the physics behind the devices.

In this paper, we adopt a Verilog-A based approach^10–17 that is capable of performing electro-optical co-simulation for neuromorphic architectures. Verilog-A is a time-domain modeling language for analog circuits. It inherently supports the electrical components and is easily scalable to perform simulations with larger networks. Representing photonic devices would require compact models of each device presented in this paper. The simulations can be performed in the general Verilog-A supported electronic design automation (EDA) tool, such as Cadence Virtuoso, which enables AC, DC, time-domain, and stability analysis on the same platform. Verilog-A is popularly used for simulating CMOS electrical models and, thus, opens the scope of CMOS-enabled neuromorphic electro-optical architecture simulations. Furthermore, this approach can be easily extended to simulate larger networks with multiple neurons and layers, such as convolutional neural networks (CNNs),¹⁸ and study hierarchical and complex photonic networks.^19,20

Artificial neural networks (ANNs) are widely used in modeling non-linear processing and performing various machine learning tasks, such as face recognition, pattern detection, anomaly detection, natural language processing (NLPs), etc.²¹ The basic building block of an ANN is a neuron. Each neuron in a network performs two operations: a linear operation consisting of a weighted sum, i.e., element-wise multiplication of inputs with weights before being summed and followed by a non-linear unit that applies an activation function to the weighted sum, yielding the output of the neuron. The neuron’s output becomes an input to the neurons in the next layer or to itself in some networks with recurrent connections. Therefore, implementing an ANN on hardware requires the isomorphism of the neuron on the hardware substrate. The physical dynamics (how system behavior changes over time) governing the hardware system should perform the neuron computations to represent a network correctly.

In photonics, the integration of different optical devices and electrical components provides a way to represent an artificial neuron. While there are several ways to achieve neuron isomorphism in photonics, we focus on photonic neurons that are implemented in silicon photonics using multi-wavelength tunable filters as synaptic weights and electro-optic nonlinearity. Specifically, the building blocks of such an integrated photonic neuron are shown in Fig. 1, including microring resonator weight bank, a balanced photodetector for weighted summing, and a microring modulator for the non-linear activation function.^8,22 The input values to the network are encoded as different wavelength laser powers. The weighting operation from one layer to another is performed by using tunable add–drop microring resonators. Each microring of the weight bank on a photonic neural network acts as a weight. Different weight values are encoded by shifting the resonance wavelength of the ring resonator. The shift in resonance wavelength is done by changing refractive index n_eff of the silicon microring resonator by applying an external effect, such as thermo-optic effect, plasma dispersion effect, absorption effects, etc. The summation operation is performed by photodetectors that output the current proportional to the sum of different wavelength laser powers. Having photodetectors in a balanced (push–pull or series) configuration and enabling microring resonators at critical coupling provides a way to encode positive and negative weight values.²⁰ The output current from balanced photodetectors is converted to a voltage by using a transimpedance amplifier (TIA). The non-linear activation function is implemented by a PN/PIN junction-based microring modulator that is driven by the output voltage from the TIA. The transfer functions of the microring modulator can represent the predominant activation functions used in neural network, such as rectified linear units (ReLU), sigmoidal, quadratic, etc.²² 

The methodology of our work is based on representing optical signal as an analytic signal by using electric field representation of light.^10,11 As Verilog-A does not support complex numbers, the optical signals propagate as real and imaginary parts of the electric field carrying information about the magnitude and phase of the electric field. The information about the wavelength of the light can be represented in two ways.^10,11,23 The first method^10,11 is based on the oscillating electric field signal where the frequency of the oscillation defines the wavelength of the light. It includes defining light source’s reference frequency (THz) and working within the bandwidth of interest. Depending on the application, the optical bandwidth can range from 100 GHz to a few THz above the reference frequency. However, working at higher frequencies (⁠$>2$ THz) requires minimal time and takes an extremely long time to run the simulation. To counter this problem, the second method²³ includes defining the wavelength of the light as an explicit signal that passes across all the devices. This enables simulating a wide bandwidth range of frequency and faster run-time. In the case of photonic neuromorphic architectures, the bandwidth of interest is within the free spectral range (FSR) of the microring resonator, which is of the order of 2–3 THz, thus making the second method suitable for the simulations. All the simulations performed in this paper are based on the second approach, where the wavelength of the light is defined as a separate signal from the light source, along with real and imaginary parts of the electric field.

Devices such as microring resonators and Mach-Zehnder interferometer (MZI) are modeled by combining different essential components such as waveguides, couplers, phase shifters, etc. Thus, we do not require any explicit model to define composite devices.¹⁰ 

The fundamental photonic devices in neuromorphic computing are lasers, waveguides, couplers, splitters, combiners, and photodetectors. The Verilog-A models for each device are based on the circuit models^{10,12–14,24} that are described next.

The laser is modeled as a continuous wave (CW) light source.^10,11 The laser model is defined by the power and wavelength of laser light, as shown in Fig. 2(a). The voltage value P at the input port represents the power of the light in Watts and is proportional to the magnitude of the electric field amplitude. The wavelength of the light is defined as an independent parameter that can be adjusted according to the simulation. The output of the laser includes three signals: (1) the real part of the electric field (E_real), (2) the imaginary part of the electric field (E_img), and (3) the wavelength of the light (E_λ). The magnitude and phase of the electric field keep changing while propagating through various optical devices according to the physical equations governing the devices, thus changing the real and imaginary electric field signals. However, the wavelength of the signal remains conserved along with all the devices. The initial phase θ_phase(0) of the electric field at the source is defined as zero at the start of the simulation,¹⁰ 

The straight waveguide is modeled by incorporating three significant effects, namely, loss due to propagation, phase shift, and time delay.^10,11 The output electric field signal from the waveguide is given by^10,25

where E_in is the input electric field signal to the waveguide, L is the length of the waveguide, c is the speed of light, n_g is the group index, α is the attenuation parameter, ϕ is the phase of the electric field signal, and ω is the angular frequency of the light signal. The implementation of complex electric field equations for the Verilog-A waveguide model is attached in the supplementary information.

Phase shifters are modeled as a waveguide where the effective refractive index n_eff and attenuation α_eff depend on the system’s perturbation, which can be due to thermo-optic effect, plasma-dispersion effect, carrier-depletion, etc. In that case, the output of the electric field from the phase shifter becomes¹⁰ 

where n_eff and α_eff can be modeled by extracting polynomial-fit coefficients^26,27 from experimental or device level simulation data such that n_eff = a₀ + a₁V + a₂V² + a₃V³ + a₄V⁴ and α_eff = b₀ + b₁V + b₂V² + b₃V³ + b₄V⁴, with V being the voltage applied across the phase-shifter.

The coupler is modeled as a unidirectional point coupler with two input signals and two output signals as given in Fig. 2(c) and defined as¹¹ 

where t is the coupling coefficient, $Ein1$ and $Ein2$ are the input electric fields to the coupler, and $Eout1$ and $Eout2$ are the output electric fields from the coupler.

The PN junction diode describes the circuit model of the photodetector with capacitive and resistive elements as shown in Fig. 2(d).¹¹ The diode works in the reverse bias condition, and the output current is proportional to the incoming power of laser light. The current through the photodetector diode is defined by using the Shockey diode equation,

where I_s is the saturation current; V_PDbias is the voltage applied across the photodetector; and v_t is the thermal voltage defined by kT/q in which k is Boltzmann’s constant, T is the temperature, and q is the charge of an electron. The current through resistive and capacitive branches is defined by

where R is the responsivity of the photodetector in A/W, E_in is the electric field of the incoming light, C is the capacitance, and V_PDbias is the bias voltage applied across the photodetector. The present model of the photodiode is sensitive to all-optical wavelengths. Still, a filter can be implemented in the Verilog-A code to work in a specific bandwidth of wavelengths.

Microring resonator is modeled by combining discrete components that include couplers, waveguides, and phase shifters.^10,11 Figures 3 and 4 show the circuit schematic and the corresponding spectrum of all-pass and add–drop ring resonators implemented in Verilog-A by combining couplers and waveguides, respectively.^10,11 The spectrum is obtained by performing a DC sweep simulation over the wavelength of the light source. The spectrum shows the expected Lorentzian output of the microring resonator for both thru and drop ports.

The photonic device model parameters, such as effective refractive index n_eff, attenuation co-efficient α_eff, Δn_eff, etc., are extracted by characterizing each of the on-chip device fabricated at Advanced Micro Foundry (AMF). This section provides information about the methods used to extract and tune the parameters from the measured data. In addition to it, Verilog-A simulations are compared with the experimental results by incorporating the extracted parameters.

One important component in photonic neuromorphic circuits is the weight bank, which includes an array of microring resonators (MRRs) with a common bus waveguide passing through each microring’s thru and drop ports. Figure 5(a) shows the microscopic image of an on-chip weight bank with five microrings with different radii, DC pads, and waveguides. Figure 5(b) shows the experimental setup to obtain the spectrum of the weight bank, which includes the mounted chip, DC probe for electrical probing, thermal controller to stabilize the chip temperature, and V-Groove for optical probing.

In Verilog-A, the N-MRR weight bank is created by combining the N number of microrings in series with common thru and drop ports bus waveguide. The spacing between microrings is defined by adding a straight waveguide in-between the consecutive rings at the bus waveguide as shown in Fig. 5(c). In the present models, the coupling coefficient t of the coupler is defined as a constant value that can also be modified as a polynomial function of the gap in between the ring and bus waveguide based on the device simulation or measured data. Figure 5(d) shows the comparison between the measured spectrum gathered by using an optical spectrum analyzer (Aragon BOSA 400) and the Verilog-A spectrum obtained at the thru port of the weight bank. As the spectrum wavelength is quite broad, the change in the refractive index of a silicon waveguide with wavelength n_eff(λ) has also been considered in the Verilog-A model as given by²⁸ 

where λ is the wavelength of light, λ₀ is the reference wavelength, and n₀ is the refractive index of silicon at the reference wavelength. Table I shows the values of fabricated radius based on the layout design and the values used in Verilog-A simulations. The radii of the last two rings have been modified slightly in Verilog-A so that the resonances match the experimental data. The slight change in the radii can be mainly due to the fabrication variations.

Each microring of the weight bank represents a weight, and the weighting mechanism is performed by changing the resonance wavelength. The method used to change the resonance of a ring in the weight bank is a thermo-optic effect that includes changing the temperature of the waveguide to affect the refractive index of the silicon, which leads to a change in the resonance wavelength according to the following equation:²⁹ 

where n_eff(T) is the effective refractive index that is a function of the temperature, L = 2πR is the length of the MRR waveguide of radius R, and m is the order of resonant mode.

An n-doped heater embedded in the microring alongside the waveguide is used to change the temperature of the silicon. The n-doped heater contains the interplay of thermal, optical, and photoelectric effects. The temperature affects the silicon’s refractive index, which changes the transmission at the laser wavelength. The change in transmission results in a change in photoelectric absorption, which affects the resistance of the heater. The difference in the resistance changes the power applied at a given applied current, which changes the temperature.^30,31 Thus, to model an n-doped heater accurately, each of these effects needs to be captured correctly. In our approach, the n-doped phase shifter is modeled by finding Δn_eff with respect to the applied current. Equation (11) provides a direct relationship between Δn_eff and Δλ_res as follows:

Thus, extracting Δλ_res with respect to the current applied to the n-doped heater and performing a polynomial fit on it will provide information about Δn_eff.

Figure 6 shows the experimental spectrum obtained by sweeping current (0–1 mA) to the first n-doped heater ring with radius of 8 μm. The resonance wavelength of the MRR changes with the applied current due to the thermo-optic effect.

The current sources (Keithleys 2600) are used to sweep current and are controlled digitally using the Lightlab package.³² The resonance wavelengths were extracted for the first ring, and the change in resonance wavelength is plotted against the swept current as shown in Fig. 7(a). Using Eq. (12), Δn_eff is found with a fourth order polynomial fit as follows:

such that

where I_heat is the current applied to the MRR and the extracted polynomial fit coefficients n₀, a₀, a₁, a₂, a₃, and a₄ are given in Table II and the fitted plot is shown in Fig. 7(b).

The equation n_eff(I_heat) is implemented in the Verilog-A code of a thermal phase shifter. The n-doped heater MRR is constructed by combining a coupler and a thermal phase shifter. Figure 7(c) shows the comparison between the simulation and the measured spectra of n-doped heater MRR. The resonance wavelengths in the Verilog-A spectrum for the n-doped heater match well with the results from the measured data as shown in Fig. 7(d).

Using the same methodology as given above, the PN junction MRR modulator is constructed by combining the coupler and PN junction phase shifter.^10,11,33 The electrical circuit model of the PN junction phase shifter includes the PN junction diode with the capacitive branch such that¹¹ 

where I_res represents the Shockley diode equation and I_cap is the current through the capacitive branch.³⁴ The capacitance value is based on the circuit model of a ring modulator³⁵ and can also be modeled by a polynomial fit based on the device simulation³⁶ performed in software, such as Ansys Lumerical MODE and CHARGE.

The on-chip MRR modulator exploits the plasma dispersion effect to modulate the laser light. In plasma dispersion effect, the concentration of the carriers is changed by either injecting or removing the carriers, which results in changing the refractive index of the silicon and the absorption.²⁸ The change in refractive index and absorption are phenomenologically described by the following equations:^28,37

where ΔN and ΔP are the carrier densities of electrons and holes, respectively. Thus, to model a PN junction phase shifter, the change in the refractive index and absorption needs to be derived.

To model n_eff and α_eff, the on-chip PN MRR modulator is characterized by measuring the change in resonance wavelength and absorption with respect to the applied bias voltage to the ring modulator. Figure 8(a) shows the microscopic image of the on-chip PN MRR modulator, and Fig. 8(b) shows the experimental setup including the chip, RF probe for electrical probing, and V-Groove for optical probing. Figures 8(c) and 8(e) show the observed spectrum obtained for the modulator operated in reverse and forward biases, respectively. The shift in resonance and absorption is higher in forward biasing compared to the reverse bias condition for a given bias voltage change. Equation (12) is used to find the difference in the refractive index using a shift in resonance wavelength at different bias voltages. A fourth order polynomial fit is done as shown in Figs. 8(d) and 8(f) for reverse bias and forward bias, respectively, such that

where V_bias is the reverse bias voltage applied to the modulator and the constants n₀, a₀, a₁, a₂, a₃, and a₄ are given in Tables III and IV for reverse bias and forward bias, respectively.

The attenuation coefficient is also calculated by using the fourth order polynomial fit such that

where a is the single-pass amplitude transmission that includes both propagation loss in the ring and loss in the couplers;²⁹ α_eff is the power attenuation coefficient; L is the length of the phase shifter; V_bias is the bias voltage across the PN junction; and b₀, b₁, b₂, b₃, and b₄ are the polynomial fit coefficients given in Tables III and IV for reverse bias and forward bias, respectively. The fourth order polynomial fit curves are an empirically good fit on the measured data based on the r-square values.

After obtaining n_eff(V) and α_eff(V), the polynomial fit equations are implemented in the phase shifter Verilog-A model. The PN junction MRR modulator is constructed by combining coupler and PN junction phase shifter and optical spectrum simulations are performed by biasing the PN phase shifter at different voltages. Figures 9(a) and 9(c) show the comparison between the simulated and the measured spectra obtained for the PN junction modulator in reverse bias and forward bias, respectively. The resonance wavelengths for the simulations are consistent with the experimental results as shown in Figs. 9(b) and 9(d) for different biases. The loss in the simulated spectrum is deviated from the measured spectrum because the measured spectrum is not symmetric about the resonance wavelength. The asymmetry in the ring modulator spectrum can occur due to several reasons, such as nonlinear interactions,³⁸ optical bistability,³⁹ fabrication imperfections in the directional coupler of the ring, or the back-reflections from outside of the ring.⁴⁰ The present Verilog-A models assume the ideal conditions and do not account for these imperfections and interactions.

The photonic neuron simulation includes implementing the activation functions and integration of laser, microring resonator, photodetector, microring modulator, and electrical components on the same platform. Verilog-A supports the integration of these devices for neuromorphic simulations.

The modulator neuron consists of two photodetectors (PDs) in a balanced configuration connected electrically to the driving input of the MRR modulator as shown in Fig. 10. The two PDs are reversed biased and sense two multiplexed incoming signals (P+ and P−). The balanced PD (BPD) outputs a current i_BPD that is proportional to the difference of the power of the incoming signals such that

where i₊ is the photocurrent from the positive PD (proportional to P+) that adds to the injected current and i₋ (proportional to P−) is the photocurrent from the negative PD that diverts the current coming from the bias (i_bias). The output photocurrent from BPD combines with the incoming bias current i_bias. The total current (i₊, i₋, i_bias) passes through the resistor-based transimpedance amplifier (TIA). The TIA output voltage drives the PN based MRR modulator that exploits the plasma-dispersion effect to change the refractive index of the silicon microring, which, in turn, results in a nonlinear transfer on the power of the input pump laser (P_pump). The plasma dispersion effect provides high-speed modulation of the MRR resonance, but it is a weak effect. In some cases, the resonance of the MRR needs to be significantly tuned so that the transfer function of the MRR modulator corresponds to the required activation function. A heater element (n-doped heater or a resister) is also fabricated alongside the MRR that exploits the thermo-optic effect to tune the resonance wavelength. The thermo-optic effect is a strong effect in changing the resonance wavelength, but it has a slow modulation response. Therefore, the plasma dispersion effect is used in the case of high-speed tuning, and the thermo-optic effect is used in the case of significant resonance wavelength tuning.

Figure 11 shows the wavelength spectra of the MRR modulator at different combinations of heater current (i_heat) and bias current/voltage (i_bias/V_bias) for the PN junction modulator simulated in Verilog-A. The balanced photodetector has two laser inputs (IN+ and IN−) corresponding to the laser connected to the upper PD and the lower PD, respectively. The different colors of the spectra correspond to the cases when only one of the lasers is on (IN + or IN−), both lasers are on, and both lasers are off. The heat current shifts the resonance to the right, whereas the bias current/voltage shifts the resonance to the left, thus enabling both red and blue shifts. It can also be noted that the change in resonance is much higher in the case of heater current i_h than the bias voltage V_bias.

The photonic neuron can implement several nonlinear activation functions that are important for different machine learning applications. The transfer functions of an MRR modulator are used to implement other different activation functions. Different activation functions are obtained by biasing the modulator at different voltages. Experimentally, the transfer functions of a PN modulator are obtained by sweeping the bias voltage at a constant input wavelength of laser light. Figure 12 shows the comparison between the simulated and the experimental transfer functions obtained for the PN junction based MRR modulator. The simulation and experimental data can be seen well aligned for quadratic and sigmoid shape transfer functions.

The implementation of the activation function is also tested on a time-varying input signal. A CW laser source is used to input the signal in one of the positive PD of the BPD, while the negative PD is set to be zero laser power. Figure 13 shows the response of the modulator neuron to the time-varying input carrier signal. For the rectified output [Fig. 13(d)], the modulator is biased such that the transfer function corresponds to the sigmoidal shape activation function (Fig. 13). The steepness of the linear region of the transfer function can be changed in many ways, such as increasing absorption by operating the modulator in forward bias, designing application-specific MRR, etc. The position of the linear region can also be changed by biasing at a different region of the Lorentzian peak. Similarly, for the quadratic output [Fig. 13(e)], the modulator is biased to obtain the quadratic transfer function [Fig. 13(b)]. Both the transfer functions are biased at the starting point (x = 0.0) of the x axis in Figs. 13(a) and 13(b) to obtain Figs. 13(d) and 13(e), respectively.

A single photonic neuron is simulated in Verilog-A to implement input, linear weighting, and non-linear activation function to show fully integrated photonic neuron implementation. Figure 14 shows the schematic of the single photonic neuron with optical and electrical parasitics. The laser power (P_in) represents the input (x_i) to the neuron, the current (I_heat) to the n-doped heater represents the weight (w_i), balanced photodetectors perform the summation operation, and the transfer function of the modulator represents an activation function (σ). R_TIA acts as a passive transimpedance amplifier and the circuit parameters such as parasitic resistance (R_par) are also included based on the literature.³⁵ The output power of the laser is set constant by providing a constant voltage to the laser. The wavelength of the laser is set according to the resonance wavelength of n-doped MRR. The weighting mechanism [−1,1] is performed by applying current I_heat to the n-doped MRR. The modulator is biased by V_bias voltage so that it operates in the required region of the transfer function. Sweeping the current at the n-doped heater changes the resonance wavelength of the MRR, which results in changing the voltage at the modulator as shown in Fig. 14. The voltage at the input of the modulator drives the modulator to implement the non-linear activation function on the input power of the modulator laser. By changing the bias voltage, the modulator can be operated to implement different transfer functions corresponding to the neural network’s activation functions.

In this paper, we adopted a Verilog-A based approach for simulating neuromorphic architectures by modeling lasers, waveguides, couplers, photodetectors, and phase shifters. The parameters for the photonic devices are extracted from the experiments and implemented in the Verilog-A models. The simulation results for the composite devices, such as MRR weight bank and MRR modulator, are found to be consistent with the experimental results, thus validating the Verilog-A models. We have also shown the characteristics of the MRR modulator by showing the implementation of transfer functions and response to different bias voltages. Furthermore, the electrical circuit parameters and parasitics are also included in the simulation of the single photonic neuron. The work done in this paper would significantly increase the capabilities to perform electro-optical co-simulations on the same platform primarily dedicated toward photonic neuromorphic computing architectures.

This work can be extended to perform large-scale neuromorphic photonic circuit simulations. A Python package can be made to automate the netlist generation and interface it with general python-based machine learning training packages, such as Tensorflow,⁴¹ Keras,⁴² and ONNX.⁴³ 

See the supplementary material for implementing the complex electric field equations in the waveguide Verilog-A model.

The authors would like to thank CMC Microsystems for facilitating access to the CAD tools. This work was supported, in part, by the Natural Sciences and Engineering Research Council of Canada (NSERC).

The authors have no conflicts to disclose.

The Verilog-A codes for the photonic devices and the simulation data are available from the corresponding author upon reasonable request.