Parity-time symmetric systems with memory

Classical open systems with balanced gain and loss, i.e. parity-time ($\mathcal{PT}$) symmetric systems, have attracted tremendous attention over the past decade. Their exotic properties arise from exceptional point (EP) degeneracies of non-Hermitian Hamiltonians that govern their dynamics. In recent years, increasingly sophisticated models of $\mathcal{PT}$-symmetric systems with time-periodic (Floquet) driving, time-periodic gain and loss, and time-delayed coupling have been investigated, and such systems have been realized across numerous platforms comprising optics, acoustics, mechanical oscillators, optomechanics, and electrical circuits. Here, we introduce a $\mathcal{PT}$-symmetric (balanced gain and loss) system with memory, and investigate its dynamics analytically and numerically. Our model consists of two coupled $LC$ oscillators with positive and negative resistance, respectively. We introduce memory by replacing either the resistor with a memristor, or the coupling inductor with a meminductor, and investigate the circuit energy dynamics as characterized by $\mathcal{PT}$-symmetric or $\mathcal{PT}$-symmetry broken phases. Due to the resulting nonlinearity, we find that energy dynamics depend on the sign and strength of initial voltages and currents, as well as the distribution of initial circuit energy across its different components. Surprisingly, at strong inputs, the system exhibits self-organized Floquet dynamics, including $\mathcal{PT}$-symmetry broken phase at vanishingly small dissipation strength. Our results indicate that $\mathcal{PT}$-symmetric systems with memory show a rich landscape.


I. INTRODUCTION
Over the past decade, open systems with balanced, spatially separated gain and loss have become a rich area of research. They are described by a special class of non-Hermitian Hamiltonians that are invariant under spaceand time-reflections, i.e. parity-time (PT ) symmetric Hamiltonians [1]. More than two decades ago, Bender and coworkers first introduced a broad class of such continuum Hamiltonians on an infinite line [2], and showed that, in spite of their non-Hermitian nature, they have purely real spectra when the non-Hermiticity is small. The initial, theoretical studies of PT -symmetric Hamiltonians were focused on developing a complex extension of quantum theory [3][4][5]. Over the past decade, however, it has become clear that such Hamiltonians describe the dynamics of classical energy density within different parts of a system, in the presence of localized sources or sinks [6][7][8][9]. A non-Hermitian Hamiltonian H(γ) is called PT symmetric if it commutes with an antilinear operator PT , where P is the parity operator satisfying P 2 = 1, and T = U K where U is a unitary operator and K denotes the (antilinear) complex conjugation operation. At γ = 0, H is Hermitian, has real eigenvalues and a complete set of orthonormal eigenvectors that lead to a unitary time evolution. At small non-Hermiticity γ < γ PT , the spectrum of H(γ) is purely real, but the non-orthogonal eigenvectors generate a nonunitary, bounded, oscillatory-in-time dynamics. At a critical gain-loss strength γ = γ PT , eigenvalues of H become degenerate, as do the corresponding eigenvectors. Such degeneracies, where the non-orthogonal eigenvectors of H(γ PT ) do not span the space, are called exceptional points (EPs) [10][11][12]. Beyond the EP, eigenvalues of H(γ) occur in complex conjugate pairs. Due to the an-tilinear nature of PT operator that commutes with H, an eigenstate |ψ α of H with eigenvalue α is a simultaneous eigenstate of the PT operator with eigenvalue +1 if and only if α is real; when α is complex, PT |ψ α gives rise to the eigenvector corresponding to the complex conjugate eigenvalue K α = * α . The transition from a PT -symmetric region (real spectrum) to the PT -symmetry broken region (complex conjugate spectrum) across the EP has been extensively studied in classical wave systems where both gain and loss are readily implemented. Realizations include coupled optical waveguides [13], fiber loops [14], microring resonators [15], acoustic setups [16], coupled mechanical oscillators [17], and coupled electrical circuits [18,19]. Due to quantum fluctuations associated with a linear gain [20], balanced gain and loss configurations are not possible at a quantum level [21]. However, EP degeneracies are also present in dissipative systems with mode-selective losses, thus extending the ideas of PTsymmetry into the quantum domain, where the passive, PT -symmetric systems [22,23] have been realized with tabletop [24] and integrated quantum photonics, ultracold atoms [25], a single NV center in diamond [26], and a single superconducting qubit [27]. Most of these systems are modeled with a static, PT -symmetric Hamiltonian whose eigenvalues and eigenvectors determine the PT -phase diagram of the system. This landscape is dramatically transformed when one considers PTsymmetric Hamiltonians that are periodic in time with period T [28,29]. In this case, the eigenvalues α (t) of the instantaneous Hamiltonian do not govern the system dynamics; instead, the PT -phase diagram is determined by an equivalent static Hamiltonian H F called the Floquet Hamiltonian [25,30,31]. Another level of complexity is added when we consider systems described by a nonlin-ear, PT -symmetric Schrodinger equation [32]. However, in almost all cases [33] the dynamics are Markovian, i.e. the state of the system at the next instance depends only on its state at present, but not on its history.
Here, we introduce PT -symmetric systems with memory. Formally, this is to be achieved by making either the gain-loss or the Hermitian part of the Hamiltonian dependent on the history of the system. We use coupled, active and lossy LC circuits as the model [18,19,30], because resistive and inductive elements with memory, i.e. memristors and meminductors, are well understood. The plan of the paper is as follows. In Sec. II we review properties of a lossy RLC circuit inductively coupled to an LC circuit with negative resistance −R. Modeling and numerical results for a system where the resistance is replaced by a memristor (and the negative resistance is matched) are presented in Sec. III and Sec. IV. Section V has corresponding results for a system where the coupling between the two LC circuits is mediated by a meminductor instead of a regular inductor. We conclude with discussion in Sec. VI.

