
==== Front
Proc Natl Acad Sci U S A
Proc Natl Acad Sci U S A
PNAS
Proceedings of the National Academy of Sciences of the United States of America
0027-8424
1091-6490
National Academy of Sciences

38470918
202314995
10.1073/pnas.2314995121
research-articleResearch Articlebiophys-physBiophysics and Computational Biology0
Physical Sciences
Biophysics and Computational Biology
Collective neural network behavior in a dynamically driven disordered system of superconducting loops
Goteti Uday S. ugoteti@ucsd.edu
a 1 https://orcid.org/0000-0002-9430-451X

Cybart Shane A. b https://orcid.org/0000-0001-9890-6537

Dynes Robert C. rdynes@ucsd.edu
a 1 https://orcid.org/0000-0001-6740-9677

aDepartment of Physics, University of California, San Diego, CA 92093
bDepartment of Electrical and Computer Engineering, University of California, Riverside, CA 92521
1To whom correspondence may be addressed. Email: ugoteti@ucsd.edu or rdynes@ucsd.edu.
Contributed by Robert C. Dynes; received August 29, 2023; accepted February 12, 2024; reviewed by Steven A. Kivelson and Terrence J. Sejnowski

12 3 2024
19 3 2024
12 9 2024
121 12 e231499512129 8 2023
12 2 2024
Copyright © 2024 the Author(s). Published by PNAS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Significance

In this paper, we describe a very different architecture for complex computations. It is a step in the direction of modeling and experimentally demonstrating a highly interactive system which we think more represents the functioning of our brains. We bring together the properties of superconducting loops linked together with Josephson junctions and the physics of disordered structures. The information is transmitted via superconducting magnetic vortices transmitted between loops through junctions and the delocalized long-range coherence is assured by the superconducting properties. The disorder of the constituent loops offers an exponentially increasing density of configurational states resulting in static and dynamic memory. We have shown both by modeling and experiment that this circuit has intrinsic associated memory: a prerequisite for cognition.

Collective properties of complex systems composed of many interacting components such as neurons in our brain can be modeled by artificial networks based on disordered systems. We show that a disordered neural network of superconducting loops with Josephson junctions can exhibit computational properties like categorization and associative memory in the time evolution of its state in response to information from external excitations. Superconducting loops can trap multiples of fluxons in many discrete memory configurations defined by the local free energy minima in the configuration space of all possible states. A memory state can be updated by exciting the Josephson junctions to fire or allow the movement of fluxons through the network as the current through them surpasses their critical current thresholds. Simulations performed with a lumped element circuit model of a 4-loop network show that information written through excitations is translated into stable states of trapped flux and their time evolution. Experimental implementation on a high-Tc superconductor YBCO-based 4-loop network shows dynamically stable flux flow in each pathway characterized by the correlations between junction firing statistics. Neural network behavior is observed as energy barriers separating state categories in simulations in response to multiple excitations, and experimentally as junction responses characterizing different flux flow patterns in the network. The state categories that produce these patterns have different temporal stabilities relative to each other and the excitations. This provides strong evidence for time-dependent (short-to-long-term) memories, that are dependent on the geometrical and junction parameters of the loops, as described with a network model.

brain-like computation
complex systems
superconductivity
disordered systems
Josephson junctions
DOE | Office of Science (SC) 100006132 DE-SC0019273 Uday GotetiRobert C. Dynes DOE | National Nuclear Security Administration (NNSA) 100006168 DE-NA0004106 Shane Cybart DOD | USAF | AMC | Air Force Office of Scientific Research (AFOSR) 100000181 FA9550-20-1-0144 Shane Cybart
==== Body
pmcEmergent behavior exhibited by complex systems can be statistically evaluated as the effects of collective spontaneous activity of a large number of its individual constituents and time-dependent interactions between them. This behavior is generally observed at the system level in networks of coupled nonlinear elements, e.g., neurons in our brain, despite the differences in many of the details of the individual components and interactions that compose them. A Hopfield network (1) models such behaviors based on disordered physical systems, e.g., spin glasses that are composed of a network of coupled individual spins. These systems allow an almost continuum of distinct energy levels to describe computational properties in time as due to the phase-space flow of the state of the network in response to external excitations, as it dynamically equilibrates with its environment and converges onto the nearest local free energy minimum. A physical manifestation of a disordered network of elements with short-to-long-range interactions can be operated to store and retrieve information through excitations and the generated responses in time. The inclusion of feedback mechanisms to actively reconfigure the energy landscape of the state space enables them to be implemented in neuromorphic computing and to model behaviors of complex systems.

In this paper, we show a disordered neural network of superconducting loops with Josephson junctions and demonstrate its operation both experimentally and in simulations. Each superconducting loop can trap multiples of magnetic flux quanta (Φ0 = 2.067×1015 T/m2) or fluxons in the form of either clockwise or counterclockwise circulating supercurrents around the loop. A multi-loop disordered system, schematically shown in Fig. 1A, can host a large number of distinct flux configurations that are either statically or dynamically stable defining a discrete state space. A subset of these can be configured to be local energy minima or attractors. Each state defines an interaction strength between pairs of different Josephson junctions and any state transition triggers a response that can be measured by the voltage across them. Active external locally applied currents or flux stream excitations drive the system through different states and into the nearest local minimum after they are turned off. We show that a small-scale network of four loops can be simulated using a lumped circuit model with Josephson junctions and inductors to reproduce experimental observations obtained from an equivalent YBCO-based loop network. The following section discusses, through simulations, the operation of a superconducting disordered system and demonstrates addressing the memory states and re-configuring the state space. Then, a network model is presented to relate the simulation and experimental results presented in the following section. The dynamic evolution of the state obtained from simulations statistically correlates external excitations to locally measured responses across junctions observed in experiments. Applying multiple local excitations leads to neural network behavior with memory states defining dominant flux flow (information flow) pathways through it. These flux pathways exhibit wide-ranging temporal stability with dynamically changing external excitations depending on the loop geometry and coupling from the relatively stable trapped flux during the dynamic flow.

