Quantum Memristors with Superconducting Circuits

Memristors are resistive elements retaining information of their past dynamics. They have garnered substantial interest due to their potential for representing a paradigm change in electronics, information processing and unconventional computing. Given the advent of quantum technologies, a design for a quantum memristor with superconducting circuits may be envisaged. Along these lines, we introduce such a quantum device whose memristive behavior arises from quasiparticle-induced tunneling when supercurrents are cancelled. For realistic parameters, we find that the relevant hysteretic behavior may be observed using current state-of-the-art measurements of the phase-driven tunneling current. Finally, we develop suitable methods to quantify memory retention in the system.

Circuit elements that intrinsically carry a recollection of their past evolution [1][2][3] promise to bring forth novel architectural solutions in information processing and unconventional computing [4] due to their passive storage capabilities.These history-dependent circuit elements can be both dissipative and non-dissipative, such as memcapacitors and meminductors [2,5], or just dissipative, such as memristors.Classical memristors [6][7][8][9] are elements whose operational definition relates the voltage V and the current I, complemented with an update of one or more internal state variable(s) x carrying information of the electrical history of the system.For a voltagecontrolled memristor The memductance (memory conductance) G depends on both the instantaneous input voltage V and the state variable x, which tracks the past memristor configuration via the update function f .Such dynamics leads to the characteristic pinched hysteresis loops under periodic driving [3,6,7,10,11], a strictly non-linear conductive effect showcasing zero-energy information storage [1].
Even though both the quantization of superconducting circuits [12] and applications of memristors are well established techniques, memristive operation in the realm of quantum dynamics is a largely unexplored area.From an intuitive point of view, the combination of powerful memristive concepts with quantum resources, such as superposition and entanglement, promises groundbreaking advances in information and communication sciences.With this motivation in mind, the idea of a quantum memristor was recently defined in Ref. [13] by introducing * mikel.sanz@ehu.eusthe fundamental components for engineering memristive behavior in quantum systems.However, superconducting circuits naturally include memristive elements in Josephson junctions, a feature exploited in a recently proposed classical superconducting memristor design [14].While this conductance asymmetric superconducting quantum interference device (CA-SQUID) design was able to produce hysteretic behavior [14], it did not include the quantum features of the circuit, including the dissipative origins of the memory or its measurement and quantification.These features are of utmost importance, as the operation of the design is based on quasiparticle tunneling, whose control and measurement have recently seen significant strives forward [15,16].Indeed, to our knowledge, up to now no experimental work has studied the hysteretic IV-characteristics of such systems.In our opinion, the reasons for this are two-fold, namely, 1) the pinched hysteresis loops were only recently predicted to exist for such systems in Ref. [14] with the use of the aforementioned CA-SQUID and a proper selection of parameters, and 2) the experimental apparatus required to control and measure quasiparticle excitations with high accuracy is just beginning to emerge (see Refs. [15,16]).
In this Article, we show that a suitably designed superconducting quantum circuit element with an external phase bias serves as a prototypical quantum memristor via low-energy quasiparticle tunneling.To this end, we describe the device in a fully quantum-mechanical fashion.We apply an ensemble interpretation of the system input and output, while the average superconducting phase difference stores information of the past dynamics.We study the hysteretic signature in a regime achievable with recent quantum nondemolition projective measurements [16], and construct a memory quantifier related to the accumulation of internal state change.Finally, we discuss the quantumness of our proposal, comparing it with Ref. [14].Our proposal represents, to our knowledge, the first design of a superconducting quantum memristor from fundamental principles, exploiting quasiparticle tunneling in memristive quantum information processing.
The envisioned device has the rf SQUID design shown in Fig. 1(a).It consists of a superconducting loop with inductance L, which is interrupted by a dc SQUID with negligible loop inductance acting as an effective fluxtunable Josephson junction.The dc SQUID junctions are made from different materials so that they have the same critical current but a different normal conductance [14].In this way, the effective critical current of the dc SQUID can be completely suppressed by a bias flux of half a flux quantum, Φ 0 /2, threading its loop [14].Finally, we also apply a bias flux Φ d to the rf SQUID loop, resulting in the phase bias The total Hamiltonian of this device is the sum of the system Hamiltonian ĤS , a term for the quasiparticle degree of freedom, and a total tunneling term.The latter includes quasiparticle contributions but, due to the vanishing effective critical current (note that these contributions would yield a renormalization of the qubit frequency in the low-energy regime considered in Ref. [17]), neither pair contributions, nor the Josephson counterterm [17,19,29].Under these conditions, the system Hamiltonian takes the harmonic form where n and φ are the Cooper-pair counting and phase difference operators of the effective junction, respectively.We define the capacitive energy scale E C = 2e 2 /C d with the intrinsic junction capacitance C d and the inductive energy scale E L = (1/L)(Φ 0 /(2π)) 2 .Regarding the dc SQUID, we assume the limit of strong conductance asymmetry needed for the effective junction picture due to the inclusion of quasiparticle excitations (see Supplementary Information).In this limit, the dissipative flow is through the physical junction with a smaller superconducting gap while the junction with a larger gap functions as a shunt for the total Josephson current through the SQUID.Furthermore, we demand that the phase bias is changed adiabatically, i.e., sufficiently slowly to avoid the generation of quasiparticles.Finally, our device operates in the low-energy regime ω 10 , δE 2∆, where ω 10 = √ 2E C E L / is the system transition frequency and δE is the characteristic energy of the quasiparticles above the gap ∆.Even though the Hamiltonian does not warrant operation as a qubit due to the lack of sufficient anharmonicity, the system dynamics is confined to the two lowest eigenstates of Eq. ( 2) when the aforementioned assumptions are complemented with operation in the highfrequency regime ω 10 δE.In this regime, there exist no quasiparticles with sufficiently high energy to excite the system.We emphasize that the slow biasing and high frequency assumptions utilized in this article are not contradictory.The former refers to suppressing unwanted generation of quasiparticles due to the biasing field [20] while the latter refers to a condition on the quasiparticle bath.
The two-level master equation describing the quasiparticle-induced decay takes the Lindblad form [19] ∂ t ρ = −i/ [ ĤS , ρ] + D{ρ} for the system density ρ, with D{ρ} the corresponding Lindbladian dissipator.Note that the master equation assumes adiabatic steering, and employs the Born-Markov and secular approximations.We omit the quasiparticle-induced average frequency shift and the pure dephasing channel.See Supplemental Material for the estimation of these effects.In the low-energy limit, the decay rate factorizes into in the lowest order in ω 10 /∆.Here, {|0 , |1 } are the lowest energy eigenstates of ĤS and the quasiparticle spectral density S qp (ω) now depends on the distribution function which may, in general, includes both equilibrium and non-equilibrium contributions.Note that the decay rate in Eq. (3) stems from the sin φ/2 dependence of the quasiparticle-system coupling and is crucial to the memristive behavior detailed in the following section.By using the properties of displaced number states (see Supplementary Information), the squared inner products in Eq. ( 3) have a convenient cosine form valid for any pair of Fock states {|n , |m }, with Here, g 0 = [E C /(32E L )] 1/4 and L y x denotes an associ-ated Laguerre polynomial.Notably, the sign of the cosine term in Eq. ( 4) depends on the parity difference between the states involved.While this potentially provides insight into interesting phenomena when multiple decay channels are involved [21,22], we concentrate on the two-level process and leave such considerations for future studies.To understand how memristive behavior emerges from quasiparticle tunneling, we study the charge flow in the device.Let â be the annihilation operator for a harmonic excitation in the system.This allows us to write φ = 2g 0 (â + â † ) + ϕ d (t) and n = i(â † − â)/(4g 0 ), and denote by φind = φ − ϕ d the operator for the phase over the rf SQUID loop inductance.The directional convention for the superconducting phase differences and the different currents are presented in Fig. 1