II. COUPLED LC CIRCUITS WITH PT SYMMETRY
We start with the review of an electrical PT -symmetric dimer [18,19,30,31]. Let us consider two identical LC circuits, with effective, parallel resistors ±R respectively, that are connected with a coupling inductor L c as shown in Fig. 1a. When the two circuits are uncoupled, the energy in the standard RLC circuit undergoes overdamped or underdamped decay, while the energy in the −RLC circuit, with negative effective resistance, grows with time reflecting a time-reversed dynamics. The state of the coupled system is characterized by a real vector is the voltage across the first (second) capacitor, I 1(2) (t) is the current through the first (second) inductor, and I c (t) is the current flowing through the coupling inductor L c . Its equation of motion, determined by the Kirchhoff laws, can be written as i∂ t |φ(t) = M |φ(t) where the purely imaginary, 5 × 5 matrix M of rank 4 is given by To map the Kirchhoff-law equations into the Schrodinger-like equation, we note that the circuit energy is given by E(t) = φ(t)|A|φ(t) where A = diag(C, C, L, L, L c )/2 is a diagonal matrix. Therefore, we consider the dynamics of energy density by defining |ψ = A 1/2 |φ such that the norm of the state |ψ is the total energy in the circuit. This mapping transforms the Kirchoff-law equations into where the purely imaginary, rank-4 Hamiltonian H eff = A 1/2 M A −1/2 is given by Here, ω 0 = 1/ √ LC is the fundamental frequency of a single LC circuit, µ = L/L c is the dimensionless inductive coupling between the lossy and active oscillators, and γ = 1/RC = Γω 0 is the gain-loss rate for the circuit. When γ = 0, the Hamiltonian H eff is Hermitian, and the norm of the state, i.e. the total energy in the circuit, is conserved. We also note that H eff (γ) is PT -symmetric with respect to where σ x is the standard Pauli matrix and 1 k is a k × k identity matrix. H eff also satisfies ΠH eff = −H eff Π where Π = PU . This chiral symmetry is responsible for its spectrum consisting of a trivial zero and two pairs of particle-hole symmetric eigenvalues [34] given by (5) It is easy to check that α are purely real for γ ≤ γ PT = ω 0 ( 1 + 2µ 2 − 1), and become purely imaginary when γ > γ c = ω 0 ( 1 + 2µ 2 + 1). Figure 1b shows the flow of the four eigenvalues in units of ω 0 as a function of the gain-loss strength γ/γ PT for a strongly coupled LC circuit with µ = 1.
When the dissipation γ is time-dependent, the electrical-circuit dynamics still maps onto a Hamiltonian given by H eff (t) = A 1/2 M (t)A −1/2 . On the other hand, if a conservative circuit element varies with time, the change-of-basis matrix A 1/2 (t) leads to a new Hamiltonian [31] Whether the additional gauge term in Eq.(6) commutes with the PT operator depends on the functional timedependence of matrix A(t). However, (with a little license with notation) we will continue to call such system a PT -symmetric dimer. In the following sections, we will consider models where the time-dependence arises through memory in either the gain-loss strength γ or the inductive coupling µ.

