# Phenomenological Model of
Superconducting Optoelectronic Loop Neurons

Jeffrey M. Shainline∗, Bryce A. Primavera, and Saeed Khan
National Institute of Standards and Technology
325 Broadway, Boulder, CO, USA, 80305
<EMAIL_ADDRESS>
(October 17, 2022)

###### Abstract

Superconducting optoelectronic loop neurons are a class of circuits
potentially conducive to networks for large-scale artificial cognition. These
circuits employ superconducting components including single-photon detectors,
Josephson junctions, and transformers to achieve neuromorphic functions. To
date, all simulations of loop neurons have used first-principles circuit
analysis to model the behavior of synapses, dendrites, and neurons. These
circuit models are computationally inefficient and leave opaque the
relationship between loop neurons and other complex systems. Here we introduce
a modeling framework that captures the behavior of the relevant synaptic,
dendritic, and neuronal circuits at a phenomenological level without resorting
to full circuit equations. Within this compact model, each dendrite is
discovered to obey a single nonlinear leaky-integrator ordinary differential
equation, while a neuron is modeled as a dendrite with a thresholding element
and an additional feedback mechanism for establishing a refractory period. A
synapse is modeled as a single-photon detector coupled to a dendrite, where
the response of the single-photon detector follows a closed-form expression.
We quantify the accuracy of the phenomenological model relative to circuit
simulations and find that the approach reduces computational time by a factor
of ten thousand while maintaining accuracy of one part in ten thousand. We
demonstrate the use of the model with several basic examples. The net increase
in computational efficiency enables future simulation of large networks, while
the formulation provides a connection to a large body of work in applied
mathematics, computational neuroscience, and physical systems such as spin
glasses.

##  1 Introduction

A technological platform capable of realizing networks at the scale and
complexity of the brains of intelligent organisms would be a tool of supreme
scientific utility. Neuromorphic hardware based on conventional silicon
microelectronics has a great deal to offer in this regard [1, 2, 3, 4, 5]. Yet
challenges remain, primarily concerning bottlenecks in the shared
communication infrastructure that must be employed to emulate the connectivity
of biological neurons. Alternative hardware may bring new benefits, and we
have argued elsewhere for the advantages of a superconducting optoelectronic
approach [6, 7, 8, 9, 10, 11]. In brief, light for communication enables high
fan-out with low latency across spatial scales from a chip to a many-wafer
system. Superconducting electronics provide single-photon detection coupled to
high-speed, low-energy analog neuromorphic computational primitives.

While the components of these superconducting optoelectronic networks (SOENs)
have been demonstrated [12, 13, 14, 15], full neurons have not. Prior to
undertaking the effort and expense of realizing the requisite semiconductor-
superconductor-photonic fabrication process, it is prudent to gain confidence
that SOENs are indeed ripe for further investigation. This confidence can be
gained through simulations of device, circuit, and system behavior using
numerical simulations on digital computers. The constituent devices are most
commonly modeled with circuit simulations carried out on a picosecond time
scale to accurately capture the dynamics of Josephson junctions. Simulation of
networks of large numbers of these neurons becomes computationally intensive.
From the perspective of the neural system, the picosecond dynamics of the JJs
are not of primary interest, and one would prefer to treat each synapse,
dendrite, and neuron as an input-output device with a model that accurately
captures the circuit dynamics on the nanosecond to microsecond time scales
while not explicitly treating the picosecond behavior of the underlying
circuit elements.

Here we introduce a phenomenological model of loop neurons and their
constitutive elements that accurately captures the transfer characteristics of
the circuits without solving the underlying circuit equations. Dendrites are
revealed to be central to the system. Each dendrite is treated with a single,
first-order ordinary differential equation (ODE) that describes the output of
the element as a function of its instantaneous inputs and internal state.
These equations take the form of a leaky integrator with a nonlinear driving
term. The input to each dendrite is flux, and the output is an integrated
current, which is coupled through a transformer into another dendrite.
Synapses are dendrites with a closed-form expression for the input flux
following each synapse event. The soma of a neuron is also modeled with the
same equations as a dendrite with two modifications. First, when the output
current reaches a specified level, threshold is reached, and an electrical
signal is sent to a transmitter circuit, which produces a pulse of light.
Second, a refractory dendrite triggered off this output is coupled back to the
soma inhibitively to achieve a refractory period.

By working at the phenomenological level, the time to simulate single
dendrites is decreased by a factor of ten thousand while maintaining an
accuracy of one part in ten thousand. The speed advantage grows with the size
of the system being simulated and the duration of the simulation. The
functional form of the speedup with size and duration has not been fully
investigated, as solving the full system of circuit equations becomes very
time consuming for even small systems. In addition, the model makes
transparent the qualitative behavior of all components of the system,
including the similarities to biological neurons as well as other physical
systems such as Ising models, spin glasses, and pulse-coupled oscillators.

We begin by motivating the form of the model based on circuit considerations.
We then describe the means by which the form of the driving term in the leaky
integrator is obtained. Error is quantified by comparison with full circuit
simulations, and convergence is investigated as a function of time step size.
Numerical examples of dendrites, synapses, and neurons are presented. An
example of a neuron with a dendritic arbor performing an image-classification
task is given to illustrate the utility of the model. Further extensions to
enable theoretical treatment of very large systems are discussed.

##  2 Overview of Loop Neurons

Figure 1: Schematic of a loop neuron with an elaborate dendritic tree. The
complex structure consists of excitatory and inhibitory synapses
($\mathsf{S_{e}}$ and $\mathsf{S_{i}}$) that feed into dendrites
($\mathsf{D}$). Each dendrite performs computations on the inputs and
communicates the result to other dendrites for further processing or on to the
cell body of the neuron ($\mathsf{N}$). The neuron cell body acts as the final
thresholding stage, and when its threshold is reached, light is produced by
the transmitter ($\mathsf{T}$), which is routed to downstream synaptic
connections via photonic waveguides.

Research in superconducting optoelectronic networks of loop neurons aspires to
realize artificial neural systems with scale and complexity comparable to the
human brain. We have introduced the concepts of loop neurons in a number of
papers [6, 7, 8, 9, 10, 11], and we have demonstrated many of the principles
experimentally [12, 16, 13, 14, 15]. A schematic diagram of a loop neuron is
shown in Fig. 1, depicting the complex dendritic tree that appears central to
the computations of loop neurons [9, 17]. In these neurons, integration,
synaptic plasticity, and dendritic processing are implemented with inductively
coupled loops of supercurrent. It is due to the prominent role of
superconducting storage loops that we refer to devices of this type as loop
neurons.

Figure 2: Circuit diagrams. (a) A two-junction SQUID that forms the dendritic
receiving loop. (b) A dendrite with receiving (R) and integrating (I) loop.
The integrating loops of two other dendrites are input to a collection coil
(C) that delivers flux to the receiving loop. (c) A synapse formed with an SPD
input to a dendrite. (d) A soma realized as a dendrite with thresholding
component in the integration loop and initial stage of the transmitter (T).

Operation of loop neurons is as follows. Photons from upstream neurons are
received by a superconducting single-photon detector (SPD) at each synapse.
Using a superconducting circuit comprising two Josephson junctions (JJs)
coupled to the SPD, synaptic detection events are converted into an integrated
supercurrent which is stored in a superconducting loop. The amount of current
added to the integration loop during a photon detection event is determined by
the synaptic weight. The synaptic weight is dynamically adjusted by another
circuit combining SPDs and JJs, and all involved circuits are analog. When the
integrated current of a given neuron reaches a (dynamically variable)
threshold, an amplification cascade begins, resulting in the production of
light from a waveguide-integrated semiconductor light emitter. The photons
thus produced fan out through a network of dielectric waveguides and arrive at
the synaptic terminals of other neurons where the process repeats.

The core active component of loop neurons is a circuit known as a
superconducting quantum interference device (SQUID), which comprises two JJs
in parallel. A circuit diagram is shown in Fig. 2(a). The SQUID is a three-
terminal device with a bias, ground, and an input which couples flux into the
loop formed by the two JJs and inductors. For the present purpose, the flux
input is the active signal. When the SQUID is current-biased below the
critical current of the two JJs and the flux input is below a bias-dependent
threshold, the SQUID remains superconducting, and the voltage across the
device is zero. When the applied flux exceeds the threshold at a given bias
point, the JJs will begin to emit a series of voltage pulses known as fluxons,
and the time-averaged voltage across the device will become non-zero.

To form a dendrite from a SQUID, we first ensure that it is biased below the
critical current of the JJs so it is quiescent when no flux is applied. The
output of the SQUID is captured by an $L$-$R$ loop that performs current
integration with a leak [Fig. 2(b)]. When sufficient flux is input to the
SQUID to drive the JJs to begin producing voltage pulses, the pulses drive
current into the $L$-$R$ loop, these pulses are summed in the inductor, and
the accumulated signal leaks with the $L/R$ time constant of the loop. This
configuration of a SQUID coupled to an $L$-$R$ loop is referred to as a
dendrite, the SQUID that receives input flux is referred to as the receiving
(R) loop, and the $L$-$R$ component is referred to as the integration (I)
loop. The current stored in the integration loop produces the signal that will
be communicated to other dendrites. To receive signal from multiple input
dendrites, a passive collection coil (C) is used. A circuit diagram of a
dendrite with two inputs, a collection coil, a receiving loop, and an
integrating loop coupled to an output is shown in Fig. 2(b). All coupling
between dendrites is through magnetic flux communicated through mutual
inductors. The use of transformers for this purpose mitigates cross talk and
enables high fan-in [17].

To form a synapse from a dendrite, we attach an SPD to the flux input to the
receiving loop. This circuit is shown in Fig. 2(c). When an SPD detects a
photon, it rapidly switches from zero resistance to a large resistance,
diverting the bias current to the other branch of the circuit. This current is
coupled into the receiving loop of the dendrite as flux, driving the dendrite
above threshold and adding current to the dendrite’s integration loop.

To form a neuron from a dendrite, two modifications are required. First, the
integration loop must be equipped with a thresholding element that drives a
transmitter circuit to produce light when the integrated signal reaches this
threshold. This thresholding element coupled to the transmitter is shown in
Fig. 2(d); it is referred to as a tron and is described in more detail in Sec.
5. Second, an additional dendrite is attached to the neuron that provides
negative feedback to achieve a refractory period (not shown in the circuit
diagram for simplicity). The refractory period is a brief period of quiescence
following a neuronal spike event. When the neuron’s integration loop reaches
threshold, this refractory dendrite is driven, accumulates signal in its
integration loop, and this signal suppresses the state of the neuron’s
receiving loop. The time constant of the refractory dendrite establishes the
refractory period of the neuron.

To summarize, a dendrite is a SQUID with a flux input adding current to a
leaky integration loop. A synapse is a single-photon detector coupled to a
dendrite. The soma of a neuron is a dendrite with a thresholding element in
the integration loop as well as a second dendrite that provides feedback for
refraction. To construct full loop neurons, many synapses are coupled into a
dendritic arbor which feeds forward into the soma. The output of the soma is
an optical pulse that couples light to a network of waveguides and delivers
faint photonic signals to downstream synapses where they are received with
single-photon detectors. This construction is illustrated schematically in
Fig. 1, and the basic circuits are shown in Fig. 2. The current in a dendritic
integration loop is analogous to the membrane potential of a biological
dendrite [18, 19], and these signals are the principal dynamical variables of
the system. Inhibitory synapses can be achieved through mutual inductors with
the opposite sign of coupling. Complex arbors with multiple levels of
dendritic hierarchy can be implemented to perform various computations [9] as
well as to facilitate a high degree of fan-in [17].

We see that dendrites are central to loop neurons. The states of current in
all dendritic integration loops specifies the state of the system. A
phenomenological model of loop neurons must therefore capture the temporal
evolution of these currents as well as the coupling between dendrites. With
these concepts in mind we proceed to construct the model.

##  3 Dendrite Model

Figure 3: The response of a SQUID. (a) Time trace of a SQUID biased in the
voltage state. Fluxon peaks are marked with crosses of different colors for
the two JJs. The time average is also shown as calculated by taking the
average of the time trace between two fluxons produced by the same junction.
(b) The time-averaged voltage across the SQUID as a function of applied flux,
$\phi$, normalized to the magnetic flux quantum, $\Phi_{0}$, for several
values of normalized bias current, $i_{\mathrm{b}}$.

As mentioned in Sec. 2, a SQUID is the primary active element of a dendrite.
To motivate the phenomenological dendrite model we require a quantitative
understanding of SQUID operation. The two-junction SQUID of Fig. 2(a) with
symmetrical inductances was modeled using a first-principles circuit model
[20], and the results are shown in Fig. 3. Figure 3(a) shows a time trace of
the voltage across the SQUID when it is in the voltage state. The peaks
corresponding to fluxon production are evident, and the time-averaged voltage
is also shown. While the voltage is a rapidly varying function of time on the
picosecond scale, the time-averaged voltage is steady.

This behavior is analyzed systematically in Fig. 3(b). The time-averaged
voltage of a symmetric SQUID is plotted as a function of the applied flux for
several values of the bias current, which has been normalized to the critical
current of a single junction ($i_{\mathrm{b}}=I_{\mathrm{b}}/I_{c}$). The
time-averaged voltage, plotted on the left $y$-axis, is proportional to the
rate of flux-quantum production, plotted on the right $y$-axis. The
relationship results from the single-valued nature of the superconducting wave
function around the closed SQUID loop, which requires that

$\int_{0}^{t_{\mathrm{fq}}}V_{\mathrm{sq}}(t)\,dt=\Phi_{0}.$ (1)

Equation 1 informs us that the time required to produce a single fluxon,
$t_{\mathrm{fq}}$, is related to the voltage across the squid,
$V_{\mathrm{sq}}$. For constant voltage, we obtain
$r_{\mathrm{fq}}=1/t_{\mathrm{fq}}=V_{\mathrm{sq}}/\Phi_{0}$. When viewed over
time scales appreciably longer than $t_{\mathrm{fq}}$, it makes sense to speak
of a rate of flux-quantum production, $r_{\mathrm{fq}}$, and this is the first
element of our model: when driven to the active state, a dendrite will begin
to produce fluxons, which carry current, and we can track this current by
monitoring the rate of fluxon production while ignoring the picosecond
dynamics by which the JJs produce the fluxons.

Several features of Fig. 3(b) are pertinent to the present study. First, for a
given value of $i_{\mathrm{b}}$, a finite value of flux is required before the
SQUID enters the voltage state. This provides a non-zero threshold that is
relevant to dendritic computation. This threshold can be adjusted with the
bias current. Second, the response is periodic in applied flux, with the
period being $\Phi_{0}/2$, where $\Phi_{0}=h/2e\approx 2\times
10^{-15}\mathrm{V}\cdot\mathrm{s}=2\,\mathrm{mV}\cdot\,\mathrm{ps}$ is the
magnetic flux quantum. To maintain a monotonic response, the applied flux must
be limited to this value [17]. Third, if the inductors $L_{1}$ and $L_{2}$ are
equal, the response of the SQUID is symmetric about $\Phi=0$. These features
will be discussed further as the study proceeds.

The second element of the model captures the integration and leak of the
current generated when the SQUID is driven above threshold. These behaviors
are accomplished by adding an $L$-$R$ loop to the output of the SQUID, as
shown by the integration loop labeled I in Fig. 2(b). The current integrated
in this branch of the circuit is the quantity of interest for the dendrite. It
is this quantitiy that will couple to other dendrites or the neuron cell body,
and it is this quantity we wish to track with our phenomenological model.
Dendrites comprising a receiving loop coupled to an integrating loop [Fig.
2(b)] are referred to as RI dendrites. We know from elementary circuit theory
that the $L$-$R$ loop will result in exponential decay of signal with a time
constant of the dendritic integration loop given by
$\tau_{\mathrm{di}}=L_{\mathrm{di}}/R_{\mathrm{di}}$.

