Dynamic Processes in a Superconducting Adiabatic Neuron with Non-Shunted Josephson Contacts

Abstract

We investigated the dynamic processes in a superconducting neuron based on Josephson contacts without resistive shunting (S_C-neuron). Such a cell is a key element of perceptron-type neural networks that operate in both classical and quantum modes. The analysis of the obtained results allowed us to find the mode when the transfer characteristic of the element implements the “sigmoid” activation function. The numerical approach to the analysis of the equations of motion and the Monte Carlo method revealed the influence of inertia (capacitances), dissipation, and temperature on the dynamic characteristics of the neuron.

1. Introduction

Theoretical and practical research of artificial intelligence systems, methods of machine learning, and artificial neural networks are being actively developed in recent years [1,2]. The main goal of artificial neural networks is to effectively use the features of the human brain such as the learning ability, the ability to model complex separating surfaces in a multidimensional feature space, distributed memory in order to solve computational problems.

The search for the optimal element base and architecture of artificial neural networks continues [3,4,5,6]. In particular, computing accelerators—hardware implementations of deep neural networks—show high efficiency. The key parameters of such devices are performance and energy efficiency. As a commercially available example, we can mention Google’s tensor processor-based supercomputers (TPU) provided in the cloud, whose performance in 2021 exceeded exaflops (10¹⁸ floating-point operations per second).

The first implementations of artificial neural networks were based on semiconductor CMOS technology. Nevertheless, despite a number of advantages of this technique [4], this technology has a number of disadvantages for a large number of neurons and synapses, the main of which are high power consumption, limited clock frequencies, and the complexity of the aforementioned base cells. The limitations imposed by the dissipative nature of CMOS technology make the hardware acceleration of both large convolutional neural networks of deep learning and the increasingly popular transformer networks extremely difficult to implement [5]. On the other hand, an alternative and actively developing approach for the implementation of artificial neural networks (ANN) is the development of energy-efficient superconducting solutions [1,7,8,9,10,11].

Initially, the concept of a superconducting neural cell for spike neural networks [7,8] was based on two Josephson transmission lines (JTL) connected via resistors, and where the magnitude of the conductivity was regarded as a synaptic weight. In this case, pulses of current and voltage in JTL, which are related to the motion of specific quasiparticles—fluxons, are considered as neuroimpulses. It worth mentioning that these were ones of the first experimental results [9] on the study of operation of an integrated superconducting neural network. We can emphasize several main branches of the further cell design evolving [10,11,12,13,14,15] such as using the Josephson comparator as a neuron with a pseudosigmoidal activation function, replacing resistors in synaptic (connecting) elements with tunable Josephson junctions, and implementation of the principles of adiabatic superconducting logic (ASL) [16,17,18]. The latter helped to overcome the fundamental limitation on energy efficiency [14]. When using ASL schemes, the functioning of the ANN is based on handling the information given not in the form of the presence or absence of a quant, but as the magnitude and direction of circulating superconducting currents. The performance of the hardware implementation of a superconducting neural processor while executing the test on examples of standard configurations of neural networks exceeds the performance of a semiconductor analog (TPU) by 23 times on average [1]. These indicators were demonstrated when using memory with a bandwidth (300 GB/s) and a typical clock frequency of a superconducting processor (52.6 GHz) [1]. A number of examples per learning iteration was set according to the capacity of the on-chip memory, which made it possible to avoid addressing the external memory inside the loop. At the same time, the energy consumption is 490 times less, so even with taking into account necessity of cooling from the room temperature level assuming the required additional power at the rate of 400 W per W, the superconducting system is more energy efficient.

Unfortunately, the existing prototypes of adiabatic superconducting network cells are quite large in size [19,20,21]. By eliminating the resistive shunting of Josephson junctions from the proposed schemes, we will significantly reduce their size and energy dissipation per operation. At the same time, the influence of “inertia” (associated with the capacity of contacts which became significant) on all dynamic processes occurring in the neuron will become critically important.

In this work, we have developed the design of the basic cell of the neural network named the S_C-neuron, which is based on Josephson junctions without resistive shunting.

