Solitonic Josephson-based meminductive systems

Memristors, memcapacitors, and meminductors, collectively called memelements, represent an innovative generation of circuit elements whose properties depend on the state and history of the system. The hysteretic behavior of one of their constituent variables, under the effect of an external time-dependent perturbation, is their distinctive fingerprint. In turn, this feature endows them with the ability to both store and process information on the same physical location, a property that is expected to benefit many applications ranging from unconventional computing to adaptive electronics to robotics, to name just a few. For all these types of applications, it is important to find appropriate memelements that combine a wide range of memory states (multi-state memory), long memory retention times, and protection against unavoidable noise. Although several physical systems belong to the general class of memelements, few of them combine all of these important physical features in a single component. Her we demonstrate theoretically a superconducting memory structure based on solitonic long Josephson junctions (LJJs). We show that the Josephson critical current of the junction behaves hysteretically as an external magnetic field is properly swept. According to the hysteretic path displayed by the critical current, a LJJ can be used as a multi-state memory, with a controllable number of available states. In addition, since solitons are at the core of its operation, this system provides an intrinsic topological protection against external perturbations. Solitonic Josephson-based memelements may find applications as memories, and in other emerging areas such as memcomputing, i.e., computing directly in/by the memory.

Memristors, memcapacitors, and meminductors, collectively called memelements, represent an innovative generation of circuit elements whose properties depend on the state and history of the system [1].The hysteretic behavior of one of their constituent variables, under the effect of an external time-dependent perturbation, is their distinctive fingerprint.In turn, this feature endows them with the ability to both store and process information on the same physical location, a property that is expected to benefit many applications ranging from unconventional computing to adaptive electronics to robotics, to name just a few [2,3].For all these types of applications, it is important to find appropriate memelements that combine a wide range of memory states (multi-state memory), long memory retention times, and protection against unavoidable noise.Although several physical systems belong to the general class of memelements, few of them combine all of these important physical features in a single component.Here, we demonstrate theoretically a superconducting memory structure based on solitonic long Josephson junctions (LJJs).We show that the Josephson critical current of the junction behaves hysteretically as an external magnetic field is properly swept.According to the hysteretic path displayed by the critical current, a LJJ can be used as a multi-state memory, with a controllable number of available states.In addition, since solitons are at the core of its operation, this system provides an intrinsic topological protection against external perturbations.Solitonic Josephson-based memelements may find applications as memories, and in other emerging areas such as memcomputing [3,4], i.e., computing directly in/by the memory.
Circuit elements, specifically, resistors, capacitors, and inductors with memory [5][6][7][8][9], i.e., elements with characteristics that depend on the past states through which the system has evolved, have recently received increasing attention.Beyond the obvious applications in storing information, these elements can be combined in complex circuits to perform logic [10] and unconventional computing operations [2,4,[11][12][13][14] in massive parallel schemes [3], and in the same physical location where storing occurs.Superconducting circuits that store and manipulate information are particularly appealing in view of their low-energy operation.Among these, a superconducting tunnel junction-based memristor was recently suggested [15,16].However, this type of element does not feature controllable multiple states that can be easily pro-tected against unavoidable noise, due to a stochastic drift of the memory [14].
Our proposal instead is based on a long rectangular tunnel Josephson junction (LJJ) subject to a suitable periodical driving.A tunnel Josephson junction is a quantum device formed by sandwiching a thin insulating layer between two superconducting electrodes, and "long" refers to the physical length of the junction (L) which is supposed to exceed the Josephson penetration depth (λ J ).A scheme of a LJJ with an in-plane magnetic field (H ext ) is shown in Figure 1a.A LJJ is the prototypical system to investigate solitons in a fully solid-state environment, and the history-dependent behavior that we envision stems from how solitons rearrange their configuration along the junction under the effect of an external magnetic field.
The phase dynamics of a LJJ is described by the sine-Gordon equation [17][18][19]: Above, ϕ is the macroscopic quantum phase difference between the superconductors, α denotes the intensity of the damping effect, x is the spatial coordinate along the junction, and t is the time (see SI).The boundary conditions of equation (1) read where H(t) is the normalized time-dependent external magnetic field, and L = L/λ J is the normalized length of the junction.By varying H(t), the phase ϕ evolves according to equations ( 1) and ( 2).For a spatially homogeneous supercurrent density, the Josephson critical current I m s (t) of the junction shows a "Fraunhofer-like" diffraction pattern consisting of overlapping lobes as the magnetic field is increased, and described by the following equation [20][21][22]: where I c is the zero-field, zero-temperature junction critical current.This behavior is shown in Figure 2a as the driving magnetic field is swept "forward" from zero.A diffraction lobe corresponds to a specific number of solitons present along the junction [22].When the external magnetic field penetrates the junction edges it induces Josephson vortices .The length and the width of the junction are L > λ J and W λ J , respectively, where λ J is the Josephson penetration depth.A LJJ excited by a magnetic flux falls into the category of field-controlled meminductive systems, since the input and output variables are the applied magnetic field and the Josephson critical current, respectively.The symbol used to represent the solitonic Josephson-based meminductive system (SJMS) is shown.Fluxons (Φ 0 ) within the junction surrounded by supercurrent loops are also represented.b, Schematic of a possible memory drive formed by an ensemble of SJMSs.The core of the device is a LJJ excited by an in-plane magnetic field, with specific read-out electronics for the critical current.As an example, we display here a junction with length L = L/λ J = 10 by which a 4-state memory element can be defined.These distinct states are labelled by the number of solitons arranged along the junction.The peaks in the dϕ/dx curves (see SI) and the number of loops of Josephson current surrounding the fluxons are indicated as well.
along the weak-link, according to the nonlinearity of equation (1).These vortices, i.e., solitons, are induced by persistent supercurrent loops carrying a quantum of magnetic flux, Φ 0 .The critical current, and the resulting patterns as the driving field is swept, are the physical quantities on which we focus since they can be measured with conventional techniques.In all forthcoming calculations we use parameters typical of Nb/AlOx/Nb tunnel junctions as the ideal materials combination to implement solitonic Josephson-based meminductive structures.
Figure 2b shows the diffraction pattern of the critical current when the magnetic field direction is reversed.The resulting "backward" diffraction pattern markedly differs from the forward pattern shown in Figure 2a.For a given magnetic field H, the current state in which the system is found depends on the field history.This is a remarkable feature of the dissipative solitonic dynamics described by equation (1).Different current states correspond to different numbers of solitons arranged along the junction, and the transition from a diffraction lobe to another corresponds to the injection, or the ejection, of solitons [22].As in any dissipative dynamics, the state of the system is not only determined by the value of the drive but it also depends on the path followed by the system.This induces the forward-backward asymmetry, and the hysteretic diffraction patterns shown in Figures 2a,b.
Figures 2b-d display the forward-backward diffraction patterns as a function of the junction length.Specifically, by in-creasing the length, the number of lobes forming the pattern grows, and the hysteretic asymmetry between forward and backward patterns is enhanced.Notably, L can be tuned as well by changing the junction operation temperature (T ) owing to the temperature dependence of λ J (T ) (see SI).
The presence of both the hysteretic behavior of the critical current and highly-distinguishable current states suggests possible applications of the LJJ.For instance, this device can be used as a field-controlled memelement [1,5,9], in which the time-dependent input/output related variables are the external magnetic field H(t) and the Josephson critical current I m s (t), respectively.We envisage here a memelement with distinct memory states which make use of the lobes of the forward/backward diffraction patterns.For a given applied magnetic field, the memelement state is determined by the value of the critical current, the latter keeping track of the field history, and pointing to a specific number of solitons present in the junction.Since the critical supercurrent and the magnetic field are the variables yielding the history-dependent behavior, our junction can be regarded as a meminductive system [1,23], specifically, a field-controlled solitonic Josephson-based meminductive system (SJMS).
More generally, the LJJ can be thought as a multi-state memory in which each memory state is represented by a specific diffraction lobe, and labeled by the number of excited solitons (see Figure 1b).For example, by referring to the diffraction patterns shown in Figures 2a,b, three backward lobes can be easily recognized within the range H ∈ [0, 2] in clear contrast to one single forward lobe, by which a 4-state memory could be built.
On general grounds, a good memelement has to read/write in short times, and has to be sufficiently robust against external fluctuations (noise) that tend to destroy the stored information.On the one hand, reading the state of the SJMS, namely the critical current I m s (H), can be performed by conventional and well-established techniques without altering the memelement state.On the other hand, the writing process of each memory state depends on the operating frequency (ω H ) of the magnetic field, and on the ability of the system to follow a fast periodic driving.To quantify the LJJ memdevice performance as the driving frequency and the temperature are changed we make use of a figure of merit defined by the difference between the forward and backward critical currents, , where H i is the magnetic field at the midpoint of the i-th backward diffraction lobe, as shown in Figure 3a for i = 1, 2, 3.For large δ I i one can safely distinguish distinct memory states, namely, the current states.Furthermore, to further characterize our memdevice we have included a Gaussian thermal fluctuation term in equation ( 1) (see SI) thereby making the SJMS a stochastic memory element [1,14,15,[24][25][26] whereas a noiseless driving field source was considered.The relevant supercurrent differences (δ I i ) are then calculated between the averaged diffraction patterns.
Figure 3b shows δ I i as a function of the driving frequency ω H , for T = 1.2 K.The memory states defined in Figure 3a are stable up to a driving frequency ω H ∼ 0.5GHz.At higher frequencies, i.e., for ω H 1GHz, the system is not able to respond anymore to the fast driving.In this region of frequencies, δ I i tends to increase (see SI), the diffraction patterns are not stable, and therefore cannot be used to safely distinguish the memory states.
As expected, due to its topological nature the LJJ memory s (H i ), computed by averaging over N exp = 100 numerical realizations of the Josephson critical current, as a function of the driving frequency ω H for T = 1.2 K.The memory states are stable up to ω H ∼ 0.5GHz.At higher frequencies, i.e., ω H 1GHz, the system is no more able to respond to the fast driving.
shows remarkable robustness against thermal disturbances: being a soliton-based memelement, it is intrinsically protected against small fluctuations.Indeed, the states of the memory are associated to the number of solitons present in the LJJ and, therefore, are quantized [22].The creation of a soliton is a macroscopic quantum phenomenon involving crossing of a potential barrier [22].Far away from the superconducting critical temperature (T c ), the presence of an energy barrier in a damped dynamics prevents noise-induced state degradations, i.e., the so-called "stochastic catastrophe" [14].thermal fluctuations, as the driving frequency is set to ω H ∼ 0.04GHz.Specifically, here we show how the temperature affects the forward (Figure 4a) and backward (Figure 4b) diffraction patterns.In particular, by increasing the temperature leads to a smoothing of the interference patterns with broadened transitions between lobes due to noise-induced creation or destruction of solitons.Nevertheless, the memory states tend to degrade only for somewhat high temperatures approaching T c (see the results for T > 4.2 K in Figure 4a and  b).
Finally, the stability of our Josephson-based memory as the temperature is changed is quantified in Figure 4c.In particular, the memory states turn out to be stable against large temperature variations, i.e., δ I i is roughly constant as long as T 4.2 K.For higher temperatures, the average forward/backward diffraction patterns tend to superimpose so that δ I i vanishes with the following suppression of the memory states at the critical temperature.
In summary, we have suggested long Josephson junctions excited by an external magnetic field as prototypical multistate superconducting memories.Our proposal for a memory element is based on the characteristic hysteretic behavior of the critical supercurrent as the driving field is swept.The resulting memelement realizes a multi-state memory with a number of states controllable via the effective length of the junction.The solitonic nature at the origin of the critical current hysteresis makes these memory states stable and robust against thermal fluctuations.Our memory scheme represent the first endeavor to combine superconductivity and solitons physics in one single memelement, and could find potential application in various emerging areas such as logic in memory and unconventional computing [3,4].

THE SINE-GORDON EQUATION AND ITS SOLUTIONS
In Fig. 5(a), a long and narrow Nb/AlO/Nb Josephson junction (JJ) is represented.The electrodynamics of a long JJ (LJJ) is usually described by a partial differential equation for the order parameter ϕ, namely, the phase difference between the wavefunctions describing the carriers in the superconducting electrodes.In normalized units, the perturbed sine-Gordon (SG) equation reads [1-5] with boundary conditions taking into account the normalized external magnetic field In equation ( 4), space and time variables are normalized to the Josephson penetration depth λ J and the inverse of the Josephson plasma frequency ω p , respectively.They read where R and C are the total resistance and capacitance of the JJ, Φ 0 = h/2e 2.067 × 10 −15 Wb is the magnetic flux quantum, µ 0 is the vacuum permeability, I c and J c = I c /A are the critical current and the critical current area density (A being the junction area).Moreover, t d = λ 1 + λ 2 + d is the effective magnetic thickness, λ i being the London penetration depth of the superconductor S i and d the interlayer thickness.If λ i exceeds the thickness t i of the i-th superconductor, the effective magnetic thickness has to be replaced by td = λ 1 tanh (t 1 /2λ 1 ) + λ 2 tanh (t 2 /2λ 2 ) + d.
The Josephson penetration depth represents the length scale of the system, so that a JJ is regarded as long and narrow if the length and the width of the junction are L > λ J and W λ J , respectively.In normalized unit, the linear dimensions of the junction read L = L/λ J > 1 and W = W/λ J 1. Moreover in equation ( 4), α = (ω p RC) −1 is the damping parameter.
The SG equation admits traveling wave solutions, called solitons [6].In the SG framework, a soliton is often referred to as a kink.For the unperturbed SG equation, i.e., α = 0 in equation ( 4), solitons have the simple analytical expression [2] where the sign ± is the polarity of the soliton (specifically, the minus sign defines an antisoliton) and u is the Swiharts velocity [2], namely, the largest group propagation velocity of the linear electromagnetic waves in long junctions.Specifically, the phase of a soliton (antisoliton) twists from 0 to 2π (from 2π to 0).Alternatively, ϕ/2π has a topological charge +1 for each soliton and -1 for each antisoliton.Moreover, a SG soliton has a well defined physical meaning in the LJJ framework, since it carries a quantum of magnetic flux Φ 0 , induced by a supercurrent loop surrounding it, with the local magnetic field perpendicularly oriented with respect to the junction length [7].Thus, a soliton is usually referred to as a fluxon, or a Josephson vortex, in the context of LJJs.
In equation ( 5), the normalized external magnetic field is , where H ext (t) is the non-normalized external magnetic field lying parallel to a symmetry axes of the junction and along y, see Fig. 5a.For the numerical simulation, we have modeled H(t) as a staircase function formed by steps with "treads" deep ∆t H and "risers" high ∆H [see Fig. 5(b)] and is ramped up from zero to H max , then reduced to −H max and subsequently raised again to zero to perform a double-swept drive.Accordingly, the driving period T H = 4(H max /∆H)∆t H and the frequency ω * H = 1/T H are defined.

THE CRITICAL CURRENT DIFFRACTION PATTERNS
The ϕ-dependent supercurrent as a function of the external magnetic field H can be expressed as where J s (x, y) is the supercurrent density per unit area, and J c (x, y) is the Josephson critical current density.We denote with i c (x) the J c (x, y) integral in the direction of the magnetic field so that the Josephson current becomes In equation (11), ϕ(x) is the phase difference induced by the applied magnetic field H ext .In fact, ϕ depends on the local magnetic field H y (x) through the equations [2] The latter equation comes from the condition W λ J , so that ϕ (x, y) ≡ ϕ (x).
For a short rectangular JJ ( L λ J and W λ J ) the external magnetic field fully penetrates the junction and is spatially homogeneous along it, namely, H y (x) ≡ H ext , so that, according to equation ( 12), the phase is just linearly increasing in the x direction, Accordingly, equation (11) becomes The Josephson critical current is the amplitude of the last integral, that is independent of any phase factor ϕ 0 .By assuming a uniform supercurrent area density within the junction, i.e., J c (x, y) ≡ J c , for 0 ≤ x ≤ L and 0 ≤ y ≤ W, and zero elsewhere, we obtain i c (x) ≡ i c = J c W, according to equation (10).Therefore, equation ( 15) becomes I m s (Φ) where k = 2π µ 0 t d Φ 0 H ext , Φ is the magnetic flux through the effective magnetic area (t d L), and c = i c L = J c LW = J c A.
The long junction case markedly differs with respect to the short case, since both the penetrating external field and the self-field generated by the Josephson current have to be considered, so that ϕ(x) nonlinearly changes along the junction according to equations (4)- (5).Therefore, in normalized units, the maximum value of the Josephson current can be written as [5,8] It only remains to include in equation ( 18) the proper phase difference ϕ(x,t) for a driven LJJ given by solving equations (4)-( 5).The magnetic field dependence of I m s results in "Fraunhofer-like" diffraction patterns [5,9,10].While in the short junction limit [2], different diffraction lobes are well separated, here we observe the overlapping of the lobes.The transitions between these lobes are usually discontinuous.These patterns can be explained in terms of solitons entering the JJ.
Each lobe corresponds to a state with a fixed number of solitons.When the magnetic field increases, the configuration with more solitons is energetically favorable and, thus, the system jumps from a metastable state to a more stable state with more solitons.In the region of H values in which the diffraction lobes overlap, several solutions with different number of solitons may co-exist [9,10].Therefore, the system stays in the present configuration until the following one is energetically more stable.
To further explore the behavior of a magnetically driven LJJ, we have implemented a double-swept drive.The forward, i.e., with H increasing, and the backward, i.e., with H decreasing, patterns are significantly different.For a given value of the magnetic field, the critical currents in the backward and forward evolutions differ and the system is found in a different diffraction lobe.We can associate the forward and backward stable states (at fixed H) with a different number of solitons in the junction.Interestingly, the overall effect is a hysteric behavior in the critical current.
The diffraction patterns of the Josephson critical current, as the driving field is first increased (forward plot) and then reduced (backward plot), are shown in Fig. 6 for several JJ normalized lengths.

THE MULTISTATE STRUCTURE AS A FUNCTION OF THE JUNCTION LENGTH
Forward-backward differences in the hysteretic behavior of the critical current are strongly evident for |H| H c = 2, see Fig. 6.In the forward pattern, the first lobe corresponds to the Meissner state, i.e., zero solitons in the junction, whereas by exceeding the threshold value H c the second lobe begins and solitons in the form of magnetic fluxons penetrate into the junction.This value of the critical field characterizes the diffraction patterns of the Josephson critical current in both overlap and inline LJJs [2,11,12].For H > 0, the backward dynamics is described by N-solitons solutions, with N ≥ 1.The amount of solitons exited depends on both the field intensity and the length of the junction.The many-soliton backward solutions suggest applications of this system as multistate memories, in which each state is clearly indicated by drastic suppressions of the critical current I b s with respect to I f s .To quantify the LJJ memory-device performance we make use of a figure of merit defined by the difference between the critical currents, where H i is the magnetic field at the midpoint of the i-th backward diffraction lobe.For large δ I i one can safely distinguish distinct memory states (MSs), namely, the current states.For instance, by focusing on the panels c,d, and e of Fig. 6, we observe that, in the range H ∈ [0 − H c ], • for L = 6, only one MS is clearly available, with a current difference δ I 1 (H 1 1) ∼ 0.5, see Fig. 6c; • for L = 8, two MSs can be defined, with δ I 1 (H 1 0.5) ∼ 0.5 and δ I 2 (H 2 1.5) ∼ 0.6, see Fig. 6d; • for L = 10, three MSs can be defined, with δ I 1 (H 1 0.32) ∼ 0.4, δ I 2 (H 2 1) ∼ 0.6, and δ I 3 (H 3 1.75) ∼ 0.6, see Figs. 6e.
Finally, junctions with different lengths are characterized by different numbers of distinct available MSs, each of them corresponding to a specific amount of solitons arranged along the junction.
Moreover, for a fixed effective junction length L, the normalized length L(T ) = L/λ J (T ) and, therefore, the amount of MSs of the memdevice can be controlled by changing the temperature T of the system (see below).

PHYSICAL QUANTITIES
To give a realistic estimate of the physical quantities used in the computations, both the superconductors and the insulator making the junction (according to which distinctive values of resistance per area R a and specific capacitance C s of the junction result), and the normalized length of the device have to be chosen.Therefore, let us set a Nb/AlO/Nb junction, characterized by R a = 50 Ω µm 2 = 5 × 10 −11 Ω m 2 and C s = 50 f F µm 2 = 5 × 10 −2 F m 2 , and a length-to-Josephsonpenetration-depth ratio equal to L = 10.In the low temperature regime, the critical current is The main physical quantities to fully describe the system are: • effective magnetic thickness t d = 2λ 0 L + d ∼ 161nm, the London penetration depth of a Nb thin film being λ 0 L ∼ 80nm and setting d = 1nm; • Josephson penetration depth • Linear dimensions L = 10λ J = 60µm and W = 1µm; • Area A = WL = 60µm 2 = 6 10 −11 m 2 ; • Critical current • Plasma frequency • Damping parameter α = 1 ω p RC ∼ 0.24; THERMAL EFFECTS Some of the quantities introduced so far have an explicit dependence on the temperature.In particular, for identical superconductors, [2] • the effective magnetic thickness t d (T ) depends on T through the London penetration depth λ i (T ) = 4 , and • the Josephson critical current I c (T ) depends on T given by the Ambegaokar and Baratoff formula [2] where ∆(T ) is the BCS gap of the superconductors.
Concerning the normalized length, since as the temperature is increased to T * ∼ 0.8T c , being λ J (T * )/λ J (0) ∼ 1.3, the corresponding normalized length becomes L(T * ) ∼ 8.The temperature of the bath influences also the dynamics of the junction.In order to take into account the thermal fluctuations on the phase dynamics, a noise current i T has to be included into the perturbed SG equation The normalized thermal current i T (x,t) is characterized by the well-known statistical properties of a Gaussian random process For a LJJ, the thermal noise amplitude reads [13] The behaviors of the normalized Josephson critical current I c , the plasma frequency ω p , the damping parameter α, the Josephson penetration depth λ J , the normalized length L, and the thermal noise amplitude γ as a function of the temperature, for T c = 9.2K, are shown in Fig. 7.
Moreover, the values of the thermal noise amplitude, the damping parameter, and the plasma frequency in correspondence of few specific temperatures (used to obtain the results shown in the manuscript) are listed in Table I.

THE FREQUENCY RESPONSE
Independently of the physical mechanism defining the state of the device, the memelement response is usually strongly dependent of the frequency of the input drive [14].At low frequencies, the system has enough time to adjust its state to the instant value of the drive, so that the device non-linearly behaves and a hysteretic evolution results.Conversely, at high frequencies, there is not enough time for any change during an oscillation period of the drive.Therefore, we explore the effects of variations of the driving frequency on the behavior of our device.However, our system is an example of a memory that benefits from, and properly work only in, the presence of noise.To discuss this point, we compare results obtained in both the deterministic and stochastic approaches, by taking into account several temperatures.
To quantify the LJJ-based memory performances as the driving frequency is changed, we use the distances δ I i defined in equation (19).Specifically, Fig. 8a shows the midpoint values H i of the backward diffraction lobes within the field range H ∈ [0 − 2], and the distances δ I i (i = 1, 2, 3) for these fields, for a junction with L = 10.
First, we analyse the device performance in absence of thermal noise.The behavior of δ I i (i = 1, 2, 3) as the driving frequency ω H is changed is shown in Fig. 8b for the deterministic case, i.e., no noise source is considered in the model.For the sake of clarity, we define in Fig. 8b two threshold frequencies, ω 1 and ω 2 , and examine the results in different frequency ranges: • for ω H ω 1 the distances δ I i approach constant values, inasmuch steady diffraction patterns are obtained.For these frequencies, the logic states are definitively robust against frequency variations; • in the range ω H ∈ [ω 1 − ω 2 ] the values of δ I i , and accordingly the amount of the logic states, significantly deviate from the steady ones; • for ω H ω 2 the system is not able to respond to extremely high driving oscillations, so that the backward patterns are highly disordered, despite the fact that δ I i → 1, and therefore the logical states cannot be safely distinguished.
As discussed above, realistic devices are subject to thermal noise.Far away from the critical temperature, the addition of the thermal noise has the effect to stabilise the dynamics and, therefore, to access to higher driving frequencies.
FIG. 8. a, Forward and backward diffraction patterns for H ∈ [0 − 2] and L = 10.The backward pattern is composed by three lobes in the place of the large lobe of the forward one.The value of the magnetic field H i in the center of each lobe and the differences δ I i (i = 1, 2, 3), see equation (19), for H ≡ H i are also shown.The values of δ I i are used to check the behavior of the logic states of the LJJ-based memory against frequency variations.b, δ I i (i = 1, 2, 3) as a function of the driving frequency ω H = ω p /T H in absence of thermal noise.As ω H reduces, the diffraction patterns tend to become steady and δ I i approach constant values.Specifically, by defining two threshold values, ω 1 and ω 2 , the behavior of the device in different ranges of frequencies can be discussed: i) for ω H ω 2 the system is not able to respond to extremely high driving frequency oscillations; ii) in the range ω H ∈ [ω 1 − ω 2 ] the memory cannot safely provide three logic states ; iii) for ω H ω 1 the distances δ I i approach constant values, and, in spite of frequency variations, the system provides three distinct states.c, d, e, and f, distances δ I i (i = 1, 2, 3), see equation (25), as a function of ω H for T = {0.02,0.3, 1.2, 4.2}K, respectively.The legend in panel c refers to these panels too.8a for L = 10 and H i = 0.32, 1.0, 1.75 with i = 1, 2, 3 (panels a, b, and c, respectively).The magnetic field value H i is chosen in the midpoint of the i-th backward diffraction lobe, so that the robustness against small field fluctuations is ensured.These graphs are obtained by setting H(t ≥ t i ) = H i , where t i is the time for the magnetic field H(t) to reach the value H i during the backward sweep.The diffraction patterns are computed for T = 1.2K.In spite of the thermal fluctuations, as the magnetic field is set to H(t) = H i , the δ I i values are practically constant in time, i.e., I m s (t ≥ t i ) ∼ I m s (t i ).These results are obtained by setting ∆t H = 4, ∆H = 0.005, and H max = 5, so that the non-normalized driving frequency is ω H = ω p /T H 0.1GHz, where ω p (T = 1.2K) 1.634THz (see Table I).
In the noisy approach, the distances s /I c are computed by averaging over the total number of numerical realizations, N exp = 100, the normalized critical currents as the magnetic field H is swept forward and backward, respectively, when the thermal fluctuations are included in the SG model.
For T = 0.02K, 0.3K, and 1.2K (see Figs. 8c, d, and e, respectively) the values of δ I i are roughly constant and the logic states of the device are definitively stable up to ω H ∼ 0.1GHz.Conversely, for higher frequencies, the inability of the system to adjust its state to rapid changes in the magnetic bias comes to light.
For T = 4.2K, i.e., the liquid helium temperature, the frequency behavior significantly changes, see Fig. 8f.In fact, δ I i (i = 1, 2, 3) approach the values obtained for low temperatures only for ω H ∼ 0.1GHz.Conversely, for lower frequencies the thermal fluctuations have enough time to guide the evolution of the system, so that the state of the system is set by noiseinduced transitions.Therefore, the backward and forward patterns tend to superimpose and δ I i → 0.
Moreover, we verify if the system is able to provide information-storage times longer than any practical reading times, so that it works as a non-volatile memory [23].To this end, we show in Fig. 9 the normalized critical currents I m s /I c and driving field H as a function of the normalized time t for the states defined in Fig. 8a, for T = 1.2K.Specifically, results in panels a, b, and c of Fig. 9 are obtained by freezing the magnetic field to H(t ≥ t i ) = H i with H i = 0.32, 1.0, and 1.75, respectively, t i being the time for the magnetic field to reach the value H i during the backward sweep.In spite of the thermal fluctuations, as the magnetic field is set to H(t) = H i , the critical current is roughly constant, i.e., I m s (t ≥ t i ) ∼ I m s (t i ), so that steady logic states are established.