We can now write down a postulated expression for the signal $s$ stored in the
integration loop of a dendrite:

$\beta\frac{d\,s}{d\tau}=r\left(\phi,s;i_{\mathrm{b}}\right)-\alpha\,s.$ (2)

Equation 2 states that the signal $s$ grows in time due to the driving term,
which is the rate of flux quantum production, denoted by $r$. This function
$r$ depends on the applied flux to the receiving loop of the dendrite,
$\phi=\Phi/\Phi_{0}$, as well as of the signal present in the integration
loop, $s$. The rate also depends on the bias current, $i_{\mathrm{b}}$, as a
parameter that throughout this work is assumed to be held constant over times
much longer than the inter-fluxon interval. In Eq. 2 we have formulated the
model in dimensionless units, where $s=I_{\mathrm{di}}/I_{c}$, and
$I_{\mathrm{di}}$ is the current present in the dendritic integration loop.
$\beta=2\,\pi\,I_{c}\,L_{\mathrm{di}}/\Phi_{0}$ is a dimensionless parameter
that quantifies the inductance of the loop, and
$\alpha=R_{\mathrm{di}}/r_{\mathrm{j}}$, where $r_{\mathrm{j}}$ is the shunt
resistance of each JJ in the resistively and capacitively shunted junction
model [21, 22, 23, 20]. The signal $s$ decays at a rate related to $\alpha$
and $\beta$. Specifically, the time constant for decay is given by
$\tau_{\mathrm{di}}=L_{\mathrm{di}}/R_{\mathrm{di}}=\beta/\omega_{c}\alpha$,
where $\omega_{c}$ is the Josephson characteristic frequency discussed in
Appendix A, and we include the subscript on $\tau_{\mathrm{di}}$ to refer to
the dendritic integration loop and distinguish that quantity from the
dimensionless time variable entering Eq. 2. Equation 2 is a leaky integrator
ODE. The drive term is the rate of flux quantum production, and the leak term
gives simple exponential decay, as expected from an $L$-$R$ circuit.

Figure 4: Rate arrays for the dendrite with a receiving and integrating loop
for three values of $i_{\mathrm{b}}$. (a) The normalized bias current
$i_{\mathrm{b}}=1.35$. (b) $i_{\mathrm{b}}=1.65$. (c) $i_{\mathrm{b}}=1.95$.

Dendrites are coupled to each other through flux. The coupling flux from
dendrites indexed by $j$ to dendrite $i$ is given by

$\phi_{i}=\sum_{j=1}^{n}\,J_{ij}\,s_{j},$ (3)

where $J_{ij}$ is a coupling term proportional to the mutual inductance that
includes contributions from all the transformers present on the collection
coil in Fig. 2(b). Equation 3 shows that coupling between dendrites is due to
the signal in the integration loop of one dendrite being communicated as flux
into the receiving loop of a subsequent dendrite. The signal in the subsequent
dendrite is then obtained through the evolution of Eq. 2 with the flux from
the first dendrite providing the flux $\phi$ in the function $r$ and the
signal from the second dendrite providing the $s$ term. The $i_{\mathrm{b}}$
term entering $r$ refers to the bias on the second dendrite and is treated
here as a parameter rather than a dynamical variable. Further details
regarding the derivation of these expressions is given in Appendix C.

Equations 2 and 3 constitute the phenomenological model of a dendrite. A
neuron or network can be simulated by solving these coupled equations for all
dendrites in the system. However, we have not specified the form for the rate
function, $r\left(\phi,s;i_{\mathrm{b}}\right)$, which is central to the
model.

We have arrived at Eq. 2 as a postulate; this expression is not directly
obtained from the underlying circuit equations. The postulate is that there
will be a function $r\left(\phi,s;i_{b}\right)$ such that Eq. 2 provides an
accurate account of the signal present in a dendrite’s integration loop under
the circumstances of interest for loop neurons, provided we only inquire about
the signal over time scales appreciably longer than the inter-fluxon interval,
which is on the order of 10 ps. We aim to interrogate dendrites on time scales
of 100 ps or longer, with neuron and network activity of interest on time
scales from nanoseconds to the longest timescales that can be simulated under
the limits of computational resources.

Our procedure for obtaining $r\left(\phi,s;i_{\mathrm{b}}\right)$ is as
follows. A dendrite with a SQUID as a receiving loop and an $L$-$R$ branch as
an integrating loop is numerically modeled with the circuit equations given in
Appendix A. A constant value of flux is applied to the receiving loop, and the
rate of flux-quantum production is monitored as a function of time while
current accumulates in the integration loop. For these simulations, the
resistance of the integration loop is set to zero. This procedure is repeated
for many values of $\phi$ and $i_{\mathrm{b}}$ to obtain what we refer to as
“rate arrays”, which are shown in Fig. 4, where
$r\left(\phi,s;i_{\mathrm{b}}\right)$ is plotted as a function of $\phi$ and
$s$ for three values of $i_{\mathrm{b}}$. Here we work in dimensionless units,
so the units of $r$ are fluxons generated per unit of dimensionless time,
$\tau$, which is related to the JJ characteristic frequency (Appendix A). It
can be seen that the value of $r\left(\phi,s;i_{\mathrm{b}}\right)$ is
monotonically increasing with $\phi$ over the range considered here, while
accumulation of $s$ decreases the rate of flux-quantum production. This
decrease is because addition of current to the integration loop diverts the
bias away from the SQUID, so the voltage is decreased, and the rate is reduced
in accordance with Eq. 1. When sufficient signal is accumulated in the
integration loop, the rate of flux quantum production drops to zero, and we
say the loop is saturated.

For a small value of $i_{\mathrm{b}}$ [Fig. 4(a)] a large amount of flux is
required to drive the dendrite above threshold to the active state, and a
small signal $s$ is present at saturation. As $i_{\mathrm{b}}$ is increased
[Figs. 4(b) and (c)] the threshold is reduced, and the saturation level is
increased. The qualitative shape of $r$ for different values of
$i_{\mathrm{b}}$ is similar, and the surfaces for smaller $i_{\mathrm{b}}$ are
seen to fit inside those for larger $i_{\mathrm{b}}$. For this reason we refer
to these surfaces as $r$-shells, with each shell corresponding to a fixed
value of $i_{\mathrm{b}}$. Network adaptation, training, and learning consist
in finding the values of $i_{\mathrm{b}}$ that achieve the desired input-
output relationships for a given computational or control task.

While Fig. 4 shows $r$ for $0\leq\phi\leq 1/2$, these surfaces are symmetric
about zero for the circuit of Fig. 2(b) with symmetric inductances in the
SQUID [as is also the case in the SQUID response of Fig. 3(b)], and they are
periodic in normalized flux $\phi$ with period of unity (period of $\Phi_{0}$
in SI units). One consequence of this response is that it may be necessary to
restrict the applied flux to $\phi\leq 1/2$ to retain a monotonic response
[17], and if inhibition is applied to a dendrite ($\phi<0$) it may be
necessary to limit this flux below the threshold for activity so that
inhibition does not drive the dendrite to the active state. These concepts are
discussed further in Sec. 6.

Figure 5: Response parameters of the RI dendrite extracted from the rate
arrays. (a) The thresholds to positive and negative applied flux as a function
of the dendritic bias current. (b) The value of saturation signal as a
function of the normalized flux applied to the receiving loop for several
values of the normalized bias current, $i_{\mathrm{b}}$.

From the rate arrays we obtain several quantities of interest that facilitate
dendrite design. In particular, the value of $\phi$ for which $r$ becomes non-
zero as a function of $i_{\mathrm{b}}$ determines the dendrite threshold. This
function is shown in Fig. 5(a) for both positive and negative flux. These
traces can be leveraged to anticipate when a dendrite will become active as
well as the maximum value of inhibitory flux that can be applied before
activity results, which is necessary in design of the refractory dendrite that
quenches activity of a soma. The saturation values of $s$ as a function of
applied flux and bias current are also important for anticipating the maximum
value of signal that a dendrite will accumulate, and therefore the maximum
value of flux it will couple into another dendrite. We refer to this value as
$s_{\mathrm{max}}$, and Fig. 5(b) shows its behavior as a function of both
positive and negative applied flux ($s_{\mathrm{max}}^{+/-}$). The inset shows
the value versus $i_{\mathrm{b}}$ for the maximum absolute value of applied
flux, $|\phi|=1/2$.

Obtaining the rate arrays is somewhat computationally intensive, requiring a
large number of numerical simulations of the dendrite circuit. However, these
computations only need to be carried out once. The rate arrays are then
accessed by the phenomenological model as a look-up table, and the same arrays
are used for all values of the integration loop inductance and leak rate.
Using the same arrays for all values of inductance is possible because in all
cases the inductance of the integration loop is much larger than that of the
receiving loop, so the fraction of current bias that initially goes into the
integration branch of the circuit is small in all cases of interest here. If
properties of the receiving loop—the SQUID—are changed, new rate arrays will
need to be calculated. Examples include changing the value of
$\beta_{\mathrm{c}}$ of the JJs or changing the asymmetry of the inductances.
However, there appears to be little advantage from employing SQUIDs with a
variety of $\beta_{\mathrm{c}}$ values. Leveraging different asymmetry designs
may be fruitful, but a small number of variants will likely suffice. We
anticipate a symmetric SQUID and a highly asymmetric SQUID will be the two
cases of primary interest.

Figure 6: Comparing the phenomenological model to circuit equations in the
case of a linearly ramping drive. The time step for the phenomenological model
was 100 ps. (a) The flux drive signal. (b) The smallest integration loop with
$\beta/2\pi=10^{2}$. (c) Integration loop with $\beta/2\pi=10^{3}$. (d)
$\beta/2\pi=10^{4}$. In (a)-(d) the black traces show the values for the
circuit simulations while the dashed colored traces show the values for the
phenomenological model. The error of the drive signal was
$\chi^{2}_{\mathrm{drive}}<5\times 10^{-11}$ for all cases shown here.

Equipped with the phenomenological model (Eqs. 2 and 3) as well as the rate
arrays, we can compare the model performance to explicit solution of the
circuit equations. Throughout this work, we numerically solve systems
described by Eqs. 2 and 3 using a forward Euler method to facilitate
computational speed. In discrete form, Eq. 2 becomes

$s_{\tau+1}=s_{\tau}\left(1-\Delta\tau\frac{\alpha}{\beta}\right)+\frac{\Delta\tau}{\beta}r(\phi,s;i_{\mathrm{b}}).$
(4)

We solve the circuit equations given in Appendix A using a Runge-Kutta method
implemented with the Python SciPy function solve_ivp for initial-value
problems, which uses an adaptive time mesh.

As a first test, Fig. 6 compares the phenomenological model to the circuit
equations in the presence of a linear ramp function as the applied flux. The
form of the applied flux into the receiving loop ($\phi$) is shown in Fig.
6(a). The rate arrays are obtained only for discrete values of applied flux
(200 values between 0 and $1/2$), so this applied drive is relevant to
quantify performance when the applied flux takes intermediate values between
those explicitly present in the rate arrays. No interpolation was used. At
each time step, the three inputs to the rate array ($\phi$, $s$, and
$i_{\mathrm{b}}$) are rounded to the closest value for which the rate array
has been evaluated, and that value of $r$ is used in Eq. 4 at that time step.
The dendrite signals are shown in Fig. 6(b)-(d) for three values of the
inductance parameter, $\beta$, from $\beta/2\pi=10^{2}-10^{4}$. For each value
of $\beta$, four time constants are considered from 10 ns to 1.25 µs. The
value of $\beta/2\pi$ quantifies the number of fluxons that can be
accommodated in the storage loop, and the time constant is converted to the
dimensionless parameter $\alpha$ for simulation in the dimensionless model.
All data from the phenomenological model in Fig. 6 were obtained using a time
step of 100 ps (converted to dimensionless time for Eq. 4). The circuit model
time traces are shown as black lines, while the traces obtained with the
phenomenological model are shown as various dashed colored lines.

Figure 7: Comparing the phenomenological model to circuit equations in the
case of a random sequence of square pulses. The time step for the
phenomenological model was 100 ps. (a) The flux drive signal. (b) The smallest
integration loop with $\beta/2\pi=10^{2}$. (c) Integration loop with
$\beta/2\pi=10^{3}$. (d) $\beta/2\pi=10^{4}$. The error of the drive signal
was $\chi^{2}_{\mathrm{drive}}=5.03\times 10^{-5}$ for all cases shown here.

Accuracy of the phenomenological model is quantified with a $\chi^{2}$ of the
form

$\chi^{2}=\frac{\sum_{i=1}^{n_{t}^{\mathrm{soen}}-1}\left|s_{\mathrm{soen}}\left(t_{i}^{\mathrm{soen}}\right)-s_{\mathrm{ode}}^{\mathrm{interp}}\left(t_{i}^{\mathrm{soen}}\right)\right|^{2}\Delta
t_{i}^{\mathrm{soen}}}{\sum_{i=1}^{n_{t}^{\mathrm{ode}}-1}\left|s_{\mathrm{ode}}\left(t_{i}^{\mathrm{ode}}\right)\right|^{2}\Delta
t_{i}^{\mathrm{ode}}}.$ (5)

The subscript or superscript “soen” refers to the phenomenological model,
while “ode” refers to the first-principles circuit model of Appendix A. In Eq.
5, $n_{t}^{\mathrm{soen}}$ is the number of time steps in a given SOEN
simulation, $s_{\mathrm{soen}}$ is the state quantity obtained with the SOEN
phenomenological model, $t_{i}^{\mathrm{soen}}$ is the time of the given
simulation at time step $i$, $s_{\mathrm{ode}}^{\mathrm{interp}}$ is the state
quantity obtained with the ODE circuit model solved with an adaptive time mesh
and interpolated to the coarser phenomenological model time mesh, and $\Delta
t_{i}^{\mathrm{soen}}=t_{i+1}^{\mathrm{soen}}-t_{i}^{\mathrm{soen}}$. In the
denominator, $n_{t}^{\mathrm{ode}}$ is the number of time steps in the ODE
solution with the adaptive mesh, $s_{\mathrm{ode}}(t_{i}^{\mathrm{ode}})$ is
the ODE solution on the non-uniform time grid without interpolation, and
$\Delta t_{i}^{\mathrm{ode}}=t_{i+1}^{\mathrm{ode}}-t_{i}^{\mathrm{ode}}$. In
all cases shown in Fig. 6, the value of $\chi^{2}$ is less than $7\times
10^{-5}$, and the time required to run the ODE model exceeded that to run the
SOEN model by a factor of more than one thousand. The accuracy is worst for
the smallest storage loop ($\beta/2\pi=10^{2}$) and the shortest time constant
($\tau_{\mathrm{di}}=10$ ns). All values of $\chi^{2}$ and ratios of
simulation times are given in Table 1 in Appendix D.

As another test case, a series of square flux pulses was input to a dendrite
for the same values of $\beta$ and $\tau_{\mathrm{di}}$, as shown in Fig. 7.
The heights of these pulses were drawn randomly from the interval between the
flux threshold [Fig. 5(a)] and the maximum value of 1/2. The pulse durations
were drawn randomly from the interval between 5 ns and 40 ns, while the pauses
between pulses were drawn randomly from the interval between 10 ns and 100 ns.
In Fig. 7, 10 such pulses were applied. The largest value of $\chi^{2}$ was
$3.32\times 10^{-4}$, again in the case of smallest $\beta$ and
$\tau_{\mathrm{di}}$, while the time to complete the ODE simulations was
longer than that for the SOEN simulations by at least $10^{3}$ in all cases.
Complete data for the square pulses is given in Table 2 in Appendix D.

Figure 8: Similar comparison as Fig. 7 except the time step was 1 ns. The
error in reproducing the drive signal was
$\chi^{2}_{\mathrm{drive}}=7.35\times 10^{-3}$. Figure 9: The same simulations
as Fig, 8 except zoomed on the time window around the first input pulse.