Fig. 1. (A) Superconducting loop disordered neural network with multiple spiking inputs and outputs. A typical trapped flux configuration with circulating currents around individual loops is shown as n1Φ0 and n2Φ0 and around multiple loops shown as n345Φ0. Each loop can be treated in simulations using an equivalent lumped element circuit model of inductors and Josephson junctions as shown. Feedback currents are useful to reconfigure the mapping between input signals and the resulting memory states. (B) Configuration space with disordered energy landscape with coordinates defined by all the possible trapped flux states. Stable memory configurations are given by the local energy minima. Any of the trapped flux states can be dynamically reconfigured to be local minima with suitable excitations (inputs or feedback). (C) Model of a 4-loop network showing possible ways (arrows) to excite the system to access various flux configurations. (D) An equivalent lumped circuit model of the 4-loop network was implemented to generate the simulation results discussed in the article. (E) YBCO-based 4-loop network with Josephson junctions defined by focused He-ion beam irradiation-induced insulating tunnel barriers labeled by white lines. (F) The number of flux configurations and the energy calculated from circulating currents from the circuit model show the density of states from zero trapped flux to maximum filled flux state for the 3-loop memory network.

Flux Configurations

The number of all possible trapped flux configurations in a superconducting disordered system is determined by its geometrical or junction parameters and limited by flux quantization. A closed superconducting loop can sustain persistent circulating currents that aggregate to an integer number of fluxons in the absence of any external excitations. Josephson junctions are the weak links through which fluxons can enter or exit a loop generating a voltage response in time (i.e., ∫0tVdt=nΦ0) when an excitation current together with circulating current from previously trapped flux exceeds their critical current barriers. The junctions in a disordered loop system can hence be treated as equivalent to the active neuron-like components of the network interacting through the trapped flux overlapping them. Nonhysteretic (damped) junctions are considered in simulations and experiments that can show single fluxons traversing them in the form of discrete voltage spikes. An integer number of fluxons can be trapped as circulating currents around individual loops (e.g., n1Φ0 and n2Φ0 around loops 1 and 2 in Fig. 1A) that constrain the interactions to a short range, i.e., between the junctions in that loop. Due to macroscopic coherence across the superconducting system, the network allows long-range direct interactions between junctions through trapped flux where circulating currents can encompass several loops (e.g., shown in Fig. 1A as n345Φ0 involving loops 3, 4, and 5). Flux configurations are typically composed of a combination of interactions ranging from short-to-long length scales in a network where interactions (both direct and indirect) between any pair of junctions can be described as due to a cumulative effect of the flux trapped in each closed-loop pathway encompassing either of them. In the rest of the article, the flux configuration or the memory state of the network is defined collectively by the integer number of trapped fluxons nk for each of the circulating current paths k (e.g., k=1,2,45,367, etc.). nk is considered to be positive for clockwise circulating currents. These configurations can be represented as co-ordinates with their respective energies in an appropriate multi-dimensional configuration space C as shown in Fig. 1B. Interaction between any pair of junctions can be controlled by traversing the configuration space to a desired state through external excitations that define the dimensionality of C.

Fluxon dynamics through the loop network can be modeled using lumped circuit analysis similar to ref. 2. Its motion through superconducting loops with junctions is comparable to angular strain propagation through torsion bars with suspended pendulums. Pendulums are analogous to Josephson junctions and are treated using their circuit equivalents described by the resistively and capacitively shunted junction (RCSJ) model. The superconducting paths between the junctions similar to the torsion bar connecting the pendulums are treated using their equivalent inductances as shown in Fig. 1A, so the angular strain in the torsion bar can be compared to the current through the inductor. The dynamic disordered system therefore behaves as a multi-dimensional array of coupled oscillators interacting over a range of lengths and time scales.

Each circulating current path k can individually accommodate a maximum trapped flux with an integer value Nk=LkIkjΦ0, where Lk is the total inductance of the circulating current path and Ikj is the critical current of the junction j where flux traverses out of that path, which is dependent on its interactions with all the junctions in the network. Effects of state fluctuations from thermal noise and quantum tunneling below junction critical currents similar to ref. 3 are neglected in the simulation model as they are expected to be negligibly small. The simulations are still consistent with experiments for statistical analysis of fluxon flow with long-time averaging as described by the neural network model in later sections. The energy of the network (flux) configuration as a function of time can be calculated from the circuit equivalent as the sum of inductive and Josephson potential energies stored in all the circuit components in the network due to cumulative currents across each element from the trapped flux and external excitations. Excitations can be considered as flux flow (voltage spike trains) into the network through Josephson junctions or as continuous current signals as shown in Fig. 1C. The flux flow excitation inputs generate a response that can be recorded as flow through output junctions, with the disordered superconducting network mimicking a spiking neural network with active currents acting as feedback signals (Fig. 1A). These feedback allow us to reconfigure the mapping between input signals and the resulting memory configurations. When driven with continuously changing excitations, the evolution of the system through the dynamic configuration space is observable through the time-dependent voltage signals across different output junctions.

Flux Configuration Mapping in a 4-Loop Network.

A simple network of 4 loops with 6 junctions is studied using the circuit model shown in Fig. 1D and the experimentally equivalent YBCO-based network with Josephson junctions defined by a focused He-ion beam induced tunnel barrier similar to ref. 4 shown in Fig. 1E. Loops 1 to 3 are constructed to independently accommodate memory states in the form of multiples of fluxons, i.e., Nk(=1,2,3)>1. Loop 4 is constructed as flux input loop such that Nk(=4)<1 and any fluxon introduced through junction J1 is stabilized across the three memory loops. To allow a rich spectrum of states disorder is important in the function of the circuit and is introduced by randomizing the network parameters such as the loop inductances, junction critical currents, and the network topology. The input current I1 induces a voltage V1 across junction J1 and the output flux stream is measured as voltage V2 across J6 representing a single channel of the neural network. A second input current I2 is used as feedback to reconfigure the state space relative to the input.

Considering the 3 memory loops, a flux configuration is represented collectively by the respective circulating currents in 7 different closed-loop pathways as C=[n1n2n3n12n13n23n123]. The number of such states as a function of their energy is shown in Fig. 1F including all possible configurations of trapped flux, calculated by iterating through combinations of nk up to their maximum Nk that can coexist without surpassing critical currents of any junctions. These configurations comprise most of the available memory states of this 4-loop network. In case of symmetry or order in the individual loops or the coupling between them, iterating through all the possibilities of k results in counting the same state multiple times (e.g., when the same set of circulating currents can be counted twice in n12 and n1+n2). If the network is completely disordered, the number of possible states increases exponentially with the increasing number of loops in the network and theoretically any of these states can be configured to be meta-stable local energy minima with a complex combination of time-dependent input and feedback signals, as demonstrated in the sections below.