(b). The average charging current Îch and the inductive current
Îind can be rigorously derived (see Supplementary Information) to obtain, by current conservation, the average quasiparticle current through the effective junction.The result is Îqp = 2eTr{ D{ρ}n} = Γ 1→0 (−e) n , which corresponds to the dissipative current induced by the interaction with the quasiparticle bath represented by the dissipator D{ρ}.Using V = −2e n /C d , the average quasiparticle current is determined by where we have preemptively written the effective conductance as a function of the selected memory variable φ , input V , and time t.Solving for the dynamics, we obtain where the average inductive phase difference only requires knowledge of the input via and we denoted the initial system coherence in the energy eigenbasis by ρ 01 (0) = 0|ρ|1 | t=0 .The memory variable update function in ∂ t φ = f [ V , t] only depends on the input and time, and has the explicit form Equations ( 5)- (7) indicate that a simple superconducting device operates as a voltage-controlled quantum memristor when the average voltage over a tunneling element is interpreted as the system input, the average quasiparticle tunneling current as the output, and the average superconducting phase difference as the memory retention variable.The quasiparticle conductance acts as the memductance corresponding to the memory-dependent average current response.It should be noted that physically speaking our device is considered a flux-controlled memristive device as it includes non-zero capacitive and inductive elements [2] while having no external capacitive coupling.However, only considering the quasiparticle contribution to the current and studying the abovementioned equations allows us to define the device as a voltage-controlled memristor from an operational pointof-view.
The operation of the constructed memristor is of ensemble nature, that is, the system input and output are quantum averages obtained from the measurement record of the corresponding observables.Experimental input consists of initialization and a slowly oscillating flux bias applied to the rf SQUID loop.In this way, one obtains in-dependently generated records which, consequently, have a complex correlation exhibiting memory features via Eq.( 5).In fact, the selected system input is not independent of the decay, but experiences a memory-dependent damping which allows one to self-consistently solve the fundamental equations above.One such solution is identifiable as mimicking the operation of the classical superconducting memristor [14], in which the memory is fully stored in the phase bias.It is obtained in the weak-damping limit by initializing the system with V | t=0 = V 0 and φ | t=0 = ϕ d (0), and by assuming a resonant sinusoidal phase bias , where the update is given by the classical Josephson relation ∂ t ϕ d (t) = 2e/ V .The solution embodies the two implicit assumptions for the classical memristor: (1) the rf SQUID loop has a negligible inductance, and (2) the internal dynamics is negligibly affected by the same dissipation that produces the output.
As a first step, we need to verify whether the abovedescribed classical-limit solution is consistent with the semiclassical results of Ref. 14.In Fig 2, one clearly sees that we observe the hysteretic current-voltage characteristic curves as required for a memristive element.In other words, a proper choice of the sinusoidal drive allows for tunable finite-area pinched loops [3].Employing the system parameters from Fig. 2, the above weakdamping solution is accurately numerically retrieved with S qp (ω 10 ) = 10 −4 ω 10 over multiple oscillation periods.This corresponds to a minimum relaxation time during the driving period of min(T 1 ) ∝ 1 µs relevant to the current state-of-the-art experimental setups [15,16].While those setups consider a different type of system, the fluxonium, very little experimental work has been able to reach the regime in which quasiparticle-induced relaxation is observable and, consequently, we use these references for initial comparison.Even though Îqp ∝ S qp (ω 10 ), the magnitude-scaled hysteresis curve is robust against decreasing the minimum T 1 -time by 2 orders of magnitude (see Supplementary Information).Beyond this, the input and output values are subject to noticeable decay.We show the parametric dependence of the average voltage and quasiparticle current in Fig. 2 for S qp (ω 10 ) = ω 10 corresponding to min(T 1 ) ∝ 100 ps.The hysteresis curve starts from a point in the weak-damping trajectory due to the identical initialization, and it is followed by a reduction of the area with time.The time evolution in Fig. 2 shows a gradual decay in the voltage and current amplitudes.Note that the system is operated in the phase regime of almost negligible loop inductance.This allows for a feasible resonant phase biasing frequency ω 10 /2π ≈ 45 GHz, achieved while ensur- ing sufficient adiabaticity max(α rs ) ≈ 0.15 (see Supplementary Information), necessary for the master-equation treatment employed for the quasiparticle bath.The initialization of the system plays a crucial role in the operation and does not simply determine the initial position in the parametric curve.Figure 3 shows the hysteresis curves for three different initializations, assuming the weak-damping limit and a resonant sinusoidal drive protocol.These curves can be interpreted by studying the time symmetry of the quasiparticle current between two consecutive crossings of the zero-energy point and indicate a tunable landscape of hysteretic behavior (see Supplementary Information).
To quantify the non-Markovian [23] character of the device, we consider the area enclosed by a hysteresis loop in the current-voltage plane as a memory measurement.This interpretation is founded in the observation that the absence of area correlates with purely time-local current response.In other words, a nonlinear conductance cannot produce a non-zero area since it depends only on the instantaneous value of the voltage.The memory quantifier for the kth traversed loop FIG. 3. Hysteresis curves with resonant sinusoidal phase bias in the weak-damping limit.System initialized such that The system parameters are the same as in Fig. 2 with Sqp(ω10) = 10 −4 ω10.
where t k c fulfills V | t=t k c = 0 for each k.This quantifier stores the evolution information of corresponds to the response specific to the selected memory variable, and the second term to the explicit time dependence of the memductance not included in the internal memory variable.However, it is in principle always possible to redefine the memory variable to absorb the explicit time dependence in the memductance, so that F m (t) = F 0 (t).The two expressions in Eq. ( 9) imply that the quantifier corresponds to a time-dependent weighted record of the change in the memory variable ∂ t φ or the instantaneous distance from its initial value If the conductance is a non-linear function of only the instantaneous input, F m (t) vanishes in integration due to input periodicity.See Supplemental Information for the decay of the quantifier as well as its response to different initializations.
Finally, our quantum memristor is formulated in the ideal case of zero leakage supercurrent.Adding a nonzero pair-tunnelling term, not only modifies the energy and state structure, but inflicts a Josephson tunneling current which may disrupt the operation.While there can be multiple factors contributing to the leakage supercurrent, such as magnetic flux noise, the primary experimental factor to tackle is possibly the critical current imbalance in the SQUID.The state-of-the-art critical current suppression factor based purely on fabrication techniques is ∼10 −2 while the balanced SQUID [24] promises a factor of ∼10 −3 -10 −4 , for a maximum critical current of 30 nA.In terms of the Hamiltonian, this implies that the imbalance term is 10 −1 -10 −3 times the charging energy scale used here.In addition, our formulation assumes only the quasiparticle decay channel and omits other natural loss channels (dielectric, inductive, radiative).Recent experimental work has studied quasiparticle-limited relaxation and shown significant progress in suppressing the additional decay channels, modifying the quasiparticle population through different means, and discerning between the different decay mechanisms [15,16].
Let us finish with a brief discussion about the quantumness of the system, as well as the role of superposition and entanglement.The dynamics of the quantum memristor described above is purely quantum, in the strict sense that the evolution cannot be emulated by a classical channel [25,26].This is not surprising, since the quantumness of our design refers to the full dissipative treatment (as an open quantum system) of the quasiparticle bath leading to memristive features in the expectation values of quantum observables.Therefore, superposition plays the same role as in any other quantum system.With respect to the entanglement, coupling two of these quantum memristors is a natural and relevant question after showing the dynamics of a single device, but beyond the scope of this manuscript.
In conclusion, we have demonstrated a prototype design for a quantum memristor in a superconducting circuit relying on quasiparticle tunneling.The pinched hysteretic behavior of the average quasiparticle current is a clear signature of conductance beyond typical nonlinearity, and modified by both the characteristics of the circuit and the quasiparticle bath.The measurement resolution can potentially be varied by tuning the nonequilibrium quasiparticle population, by just using the state-of-the-art injection and trapping methods [16] during the lifetime of the quasiparticles.Our work paves the way for the engineering of on-demand quantum non-Markovianity using the superconducting quantum memristor as a building block.Furthermore, we may consider possible applications such as the codification of Quantum Machine Learning protocols [27,28] and neuromorphic quantum computing [13].
While it is well-known that a SQUID behaves as if it was a single effective junction with a tunable Josephson energy when pair tunneling is discussed, additional conductance requirements are set by the inclusion of quasiparticle excitations.Classically speaking, a nonvanishing phase-dependent conductance can be achieved in a conductance-asymmetric SQUID, while assuming that the physical junctions have equivalent critical currents [14].Summation of dissipative currents yields a representation as an effective junction with an effective leakage conductance G eff = G 1 + G 2 , where G i is the leakage conductance of the ith physical junction, and an asymmetry term for the phase-dependent current of To use the effective junction picture in the main text in association with quasiparticle tunneling, we assume strong asymmetry G 1 G 2 such that S asym ≈ 1.Hence, the dissipative flow is effectively only through a single junction, while the total pair current can still be cancelled.This allows for the effective junction picture to be used for describing quasiparticle tunneling.Since the Ambegaokar-Baratoff relation implies for equivalent critical currents that G 1 /G 2 = ∆ 2 /∆ 1 , with ∆ i is the superconducting gap of the electrodes of the ith physical junction of the SQUID, this assumption can be enforced by demanding ∆ 2 ∆ 1 .Weakening the assumption would require the degrees of freedom of the individual physical qubits to be included in the Hamiltonian in Eq. ( 2) of the manuscript.This means that the quasiparticle tunneling Hamiltonian would have to be given for each physical junction separately, resulting in two separate decay channels, which allows for the weaker conductance discrepancy to manifest.We leave such considerations for future work and assume the limit of strong conductance asymmetry throughout the main text.Adiabaticity.-Inorder to evaluate the adiabaticity of our phase-driven system described by the Hamiltonian in Eq. ( 1) of the main text, we calculate the instantaneous adiabatic parameter for the dynamics confined to the two lowest energy levels.The instantaneous eigenstates of the abovementioned Hamiltonian are the well-known instantaneous number states where 4 and H n is the nth Hermite polynomial.By using known properties for the Hermite polynomials, we obtain two useful identities By taking a time-derivative of the instantaneous eigenstate in Eq. (B1) and applying the identities given by Eqs.(B2) and (B3), we obtain (B4) where the remaining time-derivative depends on the details of the external drive protocol.