III. MEMRISTIVE PT -SYMMETRIC MODEL
A memristor is a resistor whose resistance R(x) depends on an internal, dimensionless state variable x [35].  The equation of motion for the state variable x(t), in turn, is determined by the underlying microscopic model for the memristor. A memristor (memory resistor) was postulated by L. Chua more than half a century ago based on symmetry arguments [36]. It was realized just over a decade ago in a thin-film device with one monolayer of TiO 2 and its oxygen-vacancy doped counterpart, TiO 2−δ [37,38]. Since then it has become clear that memristive systems, with their pinched-loop hysteresis signature [39] in the current-voltage characteristics, manifest themselves in semiconductor thin films, thermistors [40], as well as ion channels [41] in biological membranes [42]. Here, we will consider the simplest model for its internal state variable [37,43]. The TiO 2 /TiO 2−δ thin film with size D (D ∼ 5 nm) can be modeled as two resistors in series, where the size of the doped part is given by w = xD, and the size of the undoped region is D − w = (1 − x)D (Fig. 1c). The resistance of this two-terminal passive device is given by where R on (∼1 kΩ) is the resistance of the device if it is entirely doped and R off ∼ 10 2 R on R on is the resistance of the insulating TiO 2 device. When voltage is applied to such a device, in addition to conduction electrons, the charge +2 oxygen vacancies also move. The effect of their motion is amplified due to the two-monolayer thickness of TiO 2 /TiO 2−δ , and it determines the fractional width x(t) of the doped region. By equating the rate of change of x(t) with the drift velocity of the oxygen vacancy dopants, we get Here I is the electronic current through the memristor, Q 0 = D 2 /µ D R on is the characteristic charge-scale for the memristor, the dopant mobility µ D (∼ 10 −10 cm 2 V −1 s −1 ) is 10-12 orders of magnitude smaller than the corresponding electron mobilities, and η = ±1 signifies whether the doped region is shrinking or growing, i.e. polarity of the memristor. The window function F (x) in Eq.(8) suppresses dopant mobility when the interface between undoped and doped regions approaches the device boundaries, i.e. x → 0 or x → 1. We use a family of window functions [43]. Since x = {0, 1} are fixed points of Eq.(8), if the time-evolved state variable x(t) reaches either fixed point, it has no further dynamics. However, the amount of charge required to change the state-variable value from x = δ 1 1 to . Thus, starting from an x ∈ (0, 1) it is impossible to reach the boundaries in finite amount of time [44]. In calculations, if x(t 0 +δt) reaches or exceeds the fixed-point values while x(t 0 ) ∈ (0, 1), either a smaller time-step ∆t is chosen or the updated x-value is shifted to just inside the boundaries to circumvent this numerical artifact.
We now consider a PT -symmetric dimer with a memristor R(x) and a balanced gain-resistor that matches the instantaneous dissipation in the memristor. Such a setup is experimentally feasible with lock-in amplifiers and synthetic elements. The circuit dynamics now are described by two coupled, nonlinear, first-order equations, where γ(x) = 1/R(x)C denotes the variable gain-loss strength, and 1|ψ(t) is the first element of the energydensity state vector |ψ(t) . Note that due to nonlinearities, the system dynamics depend on the initial statenorm or equivalently, the initial energy in the system. Therefore, the fate of the PT -symmetric system with memory cannot be analyzed in terms of its state-variable dependent Hamiltonian H eff (γ(x)). Instead, we track the time-dependent energy E(t) = ψ(t)|ψ(t) in the circuit. As in the standard PT -symmetric case, this circuit energy E(t) shows bounded, periodic, oscillatory behavior for some parameter regime, while for others, it shows non-periodic, divergent behavior. As we will demonstrate below, this behavior can sensitivity depend on the initial state |ψ(0) as well. To characterize these two trends, we calculate the long-time amplification rate, When the system is in the PT -symmetric phase, the dynamics are oscillatory and therefore Λ amp = 0. In the PT -broken phase, due to the presence of amplifying modes, the energy grows exponentially with time and Eq.(11) characterizes its growth rate. Since the energy E(t) is obtained numerically, in practice, we must choose τ to be larger than any other relevant time-scale in the problem. For a PT symmetric circuit without memory, the fast timescales are given by T 0 = 2π/ω 0 and T 0 /µ 2 T 0 . The longest timescale is given by the inverse of the smallest eigenvalue difference, Eq.(5), and it diverges as one approaches the static PT -threshold γ PT . Thus, one needs data at arbitrarily long time-scales to distinguish a system in the PT -symmetric phase from one in the PT -symmetry broken phase. For a memristor, the statevariable dynamics timescale depends on the dopant drift velocity or, equivalently, the applied voltage strength v 0 , and is given by [43]. Equating these two time scales, we obtain the characteristic voltage scale v 0 = D 2 ω 0 /2πµ D ; this voltage gives rise to a linear drift of size D in one oscillation of the LC circuit. Note that since the memristor value is confined between R on ≤ R(x) ≤ R off , the dissipation strength γ(x) is also bounded between γ off and γ on . If both strengths are below γ PT , the time-averaged dissipation, defined as will be smaller than the threshold for any t. Then, in the absence of any Floquet resonances [28,29], the system dynamics will be bounded and oscillatory. Similarly, if both are above γ PT , so is the time-averaged dissipation, and with the same caveat, the system will be in the PT -symmetry broken region, as indicated by a divergent circuit energy growth. Therefore, for simulations, we choose 0 < γ off ≤ γ PT and γ PT ≤ γ on ≤ γ c . Figure 2 shows the typical energy dynamics E(t) for a memristive dimer as a function of γ off . The circuit parameters are µ = 0.3, γ on = 2γ PT , and x(0) = 0.5. The initial state vector |ψ(0) of the system is given by At low dissipation strength, γ off = {0.25, 0.5}γ PT , the scaled circuit energy shows bounded oscillations whose period increases as the dissipation is increased. At large dissipations, γ off = {0.75, 1}γ PT , the scaled circuit energy shows exponential growth that is characteristic of a PT -symmetry broken phase. The inset shows the corresponding dynamics of the internal state variable x(t).
In the PT -symmetric phase, the doped fraction x(t) undergoes small oscillations around its mean value, whereas in the PT broken phase, it reaches its extremal values, thus driving the dissipation between γ on and γ off over the timescale T 0 . The results in Fig. 2 are, for the most part, expected even for a PT -symmetric dimer without memory. In the following paragraphs, we will show the unique features that arise from the memristive nature of dissipation and the resultant non-linearity in this system. Figure 3a shows the temporal evolution of circuit energy for four initial states |ψ(0) = {20, 40, 60, 80}|ψ 1 , with circuit parameters µ = 0.3, γ on = 2γ PT , γ off = 0.3γ PT , and x(0) = 0.5. At small initial energies E(0) = ψ(0)|ψ(0) , the dynamics are oscillatory with a period that increases with increasing energy; this behavior changes over to exponential growth when |ψ(0) = 80|ψ 1 . In a linear PT -symmetric dimer -static or  Floquet -all of these states are equivalent to |ψ 1 and will exhibit identical results for the scaled circuit energy E(t)/E(0). Instead, due to the nonlinearity introduced by dynamics of the internal state variable, now the fate of the system depends on the initial circuit energy. Figure 3b shows the diverse energy dynamics that occur for four initial states given by |ψ(0) = ±|ψ 1 , ±80|ψ 1 ; the initial doped fraction is x(0) = 0.9 and rest of the circuit parameters are the same as in Fig. 3a. We see that the fate of the scaled circuit energy depends not only on the initial energy, but also on the sign of the initial voltage V 1 (0), or equivalently, the polarity η = ±1 of the memristor. This can also be interpreted as the phase of the initial state |ψ(0) . (Recall that the phase of the purely real state-vector |ψ(t) is restricted to 0 or π). Specifically, at small initial energies, i.e. V 1 (0) = ±0.5v 0 , the scaled energy shows an exponential-growth transient followed by an oscillatory behavior that persists at long times t/T 0 ∼ 10 3 (not shown). At higher initial energy, the system starting in the initial state |ψ(0) = −80|ψ 1 shows a clear oscillatory energy dynamics with minimal net amplification, whereas with initial state +80|ψ 1 , the scaled circuit energy undergoes exponential amplification. Figure 3c shows results for a moderately coupled dimer with µ = 0.5, γ on = 2γ PT , γ off = 0.3γ PT , x(0) = 0.9, and four initial states given by |ψ(0) = {−0.2, −26.8, −53.4, −80}|ψ 1 . As the negative prefactor of |ψ 1 is increased in magnitude, we see that the scaled circuit energy E(t)/E(0) settles into oscillations about a mean that is monotonically suppressed. When the initial circuit energy E(0) is small, V 1 (0)/v 0 = −0.1, the scaled energy starts oscillating after an exponential growth transient similar to that seen in panel b for small initial voltages. It is important to note that this behavior also contrasts with Fig. 3a, where the largest V 1 (0) value resulted in a PT symmetry broken phase.
Lastly, in Fig. 3d, we explore the behavior of the circuit energy for four initial states given by 75|ψ 1 , |ψ 2 = 75P|ψ 1 , |ψ 3 = L/2[0, 0, I 1 (0), 0, 0] T , and |ψ 4 = P|ψ 4 . The initial voltages V 1 (0) = V 2 (0) = 35v 0 or initial currents I 1 (0) = I 2 (0) are chosen such that all states have the same initial energy (or state norm); the rest of the circuit parameters are the same as those for Fig. 3a. We see that when the system starts with nonzero voltage on the dissipative LC circuit (75|ψ 1 ), after a growth transient, the scaled circuit energy oscillates with a mean that is below the initial circuit energy. In contrast, when the system starts with nonzero voltage on the amplifying LC circuit (|ψ 2 = 75P|ψ 1 ), the scaled energy diverges exponentially indicating a PT -symmetry broken state. In contrast, when we start with a nonzero inductor current in the dissipative LC circuit (ψ 3 ), the scaled energy shows exponential growth, while switching the nonzero initial current to the amplifying LC circuit (|ψ 4 = P|ψ 3 ) leads to stable oscillatory dynamics for the scaled circuit energy. These properties are dramatically different from those of a traditional, memory-less PT -symmetric system. In the latter, different distributions of the initial energy density only introduce temporal shift in the dynamics of scaled circuit energy.
The insets in Fig. 3 show the evolution of the average dissipation strength γ(t) relative to the PT -breaking threshold strength γ PT in a PT -symmetric dimer with no memory. A common pattern observed is that when their ratio is less than unity at long times, the scaled circuit energy shows oscillatory behavior, whereas if the ratio exceeds unity at long times, the scaled circuit energy grows exponentially. In the following section, we will show that this naive expectation is wrong.