MEMORY DEVICES AND THE RESPONSE FUNCTION
The properties of a memory-element (memelement for short) depend on the state and the history of the system [24].Specifically, in ideal memristors, memcapacitors and meminductors [24] • (ideal memristor) the resistance depends only on the charge that flows in the system (or on the history of the voltage); • (ideal memcapacitor) the capacitance depends only on the history of the charge stored on its plates (or the history of the voltage across it); • (ideal meminductor) the inductance only depends on the history of the current that flows through it (or the history of the flux).
These definitions can be generalized by invoking a general non-linear, memory-dependent response function g [25] y(t) = g(x, u,t)u(t) where • g(x, u,t) is the response function; • u(t) is the input signal; • y(t) is the output signal; • f (x, u,t) vector function of internal state variables; • x vector of internal state variables.
Generally, in real systems ideal memdevices are usually rare, so that the relation between current and voltage defines a memristive system (i.e., the resistance depends on both the charge and other internal variables of the system), while the relation between charge and voltage specifies a memcapacitive system, and the flux-current relation gives rise to a meminductive system [24].
A distinctive signature of memory devices is the presence of a hysteresis loop in the behavior of the output y(t) and/or the response function g(t) as a function of the input u(t) [24].The features of the hysteresis loop depend on the properties of both the system and the input u(t), such as its amplitude and frequency.Hysteresis loops can be pinched, when the loop passes through the origin (y is zero whenever u is zero and vice versa).Moreover, a pinched hysteresis can be selfcrossing [24], i.e., with the crossing between opposite direction branches of the loop, or not self-crossing.
In contrast with usual memelements defined by equations (26), the behavior of our LJJ-based memory-device is not directly stated in the form of a relation between I m s (t) and H(t) through a response function g(ϕ, H,t).In other words, our memdevice is not described by a current-field expression such as I m s (t) = g(ϕ, H,t)H(t).In fact, in normalized units, the critical current reads The internal state variable of this field-controlled memelement is the phase difference ϕ(x,t), whose dynamics is ruled by equations ( 4)- (5).
We observe that a relation including a response functional comes to light by first-order expanding the cos ϕ(x,t) term in equation ( 27) around the junction edge x = 0, that is by ignoring the non-linearity of the problem, cos ϕ(x,t) ∼ x=0 cos ϕ(0,t) − sin ϕ(0,t) dϕ(x,t) dx 0 x. (28) where the field-dependence of the phase dynamics is stressed.