To push the limits of speed, the same linear ramp and square pulse cases were
carried out using a time step of 1 ns in the phenomenological model. The
results for the square-pulse drive are shown in Fig. 8 for the same values of
$\beta$ and $\tau_{\mathrm{di}}$. In all cases, the agreement between the two
models is good, except for the smallest inductance of $\beta/2\pi=10^{2}$.
Poor performance in these cases results from the fact that the term
$\Delta\tau\frac{\alpha}{\beta}$ in Eq. 4 becomes large, leading to numerical
instability. More detail is shown in Fig. 9 where a temporal zoom on just the
first input square pulse is presented. The numerical performance is seen to be
unacceptable for the smallest value of $\beta$ and only marginally acceptable
for values above that.

Figure 10: Quantification of the phenomenological model versus the circuit
model. For the dendrite considered here, $\tau_{\mathrm{di}}=250$ ns and
$i_{\mathrm{b}}=1.702$. These simulations included 160 random square pulses,
and the circuit model was set to converge to a relative and absolute accuracy
of $10^{-5}$ in the solve_ivp function. (a) The values of $\chi^{2}$ as a
function of time step $dt$ for three values of $\beta$. (b) The ratio of the
time required to run the simulations for the two numerical approaches as a
function of time step $dt$.

To determine the largest acceptable time step, square-pulse-input simulations
were conducted for numbers of input pulses ranging from 10 to 160 with time
steps ranging from 10 ps to 10 ns for the same values of $\beta$ as above. To
keep the computation time manageable, only $\tau_{\mathrm{di}}=250$ ns was
considered, but individual test cases were conducted with other values to
ensure the conclusions were insensitive to this number. The results of the
$\chi^{2}$ values for each simulation as well as the ratio of computational
run times were stored, and the results for the case of 160 random square pulse
inputs is shown in Fig. 10. Further data on the other cases is given in Fig.
27 in Appendix D. The $\chi^{2}$ values again show that the lowest value of
$\beta$ is the most difficult to match, and the $\chi^{2}$ converges close to
$10^{-4}$ with a time step of 200 ps. All other $\beta$ values approach
$\chi^{2}=10^{-5}$ with this same time step. The simulated time interval in
this case was 12 µs, and with a time step of 200 ps it took over ten thousand
times longer to solve the system of circuit ODEs as to step through the
phenomenological model. Based on this analysis, a time step of 200 ps appears
an optimal compromise between speed and accuracy. This time step is likely the
convergence point in these simulations as 200 ps is the rise and fall time of
the square pulses. To simulate 12 µs of activity of a single dendrite with the
phenomenological model required 1.02 s, while simulating the same time
interval with the ODE model required 2.78 hours. Attempting to simulate longer
time intervals with the circuit model becomes impracticable.

Figure 11: The RTTI dendrite. (a) Circuit diagram showing the addition of two
transmission loops. (b) An example $r$-shell for the RTTI dendrite.

Before proceeding to treat synaptic inputs to dendrites, we note an extension
of the model to dendrites with additional circuit complexity. The response
functions of the two-JJ dendrite is useful for many operations, but in certain
cases it is desirable for the value of the saturation current to be less
sensitive to the applied flux. This can be accomplished by adding a Josephson
transmission line between the receiving and integrating loops of the dendrite,
as shown in Fig. 11(a). Due to the presence of two transfer loops between the
receiving and integrating loops, we refer to this as the RTTI dendrite
(receive-transfer-transfer-integrate). An example rate array for the RTTI
dendrite is shown in Fig. 11(b). While it is similar to the RI dendrite in
that it has a threshold of applied flux and decreases with integrated $s$, the
overall shape is more complicated, with appreciable structure close to
saturation. The details of this structure can be used for various
computations. In particular, the fact that the response is relatively flat for
certain ranges of $\phi_{\mathrm{r}}$ can be used to achieve a digital
response in which the output is insensitive to the exact input value of
$\phi_{\mathrm{r}}$. This attribute will be exploited below to achieve the OR
gate.

Figure 12: The values of flux thresholds (a) and saturation levels (b) for the
RTTI dendrite extracted from the rate arrays, analogous to these quantities
shown in Fig. 5 for the RI dendrite.

As for the case of the RI dendrite, we extract the thresholds and saturation
values from the rate arrays. These functions are shown in Fig. 12. In Fig.
12(a) it can be observed that the thresholds for positive and negative applied
flux are now slightly asymmetric due to the breaking of the symmetry of the
circuit with the low-inductance transfer loop to the right of the receiving
loop. More importantly, the values of $s_{\mathrm{max}}^{+/-}$ shown in Fig.
12(b) are seen to have several regions of nearly flat response where the
saturation level of the dendritic integration loop does not depend on applied
flux. Such a response is useful for obtaining digital behaviors, as we show
below when discussing logic gates in Sec. 4. The peak value of
$s_{\mathrm{max}}^{+\\-}$ as a function of $i_{\mathrm{b}}$ is shown in the
inset to Fig. 12(a).

To simulate the RTTI dendrite with a conventional model, nine coupled ODEs are
used, making the simulations even slower than the RI dendrite with its five
coupled equations. With the phenomenological model, the RTTI dendrite reduces
to a single ODE, just as the RI dendrite. The only difference is the form of
the $r$-shells that provide the driving term.

##  4 Synapses

As described in Secs. 2 and 3, the dendritic receiving loop is a SQUID. In
many applications in science and technology, a SQUID is used as a measurement
device with unmatched sensitivity for detecting magnetic flux. Based on the
response curves of Fig. 3(b), it can be seen that a SQUID can be used as a
flux-to-voltage transducer. Such a device can also be used to measure very low
current levels when the current is coupled into the SQUID as flux using a
transformer. It is in this mode of operation that we convert a dendrite into a
synapse.

When an SPD detects a photon, a current pulse is diverted out of the circuit,
and that current returns with an $L/R$ time of around 35 ns, which sets the
recovery time and maximum rate of synapse events. To couple an SPD into a
dendrite to form a synapse, this current output can be coupled into the SQUID
that forms the receiving loop through a transformer, as shown in Fig. 2(c).
The equations of motion for the SPD can be solved exactly under a simple model
[24] to obtain the currents as a function of time. In this model, the
resistance of the SPD is zero until a photon is detected, at which point it
switches to a finite resistance for a finite duration. We are interested in
the flux applied to a dendrite by the SPD receiver circuit; this flux is the
product of the current $I_{2}$ with the mutual inductance of the transformer,
which is a quantity we can select in design. Therefore, obtaining the currents
is sufficient to obtain the flux, which we write as

$\begin{split}\phi_{\mathrm{r}}&=\phi_{\mathrm{peak}}\,\left(1-\frac{\tau_{\mathrm{rise}}}{\tau_{\mathrm{fall}}}\right)\\\
&\times\begin{cases}\left[1-e^{-t/\tau_{\mathrm{rise}}}\right],\text{
for}\,0\leq t\leq t_{0}\\\
\left[1-e^{-t_{0}/\tau_{\mathrm{rise}}}\right]\,e^{-(t-t_{0})/\tau_{\mathrm{fall}}},\text{
for}\,t>t_{0}.\end{cases}\end{split}$ (6)

In Eq. 6 it is assumed a synapse event occurred at time $t=0$, which is the
time at which the SPD entered the resistive state, and $t_{0}$ is the duration
for which the SPD stays in the resistive state following the absorption of a
photon, which we take to be 200 ps, following the model of Ref. 24. This
duration is commensurate with the 200 ps time step that we found to be optimal
in Sec. 3. The resistance $r_{1}$ [Fig. 2(c)] is the resistance of the SPD in
the time window $t_{0}$ after absorbing a photon, while $r_{2}$ is a fixed
resistance chosen to obtain a sufficiently long $L/R$ recovery time to prevent
the SPD from latching [25]. The total inductance of the circuit is
$L_{\mathrm{tot}}=L_{1}+L_{2}+L_{3}$, which includes the kinetic inductance of
the SPD itself ($L_{1}$), the inductance of the transformer input coil
($L_{3}$), and any additional parasitic inductance in the circuit $L_{2}$. In
terms of these parameters,
$\tau_{\mathrm{rise}}=L_{\mathrm{tot}}/(r_{1}+r_{2})$ and
$\tau_{\mathrm{fall}}=L_{\mathrm{tot}}/r_{2}$. Typically
$\tau_{\mathrm{rise}}$ is a few tens of picoseconds and $\tau_{\mathrm{fall}}$
is 30 ns - 50 ns [15]. The quantity $\phi_{\mathrm{peak}}$ can be chosen in
design and set by the transformer between the SPD and the dendrite. Throughout
this work we take $\phi_{\mathrm{peak}}=0.5$ in units of $\Phi_{0}$.

Figure 13: Output of the phenomenological model for the case of a synapse. (a)
Flux output from the SPD. (b) Signal $s$ in the dendritic integration loop for
the case of $\beta/2\pi=10^{2}$. The solid lines result from a time step of
100 ps, while the dashed lines used 1 ns. (c)-(e) $\beta/2\pi=10^{3}$,
$10^{4}$, and $10^{5}$, respectively.

We now consider several examples of the phenomenological model with the output
of an SPD coupled as input to a dendrite. In Fig. 13 we show the response of
the dendrite to a single synaptic pulse. The time course of the total applied
flux to the receiving loop, $\phi$, is shown in Fig. 13(a), and four values of
$\beta$ are shown in Fig. 13(b)-(e), each for four values of
$\tau_{\mathrm{di}}$. In each case, the dashed line corresponds to using a
time step of 1 ns in the phenomenological model, while the solid line
corresponds to 100 ps. Again we see that the worst performance occurs with
small $\beta$ and small $\tau_{\mathrm{di}}$.

Figure 14: Similar data to Fig. 13 except signals are in response to a train
of five input synapse events.

Similar data is shown in Fig. 14 for an arbitrary sequence of five synapse
events. Dendritic integration is evident, and the effects of $\beta$ and
$\tau_{\mathrm{di}}$ can be seen: for smaller $\beta$ and larger
$\tau_{\mathrm{di}}$, the integration loop reaches saturation after a small
number of input pulses, while for larger $\beta$ and/or smaller
$\tau_{\mathrm{di}}$, the integrated signal does not reach saturation even by
the end of the train of five pulses.

Using synapses and dendrites to perform various temporal extensions of basic
logic operations is a primary function of biological neural computation. Here
we consider four basic logic operations: AND, OR, AND-NOT, and XOR. Each of
these operations is accomplished with a single dendrite receiving input from
two synapses.

Figure 15: Operation of the AND gate. (a) Input flux from the two synapses
separately. (b) Combined input flux. (c) Dendrite output signal. (d) Analysis
of the gate response as a function of time delay between the two synapse
events for several values of $i_{\mathrm{b}}$.

The AND operation is shown in Fig. 15. The flux from the two synapses
independently is shown in Fig. 15(a), their sum is shown in Fig. 15(b), and
the signal in the dendritic integration loop is shown in Fig. 15(c). In Figs.
15-18 the flux from a single SPD to the dendrite is referred to as
$\phi_{\mathrm{spd}}$, while the combined flux from both SPDs to the receiving
loop of the dendrite is labeled $\phi_{\mathrm{r}}$. To achieve the AND
response, the dendrite is configured such that the flux from a single synapse
is insufficient to drive the dendrite above threshold, but the flux from two
synapses together is adequate, and signal is added to the integration loop
only when the pulses are coincident.

Because the synaptic flux and the dendritic signal are both functions of time,
it is important to know how the accumulated signal will vary with the delay
between the arrival of the two synapse events. This data is shown in Fig.
15(d) for multiple values of the bias to the dendrite, $i_{\mathrm{b}}$. Such
a circuit can be used as a temporal coincidence detector, as described in Ref.
9.

Figure 16: Operation of the OR gate. (a) Input flux from the two synapses
separately. (b) Combined input flux. (c) Dendrite output signal. (d) Analysis
of the gate response as a function of time delay between the two synapse
events for several values of $i_{\mathrm{b}}$.

The OR operation is shown in Fig. 16. The challenge with OR is that we would
like the amplitude of the signal in the dendritic integration loop to be
identical whether one or both synapses are active. Because both the rate of
fluxon production and the value of saturation depend on the applied flux for
the RI dendrite, this is difficult to accomplish. However, the flat response
of the RTTI dendrite can generate this input-output relationship, and that
type of dendrite has been used to generate Fig. 16. In Fig, 16(b) we can see
that when both synapses are active the applied flux to the receiving loop has
twice the amplitude as when only a single synapse is active, yet the
integrated signal in Fig. 16(c) is nearly identical whether one or both
synapses are active. The response is possible due to the flat-top behavior of
the RTTI dendrite. We have biased the dendrite at a point such that when the
flux from one synapse is present the dendrite just enters the plateau of the
response, and when the flux from both is present it remains on the plateau.

The response as a function of the delay between the two synapse events is
shown in Fig. 16(d). For lower bias points the flat-top response only occurs
for a relatively short coincidence window, providing a response more like the
AND gate. Yet for higher bias points the response is very flat, independent of
temporal delay between the synapse events, just as we desire for the OR gate.
We will see in Sec. 5 that this response is useful in achieving the desired
classification task.

Figure 17: Operation of the AND-NOT gate. (a) Input flux from the two synapses
separately. (b) Combined input flux. (c) Dendrite output signal. (d) Analysis
of the gate response as a function of time delay between the two synapse
events for several values of $i_{\mathrm{b}}$.

The AND-NOT operation is shown in Fig. 17. This logic gate produces a logical
one at the output if and only if the first input is a logical one and the
second input is a logical zero. This is distinct from the NAND gate. For AND-
NOT we return to the RI dendrite. In this case, synapse one is excitatory,
while synapse two is inhibitory with the opposite sign of flux applied as well
as a reduced amplitude. Fig. 17(c) shows that the dendrite only becomes active
when synapse one receives a pulse in the absence of synapse two receiving a
pulse. Extension to various temporal delays in Fig. 17(d) shows the effect of
the asymmetric coupling of the two synapses.

Figure 18: Operation of the XOR gate. (a) Input flux from the two synapses
separately. (b) Combined input flux. (c) Dendrite output signal. (d) Analysis
of the gate response as a function of time delay between the two synapse
events for several values of $i_{\mathrm{b}}$. Here the time axis is narrower
than Figs. 15, 16, and 17. The inset shows the full time axis of the other
logic gate figures.

The XOR operation is shown in Fig. 18. XOR is accomplished in much the same
manner as AND-NOT, except the coupling between the two synapses is equal.
Thus, when either synapse is active in the absence of the other the dendrite
becomes active. In this instance we are making use of the fact that the RI
dendrite response is symmetric with respect to positive or negative flux. The
response as a function of relative delay between the synapse events is shown
in Fig. 18(d), and it is evident that cancellation of the two pulses requires
relatively high timing precision. The response over longer delays is shown as
an inset. The requisite timing precision can be adjusted with the $L/R$ time
of the SPD. This implementation of XOR in a single dendrite is analogous to
that observed in human pyramidal neurons in cortical layers two and three
[26].

Figure 19: Transfer functions for the RI dendrite. (a) Peak of integrated
signal as a function of the number of synapse events in a burst for several
values of $i_{\mathrm{b}}$ and frequency of burst inputs. (b) Peak of
integrated signal as a function of the frequency of input pulse train for
several values of $i_{\mathrm{b}}$ and number of pulses in the train.

