A scalable kinetic Monte Carlo platform for charge transport dynamics in polymer-based memristive systems

Abstract

Polymer assisted ion transport plays a central role in both energy storage technologies and emerging neuromorphic computing devices. Accurately modeling how ions move is crucial for understanding the behavior of batteries and memristors, yet it remains difficult due to the combined effects of drift, diffusion, and electrostatic interactions, along with the limits of continuum models and molecular dynamics. These issues are especially important in the context of the climate and energy crisis, where high performance, low carbon technologies rely on well optimized ion conducting materials and devices. In this work, we present a scalable and flexible stochastic simulation platform based on Markov chain Monte Carlo methods to model ion migration in solid state systems. The platform uses a vectorized, rail based description of device geometry, which allows fast simulations of lateral ion transport and space charge effects while keeping the inherently random nature of ion hopping. It supports a wide range of material systems and can incorporate experimental parameters without requiring changes to the code. We also introduce an implementation optimized for highly energy efficient GPUs, which boosts performance while lowering the carbon footprint of the simulations. Validation with polymer based memristive devices shows that the simulator captures key behaviors such as relaxation decay, current voltage hysteresis, and learning and forgetting dynamics. By combining computational efficiency with relevant mesoscale physics, this platform offers a practical and versatile tool for exploring ion driven processes in energy storage and neuromorphic devices, supporting exploratory and applied research.

---

1 Introduction: the need for scalable, general-purpose ion transport simulation

Currently, scientists are facing interconnected global challenges, including the climate and energy crisis, the electrification of society, and the reduction of CO₂ emissions. At the same time, the rapid expansion of artificial intelligence (AI) and big data has led to algorithmic models of unprecedented complexity. Specifically focusing on this latter challenge, conventional von Neumann architectures suffer from intrinsic bottlenecks (i.e. energy consumption, heat dissipation, and data transfer) that limit their scalability for large scale, parallel workloads.^1,2 Although alternative non-von Neumann approaches have been explored, they remain insufficient for current computing requirements, especially since we need to simultaneously reduce our collective energy consumption. This motivates the search for new paradigms such as in-memory and neuromorphic computing.³

Neuromorphic computing aims to replicate the brain's energy efficient way of processing information by using synthetic materials that emulate neuron like behavior. In recent years, proposals from chemistry, physics, electronics, and materials science have advanced the design of artificial neurons and neural networks.

From the perspective of materials design, ion migration plays a central role in both neuromorphic computing and energy technologies. In batteries and supercapacitors, this phenomenon governs charging rates, cycle life, and processing performance.^4,5 In information technologies, ion migration has become one of the most effective ways to achieve memristive behavior and neuromorphic devices, since controlled ion movement allows for non volatile memory and synaptic plasticity.^6–8 Unlike the relatively free movement of ions in liquids, transport in solid materials occurs through hopping between vacancies, interstitial sites, or disordered paths within crystalline or amorphous structures. Gaining insight into these processes is therefore key to enhancing the performance of both energy storage systems and neuromorphic devices.^9–11

Ion migration in solid state materials results from the combined effects of drift, diffusion, and electrostatic interactions. When an external bias is applied, drift becomes the dominant process, with mobile ions responding to the electric field. Cations move toward the negative electrode, while anions are pushed toward the positive one. As ions accumulate at interfaces, repulsion between like charges increases and creates an internal electric field that counteracts the applied bias. This internal field can weaken the overall driving force for drift, encourage the buildup of space charge, and eventually shift transport toward a diffusion dominated regime.

When the bias is removed, transport is governed mainly by diffusion as ions relax from the high density interfacial regions toward a more uniform distribution. However, Coulombic effects persist: densely packed ions still resist further compaction, slowing the return to equilibrium.¹² These dynamics, marked by a rapid accumulation driven by drift and followed by a slower relaxation shaped by electrostatic repulsion, are commonly observed in memristive devices, perovskite films, and other ionic conductors.^9–11 They have a strong influence on macroscopic properties such as hysteresis in current voltage curves, charge retention times, and the overall efficiency of ion transport.^13,14

However, simulating ion transport remains challenging: Continuum methods, such as drift and diffusion models, are difficult to use routinely because they rely on a large number of parameters that are often hard to determine.¹⁵ In addition, several ion related processes remain poorly understood and are therefore rarely included, especially those involving interactions with transport layers and ion penetration effects.¹⁵ Molecular dynamics can describe atomic scale behavior in detail, but it quickly becomes impractical for large systems or long time windows. In most cases, simulations are limited to a few nanoseconds, while the transport is often governed by rare, long range hopping events that can extend into the microsecond regime.¹⁶

Alternatively, modeling tools like Molecular Dynamics are employed to study ion movement in nanoscale fluid systems, such as investigating ion drift in nanochannel water flows under external electric fields for desalination.¹⁷ By contrast, when it comes to simulating kinetic processes in solid state materials, advanced tools such as kMCpy, a Python package for kinetic Monte Carlo simulations (KMC), have been developed to study ionic transport in crystalline systems. Unlike molecular dynamics, KMC methods are well suited for reaching much longer timescales, up to microseconds or even milliseconds, as well as larger length scales on the order of micrometers. This makes them particularly useful for capturing ion transport behavior in crystalline materials.¹⁶