(B6)
In order for the master equation approach in the main text to be valid, Landau-Zener transitions must be suppressed at all times, such that any population transfer is dissipation-induced and occurs between instantaneous eigenstates.This is guaranteed by operating in the regime max(α rs ) 1.
Quasiparticle-induced average shift of the system transition frequency.-Theexistence of quasiparticles induces an average frequency shift for the system that can be attributed to two different mechanisms [19,29]: the quasiparticle renormalization of the Josephson energy and the quasiparticle-mediated virtual transitions between different energy levels.In general, each physical junction generates different shift terms, since nonvanishing total quasiparticle current requires conductance asymmetry (see previous section) and, hence, different superconducting gaps for the electrodes of the junctions (∆ 1 = ∆ 2 ).However, our assumption of strong conductance asymmetry (∆ 1 ∆ 2 ) implies that the frequency shift is dominated by the junction with large quasiparticle flow, so we denote the effective gap by ∆ ≈ ∆ 1 .The renormalization term comprises of contributions from the pair tunneling and Josephson counterterms ( ω 10 , δE 2∆), as well as the terms in the quasiparticle tunneling Hamiltonian that do not contribute to pure dephasing and relaxation [19].Since these contributions are renormalized Josephson energy terms for the individual junctions, they vanish for our effective junction.For a nonvanishing Josephson term, the quasiparticle contribution can be approximated knowing the CP-density-normalized quasiparticle density x qp and the energy mode occupation of the quasiparticles at the gap x A qp = f E,qp (∆), where f E,qp is the energy mode of the lead quasiparticle distribution function [29].
By using the assumptions detailed in the main manuscript, the principal contribution from the virtual transitions to the ith energy level of the system is [29] where |i is the ith eigenstate with energy E i , ω ki = E k − E i , and F qp (ω) is an expression involving complex nested integration of the quasiparticle distribution in energy space [see Appendix A of Ref. [29]].As it is apparent from Eq. (B7), each energy shift generally accounts for virtual transitions to and from each other state in the Hilbert space.Our system lacks anharmonicity and, hence, we cannot restrict virtual occupation to the two lowest levels.However, we operate in the phase regime where 1 and, hence, Eqs. ( 4) and ( 5) in the main text imply that the largest terms in δE i,qp correspond to virtual transitions between i and its energetically nearest levels.Thus, the frequency shift becomes , where the non-nearest-neighbour terms are of the order g 4 0 by the constuction of the inner product in Eq. (B7), whilst the nearest-neighbour terms scale as e −g 2 0 g 2 0 , which was expanded to obtain the principal term in Eq. (B8).Notably, the lowest-order contributions stem from the virtual transitions 1 ↔ 0, 0 ↔ 1, and 2 ↔ 1. Making use of the definition of quasiparticle impedance and assuming the high-frequency limit, we finally obtain where g T ∝ G eff is the effective junction conductance and e is the elementary charge.This term is generally time-dependent, due to the phase steering.We assume that the physical parameters are selected such that δω 10,qp ω 10 , so that the frequency shift can be omitted.As evident from Eq. (B9), the validity of this assumption is determined by the details of the execution of the low-energy limit given in the main text hindering it and the execution of the phase-regime as well as typical low quasiparticle densities supporting it.Note that in thermal equilibrium and at low temperatures T ∆, it is straightforward to support our statement that the frequency shift is dominated by the junction with large dissipative flow (∆ 1 ∆ 2 ) using Eq.(B9), since then δω 10,qp ∝ e −∆/k B T .This implies that δω 1  10,qp /δω 2 10,qp ∝ e (∆2−∆1)/k B T 1, where the superscript indicates the physical junction.
Pure dephasing.-Forour system, whose memristive operation relies solely on a weak but nonvanishing dissipative current generated by the quasiparticle-induced decay, pure dephasing adds another decoherence channel and, might potentially destroy the memristive behavior.Its effect would be to both distort the hysteresis loops by adding another time-dependent contribution to the quasiparticle current and to increase the memorydependent damping of the average voltage by contributing to the total dephasing.Even though our SQUID is necessarily asymmetric as explained earlier, the individual junctions are assumed symmetric throughout this work.Thus, one would estimate pure dephasing using the self-consistent rate Γ φ in Eq. ( 28) of Ref. [19] to avoid the issue of logarithmic divergence in the lowest-order perturbative tunneling treatment.However, the self-consistent rate scales with the system energy as |A d s | 2 where A d s = ( 1| sin φ/2|1 − 0| sin φ/2|0 )/2 and the inner products can be calculated using the identities in Eqs.(C1) and (C2).This yields |A d s | 2 = g 4 0 e −g 2 0 sin 2 (ϕ d /2)/4, which is of the order g 4 0 , while the relaxation rate scales as g 2 0 .A full comparison of the relaxation and pure dephasing rates would require knowledge of the quasiparticle distribution to calculate the quasiparticle spectral density in Eq. ( 3) of the main text, as well as to iteratively solve the self-consistent expression for the pure dephasing rate.As an example, assuming Γ φ δE, the scaling becomes , which decreases faster than Γ 1→0 in small g 0 .In this work, we assume that the quasiparticle distribution can be established in a manner that the beneficial difference in scaling in the phase regime allows for the pure dephasing process to be neglected.
Finally, we remind the reader that our construction of the relaxation rate, as well as the considerations on the pure dephasing rate above, do not only exploit the adiabatic assumption but also that of low charateristic quasiparticle energies ( ω 10 , δE 2∆).This implies that the secondary terms in the rate equations in Ref. [19] can typically be omitted, since they are negligible in comparison to the primary terms we use throughout this work.However, our system crosses points during the phase-driving process in which the primary terms vanish and the secondary terms become the dominating contributions, ensuring that the total rates are nonvanishing.For example, the secondary term in the relaxation rate is proportional to e −g 2 0 g 2 0 (1 − cos ϕ d )/2, which reaches its maximum value at the point of vanishing primary term proportional to e −g 2 0 g 2 0 (1 + cos ϕ d )/2, that is, at ϕ = π + 2πn, n ∈ Z.Even though the secondary terms dominate near these points, we assume that their contribution to the total dissipative flow during the full driv-ing cycle is negligible due to the aforementioned assumptions.If low characteristic quasiparticle energies cannot be guaranteed, the full rate equations should be applied to study the potentially modified memristive function.
Appendix C: Matrix element of sin φ/2 in the quasiparticle decay rate By using φ = 2g 0 (â+â † )+ϕ d (t), where â is the bosonic annihilation operator for the harmonic system, the sinusoidal phase term can be written as where D(ig 0 ) = exp[ig 0 (â + â † )] is the displacement operator with ig 0 = i[E C /(32E L )] 1/4 the phase space displacement.General properties for displaced number states assert that, for any pair of eigenstates of the harmonic oscillator |n and |m , we have the identity [30] m|  3) and (4) in the main text.It should be noted that a similar calculation has been performed in Ref. [29] for the pair (|n , |0 ) and the results presented in the main text are a generalization to any pair of Fock states.The squared inner products corresponding to the transitions between the three lowest energy eigenstates are (C4) Appendix D: Derivation of the charging, inductive, and quasiparticle currents The average charging current for the junction is given by where we assumed a Lindblad-form evolution ∂ t ρ = −i/ [ ĤS , ρ] + D{ρ} for the system density ρ, and used the properties of the bosonic operators to rewrite the commutator [n, ĤS ] = −i ω 10 (â † + â)/(4g 0 ), after applying Ehrenfest theorem.By employing the operator for the current through the inductive element Îind = −2e/ E L ( φ − ϕ d ), the average inductive current is The same decay that enables hysteresis is responsible for the gradual decrease of the quasiparticle current due to dephasing.Observation of hysteresis requires a sufficently large T 1 -time to be achieved.To estimate the bath-induced decay of the input and output values with respect to the scaling of the T 1 -times, we present the change in the average quasiparticle current after 10 consecutive driving cycles in Fig. 4. Since the system is initialized to the classical trajectory, Îqp | t=0 corresponds to the value that the classical memristor returns after each driving cycle.Note that the memory-dependent damping in the input V determines the decay and, hence, has the same scaling as the output.As briefly mentioned in the main text, the magnitude-scaled parametric hysteresis curves do not noticeably decay after 10 cycles for the orders of magnitude of the T 1 -times between 1 ms and 10 µs.This is noticeable in Fig. 4, where the quasiparticle current after 10 cycles does not return to its original value when log[S qp (ω 10 )/ω 10 ] > −2.
Appendix F: Effect of initialization on the hysteretic behavior and memory quantifier Here, we provide a qualitative interpretation for the varied hysteretic curves displayed in Fig. 3 of the main text.As in Fig. 2 of the main text, initially assuming a nonzero voltage over the SQUID and a negligi-ble inductive phase difference yields hysteresis with πrotational symmetry in the current-voltage plane (blue curve).On the other hand, assuming vanishing initial voltage and nonzero inductive phase difference destroys the hysteretic behavior (red curve).Finally, assuming that both are nonvanishing allows for voltage-asymmetric hysteresis (black curve).Each case can be seen as an example of the time-symmetry of the quasiparticle current between two consecutive crossings of the zero-energy point.With sinusoidal driving, both the voltage and the memristance are time symmetric, whereas the symmetry of their product, the quasiparticle current, depends on the specific initialization.Hence, tuning the initial conditions yields a variety of different memristic switchings which can, in general, be made voltage asymmetric.The parametric current-voltage characteristics can also be tuned by resorting to different bias protocols which may help to achieve desired behavior.
The memory quantifier for Fig. 2 in the main text decays pairwise linearly from N 1 m /N cl m = 0.9960 to N 19  m /N cl m = 0.8511, where N cl m is the value for the classically initialized weak-damping solution, during the simulation.In the weak-damping regime, the decay is negligible.For the different initializations in Fig. 3 of the manuscript, the weak-damping solutions yield N m /N cl m = 1 (blue), N m /N cl m ∝ 10 −7 (red), and N m /N cl m = 1.3409 (black).Curiously the last case yields the same value for both asymmetric loops.
Appendix G: Area of a hysteresis loop as a memory measurement For a classical voltage-controlled memristive system, the current response is generally given by [2] I(t) = G(γ(t), V (t), t)V (t), where G(γ, V (t), t) is the memductance, which depends on the instantaneous value of the input voltage V (t), the accumulated value of the memory variable γ(t), and the time t in a parametric manner.The memory variable update is defined by where the update function f depends on the past evolution, the input and any external parametric driving.Notice that the general definitions used above allow for tracking a specific state variable γ rather than absorbing the total response into a generalized memory variable γ x such that G(γ x ) = G(γ, V, t) where γ x = f x (γ x , V, t) now updates the new variable.This absorption can always be done, since memductance is a system property depending on the input and external driving via the state variable.Define C as any closed curve in the (I, V )-space, hence the area enclosed by the curve is determined by Green's theorem as (G3) where the explicit dependences given by the response in Eq. (G1) are omitted for clarity.Transforming to temporal space, the integration takes the form where we have fixed the loop to begin at time t = 0 and end at t = T when the initial point in the (I, V )space is reached again.Assuming pinched hysteresis, the periodicity condition is conveniently written as V (nT ) = 0, n ∈ Z. Due to this condition, the second term in Eq. (G4) yields implying that any instantaneous voltage dependence in the current response does not contribute to the area.By making use of Eq. (G2), the area is (G6) where the first term in the integrand corresponds to the response related to the selected memory variable and the second term accounts for any remaining explicit time dependence of the conductance.Using the notation above, any non-linear conductor G = G(V ) can be defined via γ = V such that f (γ) = ∂ t V .Equation (G6) then yields which is zero, A = 0, due to the periodicity.In other words, purely non-linear conductance does not produce hysteresis, and cannot generate a non-zero loop area.The memory-variable related term in Eq. (G6) can be rewritten as where ∆γ(t) = γ(t) − γ(0) is the change in the memory variable from its initial value, we executed integration by parts after the first equality, and the first term after the first equality vanishes due to the periodicity condition.
Applying the results in Eqs.(G6) and (G7) yields the definition of the loop area as a memory quantifier given in the main text for our superconducting circuit.

FIG. 1 .
FIG. 1. Superconducting quantum memristor.(a) Schematic representation of the superconducting quantum memristor.The green and red strips represent junctions with different normal conductances.(b) Current diagram using effective circuit elements corresponding to the total loop inductance, charge retention of the SQUID, and the quasiparticle tunneling through it.Let us remark that, as the capacitative part has been explicitly separated in (b), the cross-notation does not refer here to the entire (effective) Josephson junction, but to the quasiparticle and phase-dependent dissipative current contributions.
Appendix B: Estimation of adiabaticity, quasiparticle-induced average frequency shift, and pure dephasing

FIG. 4 .
FIG.4.Expectation value of the quasiparticle current after 10 resonant sinusoidal driving cycles with respect to the spectral density of the quasiparticle bath.System parameters and initialization as in Fig.2of the main text.