Experimentally, a transition from any given memory state CM1 to another state CM2 can be initiated by applying multiple different combinations of time-dependent external excitation signals to traverse the dynamic configuration space followed by a relaxation period. The flux configuration at any instant of time can be known with the knowledge of the initial memory state and the flux flow through all the junctions entering or leaving different loops during state evolution. In all the simulations discussed in the article, various states of the 4-loop network are addressed using different excitation signals starting from an initial configuration of zero flux trapped in the system.

Magnetic field excitation.

The flux states of the disordered network can be uniquely addressed by applying a uniform magnetic field across the entire network and allowing it to reach equilibrium before the excitation is turned off. The system will then find the nearest stable local energy minimum. The simulation model shown in Fig. 1D is subjected to a time-dependent uniform magnetic field through a pulsed current with a short rise/fall time of 1 ps across an inductively coupled coil to the outer inductors of the 4-loop network as illustrated in Fig. 2A. At the onset of the magnetic field pulse, flux begins to enter the network through Josephson junctions as the state of the system evolves in time and approaches equilibrium with the external field. If the pulse width of the field is longer than the critical time period required to reach this equilibrium state, the final configuration reached after relaxation is independent of further increase in pulse width as shown in Fig. 2B. This critical width represents an important timescale during which the effects of dynamics of evolution of the flux state and the excitations cannot be completely separated in the statistical analysis of the observed flux flow responses. Its value is specific to the input excitations and is related to the initial state, the applied field amplitude, and the loop parameters as described using the network model. For the simulated network, the critical excitation time periods for different configurations are found to be in the range of 0 to 500 ps in simulations for different applied fields starting with a zero-flux state. Therefore, a pulse width of 1 ns is used to map different field strengths to flux configurations as shown in Fig. 2C–E. Fig. 2C shows the flux movement in both directions through different junctions in the network during excitation (i.e., from 0 to 400 ps) with an equivalent field of 500ΦΦ0 followed by relaxation (from 1 to 1.2 ns). The relaxation time is similarly related to the loop parameters. Voltage spikes equivalent to individual fluxons can be observed traversing some junctions (e.g., J2, J6) while a stable flux configuration is reached.

Fig. 2. (A) Schematic illustrating simulations of the uniform magnetic field across the 4-loop network described in Fig. 1C–E. A magnetic field pulse is applied with a pulse width of 1ns and different pulse heights and then the system is allowed to relax to realize different meta-stable trapped flux configurations. (B) Simulation of time-evolution of the state of the system due to magnetic field excitation. The relaxed state is recorded with pico-second increments to the magnetic field pulse width to map the time evolution for 5 different pulse amplitudes. This state is independent of time after a critical excitation time period (pulse width) dependent on the loop parameters of the involved trapped flux configurations. (C) The voltage across 5 junctions in the 3-loop memory network of Fig. 1C–E showing multiples of quantized flux pulses entering and leaving the network through each of the junctions for the magnetic field of 500ΦΦ0 with a pulse width of 1ns. There is no flux motion after the critical excitation time period. (D) The energy of the state of the network in Fig. 1C–E with the respective magnetic field for a fixed pulse width of 1ns at different times during the relaxation process. Discrete states reached by 1.5 ns (500 ps of relaxation) represent the local minima in configuration space. (E) The energy of the relaxed state is shown to demonstrate reconfiguration of the state-space with different constant feedback currents at I2. Feedback current has the effect of tilting the state space to access different memory states.

The energy of the state at different times during the relaxation period following the excitation pulse is shown for different magnetic fields in Fig. 2D. The state space is symmetric with respect to the sign of the applied field, although the local magnetic moments are expected to be of different directions aligned to the fields. Discrete energy levels of relaxed flux configurations can be seen beginning after 100 ps, and the network is completely relaxed at 500 ps for the entire range of field amplitudes. These states are retained indefinitely in the absence of further positive or negative field excitation. Each state can be addressed uniquely with field amplitudes extending over a range, e.g., the state at ≈8eV can be accessed with multiple different field amplitudes in the range of 230ΦΦ0 to 280ΦΦ0. Other states at ≈6 eV, 9 eV, and 11 eV partially overlap with this range of fields, but this overlap is not exact, i.e., any given applied excitation maps to a unique state. These states are adjacent to each other resulting in some uncertainty in the final state in the presence of noise in the excitation. The system reaches a maximum filled flux state along this coordinate axis at a field amplitude of ≈450ΦΦ0 and is unchanged with further increase in the field.

The mapping between the fields and the resulting memory states can be reconfigured through a feedback current applied to the system (shown as I2 in Figs. 1 D and E and 2A). Application of a constant feedback current has the effect of tilting the state-space mapping between the inputs and the final memory states, shown in Fig. 2E for I2 = 50μA, 100μA and 200μA. States that are not stable previously can be addressed after reconfiguration and are retained even when the feedback current is reduced to zero if the system is allowed to reach equilibrium. In a network with multiple excitations as in Fig. 1C, some of the available excitation channels can be used as feedback to give access to a large number of memory states shown in Fig. 1F. A uniform excitation across the entire system in equilibrium conditions therefore shows a simplified time-independent operation of the disordered network where continuous input information can be mapped onto unique discrete states.

Local spiking signal excitation.