Building on these ideas, the simulator presented here is fully vectorized to improve efficiency and introduces a rail based abstraction of device geometry based on a dimensionality reduction. In this framework, ions move stochastically along predefined one dimensional paths, referred to as rails, with each rail containing a fixed number of mobile ions. Ion migration, which in real materials can follow curved paths or depend strongly on local morphology, is projected onto these independent one-dimensional rails. This simplification reformulates the previous drift–diffusion approaches but using a stochastic hopping formalism that preserves mesoscale transport statistics while reducing computational complexity. Motion between rails or in the vertical direction is intentionally excluded, which reduces the dimensionality of the problem in a controlled way. This approach preserves device level transport statistics while allowing simulations to scale efficiently. Increasing the number of rails enhances statistical sampling without changing the underlying stochastic rules.

Experimental data can be imported directly through external files, available in both MATLAB and Python. Thanks to flexible rate definitions and a modular structure, it can be adapted to a wide range of materials and transport mechanisms. Validation using memristive device case studies shows that the model captures key experimental trends and offers microscopic insight into resistive switching and ion migration. By striking a balance between efficiency, scalability, and reproducibility, this tool supports the study of ion driven processes in energy storage, non volatile memory, and neuromorphic computing. It also contributes to ongoing efforts toward frugal modeling in molecular electronics.¹⁸

2 Previous theoretical background

At the microscopic level, as mentioned earlier, ionic transport in solids arises from several mechanisms that often coexist and influence each other. When an electric field is applied, drift tends to dominate, pushing ions in a specific direction and causing them to accumulate at interfaces. This accumulation can have a strong impact on device performance and stability.^19,20 Drift is especially sensitive to the enhanced ion mobility that occurs at the very high field strengths typical of nanoscale devices.²⁰ At the same time, hopping, or jump diffusion, allows ions to thermally activate and move between discrete lattice sites through vacancies or interstitials. This process depends on temperature, lattice symmetry, and the local chemical environment, and it can show correlated behavior in disordered regions.^21,22 Essentially, the local energy landscape around a migrating ion governs how easily it can move.²¹ Trapping and detrapping add another layer of complexity. Ions can become temporarily immobilized at deep energy sites, such as defects or grain boundaries, before being released through thermal or electrical activation.^21,23 Finally, in the absence of external fields, diffusion becomes the dominant mechanism, driving ions from regions of high to low concentration according to Fick's laws.^23,24 All of these mechanisms interact to determine the effective ionic conductivity and transport anisotropy. This highlights the importance of simulation approaches that can capture their interplay while remaining computationally practical.

Modeling ionic motion in solid state and membrane systems has traditionally relied on coupled differential equations, most notably the drift-diffusion or Nernst–Planck formulations, sometimes combined with Poisson's equation to account for electrostatic effects. These continuum approaches are computationally efficient, but solving the equations often faces challenges with numerical stability and convergence. They also struggle to capture correlated or collective migration events that play a critical role in device performance,^16,25 since they rely on macroscopic parameters and largely ignore the stochastic nature of ionic random walks.²⁵ These methods become particularly limited when the migration path is curved, such as the route that anion vacancies follow in the perovskite structure,²⁰ or when ion transport shows non linear mobility under very high electric field strengths (e.g., MV cm⁻¹).²⁰ In these cases, the standard linear assumptions of the models no longer hold, because the effect of the applied electric field changes along the curved trajectory of the migrating anion.²⁰

2.1 Differences between charge transport and ion migration in modeling terms

Considering the similarities between polymer-based memristive systems and the classical structure and composition of Light Emitting Electrochemical Cells (LECs), one can draw parallels that help us understand the models developed for the former. This approach allows the extraction of electronic properties that yield neuromorphic behavior.

A key aspect of modeling LECs is the distinction between electronic charge transport and ion migration. The interaction between these two phenomena has historically been challenging to describe, leading to the development of two principal models: the electrodynamic (ED) model and the electrochemical doping (ECD) model.

The ED model, introduced by de Mello and colleagues,^26,27 attributes enhanced charge injection to the formation of electric double layers at the electrodes, which screen the active layer from the external field. In this framework, steady state currents are dominated by electronic diffusion, with negligible drift, and no internal p–n junction is predicted within the bulk material. In contrast, the ECD model, developed by Pei and co-workers,^28–30 proposes that injected charges trigger bulk doping of the conjugated polymer via ion migration. This leads to the formation of expanding p- and n-type regions that eventually meet to create a dynamic p–n junction, where light emission occurs. Evidence supporting this model comes from spatially resolved photoluminescence quenching, scanning Kelvin probe microscopy (SKPM) potential profiling, and transient conductivity measurements.

To connect these perspectives, unified drift-diffusion-Poisson models³¹ have been introduced, which describe electrons, holes, and ions simultaneously through coupled non-linear partial differential equations. While conceptually powerful, these models require careful numerical treatment to maintain stability and convergence, and their complexity grows rapidly when including realistic physical factors such as injection barriers, disorder, recombination, and morphology dependent mobilities are included. As a result, they can become computationally expensive and are less suited for rapid parameter sweeps or adaptation to new experimental conditions.

3 KMC methodology for ion transport

