Monte Carlo simulations of the switching processes in the superconducting quantron-based neuron

We study the response function of a superconducting single-contact interferometer (quantron or rf-SQUID), which is the building block for adiabatic neuron with a nonlinear transfer characteristic. It is shown that the intrinsic shunt capacitance of the junction leads to distortions in the response function. At the same time the contact resistance leads to the suppression of these vibrations. The response function of the rf-SQUID at finite temperature is calculated by Monte Carlo method.


Introduction
Superconducting interferometers with embedded Josephson junctions are widely used in quantum electronics and computer engineering. One of the modern applications of the rf-SQUID (or quantron) is to create quantum neurons elements which may be considered as basic of neural networks. Recently, an inertialess model model of a neuron has been developed [1][2][3][4]. Currently, there is no mathematical neuron model that takes into account both the intrinsic contact capacity ("inertia") and dissipative processes in the Josephson junction.
Despite the quantum nature of the rf-SQUID, there is a range of parameters where it is possible to describe the operation of a neuron in a classical way. In this paper, we study the switching process in the rf-SQUID taking into account the shunting capacity and dissipation in the junction. We are interested in the presence of bistable rf-SQUID states and transfers between them under the action of an external magnetic flux. Numerical simulations of the rf-SQUID dynamics in the framework of a resistive model have been performed. It's assumed that initially, before the external magnetic flux through the loop of the rf-SQUID is turned on, it is in thermodynamical equilibrium with the medium. An ensemble of the junction "copies" was introduced to take this effect into account. The switching process at finite temperatures was modeled using the Monte Carlo method.

The model and basic equations
The system under the study (see figure 1) consists of a superconducting loop interrupted by a single weak link through which the Josephson current can flow where I C is the value of maximum amplitude of a superconducting current, and the ϕ is the charge-invariant phase. The voltage V across the junction follow from Josephson equation There is always a shunting capacitance, C, and normal current pass across the junction. Taking into account the capacitance and contact resistance of the junction we can write the total current through the junction in the form The main characteristic of a neuron is the nonlinear dependence of the output variable ϕ versus the external magnetic flux Φ ex . If the rf-SQUID's loop is treated by the Φ ex flux, the resulting current is described by the expression where L q is the loop's inductance.
In result of combination the expressions (3) and (4) we obtain a differential equation in the normalized variables that determines the relationship between the normalized phase ϕ(t) and the normalized control parameter ϕ ex = 2π Φex is the characteristic Josephson frequency, γ = ω p /ω d is the attenuation parameter. If we neglect the dissipation the dynamic of the system my be described by a Hamiltonian where the dimensionless momentum p is measured in the values is the measure of inertia, E J is the energy of the Josephson coupling on which the Hamiltonian is normalized.

The switching process of the rf-SQUID
Let's start with the case of a inertialess and non-dissipative rf-SQUID. In this case the Eq. (4) reduces to ϕ + l q sin ϕ = ϕ ex , where l q = 2πI C Lq Φ 0 = L q E J is a dimensionless coefficient that characterizes the degree of nonlinearity of the system. Notice, that the equation we can treat as minimum of the potential U (ϕ, ϕ ex ) . Despite the fact that in practice it is easier to measure the directly connected with current flux through the loop, this dependence can be considered as an activation function of the neuron.
The rf-SQUID's activation functions for different values of the parameter l q (figure 2) were obtained by numerically solving the transcendental Eq. (7) for a set of points ϕ ex in the range from 0 to 2π. For l q > 1, the function is multi valued (figure 2 b), and different colors of the point define the root search area. The time dependent bias turns the static activation function of the rf-SQUID into a timedependent one. It worth to remember that the inverse plasma frequency ω p was chosen as the time unit. In this case, the value of ϕ depends on both the control parameter ϕ ex and its history of changes according to Eq. (5). The case when the control function has the form of a smooth step is being considered, where ϕ ex (0) is the initial value of the external magnetic flux and the parameter τ defines the smoothness of the curve. The dependence of ϕ(ϕ ex ) is obtained by solving the differential Eq. (5) on the interval of the characteristic switching time τ and parametric plotting of the curve with respect to time. For this purpose an algorithm was written and implemented, which resulted in a function that completely coincides with the corresponding dependence obtained by applying the system of Hamilton equations to the expression (6). As it's shown in figure 3, the output variable ϕ oscillates after switching on the characteristic time ω −1 p .  Figure 3. The function ϕ(ϕ ex ), which is obtained using a system of Hamilton equations for a single virtual particle with the Hamiltonian Eq. (6) for parameter l q = 0.95 (a) and l q = 3, τ = 1000, γ = 0 (b). Initial conditions are defined as: ϕ(0) = 0, ϕ (0) = 0.

The Monte Carlo method
With the presence of the capacity, the dynamics of the rf-SQUID can be realized as the dynamics of a particle with the coordinate ϕ and momentum p, which obey the system of Hamilton equations.
Now we want to take into account the termodynamical fluctuations of the initial coordinate ϕ and momentum p. Let the junction be in the thermodynamical equilibrium with the medium at a finite temperature of T . Suitable parameters of the Josephson contact were selected for modeling. At the initial moment there were generated 100 000 copies of virtual particles, which coordinates in phase space obey the Gibbs distribution ρ(ϕ, p) = 1 Z e − H(ϕ,p) kT , where Z is the partition function and k is the Boltzmann constant. The coordinate ϕ and the momentum p of each copy are the initial values for the equations of motion, which have been solved by the fourth-order Runge-Kutta method. As can be seen from the figure figure 4, that the initial distribution, which approximately looks like a gaussian distribution, evolves into the shape of a "necklace". This is due to the fact that particles acquire CSP2020 Journal of Physics: Conference Series 1740 (2021) 012063 IOP Publishing doi:10.1088/1742-6596/1740/1/012063 5 energy due to adiabatic changes in the potential and accumulate in the phase space to an ellipse with the corresponding energy. As can be seen from the figure 5 (a), the average coordinate over the ensemble oscillates with a lower amplitude for the non-dissipative case than the coordinate of the single trajectory ( figure 3). It can be noticed that the energy of the system increases as a result of switching, therefore the oscillations occur. The number of oscillations depends on the switching time as expected. The presence of dissipation leads to localization of the distribution of states, so the function ϕ(ϕ ex ) is similar as for a single trajectory. As a result of numerical simulation, it was found that the final distribution of states depends not only on the switching time τ (as τ increases, it becomes more symmetrical and occupies a smaller area), but also on ϕ ex (0).

Conclusion
It is shown that the capacity and dissipation can affect qualitatively the dynamic of the rf-SQUID in external magnetic flux. It is found that these effects lead to oscillations of the transition phase characteristic and to an unusual distribution function of the phase in the case of a finite temperature. These properties must be taken into account when designing neurons based on Josephson elements.