IV. PT SYSTEM WITH SELF-ORGANIZED FLOQUET DYNAMICS
In realistic memristors, the resistance of the undoped region can be orders of magnitude higher than that of the doped regions. When the initial state voltage is sufficiently high, the internal state variable x(t) switches periodically between 0 and 1 at frequency ω 0 as shown in the inset in Fig. 2. It implies that the gain-loss strength in the PT -symmetric dimer effectively switches on and off, with γ off γ on , i.e. the system behaves as if it has a time-periodic gain-loss strength whose periodicity is given by T 0 = 2π/ω 0 .
It is known that in memory-less PT -symmetric systems with periodic non-Hermiticity, the landscape of exceptional points that separate the PT -symmetric phase from the PT -symmetry broken phase is dramatically altered relative to its static limit [28,29]. In particular, the PT -broken phase occurs at vanishingly small gain-loss strengths when the modulation frequency ω 0 is an odd sub-harmonic of the Hermitian energy gap ∆ ≡ ω 0 ( 1 + 2µ 2 − 1) of the system [25,30,31], where n ≥ 0. This analogy suggests that at µ n = (2n + 3)(2n + 1)/2, the memristive PT -symmetric system might show similar properties. The strongest PTsymmetry broken region should then occur at the first resonance µ 0 = 3/2 ≈ 1.225 for a vanishingly small gain-loss strength γ off γ on γ PT . We therefore carry out simulations with Eqs.(9)-(10) for circuit parameters γ off = 0.01γ PT , 0.1 ≤ γ on /γ PT ≤ 0.9; initial doped fraction 0.1 ≤ x 0 ≤ 0.9; and different initial states and memristor polarities. These initial conditions ensure that time-averaged dissipation remains below the static PT -symmetry breaking threshold at all times, i.e. γ(t) < γ PT . Figure 4a shows the typical plot for the amplification factor Λ amp , Eq.(11), in the µ − γ on plane. It is obtained with |ψ(0) = 70|ψ 1 , x 0 = 0.85, and τ = 100T 0 . We see that the system is in the PT -symmetric state (Λ amp = 0) for most of the region except in the vicinity of µ 0 = 3/2 = 1.225, where the amplification factor is positive and grows with increasing γ on . These qualitative features remain the same for different initial circuit energies, state variable values x 0 , and memristor polarities. Figure 4b shows the temporal evolution of the scaled circuit energy for |ψ(0) = 70|ψ 1 , γ on = 0.4γ PT , and three different couplings marked by squares in Fig. 4a. At µ = 1.1, the system is in the PT -symmetric phase with an approximately constant total energy. At slightly higher coupling µ ≈ µ 0 , the scaled energy shows a clear exponential growth indicating a PT -symmetric broken regime. As the coupling is increased further to µ = 1.3, the system again enters the PT -symmetric phase, albeit with an enhanced mean for the scaled circuit energy. The inset in Fig. 4b shows the dynamics of the internal state variable x(t) at the three coupling values. In the PT -symmetric region (µ = 1.1), the doped fraction oscillates without reaching both extrema; at µ = 1.225, this changes to square-wave oscillations between the two extema; and at µ = 1.3, the system is back in the PTsymmetric phase and the x(t) oscillation range is reduced again. Figure 4c shows corresponding time-averaged gain-loss strengths γ(t). We see that at µ = 1.1 and µ = 1.3, the  average γ saturates to ∼ 5% of the static threshold value, whereas for µ = 1.225, it saturates towards ∼ 20% of the static threshold. In all cases, however, it is significantly smaller than unity. This emergence of a positive amplification factor at specific couplings and very small gainloss strengths is a key hallmark of Floquet PT -symmetry breaking phenomenon [25,[28][29][30][31]. It can be qualitatively understood as follows: the energy oscillates between the lossy LC unit and the active LC circuit at frequency ∆(µ) that is governed by the dimensionless coupling. If it is on the gain unit when x = 1 and moves to the loss unit when x = 0, the circuit energy will continue to amplify even if the gain-loss strengths γ off and γ on are small relative to γ PT . However, this synchronization requires the memristor switching frequency to match the energy sloshing frequency, i.e. µ = 3/2. Our results for the PT -broken phase at small time-averaged γ in the vicinity of the primary resonance remain qualitatively unchanged irrespective of initial circuit parameters, the initial state energy, the memristor polarity, or a different distribution of energy density across different elements in the initial state. As is expected, the nonzero amplification factor increases with γ on , increases with E(0), and decreases with initial fraction of the doped region x(0). We also note that numerically, we do not find any evidence for PT -broken region at higher values of µ, suggesting the absence of such effect at odd sub-harmonics.
The wide range of dynamics displayed by the circuit energy in a memristive PT -symmetric dimer raises several considerations. First, the approximation of a constant negative resistor −R or memristor −R(x) breaks down at sufficiently large net amplification E(t)/E(0), just as the approximation of a constant gain coefficient breaks down in the optical domain. Therefore, in reality, the ex-ponential growth in the circuit energy will saturate and reach a steady-state value in the PT -symmetry broken region, just as it does in the optical domain. Second, due to the very large parameter space in Eqs.(9)-(10), gaining a global understanding of whether the trajectories of |ψ(t) in the four-dimensional space R 4 exhibit closed orbits, fixed points, open diverging orbits, or chaotic behavior is difficult. On the other hand, with increase in the studies of PT -symmetric electrical lattice models [45] and the easy availability of memristor emulators [46], experimental investigation of these systems seems highly feasible.

V. PT -SYMMETRIC DIMER WITH MEMINDUCTIVE COUPLING
A meminductor (memory inductor) is a two-terminal passive device whose inductance L c (y) depends on a dimensionless state variable y whose dynamics, in turn, are governed by the current I c flowing through the inductor [47]. Such a device shows a pinched hysteresis loop in the plane spanned by the current I c and the timeintegral of the voltage (called the flux) φ across the device [48]. The simplest, intuitive model of a meminductor is a solenoid with a ferromagnetic rod that can move in and out of its core [49]. However, the motion of the rod depends on the current I c through a second-order derivative, i.e. d 2 y/dt 2 ∝ I c [49], and sets it apart from the viscous, drift-velocity model for the internal state variable of a memristor, Eq. (8).
To circumvent this distinction, we only consider meminductors in which the internal state variable y(t) has a viscous, drift dynamics [47], i.e. φ = L c (y)I c , where the function f (y) depends on the physical realization of a meminductor and corresponding system parameters. Inspired by the ferromagnetic-rod example, we consider L c (y) = yL > + (1 − y)L < where y(t) is the fractional size of the "effective magnetic core" (Fig. 1d). One microscopic mechanism for generating a current-induced spin polarization (or magnetization) is the spin Hall magnetoresistance effect [50,51], that enabled the realization of the meminductor in a platinum yttrium-iron-garnet (Pt/YIG) hybrid structure [52]. We note that since Eq. (16) relates the flux to the current, the Kirchhoff-law equation for the current through the inductor is modified to where ∆L = L > − L < is the maximum change in the inductance. It is worth its while to emphasize that the second, I c -dependent term in Eq. (18) is absent for an inductor without memory. Its presence changes the sign of the gauge term introduced by a time-dependent change of basis, Eq. (6), giving rise to a new, effective Hamiltonian H eff = H eff (µ(y)) − (i/2)∂ t ln A(y). Thus, to investigate energy dynamics in a meminductive dimer, we solve the following set of coupled, nonlinear equations, where η is the polarity of the meminductor, F (y) is a window function that suppresses the change in y(t) near its fixed points, and Q c is the material-dependent characteristic charge that generates sufficient spin accumulation to change y from zero to unity. The corresponding characteristic current-scale for the coupling meminductor is given by i 0 = ω 0 Q c . This phenomenological drift model, Eq.(20), produces key meminductor features such as a pinched loop hysteresis in the φ − I c plot for an alternating current input [47,52]. For a pair of ±RLC circuits with variable coupling inductance L c , the dimensionless gain-loss strength is given by Γ = γ/ω 0 = L/CR 2 . The threshold coupling at which the memory-less PT -symmetric dimer transitions from PT -broken region at µ = 0 to a PT -symmetric region is given by When the meminductance value changes from maximum to minimum, the coupling increases potentially pushing the system into the PT symmetric region; on the flip side, the reduced coupling may push the system into the PTbroken region. We investigate the behavior of the system for different meminductor strengths. The dimensionless gain-loss strength in each circuit is Γ = 0.5, the initial magnetic-core fraction is given by y(0) = 0.5, and the initial state vector |φ(0) is given by Note that we specify |φ(0) instead of the energy-density state vector |ψ(0) because the latter also depends upon y(0), i.e. the initial coupling meminductor value. Figure 5a shows the transition from a PT -broken phase to a PT -symmetric phase that occurs when the coupling strength is increased; the vertical axis is logarithmic. Figure 5b shows results for µ > = 1.3µ PT and µ < = {1.4, 1.5}µ PT . When the initial state is |χ 1 , as is seen in Fig. 5a, the system is in the PT -symmetric phase with bounded oscillations for the scaled circuit energy. In contrast, when the initial meminductor current is reversed, i.e the initial state is given by −|χ 1 , the system goes into a PT -symmetry broken state with exponential growth for the scaled circuit energy. The inset shows dynamics of the internal state variable y(t). The key difference between oscillatory behavior and exponential growth is that for the latter, y(t) switches between the two extrema in an almost square-wave fashion. Although we have shown only two instances of this highly unusual behavior, it is generically found over wide parameter ranges. It is solely due to the internal dynamics, or memory, of the coupling inductor that the fate of the system depends on the sign of the initial state, with states ±|χ 1 leading to PT -symmetric and PT -broken states, respectively. Figure 5c shows a typical instance where the circuit energy dynamics are stabilized by increasing the initial meminductor current. The circuit parameters are Γ = 0.5, y(0) = 0.5, and (µ > , µ < ) = (1.1, 1.3)µ PT . We see that the exponential growth for the initial state |χ 1 changes to an oscillatory energy dynamics for the state 100|χ 1 . This qualitative change in the dynamics is reflected in the y(t) dynamics (inset), where a square-wave modulation between extremum values corresponds to the PTsymmetry broken state whereas oscillations that do not reach both extrema show PT -symmetric phase. Another example of stabilization is shown in Fig. 5d; here, the circuit parameters are (µ > , µ < ) = (1.5, 1.7)µ PT . For initial state |χ 1 , the scaled circuit energy decays before settling into a constant-amplitude oscillatory behavior. When the initial state is changed to 50|χ 1 , that decay is arrested and the system settles into an oscillatory dynamics. The inset shows that both amplifying and decaying cases have internal-state variable y(t) that switches between the two extrema.
Lastly, we obtain the landscape of PT -symmetric and PT -broken phases via the amplification factor Λ amp in the µ < − µ > plane for initial states ±|χ 1 . The relevant circuit parameters are Γ = 0.5 and y(0) = 0.5. The top plane in Fig. 5e shows that at small couplings µ <,> /µ PT 1, the system is in the PT -broken phase. As both couplings are increased, emergence of PT -  5}µPT; Γ = 0.5, y(0) = 0.5, and initial state is |χ1 . (b) With initial state −|χ1 , the scaled circuit energy dynamics changes from bounded oscillatory behavior to an exponential growth; Γ = 0.5, y(0) = 0.5, and µ> = 1.3µPT. Inset: y(t) oscillates in a square-wave fashion between its extremum values in the PT -broken phase, whereas it does not fall below 0.5 in the PT -symmetric state. The fast oscillations with period T0 = 2π/ω0 on top of the slow dynamics are present in both phases. (c) Exponential growth of the scaled energy for initial state |χ1 is stabilized to oscillatory behavior when the initial state is changed to 100|χ1 , i.e. the initial meminductor current is increased 100-fold; Γ = 0.5, y(0) = 0.5, and (µ>, µ<) = (1.1, 1.3)µPT. (d) For (µ>, µ<) = (1.5, 1.7)µPT, and initial state |χ1 , the circuit energy decays before settling into bounded, oscillatory behavior. This decay is arrested for an initial state 50χ1 . The insets in (c)-(d) show snapshots of the y(t) dynamics. (e) Amplification factor Λamp(µ>, µ<) in the lower-half plane (µ< greater than µ>) for initial states ±|χ shows differences seen in (b). symmetric phase is signaled by Λ amp = 0. However, this trend is not monotonic. When the weakest coupling µ > exceeds µ PT , Eq.(21), the system shows a re-entrant PTsymmetry broken phase. The bottom plane shows corresponding results for a system with initial state −|χ 1 . We see that in the region µ <,> /µ PT 1, changing the sign of the initial meminductor current does not change the fate of the system. On the other hand, some regions at moderate coupling values exhibit a change from exponential to oscillatory circuit energy dynamics. Since these results depend upon the internal state-variable value y(0) and the initial circuit energy, Fig. 5e provides only a glimpse of the rich and diverse amplification-factor landscape for a meminductive PT -symmetric dimer.
The results in this section, too, open up many more questions than they answer. The nonlinear, initial-state sign and strength dependent dynamics of such a system, combined with the multi-dimensional relevant parameter space, make it hard to obtain significant analytical results or straightforward insights into the long-term dynamics. Our results suggest that a dynamical systems theory approach may be required to understand the longterm temporal behavior of PT -symmetric systems with memory.