The state evolution shown for the superconducting system subjected to uniform excitations describes conceptually equivalent dynamics when compared to various experiments performed on other disordered systems, characterized by specific heat capacities of glasses (3) or magnetic susceptibilities of spin glasses (5). However, the networks described in this work allow us to generate local excitations to disordered systems that represent inputs to some neurons to generate local responses to be measured at other output neurons in a neural network. The 4-loop system of Fig. 1 D and E shows a simplified version with the junctions at the outer-most loop, i.e., J1, J3, and J6 as either inputs or outputs for local excitations or responses. Currents I1 and I2 generate flux flow excitations locally that can be used for either input or feedback. An external current pulse excitation I1 is applied in simulations to induce an input flux flow at a constant rate into the network locally through the junction J1 for a fixed time period as described by Fig. 3A. Flux flow rate can be quantified from the number of discrete fluxons traversing the junction over a fixed time period, characterized by the constant average frequency or the voltage across it as V1Φ0 from the AC Josephson effect. Similar to uniform fields, a sufficiently long period of local flux flow excitation longer than a critical time period allowed the network to reach a dynamic equilibrium. A period of 2 ns is observed to be sufficiently long for the network to reach this state giving a steady flux flow into or out of the network through all the junctions on the outer loop. Therefore, a constant flux flow excitation is applied for 2 ns and the network is similarly allowed 500 ps to completely relax to respective local energy minima for different local flux flow inputs through V1. The flux states shown by their respective potential energy stored are mapped to the input flux flow rates at different instances during the relaxation period using simulations in Fig. 3B. The state space is qualitatively similar to that of magnetic field excitation but with a larger number of accessible states and a stronger overlap. The final states are dependent on both the flux flow rate and the time period of excitation, as each fluxon during the flow induces a state transition. The configuration map in Fig. 3B is specific to the chosen time period of 2 ns.

Fig. 3. (A) Schematic illustrating simulations of spiking excitation induced at J1 in the network described in Fig. 1C–E. An input flux flow at a constant rate is applied for a fixed pulse duration of t1 = 1 ns for different pulse heights I1 and then the system is allowed to relax to realize different meta-stable trapped flux configurations. (B) The energy of the state of the network in Fig. 1C–E with the respective excitation voltage V1 (i.e., average flow rate of V1Φ0) for a fixed pulse width of t1 = 1 ns observed at different time intervals (i.e., at t2 = 1.01 ns, t3 = 1.05 ns and t4 = 1.5 ns) during the relaxation process. Discrete states reached by 1.5 ns represent the local minima in configuration space. (C) The energy of the relaxed state is shown to demonstrate reconfiguration of the state-space with active constant feedback I2=700μA showing neural network behaviors such as categorization and associative memory. (D) Simulation of time-evolution of the state of the system due to local excitations (i.e., I1=1 mA and I2=0). The state is recorded after relaxation as a function of pulse width t1 with picosecond increments to map the time evolution. (E) Firing probability p1 of junction J6 with respect to J1 as a function of integration time T during the time evolution of the state shown in Fig. 3D. (F) Simulation of time-evolution of the state of the system due to local excitations (i.e., I1=0.9 mA and I2=0.7 mA). The state is recorded after relaxation as a function of pulse width t1 with picosecond increments to map the time evolution. (G) Firing probability p1 of junction J6 relative to J1 shown as a function of integration time T during the time evolution of the state shown in Fig. 3F. (H) Simulation of time-evolution (i.e., the relaxed state after picosecond increments) of the state of the system due to dynamically varying input excitation, i.e., I1 is linearly varied from −1.5 to 0 mA in an interval of 6 ns with constant feedback current I2=0.7 mA. With time-dependent input, the trapped flux configurations that generate the output flow have different temporal stability.

Current through I2 can be used as feedback to reconfigure the mapping between input flux flow and the final flux configurations. Simultaneously applying multiple different local excitations (input and feedback) results in neural network behavior exhibiting categorization and associative memory. An example considering a constant feedback current of I2=700μA is observed to introduce asymmetry to the state space resulting in the separation of states into categories seen as energy contours in Fig. 3C. Similar behavior is observed for other values of I2 as shown in the experimental results. State categories can be observed while the excitation currents are on, within regions of V1 between 0 to 0.3 mV, −0.5 to −0.25 mV, etc., and are retained within the same regions even when one or both the currents are turned off. Specifically, the network is in a steady flow state after the input I1 is turned off with I2 turned on and the state categories that were originally created with two excitations are retained with energy gaps separating them. A gap of ≈15 eV separation can be observed at V1 of −0.1 mV, and similarly, a gap of ≈25 eV is seen at V1 of 0.25 mV. We will find that these state categories represent different flow patterns of flux through the network in experimental observations. Even though the output flow measured after I1 is off is expected to be different, this behavior is characteristic of associative memory between the two inputs I1 and I2. This is because the output flux flow measured before and after I1 is turned off is statistically equivalent (caused by the same category of states), and the flux flow pattern is effectively unchanged as discussed in the sections below.

The discrete states representing uniquely addressable flux configurations in Fig. 2 D and E or the strongly overlapping configurations in Fig. 3 B and C are the quiescent states of the network reached after complete relaxation. Therefore, information can only be retrieved through a detailed nonintrusive scan of trapped flux across the entire network. Alternatively, the memory states can be addressed by locally induced flux flow excitations and measuring the resulting flux flow responses, which are a direct consequence of the interactions between the input and the output junctions defined by the underlying trapped flux. This is a destructive readout mechanism as any excitation to the system results in the perturbation of previously written memory. A constant flux flow is induced at J1 using current pulses at I1 of increasing excitation pulse widths with pico-second increments, each followed by a relaxation period of 500 ps. The state evolution in time with each additional input fluxon, obtained from such simulations is shown in Fig. 3D. With an initial state of zero flux trapped, the state traverses through the configuration space until it reaches a dynamic equilibrium by 400 ps as it converges on to a small region observed between 2 and 4 eV. The network appears to mimic stochastic fluctuations in this region at short time scales although these fluctuations are observed to be deterministic and re-traceable. This behavior is seen in the absence of noise due to the nonlinearity of the Josephson junctions and the disordered coupling between them. In dynamic equilibrium, the interaction strength between the input junction J1 and the output junction J6 is dependent on a localized category of states around the stable local minima that can be measured by statistical correlations between the switching activities of the two junctions. Therefore the approximate address of the memories in the vicinity of these states can be experimentally measured from the relative flux flow rates, observed over a period to obtain the firing probability of a junction relative to another. Firing probability can be quantified as the number of fluxons exiting J6 due to a fixed number of fluxons entering J1 integrated over a time period T i.e., p1=∫0TV2dt∫0TV1dt. p1 is shown as a function of integration times starting from t=0 during the time-evolution of the state (Fig. 3D) in Fig. 3E. This quantity can also be considered as an equivalent to the synaptic weight between the two junctions and reaches a constant value after an integration time longer than the critical excitation time period. It is independent of time sufficiently long after the system reaches a dynamic steady state, similar to the time-independence observed in static equilibrium during uniform field excitations shown in Fig. 2. Another example of dynamic evolution of the state that results in a different firing probability between J1 and J6 with both I1 and I2 on is shown in Fig. 3 F and G. When input flux flow is constant, the dynamic memory state of the system within its vicinity in C is permanently stable. However, a much more general description of a neural network includes networks subjected to time-varying input and feedback excitations. Such an example is shown in Fig. 3H with a constant I2 but linearly increasing input current I1 over a period of 6 ns. The state of the system follows the variations in input as it evolves, but each of the underlying trapped flux states has different temporal stability relative to the changing input excitation resulting in shorter and longer-term dynamic memories. e.g., the trapped flux state between 2.8 and 3.2 ns is unchanged even as the excitation is varied. A neural network model presented below describes the time-dependent memory and the respective information flow patterns that are useful to analyze experimental results on the 4-loop network.