KMC has emerged as a powerful tool for connecting microscopic realism with computational efficiency in ion transport simulations. By reducing atomic dynamics to stochastic, event based hops, with rates derived from transition state theory or first-principles calculations, KMC can capture rare events and long timescale processes that are beyond the reach of conventional molecular dynamics. Its formulation as a Markov chain, in which transitions depend only on the current state, allows exploration of complex energy landscapes through memoryless random walks. This makes KMC particularly well suited for studying resistive switching, filament formation, and diffusion driven dynamics in realistic devices and enables simulation times far longer than those achievable with standard atomistic techniques.^16,25,32,33

Implementing KMC requires a discrete set of ion positions representing energetically distinct sites in the host material, and the migration barriers that separate them. These barriers, typically obtained from first-principles methods (e.g., density functional theory) or the nudged elastic band technique, determine how easily ions can move. The transition rate for each event is calculated using an Arrhenius-type expression:

where ν* is the attempt frequency, E_m the migration barrier, k_B Boltzmann's constant, and T the temperature. In crystalline materials, cluster expansion techniques can further refine site energies and barriers by explicitly accounting for the local chemical environment.^16,25

At each KMC step, all possible migration events, such as hopping, trapping, or detrapping, are listed with their rates. After calculating the total rate, one event is randomly selected according to its probability, and the simulation time is incremented by a stochastic interval drawn from an exponential distribution. By repeating this process over millions of steps, KMC produces statistically robust transport properties, including diffusivity and conductivity, even under conditions far from equilibrium.^16,34 Unlike classical Monte Carlo, which samples equilibrium states, KMC captures time-dependent dynamics directly, making it especially effective for studying ion migration on experimentally relevant timescales.^34,35

Despite these advantages, most existing KMC implementations are specialized and difficult to generalize, which limits reproducibility and broader applicability. Many lack modularity, making it challenging to adapt them to different materials, such as amorphous systems, composites, or defect rich structures.¹⁶ As materials science increasingly relies on high throughput and multiscale modeling, there is a growing need for frameworks that can interface seamlessly with first principles calculations while remaining computationally efficient. Rising computational demands also emphasize the importance of sustainable algorithms: efficient, low overhead simulations not only reduce environmental impact but also make advanced modeling capabilities more widely accessible.

To address these challenges, we introduce a scalable stochastic simulation platform for modeling ion migration in solid state systems. Building on KMC and Markov chain concepts and inspired by our previous work on lanthanide molecular nanomagnets,³⁶ the framework uses a probabilistic, event based methodology adapted to a different class of materials. It integrates systematic event generation, efficient rate calculation, and modular event handling to capture stochastic ion dynamics across a wide range of systems. While it does not reproduce exact atomic trajectories, the approach balances physical realism with computational efficiency, avoids the overhead of coupled PDE solvers, and reveals clear qualitative trends.

3.1 Mapping to the code

To include electrostatic interactions, ion migration is described by probabilities that depend on the position of each ion and change over time according to the local effective fields. At each step, the field is calculated as the sum of the applied bias and the Coulombic contributions from nearby ions. In dense regions, repulsion reduces the net field and lowers the probability of ions moving into that area. From this effective field, the energy cost of a hop is obtained, and the hopping probability is assigned using a Boltzmann factor:

Forward and backward probabilities are normalized to maintain detailed balance in the absence of bias. This allows both drift and diffusion to appear naturally from the same rules, while space charge effects arise from the changes in local fields. When implemented in a vectorized and parallelized way, the simulator can efficiently reproduce the combined effects of drift, diffusion, and Coulomb interactions without the need to explicitly solve Poisson's equation.

Considering the previous ideas, the simulator represents the device as a two-dimensional lattice of size N_y × L_x. The positions of the ions are stored in a matrix , where each element X_ij gives the column position of ion j along rail i. Each row represents a single horizontal path, which prevents ions from moving between rows and allows thousands of ions to be tracked efficiently. Ion motion is determined by an energy difference that depends on the local electric field:

ΔE = t_devq_ionΔV_netd_hop (3)

where ΔV_net accounts for the applied bias and local ion accumulation, d_hop is the hopping distance, and t_dev is the device thickness. To estimate the local electric field, the lattice is divided into ten vertical zones. The local effective voltage in each zone is computed from the difference in ion densities between neighboring zones, scaled by a fixed proportionality factor. This defines a local electric field as the sum of the externally applied bias and an effective Coulombic contribution arising from nearby ionic accumulation. All physical constants entering this electrostatic correction are grouped into a single effective scaling parameter, ensuring dimensional consistency while avoiding the explicit solution of the Poisson equation, as it was mentioned before. Within each zone, the Boltzmann factor determines the difference in probability between hopping in one direction and hopping in the opposite direction. This is then combined with a baseline probability P₀ to define the probabilities of moving left (P_←) and moving right (P_→) in each zone.

Moves are proposed for all ions at the same time at each KMC step, but they are accepted only if the target site is empty. A move is valid only if the target position is not already occupied (Δ_ij > 1), preventing over occupation and reproducing physical space charge limitations. This rule, together with the local electric field that includes Coulomb repulsion, ensures that drift and diffusion emerge naturally while reducing ion mobility in dense regions, without requiring a Poisson solver. To connect the motion of individual ions with measurable properties, the fraction of ions in the first and last sections of the lattice, corresponding to the regions near the electrodes that dominate transport, is tracked. Normalizing these numbers against the total ion count and applying the time vector produces trajectories that can be directly related to device conductance and current under positive or negative bias (see SI, Section S2).