In addition to discrete responses to few synapse events, the phenomenological
model can help us understand the transfer functions performed by synapses
coupled to dendrites when many pulses are incident. In Fig. 19 we analyze
these responses for various bias conditions, revealing the role of
$i_{\mathrm{b}}$ for sculpting dendritic behavior. Figure 19(a) shows the peak
of the accumulated signal in the dendritic integration loop as a function of
the number of synapse events input in a burst. Traces from the same color
family correspond to the same value of $i_{\mathrm{b}}$. Lighter traces
correspond to pulses input at a slower rate relative to the time constant of
the dendrite, while darker traces show the response for higher-frequency
inputs. The response is unsurprisingly larger for higher input frequencies and
can be broadly adjusted with the bias current. The sigmoid-like response with
threshold and saturation has been useful in various machine learning
applications. Similar responses are seen in Fig. 19(b) where the frequency of
the burst is now the quasi-continuous variable on the $x$-axis, and various
curves have different numbers of input pulses. Similar thresholding and
saturation behavior is evident, and the same qualitative behavior has been
demonstrated in experiments of related circuits [15].

If generated with conventional circuit models, the plots of Fig. 19 would have
required so much computational time as to be impracticable with the computer
used for this study. With the phenomenological model they are produced in
minutes so that many scenarios may be investigated and insight across
parameter space can be quickly gained.

##  5 Neurons

Whereas a synapse was shown to be a formed where a single-photon detector
couples to a dendrite, a loop neuron can be as simple as a dendrite with a
transmitter circuit providing output. Such a soma circuit is shown in Fig.
2(d). While all dendrites have a bias-dependent threshold for activation in
the receiving loop below which applied flux elicits no response, the soma is
distinguished in that it also has a threshold in the integration loop at its
output. Embedded in the integration loop is a superconducting thresholding
element referred to as a tron that switches from a zero-resistance state below
threshold to a high-resistance state when the current in the integration loop
reaches a critical value. Such comparators have been demonstrated as an
interface between superconductor and semiconductor electronics [27, 14]. In
the present case, accumulation of sufficient signal in the soma’s integration
loop will drive the tron to the voltage state, at which point a semiconductor-
based transmitter circuit will generate photons from a light-emitting diode
(LED). This is a spike event or action potential of a superconducting
optoelectronic loop neuron. When this threshold is reached and the tron
switches, the current in the somatic integration loop is purged, and
integration must begin again. In most dendrites in the system, any accumulated
signal in the integration loop is immediately coupled as flux into other
dendrites. However, the soma has no output other than the transmitter circuit,
so its state is only communicated to other elements of the network when
threshold is reached and light is produced.

Several approaches to transmitter circuits for loop neurons have been explored
[8], and here we base the phenomenological model on a circuit concept that
leverages MOSFETs in an essentially digital configuration. Upon reaching
threshold, activity is transferred to these semiconductor circuits, and to
understand their operation we must work backward from the light source itself.
To maintain “brevity” in the main thread of this article, we relegate the
details to appendices. Appendix E describes a rate-equation model to capture
the primary behavior of the light emitters, while Appendix F treats the
transmitter circuit from the thresholding element through the LED. Because the
subject of this work is to establish a phenomenological model, here we simply
extract the key message from those appendices: when a soma reaches threshold,
photons are produced with a few nanosecond delay followed by an exponential
probability distribution. The distribution is obtained numerically through
simulation of the transmitter circuit combined with the source rate equation
model. In the reduced phenomenological model that is the subject of this work,
each time a somatic integration current reaches threshold, a specified number
of photon-production times are drawn randomly from this distribution of a
delay followed by an exponential decay. These photon times are then assigned
randomly to the neuron’s downstream synapses and added to a list of times of
input synapse events. With this approach, the number of ODEs which must be
treated remains unchanged, and the added computational burden of treating the
transmitter circuit and sources is reduced to dealing with random number
generation only on the time steps at which the soma’s threshold is reached.
The indefatigable reader is directed to Appendices E and F for more details.

Figure 20: Schematic of a monosynaptic point neuron with excitatory synapse
($\mathsf{S_{e}}$), neuron cell body ($\mathsf{N}$), transmitter
($\mathsf{T}$) and refractory dendrite ($\mathsf{D_{r}}$).

In addition to this output thresholding behavior, a loop neuron requires a
means of establishing a refractory period so that spike events are indeed
discrete and their rate does not exceed the roughly 20 MHz at which the
synaptic SPDs can respond. A refractory period is accomplished by adding a
second dendrite to the soma. A schematic of this minimal point neuron is shown
in Fig. 20, where thie refractory dendrite is labeled $\mathsf{D_{r}}$. This
refractory dendrite is driven by a flux pulse identical in form to the SPD
response (Eq. 6) when the soma reaches threshold, and the signal accumulated
in the refractory integration loop is coupled back to the soma’s receiving
loop as inhibition. The time constant of the refractory dendrite’s integration
loop therefore sets the refractory period and maximum neuronal firing rate. To
achieve robust refraction, the coupling from the refractory dendrite back to
the soma is designed so that the inhibitory flux is as strong as possible
without driving activity on the negative-flux side of the $r$-shell. This can
be accomplished with knowledge of the values of $s_{\mathrm{max}}$ for the
refractory dendrite at its bias point as well as the values of threshold flux
$\phi_{\mathrm{th}}^{+}$ and $\phi_{\mathrm{th}}^{-}$. The refractory dendrite
is chosen to have a relatively small value of $\beta$ so that it saturates
quickly at $s_{\mathrm{max}}$ each time the soma fires. The coupling from the
refractory to the somatic dendrite is then given by
$J=(\phi_{\mathrm{th}}^{+}\,-\,\phi_{\mathrm{th}}^{-})/s_{\mathrm{max}}$.

The role of the refractory dendrite is to inhibit the neuron for a brief
duration after it spikes. Yet similar inhibition in response to activity is
often desirable on longer timescales. Homeostatic plasticity [28, 29, 30] is
one means by which the activity of neurons can be maintained in a useful
dynamic range. This self-regulatory behavior can be accomplished with feedback
analogous to refraction if further dendrites receive input when the soma fires
and couple inhibitively to the soma’s receiving loop. Such a homeostatic
dendrite is likely to have a longer time constant and also a larger
integration loop inductance than the refractory dendrite. In this way, each
time the neuron fires, a small inhibitory signal is fed back to the soma,
effectively increasing the threshold for the next firing. This signal will
decay over a longer time period, on the order of many interspike intervals.
Such a plasticity mechanism can be treated with the same phenomenological
model presented here.

Figure 21: Neuron firing rate versus applied flux. (a) Comparing several
values of the neuronal integration loop inductance parameter,
$\beta_{\mathrm{ni}}$. (b) Comparing several values of neuronal integration
loop time constant and bias current. For $i_{\mathrm{b}}=1.6$ and 1.8,
$\beta_{\mathrm{ni}}/2\pi=5\times 10^{3}$. For $i_{\mathrm{b}}=1.4$,
$\beta_{\mathrm{ni}}/2\pi=1.6\times 10^{4}$.

To begin analysis of neurons with this framework, we consider a neuron with no
dendritic tree in which all inputs are received directly at the soma. Such a
neuron is referred to as a point neuron. As a first means of gaining intuition
about the design of a point neuron, we consider the spiking behavior of a
simple loop neuron with steady-state flux coupled directly into the soma and a
single refractory dendrite to induce pulsatile behavior at the transmitter. In
Fig. 21 we show the rate of neuron spiking as a function of the steady-state
applied flux, $\phi$. Figure 21(a) shows the response for five values of the
neuronal integration loop inductance parameter, $\beta_{\mathrm{ni}}$. In all
cases the time constant of the refractory integration loop was 50 ns and the
value of the somatic integration loop threshold was set to
$0.3\,s_{\mathrm{max}}$. Smaller values of $\beta_{\mathrm{ni}}$ lead to
rapidly refilling of the integrated signal after it is purged upon reaching
threshold, while larger values require more time to accumulate signal back to
threshold, leading to slower firing activity for the same flux drive.

By contrast, in Fig. 21(b) the somatic integration loop inductance was fixed
at an intermediate value, while the operating bias point and somatic
integration loop time constant were adjusted. Lower biasing requires higher
input flux to achieve threshold and also leads to slower firing rates, while
higher bias reduces threshold and increases the output rate. Such a control
parameter is the primary means by which the soma can be reconfigured
dynamically during network operation, either in an unsupervised manner by
activity within the system or in a supervised manner by a user or control
system. The time constant cannot be changed dynamically, as it is set by the
resistance and inductance of the loop, which are fixed in fabrication.
However, the time constant can have an appreciable effect on the transfer
characteristics, as is most pronounced for the low-bias case where reducing
the time constant to 50 ns dramatically increases the flux threshold and
reduces the output rate at a given value of applied flux. All of these
parameters must be considered when designing neurons to play various roles in
a network.

As in the transfer functions of Fig. 19, the responses in Fig. 21 would have
been quite cumbersome to obtain with the conventional circuit equation
approach. Even at this basic level, the phenomenological model provides
utility in guiding the design principles of loop neurons across a broad
parameter space.

Figure 22: Schematic of the nine pixel classifier. The dendritic tree is
illustrated as are the three classes of inputs, $z$, $v$, and $n$.

We now consider an example of employing the phenomenological model to treat a
loop neuron with a more elaborate dendritic arbor to solve the standard
benchmark problem of the nine-pixel image classifier [31]. Figure 22 shows the
nine-pixel input and the manner in which the pixels can be argued to depict
three letters: ‘$z$’, ‘$v$’, and ‘$n$’. The problem as formulated for small
neural networks is to train the network to identify the image as the correct
letter even when the state of any one pixel is allowed to switch. This problem
is often solved with small neural networks trained with conventional
techniques such as backpropagation. Here we solve the problem with a single
neuron with a small dendritic tree, showing that such a tree can serve a
similar role to a feed-forward neural network [32] and to illustrate how even
relatively simple loop neurons become sophisticated computational processors,
as inspired by their biological counterparts [33]. We also achieve the desired
neural operation through inspection without a training algorithm, in the
spirit of biological vision systems that “hard-code” basic computational
primitives, such as Gabor filters matched to common spatial frequencies in
natural scenes [34].

The arbor that solves the problem is shown in Fig. 22. Each pixel is input to
a synapse as a single spike event: a spike event occurs at synapse $i$ if
pixel $i$ is active in the image under consideration, and no spike event is
input at that synapse if the corresponding pixel is not active. The various
dendrites are labeled by the logic operations they perform (Sec. 4). The
knowledge accessible at various stages of the tree is written above the
dendrites whose signal represents that information. The basic reasoning that
solves the problem can followed by stepping through the tree. For example,
starting at the top left, if pixel two or five is active but pixels four and
six are not both active, the letter cannot be $v$. Similar reasoning can be
applied to rule out $n$. If it is not $v$ and it is not $n$, it must be $z$.
At the last stage of the tree before the soma, three different dendrites know
whether their letter is present. By assigning these three dendrites different
time constants ($\alpha$ parameters in Eq. 2), they each evoke different
numbers of pulses from the neuron. An output of one spike event from the
neuron informs us a $z$ has been presented, two correspond to $v$, and three
to $n$. The neuron under consideration is shown to accomplish the task even in
the presence of 1 ns timing jitter due to the transmitter circuit and quantum
dot light source (Appendices E and F).

Figure 23: Time traces for the nine-pixel classifier neuron for the case of
$n_{4}$, the image of $n$ with pixel four switched. (a) The input synaptic
flux from all four active synapses. The inset shows higher temporal resolution
to display the timing jitter due to the model of the quantum-dot light source.
(b) The total applied flux to the neuronal receiving loop, including the
refractory contribution. (c) The integrated current in the neuronal
integration loop. (d) The flux generated by the downstream SPD receiving
photons upon neuronal firing.

Time traces from a single classification instance ($n_{4}$) are shown in Fig.
23. The four input synapse events are shown in Fig. 23(a), and the timing
jitter due to the transmitter circuit and light source is evident in the
inset. The total flux input to the neuronal receiving loop is shown in Fig.
23(b), as are the thresholds for activity due to positive and negative flux.
The accumulated signal in the neuronal integration loop is shown in Fig.
23(c), along with the value of the neuronal integration loop threshold. The
times of neuronal spike events are shown as black crosses. In this model, the
neuron is given one output synapse, and the flux present at the output of that
synapse due to photon detection events from neuronal firing are shown in Fig.
23(d).

Figure 24: Output of the nine-pixel classifier neuron versus input image. (a)
The peak of the integrated signal in the final dendrites that represent the
three letters. (b) The number of spikes out as a function of the input image.

The full set of 30 inputs is shown in Fig. 31 in Appendix G. The results of
all instances of the classification task are shown in Fig. 24. The maximum
level of integrated signal in all three final dendrites is shown in Fig. 24(a)
for all thirty presented images. The number of spikes produced by the neuron
is shown in Fig. 24(b). It can be seen that no signal is generated in any of
the dendrites except the one corresponding to the presented letter, even in
the presence of one-pixel noise, and the neuron correctly produces the number
of spike outputs to communicate the result of its calculation to other neurons
that may be present in the network.

##  6 Discussion

We have seen that the computational circuitry of loop neurons consists of a
network of interacting dendrites. Within this framework, a dendrite is a SQUID
with output coupled to a current-storage loop. The bias point of the SQUID
determines the threshold input flux required to initiate activity, and the
storage capacity of the output loop provides a saturating nonlinearity. The
various circuit parameters provide means to adjust the response
characteristics across broad operating ranges. Learning and homeostatic
adaptation can be accomplished with similar circuits providing feedback,
dynamically adjusting bias points through coupled flux that can be stored
perpetually and adjusted in small or large increments. The final stage of
dendritic processing culminates in the soma, the neuron cell body. When this
final dendrite reaches threshold, something different happens with regard to
physical hardware: an amplifier drives a semiconductor light source. The
threshold signal is the final stage of the computational process; the
production of light is a binary action potential, and the physical
transduction to photons is chosen to enable communication to many destinations
across length scales that result in wiring parasitics that render electronic
signals impracticable. A neuron is a network of coupled dendrites with local
electronic communication. A spiking neural network extends interconnectivity
to longer distances with optical communication.

Simulating loop neurons with the phenomenological model brings significant
speed increases relative to the conventional method of simulating the first-
principles circuit model. The speed of simulation can increase by a factor of
ten thousand while retaining a $\chi^{2}$ of $10^{-4}$ when simulating a
single dendrite for 12 µs. When simulating a large network, the circuit-model
approach is impracticable due to the time required to run calculations over
even a relatively short duration of network activity. The phenomenological
model thus enables design of more components, exploration of a wider range of
parameter space, and treatment of larger networks than would otherwise be
possible. Still, simulating 12 µs of activity from a single dendrite requires
around 1 s of CPU time. A physically instantiated network on a small chip
(roughly 1 cm $\times$ 1 cm) would comprise roughly 10,000 neurons and half a
million synapses and dendrites. Linearly extrapolating the time requirements,
it would take over a year to simulate 1 ms of activity of this network on the
same workstation used in this work. Two conclusions can be drawn: 1) the
phenomenological model must be used as a tool to facilitate further increases
in speed of design, and 2) constructing the physical hardware of loop neurons
could bring a computational platform of immense power and utility. Regarding
the first point, by using the phenomenological model to design a variety of
dendrites, coupled-dendrite computational blocks, and a variety of neurons,
libraries of optimized components can be assembled. The model can be used to
efficiently design these components and reduce their behavior to stored
transfer functions to achieve a further level of abstraction that can be used
to model behavior of larger networks without explicitly stepping through the
ODEs, as was done for the logic gates, burst, and rate transfer functions in
Sec. 3. Regarding the second point, it would take a year on a workstation with
a 3.7 GHz processor to simulate one millisecond of a network that would fit on
a chip of similar size and power consumption to the processor, including
cryogenic cooling. The increase in computational speed that would result from
realizing SOEN hardware would be 32 billion for even a modest network. SOENs
therefore appear to offer a path to neuromorphic supercomputing.

