Intrinsic gradient oxygen-driven second-order memristors for continual reinforcement learning

Domain: chemistry
Published: 2026-03-03
Editorial outcome: minor_revision

Abstract
The intrinsic gradient in organisms governs long-term information processing and continual learning in dynamic environments. Reinforcement learning (RL) aims to emulate such temporal regulation to enhance learning efficiency and adaptability. However, the absence of intrinsic gradient construction in existing memristive devices leads to stochastic and abrupt state changes, disrupting the generation of temporally correlated internal states critical for continual RL. In this work, a second-order memristor incorporating a stable intrinsic oxygen gradient was designed via a molecular-coordinated layer, which enables a prolonged dynamic barrier evolution (>102 s). This slow dynamic response facilitates balanced oxygen ion migration and diffusion under unipolar spike stimulation, resulting in a significant conductance modulation (ΔG = –98.1%). These temporally adaptive conductance states were quantitatively mapped to learning rates in the RL algorithm, allowing the learning task timescale to co-evolve with device dynamics. Compared with conventional strategies, intrinsic-gradient-driven modulation reduced training iterations by 68.75% and 35.65% in static and dynamic environments, respectively. These findings underscore the potential of slow-dynamic second-order memristors as physically grounded time-adaptive units bridging device dynamics and algorithmic learning in neuromorphic systems.

Introduction
Continual learning, a fundamental cognitive ability in humans, enables the acquisition, refinement, and transfer of knowledge and experiences throughout an individual's lifespan.[1] This sophisticated process is underpinned by intricate in vivo electrophysiological mechanisms that establish gradients across intra- and extracellular membranes, thereby facilitating memory consolidation and retrieval processes. The significance of continual learning extends to artificial intelligent systems operating in dynamic environments, where the capacity to process continuous information streams is paramount.[2,3] Reinforcement learning (RL) algorithms are considered one of the most promising candidates for continual learning realization in computational systems.[4-6]

The emergence of solid-state and fluidic ionic memristors that emulate biological ion-based signal processing has further advanced the potential of RL in this domain.[7-14] However, persistent challenges remain, particularly in handling non-stationary data distributions that often result in convergence instability, non-monotonic reward accumulation, and redundant learning in continual RL systems. At the device level, the inherent stochasticity and fluctuations in ionic redistribution gradients often result in abrupt switching and unstable conductance, hindering the retention of temporally correlated internal states.[15,16]

In contrast, second-order memristors demonstrate distinct advantages through their biorealistic gradient modulation capabilities, mediated by multiple state variables.[17,18] These devices exhibit enhanced spatiotemporal plasticity compared to first-order memristors, owing to their time-dependent spontaneous decay processes.[19,20] However, the rapid diminishing of these state variables, like ultrafast temperature decrease or metal ion diffusion occurring within hundreds of nanoseconds, leads to a fast decay. Though advantageous for high-speed operations, such rapid decay hinders the retention of temporal conductance correlations, unlike nonvolatile memristors.[21-25] Consequently, second-order memristors struggle to process continual information over extended timescales, especially in RL tasks that require stable memory traces for long-term strategy formulation.[26] To enable continual RL with second-order memristors, it is crucial to construct long-term state variable changes to achieve a stable and continuous conductance modulation process.

Notably, the dynamic barrier height (Φ) state variable exhibits relatively slow decay, spanning milliseconds to seconds. This behavior arises from interfacial barrier recovery processes, which are predominantly governed by oxygen concentration gradients accumulated during historical operations. While pronounced gradients facilitate second-order behavior, they simultaneously accelerate recovery kinetics, thereby reducing the temporal window for interfacial barrier evolution. This accelerated recovery undermines the temporal consistency essential for continual RL, resulting in volatile conductance states that compromise long-term stability.[27-29] Although intrinsic concentration gradients mitigate abrupt oxygen redistribution and prolong barrier-height modulation, their long-term stability remains limited by operation-induced deformation in the absence of a stabilizing mechanism.[30-33] Therefore, stabilizing intrinsic concentration gradients to extend operation temporal window of Φ is crucial for achieving longer gradient-based conductance change, which stabilizes the time-varying state relevance and enhances the feasibility of continual RL. Through biomimetic gradient modulation, second-order memristors enable efficient execution of long-term tasks with reduced training iterations, highlighting their promise in next-generation continual learning systems.

In this work, the second-order memristor with ITO/zinc-porphyrin (ZnTPP)/ALD-AlO x /Al was fabricated. The incorporation of the ZnTPP layer was designed to establish stable intrinsic oxygen gradients, thereby prolonging dynamic interfacial barrier modulation (>102 s) and enabling long-term spontaneous decay characteristics. Cooperating with the molecular-coordinated transport process, the unipolar spike voltage-dependent plasticity (U-SVDP) was utilized to introduce the different degrees of oxygen ion gradient through versatile voltage amplitudes. The sustained barrier modulation facilitated the reconfiguration of oxygen ion migration and diffusion through U-SVDP, yielding 40 distinct pseudo-nonvolatile (PNV) conductance states. This intrinsic-gradient driven approach enabled dynamic learning rate updates (from ca. 0.14 to 0.77) for RL tasks, achieving a 68.75% reduction in training iterations. Moreover, the required convergence episodes were reduced by approximately 35.65% compared to conventional methods in the dynamic environments.

Results
Design of the intrinsic oxygen gradient for second-order memristor

Fig. 1 illustrates the conceptual and experimental framework underlying the intrinsic gradient-driven learning paradigm. As shown in Fig. 1a, the design is inspired by biological systems, where intrinsic gradients, such as ion concentration differences across cell membranes, govern long-term adaptation and signal regulation. The inset highlights the directional concentration change from the intracellular to the extracellular, representing the natural driving force behind gradient-based modulation. Based on this principle, Fig. 1b illustrates an RL agent whose learning rate is adaptively modulated over time in response to dynamic environmental conditions. This modulation is driven by the physical state evolution of the underlying memristive device. The device configuration, shown in Fig. 1c, comprises an ITO/ZnTPP/ALD-AlO x /Al stacking, which supports the formation of an intrinsic oxygen gradient. As detailed in Fig.1d, this gradient is induced by the molecular coordination effect of ZnTPP and is spatially fixed through a consecutive atomic layer deposition (ALD) process (Fig. 1e and Methods). The resulting intrinsic gradient significantly prolongs the dynamic barrier recovery process (Fig. 1f), leading to a slow second-order conductance decay behavior with characteristic timescales exceeding 102 seconds (Fig. 1g). Such slow, device-driven dynamics can be further modulated through unipolar spike voltage amplitude, as illustrated in Fig. 1h, enabling the generation of well-separated pseudo-nonvolatile conductance states. These states offer temporally rich and non-stationary input signals to the agent, which are essential for long-timescale adaptation in continual learning scenarios.

Firstly, the critical role of the molecular-coordinated layer in the device structure was discussed. In general, to ensure uniform growth during the subsequent ALD process, the substrate surface must possess sufficient hydrophilicity and specific functional groups to provide reactive sites and facilitate the adsorption and reaction of precursor molecules.[34,35] However, ZnTPP lacked the hydrophilic functional groups necessary for bonding with the precursors. Thus an in-depth investigation into the ALD growth mechanism was performed. For XPS depth profiling, a 3 nm-AlO x film was selected to deduce the chemisorption sites of ZnTPP through binding energy shifts (Supplementary Fig. 4a-b). As shown in Supplementary Fig. 5a, the Zn 2p core level spectra were located at 1044.85 eV and 1021.85 eV, respectively. Following ion beam sputtering, a 0.80 eV red shift could be observed. It showed that electron transfer occurred between the Zn element and the precursors, leading to changes in the chemical bond state. In contrast to the original energy spectra of Al 2p at 74.60 eV, the etched Al 2p energy spectra in Supplementary Fig. 5b were located at 74.90 eV. These results signified that AlO x with a higher degree of oxidation was generated in the filling stage.[36] In Supplementary Fig. 5c, the 532.50 eV and 531.30 eV of O 1s spectra corresponded to the binding energy of Zn-OH and Al-O could be observed on the unetched film, respectively. After etching, the asymmetric peaks of O 1s spectra located at 533.1 eV and 531.74 eV were attributed to H-O-H and Al-O bonds. These findings suggested that the ZnTPP film layer mildly adsorbed trace amounts of H 2 O molecules from the environment. At elevated temperatures, hydrogen bonds were disrupted, allowing Zn atoms to coordinate with hydroxyl groups, leading to the formation of zinc-based hydroxides, which acted as chemisorption sites.[37] Moreover, the rough ITO substrate introduced additional reactive sites (Supplementary Fig. 4c-g), facilitating increased oxygen incorporation during the initial filling stage, leading to a region with higher oxygen content. The calculated O/Al ratios in the ALD-AlO x films were 1.33 and 1.15 during the filling and conformal growth stages, respectively, remaining stable throughout the 10 nm growth (Supplementary Note 1, Table 1 and Supplementary Fig.6a-b). Compared to conventional substrates (e.g., silicon), the ZnTPP film promoted a higher oxygen content at the ZnTPP/ALD-AlO x interface, facilitating ion redistribution under an external electric field (Supplementary Fig. 6c-d).[38]

Intrinsic oxygen gradient-based second-order memristive characteristic and U-SVDP

As shown in Fig. 2a, the device exhibited typical current-voltage (I-V) characteristics under continuous sweeps of 0→6→0 V and 0→9→0 V, with clockwise current decrease and counterclockwise increase, respectively. An external bias was applied to the bottom ITO electrode, while the top Al electrode was grounded. The device uniformity was counted in Supplementary Fig. 7-8. These results indicate that excitation/inhibition behavior can be tuned by adjusting the voltage sweep range. A voltage amplitude-dependent transition was observed between 5-10 V, with 7 V identified as the critical point for switching between inhibition and excitation (Fig. 2b and Supplementary Fig. 9). This suggested different conduction mechanisms dominate at low voltage (V Low =5/6/7 V) and high voltage (V High =8/9/10 V), respectively. Additionally, a typical hysteresis curve appeared under bipolar operation (Supplementary Fig. 10), confirming the basic memristive behavior.

To unveil the voltage amplitude-dependent conduction mechanism, I-V curves fitting were performed, as shown in Fig. 2c-d. When the ITO electrode was forward-biased in the low voltage region (0-3 V), a positive slope ln(I) versus V 1/2 fitting curve could be fitted, manifesting the Schottky emission dominated carrier transport, the majority of charges primarily need to overcome the Schottky barrier generated by the Al/AlO x metal-semiconductor contact in the initial stage.[39] As the voltage increased from 3 to 6 V, distinct positive and negative slopes were identified in the ln(I/V 2 ) versus V −1 plot. The tunneling junction further transitioned to triangular, indicating a change from the direct tunneling (DT) to Fowler-Nordheim tunneling (F-N-T).[40] The larger slope observed during the backward sweep indicated a higher energy barrier, which caused the current to gradually decrease with successive scans (Supplementary Fig. 11).[41] Only a bit of oxygen ions (O 2− ) could be extracted from AlO x and bind to ZnTPP molecular through Zn-O coordination, exhibiting shallow trap and transport (Supplementary Fig. 12). When under the large voltage (V High ≥8 V), a large number of O 2− would penetrate into the ZnTPP layer from the AlO x layer, and successively hop from one ZnTPP molecule to another through binding and dissociation (Supplementary Fig. 13).[42] Therefore, the field-dependent ZnTPP coordination layer acts as a stabilizing electronic foundation, enabling a deterministic coupling between ionic drift and interfacial-barrier recovery. This conclusion is supported by thickness-dependent modulation of both the ZnTPP and ALD-AlO x layers and by comparative single layer and stacking-order control experiments (Supplementary Figs. 14-17).

The spatial-temporal dependence of ion dynamics between synapses is fundamental to the brain for information processing and communication.[43] Supplementary Fig.18a depicted the configuration relationship between the biological synapse and the memristor. To verify the impact of voltage amplitude on internal dynamics, synapse characteristics were examined. Supplementary Fig.18b exhibited clear inhibitory/excitatory transient conductance responses with demarcation under a pulse (V Low =5/6/7 V and V High =8/9/10 V) with a duration time (t d ) of 10 ms (Supplementary Fig.18c-e). The strong electric field was derived from the extremely low diffusion coefficient and high activation energy requirement of O 2− in AlO x .[44] In Supplementary Fig.18f, the continuous conductance changes under V Low or V High suggested a progressive hindrance or transfer process of O 2–. And the pulse interval (Δt)-based spike-rate-dependent plasticity (SRDP) was demonstrated as well (Supplementary Fig.18g). The higher frequency pulses of V High or V Low showed a more pronounced excitation/inhibition trend, which was determined by the dominant conduction process. Unlike the static storage characteristic, the spontaneous decay process could be achieved with the appropriate amount of O 2− driven by a 9 V pulse (Fig. 2e). This process could be well-fitted by a double exponential decay function as follow:

I(t)=I 0 +A 1 e(−t/τ 1 ) +A 2 e(−t/τ 2 ) (1)

where I 0 is the current offset, A 1 and A 2 are the constants, τ 1 and τ 2 are the decay time constants. The time constants of τ 1 (3.93 s) and τ 2 (53.60 s) obtained from the fitting result correspond to the faster and slower oxygen diffusion-back process, respectively.[29] This process demonstrated the second-order characteristic.[27] The reproducibility and pulse-width-dependent tunability of this second-order memristive behavior are confirmed in Supplementary Figs. 19-20. History-dependent SRDP was investigated (Fig. 2f), revealing frequency-dependent O 2– migration-diffusion behaviors. Under 50 Hz pulses, conductance surged due to O 2– drift against the intrinsic gradient, while 9.09 Hz pulses led to a diffusion-dominated conductance decrease. At 0.99 Hz, a balance between diffusion and drift resulted in negligible conductance changes. Reapplying 9.09 Hz pulses induced drift-driven O 2– accumulation, which enhanced the conductance. In contrast, V Low pulses consistently suppressed conductance changes, as the low electric field hindered electrons from overcoming the potential barrier, impeding O 2– accumulation and diffusion. By utilizing the single channel inputs with relative time difference (Supplementary Fig.21), the anti-Hebbian learning rule was implemented as a function of the Δt, as shown in Fig. 2g. This behavior can be explained by the inability of conductance accumulation to counteract decay effects when larger Δt pair-pulses were applied. Meanwhile, due to the non-overlapping pulse scheme, conductance was more dependent on history rather than the instantaneous sum of pulse amplitudes.[17]

By utilizing molecules to coordinate ion transport, O 2− could be controlled in conjunction with the intrinsic gradient in the ALD-AlO x layer.[45] As shown in Fig. 2h and i, the unipolar-spike-voltage-dependent plasticity (U-SVDP) was implemented using a pair of pulses with the same polarity but different amplitudes (V High and V Low ).[46] The conductance state set by V High could be modulated by V Low , rather than relying solely on O 2− diffusion (Fig. 2j). To evaluate the degree of modulation of conductance states based on the U-SVDP characteristic, ΔG U-SVDP was defined as follows (Fig. 2k):

ΔG U-SVDP =(G 2 −G 1 )/G 1 ×100% (2)

Where G 1 and G 2 were the conductance states before and after the U-SVDP modulation, respectively. There existed an enhanced inhibition process when V Low were 6 V and 7 V. The V High established a pronounced O 2− concentration gradient and metastable migration state, upon which the subsequent V Low (Δt ≤100 ms) facilitate directional migration and accumulation of O 2− along the pre-established low-barrier pathways. This resulted in a enhancement of conductance modulation, embodying a characteristic cooperative migration behavior. When the U-SVDP operation resembled the read operation (e.g., 9 V/2 V), ΔG U-SVDP followed an exponential decay process, as shown in Supplementary Fig. 22. This unipolar modulation synchronizes the fast ionic drift with the slow interfacial-barrier relaxation under the guidance of the intrinsic oxygen concentration gradient, forming a dual-variable interaction that dictates the device’s temporal conductance dynamics and prepares the ground for higher-order metastable behavior.

Intrinsic gradient-driven oxygen dynamics and physical verification

To shed light on the internal dynamics of the second-order memristor, the memristive mechanism schematic was proposed. For fresh devices, the O 2− concentration gradually decreased from the ZnTPP interface to the AlO x layer (Fig. 3a). Upon applying V High , O 2− overcame the Al-O Coulombic force and migrated via the ZnTPP network, establishing the interfacial electric field (𝐄⃑ in ) (Fig. 3b-c). The hopping transport process of the O 2− in the ZnTPP matrix were demonstrated in our previous work.[42] Multiple oxygen-deficient regions generated oxygen vacancies (V o ) filaments in the oxide layer, while the oxygen-rich region at the ZnTPP/ALD-AlO x interface acted as a buffer layer to suppress bulk filament growth.[47,48] After removing V High , part of the O 2− diffused back along the intrinsic gradient. Retention tests indicated that recovery time increased with stronger extra gradients induced by longer V High durations (Supplementary Fig. 23). The direction of 𝐄⃑ in was opposite to the applied external field (𝐄⃑ ext ), such that a V Low could only drive limited O 2− migration. This allowed partial filling of V o sites and dissolved filaments (Fig. 3d). Subsequently, the 𝐄⃑ in would further enhance as the progressive accumulation of O 2− at the interface (Fig. 3e).