On the basis of the Hamiltonian formalism and using the Monte Carlo method, the analysis of dynamic processes affecting the neuron’s activation function calculated in one clock cycle is carried out, and the nature of switching is studied taking into account the inertia (the capacitance in the circuit), dissipation, and the temperature. We carried out the selection and optimization of the parameters of a perceptron-type superconducting circuit for the implementation of the “sigmoid” activation function, which is most convenient for neural network training algorithms and solving problems of pattern recognition and images.

2. Neuron Model and Research Methods

The most important characteristic of the cells of superconducting adiabatic ANNs is the flux-to-current transformation. This is the dependence of the output current, $I_{out}$, versus input signal (applied magnetic flux, ${\Phi }_{in}$). In Figure 1a, we schematically showed the corresponding normalized values. Earlier in [18,19,20,21,22], it was demonstrated how it is necessary to change the quantron to obtain a computationally effective transfer characteristic (activation function). In this case, the quantron is actually a single-contact Josephson superconducting interferometer for which it is assumed that the critical current of its Josephson contact can change under the influence of an external signal. The inertia-less mode of operation of this device was studied.

To study the flux-to-current transformation, we will consider the S_C-neuron shown in Figure 1a. The idea of the proposed circuit is to combine the nonlinearity of the Josephson junction (quantum part) with the linear flux-to-current transformation of a conventional inductance in order to obtain a transfer characteristic close to the mathematical sigmoidal function.

We assume that during operation, a magnetic flux pulse, ${\phi }_{in}(t)=\frac{2\pi {\Phi }_{in}(t)}{{\Phi }_0}$, is supplied to the input of the Sc neuron, which we normalize to the magnetic flux quantum ${\Phi }_0$:

$${\phi }_{in}(t)=A_{in}\left(\frac{1}{1+\exp [-2D(t-t_1)]}+\frac{1}{1+\exp [2D(t-t_2)]}\right)-A_{in}.$$

(1)

In this case, the amplitude, A_in, and the duration of the action, $\tau =t_2-t_1$ ($t_2>t_1$), are selected on the basis of considerations for the system to go through the entire cycle of evolution (activation function). In addition, the rate of rise/fall of the control signal is determined by the parameter D in Expression (1), see Figure 1b.

The flux-to-current transformation of the system in Figure 1a can be represented as a joint solution of the Kirchhoff equations and the phase balance condition in superconducting circuits:

$$\left\{\begin{array}{l}\phi +l\cdot i=\raisebox{1ex}{${\phi }_{in}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.+l_{out}{i}_{out}\\ \phi +l\cdot i={\phi }_{in}+l_a{i}_a\\ i+i_{out}+i_a=0\end{array}\right.,$$

(2)

where $\phi$ is the phase drop across the Josephson contact, the output current $i_{out}=\frac{I_{out}}{I_c}$ is normalized to the critical current $I_c$, and inductances, $l,\hspace{0.33em}{l}_a,l_{out}$, are normalized to $\frac{2\pi I_c}{{\Phi }_0}$, as shown in Figure 1a.

From the first equation in system (2), it is easy to obtain the dependence of the output current $i_{out}$ on the input flux ${\phi }_{in}$, defined by Expression (1), and the current $i$ at the Josephson junction:

$$i_{out}=\frac{{\phi }_{in}-2l_{a}i}{2(l_a+l_{out})},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}i=b{\phi }_{in}-a\phi ,$$

(3)

where the parameters are defined as $a=\frac{l_a+l_{out}}{l\cdot l_a+l_{out}\cdot \left(l+l_a\right)}$ and $b=\frac{l_a+2l_{out}}{2\left(l\cdot l_a+l_{out}\cdot \left(l+l_a\right)\right)}$.

In this case, within the framework of the “resistive—capacitance model” for a Josephson junction (RCSJ-model), by introducing the Stewart-McCumber parameter $\beta =\frac{2\pi R^2{I}_{c}C}{{\Phi }_0}$ and the characteristic Josephson junction frequency ${\omega }_c=\frac{2eR{I}_c}{\hslash }$ (here R and C are the characteristic resistance and capacitance of the junction), we can write the equation that determines the dynamics of the phase at a given current:

$$\frac{\beta }{{\omega }_c^2}\ddot{\phi }+\frac{1}{{\omega }_c}\dot{\phi }+\sin \phi =\frac{{\phi }_{in}\cdot \left(l_a+2l_{out}\right)-2\phi \cdot \left(l_a+l_{out}\right)}{2\left(l\cdot l_a+l_{out}\cdot \left(l+l_a\right)\right)}.$$

(4)

In the case of contact with low losses (taking into account only conservative forces), we have:

$$\frac{1}{{\omega }_c}=\frac{\hslash }{2eR{I}_c}\to 0, \frac{\beta }{{\omega }_c^2}=\frac{\hslash C}{2e{I}_c}=\frac{1}{{\omega }_p^2},$$

where ω_p is the plasma frequency of the Josephson junction. Therefore, for a system without dissipation, the following equation can be used:

$$\frac{1}{{\omega }_p^2}\ddot{\phi }+\sin \phi ={\phi }_{in}\cdot b-\phi \cdot a.$$

(5)

In this case, it will be convenient to rescale the time according to $\tilde{t}\to {\omega }_{p}t$. Note that the equation for the inertialess model of a neuron is obtained when $C\to 0$, i.e., for $\frac{1}{{\omega }_p^2}\to 0$.

In the context of the Josephson junction working as a part of a S_C-neuron, the inverse characteristic Josephson frequency, $1/{\omega }_c$, determines the characteristic relaxation time of the neuron, and the ${\omega }_p$ parameter is responsible for the oscillatory mode at low dissipation. This will be clearly seen in the next section.

To describe the S_C-neuron (see Figure 1a), we will use the Hamiltonian formalism. The Equation (5) can be written as a system of two first-order equations:

$$\left\{\begin{array}{l}\dot{\phi }=\frac{\partial H(p,\phi )}{\partial p}=p,\\ \dot{p}=-\frac{\partial H(p,\phi )}{\partial \phi }=-a\hspace{0.33em}\phi +b\hspace{0.33em}{\phi }_{in}-\sin \phi ,\end{array}\right.$$

(6)

where the Hamilton function is introduced:

$$H\left(\phi ,p\right)=\frac{p^2}{2}+U(\phi ,{\phi }_{in}).$$

(7)

Note that the system under consideration becomes similar to the motion of a classical particle with unit mass and momentum $p=\dot{\phi }$ in the potential:

$$U\left(\phi ,{\phi }_{in}\right)=\frac{{\left(b{\phi }_{in}-a\phi \right)}^2}{2a}+1-\cos \phi .$$

(8)

The symmetry of the potential in (7), and, consequently, the dynamic processes in the system, significantly depends on the inductances included in the parameters a and b of the system under consideration. For a shunted Josephson junction [21], the following relation should be satisfied: $l_a=1+l$. Having fixed this condition, we analyzed the changes in the potential during the operation of the S_C-neuron. Note that the input signal smoothly changes the potential of the system between two stable equilibrium positions near ϕ = 0 and φ = 2π. As can be seen from Figure 2, the behavior of the potential significantly depends on the inductance parameter l of the circuit. At the same time, two fundamentally different modes can be distinguished: “neural mode”, when ϕ(ϕ_in) is uniquely defined, and “logical element mode”, when there is an ambiguity area for the dependence of the junction phase on the input flux. For the neural mode, there is a symmetric adiabatic and reversible transition of the system between the initial and final states: $\phi =0\iff 2\pi$. At the same time, the potential energy always has one minimum. In the “logical element mode”, the potential energy becomes double-well (bistability arises, as in the cells of ASL-circuits) near $\phi =\pi$ (see Figure 2b).

The range of parameters for the existence of a logical mode is not difficult to obtain from the condition:

$$\phi +\frac{1}{a}\sin \phi =\frac{b}{a}{\phi }_{in},$$

Whatwhich follows from Equation (5), when $C\to 0$. This equation is transcendental and will have a pair of solutions if the coefficient at sinϕ is greater than one. Therefore, we have the condition

$$\frac{1}{a}=l+\frac{l_a\cdot l_{out}}{l_a+l_{out}}=\left\{l_a=l+1\right\}>1, l^2+2l_{out}\cdot l-1>0.$$

This inequality will hold (the potential will have more than one extremum) and have physical meaning only if $l>\sqrt{l_{out}{}^2+1}-l_{out}$.

Dissipative effects in the system, according to [23], are taken into account by adding non-conservative forces to Equation (6):

$$\left\{\begin{array}{l}\dot{\phi }=p,\\ \dot{p}=-a\hspace{0.33em}\phi +b\hspace{0.33em}{\phi }_{in}-\sin \phi -\gamma p,\end{array}\right..$$

(9)

where $\gamma =\raisebox{1ex}{${\omega }_p$}\!\left/ \!\raisebox{-1ex}{${\omega }_c$}\right.$ is the damping parameter.

Thus, the established relationships between the inductances of the circuit show how to control the potential relief, providing bistable states of the neuron. If the capacitance and resistance can be neglected, then the inertia-less mode of the neuron takes place, when the output current immediately follows the control flux [18,19,20,21,22]. The presence of capacitance will lead to nontrivial dynamics in the phase space, which will also depend on the initial conditions of the phase and its velocity. In this case, the role of the initial temperature-dependent distribution of the phase and its velocity becomes essential: the distribution determines the available regions in the phase space of states and their final distribution that arises after the action of the input flux. In the next section, we will simulate the processes of switching a neuron under the action of these factors.

3. Results of the Analysis of Dynamic Processes

3.1. Excluding Dissipation

First, we investigated numerically dynamic processes in the system described by Equation (9) for the case $T\to 0$. Let at the initial moment of time that the system was localized at the minimum of the potential, which is equivalent to considering the dynamics of the behavior of a particle with the initial conditions: $\phi (t=0)=0;\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}p(t=0)=0$. We have studied the behavior of the flux-to-current transformation and the phase portrait of the system for two typical modes of operation (neuro- and logic-mode).

When changing the parameters of the system (values of inductances), it was found that l is responsible for the amplitude of the current value at the output, $i_{out}$. At the same time, the general sigmoidal form of the transfer function is preserved, if l is greater, then the lower the value of the current at the output. The values of the shoulder inductances and the inductive coupling coefficients affect the “delay” of the supplied flux at the device input. The quantity l_a is responsible for the tilt of the symmetry axis of the sigmoid relative to the horizontal axis [21].

In an adiabatically smooth switching mode (when the particle momentum is less than the rate of change of the input signal, for example, D = 0.008 in Figure 3) for l < 0.5, a complete coincidence of the two evolutionary paths of the system is achieved with an increase in the rising and falling phases: “↑” ($\phi =0\to 2\pi$) and “↓” ($\phi =2\pi \to 0$) of the input signal. This is clearly seen in Figure 3a, when the flux-to-current transformation has a sigmoidal dependence, and the phase trajectories (black curves in Figure 3b) have a symmetrical shape with the rise/fall of the input flux. With an increase in the inductance value $l\to \sqrt{l_{out}{}^2+1}-l_{out}$ oscillations appear on the transfer characteristics at a frequency close to ${\omega }_p=\sqrt{\frac{2e{I}_c}{\hslash C}}$. This occurs due to imparting an additional impulse to the particle (destruction of the adiabatic approximation) upon passing the value ϕ = ϕ_in = π (Figure 2b). Under this condition, there is a sharp increase/decrease in the signal amplitude caused by external influence. This is due to the fact that the system does not have time to relax and cannot reach the minimum of its potential energy during the impact change. It is worth noting that for the neural mode of operation $l<\sqrt{l_{out}{}^2+1}-l_{out}$, the transfer characteristics averaged over plasma oscillations $\phi =0\iff 2\pi$ coincide and have a pronounced sigmoidal character. A different type of dependencies is manifested for the logic element mode, when the flux-to-current dependencies averaged over the plasma frequency have a hysteresis character (solid and dotted black curves in Figure 3c). At the same time, on the phase diagrams (Figure 3d), it can be seen that the particle completes movement near the equilibrium positions: ϕ = 0 (the red trajectory in Figure 3d for the rising phase of the control signal “↑” at $\phi =0\to \phi =2\pi$) and ϕ = 2π (the orange trajectory in Figure 3d) for the falling phase of the control signal “↓” at $\phi =2\pi \to \phi =0$).

For a more detailed study of the dynamic processes associated with inertia in the system, we investigated the behavior of integral characteristics, such as the amplitudes of A/I_c current oscillations at the output and their period $T_p=\frac{2\pi }{{\omega }_p}$ at the different values of inductors and capacitance. The whole process of evolution was divided into two parts: the rise of the input signal, $\phi =0\to \phi =2\pi$ and its fall, $\phi =2\pi \to \phi =0$. As can be seen from Figure 3a,c current oscillations in the course of evolution reach stationary values A_↑ and A_↓ (solid and dashed curves in Figure 4a, respectively), except for the transition region (hysteresis loop in Figure 3c). Dependencies in Figure 4 show that at low values of the inductance l in the circuit (up to l < 0.5), the system evolves adiabatically slowly between the equilibrium positions $\phi =0\iff \phi =2\pi$. In this case, the amplitude of oscillations on the transfer function is equal to zero even with an increase in inertia (see curves for different values of the Stewart-McCumber parameter in Figure 4a). The contribution of kinetic energy, and, consequently, an increase in the amplitude and period of oscillations is observed for the logic element mode. This is manifested by the fact that the amplitudes of plasma oscillations are correlated, but in absolute values A_↓ and Tp↓ are 1.2 times higher than A_↑ and Tp↑. At the same time, with an increase in the capacity (Stewart-McCumber parameter), there is an increase in the amplitude of oscillations during the “reverse way” ($\phi =2\pi \to \phi =0$) and the period of oscillations, Tp, which is associated with an increase in the effective “mass” for the system.

3.2. Influence of the Dissipation

Undoubtedly, an important point in studying the principles of operation of an artificial S_C-neuron is taking into account the processes associated with energy dissipation due to the appearance of voltages on the resistive element and the destruction of the superconducting state. We have investigated the dynamic processes in the system based on the solution of the system of Equation (9) for two operating modes: the neuro-mode at l = 0.22 (Figure 5a) and the bistable logical element mode at l = 2.5 (Figure 5b). The presence of dissipation in the system is associated with a nonzero rate of change in the phase of the Josephson junction under the action of a current pulse. This leads to the suppression of plasma oscillations on the transfer characteristics $i_{out}({\phi }_{in})$, as well as the establishment of a stationary value of the output current for the logic element mode (Figure 5b). Note that the sigmoidal form of the activation function is preserved, which indicates the possibility of creating a superconducting neuron without resistive shunting (with the large Stewart-McCumber parameter). We calculated the dissipation energy $E_d=\int {\left(u_{JJ}(t)\right)}^{2}dt$, where $u_{JJ}(t)=p{(t)}^2/2$ is the voltage across the Josephson junction under prolonged external action ${\phi }_{in}(t)$. The inset in Figure 5a shows that an increase in the Stewart-McCumber parameter β leads to an increase in energy losses in the system. However, we can choose the parameters of the circuit so as to suppress unwanted processes of inertia (plasma oscillations) and at the same time reduce the geometric dimensions of the artificial neuron circuit (by reducing shunting), which is important for integrating it onto a chip. For typical parameters of the neuron circuit (I_c = 100 mA, R = 1 Ω, C = 10 pF), we estimated the operation time: τ = 10 ns (ω_p~10¹²c⁻¹). The dissipation energy in the neuro-mode (l = 0.22) E_d = 1.6 aJ; for the logical mode (l = 2.5) E_d = 75 aJ at the value of the Stewart-McCumber parameter β = 5.

3.3. Influence of Temperature on the Dynamics of a Neuron

Let us now analyze the influence of the initial distribution of the phase and its velocity on the behavior of the trajectories of the system in the phase space of states and on their final distribution that arises after the action of the control flux. In other words, we are interested in the effect of thermal “blur” on the output characteristic of the device.

We carried out a numerical solution of the system of Equation (9) by the Monte Carlo method to take into account the effect of temperature on the dynamic characteristics and phase portraits. We believe that at the initial moment of time the S_C-neuron is in thermodynamic equilibrium with the external environment at a final temperature T. In this case, the numerical method consists in generating random initial states of system (9) according to the Gibbs distribution law:

$$\rho (\phi ,p,t=0)=\frac{1}{Z}{e}^{-\frac{E_{J}H(\phi ,p,t=0)}{kT}},$$

(10)

where $H(\phi ,p,t=0)$ is the dimensionless Hamiltonian defined by expression (7), Z is the partition function, k is the Boltzmann constant, and $E_J=\frac{\hslash I_c}{2e}$ is the characteristic Josephson energy. The evolution of a set of states in time gives an idea of the density distribution function at each of these moments.

We investigated the dynamics of 10⁴ particles for switching a neuron at a characteristic temperature of 4 K. The initial conditions are specified according to distribution (10) for the case of a neuro mode for a potential with one local minimum (l = 0.22, solid black curve in Figure 2) and a sigmoidal transfer characteristic (see Figure 6a). In Figure 6 we observe that all particles move adiabatically slowly $\phi =0\iff \phi =2\pi$ (the particle velocity is less than the rate of potential change), which can be seen in Figure 6b at time t_m (see Figure 1b). The color in Figure 6 and lines show equipotential levels, corresponding to the same energy. It was found that the form of the dynamic characteristic after averaging over the ensemble of particles remains similar to Figure 5a, which indicates the stable operation of the proposed artificial neuron at a finite temperature.

For the case of a logical element (l = 2.5), we investigated the phase trajectories of an ensemble of 10⁴ particles. As can be seen from Figure 7, all particles for the rising phase of the input signal $\phi =0\to \phi =2\pi$ evolve in a correlated manner and acquire a final momentum, moving in a circle in Figure 7b. Further, due to giving the system an additional pulse at the signal decay interval (Figure 7c), a rearrangement of particles is observed. Some of them evolve to the initial state $\phi =0$ at the time of the termination of one pulse cycle, and the second group of particles continues to move with a finite momentum (Figure 7d).

Thus, taking into account the capacity, the final state of the system depends on the ratio between the effective depths of the wells and the temperature. In this case, taking into account the effects of dissipation leads to the localization of the distribution of states (the mode of attractors), similar to what was shown in Figure 5 for one path. In addition, as a result of a numerical calculation, it was found that the final distribution of states depends on the switching time; in this case, with an increase in the duration of the input signal (a decrease in the parameter D in Expression (1)), the density of states becomes more localized (occupies a smaller area). In the adiabatic switching mode (when the particle velocity is less than the rate of potential change), the distribution function becomes more and more symmetric with respect to the initial distribution, and the finite temperature leads to smearing.

4. Conclusions

We have studied the dynamics of a superconducting neuron based on non-shunted Josephson contacts (S_C-neuron) in different situations. Earlier, the functioning of a superconducting neuron was discussed only within the framework of the inertialess model, in which the output current instantly follows the change in the control magnetic flux (instantaneous flux-to-current transformation). As it was shown in this paper, if we take into account the capacity of the Josephson contact, this leads to a qualitative change in the switching mode: i) if the element in “logical mode” was initially in the equilibrium position and its state was represented by a point in the phase space (phase-conjugate momentum); ii) then after half of the switching cycle, its states are described by the dots on a closed curve. After single cyclic switching, the element does not return to its original state, which indicates its hysteretic behavior. This property can be used to remember the history of the process. We also investigated the effect of dissipation on the switching process and obtained a picture of the formation of an attractor in the phase space for the realistic system parameters. An interesting prediction is related to the effect of the initial thermal “blurriness” on the states of the system. As it turned out, the initial heat spot in the phase space undergoes a complex transformation after a switching cycle.

We believe that the newfound features of the switching process can be attributed to an additional appealing step to design superconducting neurons and neuromorphic networks. The results from the analysis estimates of the performance and energy efficiency of the analog part of hybrid analog-digital “adiabatic” superconducting circuits allow us to raise the question of the prospects of their application in hardware accelerators of neural networks. The utilization of verified technologies for creating basic components of a universal set of logic elements (with reasonable values of the capacity of the Josephson junctions) allows us to count on the creation of a neuroprocessor with target performance indicators of the order of 10⁻¹¹J and energy consumption at the level of 10⁻¹⁸J per operation of “calculating” the activation function. Examples of practical tasks with high requirements for performance, energy efficiency, and the accuracy of synchronization of the hardware platform that require the development of such devices are listed below [24,25]:

• Internet search using neural networks-transformers that places extremely high requirements on the hardware platform;
• Object detection in aerial and satellite images;
• Analysis of network traffic in order to ensure network information security and neutralize Internet fraud.