Fig. 1 illustrates the physical abstraction underlying the kinetic Monte Carlo framework used in this work. The polymeric layer is discretized into a set of energetically equivalent trapping sites arranged along independent one-dimensional rails, with ion motion restricted to site-to-site hopping along each rail. Under an applied bias, mobile cations experience a net drift toward the negatively biased electrode, while diffusive relaxation is incorporated by the stochastic nature of the hopping process. As ions accumulate near the interfaces, Coulombic repulsion progressively screens the applied field, lowering the probability of further hops into densely occupied areas.

Fig. 1 Schematic illustration of the kinetic Monte Carlo model used to describe ion transport within the polymeric memristive layer confined between two electrodes (device length = 200 sites). The active layer is modeled as a collection of independent one dimensional transport paths made up of discrete hopping sites. Positively charged mobile ions (blue circles marked with “+”) move only along their assigned path, hopping thermally between neighboring vacant sites. Lateral motion and coupling between adjacent paths are not included. The applied bias between the left (negative) and right (positive) electrodes is distributed across the total number of sites, defining the local potential energy that drives ionic drift, while the stochastic hopping mechanism captures diffusion. Electrostatic interactions are taken into account through density dependent corrections to the local hopping probabilities, which regulate ion redistribution, particularly in the vicinity of the electrodes.

As indicated above, to benchmark the platform, we focused on memristive materials, particularly polymer-based organic devices where resistive switching involves a complex interplay of ionic and electronic processes. In these systems, ion migration couples with charge trapping, interfacial barrier modulation, filament formation, and redox reactions, all of which shape switching kinetics and stability.^37–41 Their behavior, involving multiple mechanisms, makes them an ideal test case for validating the simulator against experimental results, ensuring that it reproduces both qualitative trends and quantitative transport properties.

The rail-based abstraction introduces several simplifying assumptions that define the scope of applicability of the model. First, ion migration is confined to one dimensional motion along independent rails. Lateral diffusion and correlations between neighboring rails are not explicitly included. As a result, morphological inhomogeneities are taken into account only indirectly, through the zone based effective field correction. Electrostatic interactions are also treated in an approximate manner. Instead of solving Poisson's equation self consistently, the model uses an effective density based correction term. This approach captures space charge screening at the mesoscale, but it does not resolve local potential fluctuations at nanometer length scales. In addition, transport paths that may be curved in crystalline or perovskite systems are represented as linear effective trajectories, meaning that variations in field direction along curved paths are not explicitly described. These approximations make it possible to perform large scale stochastic simulations with very low computational cost. At the same time, they can reduce predictive accuracy in systems where lateral correlations, anisotropic migration pathways, or strongly morphology dependent transport processes play a dominant role.

4 Simulation results

The stochastic kinetic Monte Carlo framework developed in this work successfully captures the mesoscale transport behavior of polymer-based memristive devices across different operating regimes. To test the versatility of the simulator, we modeled three representative experimental protocols: relaxation decay, current–voltage hysteresis, and potentiation/depression tests.³⁷ In each case, our goal was not to replicate every microscopic detail of the device, but to capture the dominant mesoscale processes, including drift, diffusion, and space charge effects, that govern device performance. We selected geometry, ionic density, and operating conditions to reflect typical parameters for polymer composite memristive devices reported in the literature, which allowed us to compare with experimental benchmarks while keeping computational cost low.

To further improve computational and energy efficiency, we implemented an alternative version of the simulation using the JAX framework.⁴² Modern GPUs, with their massively parallel architectures, are naturally more energy efficient than CPUs for highly parallel numerical tasks. This efficiency comes from the high throughput per watt achieved by streaming multiple processing units optimized for vectorized operations, which is especially helpful for heavy linear algebra workloads. The JAX implementation closely follows the original NumPy code but takes advantage of GPU acceleration. All array operations were replaced with their JAX equivalents, enabling execution on both CPU and GPU without major structural changes. To further enhance performance, we applied just-in-time (JIT) compilation, which converts computational kernels into optimized machine code and significantly reduces Python overhead during repeated evaluations. This approach provides substantial performance gains while being consistent with the reference NumPy implementation. The quantitative comparison of CPU and GPU performance and energy metrics is provided in Tables S2 and S3 of the SI, showing that GPU execution consistently outperforms the CPU, achieving up to a 17-fold speedup and more than seven-fold energy savings for large scale simulations.

In the relaxation decay protocol, the model begins with rapid ion migration toward the electrode interfaces under an applied bias, followed by a slower redistribution modulated by repulsion once the bias is removed. Tracking individual ion positions over time reveals the evolving spatial density gradients, providing a direct view of how local electric fields weaken as space charge regions develop. Fig. 2 shows the relaxation dynamics after the initial bias is removed.

Fig. 2 Simulated decay of current at different polarization times and voltages. The plotted signal (a.u.) corresponds to the ratio between the number of ions located in the active decile and the total number of ions, computed at each simulation step. Solid curves represent simulations with varying polarization voltages (2.5–10 V) at a fixed polarization time of 4 s, whereas dotted curves correspond to varying polarization times (2–8 s) at a fixed voltage of 5 V. Increasing the polarization voltage or time leads to a higher initial ratio and a slower relaxation process, indicating stronger ionic polarization effects.