The dynamic behavior of O 2− ions was simulated using finite-element simulation, with model details provided in Supplementary Note 2 and Supplementary Fig. 24. Under applying a 100 ms 𝐄⃑ ext , regions with higher initial oxygen concentration become further enriched, while depleted regions become more deficient, forming a strongly asymmetric ionic profile that stores additional potential energy in the built-in landscape. (Fig. 3f-g and Supplementary Fig. 25e). After removing the 𝐄⃑ ext , the intrinsic gradient provides a kinetic buffering effect that slows down the subsequent relaxation (Fig. 3h and Supplementary Fig. 25a-d). As shown in Supplementary Fig. 25e-g, the conductance does not drop immediately but enters a gradual transition stage lasting for about 30 s. During this period, the steep interfacial gradient is slowly released, and the residual built-in potential continues to regulate the ionic motion, maintaining a non-equilibrium but stable transport pathway. It could be observed that the concentration evolution process was still performed even after hundreds of seconds (Supplementary Fig. 25f-k and Fig. 3i). This buffering effect effectively converts the rapid diffusion process into an operable slow dynamic that preserves partial ionic activity for an extended second-order variable operation timescale, which is different from the homogeneous distribution (Supplementary Fig. 26).

In order to intuitively elucidate the physical dynamics of O 2− in the ALD-AlO x layer, in-situ c-AFM measurements were performed (Fig. 3j). A 6 V bias applied via the ITO electrode induced uniform tunneling current (~20 pA), which rapidly returned to the initial state after bias removal (Supplementary Fig. 27a-c). When increased to 9 V bias, the significant redistribution of O 2− generated multiple V o filaments whose locations correlated with reactive sites at ZnTPP/AlO x interface (Supplementary Fig. 27d-e).[49] This indicated an interwoven, rather than sharply defined, interface region. After filament formation, applying 6 V could partially dissolve filaments, confirming U-SVDP functionality (Fig. 3k-m). The residual filaments contributed to stable PNV conductance states. Moreover, current measurements at 5 and 10 minutes after bias removal confirmed second-order behavior with a prolonged recovery time (>102 s) compared to other second-order memristor (Supplementary Fig. 27f-g and Supplementary Table 2).[27,50] This slow dynamic, governed by the intrinsic gradient, enables fine-tuned conductance modulation via U-SVDP and supports time-dependent learning processes with improved adaptability.

Reconfigured PNV conductance states achieved by bioinspired U-SVDP modulation

In biological systems, continual learning relies not only on memory retention but also on the ability to modulate responses over time. A key contributor to this dynamic adaptability is the operation of non-voltage-gated ion channels (NVGCs), as illustrated in Fig. 4a.[51,52] Unlike action potential-driven channels, NVGCs remain responsive to gradual membrane potential changes and utilize intrinsic ionic gradients to regulate ion flux across membranes. This enables fine-tuned, amplitude-dependent ion transport, allowing cells to continuously update their internal state without abrupt switching. Such a mechanism is central to continual learning and long-term adaptation. Mimicking this biological principle, memristor can employ U-SVDP to regulate O 2‒ migration-diffusion along a pre-established intrinsic concentration gradient (Fig. 4b-c). Like NVGC’s principle, U-SVDP facilitates amplitude-tunable and gradient-guided ionic behavior, enabling gradual, history-dependent conductance evolution. This allows the device to physically encode time-resolved internal states that mirror the biological temporal window for adaptation.

Inspired by the NVGCs, U-SVDP-based modulation was further explored by pulse patterning. As shown in Fig. 4d, the applied V High determined the quantity of filaments and thus the conductance state for V Low operation. Different durations of V High input (t High =1000 or 100 ms) followed by V Low modulation (t Low =10 ms) resulted in a conductance change of −2.11 μS and −2.06 μS at the 4s reading, respectively (Fig. 4e-f). This suggested that under the premise of the intrinsic gradient, the pronounced O 2− redistribution stimulated by V High further changed the oxygen concentration gradient. The 1000 ms V High induced large extra gradient changes that caused a longer diffusive process, making V Low operation obtain an unstable state. By contrast, the 100 ms V High stimulation resulted in a slower gradient recovery, allowing V Low to stabilize the decay conductance into a pseudo-nonvolatile (PNV) state. The PNV state refers to metastable, temporally distinguishable conductance plateaus that persist for a finite period after the stimulus is removal, and this state originates from the coupled dynamics of the interfacial barrier and the multi-filament network. Unlike conventional relaxation-driven decay, the PNV state represents a transiently frozen equilibrium between adjacent energy minima, arising from the asynchronous relaxation of multiple conductive filaments and the slow recovery of the interfacial barrier. However, under a lower additional gradient induced by 10 ms V High stimulation, the conductance changes (–0.90 μS and –1.04 μS, respectively) remained similar regardless of the V Low duration (Fig. 4g-h), due to the limited diffusible O 2−. To achieve a suitable additional gradient for slow gradient recovery and PNV state promotion, both t High and t Low were optimized to 100 ms.

As a small amount of O 2− driven by consecutive V Low pulses increased the interfacial barrier and partially dissolved the existing V o multi-filament (Fig. 4i). As a result, the conductance change (ΔG U-SVDP =−72.4%) after 16 successive PNV states was significantly greater than the −51.7% spontaneous decay (ΔG) (Fig. 4j). Furthermore, the stability of the intrinsic-gradient-driven U-SVDP modulation PNV states were demonstrated in Supplementary Fig. 28. After a varying number of V Low stimulations within a defined period, the PNV state remained stable (Fig. 4k). Meanwhile, a small amount of O 2− accumulated by V Low would decay under the intrinsic gradient after removing the stimulus (Supplementary Fig.29). These results indicated that the multi-filament system supported by the intrinsic gradient dissolved partially, providing a basis for achieving multi-PNV conductance states.[53] Inspired by non-voltage-gated ion channels, which rely on weak yet persistent ionic modulations to regulate long-term biological functions, the device operates beyond binary switching and supports the formation of multi-level PNV conductance states governed by intrinsic ionic dynamics. Importantly, these PNV states are not static memory states but arise from gradient-driven transient and relaxation processes, whose stability and repeatability are experimentally verified under repeated electrical cycling and thermal perturbation (Supplementary Figs. 30-34).

Intrinsic oxygen gradient-driven dynamic learning rates update for continual RL

In conventional RL, an agent aims to discover a policy that maximizes long-term reward.[54,55] Continual RL, however, requires the agent to adapt continuously without forgetting prior knowledge. Although the learning rate (α) can be manually tuned, issues like slow convergence and local optima remain. By leveraging device-level time-dependent states, learning rates can be dynamically updated to support continual adaptation. To demonstrate this, a continual RL model was applied to an optimal pathfinding task for autonomous vehicles within signal coverage (Fig. 5a-b). The Q-learning-based continual RL was illustrated in Fig. 5c and was explained in depth in Supplementary Note 3. The ΔG (−98.1%) exhibited by 40 PNV conductance states were clearly distinguishable from spontaneous decay (Fig. 5d). Their low variation confirmed the stable multi-filament system (Fig. 5e). In this system, the learning rate α corresponds to the time-dependent change in memristor weight, defined as:

α (t) = (ω t+1 − ω t ) / ω t × 100% (3)

where t is the time after V High stimulation and ω is the conductance at the corresponding time t. The ω 0 corresponding to t 0 was the initial state after the V High (t=1.135 s). Notably, t 1 (7.658 s) marked the rapid O 2− diffusion phase (τ 1 ), during which spontaneous decay yielded the largest ΔG. In contrast, ΔG U–SVDP continued to increase beyond t 1 due to ongoing O 2− redistribution triggered by V Low pulses (Fig. 5f-g). With equal time intervals between t 1 , t 2 , t 3 , and t 4 (each 60 s), the learning rates α U–SVDP and α Decay exhibited divergent temporal trends, as shown in Fig. 5k. The optimal path derived from Q-learning was shown in Fig. 5h. The total reward received and the exploration steps spent in each iteration epoch with different α s, as shown in Fig. 5i-j and Supplementary Fig. 35. When the α U-SVDP dynamic increased from 0.14218 to 0.77311, the number of iterations required for convergence decreased from 3927 to 1227, leading to a 68.75% reduction. On the contrary, more iterations was required under α Decay s modulation, as shown in Fig. 5k and Supplementary Fig. 36. The standard deviation was used to assess the stability of the convergence results as well. It was observed that the α U-SVDP -based results gradually approached the optimum and exhibited the most stable performance (Fig. 5l-m). For comparison, linear α manual modulation exhibited local optima and convergence oscillations (Supplementary Fig. 37).

To evaluate whether the learning-rate schedule reconstructed from intrinsic-gradient-driven U-SVDP conductance evolution remains effective beyond static settings, the continual RL task was extended to dynamic multi-stage environments. The U-SVDP-induced conductance evolution was fitted to a continuous nonlinear function α U-SVDP (p), which was then applied as the learning-rate schedule in the per-step Q-learning update (Supplementary Note 4 and Supplementary Fig. 38). The linear and fixed learning rate strategies were also implemented for comparison. As shown in Fig. 6a, the α U-SVDP (p), α Linear (p), and α Fixed modulation strategies were applied across four consecutive stages. In this scenario, the agent sequentially navigates through environments of increasing complexity, represented by map sizes of 50×50, 75×75, 100×100, and 200×200, to mimic the challenges of continual RL in evolving real-world tasks (as shown in Fig. 6b). The optimal trajectories illustrated that α U-SVDP enables more flexible policy adaptation as the environment scales (Fig. 6c). The convergence curves in Fig. 6d further showed that α U-SVDP consistently accelerate learning across stages, achieving shorter converged paths and enhanced stability relative to both α Linear and α Fixed (Supplementary Figs. 39-41). Compared with traditional method (linear modulation), the α U-SVDP (p) reduced the total training episodes of the agent by 35.65%. The optimal path results were compared across map sizes in detail, demonstrating the effectiveness of this device-derived learning-rate modulation (Fig. 6e-h). Beyond map-size scaling, additional evaluations were performed under progressively increasing obstacle densities, multiple reward structures, and a fully non-stationary dynamic-goal environment, all of which consistently confirmed the robustness and adaptability of the α U-SVDP modulation (Supplementary Figs. 42-53). These results underscore that the nonlinear learning-rate trajectory derived from intrinsic-gradient-mediated and slowly relaxing ionic dynamics provides an effective physical substrate for dynamic learning-rate modulation, enabling efficient stage-wise adaptation in continual RL through the close coupling between device-level processes and algorithmic behavior.

Methods
Device preparation

Devices were fabricated on commercial ITO-coated 1.5 × 1.5 cm glass substrates with a sheet resistance of ≈10 Ω sq−1.The etched ITO acts as the bottom electrode, and the area of the memristor device is 0.01 mm2, as defined by 100 μm wide ITO stripes crossed over by 100 μm wide top Al stripes. ZnTPP (5,10,15,20-tetraphenyl-21H,23H-porphyrin zinc (II)) was obtained from Sigma Aldrich without further purification. The ZnTPP film with a thickness of ≈25 nm were primarily deposited through thermal vacuum evaporation method with an evaporation growth rate of ≈0.1 Å s−1. The film thickness was measured by employing the step height down mode of AFM.

The ALD-AlO x film was deposited on ZnTPP film by ALSD using trimethylaluminum (TMA) and H 2 O as precursors in a continuously evacuated vacuum chamber. The ALD pulse lengths were set using ALD series valves (Swagelock) with a microelectronic controller to 20 ms for TMA and 15 ms for water, cycle in the order of TMA (20 ms) /N 2 (15 s)/H 2 O (15 ms)/N 2 (15 s). The 10 sccm N 2 gas was used as the carrier gas in the reaction. The growth per cycle (GPC) was around 0.8~1.2 Å. The ALD deposition cycle of AlO x on the ZnTPP film at 150℃ (shown in Fig.1f). After 83 cycles, a 10 nm AlO x film was obtained on the ZnTPP film (Supplementary Fig. 1). And the conformal growth mode based on ZnTPP was qualitative analyzed in Supplementary Fig. 2-3. Finally, a top Al electrode (≈125 nm) was prepared by the conventional evaporation process (0.5 Å s−1) with a shadow mask.

Electrical characterization

All electrical measurements were performed using a Keithley2636B and FS-Pro semiconductor parameter analyzer under ambient conditions. A defined voltage bias was applied to the ITO bottom electrode while the Al top electrode grounded.

Simulation of ionic migration under intrinsic concentration gradient

Finite-element simulations were used to model oxygen-ion migration in the ALD-AlO x layer. Ionic transport was described by a one-dimensional Nernst-Planck-Poisson framework, incorporating drift under the external electric field and diffusion regulated by the intrinsic built-in potential. The intrinsic oxygen concentration gradient was extracted from depth-resolved measurements and imposed as the zero-flux equilibrium profile. A time-dependent diffusion coefficient was introduced to capture the experimentally observed transition from fast field-assisted redistribution to slow relaxation governed by the built-in potential. The simulated interfacial oxygen-ion concentration was converted to conductance using an exponential barrier-modulation relation. All governing equations, boundary conditions, and parameter values are provided in Supplementary Note 2.

RL tasks in dynamic environments

In the dynamic multi-stage task, the agent must adapt to sequentially enlarged map sizes. A physically derived, stage-wise learning-rate schedule is used to model long-term adaptability. Under the U-SVDP scheme, the device produces a set of PNV multilevel conductance states {G i }. The instantaneous conductance evolution is quantified as:

ΔG i = G i+1 − G i (4)

These gradients are normalized into a bounded learning-rate range using:

α i = α min + (ΔG i − ΔG min ) / (ΔG max − ΔG min ) × (α max − α min ) (5)

The discrete sequence {α i }is then fitted into a smooth nonlinear curve α U-SVDP (p) (as shown in Fig. 6a) and subsequently segmented into successive blocks of episodes, ensuring gradual learning-rate evolution while avoiding excessive fluctuation at every episode.

During RL training, the learning rate for all episodes within the k-th segment is set to:

α = α (p) (6)

p ∈ Φ(k) (7)