FIG. 2 .
FIG. 2. Josephson critical current diffraction patterns.a and b, Normalized Josephson critical currents I f s /I c and I b s /I c as the driving field H is swept forward from H=0 to H=5 (right half of panel a), then backward from H=5 to H= − 5 (panel b) and again forward from H= − 5 to H=0 (left half of panel a).The inset in panel a shows one period (T H ) of the driving field.The critical current as a function of H(t) exhibits a diffraction-like pattern formed by lobes which are directly related to the number of solitons arranged along the junction.By sweeping the magnetic field forward and then backward leads to the appearance of a clear hysteretic behavior.This is a distinctive signature of any memdevice.According to this hysteretic behavior, the Josephson junction can be effectively used as a multi-state memory.For any specific range of magnetic field values, each state of the memory is represented by a forward or backward diffraction lobe, labeled by the number of excited solitons present along the junction.c, d, and e, Diffraction patterns for a few junction lengths L. The number of memory states provided by the SJMS can be changed by varying the junction length.The memory states associated with current lobes are indicated with the same notation as in Figure 1b.

FFIG. 3 .
FIG. 3. Frequency response of the memory states.a, Forward and backward diffraction patterns for H ∈ [0, 2] and L = 10.For each backward diffraction lobe, we have considered the middle magnetic field value H i , and calculated the current difference δ I i = I f s (H i ) − I b s (H i ) /I c (i = 1, 2, 3), where I f s (H i ) and I b s (H i ) are the corresponding forward and backward critical currents.b, Difference δ I i (i = 1, 2, 3) between average forward and backward diffraction patterns I f s (H i ) and I bs (H i ), computed by averaging over N exp = 100 numerical realizations of the Josephson critical current, as a function of the driving frequency ω H for T = 1.2 K.The memory states are stable up to ω H ∼ 0.5GHz.At higher frequencies, i.e., ω H 1GHz, the system is no more able to respond to the fast driving.

