Three-Dimensional Fully-Coupled Electrical and Thermal Transport Model of Dynamic Switching in Oxide Memristors

Since researchers from the HP Labs published “The missing memristor found” [1] in Nature in 2008, memristors have become a new research frontier, as they have shown many potential applications, including non-volatile memory, hybrid CMOS/memristor 3D circuits, and brain-inspired computing. Oxide-based memristors reported so far have shown a wide spectrum of complex switching current-voltage characteristics [2], ranging from bipolar nonlinear, bipolar linear, unipolar bistable, to the unipolar threshold switching. It is generally recognized that the complex switching behavior in memristors is governed by temporally and spatially intertwined electrical and thermal processes of electrons, holes, and mobile vacancies,  including field drift, Fick diffusion, Soret effect, and Joule heating.  In order to facilitate the understanding of switching physics in memristors, a high-fidelity model is necessary and should be able to simulate the entire switching spectrum. However, no existing models in literature have captured all these processes and their couplings for electrons, holes, and mobile vacancies on an equal footing in realistic 3D device structures.

               In this work, we present a fully-coupled electrical and thermal model that includes all the important processes and treats the transport of electrons, holes, and vacancies together with the Joule heating on an equal footing. Namely, we solve simultaneously the five coupled partial differential equations: the time-dependent continuity equations for electrons, holes, and vacancies, the time-dependent lattice heat equation, and the Poisson equation for all the charged species in three spatial dimensions. The current density for each species includes contributions from the field drift due to potential gradient, Fick diffusion due to concentration gradient, and the Soret effect due to temperature gradient. The heat generation in the temperature equation includes Joule heating due to all carriers. The fully-coupled equations are implemented in Charon [3], an open-source TCAD code developed at Sandia National Labs. This fully-coupled model allows us to investigate the microscopic interplay of field drift, Fickian diffusion, Soret effect, and Joule heating, and to facilitate the understanding of physical switching mechanisms in all relevant oxide memristors.

               We apply the coupled model to simulate the switching process in a 3D filamentary tantalum oxide memristor shown in Fig.1. Since the Ta electrode and conduction filament (CF) are the active regions, we solve the coupled electrical and thermal equations there, and allow the heat to conduct in the Ta₂O₅ and Pt regions. Voltage is applied to the top of Ta electrode and the bottom of CF is grounded. The top and bottom surfaces of the structure are set to 300 K. The rigid point ion model [4] is used to describe the movement of vacancies, which exponentially depends on the temperature and the electric field when the field exceeds a critical value. Thermal and electrical parameters are being calibrated with experimental data as a function of temperature and vacancy density, and reasonable constant values are currently used for the results shown here. 

               Figure 2 shows center view of the simulated time evolution of vacancy density, electron density, and lattice temperature in the Ta and CF regions during ON switching, when 1V is applied to the top of Ta electrode, and the Ta and CF regions have an initial uniform vacancy density of 1x10²¹ cm^-3 and 5x10²⁰ cm^-3 respectively. With the positive voltage applied to the Ta electrode, Joule heating causes an increase of temperature, which is the strongest near the Ta/CF interface, causing vacancies from the Ta electrode to flow into the CF, as seen at 1ms. As temperature further increases, vacancies continue filling in the CF, and by 10ms, the CF is nearly completely filled (CF is assumed to be less conductive than Ta). Time evolution of electron density shows a process similar to vacancies, indicating that vacancies indeed act as n-type mobile dopants. Radial expansion of CF is not observed here, which may be due to a combination of small radius, uniform initial vacancy density, and constant physical parameters used.        

[1] D. B. Strukov, G. S. Snider, D. R. Stewart, and R. S. Williams, Nature 453, 80-83 (2008).

[2] J. J. Yang, D. B. Strukov, and D. R. Stewart, Nature Nanotechnology 8, 13-24 (2013).

[3] G. Hennigan, D. Fixel, J. Castro, and P. Lin, “Charon Parallel Semiconductor Device Simulator”, SAND2010-3905, June 2010.

[4] N. Cabrera and N. F. Mott, Rep. Prog. Phys. 12, 163-184 (1948).

This work was supported by the Laboratory Directed Research and Development program at Sandia National Laboratories. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy's National Nuclear Security Administration under Contract DE-AC04-94AL85000.