Network Model

The simulation model evaluates a disordered network as an ideal deterministic system so that the time-dependent behavior can be accurately predicted when the details of the network parameters and the inputs are all known. Simulations of such nonlinear dynamic systems of any size require large computational resources. Experimentally, the time evolution of the disordered network can be precisely observed by tracking each discrete state transition through measurements of single fluxon movement through all the junctions in the network with a sufficiently high time resolution comparable to the scale of Josephson plasma frequency. These observations become impractical, particularly in large networks that are of interest for modeling complex systems leading toward the scale of our brain, as they are subject to continuously varying input information. Alternatively, the approximate state of the network within the vicinity of neighboring states can be experimentally observed statistically from correlations between junction activities obtained by accumulating their firing statistics over a time period T as discussed previously using simulation examples. This analysis is practically feasible and robust to noise as only some of the junctions are required to be monitored with a longer-time integration. The junctions that can be measured for flux flow and the time resolution (i.e., T) limit the precision of the observed states. This leads to multiple possible outcomes with different probabilities for similar inputs. The network model below describes the statistical relation between the input signals (including feedback), the memory state evolution, and the output signals within these limits.

In a continuously evolving network, input information is translated from external excitations to flux flow through some (input) junctions that then interact with the previously trapped flux and generate a flow across other (eventually output) junctions in the network. The dynamics of flux entering or leaving the junctions on the outer-most loop of the network (the input or the output junctions) is related to the state evolution in configuration space that can be described by conserving the number of fluxons in the network. This relation is described below by equating the rate of change of total trapped flux in the loop network (i.e., in all the sub-loops k′ contained within an outer loop k) observed at intervals of integration period T, to the flux flow rates measured across the junctions on the outer loop k during this period:[1] ∑k′⊆k(nk′(T)−nk′(0))=12π∫0T∑ikϕik(t)˙dt.

Flux flow is characterized by the voltage across the junctions, but the Josephson phase description is used since the above relation is also true for quantum fluctuations between flux states. The interaction strengths or weights between a pair of junctions on the outer loop at any time can be estimated from their statistical correlation as the probability of firing of one (output) junction due to a fluxon induced at another (input), i.e., pIO=∫0TVOdt∫0TVIdt. When the network is in a dynamic equilibrium, flux configuration is constrained to within its neighboring states and the change in the trapped flux is negligible (the left side of the above equation ≈0). This is observed in simulation examples in Fig. 3D–G and also shown in experimental results discussed later. The firing probability pIO is related to the state (energy) but is independent of time or the additional number of fluxons considered in the statistical accumulation during T. The stable flux configuration in this steady state directs the incoming flux through different pathways to the outputs as it fluctuates locally at shorter timescales than T, to generate flux (information) flow patterns of different relative strengths defined in the network topology as shown in Fig. 4A. The minimum required number of junctions to be monitored to estimate the approximate state from the flow patterns is therefore equal to the number of junctions on the outer loops when [Eq. 1] is applied to it and the loops it contains. Detailed flow patterns can be mapped across the entire network topology, by monitoring the junctions in the other loops where multiple flux flow paths cross each other.

Fig. 4. (A) Graphical illustration of the flux flow patterns between junctions on the outer loop that are defined by the underlying flux configurations that are stable over the integration period. The flux configurations are directly related to the relative junction firing probabilities. (B) A representation of distinct circulating current (trapped flux) paths in the 4-loop network that contribute to the switching activity of the output junction J6. Current paths from 1 to 4 show the contributions from the 3-loop memory while paths 5 and 6 show directly coupled paths between the input junction J1 and the output junction J6. Flux configuration of the network is composed of combinations of flux from all the paths. (C) Input and output flux streams can be measured as voltages V1 (input), V2, and V3 (outputs) across J1, J3, and J6, respectively. Output spiking probabilities are constrained by the sum rule p1+p2≈1 for long integration times in a steady state. (D) Correlation between activities of any pair of junctions can be described as due to transitions between directly coupled closed loops including both the junctions, together with transitions in closed loops encompassing either of the junctions. For the 4-loop network, the indirectly coupled loops encompassing only J6 refer to the trapped flux in the 3-loop memory network shown in Fig. 4B (1 to 4). All the loops are coupled to each other and p1 is affected by transitions between any of these paths. (E) Network topology representation of the 4 disordered loops with nodes referring to Josephson junctions highlighting the dominant flux flow pathways. Two distinct information flow paths (p1 and p2) exist from input to outputs.

In a more general case of neural network operation, input excitations are continuously changing with time, the flux configuration evolves as shown in Fig. 3H and the network deviates from dynamic equilibrium. The rate of change of trapped flux in the network is finite and sometimes comparable to the flow rate of flux on the outer-loop junctions. The firing probabilities are no longer independent of time or the number of state transitions considered. For any choice of T defined by an external clock, the states that are stable during this period resemble statistics during dynamic equilibrium with steady flux flow patterns. Short-to-long-term memories can therefore be defined relative to this clock, e.g., when the rate of change of trapped flux in loops is very small compared to the flow rate of outer-loop junctions, the states represent a long-term dynamic memory. States that are stable for a shorter duration relative to T (i.e., when the rate of trapped flux change is comparable to the output flow rates) cannot be identified from the junctions on this outer loop alone and require additional junctions with suitable shorter integration periods. The timescales of flux motion between loops that limit the precision of the evolving flow patterns depend on the junctions in the path and their interactions with other junctions or external excitations.