VI. DISCUSSION
During the past decade, open, classical systems with balanced gain and loss have seen an explosion of interest. This interest has been driven by their counterintuitive behavior, their diverse experimental realizations that span decades in relevant space and time scales, and the rich landscape of their properties that emerges when simple, canonical, static models are generalized to include timeperiodicity (Floquet), time delay, noise, correlations, and nonlinearities. Here, we have presented a new paradigm for PT symmetric systems. By introducing memory in a physically meaningful and achievable manner, we have investigated the dynamics of a PT -symmetric electric dimer. The nonlinearity introduced by the internal state variable that instills the memory, we find, leads to circuit energy dynamics that depend on both the strength and the sign of the initial state. Surprisingly, a PTsymmetric electric dimer with memristive gain and loss shows self-organized Floquet dynamics that lead to a PT -broken phase at small gain-loss strength and large coupling. Similar results are obtained when the coupling between the gain and loss LC circuits have memory. It is worthwhile to point out that the energy dynamics' sensitivity to the sign of the initial state (Fig. 5b,e) is most unusual. Typical models for the Schrodinger equation have nonlinearities that depend on the absolute value of elements in the state vector (or the wave function).
Our results expand the pool of simple, canonical, PT -symmetric models.
The notion of memory or non-Markovian behavior arises across diverse platforms, both classical and quantum. In particular, some non-Markovian aspects of quantum information and its flow between the system and its environment have been studied in dissipative [53] or PT -symmetric quantum models with a static Hamiltonian [54][55][56]. Our work, on the other hand, introduces memory into the effective non-Hermitian Hamiltonian, and leads to nonlinear and sign-dependent effects that are absent in aforementioned works. Our memory mechanism is implemented through the internal-state-variable dynamics, a non-local-in-time memory kernel, or a process tensor. Our work suggests that adding the non-Markovian aspect to PT -symmetric systems will lead to non-trivial, unanticipated results.