The simulations start from a uniform ionic distribution in which approximately one-third of the available sites are initially occupied, with ions randomly positioned across the device. During the polarization stage, a constant voltage is applied for a prescribed duration, and the curves in the figure reflect these two experimental controls: solid lines correspond to different applied voltages (2.5–10 V) at a fixed polarization time, while dotted lines correspond to different polarization durations (2–8 s) at a fixed voltage. Once this charging period ends, the voltage is removed and the system evolves freely, allowing the ionic distribution to relax toward equilibrium. Under these conditions, the fraction of ions accumulated near the electrode interface decreases exponentially with time, which indicates a relaxation process controlled by diffusion once field drift is no longer dominant. Regardless of the applied voltage or polarization time, all curves converge to a similar asymptotic value, demonstrating the system's return to a homogeneous distribution.

It is important to note that, at short enough times, the system can present short-time memory (STM), corresponding to noticeable values of current at times <10 s and applied fields lower than 5 V. Increasing the voltage or the duration of the applied field can turn this effect into long-time memory (LTM), which persists up to 45 s. Regarding the relaxation dynamics, a single polarization phase of a few volts and a few seconds still corresponds to the linear response regime. This is especially true for the time, where we see a noticeable proportionality between the number of ions accumulated in the first decile and time. At moderate voltages (<5 V), the ionic population scales linearly with initial pulse duration, confirming operation within the quasi-linear response regime. Deviations appearing at higher voltages (>7.5 V) reveal the onset of field-induced saturation and Coulombic crowding, which limit further ionic drift and enhance the consequent Coulombic effects. In this context, in the absence of strong interactions between ions, no significant deviations are expected at lower voltages or shorter polarization times. In the present model, however, both exclusion of occupied sites and a locally reduced effective electric field are explicitly included, which progressively suppress ionic drift as the ion density increases. These results show that the model accurately captures the transition from transport dominated by drift to transport controlled by diffusion, a key feature of ion-conducting polymers.

Monitoring the spatial and temporal dynamics of ionic charge carriers is central to the hysteresis simulations, where the program captures how the history of ionic accumulation shapes the width and asymmetry of the I–V loop under different sweep speeds and bias polarities.

The current–voltage characteristics depicted in Fig. 3 further validate the model's ability to capture the interplay between ionic motion and electrical response. The simulated I–V loops exhibit pronounced hysteresis, which is a well-known characteristic of memristive devices. In this case we also observe that the width and asymmetry present a strong dependence on the voltage sweep rate. At slow sweeps (<0.4 V s⁻¹), ions have sufficient time to migrate and accumulate at the interfaces, leading to strong polarization effects and broad hysteresis loops as expected. As the sweep rate increases, the ionic subsystem proves to be unable to follow the rapidly varying alternate electric field, resulting in progressively narrower and more symmetric loops. This rate dependent effect is a physical representation of the competition between electronic conduction and ionic drift over time. Again, a demonstration of the power of the simulator to reproduce memory effects, with the voltage history dependence arising intrinsically from ion redistribution.

Fig. 3 Simulated I–V hysteresis loops obtained at different voltage sweep rates (0.1–0.8 V s⁻¹). Each sweep begins from the final ionic configuration reached at the end of the preceding sweep, thereby capturing the cumulative influence of sweep history on the device response. The simulated current (a.u.) corresponds to the decile-based ionic population used in the model: during forward (positive-voltage) segments, the active population is taken as the first decile, whereas during reverse (negative-voltage) segments it is taken as the tenth decile, following the selection rule prescribed by the simulation algorithm.

Beyond continuous voltage application, the simulator is also able to capture the fine dynamic conductance modulation to resemble synaptic plasticity in typical neuromorphic materials. Fig. 4 shows the potentiation–depression behavior obtained under alternative packages of voltage pulses for two representative relaxation times. The conductance increases during positive stimulation (potentiation phase) and decreases when the bias is reversed (depression phase), forming an asymmetric cycle typical of a memristive learning process. The degree of asymmetry and memory retention specially depends on the relaxation constant: short relaxation times (τ < 1 × 10⁻³ s) produce fast conductance changes with labile persistence of the memory effect, whereas increasing relaxation times (τ > 5 × 10⁻³ s) result in slower, more persistent responses. This tunability highlights how ionic relaxation governs the balance between short- and long-term memory analogs.

Fig. 4 Simulated potentiation–depression responses for two ionic relaxation times (τ = 0.0010 s and τ = 0.0050 s). The upper panel shows the evolution of device conductance during the learning (left) and forgetting (right) phases, with shorter relaxation times producing sharper potentiation and faster decay. The lower panel shows the applied voltage waveform used to drive the dynamics.

5 Prospective extension of the simulator

Validated here for memristive materials, the platform is general and extensible to other device's phenomena.^16,18,25,36 Indeed, the model can be applied to, with the appropriate extensions, to any systems with migration of charged particles driven by an external voltage. A specific feature for the current purpose is the focus on the charge density near the electrodes, which regulates the device's memristance. In other cases, similar specific considerations would be needed. In particular:

Energy storage (supercapacitors): for this application, the model is mostly adequate as it is. The most significant change here would be the metric we use to quantify the phenomenon: for supercapacitors one would like to quantify the energy storage, which could be estimated simply by the (electric) potential energy due to the charge separation.

Optoelectronics: for this application, one would need to add a new different physical process, namely carrier recombination where this specific quantum effect accounts for the light emission of the device. This modification would involve the explicit inclusion of electrons/holes in the MC simulation and the recombination probability/quantum yield.

Interaction between spin and chirality (CISS/Electric Magneto Chiral Anisotropy): for this field one would need to consider the spin of the charge carriers, which would be electrons, and also the effect on the resistance of the combination of spin state, transport direction and chirality of the device. Indeed, some of us already present this in an upcoming manuscript.⁴³

Whereas the current version of the model does not reproduce filament formation, recovering this phenomenon would just require a small modification of the code and a moderate increase in the calculation time. The position of each charge carrier is currently registered, so it would be possible to change the behavior for each individual rail that satisfies the condition of a continuous line of charge carriers starting from the electrode. In those cases, all charge carriers except the last one should be considered as part of the electrode, meaning that the total voltage would act on a shorter device and thus produce a more intense electric potential. In turn, this would naturally generate a faster movement of the charge carriers and naturally result in the lengthening of these rails with a continuous line of charge carriers starting from the electrode, i.e. filaments.

In particular, to consider the systems where non-linear ion mobility is found (such as the perovskytes presented in ref. 20), the simulator should add an additional term on the relative energies of the successive sites which would account for the field dependence. To consider ionic transport along particular crystallographic directions we have projected all the trajectories in a single rail which accounts for a simplified version of the real system. Another simple improvement would be allowing more than one class of rails, with different proportions and hopping parameters: this would allow complex phenomena such as different characteristic times, frequencies or temperatures coexisting in the same device. Furthermore, the modeling could also be made more realistic by allowing motion between rails, although this would significantly increase computation time and affect performance.

Conclusions

The simulator predictions are consistent with key qualitative trends reported in experimental studies of polymer based memristive materials, including sensitivity to different relaxation times, the amplitude of the applied bias, and how the hysteresis loop shape depends on the voltage sweep. Since the implementation is computationally light and scales almost linearly with the number of ions and time steps, large parameter sweeps can be performed in seconds or minutes on a standard workstation. This makes it possible to explore, in a fast and efficient way, how electric field strength, morphology, and material properties jointly affect ionic behavior. As a result, the simulator offers a practical complement to more demanding numerical techniques and, in some cases, may even remove the need for them.

These results show that the KMC framework provides a physically grounded yet computationally efficient description of ionic migration in polymeric memristive materials. By relying on fitting parameters that can be obtained experimentally, the model reproduces essential observables such as relaxation decay, hysteresis modulation, and asymmetric potentiation and depression cycles through a single stochastic mechanism. The agreement between the simulated trends and reported experimental behavior confirms that the model captures the correct hierarchy of time scales governing drift, diffusion, and electrostatic screening. This makes the simulator a valuable exploratory tool for linking device architecture, operating conditions, and material properties to emergent resistive switching behavior, and for supporting the rational design of future neuromorphic components.

Within this context, the model presented in this work a significant improvement. By adopting a probabilistic description of the relevant physical events, we eliminate the need for globally convergent differential equation solvers and enable the direct use of realistic experimental descriptors within the simulations. This approach supports efficient modeling while adding flexibility for incorporating new physical mechanisms without compromising numerical stability. For example, charge trapping effects associated with sparse defects can be introduced by modifying the transition probabilities at a small number of selected sites. Using this MC model, we are able to simulate the spatiotemporal evolution of ionic migration in polymer devices under realistic operating conditions, with a level of robustness and explanatory power that was not accessible in earlier modeling approaches.

Author contributions

Gerliz M. Gutiérrez-Finol: conceptualization, methodology, software, formal analysis. Kirill Zinovjev: methodology, formal analysis. Alejandro Gaita-Ariño: supervision, project administration, validation. Salvador Cardona-Serra: supervision, project administration, validation. All authors: writing – original draft, writing – review and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

All datasets generated and analyzed in this study are openly provided in the accompanying repository. Raw and processed files required for full reproducibility are available in the “data for reproducibility” directory within https://github.com/gerlizg/CupFlow.

All simulation code used in this study is openly available in the CupFlow repository at https://github.com/gerlizg/CupFlow. The repository contains the full kinetic Monte Carlo framework, example scripts, and all configuration files required to reproduce the results presented in this work.

Supplementary information is available. See DOI: https://doi.org/10.1039/d5qm00811e.

Acknowledgements

A.G.A. has been supported by the Generalitat Valenciana (GVA) CIDEGENT/2021/018 grant. A.G.A. thanks grant PID2020-117177GB-I00 funded by MCIN/AEI/10.13039/501100011033 (co-financed by FEDER funds). This study is part of the Quantum Communication programme and was supported by grant PRTR-C17.I1 funded by MCIN/AEI/10.13039/501100011033 and European Union NextGenerationEU/PRTR, and by GVA (QMol COMCUANTICA/010).

References