Figure 4 FIG. 4 .
FIG. 4. Effects of the temperature.a and b, Average forward and backward diffraction patterns I f s /I c and I b s /I c , respectively, calculated for a few temperatures, L = 10, and ω H ∼ 0.04GHz.The patterns are computed by averaging over N exp = 100 numerical realizations of the critical current as the magnetic field is swept forward and backward when thermal fluctuations are taken into account.The legend in panel b refers to both panels.c, Differences δ I i (i = 1, 2, 3) for L = 10 and ω H ∼ 0.04GHz calculated in correspondence of the temperatures set to obtain the results shown in panels a and b.By approaching the superconducting critical temperature (T c 9.2K for a Nb/AlOx/Nb JJ) the forward and backward diffraction patterns tend to superimpose, and δ I i vanishes.

) b FIG. 5 .
FIG.5.a, An Nb/AlO/Nb long Josephson junction (LJJ) in the presence of a homogeneous external magnetic field H ext applied in the y direction.The length and the width of the junction are L > λ J and W λ J (according to the long junction regime), respectively, and A = LW is the junction area, λ J being the Josephson penetration depth.Moreover, t i and d denote the thicknesses of the i-th superconductor and the insulating interlayer, respectively.b, Numerical implementation of the magnetic field drive.Normalized staircase external magnetic field H(t), formed by small steps with height ∆H = 0.005 kept constant for time intervals ∆t H = 10 3 , with H max = 5.In the inset, a driving period T H = 4 × 10 6 of H(t) is shown.The times are normalized with respect to the inverse of the Josephson plasma frequency.