If the aspiration of loop neurons is to constitute systems of exceptional
scale and complexity, might the simplicity of the dendritic building block
limit the dynamical repertoire? We hope this model helps answer this question.
It is known that even simple systems with simple rules for propagation, such
as cellular automata, can give rise to behavior of great sophistication [35,
36]. It has been argued that similar concepts can be applied as a starting
point for physics, with the rich, natural world emerging based on the
interactions of fundamental nodes [37]. Similarly, elements as simple as two-
state spins interacting with nearest neighbors (Ising models) can give rise to
phase transitions and critical phenomena, including crucial long-range
correlations [38, 39] as well as attractor dynamics for associative memory
storage and retrieval [40, 41, 42, 43]. By comparison, the dendrites studied
here are multi-dimensional and nuanced. The nonlinear spatio-temporal
convolutions occurring in each neuron’s synapto-dendritic tree provide a deep
repository that can inform the neuron’s behavior. Transmission of action
potentials to many destinations at light speed enable complex network
topologies far beyond nearest-neighbor interactions. The use of
superconducting circuits allows signal retention across a broad range of time
scales through the choice of the $\alpha$ parameters (the $\alpha$
distribution) that specify the leak rates in Eq. 2, as has been demonstrated
experimentally [15]. The range of integration loop inductances (the $\beta$
distribution) in conjunction with the bias currents will establish the dynamic
range of the network. With near and distant connectivity as well as a wide
spectrum of dissipation times and an enormous dynamic range, SOENs are likely
to achieve the long-range spatiotemporal correlations present in the critical
states that optimize neural information processing [44, 45, 46, 47]. The
hierarchy of spatial connectivity in conjunction with the hierarchy of
information retention times and response magnitudes appears excellent for
enabling the fractal use of space and time that supports information
integration and cognition [48, 49, 50, 51, 52, 53]. Yet the simplicity of
constructing the majority of computational grey matter from similar building-
block components brings an advantage in modeling and technological
implementation.

The significance of dendrites in neural computation is well documented [54,
32, 55, 56, 57, 58, 59], yet the incorporation in neuromorphic hardware has
received proportionally less attention [60, 61, 62, 63, 64]. The dendrites
studied here may provide valuable means for implementing credit assignment in
training algorithms that utilize local information in conjunction with
population information [65, 66, 67, 68, 69], potentially leveraging much of
what is know from neuroscience about the key role of dendrites in learning
[70, 71, 72, 73, 74, 75, 76, 77] and overcoming a major obstacle to widespread
adoption of spiking neural networks for artificial intelligence. The presence
of continuous dendritic signals that are continuous in time without erasure
following a spike may provide new methods for training spiking neural networks
that are not available based on spiking activity alone. Perhaps these
continuous signals will be useful for constructing cost functions and training
networks with variants of the backpropagation technique, a feat which has been
difficult in the case of spiking neurons outside of the rate-coding domain,
and a subject which is significant in bridging machine learning and
neuroscience [78].

While the model presented here is interesting to explore, it is intended as a
tool in a larger project. The objective of such research at this stage is to
determine if SOENs are indeed as promising for future study and worthy of
appreciable investment as our hypothesis contends. We hope this model will
help answer several questions: Do the neuromorphic circuit principles
demonstrated in these superconducting optoelectronic embodiments bestow
systems with the dynamic interplay of structure and function that scaffolds
successful neural systems? Do device and circuit features such as the
achievable breadth of the $\alpha$ and $\beta$ distributions, the particular
dendritic nonlinearities, and the available forms of plasticity bring benefits
in support of cognition? Does refraction of the neuron without erasure of
information in the dendritic tree offer advantages in information processing,
lead to new forms of neural coding, or offer new means of training? How should
the structure of the dendritic trees and the structure of the network be co-
designed to achieve desired operations? The objective of the model presented
here is to enable computational studies that answer these questions.

With the phenomenological model as a guide, we can compare and contrast loop
neurons with other complex systems. Spin systems have received considerable
attention over decades due to both their simplicity as well as the emergence
of interesting phenomena such as phase transitions and spontaneous symmetry
breaking. Common spin models include the Ising model, where each spin $s$
takes a binary, scalar value of -1 or 1. Coupling is typically between nearest
neighbors and is also a scalar value, $J_{ij}$, typically with the symmetry
$J_{ij}=J_{ji}$. The $x-y$ model and the Heisenberg model extend the spin to
have two or three vector components, maintaining $|s|=1$. In the present
context, we can think of the current stored in the integration loop of
dendrite $i$, $s_{i}$, as a spin. It is more like an Ising model than an
$x$-$y$ or Heisenberg in that it is a scalar value, but in general it is
analog instead of binary and can take a large number of values between zero
and one. Coupling is through the scalar quantities $J_{ij}$ (Eq. 3) and is
manifest through transformers. In this case, there is no requirement that
$J_{ij}=J_{ji}$, and the spatial limitations regarding which spins can couple
are relaxed. Coupling is mediated by a transformer circuit rather than an
exchange interaction, so flexibility in the adjacency matrix is available. The
values of $J_{ij}$ are fixed in time, but the sign can be positive or negative
so that excitation or inhibition are possible. While the $J_{ij}$ are fixed in
circuit fabrication, the interaction between any two spins is dynamically
reconfigurable either by an external influence (experimenter) or through
network activity (plasticity) by changing applied bias currents to change the
$r$-shell on which the dendrite resides. Furthermore, the magnitude of the
spins in loop neurons is not fixed in time, but is a dynamical quantity that
grows through interactions and decays passively, following the leaky
integrator formalism. In general, the threshold for one spin to induce another
is non-zero, an important nonlinearity familiar from neuroscience.

Perhaps the most interesting connection to spin systems is to glasses in which
disorder and frustration play central roles. Frustration in spin glasses
occurs when any spin cannot achieve a configuration in which the interaction
energy to all of its neighbors is minimized. This phenomenon leads to energy
landscapes with many peaks and valleys and a large number of metastable
configurations [79]. Competing interactions are a hallmark of complex systems
and result in interesting and often computationally useful dynamics as the
system traverses the landscape of metastable states [80]. In the present case,
the dendritic tree of each loop neuron can be compared to a spin glass in
which each dendrite has multiple interactions with competing signs of
coupling. Frustration can be quantified for dendrite $i$ as

$f_{i}=\frac{\left(\sum_{j=1}^{n_{+}}J_{ij}\,s_{j}\right)\;\left(\sum_{j=1}^{n_{-}}\left|J_{ij}\right|\,s_{j}\right)}{\left(\sum_{j=1}^{n_{+}}J_{ij}\,s_{j}^{\mathrm{max}}\right)\;\left(\sum_{j=1}^{n_{-}}\left|J_{ij}\right|\,s_{j}^{\mathrm{max}}\right)},$
(7)

where the sum to $n_{+}$ runs over the excitatory inputs, the sum to $n_{-}$
runs over inhibitory inputs, and $s_{j}^{\mathrm{max}}$ is the saturation
level of the signal in dendrite $j$. This value for frustration is zero when
the sum of excitatory or inhibitory inputs is zero and reaches a maximum value
of one when all excitatory and inhibitory inputs are at their maximum values.
The frustration of neuron $p$ can be quantified by summing over its arbor,
$\mathcal{F}_{p}=\sum_{i}f_{i}$. Each neuron’s dendritic tree has qualitative
features in common with a dynamical spin glass, with randomness and transitory
inhibitory and excitatory interactions leading to ever-shifting states of
competition and frustration, as is also evident in biological neural systems
[81]. The frustration of the network can be obtained by summing once more over
neurons, $\mathscr{F}_{N}=\sum_{p}\mathcal{F}_{p}$. Future work will explore
the relationship between frustration, critical phenomena, and useful neural
computation in the context of SOENs.

Above the level of the dendritic tree, at the level of coupled neurons,
analogies to spin systems are less relevant, and comparison to pulse-coupled
oscillators and their more complex cousins, spiking neurons, are
straightforward. Superconducting optoelectronic loop neurons are spiking
neurons that produce binary, pulsatile communication signals upon being driven
to an internal threshold. The optical communication signals that facilitate
inter-neuron interaction are not strictly based on spatial location (nearest
neighbors, next-nearest neighbors, etc.), but can be engineered with great
flexibility through an adjacency matrix physically realized by photonic
interconnects. As compared to typical point neurons, a complex, active
dendritic tree appears central to the function and construction of loop
neurons [17]. Yet at the base of that tree resides a soma that sums inputs and
registers a threshold. While the signal in that final integration loop of the
soma is purged upon reaching threshold, the signals stored in the loops of the
rest of the arbor need not change upon neuronal firing. As compared to
archetypal leaky integrate-and-fire neurons, the lack of erasure following
firing is distinct. This trait allows a single synapse event to push a neuron
into a dynamical orbit and allows information retention regardless of firing
state.

In computational neuroscience a list of discrete spike times is sometimes
considered a complete description of the behavior of a neuron or population in
a given context. Such a list is available for networks of loop neurons, yet
continuous variables specifying the states of all dendritic loops can also be
used to gain more granular information about the system, much like retaining
the time-continuous values of the membrane potentials on all dendrites of all
neurons throughout the network. While loop neuron spike times, bursts, and
rates are crucial to network activity, the state of any given neuron or
population may also be specified by the signal in all loops, providing richer
depiction of state space, high temporal resolution even between spike events,
and continuous time-series representation suitable for defining distances and
overlaps between network states as well as correlation functions between
components at different times. The state of dendrite $i$ attached to neuron
$p$ is specified by $s_{i}$. This state description can be extended to neurons
as $S_{p}=[s_{1},\dots,s_{n}]_{p}$: the state of the neuron is the state of
all of its dendrites arranged as a vector. The state of the network can then
be defined in full detail as the full vector including the state of all
dendrites, $\mathcal{S}=[S_{1},\dots,S_{N}]$, or the state of the each neuron
can be reduced to an average for computational efficiency when necessary. The
distance between two states can then be defined as the Euclidean distance
between two vectors, applicable to single neurons or networks. Such a metric
is valuable in many contexts, including in quantifying sensitivity to
variation in initial conditions to determine when a network is in a chaotic
regime.

We can then discuss correlations between dendrites with themselves at
different times, with other dendrites in the same neuronal arbor, and with
other dendrites attached to different neurons. For example, a self-correlation
function of the form

$G_{i}(t,t^{\prime})=s_{i}(t)\,s_{i}(t^{\prime})$ (8)

may be useful for identifying temporal coherence related to recurrent
activity. Similarly, a cross-correlation function between two dendrites,

$G_{ij}(t,t^{\prime})=s_{i}(t)\,s_{j}(t^{\prime}),$ (9)

may be useful in identifying functional coalitions within a neuron with
dendrites $i$ and $j$ within the same arbor, or it may be useful for
identifying synchronized activity across the network when they are attached to
different neurons. Extending this concept to neurons, we can write down a
temporal correlation function of a neuron with itself,

$G_{p}(t,t^{\prime})=S_{p}(t)\,\cdot\,S_{p}(t^{\prime}),$ (10)

or with other neurons,

$G_{pq}(t,t^{\prime})=S_{p}(t)\,\cdot\,S_{q}(t^{\prime}).$ (11)

In Eqs. 10 and 11, the vector dot product can be taken to reduce the
information to a single, scalar metric for temporal correlations. The cross-
correlation function of Eq. 11 contains information about correlations across
space and time, as does the dendritic cross-correlation function of Eq. 9 on a
shorter length scale. Such functions inform us as to how inputs are
represented and stored in correlations, and from these functions we can obtain
quantities analogous to the susceptibility which provide information about
phase transitions and critical phenomena. These concepts may be further
extended to a continuum field theory, potentially conducive to formal analysis
that is gaining traction in neural-network applications [82, 83].

In spin [39], neural [84], and other complex systems, free energy is a
powerful concept capturing a system’s competing drives to minimize energy and
maximize entropy. In an elementary construction, the free energy is given by
$\mathcal{F}=\mathcal{U}-T\mathcal{S}$, with $\mathcal{U}$ the total energy of
the system, $T$ the temperature, and $\mathcal{S}$ the entropy. Within the
present framework, the energy of a single dendrite is given by
$L\,I^{2}/2\propto\beta s^{2}$, and the total energy can be obtained by
summing over dendrites. The entropy can be obtained from the logarithm of the
number of microstates at a given value of total energy, which is a well-
defined if tedious problem in combinatorics. Alternatively, an approach in
terms of a partition function may be tractable. The temperature of the system
may be conceived as the actual, physical temperature in units of kelvins, with
the $T\mathcal{S}$ term growing due to stochastic switching of JJs close to
$I_{c}$, or temperature may be engineered as a more abstract quantity
introduced via random photonic inputs to synapses. The conceptual foundation
of the loop-neuron system in terms of simple component state variables as
formulated with this phenomenological model offers multiple routes to connect
to a deeper free-energy formulation.

From this list of comparisons and extensions it is evident that loop neurons
and superconducting optoelectronic networks have potential to lead to physical
systems of immense complexity and useful computational functionality. The
phenomenological model developed here may serve in the near term to facilitate
design of useful circuits and networks for computation. In the longer term the
framework may enable more complete theoretical development of the hardware as
a complex physical system, as a tool for artificial intelligence, as a
platform for hypothesis testing in neuroscience, and as a means to devise
systems with intelligence exceeding our own.

##  Appendix A Circuit Equations for the RI Dendrite

Figure 25: Circuit diagram of the RI dendrite with relevant circuit parameters
labeled.

A circuit diagram of the RI dendrite is shown in Fig. 25. Kirchoff’s current
law gives $I_{\mathrm{b}}=I_{1}+I_{2}+I_{\mathrm{di}}$. Summing voltages
around the R loop gives

$\frac{\Phi_{0}}{2\pi}\frac{d\,\delta_{1}}{dt}+L_{1}\frac{d\,I_{1}}{dt}-L_{2}\frac{d\,I_{2}}{dt}-\frac{\Phi_{0}}{2\pi}\frac{d\,\delta_{2}}{dt}-\frac{d\,\Phi_{\mathrm{r}}}{dt}=0,$
(12)

where $\delta_{i}$ is the phase of the superconducting wave function across
the $i$th junction. We introduce the following dimensionless variables [20]:

$\begin{split}i_{x}&\equiv I_{x}/I_{c},\\\\[10.0pt] \beta_{x}&\equiv\frac{2\pi
L_{x}I_{c}}{\Phi_{0}},\\\\[10.0pt]
\phi&\equiv\Phi_{\mathrm{r}}/\Phi_{0},\\\\[10.0pt] \alpha&\equiv
r_{\mathrm{di}}/r_{\mathrm{j}},\\\\[10.0pt] \tau&\equiv\frac{2\pi
r_{\mathrm{j}}I_{c}}{\Phi_{0}}\,t\equiv\omega_{c}t.\end{split}$ (13)

$\omega_{c}$ is the characteristic Josephson frequency, and $r_{\mathrm{j}}$
is the Josephson junction shunt resistance in the resistively and capacitively
shunted junction (RCSJ) model [21, 22, 23, 20]. A dimensionless current
$i_{x}$ is defined for each branch of the circuit ($i_{\mathrm{b}}$, $i_{1}$,
$i_{2}$, and $i_{\mathrm{di}}$), while a dimensionless screening parameter,
$\beta_{x}$, is defined similarly for each inductor. We also use the notation
$\beta_{\mathrm{r}}=\beta_{1}+\beta_{2}$ and
$\bar{\beta}=\beta_{1}\beta_{2}+\beta_{1}\beta_{\mathrm{di}}+\beta_{2}\beta_{\mathrm{di}}$.
We assume all JJs have identical critial current, $I_{c}$. Algebraic
manipulation of Eq. 12 leads to

$\frac{d\,i_{1}}{d\tau}=\frac{1}{\beta_{\mathrm{r}}}\left(\frac{d\,\delta_{2}}{d\tau}-\frac{d\,\delta_{1}}{d\tau}+2\pi\frac{d\,\phi}{d\tau}\right)+\frac{\beta_{2}}{\beta_{\mathrm{r}}}\left(\frac{d\,i_{\mathrm{de}}}{d\tau}-\frac{d\,i_{\mathrm{di}}}{d\tau}\right),$
(14)