When a fluxon is induced with local excitation across a junction, it can flow through one of the possible pathways in the network while perturbing previously trapped flux with discrete state transitions. Any state transition is accompanied by a firing event across one or several junctions in that path. Each junction can be described to be in one of the two distinct states in time (i.e., 1: at a finite voltage and firing when I>Ikj and 0: at zero voltage or idle when I<Ikj) depending on if the current I through it is either above or below its critical current Ikj for junction j in loop k. Considering the junctions as neuron-like firing elements in the network allows description of associative memories such as in modern neural networks (6). Each junction is coupled to all other junctions in the network through one or more closed-loop pathways that can encompass circulating currents from trapped flux. The firing rate of a junction j is proportional to the current I through it given by the cumulative circulating currents from trapped flux nk(t) in all the loops k through it in addition to the currents from external excitations (x) as[2] Ij=∑k∈jnk(t)Φ0Lk+∑xcxj(t)Ix(t),

cxj describes a fraction of the excitation Ix through the junction j that depends on the network parameters, the flux configuration nk(t), and also the excitation currents Ix(t). While the loop network is a discrete state system, it emulates an analog neural network due to statistical analysis. Even a qualitative description of a junction state as in Eq. 2 is useful to estimate the temporal stability of a particular dynamic memory localized to a small category of states. The rate of change of trapped flux in any loop k is given by the derivative of nk(t). Therefore, when the flow patterns are stable, nk(t) is unchanged (from Eq. 1) and the flow rate through the junctions follows the excitations Ix(t) (from Eq. 2) to produce constant firing probabilities pIO. The temporal stability of any single flux flow pathway in the network depends on the relative dynamics of the first and the second terms of Eq. 2 to the current I through that junction in time intervals of T. For stable excitations, nk(t) converges to a constant value on the timescale of the loop parameters for charging (or discharging) of flux related to their time constants LkRk. Lk is the inductance of loop k and Rk is the equivalent resistance of the junctions in the normal state in that loop. Since Lk∝lkwk, where lk is the length and wk is the width of the loop branch, larger loops with longer-range spatial interaction between junctions have higher time constants. They can store longer-term memories relative to smaller loops. The temporal stabilities of different pathways of flux relative to each other depend on their length scales of interaction requiring an appropriately variable T to capture them. A complete evaluation of the time-dependent flow patterns can therefore be performed using frequency spectrum analysis of the junction flux responses relative to the excitation spectrum.

Computations can be performed with controlled excitations and generated responses. Complex computations can be separated into a series of information write and read operations in time. A read operation can be performed with suitable input excitations when the memory state is weakly perturbed within its vicinity during this time period (or the left side of Eq. 1≈0). Since the observed firing probability is independent of time and the number of firing activities accumulated during this period, a sensitive read operation can be performed with small excitations to weakly perturb the trapped flux state over long integration times. Similarly, a write operation can be performed with input flux flow to induce large changes in the flux configuration relative to feedback currents. Flow at a faster rate relative to the time constant can perform an effective write.

The model can be used to map the flux flow patterns in a dynamically evolving network in response to unsupervised excitations or to perform supervised computations by programming the network through controlled feedback as described in the 4-loop network. The firing activity of the output junction J6 (Fig. 4C) is dependent on the current through it from the trapped flux in loop pathways encompassing it, some examples of which are shown in Fig. 4B, together with the excitation currents. Loops that are labeled 1 to 4 in Fig. 4B contribute to trapped flux in the three memory loops while the examples in Figs. 4B5 and 4B6 show trapped flux current paths that encompass both the input J1 and the output J6. Assigning positive sign for flux flow from left to right aligned with clockwise circulating currents, an input flux stream induced at J1 can flow toward the other junctions on the outer loop J3 and J6 that can be observed from relative firing probabilities given by p1=∫0TV2dt∫0TV1dt and p2=∫0TV3dt∫0TV1dt as shown in Fig. 4C. In a steady flow state with stable flux paths, [Eq. 1] can be rewritten as p1+p2=1 and the flow through J6 is statistically equivalent to flow through J4 and J5. Similarly, the flow through J2 and J3 are equivalent.

Flux flow strength and its temporal stability between two junctions (i.e., J1 and J6 here) depend on direct interactions between the two junctions through closed-loop pathways in common along with indirect interactions between either of the junctions and the remaining loops or excitations in the network. The input J1 and the output J6 are directly coupled through common closed loops (e.g., Fig. 4B loops 5, 6) and indirectly coupled through other closed loops (e.g., Fig. 4B current loops 1 to 4 for junction J6) that are schematically shown in Fig. 4D. Each of these pathways represents different length and hence timescales for interactions. The state of the junction J6 and its firing probability p1 with respect to input at J1 is strongly dependent on the indirectly coupled (memory) loops (k=1,2,3) due to their larger size and flux trapping capacity. Excitation current I1 is strongly coupled to input junction J1 and I2 to output junction J6 due to their short interaction length scales. I2 can be used as feedback to guide the flux flow from input to either of the paths p1 or p2. The probability of a fluxon induced at J1 through the path p1 is higher if the flux is concentrated in loops 2 and 3 relative to loop 1 and aligned to the flow direction. When the excitations are changing at long timescales as in the experiments reported, the temporal stabilities of the flow patterns are dependent only on the loop geometries. For weak excitations i.e., at low currents at I1 and I2 relative to the circulating current from trapped flux, the flux trapped in the 3 memory loops strengthens or weakens the path p1 for output J6. p1 can be programmed by filling the loops (both directly and indirectly coupled) up to a desired memory configuration with large local current excitations over short timescales relative to the time constants of the loops. Flow patterns can be mapped along different pathways with time-dependent strengths p1(t) and p2(t) in the network topology as shown in Fig. 4E from statistical correlations between junction firing activities using a qualitative analysis of the network topology. This model can be extended to other (dynamic) disordered systems with the translation between length scales of interactions and temporal stability of the memories defined by the physical parameters.

Experimental Measurements

