Manufaturação industrial
Internet das coisas industrial | Materiais industriais | Manutenção e reparo de equipamentos | Programação industrial |
home  MfgRobots >> Manufaturação industrial >  >> Industrial materials >> Nanomateriais

Um estudo coletivo sobre modelagem e simulação de memória de acesso aleatório resistiva

Resumo


Neste trabalho, fornecemos uma discussão abrangente sobre os vários modelos propostos para o projeto e descrição de memória de acesso aleatório resistiva (RRAM), sendo uma tecnologia nascente que depende fortemente de modelos precisos para desenvolver projetos de trabalho eficientes e padronizar sua implementação em dispositivos. Esta revisão fornece informações detalhadas sobre as várias metodologias físicas consideradas para o desenvolvimento de modelos para dispositivos RRAM. Abrange todos os modelos importantes relatados até agora e elucida suas características e limitações. Vários efeitos adicionais e anomalias decorrentes do sistema memristivo foram tratados, e as soluções fornecidas pelos modelos para esses problemas também foram mostradas. Todos os conceitos fundamentais do desenvolvimento do modelo RRAM, como operação do dispositivo, dinâmica de chaveamento e relações entre corrente e tensão são abordados em detalhes neste trabalho. Modelos populares propostos por Chua, HP Labs, Yakopcic, TEAM, Stanford / ASU, Ielmini, Berco-Tseng e muitos outros foram comparados e analisados ​​extensivamente em vários parâmetros. O funcionamento e as implementações das funções da janela como Joglekar, Biolek, Prodromakis, etc. também foram apresentados e comparados. Novos conceitos de modelagem bem definidos foram discutidos, o que aumenta a aplicabilidade e precisão dos modelos. A utilização desses conceitos traz diversos aprimoramentos nos modelos existentes, os quais foram enumerados neste trabalho. Seguindo o modelo apresentado, modelos altamente precisos serão desenvolvidos, o que ajudará imensamente os futuros desenvolvedores de modelos e a comunidade de modelagem.

Histórico


Esta nova era da computação requer uma tecnologia que seja igualmente capaz de acompanhar seu crescimento. A nova tecnologia deve ser capaz de atender às demandas de desempenho aprimorado e escalonável para atender aos dispositivos do futuro. Memristors, postulado em 1971 [1] por Leon O. Chua, parece cumprir esses requisitos e lançou as bases para novas classes de dispositivos. Memristors, abreviação de "resistores de memória", são dispositivos básicos de dois terminais que lembram seu estado de resistência interna dependendo da história do estímulo de entrada fornecido. Chua imaginou que os memristores são caracterizados por uma relação entre fluxo e carga, que são as integrais de tempo da corrente e da voltagem, respectivamente.

Mais tarde, em 1976, Chua e Kang [2] generalizaram os memristores para incluí-los em uma nova classe de sistemas dinâmicos chamados sistemas memristivos. No final do século XX, o interesse por esses dispositivos diminuiu, apesar de seus muitos benefícios. Isso se deveu em parte aos avanços na tecnologia de circuitos integrados de silício. Mas com o envelhecimento das tecnologias de silício e sua incapacidade de suportar a redução, a busca por dispositivos alternativos de comutação ganhou atração no início do século XXI. Foi igualmente auxiliado pelos avanços no crescimento e caracterização de materiais em nanoescala. Isso invariavelmente leva a um progresso significativo na compreensão da comutação memristiva microscópica.

A tecnologia de memristor teve um grande avanço no ano de 2008, quando Strukov et al. [3] estabeleceu uma ligação entre a teoria e o experimento para seu TiO x dispositivos baseados em Além disso, eles obtiveram uma histerese pinçada na relação corrente-tensão, que é uma das características identificáveis ​​dos sistemas memristivos [4, 5]. Isso abriu a tecnologia de memristor para uma ampla gama de dispositivos seguindo as pegadas do metal / filme de óxido / estrutura metálica. Alguns dos tipos semelhantes de dispositivos populares foram Oxygen RRAM (OxRRAM) [6,7,8,9,10] e Conductive Bridge RAM (CBRAM) [11,12,13] entre muitos outros. Esses dispositivos são geralmente classificados com base em seu mecanismo de comutação.

Memória de acesso aleatório resistiva (RRAM)


O interesse de pesquisa nesses dispositivos emergentes aumentou porque o comportamento memristivo não volátil demonstrado pode ser aproveitado na memória não volátil. Eles estão sendo vistos como alternativas potenciais da tecnologia de memória flash. Com a computação da era atual sendo cada vez mais orientada por dados, tem havido demandas por uma tecnologia de memória que esteja mais em sintonia com os requisitos presentes e futuros. Em comparação com os vários dispositivos emergentes, os dispositivos RRAM são mais escaláveis ​​[14,15,16,17,18], têm alta densidade [19,20,21,22,23,24], consomem baixa energia [25,26,27 , 28,29], são mais rápidos [30,31,32,33], têm maior resistência e retenção [34,35,36,37] e altamente compatíveis com CMOS [38,39,40,41,42]. Dispositivos RRAM são uma das tecnologias de memória não volátil mais populares, com extenso estudo sendo realizado para entender seu mecanismo e desenvolver modelos para realizar a operação do dispositivo e projetar uma estrutura de dispositivo simples e precisa. Os dispositivos são uma estrutura simples de metal isolante de dois terminais (MIM) e alternam entre dois estados de resistência:estado de baixa resistência (LRS) e estado de alta resistência (HRS). Um LRS sugere que o dispositivo está no estado SET ou ON. Um HRS contrastante significa que o dispositivo está no estado RESET ou OFF. Por meio dessa comutação de estados de resistência no dispositivo, o bit de dados é armazenado [43,44,45]. Os dispositivos RRAM podem ser classificados em dispositivos bipolares e unipolares, dependendo da polaridade de comutação. Na comutação unipolar, os dispositivos alternam na mesma polaridade de polaridade, enquanto na comutação bipolar, é necessária a polarização de ambas as polaridades.

Diversas abordagens foram propostas para explicar o mecanismo de comutação de dispositivos RRAM, mas a mais popular e amplamente aceita, para dispositivos RRAM baseados em óxidos binários, é a formação e ruptura de filamentos condutores localizados (CF) pela deriva de íons de oxigênio / vacâncias [9, 16, 46,47,48,49]. O SET / RESET ocorre como resultado da combinação / regeneração dos íons de oxigênio / vacâncias [50,51,52]. Foi demonstrado que o desempenho dos dispositivos RRAM é fortemente afetado pela escolha da camada de óxido ativo [53,54,55]. Uma variedade de sistemas de óxido, como HfO x , TiO x , NiO x , TaO x , ZnO x , etc. [56,57,58,59,60,61,62,63,64,65,66] foram usados ​​para demonstrar o comportamento de comutação resistiva. Tem havido algumas controvérsias se os dispositivos RRAM são realmente dispositivos memristivos. Para deixar clara a posição dos dispositivos RRAM, Chua forneceu esclarecimentos de que eles são, de fato, dispositivos memristivos [67].

Importância da modelagem RRAM


Um aspecto muito importante do desenvolvimento de dispositivos eletrônicos baseados em novas tecnologias de semicondutores é o papel da modelagem. Um modelo preciso e abrangente é de suma importância para entender a operação do dispositivo, projetá-lo para desempenho ideal e verificar se ele atende às especificações exigidas. Vários modelos foram propostos com vários graus de precisão, diferentes recursos e resultados mistos. Portanto, qualquer desenvolvedor que pretenda projetar um modelo robusto e flexível para dispositivos RRAM deve ter informações sobre os métodos experimentados anteriormente e as restrições enfrentadas.

Neste trabalho, discutimos em detalhes todos os recursos e características dos vários modelos RRAM. Modelos gerais de memristor também são considerados para explicar os dispositivos RRAM [67]. Partindo do modelo Chua [1], que fornece os fundamentos dos memristores, discutimos a definição fundamental de memristores. O avanço para memristores e dispositivos RRAM fornecidos pelo modelo HP [3] é discutido em detalhes. São considerados os efeitos da deriva iônica linear, que constituem os fundamentos do mecanismo desses dispositivos, juntamente com os efeitos não lineares [46, 68, 69]. O modelo Pickett-Abdalla [70,71,72], que lançou as bases para modelos baseados na física compatíveis com o SPICE, é abordado em profundidade. Suas várias características que foram adotadas e refinadas pelo modelo Yakopcic [73, 74] também são abordadas.

Modelos que introduziram novos recursos, como efeitos de limiar [75,76,77], tendo a lacuna do filamento como a variável de estado [78,79,80,81], foram revisados. Alguns dos modelos que levam em consideração dispositivos unipolares e efeitos de temperatura [82,83,84] são revisados ​​em detalhes. Também são considerados modelos físicos [85, 86] com base na dinâmica de crescimento do dispositivo. Junto com estes, modelos considerando apenas dispositivos bipolares [87,88,89], mudança de tamanho do CF [90, 91] e muitos outros fatores [92, 93] são levados em consideração. Uma análise concisa de todos os modelos discutidos foi apresentada na Tabela 1.

Vários modelos baseados em implementações de função de janela, como Joglekar [94], Biolek [95], Benderli-Wey [96], Shin [97], Prodromakis [98, 99], etc. também foram considerados para as limitações e restrições em os vários modelos e os métodos usados ​​pelos modelos subsequentes para superá-los foram apresentados de maneira abrangente. O trabalho significativo feito por Wang e Roychowdhury [100] para melhorar a modelagem RRAM também foi revisado em profundidade, pois é um impulso considerável na direção certa para toda a comunidade de modelagem RRAM. Junto com esses exemplos, são discutidos estudos de simulação e verificação dos dispositivos em diferentes plataformas. Esta é a análise mais abrangente relacionada aos modelos RRAM e memristor no estágio atual. A descrição dos modelos foi dividida em aquelas que descrevem dispositivos bipolares e dispositivos unipolares. Os modelos de implementação de função de janela são descritos em uma seção separada.

Anteriormente, houve várias revisões sobre os mecanismos do dispositivo RRAM [46, 101,102,103,104,105], tecnologia de fabricação [106,107,108,109], pilhas de materiais [110,111,112,113] e uma discussão concisa sobre alguns dos modelos presentes naquele momento [114]. Muito recentemente Villena et al. [115] combinou a teoria de toda modelagem RRAM e propôs um modelo de otimização. Neste estudo, nos concentramos mais nas várias técnicas de modelagem, juntamente com as soluções fornecidas para várias desvantagens. Uma discussão abrangente sobre modelos de condição de contorno que podem ser classificados como modelos pseudo-compactos também foi discutida. Algumas técnicas críticas de modelagem foram investigadas neste trabalho, as quais podem ajudar significativamente os desenvolvedores de modelos. Além disso, uma discussão sobre várias técnicas de simulação e plataformas para modelos RRAM, como SPICE [116, 117], foi incluída, o que é altamente essencial. Nosso trabalho visa preencher uma lacuna significativa na comunidade de modelagem RRAM.

Modelos RRAM para dispositivos bipolares

Modelo Chua


Leon O. Chua em 1971 apresentou a ideia do memristor [1] de que ele era de fato o quarto elemento básico ao lado do resistor, capacitor e indutor. Acredita-se que as características básicas de um memristor sejam controladas por fluxo ( φ ) ou com carga controlada ( q ) e são definidos por uma relação do tipo g ( φ, q ) =0.

Chua definiu a voltagem de um memristor como [1]:
$$ v (t) =M \ esquerda (q (t) i (t) \ direita) $$ (1)
Onde
$$ M (q) =d \ varphi (q) / dq $$ (2)
A corrente fluindo através de um memristor controlado por fluxo foi formulada como 1 :
$$ i (t) =W \ left (\ varphi (t) v (t) \ right) $$ (3)
Onde
$$ W \ left (\ varphi \ right) =dq \ left (\ varphi \ right) / d \ varphi $$ (4)
Aqui, os parâmetros M ( q ) e W ( φ ) são definidos como memristância incremental e memdutância incremental, respectivamente, devido a eles terem unidades semelhantes à resistência e condutância. O φ-q curvas para os três dispositivos de memristor são mostradas na Fig. 1. Essas curvas são geradas por um circuito básico de resistor de memristor (M-R) que dá origem a três tipos de memristores. O φ-q a variação para esses dispositivos é mostrada na Fig. 1a-e, respectivamente. A Figura 1b-f representa o I-V correspondente relações dos mesmos três memristores.

a - f Carga de fluxo ( ϕ - q ) curvas obtidas a partir de três memristores diferentes [1]

As equações apresentadas acima podem ser simplificadas no seguinte [1]:
$$ v =R (w) \ vezes i $$ (5) $$ \ frac {dw} {dt} =i $$ (6)
onde w é a variável de estado do dispositivo e R uma resistência generalizada que depende do estado interno do dispositivo.

O valor da memristância incremental (memdutância) em um instante de tempo t 0 depende do tempo de integração da corrente (tensão) de memristor completa de t =- t para t = t 0 . Então, isso se traduz no fato de que enquanto um memristor atua como um resistor normal a qualquer instante t 0 , mas seus valores de resistência (condutância) dependem do histórico completo da corrente do dispositivo (tensão), daí a justificativa do nome resistor de memória.

Curiosamente, no momento da voltagem especificada do memristor v ( t ) ou atual i ( t ), o memristor se comporta como um resistor linear variável no tempo. Mas no caso em que o φ-q curva é uma linha reta, ou seja, M ( q ) =R ou W ( φ ) =G , o memristor atua como um resistor invariante no tempo linear. Portanto, um dispositivo de memristor não pode ser usado na teoria de rede linear, mas pode ser usado para definir circuitos onde o estado presente dos parâmetros é dependente dos estados anteriores.

Mais tarde, em 1976, Chua e Kang [2] generalizaram o conceito de memristor para incluir sistemas memristivos que incluem muitos sistemas dinâmicos não lineares. Foi descrito pelas equações [2]:
$$ v =R \ left (w, i \ right) \ times i $$ (7) $$ \ frac {dw} {dt} =f \ left (w, i \ right) $$ (8)
onde w é definido como um conjunto de variáveis ​​de estado, R e f são funções explícitas do tempo. Uma diferença básica entre os memristores e os sistemas memristivos é que, no último, o fluxo não é mais definido exclusivamente pela carga. Os sistemas memristivos podem ser distinguidos de um sistema dinâmico geral em que não há corrente fluindo no dispositivo quando a queda de tensão nele é zero.

As equações de memristor foram usadas razoavelmente para definir o estado variável de uma chave de limite por Chua [1], que é a primeira instância do uso de memristores na modelagem de dispositivos. A formulação do memristor por Chua estabeleceu as bases para uma nova classe de dispositivos e aplicações variadas que usam um elemento de circuito básico para armazenar dados. Este conceito básico de memristores levou ao projeto de novas arquiteturas para futuras aplicações de memória não volátil, das quais RRAM é um candidato promissor. Tem havido uma quantidade significativa de teorias explicando o funcionamento de dispositivos RRAM e modelos que os definem, que são fundamentalmente baseados no modelo de memristor.