where the relation $d/dt=\omega_{c}d/d\tau$ has been used. Summing voltages
around the I loop gives

$\frac{\Phi_{0}}{2\pi}\frac{d\,\delta_{2}}{dt}+L_{2}\frac{d\,I_{2}}{dt}-L_{\mathrm{di}}\frac{d\,I_{\mathrm{di}}}{dt}-r_{\mathrm{di}}I_{3}=0.$
(15)

Moving to dimensionless units and making use of Eq. 14 we can state the
primary equation of motion governing the circuit of Fig. 25:

$\bar{\beta}\,\frac{d\,i_{\mathrm{di}}}{d\tau}=\beta_{1}\frac{d\,\delta_{2}}{d\tau}+\beta_{2}\frac{d\,\delta_{1}}{d\tau}-2\pi\frac{d\,\phi}{d\tau}+\beta_{1}\beta_{2}\frac{d\,i_{\mathrm{b}}}{d\tau}-\alpha\beta_{\mathrm{r}}i_{\mathrm{di}}.$
(16)

Equation 16 must be solved in a system with an ODE for each of the JJs of the
form

$\beta_{c}\frac{d^{2}\,\delta_{x}}{d\tau^{2}}=i_{x}-\mathrm{sin}(\delta_{x})-\frac{d\,\delta_{x}}{d\tau}.$
(17)

Here $\beta_{c}=2\pi I_{c}r_{\mathrm{j}}^{2}c_{j}/\Phi_{0}$ is the Stewart-
McCumber parameter, with $c_{j}$ the capacitance of a JJ in the circuit. Such
an expression can be derived in the framework of the RCSJ model [21, 22, 23,
20]. This second-order ODE must be converted to two first order ODEs for each
of the JJs. The system of coupled ODEs combining Eqs. 16 and 17 thus results
in a system of five ODEs that can be solved numerically for arbitrary
parameter values, bias conditions, and flux drives. After the system of ODEs
has been integrated in time, $i_{1}$ can be obtained as

$i_{1}=\frac{1}{\beta_{\mathrm{r}}}\left(\delta_{2}-\delta_{1}+2\pi\phi\right)+\frac{\beta_{2}}{\beta_{\mathrm{r}}}\left(i_{\mathrm{de}}-i_{\mathrm{di}}\right),$
(18)

and $i_{2}=i_{\mathrm{b}}-i_{1}-i_{3}$. This numerical model has been
implemented in Python and the system of ODEs solved with the solve_ivp
function from SciPy, employing a Runge-Kutta integration method with fifth-
order accuracy and an adaptive time grid.

In constructing the phenomenological model, we have essentially replaced Eq.
16 with the leaky integrator Eq. 2 and abstracted away the driving terms into
the phenomenological rate function. By comparing Eq. 16 and 2 and noting that
in the main text we have made the replacement
$i_{\mathrm{di}}\,\rightarrow\,s$, we see that, roughly speaking,

$\beta_{1}\frac{d\,\delta_{2}}{d\tau}+\beta_{2}\frac{d\,\delta_{1}}{d\tau}-2\pi\frac{d\,\phi}{d\tau}\,\rightarrow\,r(s,\phi;i_{\mathrm{b}}),$
(19)

aside from differences in $\beta$ prefactors, and ignoring the temporal
derivative of the bias current. While the two differential equations tracking
the current in the dendritic integration loop are qualitatively similar, by
moving to the phenomenological model we no longer require tracking the time
derivatives of the phases of the JJs and the applied flux, which evolve on the
picosecond time scale, and instead require only the use of previously
calculated rate arrays as look-up tables. The demonstrated speed benefits
arise from this replacement.

##  Appendix B Circuit Equations for the RTTI Dendrite

The primary equation of motion for the RTTI dendrite is

$\beta_{\mathrm{di}}\frac{d\,i_{\mathrm{di}}}{d\tau}=\frac{d\,\delta_{4}}{d\tau}-\alpha\,i_{\mathrm{di}}.$
(20)

Equation 20 must be solved in a system with an ODE of the for given by Eq. 17
for each of the four JJs in the circuit. Upon reducing the second-order odes
to first order, the system comprises nine coupled ODEs. The remaining currents
in the circuit can be obtained from

$\begin{split}i_{4}&=\beta_{4}^{-1}\left(\delta_{3}-\delta_{4}\right)+i_{\mathrm{b3}}-i_{\mathrm{di}},\\\\[10.0pt]
i_{3}&=\bar{\beta}^{-1}\left[\beta_{1}\left(\delta_{2}-\delta_{3}\right)+\beta_{2}\left(\delta_{1}-\delta_{3}-2\pi\phi\right)+\beta_{1}\,\beta_{2}\,i_{\mathrm{b1}}\right]\\\
&+i_{\mathrm{b2}}+i_{\mathrm{b3}}-i_{4}-i_{\mathrm{di}},\\\\[10.0pt]
i_{2}&=\frac{1}{\beta_{1}+\beta_{2}}\,\left(\delta_{1}-\delta_{2}-2\pi\phi\right)\\\
&+\frac{\beta_{1}}{\beta_{1}+\beta_{2}}\,\left(i_{\mathrm{b1}}+i_{\mathrm{b2}}+i_{\mathrm{b3}}-i_{3}-i_{4}-i_{\mathrm{di}}\right),\\\\[10.0pt]
i_{1}&=i_{\mathrm{b1}}+i_{\mathrm{b2}}+i_{\mathrm{b3}}-i_{2}-i_{3}-i_{4}-i_{\mathrm{di}}.\end{split}$
(21)

##  Appendix C Obtaining the Phenomenological Model in Dimensionless Units

In this Appendix we begin with slightly different and more explicit variable
names that used in the main text. We proceed with cluttered notation through
the derivation to reduce ambiguity, and all notation will be simplified at the
end.

Figure 26: Full circuit diagram for the RI dendrite labeling additional
parameters relevant to the derivation of the model in dimensionless units.

The key dynamical variable of each dendrite is the integrated current in the
DI loop, $I_{\mathrm{di}}$ (see Fig. 26). The foundational postulate of this
work is that time-evolution of this quantity can be captured in a
phenomenological model of the form

$\frac{d\,I_{\mathrm{di}}}{dt}=I_{\mathrm{fq}}\,R_{\mathrm{fq}}\\!\left(\Phi_{\mathrm{r}},I_{\mathrm{di}};I_{\mathrm{b}}\right)-\frac{I_{\mathrm{di}}}{\tau_{\mathrm{di}}},$
(22)

where $I_{\mathrm{fq}}=\Phi_{0}/L_{\mathrm{di}}$ is the current associated
with each fluxon, $R_{\mathrm{fq}}$ is the rate of fluxon production as a
function of $\Phi_{\mathrm{r}}$, the flux applied to the receiving loop, and
$I_{\mathrm{di}}$, the current accumulated in the integration loop.
$\tau_{\mathrm{di}}=L_{\mathrm{di}}/R_{\mathrm{di}}$ is the time constant
governing the leak of the integrated signal. As in Appendix A, we convert this
equation to dimensionless units:

$\beta_{\mathrm{di}}\frac{d\,i_{\mathrm{di}}}{d\tau}=r_{\mathrm{fq}}\left(\phi,i_{\mathrm{di}};i_{\mathrm{b}}\right)-\alpha_{\mathrm{di}}\,i_{\mathrm{di}}.$
(23)

The dimensionless parameters $i_{\mathrm{di}}$, $\beta_{\mathrm{di}}$,
$\alpha_{\mathrm{di}}$, and $\tau$ have been defined in Eq. 13. The
dimensionless rate of fluxon production is given by

$r_{\mathrm{fq}}=\frac{R_{\mathrm{fq}}}{\omega_{c}/2\pi},$ (24)

with $\omega_{c}=2\pi r_{\mathrm{j}}I_{c}/\Phi_{0}$.

Each dendrite in the system will obey an equation of the form of Eq. 23, and
knowledge of the system is complete when the values of $i_{\mathrm{di}}$ for
all dendrites is obtained. However, to solve the complete system of coupled
ODEs we need to calculate how the values of $i_{\mathrm{di}}$ from input
dendrites convert to coupling flux in the receiving dendrite. Thus, the next
step is to determine the total applied flux to the DR loop of dendrite $i$ in
terms of the values of $i_{\mathrm{di}}$ present in the input dendrites, which
we index with $j$. To make space for indices, we move the part labels to
superscripts.

The total applied flux to the DR loop of dendrite $i$ is

$\Phi_{i}^{\mathrm{dr}}=M_{i}^{\mathrm{dr|dc}}\,I_{i}^{\mathrm{dc}}.$ (25)

The notation $M_{i}^{\mathrm{dr|dc}}$ refers to the mutual inductance coupling
from the DC loop to the DR loop on dendrite $i$. Explicitly,
$M_{i}^{\mathrm{dr|dc}}=k_{i}^{\mathrm{dr|dc}}\,\left(L_{i}^{\mathrm{dr}}\,L_{i}^{\mathrm{dco}}\right)^{1/2}$,
with $k_{i}^{\mathrm{dr|dc}}$ the transformer efficiency. In Eq. 25,
$I_{i}^{\mathrm{dc}}$ is the net induced current in the collection coil input
to dendrite $i$. This current obeys the differential equation

$\frac{d\,\Phi_{i}^{\mathrm{dc}}}{dt}-L_{i}^{\mathrm{dc}}\,\frac{d\,I_{i}^{\mathrm{dc}}}{dt}=0,$
(26)

which can be integrated to obtain

$I_{i}^{\mathrm{dc}}=\frac{1}{L_{i}^{\mathrm{dc}}}\,\Phi_{i}^{\mathrm{dc}},$
(27)

where
$L_{i}^{\mathrm{dc}}=L_{i}^{\mathrm{dcp}}+L_{i}^{\mathrm{dco}}+\sum_{j=1}^{n}L_{j}^{\mathrm{dc}}$,
with $L_{j}^{\mathrm{dc}}$ the inductance of the transformer input to the DC
coil connecting dendrite $j$ to dendrite $i$, $L_{i}^{\mathrm{dco}}$ the
output inductance of the DC loop, which is the input to the transformer to the
DR loop, and $L_{i}^{\mathrm{dcp}}$ any additional parasitic inductance on the
DC coil. The sum is over the $n$ inputs to the DC loop of dendrite $i$. The
total flux applied to the DC loop is

$\Phi_{i}^{\mathrm{dc}}=\sum_{j=1}^{n}M_{j}^{\mathrm{dc|do}}\,I_{j}^{\mathrm{di}},$
(28)

with
$M_{j}^{\mathrm{dc|do}}=k_{j}^{\mathrm{dc|do}}\,(L_{j}^{\mathrm{dc}}\,L_{j}^{\mathrm{do}})^{1/2}$.
Inserting Eq. 28 into Eq. 27 we have

$I_{i}^{\mathrm{dc}}=\frac{1}{L_{i}^{\mathrm{dc}}}\,\sum_{j=1}^{n}k_{j}^{\mathrm{dc|do}}\,(L_{j}^{\mathrm{dc}}\,L_{j}^{\mathrm{do}})^{1/2}\,I_{j}^{\mathrm{di}}.$
(29)

Now inserting Eq. 29 into Eq. 25 we obtain

$\Phi_{i}^{\mathrm{dr}}=k_{i}^{\mathrm{dr|dc}}\frac{(L_{i}^{\mathrm{dr}}\,L_{i}^{\mathrm{dco}})^{1/2}}{L_{i}^{\mathrm{dc}}}\sum_{j=1}^{n}k_{j}^{\mathrm{dc|do}}(L_{j}^{\mathrm{dc}}\,L_{j}^{\mathrm{do}})^{1/2}\,I_{j}^{\mathrm{di}}.$
(30)

We now define the dimensionless quantities

$\begin{split}s_{i}&\equiv
i_{i}^{\mathrm{di}}=\frac{I_{i}^{\mathrm{di}}}{Ic},\\\\[10.0pt]
\phi_{i}&\equiv\phi_{i}^{\mathrm{dr}}=\frac{\Phi_{i}^{\mathrm{dr}}}{\Phi_{0}},\\\\[10.0pt]
J_{ij}&\equiv\frac{k_{i}^{\mathrm{dr|dc}}\,k_{j}^{\mathrm{dc|do}}}{2\pi}\,\frac{(\beta_{i}^{\mathrm{dr}}\,\beta_{i}^{\mathrm{dco}}\,\beta_{j}^{\mathrm{dc}}\,\beta_{j}^{\mathrm{do}})^{1/2}}{\beta_{i}^{\mathrm{dc}}},\end{split}$
(31)

where the various $\beta$ parameters are related to the inductances as in Eq.
13. For simplicity of notation we also specify
$\beta_{i}^{\mathrm{di}}\rightarrow\beta_{i}$ and $r_{\mathrm{fq}}\rightarrow
r$, and we drop the subscript $i$ for notational simplicity. We can now recast
the phenomenological model of Eq. 23 in its final form:

$\beta\frac{d\,s}{d\tau}=r\left(\phi,s;i_{b}\right)-\alpha\,s,$ (32)

where the coupling flux from the $n$ dendrites indexed by $j$ into dendrite
$i$ is given by

$\phi_{i}=\sum_{j=1}^{n}\,J_{ij}\,s_{j}.$ (33)

Equations 32 and 33 are the complete phenomenological model of a system of
coupled superconducting loop dendrites. The $J_{ij}$ are coupling parameters
that can be chosen across a broad range in the circuit layout. To extend the
model to account for neurons, certain dendrites are assigned to be somas,
wherein a transmitter circuit is activated and light is produced when the
value of $s$ reaches a specified threshold. The times of these threshold
events are spike times and lead to synapse events on that neuron’s downstream
synapses. To extend the model to account for synapses, certain dendrites are
assigned to receive input flux from an SPD as in Eq. 6 at a time corresponding
to the input neuron’s spike times plus a delay to account for the activity of
the transmitter circuits and light sources, as discussed in Appendices E and
F.

##  Appendix D Additional Data Comparing Phenomenological and Circuit Models

Here we tabulate the $\chi^{2}$ and simulation time ratio data from Figs. 6,
7, and 8. Figure 27 shows further data from square pulse driving functions.
The data points come from cases of 10, 20, 40, 80, and 160 square pulses. The
$x$-axis is the duration of the simulated time.