The experimental data on the 4-loop network with time-dependent local excitations can be understood from this model. A fabricated YBCO loop network is shown in Fig. 1E. Measurements were made in a range of temperatures and we show the data collected at 28 K. The details of methods used for fabrication and characterization of the network including additional experimental results that are consistent with the model described in this paper are reported in ref. 7. A current excitation I1 is used to induce an input flow V1 across J1 and the output flow is measured as V2 across J2, which can be modified using feedback current I2 (Fig. 1E). The excitations are varied over long timescales with sinusoidally varying I1 and I2 at a few Hz to kHz which are several orders of magnitude slower than the Josephson plasma frequencies ranging from 109 to 1012 Hz that define the flow rate. Therefore, the network is in a dynamic steady state at each data point in all the presented data. The output firing probability is given by p1=V2V1. To scan the parameter space, the excitation currents I1 and I2 are systematically varied over a range (I1 in ±1 mA, I2 in ±100μA). The flux flow probabilities p1 are shown for different input and output flux flow rates (V1 and V2) in Fig. 5A with an experimental integration time of T≈1 ms per data point. p1 can be varied smoothly at these timescales by varying inputs I1 and I2 but any particular value of p1 can be achieved at multiple different combinations of I1 and I2. When one of the currents is large with respect to the other, p1 is the same along the line that represents a constant value of I1I2. This is not true when the values of I1 and I2 are close, causing oscillations in p1 along the transition between C1 and C8 (or between C4 and C5) in Fig. 5A.

Fig. 5. (A) Flow strength in pathway p1 in the range −1 to 1 as a function of excitation currents I1 and I2. Transitions between flow patterns are disordered reflecting the network structure. Oscillations represent differences in temporal stabilities of states discussed in the text. (B) Flow strength in pathway p1 in the range −1 to 1 across the measurement state-space of V1 and V2 shows 8 flow patterns and the transition edges between them. (C) Simulation results of the energy of the state with active excitations across the state space of V1 and V2. Energy monotonically increases with V1 and V2 due to the modulation from increasing excitation currents. (D) Simulation results of the relaxed state energy after excitations are turned off, plotted against the state space of V1 and V2 while they were active in a steady state. Saturated states in orange represent the saturated trapped flux in the 3-loop memory also shown in Fig. 6A. Useful memory states store an output firing probability in the range 0>p1>1. (E) Flux flow patterns between input and output are mapped from the experimental results on the 4-loop YBCO network (Fig. 1E) for a long integration time of T=1 ms. Different flow patterns characterize the memory state categories separated either by a change in flow direction along a path or large changes in relative strengths.

As I1 and I2 are varied to traverse the state space, transitions can occur between states separated by larger energy gaps as shown in the simulations. These transitions are characteristic of large changes to the flux configuration causing dynamically changing flux flow patterns in the network. 8 regions (labeled C1 to C8) represent different stable flux flow patterns through the network. They can be understood by re-plotting the firing probabilities in the state space of input and output flow rates (i.e., V1 and V2) as shown in Fig. 5B. The state transitions occur either along V1, V2=0 or along V1=V2 (i.e., p1=1). These transitions are therefore a result of a change in the direction of flux flow (caused by changes in trapped flux) through one or more loops as mapped in Fig. 5E using Fig. 5B. Transitions from C8 to C1 or C4 to C5 occur due to a switch in the flux flow direction through J1 or at V1=0. Similarly transitions between C2 and C3 or C6 and C7 occur due to flux flow direction change at J6 or at V2=0. The flux direction changes between C1 and C2 or C5 and C6 are at J3 along V1=V2.

The experiments can be reproduced in simulations to map the energy of the state across the state-space of V1-V2 as shown in Fig. 5C when excitation is on and in Fig. 5D when excitation is off to let the network completely relax. State transitions between 8 different flow patterns are characterized by contours in energy modeling the experiments. The state categories seen in Fig. 4C overlap with these different flow patterns. The energy of the state is modulated by the excitations in comparison to trapped flux currents due to their relatively larger contributions and is therefore increasing for higher currents (or flux flow rates) observed in Fig. 5C. When allowed to relax, the underlying stable flux configurations can be realized in Fig. 5D. It is clear from both experimental and simulation results that when the network is in a steady flow state, the observed firing probabilities in a pathway depend on the flux configuration but are independent of time or the number of firing events, consistent with the network model. The interaction strengths between junctions are independent of the scale of applied excitations or time. Exceptions to these are when the excitation strengths (i.e., flow rates) or timescales match the timescales of flux movement through loops (defined by interaction length scales). States can be programmed at these timescales. Any state can be accessed with multiple different combinations of input flux flows and feedback currents as shown in Fig. 5D but over different timescales.

From Fig. 5D, most of the stable states are confined to a narrow band in regions C1 and C5 for flux input from J1 distributed along the 3 memory loops. Both regions span the region of p1 between 0 and 1 but the flow in C5 is inverted in direction. The transition from C8 to C1 is due to the change in the circulating current direction in both loops 2 and 3 combined. Similarly, C1 to C2 is due to a flip in the trapped flux direction in loop 1. The remaining 6 regions represent a single saturated flux state that is stable over the entire region. The differences in the flow rates in these regions are due to the relative changes between I1 and I2. These saturated states on either extreme bound the range of information in excitations that can be stored in unique states. Therefore, they can be accessed with large values of one excitation (I1) with respect to another (I2) or when I1I2≈0 or ∞. Along a line of a specific value of I1I2 at the transition between C1 and C8, the strength of the pathway p1 is not independent of the scale of excitations but instead shows oscillations. These oscillations span all the flux configurations that cause the flow direction switch through loops 2 and 3. At low excitation currents, the contribution to switching from the circulating currents is larger (from Eq. 2). Therefore, the oscillation width (equivalent to wavelength) appears to decrease at higher currents. There are significant differences in the interaction length scales between J1-J3 and J1-J6 due to disorder. Therefore, the oscillations seen relative to a linearly changing excitation are a direct consequence of differences in the temporal stabilities of these different states.