FIG. 7 .
FIG. 7. Normalized Josephson critical current I c , plasma frequency ω p , damping parameter α, Josephson penetration depth λ J , normalized length L, and thermal noise amplitude γ as a function of the temperature, for T c = 9.2K.

FIG. 9 .
FIG. 9. Normalized critical current I ms /I c (left ordinate scale, full symbols) and driving field H (right ordinate scale, solid lines) as a function of the time t, normalized to ω p , for the MSs defined in Fig.8afor L = 10 and H i = 0.32, 1.0, 1.75 with i = 1, 2, 3 (panels a, b, and c, respectively).The magnetic field value H i is chosen in the midpoint of the i-th backward diffraction lobe, so that the robustness against small field fluctuations is ensured.These graphs are obtained by setting H(t ≥ t i ) = H i , where t i is the time for the magnetic field H(t) to reach the value H i during the backward sweep.The diffraction patterns are computed for T = 1.2K.In spite of the thermal fluctuations, as the magnetic field is set to H(t) = H i , the δ I i values are practically constant in time, i.e., I m s (t ≥ t i ) ∼ I m s (t i ).These results are obtained by setting ∆t H = 4, ∆H = 0.005, and H max = 5, so that the non-normalized driving frequency is ω H = ω p /T H 0.1GHz, where ω p (T = 1.2K) 1.634THz (see TableI).
are taken into account.The behaviors of δ I i (i = 1, 2, 3) as a function of ω H for T = {0.02,0.3, 1.2, 4.2}K are shown in panels c, d, e, and f of Fig. 8, respectively.The quantities I f s /I c and I b

TABLE I .
Thermal noise amplitudes, damping parameters, and plasma frequencies in correspondence of few specific temperatures used to obtain the results shown in the manuscript.