where Φ(k) denotes the k-th episode block. The same segmentation function Φ(·) and boundary values α max , α min are used to construct the linear baseline for fair comparison. This procedure ensures that the RL agent experiences a physically grounded, slowly evolving plasticity schedule consistent with the device-derived adaptation curve. Full mathematical details and implementation steps are documented in Supplementary Note 4.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper shows an interesting way to link a memristor’s slow decay behavior with the learning rate in Q-learning. By taking the decay curve from a second-order memristor and using it to set α(t), the authors manage to speed up convergence in both static and dynamic RL tasks. The experiments and modeling are solid, and the connection from the device physics to the RL performance is clear. The idea of using real device behavior to control a training parameter is nice and could be useful in neuromorphic systems. On the other hand, the RL work is done offline with data from earlier measurements, not in a real-time hardware-in-the-loop setup. That makes it harder to see how this would run in practice. Also, important RL settings like discount factor, exploration rate, and reward design aren’t fully given, and the tasks are pretty simple. It would be good to see how this works in more complex or realistic environments. Overall, the work is creative and well-executed, but it needs more details and broader testing to be fully convincing.\n- Is the tunability of the decay speed an inherent advantage of this device, and if so, how controllable and stable is this tunability across different devices or operating conditions?\n-What is the endurance (number of cycles). Can you suggest the requirements of cycling for RL in this paper? And temperature stability of the PNV states, and what is the observed variability range?\n-What were the values of the discount factor, exploration rate, and initial Q-table settings used in the RL simulations?\n-How were the reward functions and reward magnitudes determined for both the static and dynamic environments?\n-What was the exact mapping rule between the RL episode progression and the sampling points from the decay curve when generating α(t)?\n-Has the feasibility of implementing this approach in a real-time hardware-in-the-loop (HIL) setup been evaluated?\n-Beyond map size scaling, have you tested increased environment complexity through changes in obstacle density or reward structure?", "reviewer_2": "The manuscript entitled “Intrinsic Gradient Oxygen-Driven Second-Order Memristors for Continual Reinforcement Learning” proposed a second order memristor with long retention time and its application in continual reinforcement learning.\nGenerally, the manuscript would be suitable for publishing on Nature Communication after minor revision:\n1. How is the stability/durability of the device?\n2. As suggested by Fig.S25, even if the concentration is homogenous, the change of ion seems slow enough (color distribution at 1s and 50 s are obviously different) to induce the long retention time. Therefore, the intrinsic gradient may not be the key factor to achieve long retention, undermining the reliability of mechanism. To clarify this, the authors should give simulated conductance-time curves to demonstrate the sharp/smooth difference between gradient/homogeneous concentration. Simulated I-V curves may also be helpful.\n3. Would the dynamics of binding/dissociation between O2– and ZnTPP contribute to the long retention time?\n4. The new concept, pseudo-nonvolatile (PNV) states, should be further clarified. What is the difference between PNV and nonvolatile memory and volatile memory? Is there a paradigm of PNV in biological system?\nMinor issues:\n1. Page 4, Line 7, “Additionally” is redundant.\n2. Both ΔT and Δt are used to represent time interval, which could be confusing, and should be unified into Δt.\n3. In the last two graphs of fig.5e, the numbers on x-axis are overlayed.\n4. There is no Fig.3n, but was mentioned in text and the caption of Fig.3. Besides, Fig. 3j is not mentioned.\n5. Fluidic memristor as a new kind of memristor should be mentioned in the introduction part. (typical ref., science, 2023, 379, 156-161, and 161-167)"}, "author_response": "We sincerely thank the reviewers for taking the time to review our manuscript. We are very grateful for the positive assessment and constructive comments. All changes are highlighted in red in the revised manuscript and Supplementary Information.\nThe main changes we have made are:\n1) Expanded verification of device stability and durability, including long-term retention (1-21 days), 32 device-to-device variability of decay curves and PNV states, temperature-dependent stability (room temperature to 120°C), and endurance under both bipolar cycling (104 cycles) and U-SVDP cycling (5×103 cycles). These results collectively demonstrate that the intrinsic oxygen gradient remains stable and available under repeated external-field perturbations.\n2) Clarified the reinforcement learning framework and reward design through systematic evaluation across varying map sizes, obstacle densities, reward structures, and non-stationary tasks, thereby confirming that learning-rate modulation derived from intrinsic device dynamics enables consistently faster and robust convergence, even in ten-stage non-stationary scenarios.\n3) Refined physical modeling and simulations by distinguishing ionic transport under intrinsic concentration gradients from homogeneous diffusion, including time-dependent conductance analysis enabled by a time-dependent diffusion coefficient, thereby confirming that the intrinsic concentration gradient acts as the internal physical origin supporting long-timescale system dynamics.\nPlease find the enclosed point-by-point responses to the reviewers’ comments. Our responses to Reviewer 1’s comments are provided on pages 2-51, those to Reviewer 2 on pages 52-86. If any additional information is required, we will be more than happy to provide it as soon as possible.\nKind regards,\nHaifeng Ling, Linghai Xie and Wei Huang\nReviewer #1:\nComments: This paper shows an interesting way to link a memristor’s slow decay behavior with the learning rate in Q-learning. By taking the decay curve from a second-order memristor and using it to set α(t), the authors manage to speed up convergence in both static and dynamic RL tasks. The experiments and modeling are solid, and the connection from the device physics to the RL performance is clear. The idea of using real device behavior to control a training parameter is nice and could be useful in neuromorphic systems. On the other hand, the RL work is done offline with data from earlier measurements, not in a real-time hardware-in-the-loop setup. That makes it harder to see how this would run in practice. Also, important RL settings like discount factor, exploration rate, and reward design aren’t fully given, and the tasks are pretty simple. It would be good to see how this works in more complex or realistic environments. Overall, the work is creative and well-executed, but it needs more details and broader testing to be fully convincing.\nResponse: We appreciate the perspective, which involves using the intrinsic relaxation behavior of the device to directly modulate the learning rate, representing a potentially valuable approach to neuromorphic learning systems based on physical principles. In response, we have carefully addressed all specific comments by providing additional simulation results, detailed discussions on the underlying issues, and clearer figures in the manuscript and Supplementary Information. All revisions have been clearly marked in the revised manuscript, and we hope that these revisions satisfactorily address the reviewer’s concerns.\nResponse: We appreciate the reviewer for the valuable comments to improve our work.\nThe tunability of the decay speed is an inherent characteristic of the device, arising from the presence of an intrinsic concentration gradient that governs the relaxation timescale after electrical excitation. This intrinsic gradient enables the decay process to proceed through a gradual, directional evolution rather than abrupt or stochastic conductance collapse, thereby providing a physical basis for tunable relaxation dynamics.\nFirst, the device-to-device variation of the decay stability was demonstrated by statistical analysis across 32 devices. As shown in Supplementary Fig. 19a-b, the ensemble-averaged conductance decay can be well described by a bi-exponential model, yielding clearly separated fast (τ ) and slow (τ ) 1 2 1 relaxation components. Here, τ corresponds to rapid, local relaxation processes such as interfacial rearrangement or shallow trapping/detrapping, whereas τ reflects slow ion redistribution constrained by 2 the intrinsic concentration gradient. Importantly, the ratio τ /τ serves as a key metric 2 1 to quantify the relative dominance of gradient-governed diffusion over fast local relaxation. A large τ /τ ratio indicates that the relaxation dynamics are primarily 2 1 governed by slow, gradient-constrained ion migration, signifying an effective stabilization of non-equilibrium conductance states by the intrinsic gradient. In contrast, a reduced τ /τ ratio suggests that fast relaxation processes dominate, corresponding to 2 1 a diffusion-driven regime in which the intrinsic gradient no longer acts as the primary constraint. The relatively narrow distribution of the τ /τ ratio (Supplementary Fig. 19c), 2 1 confirming that both the characteristic decay speeds and their relative weighting are highly reproducible across devices, indicating that the tunability is not device-specific but intrinsic to the device architecture.\n2 1\nSecond, the controllability under different operating conditions was examined by systematically varying the programming pulse width. As shown in Supplementary Fig. 20a-b, increasing the pulse width progressively increases the number of electrically activated oxygen ions, leading to an enhanced contribution of the fast relaxation component associated with rapid diffusion. When the pulse width exceeds approximately 400 ms, the system enters a distinct regime in which excessive ionic activation perturbs the original intrinsic concentration gradient, preventing full recovery to the initial gradient profile and driving the device toward a reconfigured concentration landscape. This behavior marks the onset of a diffusion-dominated regime, in which the intrinsic gradient no longer serves as the primary constraint on the relaxation dynamics. Such operating conditions are therefore intentionally avoided in the U-SVDP scheme. Importantly, the existence of this pulse-width-dependent operational window enables deliberate adjustment of the relative contributions of fast and slow relaxation processes. This tunability allows the decay dynamics to be tailored for different task requirements and provides a physical handle to regulate the effective learning-rate strength in neuromorphic or reinforcement-learning-inspired applications, while avoiding irreversible gradient distortion.\n2 1\nThird, Supplementary Fig. 28 demonstrates that, once a stable intrinsic concentration-gradient-governed relaxation is established, a properly designed U- SVDP scheme enables reproducible modulation of pseudo-nonvolatile (PNV) states across different devices. Although the absolute conductance levels of the PNV states exhibit device-to-device variations, the temporal evolution and modulation trends remain highly consistent among devices. Importantly, each conductance change (ΔG U- SVDP ) is defined relative to the intrinsic baseline of the individual device rather than to an absolute conductance value. As a result, device-specific offsets in the PNV level do not affect the effective modulation behavior.\nThis trend-level consistency ensures that the U-SVDP-induced ΔG can be U-SVDP reliably used as a device-independent control variable, which is particularly critical for learning-rate mapping in neuromorphic or reinforcement-learning-inspired applications. In this context, the relative evolution of conductance, rather than its absolute magnitude, governs the learning dynamics, making the observed device-to-device variations in PNV states non-detrimental to system-level functionality.\nIn response to the reviewer’s concern about the decay-speed tunability, we have clarified the role of the intrinsic concentration gradient in governing the relaxation dynamics and introduced τ /τ as a quantitative metric of gradient 2 1 dominance. We further added pulse-width-dependent and device-to-device analyses to demonstrate the controllability and stability of the tunable decay behavior under the U-SVDP modulation. These revisions are included in the revised manuscript and Supplementary Figures 19-20 and 28.\nOn page 7 of the manuscript, “This process demonstrated the second-order characteristic.[27] The reproducibility and pulse-width-dependent tunability of this second-order memristive behavior are confirmed in Supplementary Figs. 19-20.”\nOn page 11 of the manuscript, “As a result, the conductance change (ΔG = U-SVDP −72.4%) after 16 successive PNV states was significantly greater than the −51.7% 6 spontaneous decay (ΔG) (Fig. 4j). Furthermore, the stability of the intrinsic-gradient-driven U-SVDP modulation PNV states were demonstrated in Supplementary Fig. 28.”\nOn page 23 of the Supplementary Information:\nSupplementary Fig. 19. Statistical analysis of second-order relaxation dynamics across 32 devices. a Mean relaxation curve averaged over 32 devices after programming with a 9 V, 100 ms voltage pulse, measured at a read voltage of 1 V. b Statistical distribution of the extracted fast (τ ) 1 and slow (τ ) bi-exponential relaxation time constants across 32 devices, presented as box plots with 2 individual data points overlaid. c Statistical distribution of the ratio τ /τ .\n2 1\nAs shown in Supplementary Fig. 19, the statistical evaluation of the decay stability across 32 devices was performed. The ensemble-averaged conductance decay can be well fitted by a bi-exponential function, yielding distinct fast (τ ) and slow (τ ) 1 2 relaxation components. Here, τ is associated with rapid local relaxation processes, 1 while τ reflects slower ion redistribution governed by the intrinsic concentration 2 gradient. The ratio τ /τ is extracted as a quantitative indicator of the relative 2 1 contribution between gradient-constrained diffusion and fast local relaxation. The relatively narrow distribution of τ /τ across devices indicates that the decay dynamics 2 1 and their relative weighting remain consistent from device to device, confirming that the relaxation behavior is intrinsic to the device architecture rather than arising from device-specific variations.\nSupplementary Fig. 20. Pulse-width-dependent relaxation dynamics and tunability of decay behavior. a Conductance decay characteristic measured after applying programming pulses with different pulse widths (10/50/100/200/300/400/500/700/1000 ms), read at 1 V. b Pulse-width- dependent evolution of the ratio τ /τ .\n2 1\nThe controllability of the decay dynamics was demonstrated under different operating conditions by systematically varying the programming pulse width. As shown in Supplementary Fig. 20a-b, increasing the pulse width activates a larger number of oxygen ions, resulting in an enhanced contribution of the fast relaxation component associated with rapid diffusion. When the pulse width exceeds approximately 400 ms, excessive ionic activation perturbs the original intrinsic concentration gradient, preventing full recovery to the initial gradient profile and driving the device toward a reconfigured concentration landscape. This behavior indicates a transition to a diffusion-dominated relaxation regime, in which the intrinsic concentration gradient no longer governs the decay dynamics. Consequently, this operating regime is excluded in the design of the U-SVDP scheme. Within the gradient-preserving pulse-width range, the relative contributions of fast and slow relaxation components can be systematically adjusted, demonstrating controllable decay behavior under pulse-width modulation.\nSupplementary Fig. 28. Statistical evaluation of U-SVDP modulation across 32 devices. a Transient conductance measured under the 100 ms V (9 V) pulse. b Transient conductance High measured under the consecutive100 ms V pulse (6 V, Δt = 5 s). c Evolution of the conductance Low states reconstructed by the V pulses (read at 1 V), corresponding to the reconfiguration of PNV Low states.\nAs shown in Supplementary Fig. 28, the device-to-device consistency of PNV state modulation was demonstrated under the U-SVDP scheme once a stable intrinsic concentration-gradient-governed relaxation is established. Although the absolute conductance levels of the PNV states vary among different devices, the temporal evolution and modulation trends remain highly consistent. Notably, the conductance change (ΔG is defined relative to the intrinsic baseline of each individual device U-SVDP) rather than an absolute conductance value. As a result, device-specific offsets in the PNV level do not affect the relative modulation behavior, confirming that the U-SVDP scheme enables reproducible and device-independent PNV state modulation across different devices.\nResponse: We appreciate the reviewer for the valuable comments. These questions address a central issue of this work, namely how cycling stability should be defined for a second-order memristive system in which the functional role arises from intrinsic-gradient-driven dynamics rather than from static memory conductance states.\nStability under conventional bipolar cycling. We first evaluate the device under conventional bipolar voltage cycling using positive and negative pulses, which is a commonly employed reliability assessment for memristive devices. Unlike nonvolatile memory applications, this test does not aim to evaluate the endurance of precisely programmed conductance states. Instead, since the device operates as a second-order memristor, the focus is on the reproducibility of the transient conductance response under repeated electrical excitation. As shown in Supplementary Fig. 30, over 104 cycles of bipolar excitation, the transient conductance response exhibits no cumulative degradation, suggesting that the device can tolerate repeated electrical perturbations while maintaining stable dynamical behavior. This measurement therefore serves as a baseline verification of device robustness, rather than a metric of memory endurance.\nSupplementary Fig. 30. Endurance characteristics evaluated using a conventional bipolar switching scheme. a Endurance performance measured under a conventional bipolar voltage switching scheme over 10,000 consecutive cycles (V = 9 V, t = 50 μs; V = ‒9 V, t = 80 positive d negative d μs, Δt = 20 μs), the interval between consecutive cycles was fixed at 150 μs. b Enlarged view of several representative cycles, illustrating the detailed conductance response to the bipolar programming pulses and the corresponding read operation (1 V).\nMore importantly, the endurance requirement in this work arises from the fundamental operating characteristics of a second-order memristor. Unlike conventional nonvolatile memristive devices, a second-order memristor does not possess well-defined static memory states that can be repeatedly programmed and retained. Instead, its functionality is encoded in transient and relaxation dynamics governed by ionic processes. Consequently, conventional endurance tests that rely on repeated programming of stable conductance states are not directly applicable. Endurance must therefore be evaluated in a manner consistent with the intrinsic operating mechanism of the device, namely the stability of the intrinsic oxygen concentration gradient that underpins the dynamic response. Based on this consideration, we designed a specialized cycling scheme to investigate whether the intrinsic gradient remains effective and directionally stable under repeated electrical disturbances, rather than evaluating the endurance of static conductance states.\nTo directly probe this aspect, we employed a U-SVDP cycling scheme\n(Supplementary Fig. 31). In this scheme, a V pulse was used to activate oxygen ions High and drive the system away from equilibrium, thereby generating a population of ions that could participate in subsequent diffusion processes. This pulse was not intended to program a stable conductance state, but rather to provide a controlled perturbation to the ionic distribution. Under short-pulse cycling (t = 50 μs), the oxygen ions activated d by V rapidly redistributed and relaxed, such that the subsequent V pulse could High Low not effectively reutilize these transiently activated ions to reconstruct the conductance state. Nevertheless, owing to the presence of the intrinsic concentration gradient, the transient response induced by V still exhibited a modulation polarity opposite to that Low induced by V .\nHigh\nIn contrast, when the pulse width of V was increased to 1 ms, a substantially High larger population of diffusible oxygen ions remained available for redistribution (Supplementary Fig. 32). In this case, the V pulse was able to actively engage these Low ions under the constraint of the intrinsic concentration gradient, leading to a more pronounced and reproducible conductance modulation with negative polarity under readout conditions. These behaviors indicate that, while the intrinsic concentration gradient dictates the directionality of the response, the availability of diffusible ions determines the magnitude and persistence of the V -induced conductance modulation.\nLow\nOverall, these observations confirm that the U-SVDP cycling scheme provides an effective means to evaluate the stability of the intrinsic concentration gradient, which constitutes the essential endurance metric for this second-order memristive system.\nSupplementary Fig. 32. Endurance performance under the U-SVDP modulation scheme. a Endurance performance evaluated using the U-SVDP modulation scheme (V = 9 V V = 6 V, High Low\n12\nt = 1 ms, Δt = 10 ms), the interval between consecutive cycles was fixed at 1 ms. b Enlarged view d of several representative cycles.\nImplications for cycling requirements in RL. The relevance and requirements of cycling endurance in RL depend strongly on the learning scenario. In short-horizon or stationary tasks, learning-rate dynamics may only need to be invoked once, placing minimal demands on cycling capability. In contrast, in continual and non-stationary RL, agents must repeatedly adapt to evolving environments, reward structures, or task objectives, requiring the learning-rate modulation to be regenerated multiple times to enable rapid re-learning followed by gradual stabilization. Similarly, in multi-stage or curriculum-based RL, each training stage benefits from a fresh learning-rate modulation cycle rather than from repeated programming of identical conductance states.\nIn this work, the learning rate is extracted from gradient-driven transient and relaxation dynamics within each cycle and is evaluated relative to the device’s instantaneous baseline rather than an absolute target conductance. Consequently, the critical endurance requirement is not ultra-high cycling capability for precise conductance programming, but the stability of the intrinsic concentration gradient and its repeated invocation as a reliable dynamical resource. From a practical perspective, this corresponds to the ability to regenerate learning-rate dynamics over hundreds to thousands of cycles, provided that the applied electrical stimulation does not induce cumulative ionic redistribution that erases the intrinsic gradient. The demonstrated robustness under both conventional bipolar cycling and asymmetric U-SVDP cycling therefore satisfies the practical cycling requirements for sustained and adaptive RL in this study.\nThermal stability of the PNV states. Within the temperature range of 40-80°C, the conductance states remain nearly unchanged, indicating robust stability of the intrinsic-gradient-governed transport behavior (Supplementary Fig. 33c). Although elevated temperature facilitates oxygen-ion activation and enhances the transient conductance response, the underlying coordination environment and gradient profile\n13\nremain largely preserved in this regime. At higher temperatures (100-120°C), the transient conductance level further increases due to enhanced ionic activity. However, the high-temperature environment progressively destabilizes the ZnTPP coordination and perturbs the multi-filament network, leading to reduced retention of the reconstructed PNV states. As a result, an increasingly pronounced conductance drift is observed at elevated temperatures (Supplementary Fig. 33a-b and Supplementary Fig. 34a-d), reflecting partial thermal redistribution of the intrinsic concentration gradient.\nSupplementary Fig. 33. Evaluation of U-SVDP-modulated PNV states at different temperatures. a,b Transient conductance responses of V (9 V, t = 100 ms) and V (6 V, t = High d Low d 100 ms, Δt = 5 s) measured at different temperatures. c Conductance evolution of the PNV states under different temperature conditions.\nIn response, we have revised the manuscript and Supplementary Information as follows:\nOn page 11 of the manuscript, “Inspired by non-voltage-gated ion channels, which rely on weak yet persistent ionic modulations to regulate long-term biological functions, the device operates beyond binary switching and supports the formation of multi-level PNV conductance states governed by intrinsic ionic dynamics. Importantly, these PNV states are not static memory states but arise from gradient-driven transient and relaxation processes, whose stability and repeatability are experimentally verified under repeated electrical cycling and thermal perturbation (Supplementary Figs. 30-34).”\nOn page 38-42 of the Supplementary Information:\nThe stability of the device under repeated electrical excitation was first examined using conventional bipolar voltage cycling with alternating positive and negative pulses. This measurement provides a baseline assessment of the robustness of the transient conductance response under repeated electrical perturbation. As shown in Supplementary Fig. 30, the transient conductance evolution remains highly reproducible over 104 consecutive cycles, with no observable cumulative degradation. The persistence of the transient response indicates that the device maintains stable dynamical behavior under repeated bipolar excitation.\nSupplementary Fig. 31. Endurance performance under the U-SVDP modulation scheme. a Endurance performance evaluated using the U-SVDP modulation scheme (V = 9 V V = 6 V, High Low t = 50 μs, Δt = 100 μs), the interval between consecutive cycles was fixed at 950 μs. b Enlarged d view of several representative cycles, illustrating the detailed conductance response to the U-SVDP\n16\nprogramming pulses and the corresponding read operation (1 V).\nTo further examine the stability of the intrinsic-gradient-governed dynamics, an asymmetric unipolar cycling scheme based on the U-SVDP scheme was employed. In this scheme, a high-amplitude pulse V was applied to activate oxygen ions High and perturb the ionic distribution away from equilibrium, followed by a low-amplitude pulse V that probes the subsequent dynamical response under gradient-constrained Low conditions.\nSupplementary Fig. 32. Endurance performance under the U-SVDP modulation scheme. a Endurance performance evaluated using the U-SVDP modulation scheme (V = 9 V V = 6 V, High Low t = 1 ms, Δt = 10 ms), the interval between consecutive cycles was fixed at 1 ms. b Enlarged view d of several representative cycles.\nUnder short-pulse excitation (t =50 μs), the oxygen ions activated by V d High rapidly redistributed and relaxed. As a result, the subsequent V pulse does not induce Low a persistent conductance reconstruction but instead produces a transient response governed by the intrinsic concentration gradient. Notably, the modulation polarity induced by V is opposite to that induced by V (Supplementary Fig. 31), indicating Low High that the response directionality is dictated by the intrinsic gradient rather than by residual field-activated ions.\nWhen the pulse width of V was increased to 1 ms, a substantially larger High population of diffusible oxygen ions remained available for redistribution (Supplementary Fig. 32). In this case, the V pulse was able to actively engage these Low ions under the constraint of the intrinsic concentration gradient, leading to a more pronounced and reproducible conductance modulation with negative polarity under readout conditions. These behaviors indicate that, while the intrinsic concentration gradient dictates the directionality of the response, the availability of diffusible ions determines the magnitude and persistence of the V -induced conductance modulation.\nLow\nOverall, these observations confirm that the U-SVDP cycling scheme provides an effective means to evaluate the stability of the intrinsic concentration gradient, which constitutes the essential endurance metric for this second-order memristive system.\nSupplementary Fig. 33. Evaluation of U-SVDP-modulated PNV states at different temperatures. a,b Transient conductance responses of V (9 V, t = 100 ms) and V (6 V, t = High d Low d 100 ms, Δt = 5 s) measured at different temperatures. c Conductance evolution of the PNV states under different temperature conditions.\nWithin the temperature range of 40-80°C, the conductance states remain nearly unchanged, indicating robust stability of the intrinsic-gradient-governed transport behavior (Supplementary Fig. 33c). Although elevated temperature facilitates oxygen-ion activation and enhances the transient conductance response, the underlying coordination environment and gradient profile remain largely preserved in this regime. At higher temperatures (100-120°C), the transient conductance level further increases due to enhanced ionic activity. However, the high-temperature environment progressively destabilizes the ZnTPP coordination and perturbs the multi-filament network, leading to reduced retention of the reconstructed PNV states. As a result, an increasingly pronounced conductance drift is observed at elevated temperatures (Supplementary Fig. 33a-b and Supplementary Fig. 34a-d), reflecting partial thermal redistribution of the intrinsic concentration gradient.\nResponse: We appreciate the reviewer for the valuable comments. In all RL simulations conducted in both the static and dynamic environments, we used a standard Q-learning formulation with identical hyperparameters unless otherwise specified. The hyperparameters used in all reinforcement-learning simulations are summarized in the table below.\nTable R1. RL parameter settings used in this work\nParameter Static environment Dynamic environment Description\nControls reward\nDiscount factor(γ) 0.9 0.9\nhorizon\nε-greedy\nExploration rate(ε) 0.1 0.1 exploration\npolicy\nStandard zero\nQ-table initialization Q(s,a) = 0 Q(s,a) = 0\ninitialization\nConsistent\nAction space 8 directions 8 directions\nacross tasks\nDetailed in\nLearning-rate Fixed per run (device- Stage-wise segmented\nSupplementary\nschedule derived) (device-derived)\nNotes 3-4\nThese settings were kept consistent across all experiments to ensure fair comparison among different learning-rate modulation strategies.\nResponse: We thank the reviewer for the valuable comments. In both the static and dynamic environments, we employed a sparse reward scheme with fixed magnitudes, determined directly by the task design. Specifically, a positive reward was assigned when the agent reached the goal, a penalty was applied upon colliding with an obstacle, and a small step cost was used for all regular transitions. The exact values used in the simulations are summarized in the table below.\nTable R2. Reward functions and magnitudes used in this work\nEvent Reward magnitude Description\nGoal reached +10 Terminal success\nCollision ‒1 Obstacle penalty\nNormal step ‒0.1 Small movement cost\nThis sparse reward formulation was applied consistently in both tasks to ensure comparability across learning-rate modulation strategies.\nResponse: We appreciate the reviewer for the valuable comments. The mapping between the reinforcement-learning (RL) episode index and the device-derived conductance states differs slightly between the static and dynamic environment settings, and we clarify the exact rules below.\n1) Static environment (single-task setting):\nIn the static environment experiments (Fig. 4 and Supplementary Fig. 35-37), the goal was to compare different fixed learning-rate values rather than to realize a time-varying schedule. In this case, a given device state ω (obtained from the intrinsic t gradient driven (IGD) and U-SVDP conductance modulation process) is first converted into a scalar learning rate α using the procedure described as\n(ω ω)\nα(t)= t+1 t ×100% (1)\nω\nt\nFor each trial, this α is then kept constant for all episodes in that training run. In other words, in the static environment, there is no episode-wise progression of α(t); instead, each conductance state ω corresponds to a single fixed learning rate that is k applied throughout the entire training process, and different runs use different α values k for comparison. Under static conditions, the IGD decay learning rate α reflects the Decay intrinsic monotonic attenuation of conductance change, which captures the slow, physically constrained relaxation dynamics of the device but suffers from a rapid decline in update magnitude at later states. In contrast, the U-SVDP-modulated learning rate α incorporates the same physical relaxation background while restoring U-SVDP conductance contrast through pulse-induced multilevel discretization. This enables the learning rate to remain within a more effective range across the entire training process.\nThe comparison therefore highlights the extent to which pulse-mediated recovery of conductance gradients contributes to stabilizing the learning-rate scale and preserving effective plasticity even when the underlying IGD relaxation becomes shallow.\nEmpirically, α consistently achieves faster convergence and lower variance U-SVDP\n22\nin path-length fluctuations compared with both α and the linear baseline. The Decay decay-only schedule α offers stable but overly conservative updates due to the Decay rapidly diminishing conductance differences in the late relaxation regime, resulting in slow convergence. The linear schedule (Supplementary Fig. 37), although straightforward to implement, exhibits sensitivity to initial parameter choices and often enters unstable update regions, leading to oscillatory convergence or premature stagnation. These results (Supplementary Table 3) demonstrated that the physically discretized, device-constrained learning-rate profile produced by U-SVDP strikes a more favorable balance between stability and adaptability, even in the absence of environmental dynamics.\nImportantly, this static-task comparison serves as a control experiment that isolates the contribution of the learning-rate schedule itself. The superior performance of α U- in such a simplified setting confirms that its advantage originates from the device- SVDP informed modulation mechanism, rather than from incidental factors such as task difficulty or exploration pressure. This provides a mechanistic basis for understanding the improvements subsequently observed in dynamic, multi-stage environments.\n2) Dynamic environment (multi-stage setting):\nIn the dynamic multi-stage environment task (Fig. 6 and Supplementary Fig. 38- 41), the environment undergoes systematic expansions from 50×50 to 75×75, 100×100, and finally 200×200. This progressive increase in state space and obstacle configuration introduces evolving exploration and planning demands, making the stability of the learning rate particularly critical. Unlike the static setting, where the learning rate remains fixed throughout each run, the dynamic environment requires a temporally evolving learning-rate schedule that can adapt to changes in task difficulty while maintaining stable convergence across transitions between stages.\nTo construct the α schedule, we first extracted 40 pseudo-nonvolatile U-SVDP conductance states from the IGD and U-SVDP scheme (one V followed by a series High of V ). These states were converted into 40 discrete learning-rate samples through the Low gradient-based mapping described in Supplementary Note 4. A smooth nonlinear\n23\nfunction α (p) was then obtained by fitting these samples over the normalized U-SVDP progress p∈[0,1].\nDirectly assigning a distinct learning rate to every episode would introduce excessive fluctuations in the temporal-difference updates and degrade convergence stability. To avoid this, the continuous nonlinear curve α (p) is discretized into a U-SVDP sequence of K monotonically evolving segments, each applied over a fixed (but unspecified) number of consecutive training episodes. Formally, let the training progress of episode 𝑒 be expressed as:\ne1\np = (2)\ne E 1\nmax\nwhere E is the total number of training episodes. The episode index 𝑒 is then mapped max to a segment index k∈{1, …, K} through a segmentation function:\nk =Φ(e) (3)\nwhere Φ(·) partitions the training horizon into K contiguous intervals of equal length.\nThe learning rate used at episode e is therefore given by:\nα(e)=α (p ) (4)\nU-SVDP k\nwhere p denotes the normalized progress associated with segment k.\nk\nFor fairness, the linear-learning-rate baseline is constructed using the same segmentation function Φ(·) and the same boundary values α , α , differing only in min max functional form:\nα (e)=α +(α α )p (5)\nLinear min max min k\nThe fixed-learning-rate baseline is expressed as:\na (e)=α (6)\nfixed const\nusing the identical update granularity.\nAcross all dynamic stages, the IGD-based nonlinear schedule provides faster convergence, reduced oscillation, and enhanced robustness. The steeper early-stage rise of α facilitates accelerated exploration, while its gradual saturation in later stages U-SVDP prevents overly aggressive updates and avoids the plateauing behavior frequently\n24\nobserved under linear scheduling. Fixed learning rates, lacking adaptive capability, perform poorly as environmental complexity increases. These results show that the device-derived nonlinear learning-rate trajectory acts as an intrinsic curriculum, enabling high plasticity in early exploration and increased stability during late-stage refinement across environments of progressively increasing complexity.\nWe thank the reviewer for pointing out the ambiguity in the original description of the device-to-learning-rate mapping. In response, we have revised Supplementary Note 4 to clarify the learning-rate construction procedure. Specifically, we now explicitly distinguish between (i) the continuous device-derived curve obtained from the IGD and U-SVDP measurements, and (ii) the segmented learning-rate schedule used during RL training. The specific modifications are as follows:\nSupplementary Note 4\nTo construct a physically grounded learning-rate schedule, we first utilized the recorded pseudo-nonvolatile conductance states obtained from the intrinsic gradient driven (IGD) and U-SVDP pulse scheme. Let {G} denote the sequence of conductance i values extracted from one high-voltage pulse followed by multiple low-voltage pulses.\nThe local conductance evolution was quantified using an absolute finite-difference operation:\nΔG = G G (7)\ni i+1 i\nThis generated a gradient profile that reflects the instantaneous adaptability of the device. The gradient sequence was then normalized into a bounded learning-rate range suitable for reinforcement learning updates:\nΔG min(ΔG)\nα =α + i ×（α -α ） (8)\ni min max(ΔG)min(ΔG) max min\nwhere α and α are predefined boundaries (e.g., 0.01 and 1.0). The set{α}, serves min max i as discrete, device-derived samples of learning plasticity. To obtain a smooth evolution profile applicable to long-horizon, multi-stage RL training, these samples were fitted to a continuous nonlinear function α (p) defined over the normalized training U-SVDP\n25\nprogress p∈[0,1] (as shown in Supplementary Fig.38).\nDirectly applying a distinct learning rate at every episode would introduce excessive fluctuations and destabilize Q-learning. Therefore, the fitted function α U- (p) was discretized into a sequence of gradually evolving learning-rate segments. SVDP Each segment was applied over a fixed (but unspecified) number of consecutive training episodes, ensuring smooth and stable learning-rate transitions while preserving the physical trend imposed by the device-derived curve.\nTo ensure fair comparison, the linear-learning-rate baseline was constructed using the same segmentation structure and the same boundary values α , α . Its min max learning-rate sequence is defined as:\nα (e)=α +(α α )p (9)\nLinear min max min k\nwhere p corresponds to the normalized progress associated with the k-th segment. The k fixed-learning-rate baseline was implemented using the same update granularity but with a constant learning rate:\na (e)=α (10)\nfixed const\nIn summary, these formulations ensure that all learning-rate strategies differ only in functional form, while maintaining identical update boundaries and identical segmentation granularity. This allows the performance difference observed in the dynamic multi-stage task to be attributed solely to the physical, device-informed nonlinear learning-rate evolution.\nResponse: We thank the reviewer for the valuable comments. While a full HIL demonstration is beyond the scope of the present study, the proposed approach is, in principle, well suited for real-time HIL integration.\nIn this work, the memristive device does not function as a static weight storage element, but instead provides a dynamic, continuously accessible learning-rate signal derived from intrinsic gradient-driven transient and relaxation behavior. Such a signal can be readily interfaced with real-time control systems through standard analog-to-digital conversion, making it compatible with online parameter modulation in reinforcement learning frameworks.\nMoreover, the demonstrated stability and repeatability of the transient dynamics under repeated electrical cycling and thermal perturbation are essential prerequisites for HIL operation, where consistent real-time feedback is required. The characteristic timescales of the device dynamics are also compatible with learning-rate updates in typical control and reinforcement learning loops. A full real-time HIL implementation, including closed-loop integration with embedded controllers and physical environments, is planned as part of future work. This will enable direct evaluation of the system-level performance and latency considerations in practical robotic or control applications.\nResponse: We appreciate the reviewer for the valuable comments to improve our work.\nFirstly, regarding the effect demonstrated in Fig. 6, increasing the map size from 50×50 to 200×200 produced a clear escalation in task difficulty. The enlarged spatial domain expanded the state space quadratically and lengthened the minimal feasible trajectories, thereby requiring the agent to operate over substantially longer exploration horizons and deeper value-propagation dependencies. The results in Fig. 6 confirmed that the IGD+U-SVDP-derived adaptive learning rate (α ) enabled consistently faster U-SVDP convergence as the navigation demands increased. However, although map-size scaling captured a form of increasing difficulty, it introduced only a limited degree of environmental non-stationarity, as the underlying reward structure and obstacle topology remained relatively smooth across scales.\nTo address the reviewer’s concerns regarding both task complexity and the sufficiency of non-stationary conditions, we therefore extended our evaluation beyond map-size scaling and systematically examined agent behavior under varying obstacle densities and reward structures.\n1) Evaluation under varying obstacle densities:\nIn the original configuration, a sparse reward paradigm was used, where positive reinforcement was provided only upon goal achievement. While suitable for small mazes, this setting becomes increasingly restrictive as environmental complexity grows. To assess the adaptability of the α under controlled complexity escalation, we U-SVDP therefore fixed the map size at 50×50 and progressively increased the obstacle density to 0.15, 0.25, 0.35, and 0.45. To ensure methodological consistency, both the α U-SVDP schedule and the linear learning-rate schedule (α ) followed the same episode-wise Linear update rules as those used in Fig. 6 and Supplementary Note 4. The discount factor (γ) and exploration rate (ε) were kept constant at 0.9 and 0.1, respectively, and the maximum episode length was set to 50×50×2 iterations to ensure a fair comparison across conditions.\nSupplementary Fig. 42. The optimal paths and convergence iterations of the α and U-SVDP α modulation schedules when the obstacle density of the map is 0.15\nLinear\nAs the obstacle density increased, the advantages of the α became more U-SVDP pronounced. As illustrated in Supplementary Figs. 43-44, at densities of 0.25 and 0.35, the α converged in only 1,236 and 592 iterations, respectively, while achieving U-SVDP shorter paths (50 and 49 steps). In contrast, the α agent generated trajectories with Linear considerably longer detour distances (3,694 and 2,300 iterations, respectively) and more curvature (51 and 54 steps), reflecting poorer adaptability to increased structural complexity.\nSupplementary Fig. 43. The optimal paths and convergence iterations of the α and α\nU-SVDP Linear\n29\nmodulation schedules when the obstacle density of the map is 0.25.\nSupplementary Fig. 44. The optimal paths and convergence iterations of the α and\nU-SVDP\nα modulation schedules when the obstacle density of the map is 0.35.\nLinear\nUnder the highest obstacle density of 0.45 (Supplementary Fig. 45), training with α remained effective, achieving convergence within ~1,000 iterations. In contrast, U-SVDP the α agent again required more iterations (6,507) to reach an equivalent optimal Linear path length (54 steps). As summarized in Supplementary Fig. 46, at all tested density levels, α consistently demonstrated faster convergence speed, higher stability, U-SVDP and stronger adaptability to environmental complexity.\nSupplementary Fig. 45. The optimal paths and convergence iterations of the α and\nU-SVDP\nα modulation schedules when the obstacle density of the map is 0.45.\nLinear\nα and α modulation strategies under different obstacle density of the map.\nU-SVDP Linear\n2) Evaluation under varying reward structures:\nTo comprehensively assess the robustness and generality of the proposed adaptive learning-rate framework, we further examined the performance of α under U-SVDP different reward structures, following its validated superiority across varying map sizes and obstacle densities. In the manuscript experiments, a sparse reward (SR) scheme was employed. Although conceptually simple, sparse reward structures impose well-known credit-assignment challenges that tend to slow learning in complex environments. In our experimental setting, this property is not a limitation; rather, it makes sparse rewards a rigorous testing platform for evaluating whether learning-rate scheduling strategies can maintain stability and adaptability in the presence of minimal feedback — which is precisely why we retained this reward setting in the final dynamic-environment task.\nTo determine whether the advantages of α depended on the reward structures, U-SVDP we introduced several alternative reward mechanisms and systematically compared the behavior of agents trained with α and with a conventional α schedule under U-SVDP Linear each setting.\n2.1) Sparse reward structure\nUnder the SR setting, where feedback was provided only through terminal\n31\noutcomes and a small step penalty, the agent trained with α exhibited markedly U-SVDP faster and more stable convergence. Across trials, α reached the optimal trajectory U-SVDP within 516 iterations and achieved a 50-step path, whereas the α schedule required Linear 3,589 iterations and converged to a longer 53-step solution (Supplementary Fig. 47).\nThe α agent also showed higher sensitivity to delayed rewards, often displaying Linear slower value propagation and increased variability during training. These results indicated that the advantages of α persisted even when intermediate feedback was U-SVDP extremely limited, confirming its robustness under SR conditions.\nSupplementary Fig. 47. Performance comparison between α and α modulation\nU-SVDP Linear\nschedules under the SR mechanism, with a map size of 50×50 and an obstacle density of 0.35.\n2.2) Distance-difference-based dense reward structure\nTo address the directional deficiency inherent in the SR structure, the distance-difference-based dense reward (DDR) structure was further introduced. The essence of this mechanism lies in transforming the spatial distance between the agent and the target into an immediate reward feedback signal, thereby providing explicit directional guidance for each decision step.\nTo verify whether the α U-SVDP agent maintains its superior performance under the DDR mechanism, a systematic comparison was conducted between the α and U-SVDP traditional α modulation schedules. The start and goal coordinates were set to (5,5) Linear and (45,45), respectively. The γ and ε were fixed at 0.9 and 0.1, respectively. As shown in Supplementary Fig. 48, α agent converged to a shorter final path (65 U-SVDP steps) while requiring markedly fewer iterations than α agent (68 steps), corresponding Linear\n32\nto improvements of approximately 4.6% in path optimality and 24.09% in convergence efficiency.\nSupplementary Fig. 48. Performance comparison between α and α modulation\nU-SVDP Linear\nschedules under the DDR mechanism, with a map size of 50×50 and an obstacle density of 0.35.\nAlthough the α agent exhibited small stepwise refinements near convergence, U-SVDP these reflected local improvements in the final trajectory rather than oscillatory instability. Unlike α agent, which settled prematurely into a stable yet suboptimal Linear plateau, α agent continued to make consistent incremental progress and ultimately U-SVDP achieved a shorter final trajectory. These results clearly demonstrate the advantages of the α modulation schedule under DDR conditions, achieving both faster U-SVDP convergence and higher final path quality while maintaining enhanced training stability.\n2.3) Multi-component dense reward structure\nTo overcome the discretization limitation inherent in the DDR-based reward introduced in the previous section, a multi-component dense reward (MCDR) structure was applied. Compared with the DDR, the key improvement of this method lies in redefining the reward as a continuous function of distance improvement, enabling the reward signal to reflect the actual magnitude of each action’s effectiveness more accurately. Because dense and continuous rewards tend to minimize performance differences in small grids, a larger 100×100 environment was selected to ensure that the MCDR evaluation remained sensitive and discriminative for comparing the two learning-rate schedules.\nSupplementary Fig. 49 illustrates the convergence behavior of the two learning-rate schedules under the MCDR configuration. In the 100×100 environment, the agent\n33\ntrained with α converged more rapidly and achieved a shorter steady-state path U-SVDP length (135 steps) than the agent using αLinear (138 steps), approaching the theoretical limit for this setting. Notably, α maintained this advantage even when the state U-SVDP space increased fourfold and the required trajectory length grew by a factor of 2.5, indicating stable performance in large, highly constrained environments.\nSupplementary Fig. 49. Performance comparison between α and α modulation\nU-SVDP Linear\nschedules under the MCDR, with a map size of 100×100 and an obstacle density of 0.40.\nThe α agent reached the target earlier due to stochastic early exploration, Linear whereas the α agent exhibited a delayed first arrival because its higher initial U-SVDP learning rate promoted broader and more systematic exploration. This early-phase delay proved beneficial in the subsequent training stage, as the more comprehensive state-space coverage enabled α to converge more rapidly and to a shorter steady-state U-SVDP path, while the linear schedule improved only marginally after its initial success.\nFrom a physical perspective, this early-stage behavior originates from the device dynamics. The built-in concentration gradient provides a strong initial driving force for ionic redistribution, producing a relatively large conductance-change rate (dG/dt) at the beginning of training. As the system undergoes slow interfacial-barrier recovery, this intrinsic driving force becomes progressively suppressed, naturally reducing the rate of conductance evolution. Importantly, within this gradually decelerated dynamical regime, the U-SVDP process continues to modulate vacancy drift through its inherent asymmetry and slow relaxation, allowing the device to sustain an effective learning-rate adaptation even when the underlying ionic motion becomes slow.\n2.4) Normalized progress reward\nThe strong normalization inherent to NPR compresses the temporal-difference signal into a narrow range, producing a weak-modulation regime that reduces sensitivity to the learning-rate schedule. Unlike SR, DDR, and MCDR, where reward gradients provide substantial modulation pressure, the NPR structure smooths the learning signal and naturally suppresses differences between α and α . The U-SVDP Linear objective of this experiment is therefore not to amplify performance gaps but to assess the stability of α under conditions where reward gradients provide minimal U-SVDP guidance.\nSupplementary Fig. 50 showed the convergence behavior on a 100×100 map with a low obstacle density of 0.1. Both schedules displayed similarly rapid improvement at the beginning because the strong normalization of the NPR structure compressed the temporal-difference signal and suppressed early differences between learning-rate schedules. As training progressed, the α agent entered a prolonged plateau around Linear 173 steps, whereas the α agent continued to refine the policy and stabilized at 170 U-SVDP steps. A further refinement to 165 steps was achieved under α , although this U-SVDP required additional iterations. This extended optimization reflected continued value-function improvement enabled by adaptive learning-rate modulation rather than instability. These observations confirmed that α maintained stable and consistent U-SVDP behavior even in weak-modulation regimes, while more substantial advantages appeared in settings where reward gradients were stronger.\nU-SVDP Linear\nschedules under the NPR mechanism, with a map size of 100×100 and an obstacle density of 0.10.\nIn summary, the α -based agent consistently exhibits superior convergence U-SVDP performance under different reward mechanisms. All comparative results and detailed design settings are illustrated in Supplementary Fig. 51 and Table R3.\nSupplementary Fig. 51. The comparison of optimal paths and convergence iterations of the\nα and α modulation schedules under different reward structure.\nU-SVDP Linear\nMap size/ Convergence Design Physical\nReward structure Reward design Optimal path\nObstacle density iteration motivation perspective\nThe slow relaxation\nof intrinsic oxygen-\nTo evaluate ionic distributions\nlearning-rate preserves learning\nSR Collision: ‒1; modulation under responsiveness\n50×50\n(γ = 0.9 Goal: +10; 50 53 516 3589 minimal guidance, under infrequent\n0.40\nε=0.1) Others: ‒0.1 where credit feedback, enabling\nassignment is α to enhance\nU-SVDP\nmost challenging. long-horizon value\npropagation when\nrewards are sparse.\nThe gradual PNV-\ndriven conductance\nTo provide modulation filters\nexplicit rapid temporal-\ndirectional difference\nDDR Collision: ‒5;\n50×50 feedback by fluctuations,\n(γ = 0.9 Move closer: 65 68 6250 7756\n0.35 shaping rewards producing smoother\nε=0.1) +0.1;\naccording to and more stable\nstepwise distance update dynamics\nimprovement. under stepwise\ndistance-based\nrewards.\ndiffusion-drift\ncoupling within the\nTo strengthen\ndevice generates\nexploration\nCollision: ‒5;\nnatural gain\nefficiency by\nMCDR Goal: +100;\ncompression,\n100×100 integrating\n(γ = 0.9 Move closer: 135 138 7389 7734\npreventing\nmultiple shaping\nε=0.1) +2Δd;\novershooting when\nterms that more\nStep penalty:\ndense and strongly\nstrongly encode\n‒0.1\nshaped reward\ntask structure.\ncomponents\namplify update\nmagnitudes.\nThe intrinsic\noxygen-gradient\nevolution provides a\nTo achieve scale- slowly varying\ninvariant feedback internal potential\nNPR\nCollision: ‒10;\nand eliminate landscape, acting as\n(γ = 0.9 165 169 6940 5671\nGoal: +50; 0.10 reward imbalance a physical stabilizer\nε=0.1)\nacross maps of that sustains scale-\ndifferent sizes. invariant learning\nefficiency when\nreward gradients are\nheavily normalized.\n38\n2.5) Dynamic Non-Stationary Task\nTo rigorously examine the adaptability of the proposed α learning-rate\nU-SVDP\nmodulation under genuinely non-stationary conditions, we designed a continuous dynamic-goal navigation task in a 100×100 environment. Unlike the earlier dynamic tasks that only enlarged the map size while keeping the underlying value landscape relatively uniform (Figure 6), this dynamic-goal configuration establishes a genuinely non-stationary learning process in which: (i) all goals share the same Q-table without reinitialization, causing value estimates to accumulate across stages and making later decisions inherently history-dependent; (ii) the target position shifts immediately after each success, producing a continually changing value landscape that forces the agent to readjust its policy after every transition; and (iii) the obstacle density is progressively increased during training (up to 0.5), so that the geometric constraints tighten at the same time as earlier Q-value biases propagate, substantially elevating the difficulty of mid-to-late stages. To prevent reward shaping from overshadowing the intrinsic learning-rate dynamics, the SR structure was adopted throughout (γ = 0.9, ε = 0.1), ensuring that adaptation is primarily driven by the modulation mechanism itself. Each of the ten goals was allocated 400 episodes, forming a continuous sequence of 4,000 episodes with no reset of either the Q-table or the exploration history.\nIn the first three stages, both α and α failed to reach the target, a natural\nLinear U-SVDP\nconsequence of beginning with an uninitialized Q-table under a SR scheme that provides almost no intermediate directional signal (Supplementary Fig. 52a-b). As the shared value function gradually accumulated structural information, the behavioral differences between the two schedules became increasingly clear. Across the ten sequential goals, the α -modulated agent exhibited smoother episode-wise\nU-SVDP\nrecovery after each goal relocation, a progressive reduction in converged path length, and a consistently declining convergence envelope even as obstacle density rose (especially the goal 4 to 7, Supplementary Figs. 52c-d and Fig. 53a-d). In contrast, the α schedule produced strong oscillations, repeated divergence-reconvergence cycles,\nLinear\nand slower re-stabilization in the mid-to-late stages where non-stationarity and environmental constraints were most severe. The accompanying bar plots further reveal\n39\nthat α achieved shorter converged paths and fewer effective iterations in most\nU-SVDP\nstages, indicating that its adaptive learning-rate dynamics were better suited to the continual reorganization required by this history-dependent, dynamically evolving environment.\nSupplementary Fig. 52. Comparison of optimal paths and convergence performance between\nα -based and α -based modulation in dynamic environments. The convergence curve\nU-SVDP Linear\nof a α -based agent and b α -based agent. c Optimal paths obtained in each stage. d\nU-SVDP Linear\nConvergence performance across stages.\nU-SVDP\ndynamic-goal tasks in a 100×100 environment.\nOverall, the continuous dynamic-goal experiments demonstrate that α\nU-SVDP\nprovides a fundamentally more stable and adaptive learning regime for long-horizon, history-dependent RL. Mechanistically, the modulation advantage originates from the coupled IGD+U-SVDP ionic dynamics. The intrinsic oxygen-concentration gradient provides a directionally biased yet kinetically constrained driving force for conductance evolution, enabling newly acquired value information to be integrated without the abrupt redistribution seen in homogeneous systems. As training progresses, the subsequent U-SVDP-governed slow relaxation further attenuates the update magnitude, forming a naturally moderated learning-rate trajectory that prevents instability under the accumulated non-stationarity across stages. This synergy yields a self-regulated learning-rate trajectory that is naturally matched to the temporal structure of continual RL—fast enough to adapt immediately after each environmental shift, yet sufficiently tempered to avoid catastrophic overwriting of accumulated knowledge.\nIn response to the reviewer’s suggestions, we have conducted additional analyses covering obstacle-density escalation, alternative reward structures, and\n41\na more challenging continual dynamic-goal task. These results have been incorporated into both the revised manuscript and the Supplementary information, providing a more comprehensive validation of the proposed α modulation under\nU-SVDP\ndiverse and increasingly non-stationary conditions.\nOn page 12-13 of the manuscript, “Compared with traditional method (linear modulation), the α (p) reduced the total training episodes of the agent by 35.65%.\nU-SVDP\nThe optimal path results were compared across map sizes in detail, demonstrating the effectiveness of this device-derived learning-rate modulation (Fig. 6e-h). Beyond map-size scaling, additional evaluations were performed under progressively increasing obstacle densities, multiple reward structures, and a fully non-stationary dynamic-goal environment, all of which consistently confirmed the robustness and adaptability of the α modulation (Supplementary Figs. 42-53). These results underscore that the\nU-SVDP\nnonlinear learning-rate trajectory derived from intrinsic-gradient-mediated and slowly relaxing ionic dynamics provides an effective physical substrate for dynamic learning-rate modulation, enabling efficient stage-wise adaptation in continual RL through the close coupling between device-level processes and algorithmic behavior.”\nOn page 55-64 of Supplementary Information:\nSupplementary Fig. 42. The optimal paths and convergence iterations of the α and α\nU-SVDP Linear\nmodulation schedules when the obstacle density of the map is 0.15.\nU-SVDP Linear\nmodulation schedules when the obstacle density of the map is 0.25.\nSupplementary Fig. 44. The optimal paths and convergence iterations of the α and α\nU-SVDP Linear\nmodulation schedules when the obstacle density of the map is 0.35.\nSupplementary Fig. 45. The optimal paths and convergence iterations of the α and α\nU-SVDP Linear\nmodulation schedules when the obstacle density of the map is 0.45.\nα and α modulation strategies under different obstacle density of the map.\nU-SVDP Linear\nEvaluation under varying obstacle densities\nTo assess the adaptability of the α under controlled complexity escalation, we\nU-SVDP\ntherefore fixed the map size at 50×50 and progressively increased the obstacle density to 0.15, 0.25, 0.35, and 0.45. The γ and ε were kept constant at 0.9 and 0.1, respectively, and the maximum episode length was set to 50×50×2 iterations to ensure a fair comparison across conditions.\nAt a low obstacle density of 0.15 (Supplementary Fig. 42), the agent under α\nU-SVDP\nrapidly identified smooth and efficient trajectories, achieving convergence with substantially fewer iterations, and reduced computational cost. Although the linear modulation eventually converged to a slightly shorter path (47 steps), it required substantially more iterations (7,194) to reach this solution, approximately 4.1 times that of the α case, highlighting its inefficiency in low-complexity environments.\nU-SVDP\nAs the obstacle density increased, the advantages of the α became more\nU-SVDP\npronounced. As illustrated in Supplementary Fig. 43-44, at densities of 0.25 and 0.35, the α converged in only 1,236 and 592 iterations, respectively, while achieving\nU-SVDP\nshorter paths (50 and 49 steps). In contrast, the α agent generated trajectories with\nLinear\nconsiderably longer detour distances (3,694 and 2,300 iterations, respectively) and more curvature (51 and 54 steps), reflecting reduced adaptability to increasing structural complexity.\nα remained effective, achieving convergence within ~1,000 iterations. In contrast,\nU-SVDP\nthe α agent again required more iterations (6,507) to reach an equivalent optimal\nLinear\npath length (54 steps). As summarized in Supplementary Fig. 46, at all tested density levels, α consistently demonstrated faster convergence speed, higher stability,\nU-SVDP\nand stronger adaptability to environmental complexity.\nSupplementary Fig. 47. Performance comparison between α and α modulation\nU-SVDP Linear\nschedules under the SR mechanism, with a map size of 50 × 50 and an obstacle density of 0.35.\nSupplementary Fig. 48. Performance comparison between α and α modulation\nU-SVDP Linear\nschedules under the DDR mechanism, with a map size of 50×50 and an obstacle density of 0.35.\nU-SVDP Linear\nschedules under the MCDR mechanism, with a map size of 100×100 and an obstacle density of 0.40.\nSupplementary Fig. 50. Performance comparison between α and α modulation\nU-SVDP Linear\nschedules under the NPR mechanism, with a map size of 100×100 and an obstacle density of 0.10.\nSupplementary Fig. 51. The comparison of optimal paths and convergence iterations of the\nα and α modulation schedules under different reward structure.\nU-SVDP Linear\nEvaluation under varying reward structures\nTo examine the generality of the α learning-rate modulation beyond\nU-SVDP\nvariations in map size and obstacle density, the performance under four representative reward structures were further evaluated, including sparse reward (SR), distance-difference reward (DDR), multi-component dense reward (MCDR), and normalized-progress reward (NPR). These settings progressively enhance directional guidance and\n46\nreward continuity, enabling a systematic assessment of how αU-SVDP interacts with different forms of temporal-difference signals.\nSparse reward (SR)\nUnder the SR setting, where feedback was provided only through terminal outcomes and a small step penalty, the agent trained with α exhibited markedly\nU-SVDP\nfaster and more stable convergence. Across trials, α reached the optimal trajectory\nU-SVDP\nwithin 516 iterations and achieved a 50-step path, whereas the α schedule required\nLinear\n3,589 iterations and converged to a longer 53-step solution (Supplementary Fig. 47).\nThe α agent also showed higher sensitivity to delayed rewards, often displaying\nLinear\nslower value propagation and increased variability during training. These results indicated that the advantages of α persisted even when intermediate feedback was\nU-SVDP\nextremely limited, confirming its robustness under SR conditions.\nDistance-difference reward (DDR)\nDDR introduces stepwise directional feedback by rewarding distance reductions.\nAs shown in Supplementary Fig. 48, α agent converged to a shorter final path (65\nU-SVDP\nsteps) while requiring markedly fewer iterations than α agent (68 steps),\nLinear\ncorresponding to improvements of approximately 4.6% in path optimality and 24.09% in convergence efficiency. The small refinements near convergence reflected incremental path improvement rather than instability, demonstrating that α\nU-SVDP\nremains effective when denser and more informative reward gradients are available.\nMulti-component dense reward (MCDR)\nTo capture richer task structure, MCDR incorporates continuous distance gains, goal bonuses, collision penalties, and step costs. Because dense shaping often minimizes performance differences in small grids, a larger 100 × 100 environment was adopted. As shown in Supplementary Fig. 49, α maintained faster convergence\nU-SVDP\nand reached a shorter steady-state path (135 steps) compared to αLinear (138 steps).\nThe early-stage delay in first-goal arrival under α resulted from its broader initial\nU-SVDP\nexploration, which led to more comprehensive value propagation and faster downstream optimization in deeper stages.\nGiven that MCDR’s shaping terms are scale-dependent, we introduced NPR to achieve scale-invariant reward magnitudes. NPR compresses temporal-difference signals, forming a weak-modulation regime where differences between learning-rate schedules naturally diminish. As shown in Supplementary Fig. 50, both schedules improved rapidly at the beginning; however, αLinear entered a persistent plateau at 173 steps, whereas α continued refining the policy to 170\nU-SVDP\nsteps, and further to 165\nsteps after extended iterations. Although the performance gap is intentionally small under NPR, α nonetheless demonstrated continued refinement and stable\nU-SVDP\nbehavior.\nAcross SR, DDR, MCDR, and NPR (Supplementary Fig. 51 and Supplementary Table 4), α consistently produced equal or better steady-state paths and more\nU-SVDP\nfavorable convergence characteristics. Dense reward settings emphasized its efficiency; sparse or normalized settings highlighted its stability. These results collectively show that α adapts robustly across a broad spectrum of reward-shaping regimes, and its\nU-SVDP\nadvantages intensify when reward gradients play an active role in value propagation.\nreward structures\nReward Reward\nDesign motivation Physical perspective\nstructure design\nConstrained IGD drift\nEvaluate learning dynamics\nSR Collision: ‒1;\nsuppresses abrupt updates; U-\nunder minimal feedback and\n(γ = 0.9 Goal: +10;\nSVDP reconstruction stabilizes\nlong credit-assignment\nε=0.1) Others: ‒0.1\nvalue propagation under\nchains.\nsparse/delayed rewards.\nPNV-driven conductance\nDDR Collision: ‒5;\nProvide stepwise directional smoothing prevents oscillatory\n(γ = 0.9 Move closer:\nguidance without heavy temporal difference updates,\nε=0.1) +0.1;\nshaping. aiding stable directional\nstepwise distance\nrefinement.\nimprovement.\nCollision: ‒5;\nDiffusion-governed\nGoal: +100;\nMCDR Capture richer task structure conductance dynamics impose\nMove closer:\n(γ = 0.9 through continuous shaping natural gain compression that\n+2Δd;\nt\nε=0.1) terms. prevents overreaction to dense\nStep penalty:\ngradients.\n‒0.1\nIGD gradient evolution acts as\nNPR Collision: ‒ Achieve scale-invariant internal normalization,\n(γ = 0.9 10;\nshaping independent of map preserving adaptive behavior\nε=0.1) Goal: +50;\nsize. even when reward magnitudes\nare compressed.\nα -based and α -based modulation in dynamic environments. The convergence curve\nU-SVDP Linear\nof a α -based agent and b α -based agent. c Optimal paths obtained in each stage. d\nU-SVDP Linear\nConvergence performance across stages.\nSupplementary Fig. 53. Optimal path results with α -modulated agents across sequential\nU-SVDP\ndynamic-goal tasks in a 100×100 environment.\n50\n“To rigorously examine the adaptability of the proposed α learning-rate\nU-SVDP\nmodulation under genuinely non-stationary conditions, we designed a continuous dynamic-goal navigation task in a 100×100 environment. Unlike the earlier dynamic tasks that only enlarged the map size while keeping the underlying value landscape relatively uniform (Figure 6), this dynamic-goal configuration establishes a genuinely non-stationary learning process in which: (i) all goals share the same Q-table without reinitialization, causing value estimates to accumulate across stages and making later decisions inherently history-dependent; (ii) the target position shifts immediately after each success, producing a continually changing value landscape that forces the agent to readjust its policy after every transition; and (iii) the obstacle density is progressively increased during training (up to 0.5), so that the geometric constraints tighten at the same time as earlier Q-value biases propagate, substantially elevating the difficulty of mid-to-late stages. To prevent reward shaping from overshadowing the intrinsic learning-rate dynamics, the SR structure was adopted throughout (γ = 0.9, ε = 0.1), ensuring that adaptation is primarily driven by the modulation mechanism itself. Each of the ten goals was allocated 400 episodes, forming a continuous sequence of 4,000 episodes with no reset of either the Q-table or the exploration history.\nIn the first three stages, both α and α failed to reach the target, a natural\nLinear U-SVDP\nconsequence of beginning with an uninitialized Q-table under a SR scheme that provides almost no intermediate directional signal (Supplementary Fig. 52a-b). As the shared value function gradually accumulated structural information, the behavioral differences between the two schedules became increasingly clear. Across the ten sequential goals, the α -modulated agent exhibited smoother episode-wise\nU-SVDP\nrecovery after each goal relocation, a progressive reduction in converged path length, and a consistently declining convergence envelope even as obstacle density rose (especially the goal 4 to 7, Supplementary Fig. 52c-d and Supplementary Fig. 53a-f). In contrast, the α schedule produced strong oscillations, repeated divergence-\nLinear\nreconvergence cycles, and slower re-stabilization in the mid-to-late stages where non-stationarity and environmental constraints were most severe.”\nThe manuscript entitled “Intrinsic Gradient Oxygen-Driven Second-Order Memristors for Continual Reinforcement Learning” proposed a second order memristor with long retention time and its application in continual reinforcement learning. Generally, the manuscript would be suitable for publishing on Nature Communication after minor revision:\nResponse:\nWe sincerely thank the reviewer for the positive and encouraging evaluation of our work. We are pleased that the reviewer considers the manuscript potentially suitable for publication after minor revision. In response, we have carefully addressed all specific comments by providing additional simulation results, detailed discussions on the underlying issues, and clearer figures in the manuscript and Supplementary Information. All revisions have been clearly marked in the revised manuscript, and we hope that these revisions satisfactorily address the reviewer’s concerns.\nResponse: We thank the reviewer for the valuable comments. The stability of the proposed second-order memristive system has been evaluated from multiple complementary aspects, including long-term storage stability, device-to-device reproducibility, and robustness under repeated electrical cycling.\nStability of long-term storage. As shown in Fig. R1,we conducted long-term storage stability tests for 1, 7 and 21 days. The results showed that the I-V characteristic of the device remained almost unchanged, indicating that its underlying transport mechanism did not change. Although the absolute current level exhibited a moderate drift, this behavior was attributed to mild self-oxidation of the AlO /ZnTPP stack, x which introduced additional oxygen ions and slightly enhanced the electronic conductivity without altering the conduction pathway. The strong consistency in curve shape confirms the robustness of the intrinsic-gradient-regulated multi-filament network and validates the long-term stability of the device.\nFig. R1. I-V characteristics measured after different storage times. a-c Typical I-V characteristics under 0→9→0 V consecutive scanning loops with storage for 1, 7, and 21 days. d-f Typical I-V characteristics under 0→6→0 V consecutive scanning loops with storage for 1, 7, and 21 days.\nSupplementary Fig. 19. Statistical analysis of second-order relaxation dynamics across 32 devices. a Mean relaxation curve averaged over 32 devices after programming with a 9 V, 100 ms voltage pulse, measured at a read voltage of 1 V. b Statistical distribution of the extracted fast (τ ) 1 and slow (τ ) bi-exponential relaxation time constants across 32 devices, presented as box plots with 2 individual data points overlaid. c Statistical distribution of the ratio τ /τ .\n2 1\ndemonstrates that, once a stable intrinsic concentration-gradient-governed relaxation is established, a properly designed U-SVDP scheme enables reproducible modulation of PNV states across different devices. Although the absolute conductance levels of the PNV states exhibit device-to-device variations, the temporal evolution and modulation trends remain highly consistent among devices. Importantly, each conductance change (ΔG ) is defined relative to the intrinsic baseline of the individual device rather U-SVDP than to an absolute conductance value. As a result, device-specific offsets in the PNV level do not affect the effective modulation behavior.\nThis trend-level consistency ensures that the U-SVDP-induced ΔG can be U-SVDP reliably used as a device-independent control variable, which is particularly critical for learning-rate mapping in neuromorphic or reinforcement-learning-inspired applications. In this context, the relative evolution of conductance, rather than its absolute magnitude, governs the learning dynamics, making the observed device-to-device variations in PNV states non-detrimental to system-level functionality.\nSupplementary Fig. 28. Statistical evaluation of U-SVDP modulation across 32 devices. a Transient conductance measured under the 100 ms V (9 V) pulse. b Transient conductance High\n55\nmeasured under the consecutive100 ms V pulse (6 V, Δt= 5 s). c Evolution of the conductance Low states reconstructed by the V pulses (read at 1 V), corresponding to the reconfiguration of PNV Low states.\nStability under conventional bipolar cycling. We first evaluate the device under conventional bipolar voltage cycling using positive and negative pulses, which is a commonly employed reliability assessment for memristive devices. Unlike nonvolatile memory applications, this test does not aim to evaluate the endurance of precisely programmed conductance states. Instead, since the device operates as a second-order memristor, the focus is on the reproducibility of the transient conductance response under repeated electrical excitation. As shown in Supplementary Fig. 30, over 104 cycles of bipolar excitation, the transient conductance response exhibits no cumulative degradation, suggesting that the device can tolerate repeated electrical perturbations while maintaining stable dynamical behavior. This measurement therefore serves as a baseline verification of device robustness, rather than a metric of memory endurance.\nSupplementary Fig. 30. Endurance characteristics evaluated using a conventional bipolar switching scheme. a Endurance performance measured under a conventional bipolar voltage switching scheme over 10,000 consecutive cycles (V = 9 V, t = 50 μs; V = ‒9 V, t = 80 positive d negative d μs, Δt = 20 μs), the interval between consecutive cycles was fixed at 150 μs. b Enlarged view of several representative cycles, illustrating the detailed conductance response to the bipolar programming pulses and the corresponding read operation (1 V).\nEvaluation of second-order memristor stability using U-SVDP cycling scheme.\nMore importantly, the endurance requirement in this work arises from the fundamental operating characteristics of a second-order memristor. Unlike conventional nonvolatile\n56\nmemristive devices, a second-order memristor does not possess well-defined static memory states that can be repeatedly programmed and retained. Instead, its functionality is encoded in transient and relaxation dynamics governed by ionic processes. Consequently, conventional endurance tests that rely on repeated programming of stable conductance states are not directly applicable. Endurance must therefore be evaluated in a manner consistent with the intrinsic operating mechanism of the device, namely the stability of the intrinsic oxygen concentration gradient that underpins the dynamic response. This consideration, we designed a specialized cycling scheme to investigate whether the intrinsic gradient remains effective and directionally stable under repeated electrical disturbances, rather than evaluating the endurance of static conductance states.\nTo directly probe this aspect, we employed a U-SVDP cycling scheme.\n(Supplementary Fig. 31). In this scheme, a V pulse was used to activate oxygen ions High and drive the system away from equilibrium, thereby generating a population of ions that could participate in subsequent diffusion processes. This pulse was not intended to program a stable conductance state, but rather to provide a controlled perturbation to the ionic distribution. Under short-pulse cycling (t =50 μs), the oxygen ions activated d by V rapidly redistributed and relaxed, such that the subsequent V pulse could High Low not effectively reutilize these transiently activated ions to reconstruct the conductance state. Nevertheless, owing to the presence of the intrinsic concentration gradient, the transient response induced by V still exhibited a modulation polarity opposite to that Low induced by V .\nHigh\nIn contrast, when the pulse width of V was increased to 1 ms, a substantially High larger population of diffusible oxygen ions remained available for redistribution (Supplementary Fig. 32). In this case, the V pulse was able to actively engage these Low ions under the constraint of the intrinsic concentration gradient, leading to a more pronounced and reproducible conductance modulation with negative polarity under readout conditions. These behaviors indicate that, while the intrinsic concentration gradient dictates the directionality of the response, the availability of diffusible ions\n57\ndetermines the magnitude and persistence of the V -induced conductance modulation.\nLow\nOverall, these observations confirm that the U-SVDP cycling scheme provides an effective means to evaluate the stability of the intrinsic concentration gradient, which constitutes the essential endurance metric for this second-order memristive system.\nSupplementary Fig. 31. Endurance performance under the U-SVDP modulation scheme. a Endurance performance evaluated using the U-SVDP modulation scheme (V = 9 V V = 6 V, High Low t = 50 μs, Δt = 100 μs), the interval between consecutive cycles was fixed at 950 μs. b Enlarged d view of several representative cycles, illustrating the detailed conductance response to the U-SVDP programming pulses and the corresponding read operation (1 V).\nSupplementary Fig. 32. Endurance performance under the U-SVDP modulation scheme. a Endurance performance evaluated using the U-SVDP modulation scheme (V = 9 V V = 6 V, High Low t = 1 ms, Δt = 10 ms), the interval between consecutive cycles was fixed at 1 ms. b Enlarged view d of several representative cycles.\nIn response to the reviewer’s concern about the stability/durability of the device. We further added device-to-device analyses to demonstrate the storage stability, decay stability, and endurance stability. These revisions are reflected in the manuscript and\nOn page 7 of the manuscript, “This process demonstrated the second-order characteristic.[27] The reproducibility and pulse-width-dependent tunability of this second-order memristive behavior are confirmed in Supplementary Figs. 19-20.”\nOn page 11 of the manuscript, “As a result, the conductance change (ΔG =−72.4%) after 16 successive PNV states was significantly greater than the −51.7% U- SVDP spontaneous decay (ΔG) (Fig. 4j). Furthermore, the stability of the intrinsic concentration-gradient-driven U-SVDP modulation PNV states were demonstrated in Supplementary Fig. 28.”\nOn page 11 of the manuscript, “Inspired by non-voltage-gated ion channels, which rely on weak yet persistent ionic modulations to regulate long-term biological functions, the device operates beyond binary switching and supports the formation of multi-level PNV conductance states governed by intrinsic ionic dynamics. Importantly, these PNV states are not static storage states but arise from gradient-driven transient and relaxation processes, whose stability and repeatability are experimentally verified under repeated electrical cycling and thermal perturbation (Supplementary Figs. 30-34).”\nSupplementary Fig. 19. Statistical analysis of second-order relaxation dynamics across 32 devices. a Mean relaxation curve averaged over 32 devices after programming with a 9 V, 100 ms voltage pulse, measured at a read voltage of 1 V. b Statistical distribution of the extracted fast (τ ) 1 and slow (τ ) bi-exponential relaxation time constants across 32 devices, presented as box plots with 2 individual data points overlaid. c Statistical distribution of the ratio τ /τ .\n2 1\nAs shown in Supplementary Fig. 19, the statistical evaluation of the decay stability across 32 devices was performed. The ensemble-averaged conductance decay can be well fitted by a bi-exponential function, yielding distinct fast (τ ) and slow (τ ) 1 2 relaxation components. Here, τ is associated with rapid local relaxation processes, 1 while τ reflects slower ion redistribution governed by the intrinsic concentration 2 gradient. The ratio τ /τ is extracted as a quantitative indicator of the relative 2 1 contribution between gradient-constrained diffusion and fast local relaxation. The relatively narrow distribution of τ /τ across devices indicates that the decay dynamics 2 1 and their relative weighting remain consistent from device to device, confirming that the relaxation behavior is intrinsic to the device architecture rather than arising from device-specific variations.\nSupplementary Fig. 30. Endurance characteristics evaluated using a conventional bipolar switching scheme. a Endurance performance measured under a conventional bipolar voltage switching scheme over 10,000 consecutive cycles (V = 9 V, t = 50 μs; V = ‒9 V, t = 80 positive d negative d μs, Δt = 20 μs), the interval between consecutive cycles was fixed at 150 μs. b Enlarged view of several representative cycles, illustrating the detailed conductance response to the bipolar programming pulses and the corresponding read operation (1 V).\nThe stability of the device under repeated electrical excitation was first examined using conventional bipolar voltage cycling with alternating positive and negative pulses. This measurement provides a baseline assessment of the robustness of the transient conductance response under repeated electrical perturbation. As shown in Supplementary Fig. 30, the transient conductance evolution remains highly reproducible over 104 consecutive cycles, with no observable cumulative degradation. The persistence of the transient response indicates that the device maintains stable dynamical behavior under repeated bipolar excitation.\nTo further examine the stability of the intrinsic-gradient-governed dynamics, an asymmetric unipolar cycling scheme based on the U-SVDP scheme was employed. In this scheme, a high-amplitude pulse V was applied to activate oxygen ions High and perturb the ionic distribution away from equilibrium, followed by a low-amplitude pulse V that probes the subsequent dynamical response under gradient-constrained Low conditions.\nUnder short-pulse excitation (t =50 μs), the oxygen ions activated by V d High rapidly redistributed and relaxed. As a result, the subsequent V pulse does not induce Low a persistent conductance reconstruction but instead produces a transient response governed by the intrinsic concentration gradient. Notably, the modulation polarity induced by V is opposite to that induced by V (Supplementary Fig. 31), indicating Low High that the response directionality is dictated by the intrinsic gradient rather than by residual field-activated ions.\nWhen the pulse width of V was increased to 1 ms, a substantially larger High population of diffusible oxygen ions remained available for redistribution (Supplementary Fig. 32). In this case, the V pulse was able to actively engage these Low ions under the constraint of the intrinsic concentration gradient, leading to a more pronounced and reproducible conductance modulation with negative polarity under readout conditions. These behaviors indicate that, while the intrinsic concentration gradient dictates the directionality of the response, the availability of diffusible ions determines the magnitude and persistence of the V -induced conductance modulation.\nLow\nOverall, these observations confirm that the U-SVDP cycling scheme provides an effective means to evaluate the stability of the intrinsic concentration gradient, which constitutes the essential endurance metric for this second-order memristive system.\nSupplementary Fig. 33. Evaluation of U-SVDP-modulated PNV states at different temperatures. a,b Transient conductance responses of V (9 V, t = 100 ms) and V (6 V, t = High d Low d 100 ms, Δt = 5 s) measured at different temperatures. c Conductance evolution of the PNV states under different temperature conditions.\nWithin the temperature range of 40-80°C, the conductance states remain nearly unchanged, indicating robust stability of the intrinsic-gradient-governed transport behavior (Supplementary Fig. 33c). Although elevated temperature facilitates oxygen-ion activation and enhances the transient conductance response, the underlying coordination environment and gradient profile remain largely preserved in this regime. At higher temperatures (100-120°C), the transient conductance level further increases due to enhanced ionic activity. However, the high-temperature environment progressively destabilizes the ZnTPP coordination and perturbs the multi-filament network, leading to reduced retention of the reconstructed PNV states. As a result, an increasingly pronounced conductance drift is observed at elevated temperatures (Supplementary Fig. 33a-b and Supplementary Fig. 34a-d), reflecting partial thermal redistribution of the intrinsic concentration gradient.\nResponse: We appreciate the reviewer for the insightful comment. The reviewer pointed out that, as suggested by Supplementary Fig. 25, even in the homogeneous case, the ionic concentration evolves slowly between 1 s and 50 s, implying that the long retention may not necessarily originate from the intrinsic gradient. We appreciate this observation and would like to clarify that the behavior shown in the updated Supplementary Fig. 26 is distinctly different from the kinetic mechanism caused by the intrinsic oxygen gradient.\nAs the reviewer noted, the homogeneous configuration indeed exhibits a long-term slow diffusion process after the removal electric field. However, this nonvolatile behavior merely reflects a passive first-order memristive relaxation, where the conductance monotonically decays toward a new equilibrium state without providing any temporal plasticity or re-programmability.[1,2] In other words, the “slowness” in the homogeneous system only extends the relaxation timescale, but it lacks a controllable dynamical component. Consequently, no effective modulation of the learning rate can be achieved (i.e., it’s hard to modulate the strong bulk conductive filament/channel), and the U-SVDP pulse scheme has no operational space under such a static relaxation process.\nAs shown in Supplementary Fig. 26a-b, the ionic migration in the homogeneous distribution is dominated by a strong local electric-field gradient once the external bias is applied (0.1-0.2 s). The initially uniform O2‒ ions configuration becomes destabilized by this abrupt field, inducing random and spatially uncorrelated ionic motion. Upon field removal (> 0.2 s), the steep concentration difference instantaneously drives a rapid\n64\nback-diffusion toward a new equilibrium state, as visualized in the simulated concentration profiles of Supplementary Fig. 26c-d and f-i. Correspondingly, the temporal conductance G(t)/G exhibits a sharp rise followed by a fast decay into a 0 steady plateau, characteristic of a first-order memristive relaxation process (Supplementary Fig. 26e). This transient behavior indicates that the homogeneous system lacks the spatial or energetic constraints required to sustain a distributed, time-correlated evolution of ionic states. Moreover, this stable state means that there are rarely diffusing O2‒ ions for the U-SVDP scheme to modulate, and PNV states cannot be generated. Consequently, the relaxation is fast, monotonic, leading to a purely non-volatile response that cannot support gradual adaptation or long-term modulation.\nSupplementary Fig. 26. Simulation of homogeneous ionic migration and corresponding conductance evolution. The simulated O2‒ distribution at a 0.1 s, b 0.2 s, c 1.0 s and d 5.0 s. e Corresponding G(t)/G curve, where the G is the resting state and G(t) is the temporal conductance. 0 0 The simulated O2‒ distribution at f 100 s, g 200 s, h 400 s and i 600 s.\nIn continual RL, the agent must balance plasticity for adapting to new stimuli with stability for retaining prior knowledge. In the homogeneous system, abrupt ionic redistribution produces an excessively high yet short-lived learning rate α(t), enabling fast adaptation but causing rapid forgetting. Conversely, the intrinsic-gradient naturally maintains a plasticity-stability equilibrium, where the built-in field preserves stored states and slow ionic relaxation introduces a gradual evolution of α(t).\nTherefore, the significance of the intrinsic gradient mechanism lies not merely in “slowing down” the relaxation, but in transforming the slow process into a programmable and reconfigurable kinetic channel with U-SVDP. The built-in potential and interfacial coordination jointly maintain a residual ionic activity that bridges transient diffusion and steady nonvolatility. This time-resolved and state-distinguishable dynamic enables continuous, physically embedded learning-rate\n68\nmodulation at the device level — providing the material foundation for hardware-sustainable continual RL.\nIn response, we have clarified this point in the revised manuscript and added additional simulation data (Figure 3, Supplementary Note 2, and Supplementary Figs. 25-26) comparing the conductance-time dynamics between the homogeneous and intrinsic-gradient configurations, highlighting their distinct relaxation characteristics and their implications for sustainable learning architectures.\nIn order to more accurately capture the underlying physical processes, we have revised the finite-element simulation model and updated the corresponding description of Supplementary Note 2 in the revised Supplementary Information.\nSupplementary Note 2 Finite-element simulation of ionic migration under intrinsic concentration gradient\nModel overview\nFinite-element simulations were conducted to investigate the migration and relaxation behavior of oxygen ions within the ALD-AlO layer. The model was based x on a coupled Nernst-Planck-Poisson framework that incorporates both the external electric field and an intrinsic built-in potential arising from the experimentally observed oxygen concentration gradient. This intrinsic potential regulates the ionic drift and diffusion processes, enabling the system to reproduce the characteristic slow relaxation behavior observed in the experiment. The governing equations are expressed as:\n∂c\n=J (11)\n∂t\nJ =D(t)cμc(ϕ +ϕ ) (12)\next int\nwhere c is the oxygen-ion concentration, D(t) is the time-dependent diffusion coefficient, μ is the ionic mobility, and ϕ and ϕ denote the externally applied and ext int intrinsic potentials, respectively. The intrinsic potential satisfies:\ndϕ D d(lnc )\nint = eq (13)\ndy μ dy\nensuring that in the absence of an external bias, the intrinsic equilibrium profile c (y) eq corresponds to a zero-flux condition (J=0). The intrinsic oxygen concentration gradient\n69\nused in the simulation was parametrized from the depth-resolved characterization of the ALD-AlO /ZnTPP interface. The equilibrium concentration was fitted by an x exponential form:\ny\nc (y)=C +(C C )exp( ) (14)\neq top bottom top λ\nwhere C =0.0538 mol/L, C =0.0322 mol/L, and λ=1.0 nm correspond to the bottom top oxygen concentrations at the bottom and top interfaces extracted from the experimental data. This measured c (y) was used as the target equilibrium distribution throughout eq all simulation stages, ensuring that the system relaxes back to the experimentally determined intrinsic state once the external bias is removed.\nTime-Dependent Diffusion Coefficient and Physical Basis\nIn a conventional homogeneous simulation, the ionic diffusion coefficient D is typically treated as a constant (e.g., ~10‒22m2/s[3]), yielding a single-rate exponential relaxation of conductance. However, such a single-timescale model could not reproduce the two-stage (fast to slow) conductance decay observed in the intrinsic-gradient devices. We therefore implemented a time-dependent diffusion coefficient D(t) to capture the transition from an initial fast ion redistribution to a subsequent slow relaxation. In this model, D(t) is expressed as a time-dependent parameter:\nt\nD(t)= D + (D D )exp( ) (15)\nslow fast slow τ\nD\nwhere D (=5×10‒20m2/s) represents the effective diffusion coefficient immediately fast after the bias removal, when the concentration gradient is steep and the field-driven redistribution dominates. D (=5×10‒21m2/s)describes the long-term diffusion regime slow governed by the built-in potential and interfacial coordination effects. The parameter τ D defines the characteristic time for the transition between the two diffusion regimes. The transition time constant τ =36 s (extracted from the Supplementary Fig 23) defines the D characteristic timescale for the reduction of effective diffusivity, consistent with the ≈30 s onset of the slow relaxation observed in the conductance decay curve (as shown in Supplementary Fig 25). This time-dependent formulation naturally describes the progressive reduction of effective ionic mobility as the system evolves toward\n70\nequilibrium, leading to the bi-exponential relaxation behavior observed experimentally.\nThe selection of dual diffusion coefficients is physically motivated by both the device’s operational environment and literature-reported oxygen ion diffusivities in amorphous Al O . Under applied bias, the strong electric field significantly enhances 2 3 ionic transport by lowering the migration barrier of O2⁻ ions. Ab initio studies have demonstrated that increasing electric field strength in alumina reduces the activation energy for oxygen diffusion, thereby accelerating ionic motion.[4] This field-driven enhancement is reflected in the choice of D =5×10‒20m2/s, a value that falls within the fast experimentally observed range for oxygen mobility in amorphous Al O under 2 3 moderate-to-strong electrical stress.[5]\nOnce the external field is removed, ionic transport is no longer driven by the applied potential but instead regulated by intrinsic factors such as the built-in oxygen concentration gradient and thermal energy.[6] In this post-bias regime, the ion migration becomes progressively restricted as the driving forces diminish and interfacial effects, including oxygen coordination and trapping at the oxide-metal interface, start to dominate. At this stage, interfacial effects such as oxygen coordination (e.g., Al–O or Zn–O bonding) and Coulombic trapping significantly reduce mobility. Reported diffusivities in dense AlO thin films can be as low as 10‒22m2/s, and amorphous x structures generally exhibit higher transport due to their disordered networks.[3,7] Accordingly, the D =5×10‒21m2/s was assigned to represent long-term constrained slow diffusion under built-in field and interfacial coordination effects.\nThe model captures two distinct diffusion regimes. The initial phase features high ionic mobility due to electric-field-assisted drift, while the later stage reflects reduced mobility governed by intrinsic material properties and interfacial coordination. This dual-regime behavior is also consistent with device-scale relaxation dynamics. Based on Fick’s law over a 10 nm oxide layer, a decay occurring over tens of seconds corresponds to an effective diffusivity near 10−20m2/s whereas a slower relaxation spanning several minutes corresponds to diffusivities around 10−21m2/s. These values closely match the chosen parameters, supporting their physical validity in describing\n71\noxygen transport within ALD-AlO layer.\nx\nMapping of ionic distribution to conductance response\nTo establish a quantitative correspondence between the simulated oxygen-ion dynamics and the experimentally measured conductance evolution, the local ionic concentration c(y,t) was converted into the temporal conductance G(t) using an exponential mapping function:\nG(t)=G exp[β(c (t)c (0))] (16)\n0 int int\nwhere G represents the initial conductance prior to the electrical stimulus, and c (t) 0 int denotes the spatially averaged oxygen-ion concentration near the active interface (typically the bottom electrode region). The proportionality coefficient β captures the sensitivity of the tunneling barrier height or defect-related conduction channel to the local ionic accumulation, reflecting the exponential dependence of conductance on oxygen vacancy density in a barrier-type memristive structure.\nThis exponential form is physically justified by the Schottky emission and tunneling mechanisms dominating the device conduction, where local ionic redistribution modulates the interfacial barrier ϕ(t) through the relation:\nϕ(t)qN (t) (17)\nO2\nwhere the q is the elementary charge (1.602×10‒19C), and N (t) is the local oxygen- O2 ion density near the interface at time t. The term qN (t) thus represents the O2 interfacial negative charge induced by accumulated O2‒ ions, which modulates the potential barrier and hence the electronic conductance. Therefore, an increase in interfacial oxygen-ion concentration effectively lowers the barrier height and increases conductance, while subsequent back-diffusion leads to a gradual restoration of the barrier and a decrease in G(t).\nIn the simulation, β was empirically adjusted to reproduce the experimentally observed ratio between high and low conductance states without affecting the intrinsic relaxation time constants. This mapping enables direct comparison between the numerically simulated ionic relaxation and the measured conductance-time (G-t) curves, validating that the two-stage (fast-slow) diffusion dynamics can quantitatively\n72\nreproduce the observed bi-exponential decay behavior.\nIn response to the reviewer’s concern about the clarity of the intrinsic concentration-gradient effect, we have refined the physical model and provided a revised interpretation of the oxygen-ion dynamics. We have also updated Figure 3 in the revised manuscript (page 9) and Supplementary Figures 25-26 (pages 33-34) to reflect these improvements.\nOn page 9 of the revised manuscript, “The dynamic behavior of O2− ions was simulated using finite-element simulation, with model details provided in Supplementary Note 2 and Supplementary Fig. 24. Under applying a 100 ms 𝐸⃑ , ext regions with higher initial oxygen concentration become further enriched, while depleted regions become more deficient, forming a strongly asymmetric ionic profile that stores additional potential energy in the built-in landscape. (Fig. 3f-g and Supplementary Fig. 25e). After removing the 𝐸⃑ , the intrinsic gradient provides a ext kinetic buffering effect that slows down the subsequent relaxation (Fig. 3h and Supplementary Fig. 25a-d). As shown in Supplementary Fig. 25f-g, the conductance does not drop immediately but enters a gradual transition stage lasting for about 30 s. During this period, the steep interfacial gradient is slowly released, and the residual built-in potential continues to regulate the ionic motion, maintaining a non-equilibrium but stable transport pathway. It could be observed that the concentration evolution process was still performed even after hundreds of seconds (Supplementary Fig. 25h-k and Fig. 3i). This buffering effect effectively converts the rapid diffusion process into a operable slow dynamic that preserves partial ionic activity for an extended second-order variable operation timescale, which is different with the homogeneous distribution (Supplementary Fig. 26).”\nOn page 34 of the Supplementary Information, “As shown in Supplementary\nFig. 26a-b, the ionic migration in the homogeneous distribution is dominated by a strong local electric-field gradient once the external bias is applied (0.1-0.2 s). The initially uniform O2‒ ions configuration becomes destabilized by this abrupt field, inducing random and spatially uncorrelated ionic motion. Upon field removal (> 0.2 s), the steep concentration difference instantaneously drives a rapid back-diffusion toward a new equilibrium state, as visualized in the simulated concentration profiles of Supplementary Fig. 26c-d and f-i. Correspondingly, the temporal conductance G(t)/G 0 exhibits a sharp rise followed by a fast decay into a steady plateau, characteristic of a\n76\nfirst-order memristive relaxation process (Supplementary Fig. 26e). This transient behavior indicates that the homogeneous system lacks the spatial or energetic constraints required to sustain a distributed, time-correlated evolution of ionic states. Moreover, this stable state means that there are rarely diffusing O2‒ ions for the U- SVDP scheme to modulate, and PNV states cannot be generated. Consequently, the relaxation is fast, monotonic, leading to a purely non-volatile response that cannot support gradual adaptation or long-term modulation.”\nReference\n[1] S. Kumar, X. Wang, J. P. Strachan, Y. Yang, W. D. Lu, Nature Reviews Materials 2022, 7, 575.\n[2] S. Kim, C. Du, P. Sheridan, W. Ma, S. Choi, W. D. Lu, Nano Letters 2015, 15, 2203.\n[3] T. Nabatame, T. Yasuda, M. Nishizawa, M. Ikeda, T. Horikawa, A. Toriumi, Japanese Journal of Applied Physics 2003, 42, 7205.\n[4] Z. An, S. Sun, B. Dong, Scientific Reports 2024, 14, 10454.\n[5] D. Leon-Chaparro, M. D. Nguyen, C. Baeumer, G. Mul, G. Katsoukis, Advanced Materials Interfaces 2025, 12, 2400846.\n[6] R. Nakamura, T. Toda, S. Tsukui, M. Tane, M. Ishimaru, T. Suzuki, H. Nakajima, Journal of Applied Physics 2014, 116.\n[7] Y. Yang, X. Zhang, L. Qin, Q. Zeng, X. Qiu, R. Huang, Nature Communications 2017, 8, 15173.\nResponse: We appreciate the reviewer for the valuable comments. Indeed, the coordination interaction between O2‒ ions and ZnTPP molecules influences the retention characteristics by introducing local binding/dissociation dynamics. In our earlier work (Adv. Mater., 2021, 33, 2104370), we systematically analyzed this effect and confirmed that reversible Zn–O coordination can stabilize conductance through partial charge transfer and localized trapping. In this framework, the ion migration rate follows an Arrhenius-type dependence on the Me–O bond energy E . The stronger Me–O coordination (larger E ) would lower the hopping frequency and enhances retention Me–O theoretically.\nIn this scenario, the coordination kinetics acts as a local trapping mechanism, giving rise to semi-permanent, quasi-nonvolatile conductance states. However, in practical operation, this coordination process behaves quasi-statically rather than purely dynamic, since the O2‒ coordination within the MTPP matrix is strongly field-dependent. Once an electrical history is established, the coordinated configuration becomes semi-permanent, and the device cannot be fully reset to its pristine state simply by applying a reversed-polarity voltage.\nOn this basis, it is precisely the quasi-permanent coordination-induced nonvolatility that allows us to concentrate on engineering the oxide layer with an intrinsic oxygen concentration gradient. In this configuration, the ZnTPP layer acts as a stable, semi-permanent conductive matrix, while the switching dynamics are governed by the oxygen migration and redistribution within the oxide layer.\nFurthermore, this quasi-permanent coordination feature provides a stable electronic background that greatly enhances the tunability of the coupled dynamics between the interfacial barrier and the filamentary conduction. By providing a consistent electronic background and suppressing random interfacial fluctuations, the ZnTPP coordination enables a deterministic coupling between fast ionic drift and slow barrier recovery, thereby facilitating the controlled implementation of the unique second-order pseudo-\n78\nnonvolatile (PNV) behavior described in this work.\nIn response, we emphasized the role of the ZnTPP coordination layer and revised the description accordingly.\nOn the revised manuscript page 6: “When under the large voltage (V ≥8 V), a High large number of O2− would penetrate into the ZnTPP layer from the AlO layer, and x successively hop from one ZnTPP molecule to another through binding and dissociation (Supplementary Fig. 13).[40] Therefore, the field-dependent ZnTPP coordination layer acts as a stabilizing electronic foundation, enabling a deterministic coupling between ionic drift and interfacial-barrier recovery.”\nResponse: We appreciate the reviewer’s valuable suggestion to clarify the definition and physical meaning of the pseudo-nonvolatile (PNV) state. As introduced in the manuscript (Fig. 4j, Fig. 5d-e), “The PNV state is characterized by conductance that remains stable and distinguishable for a short period”.\nIn general, the volatile conductance state is characterized by a rapid and reversible response to external stimuli, where the conductance rapidly returns to its initial level once the stimulus is removed, demonstrating a short-term and transient memory behavior. The nonvolatile conductance state maintains its resistance or conductance value even after the removal of external stimuli, thus enabling long-term stable information storage through permanent or quasi-permanent structural or compositional changes.\nIn comparison, the pseudo-nonvolatile (PNV) state does not simply represent a slow relaxation process. Instead, it describes a transient, metastable conductance equilibrium, which forms when the system temporarily stabilizes between adjacent energy minima. In this regime, local conductive filaments or quantum bridges reach a transiently frozen configuration that is neither completely ruptured nor fully stabilized. This configuration yields discrete and temporally stable conductance plateaus that can be clearly distinguished from continuously decaying relaxation states. Such PNV states were observed in filamentary memristors (Nat. Commun., 2023, 14, 7655; Mater. Horiz., 2020, 7, 71), conductance quantization devices (Adv. Mater., 2022, 34, 2201248), and phase transition systems (Sci. Adv., 2017, 3, e1700160), but they are rarely explicitly defined as a distinct physical regime.\nIn this work, these states originate from a dual-variable coupled dynamic system governed by the interplay between the interfacial barrier height (Φ) and the multi-conductive filament system. The interfacial barrier slowly evolves with the intrinsic oxygen-ion concentration gradient, forming an interface-dominated temporal memory,\n80\nwhereas the multiple, spatially distributed oxygen-vacancy filaments undergo quasi-stable aggregation and partial dissolution, demonstrating a bulk-dominated quasi-equilibrium. In this multi-filament system, each filament possesses a slightly different activation energy and relaxation constant, creating a distributed ensemble of conductive pathways. The coupling between the slow barrier recovery and the asynchronous relaxation of individual filaments establishes transient equilibrium states that appear as discrete yet consecutive multi-level conductance plateaus.\nIn this coupled Φ-filament dynamic framework, the unipolar spike-voltage-dependent plasticity (U-SVDP) serves as a finely tuned method that selectively regulates both the interfacial barrier and filamentary conduction under a single-polarity electric field. As shown in Fig. R2a-d, the V excitation redistributes oxygen ions and High modulates the local barrier profile without triggering irreversible filament rupture, while the subsequent V partially reconfigured conductive channels. This unipolar Low modulation enables synchronized adjustment between fast ionic drift and slow barrier relaxation, allowing the system to enter a transiently frozen equilibrium manifested as a pseudo-nonvolatile (PNV) conductance plateau. Compared with conventional bipolar pulse programming, the U-SVDP approach achieves more controllable and energy-efficient multi-level modulation by maintaining the device within the metastable regime rather than driving it into a fully nonvolatile state.\nFig. R2. Construction and unipolar short-voltage-driven reconfiguration of pseudo-\n81\nnonvolatile (PNV) states in the interfacial barrier height (Φ) and multi-conductive filament system. a The multifilament network is triggered by a brief V pulse that induces localized High oxygen-ion redistribution and partial modulation of the interfacial barrier. b After the excitation, the system relaxes into a metastable conductance state (G ) dominated by slow ionic diffusion within 0 the coupled Φ-filament dynamics. c Sequential unipolar V pulses act as mild perturbations that Low partially reconfigure the metastable conductive channels. d Establishing stable PNV conductance plateaus through unipolar modulation, corresponding to the gradual reduction of G .\nIn summary, the PNV state originates from the intrinsic oxygen-concentration gradient that mediates a dynamic equilibrium between fast ionic drift and slow interfacial-barrier recovery. The intrinsic oxygen gradient not only provides a restoring potential that counterbalances the external field but also establishes a spatially correlated feedback channel for ionic redistribution, which is indispensable for\n82\nsustaining metastable conductance equilibria. Through U-SVDP, these equilibria can be selectively reconfigured within a single polarity window, whereas they are inaccessible to homogeneous systems lacking intrinsic gradients. Consequently, our device uniquely achieves controllable, energy-efficient PNV modulation, establishing a physically grounded analogue to intermediate-term synaptic plasticity.\nAccording to the reviewers’ comments, we have made the following adjustments to the descriptions in the revised manuscript:\nOn page 8: “When the U-SVDP operation resembled the read operation (e.g., 9 V/2 V), ΔG followed an exponential decay process, as shown in Supplementary U-SVDP Fig. 22. This unipolar modulation synchronizes the fast ionic drift with the slow interfacial-barrier relaxation under the guidance of the intrinsic oxygen concentration gradient, forming a dual-variable interaction that dictates the device’s temporal conductance dynamics and prepares the ground for higher-order metastable behavior.”\nOn page 10: “The 1000 ms V induced large extra gradient changes that caused High a longer diffusive process, making V operation obtain an unstable state. By contrast, Low the 100 ms V stimulation resulted in a slower gradient recovery, allowing V to High Low stabilize the decay conductance into a pseudo-nonvolatile (PNV) state. The PNV state refers to metastable, temporally distinguishable conductance plateaus that persists for a finite period after the stimulus is removal, and this state originates from the coupled dynamics of the interfacial barrier and the multi-filament network. Unlike conventional relaxation-driven decay, the PNV state represents a transiently frozen equilibrium between adjacent energy minima, arising from the asynchronous relaxation of multiple conductive filaments and the slow recovery of the interfacial barrier.”\nMinor issues:\n1. Page 4, Line 7, “Additionally” is redundant.\n2. Both ΔT and Δt are used to represent time interval, which could be confusing, and should be unified into Δt.\n3. In the last two graphs of fig.5e, the numbers on x-axis are overlayed.\n4. There is no Fig.3n, but it was mentioned in text and in the caption of Fig. 3. Besides, Fig. 3j is not mentioned.\n5. Fluidic memristor as a new kind of memristor should be mentioned in the introduction part (typical ref., Science 2023, 379, 156-161, and 161-167)\nResponse:\n1. Page 4, Line 7, “Additionally” is redundant.\nResponse: We thank the reviewer for the valuable comments. The redundant word “Additionally” has been removed in the revised manuscript to improve clarity and conciseness.\n2. Both ΔT and Δt are used to represent time interval, which could be confusing, and should be unified into Δt.\nResponse: We appreciate the reviewer’s attention to notation consistency. All time-interval abbreviations have been unified to Δt throughout the manuscript, figures, and Supplementary Information to maintain uniform terminology.\n3. In the last two graphs of Fig. 5e, the numbers on x-axis are overlayed.\nResponse: We thank the reviewer for the valuable comments. We have carefully re-plotted Fig. 5e to improve the axis readability. The x-axis labels are now properly spaced, and the figure layout has been adjusted for clarity in the revised version. As shown below:\nU-SVDP Decay\n85\n4. There is no Fig. 3n, but it was mentioned in text and in the caption of Fig. 3. Besides, Fig. 3j is not mentioned.\nResponse: We thank the reviewer for the valuable comments. The figure number of Figure 3 in the manuscript should only be up to m, and all figure panels (a-m) have been carefully cross-checked to ensure proper correspondence between the manuscript and figure captions. The missing reference to Fig. 3j has also been added in the relevant paragraph (on Page 9).\n5. Fluidic memristor as a new kind of memristor should be mentioned in the introduction part (typical ref., Science 2023, 379, 156-161, and 161-167).\nResponse: We thank the reviewer for this valuable suggestion. The section discussing emerging memristive mechanisms in the introduction has been expanded to include fluidic memristors.\nOn page 3, “Reinforcement learning (RL) algorithms are considered one of the most promising candidates for continual learning realization in computational systems.[4-6] The emergence of solid-state and fluidic ionic memristors that emulate biological ion-based signal processing has further advanced the potential of RL in this domain.[7-14]”\nReference\n[13] T. Xiong, C. Li, X. He, B. Xie, J. Zong, Y. Jiang, W. Ma, F. Wu, J. Fei, P. Yu, L. Mao, Science 2023, 379, 156.\n[14] P. Robin, T. Emmerich, A. Ismail, A. Niguès, Y. You, G. H. Nam, A. Keerthi, A. Siria, A. K. Geim, B. Radha, L. Bocquet, Science 2023, 379, 161.\nReviewer #1:\nThe revision has been well conducted, and I recommend publishing in its current form.\nResponse: We sincerely thank the reviewer for this positive evaluation and for recommending publication of our manuscript in its current form. We are pleased that the revisions have satisfactorily addressed the previous concerns and appreciate the Reviewer’s recognition of the improvements made to the manuscript.\nReviewer #2:\nIn the revised manuscript entitled “Intrinsic Gradient Oxygen-Driven Second-Order Memristors for Continual Reinforcement Learning”, the authors further evaluated the stability and further discussed physical working principle of the device based on the original version. We do acknowledge the great efforts made by the authors to improve the quality of the manuscript that all confusing points were carefully discussed in the revised manuscript. I think this work deserve publication in Nat. Commun. and published ASAP.\nResponse: We sincerely thank the reviewer for the thorough and encouraging assessment of our revised manuscript. We appreciate the reviewer’s recognition of our additional evaluation of device stability and the expanded discussion of the physical working principles. We are particularly grateful for the acknowledgment that the previously confusing points have been carefully clarified in the revision. We thank the reviewer for the strong support of this work."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The revision has been well conducted, and I recommend publishing in its current form", "reviewer_2": "In the revised manuscript entitled “Intrinsic Gradient Oxygen-Driven Second-Order Memristors for Continual Reinforcement Learning”, the authors further evaluated the stability and further discussed physical working principle of the device based on the original version. We do acknowledge the great efforts made by the authors to improve the quality of the manuscript that all confusing points were carefully discussed in the revised manuscript. I think this work deserve publication in Nat. Commun. and published ASAP."}, "author_response": "Response: We sincerely thank the reviewer for this positive evaluation and for recommending publication of our manuscript in its current form. We are pleased that the revisions have satisfactorily addressed the previous concerns and appreciate the Reviewer’s recognition of the improvements made to the manuscript.\nResponse: We sincerely thank the reviewer for the thorough and encouraging assessment of our revised manuscript. We appreciate the reviewer’s recognition of our additional evaluation of device stability and the expanded discussion of the physical working principles. We are particularly grateful for the acknowledgment that the previously confusing points have been carefully clarified in the revision. We thank the reviewer for the strong support of this work."}]

Decision letter