Uma aplicação muito interessante do modelo fluxo-carga é o seu uso [118] para definir um RRAM unipolar e implementá-lo no SPICE. Devido à simplicidade das equações de fluxo-carga, eles podem ser facilmente integrados em simuladores de circuito com poucas modificações. O modelo SPICE foi testado contra dados experimentais de HfO 2 dispositivo RRAM unipolar baseado em A relação não linear proposta para ajustar o obtido experimentalmente normalizou q - φ os valores são dados como [118]:
$$ q \ left (\ varphi \ right) ={q} _r \ times \ min \ left (1, {\ left (\ frac {\ varphi} {\ varphi_r} \ right)} ^ n \ right) $$ (9)
Aqui, φ r é o fluxo no ponto de RESET. Quando este valor q ( φ ) = q r for obtido, o CF desaparece e a corrente associada ao CF é redefinida para 0. Isso se traduz em que o dispositivo está no HRS. Para investigar a capacidade do modelo de reproduzir as características de chaveamento unipolar do dispositivo, uma operação de varredura de polarização padrão é realizada. A voltagem aplicada ao dispositivo no estado de reinicialização é aumentada progressivamente da polarização zero até atingir o LRS e, em seguida, a polarização é varrida de volta para zero volts. A corrente LRS é modelada usando uma forma modificada da relação de corrente do modelo Chua [1], dada como [118]:
$$ i (t) =\ left \ {\ begin {array} {c} K \ sqrt {\ varphi} v (t) \ kern0.75em \ mathrm {if} \ \ varphi <{\ varphi} _r \\ {} 0 \ kern4.25em \ mathrm {if} \ \ varphi ={\ varphi} _r \ end {array} \ right. $$ (10)
Presume-se que a corrente HRS seja controlada por uma emissão termiônica, então a corrente nesse estado é modelada como:
$$ i (v) ={I} _A \ left ({e} ^ {\ frac {v} {v_A}} - 1 \ right) $$ (11)
Os efeitos de limiar também são considerados no modelo. Foi assumido que o efeito da tensão limite surge devido aos efeitos de contato. Isso pode ser levado em consideração incluindo um limite de tensão para o cálculo do fluxo nos processos SET e RESET. A corrente modificada é dada por [118]:
$$ i (t) =\ left \ {\ begin {array} {c} {I} _A \ left ({e} ^ {\ frac {v} {v_A}} - 1 \ right) \ kern2.75em \ mathrm {if} \ \ varphi <{\ varphi} _s \\ {} K \ sqrt {\ varphi} v (t) \ kern3.75em \ mathrm {if} \ \ varphi <{\ varphi} _r \ end {array }\direito. $$ (12)
Aqui, ϕ r e ϕ s são o fluxo RESET e SET, respectivamente. Essas equações podem ser implementadas em um circuito compatível com SPICE composto por uma rede de capacitores. Os resultados da implementação do SPICE foram encontrados para seguir de perto os resultados experimentais com o modelo capaz de reproduzir características de memristor quase idênticas. Ele valida o uso do modelo de fluxo-carga de Chua [1] para ser usado também na modelagem de dispositivos unipolares.

Modelo Linear de Deriva de Íons


Com uma lacuna considerável nas décadas subsequentes após a formulação do memristor por Chua, pesquisadores do HP Labs [3] em 2008 fizeram uma descoberta empolgante com relação aos dispositivos de memristor. Embora Chua tenha formulado a presença de um elemento como um memristor, não houve um circuito ou modelo realizável desenvolvido depois disso, embora vários esforços tenham sido relatados para fabricar dispositivos RRAM no início do século XXI. A equipe do HP Labs liderada por Strukov et al. [3] perceberam um sistema memristivo em nanoescala funcional onde a memristância ocorre naturalmente, onde o transporte eletrônico e iônico de estado sólido são acoplados sob uma polarização de voltagem externa. Esses sistemas mostram uma relação histérica entre as características de corrente e voltagem semelhantes a outros dispositivos eletrônicos em nanoescala, levando assim a uma compreensão fundamental dos sistemas memristivos e do projeto de sistemas semelhantes.

Um dispositivo simples de dois terminais foi relatado, onde um óxido (TiO 2 ) de espessura D foi imprensado entre dois eletrodos de Pt. Histerese I - V curvas de comutação foram comparadas com a curva simulada. Embora o mecanismo exato desses dispositivos não fosse completamente compreendido naquela época, foi um dos primeiros casos em que as memórias de comutação resistiva foram classificadas em sistemas memristivos.

Uma estrutura esquemática de dispositivo de TiO 2 com base em memristor é mostrado na Fig. 2a [3], onde existem duas resistências variáveis ​​em série, chamadas de R LIGADO que é a baixa resistência na região do semicondutor com maior concentração de dopante. Uma concentração menor de dopante aumenta a resistência da outra parte, chamada de R DESLIGADO . Relação entre a tensão aplicada v ( t ) e atual através do sistema i ( t ) devido à condutância ôhmica eletrônica e deriva iônica linear em um campo uniforme com mobilidade iônica média é dada por [3]:
$$ v (t) =\ left (\ frac {R _ {\ mathrm {ON}} w (t)} {D} + {R} _ {\ mathrm {OFF}} \ left (1- \ frac {w (t)} {D} \ right) \ right) i (t) $$ (13)
O modelo de resistor variável acoplado para um memristor é apresentado. a Um circuito equivalente simplificado que compreende um voltímetro (V) e um amperímetro (A). b , c A tensão aplicada (azul) e a corrente resultante (verde) em função do tempo t para um memristor típico também são apresentados. Em b a tensão aplicada é v 0 pecado ( v 0 t ) e a relação de resistência é R DESLIGADO / R LIGADO =160, e em c a tensão aplicada é ± v 0 pecado 2 0 t ) e R DESLIGADO / R LIGADO =380, onde ω 0 é a frequência e v 0 é a magnitude da tensão aplicada. Os números de 1 a 6 são rotulados para ondas sucessivas na voltagem aplicada e os loops correspondentes em i – v curvas. Em cada gráfico, os eixos são adimensionais, com tensão, corrente, tempo, fluxo e carga expressos em unidades de v 0 =1 V, i 0 ≡ v 0 / R LIGADO =10 mA, t 0 ≡ 2 π / ω 0 D 2 / μ v v 0 =10 m / s, v 0 t 0 e i 0 t 0 , respectivamente. O termo i 0 denota a corrente máxima possível através do dispositivo, e t 0 é o menor tempo necessário para a deriva linear de dopantes em todo o comprimento do dispositivo em um campo uniforme v 0 / D , por exemplo com D =10 nm e μ V =10 −10 cm 2 s −1 V −1 . É de notar que para os parâmetros escolhidos, a polarização aplicada nunca força nenhuma das duas regiões resistivas a entrar em colapso; por exemplo, w / D não se aproxima de zero ou um (mostrado com linhas tracejadas nos gráficos do meio em b e c ) Além disso, o tracejado i – v trama em b demonstra o colapso da histerese observado com um aumento de dez vezes na frequência de varredura. As inserções de i – v plotagens em b e c mostre que, para esses exemplos, a carga é uma função de valor único do fluxo, como deve ser em um memristor [3]

Embora a equação acima seja não linear, a resistência do dispositivo muda linearmente com a tensão aplicada v ( t ), daí a atribuição de linearidade ao modelo. Dispositivo definido por Strukov et al. [3] atua como um memristor perfeito para apenas um determinado intervalo limitado da variável de estado w . A variável de estado é definida como [3]:
$$ \ frac {dw (t)} {dt} ={\ mu} _v \ frac {R _ {\ mathrm {ON}}} {D} i (t) $$ (14)
Memristância do sistema proposto por Chua [1] na Eq. (1) é definido usando as duas Eqs acima. (13) e (14) [3]:
$$ M (q) ={R} _ {\ mathrm {OFF}} \ left (1- \ frac {\ mu_v {R} _ {\ mathrm {ON}}} {D ^ 2} q (t) \ direita) $$ (15)
Na Eq. Acima (15), o q -termo dependente é a principal contribuição para a memristância. Uma análise interessante fornecida para saber por que esse fenômeno particular foi escondido por tanto tempo é devido ao campo magnético não desempenhar um papel explícito no mecanismo. Para que um memristor seja realizado em termos simples, deve existir uma relação não linear entre as integrais de tensão e corrente.

As Eqs. (13) - (15) também incorporam os fundamentos da comutação bipolar, ou seja, o dispositivo muda de um estado para outro pela aplicação de tensão de duas polaridades. Como resultado, dispositivos mostrando histerética bipolar I - V relacionamentos podem ser modelados por essas equações e, portanto, levar à classificação de tais dispositivos como sistemas memristivos. Tal comportamento é observado em muitos sistemas de materiais, como filmes orgânicos [119,120,121,122,123], calcogenetos [124,125,126], óxidos de metal [127,128,129], óxidos dielétricos [130,131,132], perovskitas [133,134,135,136], etc. A própria equipe da HP usou um TiO 2 [3] sistema e observaram características de comutação bipolar semelhantes, com o movimento de dopante ou impureza através da região ativa como a razão para tais mudanças dramáticas na resistência. Isso é mostrado na Fig. 2b, c com a corrente mostrando queda drástica e aumento rápido com a mudança na tensão.

Fisicamente, a região ativa nestes dois dispositivos terminais opera dentro do limite, 0 a D , a espessura da camada de óxido, então a variável de estado w também é limitado entre as espessuras. A Figura 3 indica a variação de w / D com tempo para o parâmetro nunca saindo dos limites de 0 e D [3]. A mudança repentina na resistência ou a comutação é causada pelos dispositivos que atingem esses limites. Para modelar esta condição, são utilizadas condições de contorno adequadas. Certas anomalias são observadas no dispositivo especificamente nas fronteiras. Há uma mudança não constante na taxa das variáveis ​​de estado dinâmico em relação à mudança disponível. Além disso, a mobilidade iônica é significativamente menor nas fronteiras do que no meio. Isso é atribuído aos efeitos de deriva de dopante não linear nos limites. Portanto, para contabilizar adequadamente esses efeitos, as variações de certas funções de janela são usadas para definir os limites dos dispositivos. A equipe da HP propôs uma função de janela multiplicada pela variável de estado Eq. (9) dado como [3]:
$$ f (x) =\ raisebox {1ex} {$ w \ left (1-w \ right) $} \! \ left / \! \ raisebox {-1ex} {$ {D} ^ 2 $} \ right . $$ (16)
Dispositivo memristivo simulado com voltagem. a Simulação com resistência diferencial negativa dinâmica. b Simulação sem resistência diferencial negativa dinâmica. c Simulação governada por deriva iônica não linear. Nos gráficos superiores de a , b , e c , o estímulo de voltagem (azul) e a mudança correspondente na variável de estado normalizada w / D (vermelho) é plotado contra o tempo. Em todos os casos, a mudança brusca ocorre quando w / D aproxima-se dos limites em zero e um (tracejado), e o qualitativamente diferente i - v formas de histerese são devidas à dependência específica de w / D no campo elétrico perto dos limites. d Para efeito de comparação, um experimental i – v gráfico de um Pt – TiO 2 - x –Pt dispositivo é apresentado [3]

Este modelo pode ser atribuído ao lançamento das bases para futuros modelos RRAM. Também pode ser usado para dois dispositivos semicondutores terminais com histerética bipolar I - V relacionamentos. Tomando o mecanismo de um memristor como referência, vários modelos futuros para dispositivos RRAM foram desenvolvidos.

Modelo de deriva de íons não linear


O modelo de deriva linear de íons desenvolvido pela HP [3] demonstrou principalmente os efeitos de deriva linear na região em massa do dispositivo de memristor. Eles observaram alguns efeitos não lineares nas fronteiras, mas não os definiram de forma abrangente. A dependência não linear da deriva do dopante na voltagem aplicada foi observada e formulada por Yang et al. [46] em 2008. Eles propuseram uma relação corrente-tensão responsável pelos efeitos não lineares com precisão. Posteriormente, foi melhorado e adicionado por Eero Lehtonen e Mika Laiho [68].

A condução em dispositivos memristivos é controlada por uma barreira eletrônica espacialmente heterogênea de metal / óxido, conforme relatado por Yang et al. [46]. A troca é causada pelo deslocamento de vacâncias de oxigênio carregadas positivamente agindo como dopantes nativos para formar ou dissolver canais condutores através desta barreira eletrônica. A concentração de vagas é maior nas fronteiras ou interfaces metal / óxido. A comutação ON e OFF ocorreu apenas na interface superior, o que indica que o eletrodo superior atua como o eletrodo ativo.

O efeito das vacâncias de oxigênio nas características de comutação do memristor à base de óxido de titânio é mostrado na Fig. 4 [46]. As amostras tendo diferentes vacâncias de oxigênio com diferentes sequências de camada de TiO 2 mostram comutação oposta definida por suas polaridades. Além disso, a adição de vagas extras à interface superior, mostrada na Fig. 4c, altera as curvas de comutação, confirmando assim o papel dominante das interfaces não ôhmicas em dispositivos memristivos. Isso forma a base dos efeitos de não linearidade que se originam nas interfaces e governam a comutação de dispositivos.

TiO de filme fino 2 - x dispositivos com perfis de vacância de oxigênio controlados são usados ​​para verificar o mecanismo de comutação. a As amostras I e II contêm sequências de camada reversa de TiO de 15 nm 2 e TiO de 15 nm 2 - x (mais vagas) camadas. Estes mostram polaridades opostas de I-V curvas em seus estados virgens. b . As polaridades de comutação dessas duas amostras também são opostas uma à outra. c . Quando mais vagas são introduzidas adicionando uma camada de Ti de 5 nm às interfaces superiores dessas duas amostras, o I-V as curvas mudam de maneiras totalmente diferentes, confirmando o papel dominante das interfaces não ôhmicas nos dispositivos de filme fino [46]

Yang et al. [46] explicaram o fato acima de que os dispositivos memristivos agem como resistores dinâmicos que mudam seu estado de acordo com a integral de tempo da corrente ou tensão aplicada; eles falharam em fornecer um relacionamento que descreve uma variável de estado dinâmico. A relação corrente-tensão proposta pode ser descrita como [46]:
$$ I ={w} ^ n \ beta \ sinh \ left (\ alpha v \ right) + \ chi \ left ({e} ^ {\ gamma v} -1 \ right) $$ (17)
Aqui, β, γ, n e χ são constantes de ajuste. Na equação acima, o primeiro termo β sinh ( αv ) aproxima [1] o estado ON do memristor, onde os elétrons tunelizam através da barreira eletrônica residual fina. w é definido como a variável de estado do dispositivo na faixa de 0 (OFF) e 1 (ON). A segunda parte da equação aproxima o estado OFF do dispositivo com os outros parâmetros atuando como constantes de ajuste. Parâmetro n aqui atua como o parâmetro livre usado para modificar a comutação entre os estados. Durante o ajuste de n, os efeitos não lineares entram em cena. eu - V curva do dispositivo fabricado é modelada usando a Eq. (16). O melhor ajuste é obtido em 14 ≤ n ≤ 22. Isso pode ser interpretado como evidência de que a velocidade efetiva do desvio da vacância depende de uma forma altamente não linear com a tensão aplicada ao dispositivo. Portanto, a maioria dos efeitos de deriva do dopante nas fronteiras / interfaces podem ser entendidos como de natureza não linear.

Um relacionamento que descreve a dinâmica da variável de estado w neste modelo usando SPICE [116, 117] foi proposto por Lehtonen e Laiho [68]. A derivada de tempo de w foi modelado como [68]:
$$ \ frac {dw} {dt} =a \ vezes f (w) \ vezes g (v) $$ (18)
Aqui, a é uma constante, f :[0, 1] → R é uma função janela proposta e g:R → R é considerada uma função linear proposta anteriormente no modelo de deriva linear (onde R representa números reais). Os autores demonstraram a partir das soluções que para imitar o funcionamento do memristor propostas por Yang et al. [46], g ( v ) deve ser uma função não linear, ímpar e monotonicamente crescente. Uma função não linear proposta foi [68]:
$$ g(v)={v}^q $$ (19)
Here, the exponent q is used to mimic the rapid switching process. Transition between ON and OFF state in a memristor generally takes place very fast. An input voltage with a very high sweep rate is used to obtain such behavior. This is the first implementation of memristor models in the SPICE platform [116, 117].The major advantage of SPICE implementation is the ability of the model to be used in analog circuits and simulations and can be verified as fit to be circuit implementable or not. Although many improvements were made in subsequent models, this model lays the foundation for the rest of the RRAM models by accurately taking into consideration and explaining the non-linear dopant drift effects [3, 46].

Exponential Ion Drift Model


In practice, resistance switching characteristics are non-linear in nature. To analyze such exponential characteristics, Strukov et al. [69] proposed exponential ion drift model in 2009. This non-linearity caused a significant variation in retention time and write speed. Due to the exponential dependence of the switching rate for high electric field, the exponential ion drift model is generalized to explain the phenomenon by the non-linear microscopic drift of charged species in the dielectric at high field and temperature.

The major factors considered for this model are switching speed and volatility. Switching speed is the time required for the device to switch from one resistance state to the other, i.e., it can be deemed as the time required to writing the data into the memory and is denoted as τwrite . Volatility is the time required for the device to lose its resistance state, i.e., the time taken to store the data into the device before erased denoted as τstore . The ratio between τstore and τwrite derived using the Einstein-Nernst formula is given by [69]:
$$ {\tau}_{\mathrm{store}}/{\tau}_{\mathrm{write}}\sim EL\mu /D=qEL/{k}_BT $$ (20)
Here, L is the length of the device with an active doped region D and k B the Boltzmann constant. Ratio between the two parameters is approximately three orders of magnitude when considered at room temperature and reasonable bias voltages. Such a high volatility to switching speed ratio suggests a strong non-linear ionic transport due to drift-diffusion inside the device. For high-field ionic drift, the overall effect on the average drift velocity of the ions is given by the model as [69]:
$$ \nu \approx {f}_e{a}_p{e}^{-\frac{E_a}{k_BT}}\sinh \left( qE{a}_p/2{k}_BT\right) $$$$ \nu =\left\{\begin{array}{c}-\mu E,\kern0.5em E\ll {E}_0\\ {}\mu {E}_0{e}^{E/{E}_0},\kern0.5em E\sim {E}_0\end{array}\right. $$ (21)
Here, ν is the drift velocity, f e the frequency of escape attempts, T the device temperature, a p the periodicity, E a the activation energy, and E the applied electric field.

Variation of the drift velocity with the applied electric field is shown in Fig. 5 [69]. The exponential variation can be clearly seen at high applied fields which lend non-linearity to the model. There are a few shortcomings for this model which affect its accuracy and also the calculation of the average drift velocity mentioned in Eq. (20). This model is primarily suited for application to ionic crystals where the major interaction forces are the Coulomb repulsion and van-der-Waals forces. Its application for covalent crystals will affect the accuracy of calculation due to the complex interactions of electrons and ions in high electric field. Also, electrochemical diffusion reactions and redox reactions are not explained by the model [91,92,93]. This can cause significant issues in the systems where the physical switching mechanism is governed by electrochemical processes.

Nonlinear (solid) and linear (dashed) drift velocity of doubly charged oxygen vacancies along the [110] plane direction in rutile structure at room temperature [69]

Simmons Tunneling Barrier Model


Though Lehtonen and Laiho [68] first proposed SPICE-based simulations model for non-linear ion drift model as mentioned in the “Non-linear Ion Drift Model” section, but this modeling is not suitable for use in an electrical-based time domain simulation, due to the lack of proper definition of simulation parameters and equations. This situation changed with the Pickett-Adballa et al. [70,71,72] model where a new class of model based on the device physics was demonstrated, which is capable of being explained and compatible with SPICE. The equations were modified to fit the requirements for SPICE implementation.

The analysis was based on the results from a TiO2 -based memristor device [70] where the tunneling barrier width w was considered to be the dynamic state variable. This later set the precedent for one of the most popular parameters being treated as the dynamic variable in memristor systems, the other being the length of conductive filament inside the dielectric media. The deduction based on their analysis was that the dynamic behavior for on and off switching of the devices was highly non-linear and asymmetric as can be seen in Fig. 6 [70]. The explanation provided for the deduction was the exponential dependence of the drift velocity of ionized dopants on the applied current or voltage.

Dynamical behavior of the tunnel barrier width w . The evolution of the state variable w occurs as a function of time for different applied voltages for a series of a off-switching and c on-switching state tests on the same device. Legends indicate the applied external voltage. The lines are the numerical solution to the respective switching differential equations described in the text. b , d The numerical derivative w ˙ of the data in a e c plotted as a function of w for the different applied voltages. The lines are calculated from the differential equations using the measured values of w e i at each point in time. The irregularity of the calculated vs w lines in the on-switching plots is caused by the changes in the current that accompany the change in state ( is a function of two variables, w e i , and both are changing). The derivative of the state variable can be interpreted as the speed of the oxygen vacancy front. This is because the applied voltage pushes it away from or attracts it toward the top electrode [70]

The current in the device was explained based on the Simmons tunneling barrier I-V expressions [137], and based on this analysis, the dynamic state variable was determined to be the Simmons tunnel barrier width (w ) The current was given as [72]:
$$ i=\frac{j_0A}{\Delta {w}^2}\left\{{\phi}_b{e}^{-B\sqrt{\phi_b}}-\left({\phi}_b+e\left|v\right|\right){\mathrm{e}}^{-B\sqrt{\phi_b+e\left|v\right|}}\right\} $$ (22)
Onde
$$ {j}_0=\frac{e}{2\pi h},{w}_1=\frac{1.2\lambda w}{\phi_0},\Delta w={w}_2-{w}_1 $$ (23) $$ {\phi}_I={\phi}_0-\left|{v}_g\right|\left(\frac{w_1+{w}_2}{w}\right)-\left(\frac{1.15\lambda w}{\Delta w}\right)\ln \left(\frac{w_2\left(w-{w}_1\right)}{w_1\left(w-{w}_2\right)}\right) $$ (24) $$ B=\frac{4\pi \Delta w\times {10}^{-9}\sqrt{2 me}}{h} $$ (25) $$ {w}_2={w}_1+w\left(1-\frac{9.2\lambda }{\left(3{\phi}_0+4\lambda -2|{v}_g|\right)}\right) $$ (26) $$ \lambda =\frac{e.\mathit{\ln}(2)}{8\pi \varepsilon {\varepsilon}_0w\times {10}^{-9}} $$ (27)
The parameters have been adjusted here such that the barrier height φ b is in volts (not in electron volts), and the time-varying tunnel barrier width w is in nanometers. In the equations above, A is the channel area of the memristor, e is the electron charge, h is the Planck’s constant, ε is the dielectric constant, m is the mass of electron, φ 0 is a standard barrier height taken from reference [70], and v is the voltage across the tunnel barrier. B is a fitting constant. In lieu of the analytical form of the equations, they can be conveniently described and implemented in SPICE, or it can be implemented with the any SPICE compatible electrical simulator.

The dynamic state variable w varies with time as [72]:
$$ \frac{dw}{dt}={f}_1\sinh \left(\left(\frac{\mid i\mid }{i_1}\right)\exp \Big(-\exp \left(\frac{w-{a}_1}{w_c}-\frac{\mid i\mid }{b}\right)-\frac{w}{w_c}\right) $$ (28)
This is in the case of off switching state (i > 0). Whereas for on switching state (i < 0), the state variable varies as [72]:
$$ \frac{dw}{dt}=-{f}_2\sinh \left(\left(\frac{\mid i\mid }{i_2}\right)\exp \Big(-\exp \left(\frac{a_2-w}{w_c}-\frac{\mid i\mid }{b}\right)-\frac{w}{w_c}\right) $$ (29)
Here, f 1, i 1 , a 1 , b , w c , f 2 , i 2 , and a 2 are fitting parameters. The abovementioned equations are used to model the memristor on the circuit level considering the electron tunnel barrier as a voltage-dependent current source, and the conducting channel (TiO2 ) is modeled as a series resistance. The voltage drops across the tunnel barrier and the series resistance make up the complete voltage drop across the circuit.

The dynamic behavior of the device is visibly complex as it is physics-based modeling approach and has been articulated as such by the Eqs. (27) and (28). The rate of switching possibly has contributions from the nonlinear drift at high electric fields and local Joule heating of the junction speeding up the thermally activated drift of oxygen vacancies [16, 46, 82, 83]. This can be clearly seen in the case of Fig. 6a, c [70] where the nature of the curves at high electric fields is quite different to those in low fields. The switching in the device is directly affected by the width of the gap. Application of a positive bias on the top electrode increases the state variable w resulting in an exponential increase in the resistance of the device as illustrated in Fig. 6b, d [70]. An opposite phenomenon occurs when negative bias is applied on the top electrode. This signifies the bipolar nature of the switching characteristics and their dependence on the dynamic state variable w .

The SPICE simulation of the model equations is illustrated in Fig. 7 [72]. The experimental data from the fabricated device is plotted against the simulated I-V curves showing a good fit between the two. This implementation paves the way for future SPICE simulations of RRAM devices [74, 77, 81]. A possible shortcoming in this model is the lack of a boundary for the dynamic variable and a threshold voltage within which the model should work. The growth of tunneling barrier width w can possibly go to unlimited quantities owing to the lack of a bound for the same, thus creating non-realizable scenarios for the device mechanism. Many models have employed what is called a window function to define the limits for the defined dynamic state variable in the model.

Experimental data (black dots) and corresponding simulated I - V curve for the memristor (solid line) where i mem is the current through the memristor and v mem is the voltage across the entire memristor. The inset shows the externally applied voltage sweep is shown and the initial condition for w is set at 1.2 nm [72]

Yakopcic Model


Although not validated specifically for RRAM devices at the time of development, the Yakopcic model [73, 74] closely resembled a variety of RRAM devices. The model was initially tested for TiO2 systems [73], and these systems are indeed one of the most popular ones along with HfO2-based RRAM devices.

This model was based on the Pickett-Adballa model [70,71,72] using a similar state variable, but it was modified to include neuromorphic systems as well. It was one of the first models to consider the functioning of synapses into their equations. This model was verified for the device used by the HP lab team to explain the working of memristive systems.

The state variable w ( t ), a value between zero and one considered here, directly affected the current through the device and also the dynamics of the device, i.e., the resistance. The current in the device is given as [73]:
$$ I(t)=\left\{\begin{array}{c}{a}_1w(t)\sinh \left( bv(t)\right),\kern2.25em v(t)\ge 0\\ {}{a}_2w(t)\sinh \left( bv(t)\right),\kern2.25em v(t)<0\end{array}\right. $$ (30)
Two functions, namely g (v ( t )) and f ( x ( t )), are responsible for the change in the state variable. a 1 , a 2 , and b are fitting constants. Change of the state of the variable is generally governed by a threshold voltage, i.e., there is a physical change in the device structure above a certain threshold voltage. The function g (v ( t )) here models the ON and OFF voltages of the device which also takes into account the polarity of the input voltage. This results in a better fit to the experimental data in case of bipolar switching where the values of set (v p ) and reset (v n ) voltage, i.e., the thresholds are different. It is defined as [73]:
$$ g\left(v(t)\right)=\left\{\begin{array}{c}{A}_p\left({e}^{v(t)}-{e}^{v_p}\right),\kern0.5em v(t)>{v}_p\\ {}-{A}_n\left({e}^{-v(t)}-{e}^{v_n}\right),\kern0.5em v(t)<-{v}_n\\ {}\kern2.75em 0,\kern3em -{v}_n\le v(t)\le {v}_p\end{array}\right. $$ (31)
A p and A n indicate the rate of the change of state once the voltage threshold is crossed. It can be understood as the dissolution or the rupture of the filament in terms of RRAM devices. There is in-built support for threshold values in the model, which enhances its applicability.

The state change variable modeled by the function f (w ( t )) is used to define the boundaries for the variable. It explains the motion of the charge carrying particles based on the threshold values, also adding the possibility to define the motion of the particles based on the polarity of the input voltage. This basically acts as a window function which restricts the state change variable within certain boundary given as [73]:
$$ f(w)=\left\{\begin{array}{c}{e}^{-{\alpha}_p\left(w-{w}_p\right)}{f}_p\left(w,{w}_p\right),\kern0.5em w\ge {w}_p\\ {}1,\kern10em w<{w}_p\end{array}\right. $$ (32) $$ f(w)=\left\{\begin{array}{c}{e}^{\alpha_n\left(w+{w}_n-1\right)}{f}_n\left(w,{w}_n\right),\kern0.5em w\le 1-{w}_n\\ {}1,\kern10.5em w>1-{w}_n\end{array}\right. $$ (33)
Here, f p (w ,w p ) is a window function which limits the value of f (w ) to 0 when x ( t ) = 1 and v ( t ) > 0. f n (w ,w n ) is a similar window function which does not allow the value of w ( t ) to become less than zero when the current flow is reversed.

The window functions are defined as [73]:
$$ {f}_p\left(w,{w}_p\right)=\frac{w_p-w}{1-{w}_p}+1 $$ (34) $$ {f}_n\left(w,{w}_n\right)=\frac{w}{1-{w}_n} $$ (35)
The movement of dynamic state variable, in simple words, the rate of switching, is governed by a differential equation. The growth and decay of the tunneling barrier width are the defining mechanism for this particular model, and it is given by [73]:
$$ \frac{dw}{\mathrm{d}t}=g\left(v(t)\right)f\left(w(t)\right) $$ (36)
Owing to the analytical nature of the coupled equations, they can be solved using a mathematical solver such as MATLAB [138, 139]. The differential equation can also be solved in MATLAB using the in-built solvers idt () and ddt () functions, which employ the time step integration method. This particular model was simulated using the characterization data of the TiO2 memristor from HP Labs [3], and the fitting obtained was pretty good when the fitting parameters are properly calibrated.

A separate SPICE implementation of the same model was reported by Yakopcic et al. [74] which were fitted and characterized for a multitude of devices for both sinusoidal and repeated sweep inputs. The SPICE implementation revealed a good accuracy and applicability of the model at the circuit level. The model was correlated with a variety of experimental data, and low error rates of about 6% were obtained. It was one of the first SPICE implementation where the model was tested under sinusoidal as well as repetitive sweeping inputs. This helps in determining the AC behavior of the device. Along with that, very important device variability analysis is performed which defines the error tolerance in the device. Variability is an important issue, when the RRAM device is used in large systems, such as arrays. The variability analysis performed is essential in knowing until which point the system can tolerate the variability. After reaching the critical point, there is possibility of errors in device read/write.

The model was also tested for read/write operations using 256 devices, which helps determine its usability in crossbar arrays. Similarly, it can be used for neuromorphic read/write operations to test the model applicability in that system. Device variability in the model is defined with change in the device parameters. So, changing the device parameters leads to a change in the simulated device I - V which is very useful in fitting the model with the experimental data. The values of the device parameters used can help define the accepted values of the particular parameters in the real case scenario. No convergence errors were found in the 256 array system, but with new RRAM array systems reaching higher density, applicability of the model there remains a question. Higher density array systems generally pose a convergence problem in SPICE simulations, but with proper parameter definition, it can be avoided. This model can be considered a new paradigm when it comes to circuit level SPICE simulations, variability analysis, and read/write operation simulations for RRAM devices.

TEAM/VTEAM Model


Threshold Adaptive Memristor (TEAM) model [75, 76] builds based on the Simmons Tunneling Barrier model [70,71,72] (discussed in the “Simmons Tunneling Barrier Model” section) and delivers a much simpler physics-based modeling approach for memristive systems. I - V relationship in this case is not fixed and can be chosen to fit any device which provides some amount of flexibility in the model. TEAM model arose from the need of simpler analytical equations which describe the mechanism of memristive systems accurately and which take less computation time.

This model is based on the approximation of the high non-linear dependence of the memristive device current; the device can be modeled as a device with threshold currents. The results are evident in Fig. 8. As with the tunneling barrier model, the internal state derivate is dependent on the current and the state variable itself, which is the effective tunnel width. It can be modeled effectively by [76]:
$$ \frac{dw(t)}{dt}=\left\{\begin{array}{c}{k}_{\mathrm{off}}\times {\left(\frac{i(t)}{i_{\mathrm{off}}}-1\right)}^{\alpha_{\mathrm{off}}}\times {f}_{\mathrm{off}}(w),\kern0.5em 0<{i}_{\mathrm{off}}
A sinusoidal input of 1 V applied to the TEAM model using the same fitting parameters as used in Fig. 10 [76]. The values of R ON e R OFF are set as 50 Ω and 1 kΩ, and an ideal rectangular window function is applied in Eqs. (38) and (37). a I - V curve and b state variable. It is to be noted that the device is asymmetric, i.e., switching OFF is slower than switching ON [76]

Variation of the state variable with time is asymmetrical in nature, as shown in Fig. 8b. This means that the ON and OFF switching times are not equal. In the Eq. (36), i em e i desligado act as the current thresholds. Functions f em and f desligado are window functions which bound the internal state variable x ( t ) within [w em , w desligado ] Window functions are described as [76]:
$$ {f}_{\mathrm{off}}(w)=\exp \left[-\exp \left(\frac{w-{a}_{\mathrm{off}}}{w_c}\right)\right], $$ (38) $$ {f}_{\mathrm{on}}(w)=\exp \left[-\exp \left(-\frac{w-{a}_{\mathrm{on}}}{w_c}\right)\right], $$ (39)
The window functions describe the dependence of the derivative in the state variable x . They work well within the described boundaries, but the problem arises when the device goes beyond the boundaries. There are no limiting parameters here, and the window function only describes the state variable inside a particular limit. If the device goes beyond the boundaries, it can cause convergence issues with the simulator and it does not make sense for good modeling practice in case of analog devices.

I - V relationship in this model is derived from the tunneling barrier model, as discussed in the “Simmons Tunneling Barrier Model” section. Due to the non-linear nature of the tunneling current, the change in resistance varies exponentially with the state variable. So, it is assumed that any change in the tunnel barrier width changes the memristance in an exponential manner which deduces to [76]:
$$ v(t)={R}_{\mathrm{ON}}{e}^{\left(\lambda /{w}_{\mathrm{off}}-{w}_{\mathrm{on}}\right)\left(w-{w}_{\mathrm{on}}\right)}\times i(t) $$ (40)
Here, λ is a fitting parameter and R ON the equivalent effective resistance at the bounds.

I - V relationship for this model can be seen in Fig. 8a [76]. Although there is a presence of a pinched hysteresis, the form and structure of the curve are not well-defined. The model is driven with a sinusoidal input of 1 V. The verification done for this model is different from the tunneling model [70,71,72] in terms of the platform used to simulate it. The latter model uses a SPICE macro model [72] to describe the equations, but SPICE takes up a significant amount of computation time. Modeling in Verilog-A [140,141,142,143] is much more efficient, and the TEAM model [75] utilizes this functionality to model the equations presented by them.

A slightly modified version of the TEAM model with the introduction of voltage threshold levels was reported by the same group, called Voltage Threshold Adaptive Memristor model (VTEAM) [77]. Discussed TEAM model was based on threshold currents, whereas VTEAM is based on threshold voltages. The major advantages cited for using threshold voltages is that comparison with current causes performance and reliability issues if the condition is not satisfied, i.e., a low-current threshold will automatically have a low-voltage threshold as well. This might affect the overall performance of the device. Also with a threshold voltage, there is no risk with going overboard with high power and voltage destroying the device as the values are automatically controlled.

The VTEAM follows a similar concept to the TEAM model, being based on an expression of the derivative of an internal state variable. The current is dependent on the state variable itself. The only difference is inclusion of a threshold voltage. The internal state variable (w ) is defined as [77]:
$$ \frac{dw}{dt}=\left\{\begin{array}{c}{k}_{\mathrm{off}}\times {\left(\frac{v(t)}{v_{\mathrm{off}}}-1\right)}^{\alpha_{\mathrm{off}}}\times {f}_{\mathrm{off}}(w),\kern0.5em 0<{v}_{\mathrm{off}} Similar to the TEAM model, the functions f em and f desligado act as window functions which bound the internal state variable w within [w em , w desligado ] As has been assumed in the model, current varies exponentially with the internal state variable on most occasions which is defined by [77]:
$$ i(t)=\frac{e^{-\frac{\lambda }{w_{\mathrm{off}}-{w}_{\mathrm{on}}}\times \left(w-{w}_{\mathrm{on}}\right)}}{R_{ON}}\times v(t) $$ (42)
The comparative analysis of the VTEAM model with the Yakopcic model [73, 74], BCM model [99] (discussed further in this article), and the TEAM model are presented in Fig. 9 [77]. It represents the flexibility that the model possesses, as it can be tuned to fit all the three models. It shows good agreement with all the three models illustrated, respectively, in Fig. 9a–c [77]. Fundamentally, the TEAM/VTEAM models are quite generalized physics-based models. This means that with the help of fitting parameters, they can be comparable with the multitude of other models, and fit to a variety of experimental characterization data from memristive systems.

The VTEAM model is compared with previously proposed memristor models [77]. a Yakopcic model [73]. b BCM model [99]. c TEAM model [76]

Stanford/ASU Model


A physics-based model which has become very popular is the one developed by Guan et al. and Chen et al. of Stanford University and ASU, known as Stanford/ASU model [78,79,80]. This model is exclusively developed for RRAM devices, rather than a generalized one for memristive systems which was fitted for those particular devices. It included the effect of critical phenomenon of switching such as Joule heating and temperature change, which had been neglected before. The developed model was applied in the I - V switching characteristics of HfO2 RRAM [144]. Along with it, Verilog-A [79] and SPICE [81] implementations of the model are also presented.

This model is based on the growth of conductive filament. The CF growth leaves a gap with the top electrode which is called as the filament gap. This growth of the filament gap is considered as internal state variable in this case. So, the rate of filament growth and the filament gap govern the dynamics of the model. The filament growth is explained due to the movement of oxygen ions and vacancy regeneration and recombination [145]. Considering the gap value g (nominally in the range of 0–3 nm) to be the state variable, the rate of change of g is defined as [78]:
$$ \frac{dg}{dt}={\nu}_0\exp \left(\frac{-{E}_{a,m}}{k_bT}\right)\sinh \left(\frac{q{a}_h\gamma v}{L{k}_bT}\right) $$ (43)
The parameter E a is the activation energy for vacancy generation and oxygen vacancy migration in the SET and RESET processes, respectively. v is the applied voltage across the device, ν 0 the velocity containing the attempt-to-escape frequency, L the switching material thickness and a h , the hopping site distance.

A significant feature of this model is the inclusion of variations in the model caused due to the stochastic property of the ion process and the spatial variation in the gap size among multiple filaments. To account for these variations in the model, a noise signal is added to the gap distance as [78]:
$$ g\mid t+\Delta t=F\left[g|t,\frac{dg}{dt}\right]+{\delta}_g\times \overset{\sim }{X}(n)\Delta t,\kern2.25em n=\left\lfloor \frac{t}{T_{GN}}\right\rfloor $$ (44)
The variation in the gap size δ g is defined as a function of the ions’ kinetic energy and invariably on the temperature in the filament and is given as [78]:
$$ {\delta}_g(T)=\frac{\delta_g^0}{\left\{1+\exp \left[\frac{\left({T}_{\mathrm{crit}}-T\right)}{T_{\mathrm{amb}}}\right]\right\}} $$ (45)
Here, T crit is defined as a threshold temperature beyond which there is a significant change in the gap size. This can be understood as the point where the device undergoes a physical transformation such as transitioning into a SET or RESET state. In this case, threshold is considered in terms of temperature, rather than voltage or current, whatever employed in the previous models [75,76,77]. So, the equation basically depicts the resistance fluctuation that occurs when the CF temperature is increased beyond the room temperature.

Now that temperature can be considered a critical driving force in the model, a modified form of the steady-state Fourier heat flow equation is implemented in this model. Rather than considering heat flow throughout the filament, the vicinity of the tip of the filament is considered. There is a dynamic inner domain temperature T which significantly changes with change in the cell characteristics, and an outer domain remains at an ambient room temperature T amb , related as [78]:
$$ {c}_p\frac{\partial T}{\partial t}=v(t)i(t)-k\left(T-{T}_{\mathrm{amb}}\right) $$ (46)
c p is the effective heat capacitance of the inner domain, and k the effective thermal conductivity are both fitted based on the type of oxide and electrodes used in the RRAM system. RESET transition from LRS to HRS generally has higher temperature associated with it across the device, while the SET transition has a considerably lower temperature. The current inside the device is modeled using a generalized conduction mechanism where the tunneling distance and field strength have an exponential relationship. This is true in case of tunneling current conduction mechanisms such as Poole-Frenkel, Fowler-Nordhiem, trap-assisted, or direct tunneling [9, 16, 46, 49, 51, 55]; these are the mechanisms most commonly associated with RRAM systems [51, 55, 61, 66]. The current conduction is defined as [78].
$$ i\left(g,v\right)={i}_0\exp \left(\frac{-g}{g_0}\right)\sinh \left(\frac{v}{v_0}\right) $$ (47)
The advantage with a generalized current equation is that for a particular device if some other mechanism is fitting better, it can be incorporated easily by adding the required parameters and adjusting their values accordingly. I - V response of the model compared with experimental data is shown in Fig. 10. The experimental response is shown in Fig. 10a while the simulated curve is shown in Fig. 10b. Simulated transient response shows the capabilities of the model in taking variations into account. Developed model was verified using Ngspice [146] as a macrocircuit. Ngspice is an open source SPICE simulator which is quite efficient and convenient for doing DC and AC analysis. This model can be implemented in MLC memory circuits and also to verify the efficiency of programming strategies and error correction codes [78].

a Experimental and b simulated transient responses of a HfO x RRAM device to the  2.3 V 50 ns input pulses. The experimental result is reported elsewhere [144] and included here in a for convenience. c In a larger time range, the simulated transient response for the same device including the gap size and temperature is shown. Current compliance set at 200 μ A in simulation [78]

A major feature of this model is implemented in the neuromorphic systems and RRAM synaptic device design [147]. This model has been tested against a HfO x /TiO x multi-stack RRAM system [148] which is implemented in a neuromorphic system. This gives the model great flexibility and wide applications as there are only a few models that are actually applicable for neuromorphic systems. Also, the model defined for these systems has been deemed tolerant to training error caused by device variation [149]. The gradual resistance modulation which is critical to the learning process in a synaptic device can be quantified in the model [150] which marks a significant development in using RRAM synaptic stacks in neuromorphic computing systems.

Physical Electro-Thermal Model


This model is an extensive physical model which describes the bipolar operation in RRAM devices using equations closely resembling the physical mechanisms. This model was reported by Kim et al. [87], and it was verified with a tantalum pentoxide (Ta2 O 5 )-based bi-layered RRAM structure [15, 151, 152]. It makes use of the finite element solving method employed in the previous model to solve the differential equations. The major value addition by this model over the model proposed by Larentis et al. [86] was the proper description provided for the SET state in the bipolar RRAM device. The previous model was inadequate in accommodating the complete transition and explaining it properly but this model makes up for that. Also, it improved upon a physical electro-thermal model reported by Menzel et al. [153] which attempts at calculating the CF temperature precisely.

It also uses the electro-thermal physics phenomenon approach for modeling which we have seen in the previous model [86]. The major advantage with models based on this concept is their ease of use owing to the simple fundamental equations and the flexibility to employ a proper finite element method (FEM) solver to simulate the system very accurately. But a major disadvantage is that the model becomes very difficult to implement in circuit solvers based on SPICE and providing an equivalent implementation in Verilog. This is because of the lack of support in SPICE and Verilog for properly defining partial differential equations which make up for the vastness of the model. Normal ordinary differential equations and the ones which are in analytical form can be solved in circuit solvers but partial differential equations (PDE) cannot be solved.

Electro-thermal models are equally important as compared to the other physics-based models discussed before because temperature is an important factor governing the set and reset processes. Ion and vacancy migration plays a dominant role for switching mechanism [16, 46], although the governing factors are behind this process and the exact type of ions is still up for debate. So, the fact that temperature is a governing factor in this process makes these models attention worthy. Also, experiments [85, 154] in this regard suggest that there is significant change in the temperature in the CF during the switching process. Some of the previous models discussed above have neglected this effect by considering conducting filament-oxide interface to be at room temperature or by taking constant conducting filament temperature [39, 86, 88, 89, 144].

The major difference between this model and the previously discussed electro-thermal model is in the expressions used to describe the drift-diffusion process. CF is described as a doped region where the oxygen vacancies act as dopants, and the CF runs from the top electrode to the bottom electrode. This is an assumption that many models take that the CF runs from one end of the electrode to the other when the state variable is considered as the length of CF. A few models discussed previously [78, 80] have used the filament gap to the top electrode as state variable. So, the assumptions generally vary from system to system and are dependent on what mechanism is employed to describe the device.

Another assumption taken to describe the drift-diffusion of vacancy migration is that the same equation used can describe both the oxygen ions and vacancies. This is generally the case to simplify the model and reduce the complexity of the equations. The rigid point ion model by Mott and Gurney [155] is employed here to describe the process given as [87].
$$ \frac{\partial {n}_D}{\partial t}=\nabla \times \left({D}_s\nabla {n}_D-\mu v{n}_D\right)+G $$ (48)
onde D s describes the diffusion process, v gives the drift velocity of the vacancies, and G is the generation rate of vacancy or the CF growth rate which actually describes the SET process. The G term is a specialized parameter added to better describe the complete switching process [156, 157]. The parameters are defined as [87]:
$$ {D}_s=\frac{1}{2}\times {a}^2\times {f}_e\times \exp \left(-{E}_a/{k}_{\boldsymbol{B}}T\right) $$ (49) $$ v={a}_h\times f\times \exp \left(-{E}_a/{k}_BT\right)\times \sinh \left(q{a}_hE/{k}_BT\right) $$ (50) $$ G=A\times \exp \left(-\left({E}_a-q{l}_mE\right)/{k}_BT\right) $$ (51)
Here, l m is the mesh size. So, using the Eqs. (48)–(50), the oxygen vacancy transport given in Eq. (47) can be defined which contains all the factors of drift-diffusion as well as the vacancy regeneration. These equations govern the CF growth and rupture which defines the physical transformation of the device during the SET and RESET transition of the device. So, it basically acts as a dynamic internal state variable which controls the switching rate of the device.

The simulation results for the reset transition is shown in Fig. 11 [87]. Concentration of the oxygen ions is shown at different voltages in Fig. 11a [87] which invariably governs the switching in the device. The point C (3.0 V) is the point where the reset transition occurs, so the concentration of ions is also the highest at the interfaces for that voltage point as evident in Fig. 11b [87]. On similar lines, the temperature and flux are on the higher side which can be seen in Fig. 11c, d, respectively [87].

Simulation results for the reset transition of the device. a V o density (n D ) map. Calculated profiles of b n D , c T , e d y for states A (1.0 V), B (1.7 V), and C (3.0 V). The position of z  = 15 nm indicates the Ta2 O 5 /TaO x interface in the structure schematic. The shaded area shows the depleted gap, defined for n D  < 5 × 10 21 cm -3 [87]

Equations (95) and (98)mentioned further are also used in the model to describe the current conduction and the temperature change due to Joule heating in the device. The equations are simultaneously solved in COMSOL to generate the required simulated profiles. The obtained simulated profiles are compared and verified against a TaOx bi-layered RRAM system [87]. In addition to the DC I - V characteristics the model was also used to generate time-dependent reset characteristics by investigating its response to square pulses.

Huang’s Physical Model


A very comprehensive physical model of RRAM devices is developed by Huang et al. [88, 89]. Its major feature is its consideration of the multitude of factors affecting the CF dynamics in the RRAM device. This model is comprehensive in the sense that it considers both the width of CF as well as filament gap to the electrode as factors affecting the state variable dynamics. The model was validated in a TiO2 based device and also applied in a 2 × 2 RRAM array cell [88].

Covering bipolar devices primarily, it also accounts for the temperature distribution in the device with multiple heating sources. SET/RESET process is considered to be caused due to generation/recombination process of the oxygen ions (O 2− ) and oxygen vacancies (V o ) Top electrode (TE) is the active electrode and acts as an oxygen reservoir for the release or absorption of oxygen ions [88]. The CF evolution during the SET process is modeled based on the width of the CF. Growth of the CF is thought to start from the tip of the active electrode. With an increase of voltage the CF enlarges along the radius resulting in a final width of the CF as w . So, the value of w is critical to determine the LRS resistance in the SET process. Huang et al. [88] assumed that the CF grows in a symmetrical cylindrical shape which is simplifying at best. While the cylinder has been the most popular to describe the shape of the CF, it might not be the most accurate.

Rupture of the CF during the reset process is considered to start from the TE first. CF disconnects from the starting point and then dissolves internally with increase in the voltage. Distance between the tip of the CF and the active electrode layer is defined as the filament gap distance (x ) The value of x determines the resistance of HRS during the RESET process. x and dx/dt are thus critical in defining the RESET process. A very important feature of the model is that there are two parameters defining the state of the system, in place of one parameter. The parameter w acts as the state variable for the SET process and x for the RESET process. So dx/dt and dw/dt define the dynamics of the device during the SET/RESET transition. Analytical model for a RRAM cell presented by Huang et al. [88] is developed by modeling the parameters x, w and their evolving speeds.

This model also presents one of the most detailed descriptions for the processes involved behind the RESET process. The rate of the CF shortening is affected by three processes, (a) O 2− release by the electrode, (b) O 2− hopping in the oxide layer, and (c) recombination between O 2− e V o . Slowest process among the three dominates the CF reduction process which is defined by the parameter x . Speed of the processes is affected by the specific device characteristics and the oxide used.

CF reduction rate during first reset process, i.e., O 2− release by the electrode can be given as [89]:
$$ \frac{dx}{dt}=a\times f\times \exp \left(-\frac{E_i-\gamma ZeV}{k_BT}\right) $$ (52)
In case of the O 2− hopping in the oxide layer, the CF with a being the distance between two Vo , reduction rate is described by [89]:
$$ \frac{dx}{dt}=a\times f\times \exp \left(-\frac{E_h}{k_BT}\right)\sinh \left(\frac{a_h ZeE}{k_BT}\right) $$ (53)
The RESET process when dominated by the recombination between O 2− and Vo is written as [89]:
$$ \frac{dx}{dt}=a\times f\times \exp \left(-\frac{\Delta {E}_r}{k_BT}\right) $$ (54)
The value of x is fixed to x 0 after the RESET process. This invariably will act as the boundary condition for the model. But the problem here is the value and the role of x 0 is not clearly defined here. This will possibly create ambiguities while defining the states of the device or switching between two states. In the first step of the SET process which is dominated by recombination of oxygen vacancies and where a thin CF is initially grown is described by [89]:
$$ \frac{dx}{dt}=-a\times {f}_e\times \exp \left(-\frac{E_a-{\alpha}_a ZeE}{k_BT}\right) $$ (55)
Here, Z and αa are fitting parameters. In the second step, the CF grows along the radial direction of the CF is defined as [89]:
$$ \frac{dw}{dt}=\left(\Delta w+\frac{\Delta {w}^2}{2w}\right)\times {f}_e\times \exp \left(-\frac{E_a-\gamma Zev}{k_BT}\right) $$ (56)
Current flowing through the device has been taken in the model due to the hopping conduction and metallic conduction. The current in CF region can be calculated using the basic structures of Ohm’s law and Arrhenius law [158]. But the current in the gap region as a result of hopping conduction is given a little different. It is modeled as a correlation of the hopping current with the voltage and gap distance is given by [147]:
$$ i={i}_0\exp \left(-x/{x}_T\right)\sinh \left(v/{v}_T\right) $$ (57)
Temperature effects in the model are considered from the Filament Dissolution model [82, 83] discussed further in the “Filament Dissolution Model” section. Validation of the model is performed in HfOx /TiOx system [88, 89]. Transient results obtained from simulating the model are compared against the data from the device, which shows a good match as demonstrated by Huang et al. [88]. The model is also validated against devices fabricated by other groups [144, 159] and the parameters are adjusted accordingly. A pretty accurate match between the simulation and the experimental results suggests a good level of flexibility with the model. The model also demonstrates that the switching speed of the device is highly dependent on the input voltage sweep rate.

Although the model is very comprehensive and takes into account a variety of detailed processes affecting the RRAM operation; it has some critical shortcomings. A major one is the non-compatibility with SPICE or Verilog-A. Implementations in any of the circuit simulators based on these platforms has not been demonstrated which raises a question on its readiness for simulations. Also, boundary conditions and non-linear effects have not been applied in the model which leaves it open to unphysical solutions. There has been no attempt to fit a window function with the model to account for this effect. These shortcomings make the model difficult for application for simulations, but its physics give a lot of insights into the functioning of RRAM devices.

Bocquet Bipolar Model


A very interesting and unique model from Bocquet et al. [90, 92] which utilizes a physics based modeling approach to describe bipolar oxide based resistive switching memories. This was a model developed exclusively for the RRAM devices. Although a point of speculation still exists, it has been more or less accepted that the bipolar resistive switching mechanism is governed by the valence change mechanism which occurs in specific transition metal oxides and the field-assisted motion of oxygen ions O 2− [160].

This is also one of the few models that can describe electroforming process. This process basically initiates the CF growth for the first time when the device is in a pristine state. It requires significantly higher voltage as compared to the set or reset voltage because the CF formation requires an electric breakdown of the oxide and this requires higher voltage and energy. However, forming free RRAM devices have been reported [85] by adjusting the oxygen stoichiometry of the active layer. Removal of the forming process will reduce the voltage requirement of the device and make it more energy efficient.

Bocquet bipolar model uses some concepts from the Bocquet unipolar model [90] and modifies it significantly according to the bipolar switching characteristics. Major features of the model are its intrinsic simplicity in the model equations, full compatibility with SPICE based electric simulators and inclusion of voltage and time dependencies of the device. Internal state variable here is the radius of the CF which governs the switching rate. Radius of the CF varies with growth/rupture mechanism of the CF which is explained in the model with the help of local electrochemical redox processes [82, 83, 105, 161] which are dependent on the applied bias polarity. A single master equation in which both the SET and RESET processes are accounted for simultaneously is controlled by the CF radius which thus gives the switching rate of the device.

Electroforming stage is modeled using electroforming rate which describes the process of conversion of the pristine oxide into a switchable sub-oxide layer. CF radius (r CF ) varies from a minimum value of 0 to a maximum value of r CFmax . The electroforming stage is modeled as [92]:
$$ {\tau}_{\mathrm{form}}={\tau}_{\mathrm{form}0}\times {e}^{\frac{E_{a\mathrm{Form}}-q\times {\alpha}_s\times {v}_{\mathrm{Cell}}}{k_b\times T}} $$ (58) $$ \frac{d{r}_{\mathrm{CFmax}}}{dx}=\frac{r_{\mathrm{work}}-{r}_{\mathrm{CFmax}}}{\tau_{\mathrm{form}}} $$ (59)
Some of the simplifying assumptions in the model are regarding the current conduction in the LRS and HRS. During the LRS, the conduction is assumed to be Ohmic, i.e., it follows Ohm’s law. In the HRS region, the current is dominated by a leakage current in the sub-oxide region which is basically due to trap-assisted conduction, but for simplicity sake, Ohmic conduction is considered here. The SET/RESET operation in the model is described by the electrochemical redox reaction derived from the Butler-Volmer equation [162] given as [92]:
$$ {\tau}_{\mathrm{Red}}={\tau}_{\mathrm{Red}\mathrm{ox}}\times {e}^{\frac{E_a-q\times {\alpha}_s\times {V}_{\mathrm{cell}}}{k_b\times T}} $$ (60) $$ {\tau}_{Ox}={\tau}_{Redox}\times {e}^{\frac{E_a+q\times \left(1-{\alpha}_s\right)\times {V}_{\mathrm{Cell}}}{k_b\times T}} $$ (61)
Here, τRed and τOx are the reduction and oxidation reaction rates, respectively. τRedox is the effective reaction rate considering both the reduction and oxidation reactions. Above two equations are coupled together in a master equation which define the switching rate given as [92]:
$$ \frac{d{r}_{CF}}{dt}=\frac{r_{CF max}-{r}_{CF}}{\tau_{red}}-\frac{r_{CF}}{\tau_{Ox}} $$ (62)
This is quite a comprehensive model in the sense that it includes the temperature effects as well. Temperature plays a significant role in the redox reaction rates [163, 164] and thus the local temperature in the filament is a very important parameter in this regard. The basic heat equation is used in this model and modified it accordingly given as [92]:
$$ \sigma (x)\times E{(x)}^2=-k\times \frac{\partial^2T(x)}{\partial {x}^2} $$ (63) $$ T(x)={T}_{\mathrm{amb}}+\frac{v_{\mathrm{Cell}}^2}{2\times {L}_x^2\times k}\times \left(\frac{L_x^2}{4}-{x}^2\right)\times {\sigma}_{\mathrm{eq}} $$ (64) $$ T={T}_{\mathrm{amb}}+\frac{v_{\mathrm{Cell}}^2}{8\times k}\times {\sigma}_{\mathrm{eq}} $$ (65) $$ {\sigma}_{\mathrm{eq}}={\sigma}_{\mathrm{CF}}\times \frac{r_{\mathrm{CF}}^2}{r_{\mathrm{work}}^2}-{\sigma}_{Ox}\times \frac{r_{\mathrm{CF}\mathrm{max}}^2-{r}_{\mathrm{CF}}^2}{r_{\mathrm{work}}^2} $$ (66)
On the face of it, the equations seem pretty complex to evaluate. But in reality, they are analytical in nature which makes them easily solvable in a numeric solver and can be implemented in an electric simulator. This is a major advantage of this model. Almost all of the models which employ the concept of temperature change in the filament follow the basic principles of the filament dissolution model [82, 83] discussed further in the “Filament Dissolution Model” section. During set operation, the temperature rises due to the increase in the CF radius, while it falls due to a decrease in the CF radius during the reset operation. This creates a positive feedback loop between the two processes leading to a self-accelerated reaction. This forms the basis of the filament dissolution model and all models incorporating the temperature effects in the device converge on this phenomenon [82, 83, 86,87,88,89, 92].

I - V characteristics of NiO based RRAM along with simulated curve using Bocquet model is presented in Fig. 12 [91]. Figure 12a represents the set and reset transitions of the device while Fig. 12b highlights the forming process. The current conduction in the Bocquet bipolar model is treated a little differently from what we have seen from previous models [87, 88, 90]. It considers the current as a combination of contributions from three different sources. The first one is the current from the conductive area (i CF ), the second is the conduction through the switchable sub-oxide (i sub-oxide ) and then the conduction through the pristine device (i pristine ) The total current is described as [92]:
$$ {i}_{\mathrm{cell}}={i}_{\mathrm{sub}-\mathrm{oxide}}+{i}_{\mathrm{CF}}+{i}_{\mathrm{pristine}} $$ (67) $$ {i}_{\mathrm{CF}}=E\times \pi \times {\sigma}_{\mathrm{CF}}\times {r}_{\mathrm{CF}}^2 $$ (68) $$ {i}_{\mathrm{sub}-\mathrm{oxide}}=E\times \pi \times {\sigma}_{ox}\times \left({r}_{\mathrm{CF}\mathrm{max}}^2-{r}_{\mathrm{CF}}^2\right) $$ (69) $$ {i}_{\mathrm{pristine}}={S}_{\mathrm{cell}}\times {A}_e\times {E}^2\times \exp \frac{-{B}_e}{E} $$ (70) $$ {A}_e=\frac{m_e\times {q}^3}{8\pi \times h\times {m}_e^{ox}\times {\phi}_b} $$ (71)
a , b Experimental I (V) characteristics for Electroforming , Set , and Reset processes measured on a large number of memory elements to understand the device-to-device variability. The experimental device-to-device variability is accounted for in Monte Carlo simulations with a ± 5% standard deviation on parameters α and L x [21]

The parameter Be is metal-oxide barrier height (ϕ b ), dependent which is given as [92]:
$$ \mathrm{if}\ {\phi}_b\ge q{L}_xE:{B}_e=\frac{8\pi \sqrt{2{m}_e^{ox}}}{3\times h\times q}\left[{\phi}_b^{\frac{3}{2}}-{\left({\phi}_b-q{L}_xE\right)}^{\frac{3}{2}}\right] $$$$ \mathrm{otherwise}:{B}_e=\frac{8\pi \sqrt{2{m}_e^{ox}}}{3\times h\times q}\times {\phi}_b^{3/2} $$ (72)
Here, Scell is a section of the RRAM cell. A e e B e are additional parameters defined to make the equations concise. To implement the model in an electrical simulator, discrete solutions are required which are well provided in this model. This makes the model suitable for proper simulation involving electrical circuits, thus widening its use case scenario. In this model, the equations are implemented in an Eldo circuit simulator [165, 166]. Memory effect of the device was replicated in the form that for each cell of the RRAM instance during transient simulation, the previous state of the filament as well as the applied voltage are given as the present state of the device [92]. New state gets solved as a function of these new inputs and the time step considers in the transient simulation. The discrete solutions are given as [92]:
$$ {r}_{{\mathrm{CFmax}}_{i+1}}=\left({r}_{{\mathrm{CFmax}}_i}-{r}_{\mathrm{work}}\right)\times {e}^{\frac{-\Delta t}{\tau_{\mathrm{Form}}}}+{r}_{\mathrm{work}} $$ (73) $$ {r}_{{\mathrm{CF}}_{i+1}}=\left({r}_{{\mathrm{CF}}_i}-{r}_{{\mathrm{CF}\mathrm{max}}_i}\times \frac{\tau_{\mathrm{eq}}}{\tau_{\mathrm{Red}}}\right)\times {e}^{\frac{-\Delta t}{\tau_{\mathrm{eq}}}}+{r}_{{\mathrm{CF}\mathrm{max}}_i}\times \frac{\tau_{\mathrm{eq}}}{\tau_{\mathrm{Red}}} $$ (74) $$ \mathrm{where}\ {\tau}_{\mathrm{eq}}=\frac{\tau_{\mathrm{Red}}\times {\tau}_{\mathrm{Ox}}}{\tau_{\mathrm{Red}}+{\tau}_{\mathrm{Ox}}} $$ (75)
The model has been verified against electrical characterization from an HfO2 based system [167]. To better judge the model on circuit level, a 2T/1R bipolar OxRAM [168] cell was simulated using Eldo, as shown in Fig. 13a [92]. Simulation of this type helps check the stability of the model when applied to a system environment. Current variation, voltage, and CF radius (r CF ) are shown with respect to time. Voltage follows a triangular wave form, which is the input sweep. Current in the device transitions from high to low and vice-versa depend on the voltage levels. The sudden drop in the current levels, as shown in Fig. 13b [92], indicates the device transition. CF radius follows a similar path as the current which is expected behavior of the internal state variable.

Simulating the electrical characteristics of the considered 2T-1R OxRAM structure [174]

Another important feature of OxRAM that has been highlighted in the model is the soft-reset [168]. It mainly induces the dependence between resistance in HRS and the stop voltage during the preceding reset operation. This phenomenon is basically due to the incomplete destruction of the CF during the reset process. So, the CF radius and temperature decrease during this process, leading to a decrease in the reaction rate. This means a self-limited reaction rate thus getting the name soft-Reset . This model can account for the device to device variability very efficiently [169, 170]. The standard deviation obtained for the important parameters such as the length of the oxide (Lx ) is well within the accepted range, thus accounting for the variations when the materials change in different devices [167, 168].

A shortcoming in the model which can be highlighted is the lack of a voltage or current threshold. Also, it works on the simplifying assumption that the CF radius grows from one end of the top electrode to the other end of the bottom electrode. This makes the model immune to significant fluctuations if the growth of the CF is not complete, thus leaving a filament gap. There is no provision to account for the effect of the filament gap if it occurs.

Berco-Tseng Model


The proposed model and simulation approach [171,172,173,174,175] by Berco-Tseng for RRAM devices is based on describing the CF growth process. The Gibbs free energy criteria [174, 175] is used as an indicator to model the growth dynamics of the CF. Simulation approach for the forming, set and reset process in the model is based on the Metropolis Monte Carlo algorithm [174]. This approach importantly does not rely on time evolution of the CF, thus increasing the efficiency of comparison of the relative retention properties of MIM structures.

The model is quite comprehensive in terms of describing the underlying physical parameters which affect the CF kinetics in the resistive switching layer. It also introduces the concept of “hot-spots” [172,173,174] which are random localized initial clustering of oxygen vacancies which facilitate the formation of the CF. The major parameter governing the Gibbs free energy is the enthalpy of formation of an oxygen vacancy [174] is used to define the CF growth dynamics in the switching layer and integrate it into the Monte Carlo simulator. As a result, all the CF processes, namely forming, set and reset can be effectively simulated.

A monoclinic HfO2 switching layer is primarily used to implement the model. But other material stacks, such as ZrO2 , Cu-HfO2 are also studied and compared with by us. Electrical conduction in the RRAM device is modeled on the Arrhenius relation, given as [174]:
$$ {\sigma}_n={\sigma}_0\exp \left(-\frac{E_{ac}}{k_B{T}_n}\right) $$ (76)
Typical boundary condition, such as V dd at top electrode and ground at bottom electrode is applied to the device. For modeling the CF accurately, it is divided into a grid structure to discretize it, which is in line with the finite element analysis (FEA) method. The various parameters defining each grid site are its spatial coordinates (x , y ), local potential φ , temperature T , N o , N ov , trap occupancy c , electrical conductivity σ and thermal conductivity k t h . The various processes associated with the evolution of CF within the oxide layer involves generation, recombination and hopping of oxygen (O ), oxygen vacancies (OV ) and electrons.

As a result, these processes are defined in terms of probabilities in the MC simulator. The probabilities are defines according to the minimum energy criteria as discussed earlier. This approach smoothens out the iterative steps in the simulator rather than the abrupt + 1 or − 1 levels in the discrete approach. The generation (P g ), hopping (P h ), and recombination (P r ) probabilities for the oxygen species are given as [174]:
$$ {P}_{g,n}=\left(1-{C}_n\right)\exp \left(-\frac{E_a-{a}_a ZeE}{k_B{T}_n}\right) $$ (77) $$ {P}_{h,n\to m}={C}_n\left(1-{C}_m\right)\exp \left(-\frac{E_h-{a}_h ZeE}{k_B{T}_n}\right) $$ (78) $$ {P}_r={C}_n\exp \left(-\frac{\Delta {E}_r}{k_B{T}_n}\right) $$ (79)
Aqui, E is the local electric field, C n represents the ratio of N ov (density of oxygen vacancies) in the low state to the maximal one at site n (n th grid site). E a e E h are the activation energies for oxygen species generation and hopping respectively. Similarly, a a e a h are the field lowering factor for O generation and hopping.

The mechanism of the current conduction is considered to be trap-assisted tunneling [171, 173, 174] through the switching layer. The Mott variable hopping model [174, 175] is taken into account here to model the tunneling effect. Mott DC conduction considered at any two grid points m, n are given as [174]:
$$ {\sigma}_{mn}={f}_e\frac{e^2{c}_m\left(1-{c}_n\right)}{d_{mn}{k}_B{T}_n}\exp \left(-\frac{2{d}_{mn}}{\alpha}\right)\exp \left(-\frac{e\mid \nabla {\varphi}_{mn}\mid }{k_B{T}_n}\right) $$ (80)
Here, d mn is the distance between m and n , α is the typical attenuation length of the electron wave function in the trap and c being the trap occupancy.

The model is simulated and the results are compared with many different material stacks of RRAM devices from other groups as well. The Fig. 14 shows the simulated heat map of the N o (density of oxygen) concentration during the CF rupture process. The simulation result corroborates with the experimental results, showing wide applicability potential for the model.

Resistive switching layer and Ti/HfO2 interface oxygen spatial density (No) plot during the CF rupture using negative bias:a initial and b final [174]

Gonzalez-Cordero et al. Bipolar Model


It is a compact physical model proposed by Gonzalez-Cordero et al. [93] describing the working of bipolar RRAM systems. The model is unique because it considers the CF as a truncated cone, which is a significant departure from previous models considering the CF shape generally as a cylinder. Also, the model is validated by implementing it in Verilog-A which gives us a closer look into the description and simulation of RRAM devices on the circuit level using Verilog-A. The proposed model builds upon the concepts introduced in the previous Bocquet bipolar model [91, 92] and modifies it accordingly to suit the new CF shape proposed.

One of the important aspects about the model is the consideration of a truncated cone shaped CF [176,177,178,179]. Majority of the models we have encountered till now consider the CF as a symmetrical cylinder which is more of a simplifying assumption [91, 92]. This is because it has been shown that the CF can grow [39, 51, 55] from one end of either electrode to the other depending on the active electrode. So, it is quite possible that the CF in this case might not be a perfect cylinder. So, a truncated cone is equipped to account for any variability and fluctuations arising due to the shape of the CF. Shapes other than simplified geometrical shapes are not considered in the models because of algebraic complexities. In previous models [91, 92], we noticed that device to device and cycle to cycle variability’s have a significant effect on the application of particular models to devices. So, by taking a truncated cone as the CF shape provides this model more flexibility than the others.

Another significant feature is the role of temperature in the CF and the reset process. Majority of the models which describe the CF rupture due to the self-accelerated dissolution, consider that the process takes place at the CF narrowing point and temperature increases at that particular point [82,83,84, 91, 92]. This point is generally in the middle of the cylindrical CF due to its symmetry. So, when we look at it from a physical standing point, the temperature at each of the points in all the RRAMs stacked together in a circuit has to be evaluated. Realizing this from the circuital standing point and simulating thousands of devices in the circuits is a very time-consuming process and slows down the simulation. This problem can be circumvented by considering two temperatures in the CF instead of the general single temperature; this approach also keeps the simplicity of the model intact. Two temperatures represent the main CF body that is not destroyed during the reset operation and the CF narrowing. This has been implemented in this model by considering the two temperatures as the wide region and the narrow region of the truncated cone, respectively.

This model extends the previously discussed Bocquet bipolar model [91, 92] in the “Bocquet Bipolar Model” section. In the previous model, the equations were defined keeping a cylindrical CF in mind, so the equations here have been modified to account for the change in the CF shape. The truncated cone CF is described by two different radii. CF is considered to grow from the top electrode to the bottom; the interface radius with the top electrode (TE) is r CFT which is always greater than the radius of the interface with the bottom electrode (BE) r CFB . This adheres to the structure of the truncated cone. An assumption is made here that during CF rupture, height of the cone is not affected; this makes the model open to fluctuations if there is any filament gap produced due to premature growth of CF. Although a forming process is considered for the device, it is not included in the model making the model not suitable for application to devices where forming is a significant factor. A possible explanation for leaving out the forming process is to avoid adding more complexity to the model because the forming parameters have to be included in the set/reset equations as well.

Similar to the previous model, the set/reset processes are described by an electrochemical redox reaction and diffusion processes [82,83,84, 105, 161, 162] which control the growth and rupture of the CF, respectively. The reduction and oxidation reaction rates are given as [93]:
$$ \frac{1}{\tau_{{\mathrm{Red}}_{\left(T,B\right)}}}={A}_{\mathrm{Red}\mathrm{ox}}\times {e}^{-\frac{E_a-q\times {\alpha}_s\times {v}_{\left(T,B\right)}}{k_B\times {T}_{{\mathrm{CF}}_{\left(T,B\right)}}}} $$ (81) $$ \frac{1}{\tau_{{\mathrm{Ox}}_{\left(T,B\right)}}}={A}_{\mathrm{Redox}}\times {e}^{-\frac{E_a-q\times \left(1-{\alpha}_s\right)\times {v}_{\left(T,B\right)}}{k_b\times {T}_{{\mathrm{CF}}_{\left(T,B\right)}}}} $$ (82)
Velocity of the CF radius increase/decrease which is basically the switching rate is controlled by two master differential equations which include the top (r CF,T ) and bottom (r CF,B ) radius of the CF. They are given as [93]:
$$ \frac{d{r}_{CF_T}}{dt}=\frac{r_{CF_{TM}}-{r}_{CF_T}}{\tau_{Red_T}}-\frac{r_{CF_T}}{\tau_{Ox_T}} $$ (83) $$ \frac{d{r}_{CF_B}}{dt}=\frac{r_{CF_T}-{r}_{CF_B}}{\tau_{Red_B}}-\frac{r_{CF_B}}{\tau_{Ox_B}} $$ (84)
CF radius is set a boundary which defines the limits for the CF growth/rupture which is given as [93]:
$$ {r}_{CF_{Tm}}\le {r}_{CF_T}\le {r}_{CF_{TM}} $$ (85) $$ {r}_{CF_B}\le {r}_{CF_{BM}} $$ (86)
Here rCF,TM /rCF,BM is the maximum top/bottom radius that can be achieved and rCF,Tm the minimum value of the top radius. This equation indicates CF geometry following a truncated cone structure and the top radius is greater than the bottom radius at all times. The model employs a numeric solving method similar to the one used in the previous model to find the discrete solutions for the master differential equations. But the solution for this model is not found, which means it is difficult to validate the reliability of the equations.

A very interesting point here is that a separate local diffusion process is not added in this model to describe the reset process in addition to the oxidation/reduction process. Many of the previously discussed models [91, 92] have a separate equation for the diffusion. But in this model diffusion has been integrated into the Eqs. (80) and (81) for redox reactions by considering different activation energies for the reduction and oxidation rates. This has been deliberately done considering the fact that the equation used in previous models to describe the diffusion velocity is similar in structure to the redox equations. As a result, the activation energies for both are combined together to consider a single activation energy which includes diffusion as well.

The current in the RRAM device is considered to be comprised of two components, the CF current (i CF ) and the oxide current (i ox ) In CF current, contribution of the CF resistance (RCF ) and the oxide surrounding the CF included in the CF radius are given as [93]:
$$ {i}_{\mathrm{CF}}=\frac{v_{\mathrm{app}}}{R_{\mathrm{CF}}\Big\Vert {R}_{\mathrm{ox}}} $$ (87) $$ {R}_{\mathrm{CF}}=\frac{L}{\pi \times {\sigma}_{\mathrm{CF}}\times {r}_{{\mathrm{CF}}_T}\times {r}_{{\mathrm{CF}}_B}} $$ (88) $$ {R}_{\mathrm{ox}}=\left\{\begin{array}{c}\frac{L\times \beta }{2\times {\sigma}_{\mathrm{ox}}\times \pi \times {r}_{\mathrm{CF}\mathrm{ext}}\times \left({r}_{{\mathrm{CF}}_B}-{r}_{{\mathrm{CF}}_T}\right)},\kern0.75em {r}_{{\mathrm{CF}}_T}\ne {r}_{{\mathrm{CF}}_B}\\ {}\frac{L}{\sigma_{\mathrm{ox}}\times \pi \times \left({r}_{\mathrm{CF}\mathrm{ext}}^2-{r}_{{\mathrm{CF}}_T}^2\right)},\kern4.25em {r}_{{\mathrm{CF}}_T}={r}_{{\mathrm{CF}}_B}\end{array}\right. $$ (89)
The other component in the device current, i.e., the oxide current represents the current throughout the oxide accounting for the whole area of the device except the one occupied the CF. The oxide current is described as [93]:
$$ {i}_{\mathrm{ox}}=\operatorname{sign}\ \left({v}_{\mathrm{app}}\right)\times {A}_{\mathrm{HRS}}\times {S}_{\mathrm{Cell}}\times {\left\{\frac{\mid {v}_{\mathrm{app}}\mid }{L}\right\}}^{\alpha_{\mathrm{HRS}}} $$ (90)
Here, A HRS e α HRS are fitting constants. The experimental and simulated results from the model are illustrated in Fig. 15 [93]. The simulated I - V curves were tuned using device parameters to fit the reported experimental data [180]. The major device parameter which was tuned the CF radius. Simulated results are almost a perfect fit to the experiment results which highlights the flexibility of the model with the proper tuning of parameters.

a RRAM current versus voltage for different cycles; the CF features are tuned to fit different curves such as r CFTM  = 45 nm, r CFTm  = 15 nm, and r CFBM  = 7 nm for cycle 3, r CFTM  = 45 nm, r CFTm  = 25 nm, and r CFBM  = 7 nm for cycle 31, and r CFTM  = 40 nm, r CFTm  = 35 nm, and r CFBM = 7.8 nm for cycle 60; b I - V curves for the devices under consideration along with different modeled data (the CFs’ dimensions considered to fit the whole set of curves have been generated randomly). The temperature at the most representative points of the I–V curve, highlighted in a cycle 3:T SET (T,B)  = [579,656] K, T LRS+ (T,B)  = [349,368] K, T LRS-(T,B)  = [466,525] K, T RESET (T,B)  = [565,621] K; cycle 31:T SET (T,B)  = [573,669] K, T LRS+ (T,B)  = [349,368] K, T LRS- (T,B)  = [460,525] K, T RESET (T,B)  = [566,622] K; cycle 60:T SET (T,B)  = [564,685] K, T LRS+ (T,B)  = [347,367] K, T LRS- (T,B)  = [460,525] K, T RESET (T,B)  = [545,641] K [93]

When considering electrochemical redox reactions, temperature is a very important factor. As can be seen from the equations describing the set/reset process, temperature parameter is a defining factor. Modeling using this concept, temperature in the CF is described through a one-dimensional approach as given as [93]:
$$ {\sigma}_{\mathrm{CF}}\times E{(x)}^2=\frac{2{c}_h}{r_{\mathrm{CF}}(x)}\left({T}_{\mathrm{CF}}(x)-{T}_{\mathrm{ox}}\right)-k\frac{\partial^2{T}_{\mathrm{CF}}(x)}{\partial {x}^2} $$ (91)
However, this approach is not the best one for truncated cones. It leads to increased complexity and improper calculation of temperature. So, Gonzelez-Cordero et al. [93] proposed a different approach where they have developed simplified analytical equations which are suitable for simulation. In the previous Bocquet model [92] they have assumed a cylindrical uniform geometry where the calculated unique temperature is uniform throughout the CF. But to consider a more detailed physical description in this model; they have considered two temperatures in the CF. One temperature is at the hottest CF point where the reset process ruptures the CF and the other temperature is considering the main CF volume. This is a more reasonable model for a truncated cone structure as the two radii grow independently of each other depending on the oxidation/reduction processes [171, 181].

To mask this concept into a simplified analytic expression, a simplifying assumption is made. Uniform cylinders are considered to represent the truncated cone CF; one corresponding to the larger radius of the cone, i.e., r CF,T and the other to the smaller radius, i.e., r CF,B . So, the maximum temperature calculated for each cylinder considering the main CF volume (T CF,T ) and at the narrowest part (TCF,B ) is given as [93]:
$$ {T}_{{\mathrm{CF}}_{\left(T,B\right)}}={T}_{\mathrm{amb}}+\frac{\sigma_{{\mathrm{CF}}_{\left(T,B\right)}}\times {r}_{{\mathrm{CF}}_{\left(T,B\right)}}\times {E}_{\left(T,B\right)}^2}{2\times {c}_h}\left(1-\cos h\left(\frac{\alpha_{\left(T,B\right)}L}{2}\right)\right)+\frac{d{T}_{0_{\left(T,B\right)}}}{\alpha_{\left(T,B\right)}}\sin h\left(\frac{\alpha_{\left(T,B\right)}L}{2}\right) $$ (92) $$ {\alpha}_{\left(T,B\right)}=\sqrt{\frac{2h}{k\times {r}_{{\mathrm{CF}}_{\left(T,B\right)}}}} $$ (93) $$ d{T}_{0_{\left(T,B\right)}}=\frac{\sigma_{{\mathrm{CF}}_{\left(T,B\right)}}\times {r}_{{\mathrm{CF}}_{\left(T,B\right)}}\times {\xi}_{\left(T,B\right)}^2\tan h\left(\frac{\alpha_{\left(T,B\right)}L}{2}\right)}{\sqrt{2\times k\times {c}_h\times {r}_{{\mathrm{CF}}_{\left(T,B\right)}}}} $$ (94) $$ {\sigma}_{{\mathrm{CF}}_{\left(T,B\right)}}=\frac{\sigma_{{\mathrm{CF}}_0}}{1+{\alpha}_T\left({T}_{{\mathrm{CF}}_{\left(T,B\right)}}-{T}_{\mathrm{amb}}\right)} $$ (95)
Here, α T is the conductivity temperature coefficient. The model is simulated and compared with the results from the previous Bocquet bipolar model [93] on which it is based. The results compare the findings from the model considering a cylindrical CF to the one considering the truncated cone CF. There is some evidence [93] presenting a better fit with the experimental data for this particular model as compared to the previous models where cylindrical CFs are considered and also results pertaining to the cases where multiple CFs are also presented; this shows the model’s flexibility in accommodating devices where multiple CFs is existent.

Detailed simulation results presenting the variation of current, voltage and CF radius are shown in Fig. 16 [93]. The variation of the applied voltage and the device current with time is shown in Fig. 16a. The complete simulated I-V curve can be seen in Fig. 16b which follows the form of a hysteresis loop, suggesting resistive switching behavior. Variations of the CF radius with time and applied voltage are presented in Fig. 16c, d, respectively.

a Applied triangular voltage signal (blue) and corresponding device current (red) versus time. b The RRAM current against the applied voltage. c Top and bottom CF radii versus time for the devices under consideration. d Top and bottom CF radii corresponding to figure 7(c) versus applied voltage [93]

Primary aim of the model was to be simplistic enough to be implemented in electric circuit simulators. Analytical equations are properly laid out to be used in SPICE simulations and it has been represented through a 1T/1R circuit. The model is also represented through a Verilog-A representation [93] which shows its applicability in digital circuits as well.

RRAM Models Based on Unipolar Devices

Random Circuit Breaker Network Model


In 2008, Noh’s et al. [182] proposed a random circuit breaker model to explain the switching in unipolar resistive switching devices. This model evolved to clear the considerable debate regarding the switching mechanism in unipolar devices, at the early stage of the study mechanism of resistive switching memory. Some reported that switching is the result of homogeneous/non-homogeneous transition of current distribution, while some other says due to the formation and rupture of conducting filaments. A new percolation model was reported by Noh’s group [182] in this regard which was based on a network of circuit breakers with two switchable metastable states. The device used for the study is a polycrystalline TiO2 RRAM device. It shows wide distributions of SET and RESET voltage with uniform resistance change at the particular transition voltage. Conductive atomic force microscopy (C-AFM) tip was used as a top electrode, and external voltage was applied through it for the resistive switching operations.

The C-AFM images along with the switching curve are shown in Fig. 17a, b. The current map shows the state of the system in the LRS and HRS. In LRS, many conductive spots can be observed in the current map which can be considered to be conducting filaments. Corresponding to the HRS, the conducting spots are vanishing, which could be translated as the rupture of the filaments. This behavior was qualitatively described by a percolation model comprising of random circuit breakers (RCB), termed as RCB network model, and is shown in Fig. 17c, d. In this model, the RRAM system is considered as a combination of a number of circuit breaker, which can have either two resistance states, high (OFF state) or low (ON state). When the ON state circuit breaker receives the RESET voltage, it switches to the OFF state. Conversely, when the OFF state circuit breaker receives the SET voltage, it switches to the ON state.

Conducting AFM (C-AFM) topographs of TiO2 filmes. a Schematic of C-AFM measurements (b ) Typical I - V curves of TiO2 films using C-AFM tip as a top electrode. c At V tip = 8 V, mapping of the current flow through the surface just after the forming operationshows locally distributed conducting regions. d TiO 2 surface in the HRS shows locally distributed conducting regions disappear after reset operation with V tip  = 1 V [182]

This model basically laid the foundation for future percolation models used to describe RRAM switching behavior, since it is dealt with the stochastic reversible dynamic processes. Most percolation models either investigate static cluster topology problems or dynamic percolation problems. A combination of reversible and dynamic processes is quite interesting. This also enabled future model developers to account for stochastic switching in the physical equations describing unipolar RRAM devices.

Filament Dissolution Model


Filament dissolution model was proposed by Russo et al. [82,83,84] exclusively for unipolar RRAM devices, although later revisions by the same group [85, 86] made this model suitable for bipolar devices as well. This model is based on the fundamental concept of Joule heating and filament temperature change. The model primarily focuses on the RESET transition of the devices, i.e., the transition from LRS to HRS. This is because of the high resistance associated with the RESET transition and this is where the major physical transformation in the device takes place. The model is based on the concept of conductive filament ruptured or dissolute under the effect of significant temperature change [84]. This temperature change in the filament is caused due to Joule heating. The proposed filament dissolution model has been deemed as self-accelerated due to the process of the rupture of filament accelerates by itself under suitable conditions.

Major advantage of this model is that it makes use of simple well-known coupled partial differential equations which describe the various effects in the device. The model is applied on a NiO based unipolar system [82, 83] where the oxide layer is sandwiched between two Pt electrodes and the filament is considered to grow from one end of the electrode to the other. Temperature profile in the oxide layer across its geometry is considered as parabolic; meaning that the temperature in the filament is minimum at the electrodes and maximum at the middle.

The mechanism behind the filament dissolution can be explained by the basic concept of Joule heating and dissolution which acts as an activator for the CF rupture. With the application of bias across the top electrode of the device, heat is produced in the filament due to the current flowing. The temperature steadily rises with an increase in the bias and when the bias reaches a significant level called reset voltage, the temperature rises above a value called the critical temperature. At this point, the dissolution of the filament is activated and the filament gets ruptured at a very fast rate leading to the device reaching a HRS.

Filament dissolution model uses coupled partial differential equations to describe the current in the device, the temperature changes due to Joule heating and the dissolution velocity. The current conduction in the device is described by the Poisson’s equation [83]:
$$ \nabla \times \left(\frac{1}{\rho}\nabla v\right)=0 $$ (96)
Here, ρ is the resistivity of the oxide and v the electric potential developed in the device due to the application of an external bias voltage v term . The voltage bias is applied at one of the electrodes, while the other electrode is connected to ground which act as boundary conditions for the device. NiO is the switching oxide, and the CF is formed as a sub-layer comprising primarily of metal ions and oxygen vacancies. The CF is considered to have a diameter of φ d .

As a result of the potential across the electrodes, heat is produced in the device due to Joule heating which leads to an increase in the temperature inside the CF. The CF geometry is divided into a number of mesh grids, to compare the ionic motions during filament formation and dissolution. The temperature is calculated at each part of the mesh grid to describe the thermal dissolution or rupture of the CF. This effect is described by the Fourier steady state heat equation given as [83]:
$$ -\nabla \times \left(k\nabla T\right)=\rho {J}^2 $$ (97)
Where k is the thermal conductivity of the oxide layer, T is the device temperature and J is current density. Thermal conductivity and electrical resistivity values are dependent on the position they are applied in. So, ρ  = ρ CF inside the CF while it is equal to ρ OX in the oxide layer. The same analogy applies for the values of thermal conductivity as well. The temperature is considered to be equal to room temperature T 0 at the electrodes, i.e., they act as heat sinks.

The device temperature T increases up to the point of reset, where it reaches critical temperature, after which the CF dissolution takes place. As a result of the rupture of CF, the current conduction is interrupted which transitions the device into HRS or reset state. This dissolution factor is modeled as [83]:
$$ {\nu}_{\mathrm{DIS}}={\nu}_{\mathrm{DIS}-\mathrm{F}}{e}^{-\frac{E_a}{k_BT}} $$ (98)
Where, E a is the activation energy, k B is the Boltzmann constant, v DIS-F is a fitting parameter and v DIS is the velocity of the CF boundary toward the symmetry axis.

The temperature dependent resistivity in the CF can be described as [83]:
$$ {\rho}_{\mathrm{CF}}(T)={\rho}_{\mathrm{CF}-\mathrm{RT}}\left[1+c\left(T-{T}_0\right)\right] $$ (99)
where c is the experimentally calculated temperature coefficient of resistivity and ρ CF-RT the standard CF resistivity at room temperature.

The coupled equations defined in the model are self-continuously solved using the numerical solver COMSOL [183, 184] Multiphysics. This software is well suited to handle these types of simulations owing to its Multiphysics capabilities. The model here uses mechanisms from electrostatics as well as heat transfer, so to simultaneously being able to handle multiple physical phenomenon’s is a big advantage on the part of the software. The obtained simulation results are shown in Fig. 18 [185]. The results show the variation in CF temperature with the change in voltage for three samples [185]. It follows the expected pattern of the model where the CF is dissolved after reaching the reset state.

Calculated temperature maps extracted from the electro thermal simulations for applied voltages a 0.4 V, b 0.45 V, c 0.575 V, and d 0.6 V for sample P (500 °C annealed); e 0.3 V, f 0.375 V, g 0.45 V, and h 0.475 V for sample Q (400 °C annealed); eu 0.3 V, j 0.4 V, k 0.55 V, and l 0.575 V for sample R (300 °C annealed); m 0.15 V, n 0.25 V, o 0.375 V, and p 0.4 V for sample S (200 °C annealed) across the CF geometry of length t ox = 20 nm. Different bias points are considered for the devices so that the reset transition for each device can be visualized easily [185]

Filament dissolution model discussed here has been modified and presented for bipolar RRAM devices by Larentis et al. [85, 86]. It is based on the same temperature and field accelerated ion migration. The set and reset processes in the device are defined by the mechanisms of drift migration induced by local electric field, ionic/electronic conduction and Joule heating. This is a point of departure from the model for the unipolar devices where the switching mechanisms for the set state were not properly defined and understood.

As a result, an equation explaining the ion migration given in Eq. (48) is introduced in the model. It includes both drift and diffusion components which rely on ion hopping. Rate of drift and diffusion is generally governed by the external applied bias and the amount of barrier lowering caused by it which is critical for the process of ion hopping. Barrier lowering generally takes place in the direction of applied electric field F , and the ion hopping depends exponentially on this energy barrier value. This causes significant migration of ions in the direction of the electric field F . This is a very simplistic and acceptable physical approach to explain the filament formation due to the charge carriers such as metal ions or oxygen vacancies. This also explains the dependence of the filament formation on the polarity of the external bias in bipolar devices. The ion-migration factor can be defined as [86]:
$$ {J}_D={J}_{\mathrm{diff}}+{J}_{\mathrm{drift}}=-{D}_s\nabla {n}_D+\mu E{n}_D $$ (100)
Where, n D is the doping density, D s is ion diffusivity, E is applied electric field and μ is the ion mobility. The temperature activated ion diffusivity, based on the Arrhenius law and ion mobility [85] is given as [86]:
$$ {D}_s={D}_0{e}^{-\frac{E_a}{kT}} $$ (101) $$ \mu =\frac{q{D}_s}{kT} $$ (102)
Using the above equations, the continuity equation for the drift-diffusion is given as [86]:
$$ \frac{\partial {n}_D}{\partial t}=\nabla \times \left({D}_s\nabla {n}_D-\mu E{n}_D\right) $$ (103)
Another difference from the model applied for unipolar devices is the use of conductivity values rather than resistivity. Also, resistivity had a linear variance with temperature in the previous model [82, 83] while in this case conductivity is modeled to be exponentially dependent on temperature given as [86]:
$$ \sigma ={\sigma}_0{e}^{-\frac{E_{ac}}{k_BT}} $$ (104)
The current conduction defined in Eq. (95) and Joule heating described in Eq. (96) are used in the model for the bipolar devices. This suggests an assumption that the temperature profile for both types of devices follows a similar pattern. Along with it, the current conduction mechanism is also assumed to be similar. This in a sense might be an over-simplified assumption because many of the models described for bipolar resistive switching devices have been unable to be used for unipolar RRAM devices [90,91,92, 186] due to a marked mismatch in the conduction and switching mechanisms.

The differential Eqs. (95)–(98) are simulated in COMSOL Multiphysics considering a 3-D cylindrical symmetrical geometry. The oxide material system considered here for reference is HfO2 and the obtained simulation results are verified with the experimental results as shown in Fig. 19 [86]. Device shows bipolar switching characteristics owing to the extra terms added in the model equations. There is somewhat of a good match between the experimental data and the calculated results. This suggests that this model could be potentially used to model bipolar devices. It is quite comprehensive in its definition of device parameters and also agrees with experimental data.

a Measured and b calculated I - V curves, for reset (V > 0) and set (V < 0) transitions, obtained by applying triangular voltage sweeps. For preparing a reset state with variable R , the reset sweep applied to an initial set state with R  = 400 Ω, is interrupted at V stop . Then, the set sweep is applied, showing that V set increases with V stop , hence with R [86]

Bocquet Unipolar Model


This model was developed by Bocquet et al. [90] for describing both the set and reset processes in unipolar RRAM devices. It is basically a modified extension of the model proposed by Russo et al. [82, 83] in the sense that it can model both the transitions of the RRAM device while the former only considers reset transition. For set process, a local electrochemical reduction of the oxide is considered to be responsible for formation of conductive filaments. However, the reset mechanism follows the tried and tested formula for unipolar devices which considers thermally assisted destruction of the formed metallic filaments by Joule heating as the primary mechanism. Also, it has to be mentioned that the model proposes equations which are analytical in nature and can be conveniently solved in an electric circuit solver.

The CF growth and rupture in the device is described by a local redox reaction and a thermally assisted diffusion respectively [82, 83, 105, 161] for a NiO based unipolar system. The reaction velocities for the reduction and oxidation processes during CF growth are based on the Butler-Volmer equation [162]. They are given by [90]:
$$ {\nu}_{\mathrm{red}}={k}_0{e}^{-\frac{\Delta r{G}_0+2\left(1-{\alpha}_s\right)F\left(E-{E}_{eq}\right)}{P\times {T}_{\mathrm{ox}}}}\left(1-{C}_{Ni}\right) $$ (105) $$ {\nu}_{\mathrm{ox}}={k}_0{e}^{-\frac{\Delta r{G}_0-2{\alpha}_sF\left(E-{E}_{\mathrm{eq}}\right)}{P\times {T}_{\mathrm{CF}}(x)}}\times {C}_{Ni} $$ (106)
where ΔrG0 is the reaction free energy at equilibrium, E eq is the equilibrium constant, α s is the asymmetry factor, P is the ideal gas constant, k 0 kinetics constant of chemical reaction, C Ni is the dimensionless concentration of metallic species, and F is Faraday constant.

The above redox reaction describing Eqs. (104) and (105) can account for the CF growth, i.e., the SET process. But to explain the CF dissolution or the RESET process, it uses the filament dissolution model by Russo et al. [82, 83] as described further in the “Filament Dissolution Model” section. As discussed, the filament dissolution model uses the Joule heating mechanism to account for the reset transition in the device. The local diffusion velocity (νdiff ) which is the governing equation for the filament dissolution model is given as [90]:
$$ {\nu}_{\mathrm{diff}}={k}_{\mathrm{diff}}\times {e}^{-\frac{E_q}{k_b\times {T}_{\mathrm{CF}}(x)}}\times {C}_{Ni} $$ (107)
As the dissolution velocity is exponentially dependent on temperature, it gets activated only when there is significant amount of temperature. Temperature value is high only when there is a comparable amount of voltage, i.e., the reset voltage has to cross a critical value during CF dissolution. This acts as means of a self-activation voltage threshold where the voltage controls the CF dissolution.

The set of reduction, oxidation and dissolution equations are coupled together in a master equation which controls the switching rate of the device given as [90]:
$$ \frac{d{C}_{Ni}}{dt}={\nu}_{\mathrm{red}}-{\nu}_{\mathrm{ox}}-{\nu}_{\mathrm{diff}} $$ (108)
The coupled equations in this model have to be solved simultaneously and continuously due to the fact that the model relies on self-consistent kinetic equations accounting for both CF growth and destruction mechanisms. This is a key feature which has to be implemented when using simulation tools to attain numerical accuracy.

Current conduction and temperature change in the device are described using simple current and heat flow equations. Current in the device on application of a voltage v cell across the electrodes is given by [90]:
$$ {i}_0=\frac{v_{\mathrm{cell}}}{\int_0^{t_{\mathrm{ox}}}R(x) dx} $$ (109) $$ R(x)=\frac{1}{r_{CF}^2(x)\pi \left({\sigma}_{CF}(x)-{\sigma}_{\mathrm{ox}}\right)+{r}_{{\mathrm{CF}}_{\mathrm{max}}}^2\times \pi {\sigma}_{\mathrm{ox}}} $$ (110)
Current flowing in the device gives rise to a temperature due to Joule heating, and this effect is modeled as [90]:
$$ {\sigma}_{\mathrm{CF}}(x)\times E{(x)}^2=-k\frac{\partial^2{T}_{\mathrm{CF}}(x)}{\partial {x}^2}+{c}_h\frac{T_{\mathrm{CF}}(x)-{T}_{\mathrm{ox}}}{t_{\mathrm{ox}}} $$ (111)
Numerical values obtained from simulation in general profoundly interlinked between the set and reset transitions. This is because from a practical stand point the CF profile obtained after the set operation is used as the initial state to simulate the subsequent reset operation. Also, the reset current and LRS resistance depends significantly on the maximum current reached during the previous set operation [187,188,189]. This is basically due to the minimization of the CF radius which subsequently increases the resistance of the device [190].

The Bocquet unipolar model is compared against a NiO system similar to the one used in the previous Filament dissolution model [82, 83]. It is to be noted that the model is applicable for a unipolar device only. But the comparison with the NiO system is limited to a single system using a numerical solver. This is a major shortcoming in this particular model regarding the non-availability of exact experimental characteristics data comparison from other sources to calibrate the model. It means that the fitting parameters have not been tested for a variety of characterization data or other models as well. So, it is difficult to judge the accuracy and viability of the model even though it uses some interesting concepts to explain the switching process in unipolar devices.

Window Function Models


Window functions are introduced in the “Linear Ion Drift Model” section [3]. These functions are generally required to limit the values that the internal state variable can reach. The dynamics of the state variable governs the switching property of the device. So, the state variable has to be set bounds within which it can grow so that the device always remains in the permissible state and does not go out of bounds. For example, if the growth/rupture of CF is being modeled, the CF physically can only grow from one electrode end to the other. If the model growth crossed that limit, it suggests a mismatch between the physical phenomenon and the model. As a result, certain window functions [94,95,96,97,98,99] which acts as limiting functions is introduced into the model to set bounds for the device.

This is further required as near to the boundaries non-linear dopant drift effects take charge and heavily suppress the speed of the ions. As a result, it leads to a non-uniform and capricious rate of change of the state variable. To properly account for all the effects, an effective window function should have the following characteristics [98]:

  • ➢Consider and take into account the boundary conditions at the top and bottom electrodes of the device.

  • ➢Be capable of imposing non-linear drift over the entire active core of the device.

  • ➢Provide linkage between the linear and non-linear dopant drift models.

  • ➢Be scalable, meaning a range of fmax ( x ) can be obtained such that 0 ≤ fmax ( x ) ≤ 1.

  • ➢Utilize an in-built control parameter for adjusting the model.

Window functions are generally multiplied with the particular I - V equation or state variable equation and are designed in a way so as to bind the complete equation within a certain boundary. There has been various window function based models proposed. In this section, we will review the most important and popular among them. We have also reported a comparative analysis among all the window functions in Table 2.

Joglekar Window


One of the very first window functions proposed by Yogesh Joglekar and Stephen Wolf [94] is based on the linear ion drift model [3]. It was developed when memristors were still in its early stages of development after the breakthrough by the HP team proposed linear model. Window functions are aim to generalize the behavior of the model around the device boundary.

Linear ion drift model [3] in general accounted for majority of the characteristics of the memristor but where it fell short was at the physical device boundary. Behavior of the devices at their physical boundaries was much more non-linear and this was left unaccounted in linear drift model [3]. Joglekar et al. [94] window function aimed to counter this limitation and generalize the non-linear behavior. The ions mobility is significantly higher in bulk of the memristor, but when it comes closed to boundary, the speed is highly suppressed. Joglekar et al. [94] proposed a modified equation for the linear model including the window function which accounted for this effect, is given as [94]:
$$ \frac{dw}{dt}=\eta \frac{\mu_D{R}_{\mathrm{ON}}}{D}i(t)F\left(\frac{w}{D}\right) $$ (112)
This modified equation reflects the speed suppression at the edges, i.e., w~ 0 or w~D . Aqui, F ( x ) is a window function which satisfies the conditions F (0) = F (1) = 0 so that there is no drift at the boundaries. The function is also symmetric about x = ½ and increases monotonically over the interval 0 ≤ x  ≤ 1/2, 0 ≤ F ( x ) ≤ 1 = F ( x  = 1/2). The window function F ( x ) parameterized by a positive integer p, is defined as [94]:
$$ F(x)=1-{\left(2x-1\right)}^{2p} $$ (113)
Simulated results of the window function are shown in Fig. 20 [94]. The results obtained from the window function implementation, it was understood that with a significantly high value of p , the function F ( x ) provides excellent generalization of the linear drift model 3 without any of its constraints. So, the window function works best at a high value of p, where it models the linear model accurately and also accounts for the non-linear effects at the boundary. This can be particularly seen in the variation of the function. With an increase in p , F ( x ) stays constant over an increasing interval around x  = 1/2 and when r  → ∞, F ( x ) = 1 for all x except 0 and 1 which are the boundary conditions. So, Eq. (112) models the general memristive systems perfectly at these particular conditions of p .

Theoretical i -v curves for a memristor with (realistic) dopant drift modeled by window functions F p ( x ) = (1 − (2x − 1)) 2p considering p  = 1 (red solid) and p  = 10 (green dashed). An external voltage v ( t ) = 2v0 sin(ω 0 t /2) is applied. The memristor parameters are w 0 /D  = 0.5 and R OFF / R ON  = 50. The memristive behavior at p  = 10 has been enhanced. The slope of the i -v curve at a small time period is the same, R 0 −1 , in both cases whereas the slope on return sweep depends on the window function. For a large p value, the return-sweep slope is R ON −1  = 1> > R 0 −1 , corresponding to a fully doped memristor [94]

This characteristic also acts as a significant limitation for the model. On the one hand, at low values of p , the window function does not perform as per expectation. But at significantly high values of p where the non-linear effects are taken into account, the difference between the linear drift and non-linear drift in the model disappears. This means that there is no proper way to account for both the linear and non-linear drift effects at the same time, while implementing this window function. Also, the mobility at the boundaries was suppressed down to zero, which made the function to be stuck at the 0 value at the terminal states. Shortcomings of this function and the improvements made over it by the Biolek window function [95] are discussed further.

Window function implemented here should be understood as extensions to the physics based models. General limitation with the physics based models is that the models do not account for the effects at the physical boundaries of the device. Memristive devices have been found to behave differently in the bulk and the boundary of the device. So, these window functions can help to overcome these limitations in the device by setting a boundary for the model and properly accounting the boundary effects.

Joglekar et al. [94] also investigated memristive device implementation in standard fundamental circuits along with the other basic elements R, L and C. Combining these four basic elements (memristor, inductor, resistor, and capacitor), the functioning of standard circuits such as MC, MLC, etc. were studied. They came to a similar conclusion as the HP team [3] that the primary property of a memristor is the memory of the charge that has passed through it. The memristor was dimensionally characterized as a magnetic flux D 2 D onde D is the memristor size and μ is the mobility.

Biolek Window Function


Biolek window function [95] developed by Zdenek and Dalibor Biolek in 2009 was modeled on the proposed memristor equations by Strukov et al. [3] Its primary aim was to provide a marked improvement over the previous Joglekar model [94] and also provided a SPICE implementation of the Linear drift model [3]. They proposed changes to the way the window functions are defined so that a closer approximation between the model and the real circuit element can be achieved. They also reported the SPICE implementation of the linear drift model [3], which opened up the model to a wide range of circuit applications at that time.

Strukov et al. [3] first pin-pointed the pertinent problems in the Joglekar function while implementing in SPICE and then proposed improvements over it. First major problem with Joglekar function is its way of setting up of the terminal state RON and ROFF . State equation and the window function defined, respectively, in Eqs. (113) and (114) bound the value of the variable to 0 at the boundary and it is forced to hold that value. This state cannot even be changed by an external stimulus. This happens due to the HP memristor remembering the x -coordinate of the boundary between two layers and not the amount of electric charge passed through it. As a result, when a new set or reset transition is to be started from a terminal value, the device has to start from 0 and not the actual value it had in its previous state.

Second problem of the window function is noticed when the model is implemented as an actual circuit element. The circuit component exactly remembers the entire charge which is passing through it. So, in case of the Joglekar window function, to transpose the memristor from a state x 0 to x 1 , a certain amount of charge q is required. Now the same amount of charge but in the opposite polarity, i.e., –q is required to bring the memristor back from x 1 to x 0 . Thus, when a memristor is being driven by a constant current with a time interval, say t , the same time t is also required for restoring the device back to its original state.

This occurs regardless of the fact that the device could be in its terminal state all the while when the current flows. This leads to significant operating delays as documented by the SPICE simulation presented by Biolek et al. [95]. Also, when the current direction is reversed, the boundaries start to move in an opposite direction regardless of the past state, thus the state is lost along another curve.

Window function being added to a particular model should be able to enhance its accountability for most of the device characteristics and get rid of the arising discrepancies. Window function proposed by Biolek et al. [95] works in this regard, defined by a positive integer p , memristor current i , and is given as [95]:
$$ f(x)=1-{\left(x- stp\left(-i\right)\ \right)}^{2p} $$ (114) $$ stp(i)=\left\{\begin{array}{c}1\kern0.5em \mathrm{for}\ i\ge 0\\ {}0\kern0.5em \mathrm{for}\ i<0\end{array}\right. $$ (115)
If the function increases the width of the doped layer, or x  → 1, the current is positive. The function value is 0 at the boundaries. When increasing the value of p , the function yields a flat window with steep troughs to zeros at x  = 0 and x  = 1.

Figure 21 [95] shows the simulated results of the window function implementation. Tuning parameter p can be used to fine tune the function accordingly to fit to different models. The obtained results satisfy all the conclusions obtained by the linear drift model [3]. The only critical limitation is the inability to account for the hard-switching effects governed by non-linear ionic drift. This means that a symmetrical hysteresis loop obtained in the models cannot be achieved. This is an inherent disadvantage the Biolek [95] window function has in its definition. So, in hindsight none of the functions were truly applicable to the proper extent in the linear drift model [3]. It was a limitation of the time which was further accentuated by the lack of detailed information about the non-linear ionic drift.

A proposed new window function demonstrated in the model (taking p  = 2) [95]

Benderli-Way Window Function


Another window function based on the basic HP model was proposed by Benderli and Wey [96] in 2009. The end result they set out to get was similar to the Biolek function [95]. The developers wanted to develop a SPICE compatible macro model based on the HP memristor which would be suitable for applications in circuit simulations. They proposed a clipping circuit which will bind it within the constraints of the length of the device (D )

The proposed clipping circuit was comprised of four comparators whose job was to ensure that the state variable function w ( t ) did not go beyond its limits. The comparators clipped w ( t ) at its top and bottom boundaries. It basically acts as a switch in which if the comparators detect a certain value in the device they activate a switch and set the device at a particular voltage. So, when w ( t ) reached the upper boundary of the device, the device is connected to a voltage source of value D , which effectively clips w ( t ) at D . This operation occurs when the voltage bias is positive.

At negative voltage bias when w ( t ) approaches the lower boundary, the circuit is connected to ground, thus clipping w ( t ) at 0. This clip is enforced until the voltage polarity is changed suggesting the correct operation is being performed. Non-linear effects at the boundaries are modeled by a proposed window function which takes into account the non-linear dopant drift. The window function is defined as [96]:
$$ f(x)=\frac{w(t)\left(D-w(t)\right)}{D^2} $$ (116)
Also by increasing the capacitance near the boundaries, the non-linear effects could be accounted for it in the circuit. The major shortcomings of this function are its simplifying approximations and the lack of a clear description of how the linear and non-linear drift can be modeled in the circuit. Although it manages to obtain a hysteresis relation for the device, it suffers from similar limitations as the Joglekar model [94]. Additionally, lack of clear information regarding the non-linear effects was an equal deterrent to the application of the function in circuits.

Shin Window Function


Shin, Kim, and Kang [97] in 2010 tried to circumvent the issue of window functions by proposing a constitutive relationship derived from the basics of the memristors developed by Chua [1]. This is different from the previously reported window functions in the sense that they tried to model the memristors perfectly by relating charge and flux together. This was the fundamental essence of the Chua model and is a stark contrast to the linear ion drift mechanics proposed by Strukov et al. [3].

Utilizing the relationships between flux (φ ) and charge (q ) in a current-controlled and voltage-controlled memristor, the window function developed keeps the model checked within the bounds. Chua model for a current controlled memristor is defined as [97]:
$$ \frac{d\varphi}{d t}=\frac{d}{d t}\left[f(q)\right]=\left\{\frac{d}{d q}f(q)\right\}\times \frac{d q}{d t} $$ (117) $$ {v}_M=\left\{\frac{d}{dq}f(q)\right\}\times {i}_M\equiv {R}_M(q)\times {i}_M $$ (118)
Here R M (q ) is the memristance defined by a derivative of the charge-flux relationship with respect to the charge. Thus, R M (q ) = df (q )/dq defines it as a current controlled memristor.

Similar to the above equations, a voltage controlled memristor is defined in terms of the charge-flux relationship q  = g ( φ ) as [97]:
$$ \frac{d q}{d t}=\frac{d}{d t}\left[g\left(\varphi \right)\right]=\left\{\frac{d}{d\varphi}g\left(\varphi \right)\right\}\times \frac{d\varphi}{d t} $$ (119) $$ {i}_M=\left\{\frac{d}{d\varphi}g\left(\varphi \right)\right\}\times {v}_M\equiv {G}_M\left(\varphi \right)\times {v}_M $$ (120)
So, on similar terms, GM ( φ ) is a voltage-controlled memductance whose values can be calculated by measuring slopes of charge-flux relationship g ( φ )

The above written Eqs. (118) and (119) can be used in compact models for circuit implementations. But they are inadequate when it is needed to define it within a bounded resistance range. This is where window functions are written to modify the circuit parameters so that the model operates within the resistance range of R MIN e R MAX . Thus, in mathematical terms it means that R M ∈ [R MIN, R MAX ] Memristor needs to be confined within the available range of resistance so as to adhere to design requirements. When the device reaches one of its boundary values, it has to stay in that state after any excess charge or flux is applied to the device. This has to be ensured so that the device does not violate its boundary conditions under hard switching conditions.

This condition is obtained in this case by ignoring the excess input current or voltage. So, in the model it is masked by a window function H so that the value of q or φ always stays within the available range. The masked input current and input voltage in the case of a current controlled and voltage controlled memristor are respectively given by [97]:
$$ H\left({i}_M\right)=\left\{\begin{array}{c}{i}_M,\kern5.5em \mathrm{if}\ {R}_M\in \left({R}_{\operatorname{MIN},\kern0.75em }{R}_{\mathrm{MAX}}\right)\\ {}0,\kern0.5em \mathrm{else}\ \mathrm{if}\ {i}_M\ \mathrm{does}\ \mathrm{not}\ \mathrm{pass}\ \mathrm{zero}\end{array}\right. $$ (121) $$ H\left({v}_M\right)=\left\{\begin{array}{c}{v}_M,\kern5.5em if\ {R}_M\in \left({R}_{\operatorname{MIN},\kern0.75em }{R}_{\mathrm{MAX}}\right)\\ {}0,\kern0.5em \mathrm{else}\ \mathrm{if}\ {v}_M\ \mathrm{does}\ \mathrm{not}\ \mathrm{pass}\ \mathrm{zero}\end{array}\right. $$ (122)
Above equations disregard any excess input current or voltage in the model space. The boundary state of the device is held until the polarity of the input source is reversed. When the polarity is reversed, it indicates the start of a new transition; the function will force the memristor to move back into memristive region. The operation is similar to the clipping circuit proposed in the Benderli-Way function [96] discussed previously, but here, there is no requirement of a complicated comparator circuit. The major purpose of using this approach to model the devices was to remove the usage of a special window function as done previously but still be able to adhere to the boundary conditions implicitly.

But the developers were also wary of the fact that some devices may not be modeled properly by their proposed constitutive relationship, which have a non-constant dynamic state variable. This was made in order to account for the non-linear drift effects discussed and experimentally shown in the HP memristor device [3].Hence, in order to take into account, the non-linear drift effects as well, they proposed a special window function. The window function aimed to overcome the backing problem faced by the Joglekar [94] function discussed previously. It has been described as [97].
$$ W(q)=\delta +\left\{1-{\left(2z-1\right)}^{2p}\right\} $$ (123)
where δ is a non-zero positive constant (δ < < 1), z is a normalized memristive charge, and z defined as z  = (q  − q MIN )/(q MAX  − q MIN ) The non-zero value of δ ensures that the state of the model returns to the normal memristive region when the input polarity is reversed. This removes the backing problem suffered by previous window functions.

Prodromakis Window Function


Prodromakis window function was proposed in 2011 by Themis Prodromakis et al. [98] of Imperial College, London, which aimed for a simple and efficient function modeled the memristor device characteristics [191] effectively. Some of the limitations and constraints of the previous models were alleviated which made the function easy and accurate to use.

The window function is considered to be parabolic in nature. It also employs a control parameter (p ) in the exponent which provides the model with the required scalability and flexibility. It also makes the window function f ( x ) scale upwards or downwards which helps create a family of distinct curves. The function f ( x ) is given in terms of pR + as [98]:
$$ f(x)=1-{\left[{\left(x-0.5\right)}^2+0.75\right]}^p $$ (124)
Control parameter is critical in this function as it helps to remove many of the constraints and limitations of the previous functions. The function can scale upwards due to the control parameter, which suggests that f máximo ( x ) can take any value between 0 and 1 inclusive. Also a very large value of p provides a linkage with the linear dopant drift effects. A serious limitation with the previous Joglekar [94] and Biolek [95] models was that the control parameter was allowed to take only integer values. But here, p could have real values as well which added more flexibility to the model.

Results obtained from the function suggest it returning a zero value at the active bi-layer edges. The drift of the dopants is also suppressed near the metal interfaces. This accounts for the non-linear drift effects at the boundaries. Other major problem of zero value stuck at the terminal state is tackled by implementing a feedback path as suggested in Eq. (123). The function can also be adjusted for any peculiar cases such as when f máximo  > 1 with the help of a second control parameter j given as [98]:
$$ f(x)=j\left(1-{\left[{\left(x-0.5\right)}^2+0.75\right]}^p\right) $$ (125)
Hysteresis loop obtained using the window function is asymmetrical which has been explained by Prodromakis et al. [98] as a result of the different switching rates of the ON and OFF rates which is quite reasonable. The hysteresis also suggests there is no terminal state problem as highlighted in the Joglekar function. Width of the doped region does not go higher than D , and the memristance is correctly limited. In the case of reverse polarity, it does not get stuck at a zero value and does not take any error states as highlighted by the results.

Results presented by Prodromakis et al., shown in Fig. 22 [98], suggest that the function can be adjusted and scaled effectively using the two control parameters in the function. The p parameter supports lateral scaling while the j parameter supports vertical scaling. Value of p has been increased from p  = 1 to p  = 80 and the results are shown in Fig. 22, which showcase the scaling features of the window function. This window function is one of the most efficient and accurate among the others which have been reported so far.

Plot of the proposed window function f ( x ) = 1 − [(x  − 0.5)2 + 0.75] p against the normalized width of the doped region x (p is considered as a variable) [98]

Boundary Condition-Based Model (BCM) Function


This window function developed in 2012 by Fernando Corinto and Alan Ascoli [99] was aimed at improving the models proposed by Joglekar [94] and Biolek [95]. They identified possible limitations with the previous functions with respect to their exhibition of single-valuedness and multi-valuedness, respectively. Also tuning the range and the boundary conditions were not possible with the Joglekar [94] and Biolek [95] functions. This was handled by the BCM window function by deriving novel methods to propose closed-form solutions for memristor devices. Along with that they also added tuning parameters to increase the flexibility of the boundary conditions used in the models.

Design of the Joglekar function limits it to a single value of memductance-flux characteristics at all input values. Similarly, the input dependent Biolek function limits it to only multi-values of the function under sign varying input. But the BCM function allows for both single-valued and multi-valued memductance-flux characteristics under a single sign varying input. The function assumes a linear dopant drift effect, which simplifies the analytic integration as well as makes it suitable for closed form solutions under any initial condition state. But this invariably neglects the non-linear boundary effects in the device. So, on the one hand, the BCM function proposes a very simplified expression for defining the boundaries of the device but it misses out on accounting for the non-linear effects due to the simplifying assumptions.

BCM model uses tuning parameters in the window function equations in an attempt to account for the non-linear effects [3, 46, 172, 173]. But they are not accurate to the same level as implicit definition of those effects. Another assumption which further simplifies the model to allow for closed form solutions is that the ionic drift rate remains constant. The BCM model is based on a window function having unitary values for all x ( t ) ∈ (0,1), and also exhibiting the vertical transitions as [99]:
$$ 1\longrightarrow 0\left\{\begin{array}{c}\mathrm{if}\exists t\mid x(t)=1\kern0.5em \mathrm{for}\ \eta v(t)\ge -{v}_{th,1}\\ {}\mathrm{if}\exists t\mid x(t)=0\kern0.5em \mathrm{for}\ \eta v(t)\le {v}_{th,0}\end{array}\right. $$ (126)
e
$$ 0\longrightarrow 1\left\{\begin{array}{c}\mathrm{if}\exists t\mid x(t)=1\kern0.5em \mathrm{for}\ \eta v(t)<-{v}_{th,1}\\ {}\mathrm{if}\exists t\mid x(t)=0\kern0.5em \mathrm{for}\ \eta v(t)>{v}_{th,0}\end{array}\right. $$ (127)
Here, η is a linear control parameter and ∃ is a quantifier denoting “there exists” which signifies that for x ( t ) there exists exactly one solution. Values of the non-negative parameters v th,1 e v th,0 determine the occurrence of such transitions. The conditions in the first Eq. (125) are established when x ( t ) obeys the boundary conditions x  = 1 and x  = 0. But the conditions in the second Eq. (126) is established when the function x ( t ) no longer obeys these boundaries.

These dynamics are encapsulated in the mathematical description defined by the tunable boundary conditions C n (n  = 1, 2, 3) as [99]:
$$ {C}_1=\left\{\begin{array}{c}x(t)\in \left(0,1\right)\ \mathrm{or}\ \left(x(t)=0\ \mathrm{and}\ \eta v(t)>{v}_{th,0}\right)\\ {}\mathrm{or}\ \left(x(t)=1\ \mathrm{and}\ \eta v(t)<-{v}_{th,1}\right)\end{array}\right\} $$ (128) $$ {C}_2=\left\{x(t)=0\ \mathrm{and}\ \eta v(t)\le {v}_{th,0}\right\} $$ (129) $$ {C}_2=\left\{x(t)=1\ \mathrm{and}\ \eta v(t)\ge {v}_{th,1}\right\} $$ (130)
So, in case of η = + 1, v th,0 is the threshold voltage. This is the minimum value of the input needs to cross, after it enters the positive region. Similarly, v th,1 represents the negative region, for η = − 1. This means that the conditions in Eqs. (127) and (128) have to be met first, before the conditions in Eq. (129) holds. The window function F is thus defined as [99]:
$$ F\left(x,\eta v\right)=\left\{\begin{array}{c}1;\mathrm{if}\ (127)\ \mathrm{holds}\\ {}0;\mathrm{if}\ (128)\ \mathrm{or}\ (129)\ \mathrm{holds}\end{array}\right. $$ (131)
The window function qualitatively works similar to other functions. At the boundary conditions, the vertical transition from 0 to 1 or vice versa occurs depending on the polarity of the input stimulus. Thus, the input used here is sign varying in nature.

Model Verification


Several ways are there to verify the working of the presented models, in this work. Some of the implementations and verification have been included with the description of the model. Models which have been described quantitatively using mathematical equations can be verified by solving the equations in a simulator. Generally, I - V characteristics of the simulated model are compared with the corresponding experimental data from the device. This gives a fair idea on the reliability of the model. Physical models described by mathematical equations can be solved by a multitude of solvers such as MATLAB, Mathematica, COMSOL, etc. Compact models which have been translated to work in the circuit level are generally simulated using the SPICE framework. There are a variety of SPICE based simulators in the market such as HSPICE, Ngspice, etc. which could be utilized. Corresponding output characteristics can be matched with the experimental results.

Physically verifying the switching mechanism in a model is trickier. It generally involves in-situ observation of the switching process [192] which requires a lot of precision and high-end equipment. However, it is very solid evidence regarding the viability of the switching mechanism presented in the model. There have been some novel methods reported to observe the growth of the conductive filaments during the switching process. Conductive atomic force microscopy (C-AFM) [182] has been used to visualize the formation and rupture of conducting filaments. I - V switching curve is shown in Fig. 17, which is clearly shows HRS/LRS states and the corresponding state of the filaments. Electrostatic force microscopy (EFM) [193] can also be used to visualize the migration and accumulation of oxygen ions by calculating the electrostatic force between the probe and the sample. It is one kind of in-situ TEM, where the focus is primarily on the charge of the carrier. Formation of conductive channels can be observed by high resolution energy-dispersive X-ray spectroscopy (EDS) which can provide accurate detailing of the composition in the filaments. These two methods have been proven to be quite effective in verifying the physical switching mechanism and the visualization of conductive filaments in RRAM devices.

Well-Posed Memristive System Definitions


An excellent work published in 2015 by Wang and Chowdhury [100] of UC Berkeley set a new paradigm for memristor and RRAM modeling. It was a push in the right direction for the whole memristor modeling community. The major features of the work were the significant improvements made on the pre-existing models. Some tweaks were proposed in the fundamental understanding of memristor models which contribute to eradicating some of the long-standing issues which has plagued the memristor models. Also, they demonstrated implementation of the models in SPICE, Verilog-A as well as their own prototyping platform based on MATLAB called as MAPP [191].

The root cause of many of the issues affecting memristor models were improper mathematical implementation. As a result, it limited their application in a variety of simulation and design scenarios. Simulation of the models is a very critical factor in determining the applicability of the model; however, the existing models were unable to be applied in a variety of simulation studies such as DC, transient etc. This work aimed at modifying the models into a form where simulating them would be a simple task.

The common ill-posed or erroneous definitions that many of the previous models suffer are not being properly defined at all biases, outputs not being unique or continuity problems. These basic problems are to be avoided in the models for wide application. All the various problems that the authors have encountered and the improvements they have presented are discussed.

A very valid point highlighted by them is that a well-written model for a particular circuit should be able to replicate its characteristics or be valid beyond its actual boundaries as well. Even if getting outputs beyond the applicable range might not be physically possible, but in simulation environments like SPICE it is imperative that the circuits work at all level of biases. This will lead to efficient simulation and produce smooth varying outputs at all biases. Many of the models we have discussed previously have sought to ignore the operation of the model beyond the available range, leading to their incompatibility for use in circuit simulators. This requirement can also be understood by the underlying algorithms of circuit simulators such as the Newton-Ralphson (NR) algorithm [194] which is commonly used for solving non-linear equations.

The NR algorithm [194] works on the principle of applying a sequence of biases to devices so that convergence can be reached for a valid solution of the circuit. So even if the bias is physically possible or not, for the NR algorithm to find a solution the model must evaluate at all bias. This invariably means that if the NR finds a solution to a well-designed model, the input bias will be physically reasonable. In the cases where the converged solution is physically not possible, it provides insights into the problematic areas of the models and is critical to troubleshooting. Therefore, in order for the NR algorithm to work correctly and find a proper solution, models should be designed to be evaluated at all biases.

Another fundamental problem is the divide-by-zero error. Many of the models have terms such as 1/(x − a ) which cause these errors. It leads to the solutions getting unbounded and causing discrepancies. Along with that, some expressions use square root (√) with negative arguments which give rise to complex arguments. Non-real numbers are not valid arguments for models and can cause non-convergence to a valid solution in simulators.

Almost all of the models we have discussed earlier do not account for a very important aspect of device model simulations. Any mathematically viable input must produce a mathematically viable output, and the most basic among this is the DC analysis. It is commonly the crux and starting point of any analysis and a proper model should produce an accurate DC solution. A proper well-designed model should work consistently with all kind of analysis. But almost all of the models suffer from significant DC response problems. So, this is another area that needs to be improved. Wang et al. [100] also addresses the problems of the models generally faced when being defined in circuit level languages such as Verilog-A.

The crux of the improvement to the models Wang et al. [100] have proposed revolve around the correct way of modeling hysteresis itself, i.e., using internal unknowns and implicit equations. This is because the dynamics of the filament in a RRAM closely follows a hysteresis characteristic. Also, the improvements make the models simulation ready with all the major analyses like DC, AC, and transient providing acceptable results. Various techniques are also proposed to aid convergence in electric simulators including a proposed new limiting function which replaces the functionality of window functions and overcomes all their limitations.

Accurate Description of Hysteresis


Memristive systems like RRAM devices have been proven to follow the features of hysteresis [1, 2] very closely. So, tweaking the models should start from the very basic. Modeling and defining the hysteresis characteristics accurately are critical to the proper functioning of RRAM models. Devices with hysteresis do not have a simple algebraic mapping between the voltage v ( t ) and current i ( t ) A state variable s ( t ) which defines the state of the device is required for the mapping given as [100]:
$$ i(t)={f}_1\left(v(t),s(t)\right) $$ (132)
A differential equation is used to govern the dynamics of the state variable, i.e., the rate of change of the state described as [100]:
$$ \frac{d}{dt}s(t)={f}_2\left(v(t),s(t)\right) $$ (133)
The above equations serve as the model template for modeling hysteresis. Value of the s ( t ) at a particular time instant t is governed by the history or the state of v ( t ) Thus, the device is considered to be having an internal memory of the input voltage. Functions f1 and f2 are chosen accordingly to define the characteristics of the device. Choice of these functions could be termed as critical in defining the dynamic of the device.

There has been a very common thought that hysteresis shows up in transient analysis only [1,2,3]. But Wang et al. [100] have demonstrated hysteresis based on the DC solutions of the models. This has been achieved by the proper selection of the governing functions f 1 and f 2 . As discussed earlier, obtaining a DC response is imperative for the proper analysis and simulation of the devices. So, obtaining a hysteresis in a DC solution makes this concept highly noteworthy. The governing functions are defined as [100]:
$$ {f}_1\left(v(t),s(t)\right)=\frac{v(t)}{R}.\left(\tanh \left(s(t)\right)+1\right) $$ (134) $$ {f}_2\left(v(t),s(t)\right)=\frac{1}{\tau }.\left(v(t)-{s}^3(t)+s(t)\right) $$ (135)
In the function f1 , hyperbolic tangent function tanh is chosen because of its monotonically increasing properties with a range of (− 1,1). The dynamics of s ( t ) is governed by the choice of f 2 . When the value of f 2 is 0, the corresponding (v , s ) values are the DC solution of the circuit. As shown in the contour plot in the Fig. 23a, [100] the curve f 2  = 0 folds back in the middle and crosses the v  = 0 axis three times. Thus, it has three stable states or three possible values of s in the DC solution which forms the foundation of the DC analysis of hysteresis. The operation of the variable s as a hysteresis curve is shown in Fig. 23b [100]. As s modulates the current, the I - V relationship will result in a hysteresis as well. Multiple stability of the state variable and abrupt changes in the DC solutions leads to the formation of a hysteresis in DC analysis. This sets a very strong foundation for accurate and efficient modeling of RRAM devices.

f 2 function in (6) plotted in contour form and predicted s -v hysteresis curve depending on the sign of f 2 [100]

Proper Definition of Internal Unknown Variables in Verilog-A


It has been discussed several times during review of the various models, implementation of the models accurately in SPICE [116, 117] and Verilog-A [142] is critical for their acceptability. This is because SPICE is the most commonly used circuit simulation platform, and Verilog-A is the most widely used hardware description language. So, simulating in these platforms is as close as it gets to the real physical devices. A major shortcoming of the previous models is the way that internal unknown variables were handled in Verilog-A.

In a memristor model, the state variable is an internal unknown because its value gets changed with different states. Verilog-A does not have a straightforward way of handling and defining these unknowns. As a result, it can get very confusing while dealing with a constantly changing value. Wang et al. [100] proposed to declare the state variable s as a voltage or potential rather than any “real value.” Some very critical points are mentioned below while handling internal unknowns in Verilog-A.

  • Different Verilog compilers handle variables declared using “real” differently. Then, this will lead to very inconsistent results.

  • Differential equations should not be defined by using the in-built idt() function. Because this function has very inconsistent support in the compilers and causes many limitations [140, 142].

  • Time integration to obtain analytical solutions should not be coded inside the model. The process is pretty simple but has serious pitfalls as given below.

  • ➢This method makes use of “abstime” function. To define the starting point of the integration it also uses “initial_step.” These have been termed as bad modeling practices in analog simulation [140, 143].

  • ➢The internal unknowns are defined as a memory state in this method, which can create problems for periodic steady state (PSS) analysis.

  • ➢This method bypasses many of the simulators built in facilities such as the convergence aids, time step control etc.

  • ➢It can cause serious convergence issues for stiff systems due to its dependence on the Forward Euler (FE) algorithm [195].

These problems are generally a combination of bad modeling practices and the incapability of Verilog-A to handle internal unknowns efficiently. As a result, declaring s as a voltage has been demonstrated as an effective way of getting around the problem.

Developing Generic RRAM Models


Taking the previously discussed hysteresis equations as a template, Wang et al. [100] presented a generic way of developing compact models for RRAM devices. To demonstrate the RRAM model, ASU/Stanford model [78, 80] is considered.

The filament gap is used as state variable. Current across the device is considered as itb , the voltage as vtb and the unknown state variable as gap . The equations are defined based on the previously stated hysteresis template given as [100]:
$$ itb(t)={f}_1\left( vtb(t), gap(t)\right) $$ (136) $$ \frac{d}{dt} gap(t)={f}_2\left( vtb(t), gap(t)\right) $$ (137)
The above equations define the physical contexts of the RRAM device. Now, choosing the proper functions for f 1 and f 2 is critical to capture the physical properties. In most of the models, we have encountered till now the application of consistent equations for f 1 with some changes in the internal unknowns. This fact is corroborated by Wang and Roychowdhury’s work [100]. They have considered the f 1 function as [100]:
$$ {f}_1\left( vtb, gap\right)={i}_0\times \exp \left(-\frac{gap}{g_0}\right)\times \sin h\left(\frac{vtb}{V_0}\right) $$ (138)
The f 2 function is considered according to the ASU/Stanford model [78]:
$$ {f}_2\left( vtb, gap\right)={v}_0\times \exp \left(-\frac{E_a}{V_T}\right)\times \sinh \left(\frac{vtb\times \gamma \times {a}_0}{t_{ox}\times {V}_T}\right) $$ (139) $$ \gamma ={\gamma}_0-\beta \times {gap}^3 $$ (140)
The γ here is the local field enhancement factor [196] which contributes to the abrupt SET (filament growth) and gradual RESET (filament dissolution). A common property among most of the RRAM models is the fact that the sign of f 2 is same as that of –sinh (vtb ) This means in terms of the gap that it starts to decrease whenever vtb is positive and vice versa. But this growth or dissolution cannot be indefinite for numerical simulation to work. For the simulations to work in reality, they have to be bounded which has been discussed in depth in the next section.

Various methods have been proposed to account the boundary effects of the devices. It will come up short with the methods having some serious limitations. Some of the models have implemented direct “if-then-else” statements in the Verilog code [80, 81]. But the problem is that the use of “if-else” statements removes the model from the differential equations framework which is not acceptable. It also introduces hard discontinuities in the model whereas we need smooth continuous curves.

A very popular way of modeling boundary conditions is the use of window functions which we have discussed in the “Window Function Models” section. The window functions as discussed work on the principle of setting f 2  = 0 when gap = maxGap and minGap . We have seen that improvements made with the window functions make them suitable for transient simulations and they produce smooth and continuous results. But the real problem with these functions is actually deep-rooted. The problem can be understood by analyzing the sign and zero-crossings of the f 2 curve as shown in Fig. 24 reported by Wang and Chowdhury [100].

Different choices of f 2 considered in the RRAM model [100]

As is seen in Fig. 24b, the f 2  = 0 curve has three lines, V  = 0, maxGap , and minGap . Beyond the values of maxGap and minGap , there will not be any stable DC solutions so they will not show up in transient analysis. Therefore, when sweeping between those values, the transient solution will work fine with the window function multiplied to f 2 . But with the other analyses there are several problems. In DC operating point analysis, unphysical solutions can show up owing to the fact that all the lines consisting of the f 2  = 0 curves are valid. So, when the voltage is zero, any arbitrary value can be the value of gap and it no longer follows the boundaries. Hence, DC analysis is a major limitation for the window functions.

This has been tackled very efficiently by Wang et al. [100] by keeping the DC solutions in a single continuous curve and then trying to bind the value of gap . As illustrated in Fig. 24c, curves A and C contain the stable solutions and B has the unstable solutions. So, when the sweeping starts after zero, the value of gap will switch between maxGap and minGap . It is mathematically represented by introducing a couple of clipping functions F clipmin and F clipmax represented as [100]:
$$ {f}_2^{\ast}\left(\mathrm{vtb},\mathrm{gap}\right)={f}_2\left(\mathrm{vtb},\mathrm{gap}\right)+{F}_{\mathrm{clipmin}}\left(\mathrm{vtb},\mathrm{gap}\right)+{F}_{\mathrm{clipmax}}\left(\mathrm{vtb},\mathrm{gap}\right) $$ (141) $$ {F}_{\mathrm{clip}\mathrm{min}}\left(\mathrm{vtb},\mathrm{gap}\right)=\left(\mathrm{safeexp}\left({K}_{\mathrm{clip}}\times \left(\mathrm{minGap}-\mathrm{gap}\right),\mathrm{maxslope}\right)-{f}_2\left(\mathrm{vtb},\mathrm{gap}\right)\right)\times {F}_{w1}\left(\mathrm{gap}\right) $$ (142) $$ {F}_{\mathrm{clip}\mathrm{max}}\left(\mathrm{vtb},\mathrm{gap}\right)=\left(-\mathrm{safeexp}\left({K}_{\mathrm{clip}}\times \left(\mathrm{gap}-\mathrm{maxGap}\right),\mathrm{maxslope}\right)-{f}_2\left(\mathrm{vtb},\mathrm{gap}\right)\right)\times {F}_{w2}\left(\mathrm{gap}\right) $$ (143)
The functions F w1 and F w2 are smoother versions of the step functions. This serves the purpose of adhering to the property required for a clipping function to work properly while maintaining a smooth continuous curve. These are described as:
$$ {F}_{w1}\left(\mathrm{gap}\right)=\mathrm{smoothstep}\left(\mathrm{minGap}-\mathrm{gap},\mathrm{smoothing}\right) $$ (144) $$ {F}_{w2}\left(\mathrm{gap}\right)=\mathrm{smoothstep}\left(\mathrm{gap}-\mathrm{maxGap},\mathrm{smoothing}\right) $$ (145)
The functions safeexp () and smoothstep () are smoother versions of the normal variants of the exp () and step () function. They have been developed by Wang et al. [100] in their MAPP [191] platform and is available to use within the platform. The clipping functions here closely mimic the actual physical effects occurring in the device. It can be termed as a huge force which keeps the state variable gap within its acceptable physical limits.

The templates provided by Wang and Chowdhury [100] for RRAM modeling is capable of widespread applicability and can be used as an ideal platform to develop other models. It consists of accurate modeling of hysteresis, includes proper handling of internal unknowns in Verilog-A and does not need to use incompatible functions like “idt()” and “initial_step” in the differential equation framework. They also circumvent the various limitations of the window functions by the use of mathematically accurate clipping functions. The model templates support a variety of analyses such as DC, AC, transient, PSS, and homotopy [197] in Verilog-A, SPICE, and MAPP. This is a very exhaustive list of advantages which should be used for development of future models by the RRAM community.

Improving Solution Convergence


Obtaining solution convergence is one of the most important features, every RRAM model ought to possess. The convergence of the solution points to the fact that it is valid and acceptable. This has been a problem for many of the compact models describing RRAM devices. Several techniques have been proposed [100] to aid convergence of solutions in these models. The use of limiting functions compatible with SPICE is very important so that it limits the solutions whenever they cross the acceptable range. A very simple way to make sure the solutions converge is to properly scale the unknowns and variables. Proper scaling makes sure that any results obtained are defined relatively accurate to the inputs.

But the major feature that aids convergence is the use of proper limiting functions in SPICE. Generally during DC operating point analysis, the Newton Raphson (NR) [194] algorithm iterations take very large values while guessing the DC solutions. This is because of the fast-growing sinh functions used in the models. This in turn leads to large values of f 1 and f 2 . Limiting functions are the best technique to circumvent this and prevent the NR iterations from taking large input values but keeping the sinh functions intact. Presently, SPICE includes pnjlim , fetlim , and limvds as limiting functions. But these are not enough to obtain numerical accuracy for RRAM models. A new limiting function dubbed as sinhlim has been proposed for this purpose. The function is based on the original pnjlim and is given as [100]:
$$ {x}_{\mathrm{lim}}=\mathrm{sinhlim}\left({x}_{\mathrm{new}},{x}_{\mathrm{old}},k\right)=\frac{1}{k}\times \ln \left({y}_{\mathrm{lim}}+\sqrt{1+{y}_{\mathrm{lim}}^2}\right) $$ (146)
The major feature here is that the limiting function does not use very large values of x new , instead it increases the value in iterations. This function is fully compatible with SPICE and can be implemented in any SPICE compatible circuit simulator. This marks a new addition to the number of limiting functions available for use in circuit simulators apart from the ones developed decades ago.

Improving Existing Models


On the basis of the accurate generic model templates discussed in the previous section, Wang and Chowdhury [100] have proposed some improvements in some of the existing popular models. We discuss them here and present it in a concise form so that it becomes easy to understand the changes and implement it forward in future models.

The improvements have been proposed on the linear ion drift [3], non-linear ion drift [46], Yakopcic [73, 74], TEAM/VTEAM [75,76,77], and ASU [78] models. Many of the models have I-V relationships in common with other models. Therefore, the most common and important among those are considered and termed as f 1 functions as discussed earlier. The state variable equations of those models are termed as f 2 functions are corresponding improvements are proposed.

Both the f 1 and f 2 functions are generally non-linear and asymmetric. And the reported models use very discontinuous and fast-growing terms like exponential, sinh , and power (pow ) functions which results in difficulties during the convergence of the solutions. So, these can be overcome by using “smooth” and “safe” functions as proposed by Wang et al. [100]. The smooth functions are used in place of discontinuous functions. Major design criteria in the smooth functions used is a common smoothing factor which combines common elementary functions to approximate the original non-smooth functions. Smoothing factor controls the trade-off between better approximation and more smoothness. The safe functions are versions of the fast-growing functions which limit the maximum slope the functions can attain, and then linearize it to keep the slopes constant beyond it. For some functions like sqrt , log , etc., the “safe” versions clip the inputs using smoothclip so that non-valid outputs can be avoided.

In the particular f1 and f2 functions used in the models, the if-else statements are replaced by smoothswitch which removes the discontinuity of the former. The exp and sinh functions are correspondingly replaced by safeexp and safesinh . The authors have demonstrated the definition of the functions in MAPP and Verilog-A which makes it easy for future model developers to integrate it into their system. A very common problem with the f 2 functions, i.e., the state variable dynamics is the uncertainty over the range of the internal unknown. The previous models counter this by either introducing window functions that bound it within a range or do not account for the effects at all. This has been very efficiently handled by introducing self-modeled clipping functions which define the acceptable range of the internal unknown.

Another major problem which is countered is the poor DC hysteresis response of the models. There have been some attempts in the Yakopcic and TEAM/VTEAM models to model this effect by setting the value of f 2  = 0 within a certain voltage range. It has been discussed that DC hysteresis occurs due to the DC solution curve folding backwards. With the approach used in the other two models, when f 2  = 0 with the voltage close to 0, there are infinite number of solutions for the state variables within that voltage range. This is a problematic approach in those unstable DC solutions are also included here which makes the equation system ill-conditioned. This will also cause problems during DC operating point analysis and homotopy analysis. The DC solutions of the models will vary from simulator to simulator because of the manner the equations are designed. This has been very efficiently handled and circumvented by Wang and Chowdhury [100] as summarized in Tables 3 and 4.

Novel RRAM Applications


There have been several new breakthroughs with RRAM architectures and applications. Among them noteworthy in case of architectures is the use of materials such as graphene, amorphous carbon films, transition metal dichalcogenides (TMDs), black phosphorous in a RRAM device. Neuromorphic computing is a novel application scheme for RRAM devices which utilize the memory retention property to use them as synapse devices.

RRAM devices based on graphene and related materials [198] have showcased performance similar to conventional metal oxide devices. These devices are different due to their unique lattice structure and belong to a completely different family of materials. So, investigation on modeling of such devices is highly necessary. Whether conventional modeling techniques such as the ones presented in this work can be used for these devices depends on the physical phenomena governing them and the corresponding I - V charcateristics. The hypotheses presented to explain the switching in graphene oxide (GO) based devices [198] are consistent with standard bipolar RRAM switching mechanism. The absorption and creation of the conductive filaments are thought to be a result of diffusion of metallic ions from the electrode to the switching layer or transport of oxygen related carriers in the switching media. RRAM devices based on amorphous carbon [198] as the switching media are thought to operate under a similar mechanism. Owing to the similarity in the nature of the physical transport mechanisms, existing physical models can be used to explain the novel GO and amorphous carbon based devices.

Neuromorphic computing [199] is a novel architecture scheme which employs RRAM as synapse devices as its fundamental building block. It is believed that these RRAM based neuromorphic systems can replicate how our brain functions, harnessing the ability of memristors to remember their state. This enables the system to be trained for specific applications just like the human brain. With RRAM forming the crux of these systems, it is critical that device characteristics for the RRAM devices are well modeled. But modeling of RRAM based synapse devices is challenging owing to the fact that the RRAM devices used need well-defined analog behavior, which is the precondition for brain like functions. Device performance under AC stress and cycle to cycle variability are factors which affect the potentiation and depression of the synapse device. Standard models reported in this work can predict the digital behavior of the RRAM, but one may implement them for the analog behavior. Though a few models are reported [200, 201] to describe the switching mechanism in analog RRAM, but it has been difficult to mathematically define them and translate it into a compact form. However, significant research is ongoing to quantitatively describe these characteristics and translate it into a compact form to be used on the circuit level.

Conclusions


In summary, the important features of all widely accepted RRAM models have been discussed. This work fulfills the requirement of the modeling community for a unified discussion on the various RRAM models. Many of the recent models, such as Stanford/ ASU model, Gonzelez-Cordero et al. model, Prodromakis model, have provided apt explanations for RRAM processes based on the early models. Implementations of different window functions like Joglekar, Biolek, and Prodromakis have been presented and compared. Various unexplained phenomena occurring in the devices are numerically validated in the models. No one model can be deemed as the perfect one, owing to the variety of materials, fabrication processes and device operations exist in the RRAM devices. Each model has been tuned accordingly to fit the device used. Researchers are still some time away from developing a generic RRAM model owing to these factors and also due to the deficiencies in the modeling techniques. Accurate and well-defined modeling techniques have been discussed in the “Well-Posed Memristive System Definitions” section, which should act as a competent template for future model development. Combined with the detailed analysis provided for past RRAM models, this review work can potentially act as a focal point for RRAM model developers.

Nanomateriais

  1. Elatec:habilita autenticação de usuário e controle de acesso para dispositivos sem porta USB
  2. Modelagem e simulação propõe novos insights para SARS-CoV-2
  3. Características de chaveamento resistivo bipolar de dispositivos RRAM de estrutura tricamada de HfO2 / TiO2 / HfO2 em substratos revestidos com Pt e TiN fabricados por deposição de camada atômica
  4. Engenharia de processo de revestimento por imersão e otimização de desempenho para dispositivos eletrocrômicos de três estados
  5. Modelagem de Dinâmica Molecular e Simulação de Corte de Cério com Diamante
  6. Estudo de primeiros princípios sobre a estabilidade e imagem STM de Borophene
  7. Estudo in vitro da influência das nanopartículas de Au nas linhas celulares HT29 e SPEV
  8. Modelagem de Materiais de Metais em Abaqus
  9. Dispositivos de segurança e considerações
  10. O poder da modelagem e simulação 3D estimula inovações nos processos de fabricação