Figure 27: Summary of errors in comparing the phenomenological model to the circuit equations as a function of the simulated time, $T_{\mathrm{sim}}$. In these simulations, $\tau_{\mathrm{di}}=250$ µs, and the drive signals were random square-pulse inputs with 10, 20, 40, 80, and 160 pulses. (a) The values of $\chi^{2}$ for several values of time step and integration loop inductance. (b) The ratio of required time for the computation. $t_{\mathrm{soen}}$ is the time required to step through the phenomenological model, while $t_{\mathrm{ode}}$ is the time required to solve the system of circuit ODEs. The longest simulation treated $T_{\mathrm{sim}}=12.32$ µs of simulated time. Table 1: Values obtained from the simulation of Fig. 6, the case of the linear ramp input with 100 ps time step. $\beta_{\mathrm{di}}/2\pi$ | $\tau_{\mathrm{di}}$ | $\chi^{2}$ | $t_{\mathrm{ode}}/t_{\mathrm{soen}}$
---|---|---|---
$10^{2}$ | 10 ns | $6.76\times 10^{-5}$ | $1.40\times 10^{3}$
$10^{2}$ | 50 ns | $8.01\times 10^{-5}$ | $1.29\times 10^{3}$
$10^{2}$ | 250 ns | $9.95\times 10^{-5}$ | $1.17\times 10^{3}$
$10^{2}$ | 1.25 µs | $8.14\times 10^{-5}$ | $1.16\times 10^{3}$
$10^{3}$ | 10 ns | $3.83\times 10^{-6}$ | $1.78\times 10^{3}$
$10^{3}$ | 50 ns | $5.07\times 10^{-6}$ | $1.59\times 10^{3}$
$10^{3}$ | 250 ns | $1.35\times 10^{-6}$ | $1.71\times 10^{3}$
$10^{3}$ | 1.25 µs | $9.78\times 10^{-7}$ | $1.37\times 10^{3}$
$10^{4}$ | 10 ns | $2.44\times 10^{-6}$ | $3.29\times 10^{3}$
$10^{4}$ | 50 ns | $1.37\times 10^{-6}$ | $2.90\times 10^{3}$
$10^{4}$ | 250 ns | $7.26\times 10^{-7}$ | $2.48\times 10^{3}$
$10^{4}$ | 1.25 µs | $8.75\times 10^{-7}$ | $2.20\times 10^{3}$
Table 2: Values obtained from the simulation of Fig. 7, the case of the square pulse inputs with 100 ps time step. $\beta_{\mathrm{di}}/2\pi$ | $\tau_{\mathrm{di}}$ | $\chi^{2}$ | $t_{\mathrm{ode}}/t_{\mathrm{soen}}$
---|---|---|---
$10^{2}$ | 10 ns | $3.32\times 10^{-4}$ | $1.66\times 10^{3}$
$10^{2}$ | 50 ns | $1.86\times 10^{-4}$ | $1.48\times 10^{3}$
$10^{2}$ | 250 ns | $1.17\times 10^{-4}$ | $1.69\times 10^{3}$
$10^{2}$ | 1.25 µs | $7.87\times 10^{-5}$ | $1.49\times 10^{3}$
$10^{3}$ | 10 ns | $6.53\times 10^{-5}$ | $2.31\times 10^{3}$
$10^{3}$ | 50 ns | $2.20\times 10^{-5}$ | $2.00\times 10^{3}$
$10^{3}$ | 250 ns | $1.74\times 10^{-5}$ | $1.53\times 10^{3}$
$10^{3}$ | 1.25 µs | $1.82\times 10^{-5}$ | $1.62\times 10^{3}$
$10^{4}$ | 10 ns | $4.06\times 10^{-5}$ | $2.77\times 10^{3}$
$10^{4}$ | 50 ns | $1.12\times 10^{-5}$ | $2.47\times 10^{3}$
$10^{4}$ | 250 ns | $1.04\times 10^{-5}$ | $2.11\times 10^{3}$
$10^{4}$ | 1.25 µs | $9.18\times 10^{-7}$ | $1.87\times 10^{3}$
Table 3: Values obtained from the simulation of Fig. 8, the case of square pulse inputs with 1 ns time step. $\beta_{\mathrm{di}}/2\pi$ | $\tau_{\mathrm{di}}$ | $\chi^{2}$ | $t_{\mathrm{ode}}/t_{\mathrm{soen}}$
---|---|---|---
$10^{2}$ | 10 ns | $3.60\times 10^{-2}$ | $5.09\times 10^{3}$
$10^{2}$ | 50 ns | $2.34\times 10^{-2}$ | $5.43\times 10^{4}$
$10^{2}$ | 250 ns | $2.70\times 10^{-2}$ | $1.41\times 10^{4}$
$10^{2}$ | 1.25 µs | $2.45\times 10^{-2}$ | $1.48\times 10^{4}$
$10^{3}$ | 10 ns | $2.84\times 10^{-3}$ | $1.75\times 10^{4}$
$10^{3}$ | 50 ns | $9.25\times 10^{-4}$ | $1.41\times 10^{4}$
$10^{3}$ | 250 ns | $2.50\times 10^{-4}$ | $>4\times 10^{4}$
$10^{3}$ | 1.25 µs | $1.20\times 10^{-4}$ | $1.37\times 10^{4}$
$10^{4}$ | 10 ns | $2.32\times 10^{-3}$ | $2.54\times 10^{4}$
$10^{4}$ | 50 ns | $6.62\times 10^{-4}$ | $2.82\times 10^{4}$
$10^{4}$ | 250 ns | $3.07\times 10^{-4}$ | $2.14\times 10^{4}$
$10^{4}$ | 1.25 µs | $1.72\times 10^{-4}$ | $1.92\times 10^{4}$

##  Appendix E Rate Equations for the Transmitter Light Source

To accurately model a superconducting optoelectronic loop neuron, the
formalism of the computational circuitry—the dendrites—must be accompanied by
a treatment of the transmitter circuit and light source that produce photons
when a soma reaches threshold. We begin by treating the light emitters
themselves in a rate-equation framework and subsequently consider the
transmitter circuit that will provide current to the emission medium. We
consider two types of light emitters: quantum dots and silicon emissive
centers. The quantum dots we have in mind are group III-V self-assembled dots
grown by molecular beam epitaxy or metallorganic chemical vapor deposition.
Such emitters have the advantages of high efficiency and fast radiative
lifetime. However, for the application at hand scalable manufacturing and
integration with semiconductor and superconductor electronics are crucial
considerations, and III-V light emitters have a long history of difficulty
integrating with silicon. The quantum dot model is also applicable to group
IV, Ge-based quantum dots [85], which may prove adequate for the present
application. We additionally consider a rate-equation model of silicon
emissive centers, as they are easy to fabricate and integrate with
electronics, although to date their efficiency is low.

###  E.1 Quantum Dots

The areal density of quantum dots is $n_{\mathrm{D}}$ (600/µm2). We define the
concentration (number per area which can be converted to volume) in the
uncharged, unpopulated, ground state as $n_{0}$. The concentration that have
trapped a hole is $n_{1}^{\mathrm{p}}$, while the concentration that have a
trapped electron is $n_{1}^{\mathrm{n}}$. The concentration that have trapped
a hole and an electron is $n_{2}$. Assuming quantum dots are not created or
destroyed during operation we have

$n_{\mathrm{D}}=n_{0}+n_{1}^{\mathrm{p}}+n_{1}^{\mathrm{n}}+n_{2}.$ (34)

The population of quantum dots in these states is modeled by the following
rate equations:

$\begin{split}\frac{d\,n_{0}}{dt}&=-n_{0}\,p\,c_{01}^{\mathrm{p}}-n_{0}\,n\,c_{01}^{\mathrm{n}}+n_{2}\,e_{20},\\\\[10.0pt]
\frac{d\,n_{1}^{\mathrm{n}}}{dt}&=n_{0}\,n\,c_{01}^{\mathrm{n}}-n_{1}^{\mathrm{n}}\,p\,c_{12}^{\mathrm{p}},\\\\[10.0pt]
\frac{d\,n_{2}}{dt}&=n_{1}^{\mathrm{n}}\,p\,c_{12}^{\mathrm{p}}+n_{1}^{\mathrm{p}}\,n\,c_{12}^{\mathrm{n}}-n_{2}\,e_{20}.\end{split}$
(35)

In Eqs. 35, $c_{01}^{\mathrm{p}}$ is the capture coefficient for an
unpopulated dot to obtain a hole, and $c_{01}^{\mathrm{n}}$ is the capture
coefficient for an unpopulated dot to obtain an electron. The quantities
$c_{12}^{\mathrm{p}}$ and $c_{12}^{\mathrm{n}}$ are defined analogously, and
$e_{20}$ is the emission probability, related to the radiative lifetime by
$\tau_{\mathrm{rad}}=1/e_{20}$. The populations of free electrons and holes in
the semiconductor medium are given by

$\begin{split}\frac{dn}{dt}&=\frac{\eta_{\mathrm{inj}}\,I}{q\,V_{\mathrm{LED}}}-n_{0}\,n\,c_{01}^{\mathrm{n}}-n_{1}^{\mathrm{p}}\,n\,c_{12}^{\mathrm{n}}\\\
&-\frac{v_{\mathrm{s}}\,A_{\mathrm{LED}}}{V_{\mathrm{LED}}}\,n-(C_{n}\,n^{2}\,p+C_{p}\,p^{2}\,n),\\\\[10.0pt]
\frac{dp}{dt}&=\frac{\eta_{\mathrm{inj}}\,I}{q\,V_{\mathrm{LED}}}-n_{0}\,p\,c_{01}^{\mathrm{p}}-n_{1}^{\mathrm{n}}\,p\,c_{12}^{\mathrm{p}}\\\
&-\frac{v_{\mathrm{s}}\,A_{\mathrm{LED}}}{V_{\mathrm{LED}}}p-(C_{n}\,n^{2}\,p+C_{p}\,p^{2}\,n).\end{split}$
(36)

In Eqs. 36 $\eta_{\mathrm{inj}}$ is the injection efficiency, which could be
determined by a spatial carrier transport model, $v_{\mathrm{s}}$ is the
surface recombination velocity, $A_{\mathrm{LED}}$ is the surface area of the
LED, $V_{\mathrm{LED}}$ is its volume, and $C_{n}$ and $C_{p}$ are the Auger
recombination coefficients. Throughout Eqs. 35, $n_{1}^{\mathrm{p}}$ can be
eliminated with
$n_{1}^{\mathrm{p}}=n_{\mathrm{qd}}-n_{0}-n_{1}^{\mathrm{n}}-n_{2}$.

The rate of photon production from electroluminescence under current injection
is given by

$\frac{dN_{\mathrm{ph}}}{dt}=V_{\mathrm{LED}}\,e_{20}\,n_{2}(t).$ (37)

The total number of photons emitted is obtained by integrating:

$N_{\mathrm{ph}}=V_{\mathrm{LED}}\,e_{20}\,\int_{t_{0}}^{t_{f}}n_{2}(t)dt.$
(38)

The primary quantity of interest for our application is the efficiency of
light production. In the context of the present model, we can define the
efficiency as the number of photons generated divided by the number of
electron-hole pairs injected into the structure:

$\eta_{i}=\frac{N_{\mathrm{ph}}}{N_{eh}},$ (39)

where $N_{eh}$ is the total number of electron-hole pairs injected into the
intrinsic region of the diode.

To simulate the QDs we use a radiative lifetime of 1 ns, giving $e_{20}=1$
GHz. The capture coefficients are expected to be quite fast for quantum dots,
and we take $c_{01}=c_{12}=10^{-10}$m3/s. The precise values of $c_{01}$ and
$c_{12}$ have only minor impact on the numerical results as the capture rates
are always much smaller than the recombination times. For the parameters
governing non-radiative recombination we $v_{\mathrm{s}}=2.5\times 10^{3}$cm/s
[86] and $C_{n}=C_{p}=7\times 10^{-30}$cm6/s [87].

The system of ODEs given by Eqs. 35 and 36 was again solved using solve_ivp.
Complete analysis of these light sources is beyond the scope of the present
work. The important take-away is that the efficiency of the source can be
quite high, provided the number of carriers injected is commensurate with the
number of dots in the diode. For this spiking neuromorphic application, we are
interested in cases wherein an initial pulse of electrons and holes is
injected into the diode and left to decay. If the pulse of injected carriers
contains appreciably more electrons and holes than the number of quantum dots,
the excess carriers recombine non-radiatively through surfaces or Auger before
the populated dots have time to emit photons and trap additional carriers.
Thus, for a neuron to efficiently produce light upon reaching threshold, the
light source must have a number of quantum dots that is chosen in accordance
with the number of synapses that are intended to receive photons, and the
transmitter circuit that injects the light source must be designed to inject
the appropriate number of electrons and holes in a brief burst of charge. Such
a transmitter circuit is summarized in Appendix F. Time-domain simulations of
quantum dots being driven by such a driver circuit will be presented in that
Appendix. First, we briefly introduce the rate equation model of silicon
emissive centers.

###  E.2 Silicon emissive centers

We can expand the quantum dot model to treat silicon emissive centers. A
similar model was presented in Ref. 88, but in that work steady state emission
was investigated to seek a continuous-wave laser. He we require the full
dynamical equations to investigate transient behavior. The total concentration
of the particular emissive center of interest is $n_{\mathrm{D}}$. We define
the concentration (number per volume) in the uncharged, unpopulated, ground
state as $n_{0}$. The concentration that have trapped a hole is $n_{1}$. The
concentration that have trapped a hole and an electron is $n_{2}$. Following
Ref. 88 we focus on W centers [12] and assume a hole is always trapped before
an electron. Assuming emissive centers are not created or destroyed during
operation we have

$n_{\mathrm{D}}=n_{0}+n_{1}+n_{2}.$ (40)

The equations governing the population of emissive centers are

$\begin{split}\frac{dn_{0}}{dt}&=-n_{0}\,p\,c_{01}+(n_{\mathrm{D}}-n_{0}-n_{2})c_{01}k_{01}+n_{2}\,e_{20}\\\
\frac{dn_{2}}{dt}&=-n_{2}\,(e_{20}+c_{12}\,k_{12})+(n_{\mathrm{D}}-n_{0}-n_{2})\,n\,c_{12}.\end{split}$
(41)

In Eqs. 41, $c_{01}$ is the capture coefficient for an unpopulated center to
obtain a hole. As was argued in Ref. 88, W centers appear to always trap a
hole before an electron. The capture coefficient for trapping an electron
after the center has been populated with a hole is $c_{12}$. $e_{20}$ is the
emission coefficient. The emission coefficients are related to the capture
coefficients through the relations

$\begin{split}e_{10}&=c_{01}\,k_{01}=c_{01}\frac{n_{0}^{*}n_{c}^{*}}{n_{1}^{*}}\\\
e_{21}&=c_{12}\,k_{12}=c_{12}\frac{n_{1}^{2}n_{c}^{*}}{n_{2}^{*}}.\end{split}$
(42)

Asterisks refer to equilibrium values. The relations can be further specified
by

$\begin{split}k_{01}&=\frac{n_{0}^{*}n_{c}^{*}}{n_{1}^{*}}=\frac{N_{v}}{g_{n}}\mathrm{e}^{-E_{h}/kT}\\\
k_{12}&=\frac{n_{1}^{*}n_{c}^{*}}{n_{2}^{*}}=\frac{N_{c}}{g_{e}}\mathrm{e}^{(E_{e}-E_{g})/kT}.\end{split}$
(43)

Here $N_{v}$ and $N_{c}$ are the valence and conduction band densities of
states, $g_{n}$ and $g_{e}$ are the degeneracies of the states of the emissive
center, $E_{g}$ is the band gap, and $E_{h}$ and $E_{g}-E_{e}$ are the hole
and electron binding energies. The energy of the emitted photon is
$E_{e}-E_{h}$. For the case of the W center, this equals 1.018 eV. The top of
the valence band is set to zero energy. See Ref. 89, Ch. 5 for further
explanation.

In the second of Eqs. 41 the term $e_{20}+c_{12}\,k_{12}$ plays a prominent
role in determining the efficiency of emission. This term represents
competition between emission ($e_{20}$) and release of the electron in the
bound exciton back to the conduction band ($c_{12}\,k_{12}$). The ratio
$e_{20}/(c_{12}\,k_{12})$ is roughly $10^{5}$ at 4.2 K, so emission is much
more likely than exciton dissociation at this temperature.

In addition to the rate equations for the emissive center populations, the
model for Si light sources should include the effects of a population of non-
radiative centers that will inevitably be produced when the implants for the
emissive centers are performed. Currently knowledge of the number of different
types of non-radiative centers, their capture coefficients, and their
lifetimes is limited. For simplicity, we assume there is one type of dominant
non-radiative recombination center, we assume it can be populated first by an
electron or hole, we assume the capture cross section is $c_{01}$, and we
assume the secondary capture of a hole or electron is governed by $c_{12}$. We
also assume after capture of an electron, hole, or both, no carriers are ever
re-emitted to the bands, which is reasonable at the low temperatures necessary
for the superconducting electronics present in the system. The total
concentration of non-radiative defects is
$n_{\mathrm{D}}^{\mathrm{nr}}=n_{0}^{\mathrm{nr}}+n_{1}^{\mathrm{nr}}+n_{2}^{\mathrm{nr}}$.
The population rate equations for non-radiative centers then can be written
as:

$\begin{split}\frac{dn_{0}^{\mathrm{nr}}}{dt}&=-n_{0}^{\mathrm{nr}}\,c_{01}(p+n)+n_{2}^{\mathrm{nr}}\,e_{20}^{\mathrm{nr}}\\\
\frac{dn_{2}^{\mathrm{nr}}}{dt}&=-n_{2}^{\mathrm{nr}}\,e_{20}^{\mathrm{nr}}+c_{12}\,(n_{\mathrm{D}}^{\mathrm{nr}}-n_{0}^{\mathrm{nr}}-n_{2}^{\mathrm{nr}})\,(p+n).\end{split}$
(44)

The coefficients $c_{ij}$ are different here than in the emissive center
equations, but we expect them to be fast, so this is of little consequence.
The lifetime of the non-radiative transition of the center (which may be a
trap-assisted Auger process or, more likely, phonon recombination) is given by
$\tau_{\mathrm{nr}}=1/e_{20}^{\mathrm{nr}}$.

Equations 41 and 44 are coupled to the charge-carrier rate equations. For the
electron concentration we have

$\begin{split}\frac{dn}{dt}&=\frac{\eta_{\mathrm{inj}}\,I}{q\,V_{\mathrm{LED}}}-(n_{\mathrm{D}}-n_{0}-n_{2})\,n\,c_{12}+n_{2}\,c_{12}\,k_{12}\\\
&-n_{0}^{\mathrm{nr}}\,n\,c_{01}-(n_{\mathrm{D}}^{\mathrm{nr}}-n_{0}^{\mathrm{nr}}-n_{2}^{\mathrm{nr}})\,n\,c_{12}\\\
&-\frac{v_{\mathrm{s}}\,A_{\mathrm{LED}}}{V_{\mathrm{LED}}}\,n-(C_{n}\,n^{2}\,p+C_{p}\,p^{2}\,n),\end{split}$
(45)

and for the hole concentration we have

$\begin{split}\frac{dp}{dt}&=\frac{\eta_{\mathrm{inj}}\,I}{q\,V_{\mathrm{LED}}}-n_{0}\,p\,c_{01}\\\
&-n_{0}^{\mathrm{nr}}\,p\,c_{01}-(n_{\mathrm{D}}^{\mathrm{nr}}-n_{0}^{\mathrm{nr}}-n_{2}^{\mathrm{nr}})\,p\,c_{12}\\\
&-\frac{v_{\mathrm{s}}\,A_{\mathrm{LED}}}{V_{\mathrm{LED}}}p-(C_{n}\,n^{2}\,p+C_{p}\,p^{2}\,n).\end{split}$
(46)

To numerically solve Eqs. 41, 44, 45, and 46, we need to specify $c_{01}$,
$c_{12}$, and $e_{20}$. For the process of photon emission, related to
$e_{20}$, we take $e_{20}=\tau_{ec}^{-1}=1/40$ ns $=2.5\times 10^{7}$ [12,
90]. To our knowledge, $c_{01}$ and $c_{12}$ have not been measured. We take
$c_{01}=5\times 10^{-14}$m3/s and $c_{12}=1.8\times 10^{-11}$m3/s. These
numbers are motivated by a simple geometrical model as well as a hydrogenic
model of an impurity screened by the dielectric [91]. Numerical studies show
that the results are not sensitive to these values across several orders of
magnitude.

To specify the nonradiative terms in the rate equations, we take the Auger
recombination coefficients as $C_{n}=2.8\times 10^{-31}$ cm6/s and
$C_{p}=1\times 10^{-31}$ cm6/s. These are the values at 300 K [87], and they
will reduce slightly at low temperature (scaling with temperature as $T^{0.6}$
[92]), but this is not appreciable given the low level of accuracy in the
present model, so we ignore the temperature dependence of Auger and use 300 K
values. The surface recombination can be made quite low in properly treated Si
and Si-SiO2 interfaces, with values as low as 0.25 cm/s in the literature
[93]. Throughout these simulations we use $v_{\mathrm{s}}=2.5$cm/s.

Time-domain simulations of silicon emissive centers based on Eqs. 41, 44, 45,
and 46 provides a similar message to that learned for quantum dots: the
emitters will populate rapidly, but if the level of carrier injection is too
large, excess carriers will go to waste. In the case of silicon emissive
centers, loss is not due to surfaces or Auger, but rather to nonradiative
centers. The relative concentration of nonradiative centers formed during the
ion implantation used to realize W centers is likely to be a crucial quantity
determining the ultimate efficiency of such a light source. Given the
competing populations of emissive and nonradiative centers, it is not
sufficient to inject the same number of carriers as emissive centers and
desired photons; one must inject enough carriers to populate both the
nonradiative centers and the emissive centers. However, if the concentration
of nonradiative centers can be limited to something like 10 times the
concentration of emitters, total light-production efficiency may stay within a
useful range.

We next describe transmitter circuits capable of providing brief pulses of
injected current to produce the desired number of photons necessary to serve a
neuron’s downstream synaptic connections.

##  Appendix F Circuit Equations for the Transmitter Driver Circuit

Figure 28: Diagram of the transmitter circuit with relevant components
labeled. (a) The neuronal integration loop ($\mathsf{NI}$) and tron
thresholding element. (b) The superconductor-semiconductor interface where the
tron switches the first inverter. (c) The inverter pair that produces the
digital drive signal. (d) The LED driver MOSFET and the LED itself.

The transmitter circuit under consideration is shown in Fig. 28. The circuit
combines JJs, an hTron thresholding element [14], MOSFETs (PMOS and NMOS), and
an LED. The hTron serves as the interface between the superconducting and
semiconducting domains, enabling low-voltage JJs to switch a CMOS inverter.
The basic operation of the circuit is as follows. Signal moves from left to
right. When the JJ in the neuronal integration (NI) loop has added sufficient
current to drive the hTron gate above threshold, the hTron channel switches
from zero resistance to high resistance. The current bias to the hTron channel
($I_{\mathrm{b}}$) is thus shunted from the hTron channel
($I_{\mathrm{tron}}$) to the passive reset branch ($L_{\mathrm{r}}$,
$r_{\mathrm{r}}$), temporarily inducing a voltage at node $V_{1}$ sufficient
to switch the first CMOS inverter ($M_{1}$, $M_{2}$). The two inverters serve
to produce a digital signal with voltage $V_{\mathrm{dd}}$ applied to MOSFET
$M_{\mathrm{LED}}$ when the neuronal threshold is reached and the hTron
switches. This digital behavior ensures the superconducting components are
well separated from the LED, providing consistent current biasing of the LED
decoupled from the performance of the thresholding element. The LED driver
MOSFET ($M_{\mathrm{LED}}$) can then be designed with a width-to-length ratio
and doping level commensurate with the number of photons desired from the LED,
which is chosen based on the number of synaptic connections made by the
neuron. $M_{\mathrm{LED}}$ will deliver a pulse of current to the LED while
the hTron is keeping the voltage $V_{1}$ high, and when that voltage drops as
the passive reset circuit performs its operation, the current to the LED will
cease. As is shown below, the MOSFETs are active for a few nanoseconds. With
this qualitative description of the circuit operation in mind, we now present
the circuit model used in this work to provide a phenomenological treatment of
the behavior.

To treat the thresholding component where the current added to the integration
loop drives an hTron above threshold, we model the hTron as a device that
switches from a zero resistance to a high-resistance state when a certain
current is reached. Accurately modeling the electrothermal dynamics of the
hTron is difficult and will do little to improve the model for the present
purpose, so a simple resistive switch is employed.

To simulate the behavior of the MOSFETs, we use a charge-control model [94]
with current-voltage characteristics given by

$\begin{split}I_{\mathrm{ds}}(V_{\mathrm{ds}},V_{\mathrm{gs}})&=\frac{W\mu_{n}c_{\mathrm{i}}}{L}\\\
&\times\begin{cases}(V_{\mathrm{gt}}-V_{\mathrm{ds}}/2)V_{\mathrm{ds}},&\text{for
}V_{\mathrm{ds}}\leq V_{\mathrm{sat}}\\\ V_{\mathrm{gt}}^{2}/2,&\text{for
}V_{\mathrm{ds}}>V_{\mathrm{sat}}.\end{cases}\end{split}$ (47)

$V_{\mathrm{gt}}$ is the voltage above threshold:
$V_{\mathrm{gt}}=V_{\mathrm{gs}}-V_{\mathrm{t}}$. We assume the drain-source
current is zero for $V_{\mathrm{gs}}<V_{\mathrm{t}}$, i.e., we ignore
subthreshold behavior. In this model, $V_{\mathrm{sat}}=V_{\mathrm{gt}}$.

Figure 29: Time traces from the transmitter circuit. (a) The resistance of the
tron as a function of time. (b) The voltages of the MOSFETs. (c) The currents
in the tron and related circuit elements. (d) Currents through the LED. (e)
Temporal zoom of the MOSFET voltages. (f) Temporal zoom of the LED currents.
(g) Currents into the LED and capacitor when different numbers of quantum dots
are present in the light source.

We treat the LED with a conventional diode model of the form

$\begin{split}I_{\mathrm{LED}}=&\,eA\left[(D_{p}/L_{p})p_{n}+(D_{n}/L_{n})n_{p}\right]\\\
&\times\left[\mathrm{exp}\left(eV/k_{\mathrm{B}}T\right)-1\right].\end{split}$
(48)

In the present work we have considered two classes of LEDs. One is based on
silicon emissive centers [12], as have been considered in the context of SOENs
due to their unique process compatability and feasibility for enabling low-
cost manufacturing of large systems. The other is based on III-V quantum dots,
specifically with material parameters for InGaAs quantum dots on a GaAs
platform. For integration with SOENs, wafer bonding of such substrates may be
employed. The model presented here can take into account either of these light
sources. Both will obey an equation of the form given by Eq. 48, with the
diode placed in parallel with a capacitor and a resistor as shown in Fig. 28.
In reality the capacitance and parasitic shunt resistance will both be
functions of applied voltage, but for simplicity we fix these values. Based on
the source rate equation model in Appendix E we know that we would like to
have one emitter (quantum dot or emissive center) per photon to be generated
by the light source per neuronal firing pulse. We assume the emitters will
have a specified density $\rho_{\mathrm{ec}}$, so the value of the capacitance
will depend on the number of emitters through the geometry of the junction.
The number of emitters is in turn specified by the number of synaptic
connections made by the neuron. Therefore, we consider two values of
capacitance per unit area ($C_{\mathrm{a}}=1\times 10^{-7}$F/µm2 and
$C_{\mathrm{a}}=1\times 10^{-5}$F/µm2) to cover a range of values that may be
found in various devices. The value of the capacitor, $C_{\mathrm{LED}}$, is
determined by
$C_{\mathrm{LED}}=C_{\mathrm{a}}\,N_{\mathrm{ph}}/\rho_{\mathrm{ec}}$ with
$\rho_{\mathrm{ec}}$ specified as a number of emitters per unit area.

Figure 30: Time traces for the light source model. (a)-(d) The quantum-dot
model. (e)-(i) The emissive center model. (a) The current pulse injected from
the transmitter circuit. (b) The populations of quantum dots in the ground and
excited states normalized to the concentration of quantum dots in the LED
volume. (c) The concentrations of free electrons and holes normalized to the
total number injected into the LED volume. (d) The simulated
electroluminescence signal, the cumulative distribution, and a histogram of
sampled emission times. (e) The current pulse injected from the transmitter
circuit. (f) The populations of emissive centers in the ground and excited
states normalized to the concentration of emitters in the LED. (g) The
populations of non-radiative recombination centers in the ground and excited
states normalized to the concentration of non-radiative recombination centers
in the LED. (h) The concentrations of free electrons and holes. (i) The
simulated electroluminescence signal, the cumulative distribution, and a
histogram of sampled emission times.

Given these models for the hTron, MOSFETs, and LED the following circuit
equations can be derived:

$\begin{split}\frac{d\,V_{1}}{dt}&=\frac{I_{b}-I_{1}-I_{3}}{C_{\mathrm{I}}},\\\
\frac{d\,V_{2}}{dt}&=\frac{-1}{C_{\mathrm{I}}}\left[I_{\mathrm{p}}(V_{2}-V_{\mathrm{dd}},V_{1}-V_{\mathrm{dd}})+I_{\mathrm{n}}(V_{2},V_{1})\right],\\\
\frac{d\,V_{3}}{dt}&=\frac{-1}{C_{\mathrm{D}}}\left[I_{\mathrm{p}}(V_{3}-V_{\mathrm{dd}},V_{2}-V_{\mathrm{dd}})+I_{\mathrm{n}}(V_{3},V_{2})\right],\\\
\frac{d\,V_{4}}{dt}&=\frac{1}{C_{\mathrm{LED}}}\left[I_{\mathrm{n}}(V_{\mathrm{dd}}-V_{4},V_{3})-I_{\mathrm{LED}}(V_{4})-\frac{V_{4}}{r_{\mathrm{LED}}}\right],\\\
\frac{d^{2}\,I_{1}}{dt^{2}}&=\frac{I_{b}-I_{1}-I_{3}}{L_{\mathrm{t}}C_{\mathrm{I}}}-\frac{\dot{r}_{\mathrm{t}}}{L_{\mathrm{t}}}\,I_{1}-\frac{r_{\mathrm{t}}I_{1}}{L_{\mathrm{t}}}\,\frac{d\,I_{1}}{dt},\\\
\frac{d\,I_{3}}{dt}&=\frac{r_{\mathrm{t}}}{L_{\mathrm{r}}}\,I_{1}+\frac{L_{\mathrm{t}}}{L_{\mathrm{r}}}\,\frac{d\,I_{1}}{dt}-\frac{r_{\mathrm{r}}}{L_{\mathrm{r}}}\,I_{3}.\end{split}$
(49)

Here, $I_{\mathrm{p}}(V_{\mathrm{ds}},V_{\mathrm{gs}})$ and
$I_{\mathrm{n}}(V_{\mathrm{ds}},V_{\mathrm{gs}})$ are the PMOS and NMOS forms
of Eq. 47 (i.e., current flows from source to drain in PMOS and drain to
source in NMOS), $C_{\mathrm{I}}$ is the capacitance of an inverter, and
$C_{\mathrm{D}}$ is the capacitance of the MOSFET driver to the LED. The time-
dependent resistance of the hTron is $r_{\mathrm{t}}$ and its temporal
derivative is $\dot{r}_{\mathrm{t}}$.

Figure 31: All 30 inputs to the nine-pixel classifier.

Time domain circuit simulations of the transmitter model given by Eqs. 49 have
been carried out using solve_ivp. A sample of results is given in Fig. 29. The
resistance pulse of the hTron is shown in Fig. 29(a), and the relevant
voltages are shown in Fig. 29(b). The currents are shown in Fig. 29(c) and
(d). A temporal zoom of the voltage traces is shown in Fig. 29(e), where the
MOSFET threshold is shown by the grey dotted line. The gates of the first
inverter stage quickly rise above threshold, and a nearly ideal square pulse
is delivered to the gate of $M_{\mathrm{LED}}$, the driver MOSFET to the light
source, as can be seen in the temporal zoom of the LED currents of Fig. 29(f).
It is for this reason that the two-inverter digital circuit is employed. After
the voltage to the first inverter drops below threshold, $M_{\mathrm{LED}}$
quickly shuts off, and the voltage across the LED block decays with the RC
time constant of that circuit. The resistance, $r_{\mathrm{LED}}$, has been
chosen to establish a 100 ns time constant, and the actual value in a
fabricated circuit will need to be investigated.

In the temporal zooms of Fig. 29(e) and (f), the time at which the hTron
switches to the resistive state is labeled (15 ns), and the total delay of the