Dynamically stable (or unstable) states can be experimentally mapped by observing the changes in the relative flow rates observed at the output given by dV2dV1. When the network deviates from a steady flow state, the rate of change in trapped flux can be comparable in value to the observed flow rates. At long experimental integration times of T≈1 ms, most of the states are unstable and the observed differences in dV2dV1 are negligibly small. However, the states that cause changes to trapped flux over long length scales are noticeable due to the long timescales involved in flux movement. I1 is varied between ±1 mA for different constant bias currents at I2 and the relative rate of change of flow in the flux path given by dp1dt equivalent to dV2dV1 is shown for different input flows V1 in Fig. 6A. Flow patterns from C1 to C5 can be mapped to their respective relative flow rates (dV2dV1) for different feedback currents at I2. Even during seemingly abrupt transitions as observed in the state space of V1 and V2 (Fig. 5B), the transition regions are stable over a wide range of V1, labeled as C12 between C1 to C2 and as C45 between C4 and C5 as shown in Fig. 6A. The changes to the trapped flux in the three loops and the current paths induced by I1 and I2 during these transitions are mapped in Fig. 6B. First, since C1 comprises most of the stable trapped flux states with comparable short-term stabilities, a constant relative rate of flow change dV2dV1 is expected. A transition from C1 to C2 involves flipping of trapped flux direction in loop 1 as shown in Fig. 6B, during which the flow pattern is stable at p1=1. The region C12 spans all the states of loop 1 from filled to empty while flux in loops 2 and 3 is unchanged. When the excitation changes at a much slower rate compared to the timescales of transition through these states, the width of region C12 defines the timescale of stability of the flow pattern p1=1 equivalent to the length scale of loop 1. Similarly, the width of C45 is proportional to the length scales of loops 2 and 3 combined. The flow pattern of p1=0 (C45) is a longer-term memory compared to the flow pattern of p1=1 consistent with the differences in loop inductances. At smaller excitation currents, the changes to the trapped flux are larger relative to the changes to input flow V1. This is also seen in larger oscillations at lower currents in Fig. 5A. The temporal stabilities of these transition states, therefore, change relative to excitation currents and can be programmed.

Fig. 6. (A) Measurements of the relative rate of change of flow between J1 and J6 given by dV2dV1 for linearly varying input flow V1 with different constant excitations at I2. The transition regions of C12 and C45 show evidence of time-dependent memory with the transition widths proportional to their temporal stabilities. (B) Trapped flux configurations and the current flow directions for states transitioning from flux flow paths C1 to C2 and C4 to C5. Transition extends over the region of states where the moment of trapped flux in loop 1 is inverted from C1 to C2. Similarly, inversion in flux direction in loops 2 and 3 causes the transition from C4 to C5.

Summary

We report an experimental implementation of a disordered network of superconducting loops to demonstrate neural network behavior in its time evolution. Input information from flux flow excitations is translated to discretely evolving trapped flux memory states that guide the flowing flux along different pathways in the network topology. These resemble the information flow pathways formed in our brain with each having a different duration of dynamic stability. Computational properties such as associative memory and categorization of states are observed that can be programmed using feedback. This behavior is repeated in simulations performed on an equivalent circuit model. A network model is described to evaluate the results using statistical analysis of flux flow across junctions accumulated in intervals of an integration period T. This analysis can approximately estimate the state of the network from a qualitative analysis of the network topology and the loop geometries. However, it is a practical approach to perform computation operations on a large-scale network with a few local excitations by monitoring only the junctions on some outer loops. Analysis of time-dependent behavior of the experimental results using this model shows strong evidence of time-dependent (short and long-term) memories. Distinguishing features of these networks that result in this behavior are the disorder, superconductivity for macroscopic coherence (or long-range interactions), and the nature of Josephson junctions that translates excitations to oscillations in phase to cause flux flow. Variation of length scales with a distribution of small to large loops across a large-scale network as described in ref. 8, provides a brain-like architecture to perform neuromorphic computing and simulate collective behaviors of complex networks.

This work was supported as part of Quantum Materials for Energy Efficient Neuromorphic Computing, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0019273. We acknowledge the experimental work at University of California, Riverside supported by Department of Energy National Nuclear Security Agency Grant DE-NA0004106 and AirForce Office of Scientific Research grant number FA9550-20-1-0144. U.S.G and R.C.D thank Prof. Henry Abarbanel for helpful discussions. We thank the help in experimental work by Han Cai and Jay LeFebvre at University of California, Riverside.

Author contributions

U.S.G. and R.C.D. developed the theoretical model; U.S.G., R.C.D., and S.A.C. conducted experimental work and simulations; R.C.D. edited the paper; and U.S.G. wrote the paper.

Competing interests

U.S.G. and R.C.D. are inventors on a U.S. provisional patent application 63/226,743 submitted by University of California, San Diego that covers superconducting Josephson disordered neural networks.

Data, Materials, and Software Availability

All study data are included in the main text.

Reviewers: S.A.K., Stanford University; and T.J.S., Salk Institute for Biological Studies.
==== Refs
1 J. J. Hopfield, Neural networks and physical systems with emergent collective computational abilities. Proc. Natl. Acad. Sci. U.S.A. 79 , 2554–2558 (1982).6953413
2 T. Fulton, R. Dynes, P. Anderson, The flux shuttle-a Josephson junction shift register employing single flux quanta. Proc. IEEE 61 , 28–35 (1973).
3 M. Loponen, R. Dynes, V. Narayanamurti, J. Garno, Measurements of the time-dependent specific heat of amorphous materials. Phys. Rev. B 25 , 1161 (1982).
4 S. A. Cybart , Nano Josephson superconducting tunnel junctions in YBa2Cu3O7-δ directly patterned with a focused helium ion beam. Nat. Nanotechnol. 10 , 598–602 (2015).25915196
5 L. Lundgren, P. Svedlindh, P. Nordblad, O. Beckman, Dynamics of the relaxation-time spectrum in a CuMn spin-glass. Phys. Rev. Lett. 51 , 911 (1983).
6 D. Krotov, J. J. Hopfield, Dense associative memory for pattern recognition. Adv. Neural Inf. Process. Syst. 29 (2016).
7 U. S. Goteti, H. Cai, J. C. LeFebvre, S. A. Cybart, R. C. Dynes, Superconducting disordered neural networks for neuromorphic processing with fluxons. Sci. Adv. 8 , eabn4485 (2022).35452286
8 U. S. Goteti, R. C. Dynes, Superconducting neural networks with disordered Josephson junction array synaptic networks and leaky integrate-and-fire loop neurons. J. Appl. Phys. 129 , 073901 (2021).