1. R. Teja Potla, Scalable machine learning algorithms for big data analytics: Challenges and opportunities, J. Artif. Intell. Res., 2022, 2(2), 124–141.
2. S. Dhar, J. Guo, J. (Jason) Liu, S. Tripathi, U. Kurup and M. Shah, A survey of on-device machine learning: An algorithms and learning theory perspective, ACM Trans. Internet Things, 2021, 2(3), 1–49.
3. S. Zhu, T. Yu, T. Xu, H. Chen, S. Dustdar, S. Gigan, D. Gunduz, E. Hossain, Y. Jin, F. Lin, B. Liu, Z. Wan, J. Zhang, Z. Zhao, W. Zhu, Z. Chen, T. S. Durrani, H. Wang, J. Wu, T. Zhang and Y. Pan, Intelligent computing: The latest advances, challenges, and future, Intell. Comput., 2023, 2, 0006.
4. S. C. Sand, J. L. M. Rupp and B. Yildiz, A critical review on li-ion transport, chemistry and structure of ceramic-polymer composite electrolytes for solid state batteries, Chem. Soc. Rev., 2025, 54(1), 178–200.
5. A. C. Radjendirane, D. K. Maurya, J. Ren, H. Hou, H. Algadi, B. B. Xu, Z. Guo and S. Angaiah, Overview of inorganic electrolytes for all-solid-state sodium batteries, Langmuir, 2024, 40(32), 16690–16712.
6. H. Kwak, N. Kim, S. Jeon, S. Kim and J. Woo, Electrochemical random-access memory: recent advances in materials, devices, and systems towards neuromorphic computing, Nano Convergence, 2024, 11(1), 9.
7. S. Park, H.-J. Lee and H. W. Jang, Harnessing ion migration in halide perovskites: Implication and applications, Solid State Ionics, 2025, 422, 116816.
8. Y. Ji, L. Wang, Y. Long, J. Wang, H. Zheng, Z. G. Yu, Y.-W. Zhang and K.-W. Ang, Ultralow energy adaptive neuromorphic computing using reconfigurable zinc phosphorus trisulfide memristors, Nat. Commun., 2025, 16(1), 6899.
9. Z. Chen, T. Du, N. M. A. Krishnan, Y. Yue and M. M. Smedskjaer, Disorder-induced enhancement of lithium-ion transport in solid-state electrolytes, Nat. Commun., 2025, 16(1), 1057.
10. H. Yang and N. Wu, Ionic conductivity and ion transport mechanisms of solid-state lithium-ion battery electrolytes: A review, Energy Sci. Eng., 2022, 10(5), 1643–1671.
11. S. Huang, S. Li, Z. Huang, K. Zhang, W.-L. Song and S. Jiao, A review of multiscale characterization methods of ion transport in solid-state electrolytes, Chin. Chem. Lett., 2025, 110973.
12. A. O. Alvarez, M. García-Batlle, F. Lédée, E. Gros-Daillon, J. M. Guillén, J.-M. Verilhac, T. Lemercier, J. Zaccaro, L. F. Marsal, O. Almora and G. Garcia-Belmonte, Ion migration and space-charge zones in metal halide perovskites through short-circuit transient current and numerical simulations, Adv. Electron. Mater., 2024, 10(11), 2400241.
13. P. Kwantwi-Barima, S. V. B. Garimella, I. K. Attah, X. Zheng, Y. M. Ibrahim and R. D. Smith, Accumulation of large ion populations with high ion densities and effects due to space charge in traveling wave-based structures for lossless ion manipulations (SLIM) IMS-MS, J. Am. Soc. Mass Spectrom., 2024, 35(2), 365–377.
14. I. Wolaska, K. Piwowarski and J. Puton, Conservation of the charge in signal from drift tube ion mobility spectrometers, Anal. Chem., 2024, 96(43), 17337–17344.
15. M. Diethelm, T. Lukas, J. Smith, A. Dasgupta, P. Caprioglio, M. Futscher, R. Hany and H. J. Snaith, Probing ionic conductivity and electric field screening in perovskite solar cells: a novel exploration through ion drift currents, Energy Environ. Sci., 2025, 18(3), 1385–1397.
16. Z. Deng, T. P. Mishra, W. Xie, D. A. Saeed, G. S. Gautam and P. Canepa, kmcpy: A python package to simulate transport properties in solids with kinetic monte carlo, Comput. Mater. Sci., 2023, 229, 112394.
17. F. Sofos, T. Karakasidis and I. E. Sarris, Molecular dynamics simulations of ion drift in nanochannel water flow, Nanomaterials, 2020, 10(12), 2373.
18. G. M. Gutiérrez-Finol, A. Ullah, M. González-Béjar and A. Gaita-Ariño., A call for frugal modelling: two case studies involving molecular spin dynamics, Green Chem., 2025, 27(12), 3167–3177.
19. M. García-Batlle, J. M. Guillén, M. Chapran, O. Baussens, J. Zaccaro, J.-M. Verilhac, E. Gros-Daillon, A. Guerrero, O. Almora and G. Garcia-Belmonte, Coupling between ion drift and kinetics of electronic current transients in mapbbr3 single crystals, ACS Energy Lett., 2022, 7(3), 946–951.
20. D. Kemp and R. A. De Souza, Nonlinear ion mobility at high electric field strengths in the perovskites SrTiO3 and CH₃NH₃PbI₃, Phys. Rev. Mater., 2021, 5(10), 105401.
21. M. Singh, L. S. Tanwar and R. Banerjee, Impact of thermally activated ionic dynamics on the trap-mediated current-voltage characteristics of a mixed-halide hybrid perovskite, J. Mater. Chem. C Mater. Opt. Electron. Dev., 2025, 13(29), 15168–15184.
22. C. Eames, J. M. Frost, P. R. F. Barnes, B. C. O’Regan, A. Walsh and M. S. Islam, Ionic transport in hybrid lead iodide perovskite solar cells, Nat. Commun., 2015, 6(1), 7497.
23. L. Zuo, Z. Li and H. Chen, Ion migration and accumulation in halide perovskite solar cells, Chin. J. Chem., 2022, 861–876.
24. N. Yantara and N. Mathews, Toolsets for assessing ionic migration in halide perovskites, Joule, 2024, 8(5), 1239–1273.
25. P. Canepa, Pushing forward simulation techniques of ion transport in ion conductors for energy materials, ACS Mater. Au, 2023, 3(2), 75–82.
26. J. C. deMello, N. Tessler, S. C. Graham and R. H. Friend, Ionic space-charge effects in polymer light-emitting diodes, Phys. Rev. B Condens. Matter, 1998, 57(20), 12951–12963.
27. J. C. deMello, Interfacial feedback dynamics in polymer light-emitting electrochemical cells, Phys. Rev. B Condens. Matter, 2002, 66(23), 235210.
28. P. Matyba, K. Maturova, M. Kemerink, N. D. Robinson and L. Edman, The dynamic organic p-n junction, Nat. Mater., 2009, 8(8), 672–676.
29. D. L. Smith, Steady state model for polymer light-emitting electrochemical cells, J. Appl. Phys., 1997, 81(6), 2869–2880.
30. Q. Pei, G. Yu, C. Zhang, Y. Yang and A. J. Heeger, Polymer light-emitting electrochemical cells, Science, 1995, 269(5227), 1086–1088.
31. S. van Reenen, P. Matyba, A. Dzwilewski, R. A. J. Janssen, L. Edman and M. Kemerink, A unifying model for the operation of light-emitting electrochemical cells, J. Am. Chem. Soc., 2010, 132(39), 13776–13781.
32. N. Kopperberg, D. Genua Noguera and S. Menzel, 3D kinetic monte carlo-based investigation of the influence of dopants on the reliability of VCM ReRAM, J. Appl. Phys., 2025, 137(17), 175105.
33. S. Aldana, J. Jadwiszczak and H. Zhang, On the switching mechanism and optimisation of ion irradiation enabled 2D MoS2 memristors, Nanoscale, 2023, 15(13), 6408–6416.
34. M. Andersen, C. Panosetti and K. Reuter, A practical guide to surface kinetic monte carlo simulations, Front. Chem., 2019, 7, 202.
35. W. R. Saunders, J. Grant, E. H. Müller and I. Thompson, Fast electrostatic solvers for kinetic monte carlo simulations, J. Comput. Phys., 2020, 410(109379), 109379.
36. G. M. Gutiérrez-Finol, S. Giménez-Santamarina, Z. Hu, L. E. Rosaleny, S. Cardona-Serra and A. Gaita-Ariño, Lanthanide molecular nanomagnets as probabilistic bits, Npj Comput. Mater., 2023, 9(1), 196.
37. C. D. Prado-Socorro, S. Giménez-Santamarina, L. Mardegan, L. Escalera-Moreno, H. J. Bolink, S. Cardona-Serra and E. Coronado, Polymer-based composites for engineering organic memristive devices, Adv. Electron. Mater., 2022, 8(5), 2101192.
38. H. Zhang, L. Li, A. Guo, J. Li, Y.-T. Li, W. Li, M. Yi and L. Xie, P3OT-based organic polymer memristors for artificial synaptic behavior and neuromorphic computing applications, ACS Appl. Electron. Mater., 2025, 7(5), 2001–2011.
39. B. Zhang, W. Chen, J. Zeng, F. Fan, J. Gu, X. Chen, L. Yan, G. Xie, S. Liu, Q. Yan, S. J. Baik, Z.-G. Zhang, W. Chen, J. Hou, M. E. El-Khouly, Z. Zhang, G. Liu and Y. Chen, 90% yield production of polymer nano-memristor for in-memory computing, Nat. Commun., 2021, 12(1), 1984.
40. F. L. Hoch, Q. Wang, K.-G. Lim and D. K. Loke, Multifunctional organic materials, devices, and mechanisms for neuroscience, neuromorphic computing, and bioelectronics, Nanomicro Lett., 2025, 17(1), 251.
41. R. Khan, N. U. Rehman, S. Iqbal, S. Abdullaev and H. M. Aldosari, Resistive switching properties in memristors for optoelectronic synaptic memristors: Deposition techniques, key performance parameters, and applications, ACS Appl. Electron. Mater., 2024, 6(1), 73–119.
42. J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne and Q. Zhang. JAX: composable transformations of Python + NumPy programs, 2018.
43. S. Giménez-Santamarina, A. Mora Martínez, G. M. Gutiérrez-Finol and A. Gaita-Ariño. An expandable kinetic monte carlo platform for modelling electron transport through chiral molecules, arXiv, preprint, arXiv:2601.17525, 2026 DOI:10.48550/arXiv.2601.17525.

---