산업 제조
산업용 사물 인터넷 | 산업자재 | 장비 유지 보수 및 수리 | 산업 프로그래밍 |
home  MfgRobots >> 산업 제조 >  >> Industrial materials >> 나노물질

저항성 랜덤 액세스 메모리의 모델링 및 시뮬레이션에 대한 종합 연구

초록

이 작업에서 우리는 RRAM(Resistive Random Access Memory)의 설계 및 설명을 위해 제안된 다양한 모델에 대한 포괄적인 논의를 제공합니다. 초기 기술인 정확한 모델에 크게 의존하여 효율적인 작업 설계를 개발하고 장치 전반에 걸쳐 구현을 표준화합니다. 이 리뷰는 RRAM 장치용 모델 개발에 고려되는 다양한 물리적 방법론에 대한 자세한 정보를 제공합니다. 지금까지 보고된 모든 중요한 모델을 다루고 기능과 한계를 설명합니다. 멤리스틱 시스템에서 발생하는 다양한 추가 효과 및 이상 현상이 해결되었으며 이러한 문제에 대한 모델이 제공하는 솔루션도 표시되었습니다. 이 작업에서는 장치 작동, 스위칭 역학, 전류-전압 관계와 같은 RRAM 모델 개발의 모든 기본 개념을 자세히 다룹니다. Chua, HP Labs, Yakopcic, TEAM, Stanford/ASU, Ielmini, Berco-Tseng 등이 제안한 인기 모델을 다양한 매개변수에 대해 광범위하게 비교 분석했습니다. Joglekar, Biolek, Prodromakis 등과 같은 창 기능의 작동 및 구현도 제시 및 비교되었습니다. 모델의 적용 가능성과 정확성을 높이는 잘 정의된 새로운 모델링 개념이 논의되었습니다. 이러한 개념의 사용은 이 작업에서 열거한 기존 모델에서 몇 가지 개선을 가져옵니다. 제시된 템플릿에 따라 매우 정확한 모델이 개발되어 미래의 모델 개발자와 모델링 커뮤니티에 큰 도움이 될 것입니다.

<섹션 데이터-제목="배경">

배경

이 새로운 컴퓨팅 시대에는 성장에 상응할 수 있는 기술이 필요합니다. 새로운 기술은 향상된 성능의 요구를 충족하고 미래의 장치에 맞게 확장할 수 있어야 합니다. Leon O. Chua가 1971년 [1]에 가정한 Memristors는 이러한 요구 사항을 충족하고 새로운 종류의 장치를 위한 토대를 마련한 것으로 보입니다. "메모리 저항기"의 약자인 멤리스터는 제공된 입력 자극의 이력에 따라 내부 저항 상태를 기억하는 기본 2단자 장치입니다. Chua는 멤리스터가 각각 전류와 전압의 시간 적분인 플럭스와 전하 사이의 관계를 특징으로 한다는 것을 고안했습니다.

1976년 후반에 Chua와 Kang[2]은 멤리스터를 일반화하여 멤리스터 시스템이라고 하는 새로운 클래스의 동적 시스템에 포함시켰습니다. 20세기 말에 이러한 장치에 대한 관심은 많은 이점에도 불구하고 시들었습니다. 이것은 부분적으로 실리콘 집적 회로 기술의 발전 때문이었습니다. 그러나 실리콘 기술의 노화와 축소를 지원할 수 없는 능력으로 인해 21세기 초에 대체 스위칭 장치에 대한 검색이 주목을 받았습니다. 이는 나노크기 물질의 성장 및 특성화의 발전에 의해 동등하게 도움이 되었습니다. 이것은 항상 미시적 멤리스티브 스위칭을 이해하는 데 상당한 진전을 이룹니다.

Memristor 기술은 Strukov et al. [3] TiO x에 대한 이론과 실험 사이의 연결고리를 확립했습니다. 기반 장치. 또한 그들은 멤리스티브 시스템의 식별 가능한 특징 중 하나인 전류-전압 관계에서 핀치 히스테리시스를 얻었다[4, 5]. 이것은 금속/산화막/금속 구조의 발자국을 따라 다양한 장치에 멤리스터 기술을 개방했습니다. 비슷한 유형의 인기 있는 장치 중 일부는 산소 RRAM(OxRRAM)[6,7,8,9,10]과 전도성 브리지 RAM(CBRAM)[11,12,13]이었습니다. 이러한 장치는 일반적으로 스위칭 메커니즘에 따라 분류됩니다.

저항성 랜덤 액세스 메모리(RRAM)

입증된 비휘발성 멤리스티브 동작이 비휘발성 메모리에 활용될 수 있기 때문에 이러한 새로운 장치에 대한 연구 관심이 높아졌습니다. 플래시 메모리 기술의 잠재적인 대안으로 간주되고 있습니다. 현재 시대의 컴퓨팅이 점점 더 많은 데이터 중심으로 진행됨에 따라 현재와 미래의 요구 사항에 보다 잘 맞는 메모리 기술에 대한 요구가 있었습니다. 여러 새로운 장치와 비교하여 RRAM 장치는 확장성이 더 뛰어나고[14,15,16,17,18], 고밀도[19,20,21,22,23,24], 저전력 소비[25,26,27] ,28,29], 더 빠르고[30,31,32,33], 더 높은 내구성과 유지력을 갖고[34,35,36,37] CMOS 호환성이 높습니다[38,39,40,41,42]. RRAM 장치는 가장 널리 사용되는 비휘발성 메모리 기술 중 하나이며 장치 작동을 실현하고 정확하고 간단한 장치 구조를 설계하기 위한 메커니즘을 이해하고 모델을 개발하기 위해 광범위한 연구가 수행되고 있습니다. 이 장치는 단순한 2단자 금속 절연체 금속(MIM) 구조이며 두 저항 상태 LRS(저저항 상태)와 HRS(고저항 상태) 사이를 전환합니다. LRS는 장치가 SET 또는 ON 상태에 있음을 나타냅니다. 대조되는 HRS는 장치가 RESET 또는 OFF 상태에 있음을 의미합니다. 이 장치의 저항 상태 전환을 통해 데이터 비트가 저장됩니다[43,44,45]. RRAM 소자는 스위칭 극성에 따라 바이폴라 소자와 유니폴라 소자로 분류할 수 있다. 유니폴라 스위칭에서 장치는 동일한 극성 바이어스로 전환되는 반면 바이폴라 스위칭에서는 두 극성의 바이어스가 모두 필요합니다.

RRAM 장치의 스위칭 메커니즘을 설명하기 위해 여러 접근 방식이 제안되었지만 이진 산화물 기반 RRAM 장치에 대해 가장 인기 있고 널리 받아 들여지는 것은 산소 이온/공석의 드리프트에 의한 국부 전도성 필라멘트(CF)의 형성 및 파열입니다. [9, 16, 46,47,48,49]. SET/RESET은 산소 이온/공석의 조합/재생의 결과로 발생합니다[50,51,52]. RRAM 장치의 성능은 활성 산화물 층의 선택에 의해 크게 영향을 받는다는 것이 입증되었습니다[53,54,55]. HfO x와 같은 다양한 산화물 시스템 , TiO x , NiO x , TaO x , ZnO x 등 [56,57,58,59,60,61,62,63,64,65,66]이 저항성 스위칭 동작을 시연하는 데 사용되었습니다. RRAM 장치가 실제로 멤리스틱 장치인지 여부에 대해 몇 가지 논란이 있었습니다. RRAM 장치의 위치를 ​​명확하게 하기 위해 Chua는 RRAM 장치가 실제로 멤리스티브 장치라는 설명을 제공했습니다[67].

RRAM 모델링의 중요성

새로운 반도체 기술을 기반으로 한 전자소자 개발에 있어 매우 중요한 측면은 모델링의 역할이다. 정확하고 포괄적인 모델은 장치 작동을 이해하고 최적의 성능을 위해 설계하고 필요한 사양과 일치하는지 확인하는 데 가장 중요합니다. 다양한 정확도, 다양한 기능 및 혼합된 결과를 가진 여러 모델이 제안되었습니다. 따라서 RRAM 장치를 위한 강력하고 유연한 모델을 설계하려는 개발자는 이전에 시도한 방법과 직면한 제약 조건에 대한 정보를 알고 있어야 합니다.

이 작업에서 우리는 다양한 RRAM 모델의 모든 기능과 특성에 대해 자세히 논의했습니다. 일반 멤리스터 모델도 RRAM 장치를 설명하는 것으로 간주됩니다[67]. 멤리스터의 기본을 제공하는 Chua 모델[1]을 시작으로 멤리스터의 기본 정의를 논의한다. HP 모델[3]이 제공하는 멤리스터 및 RRAM 장치의 돌파구에 대해 자세히 설명합니다. 비선형 효과[46, 68, 69]와 함께 이러한 장치 메커니즘의 기본을 형성하는 선형 이온 드리프트 효과가 고려됩니다. SPICE 호환 물리학 기반 모델의 토대를 마련한 Pickett-Abdalla 모델[70,71,72]에 대해 자세히 설명합니다. Yakopcic 모델[73, 74]이 채택하고 개선한 다양한 기능도 다룹니다.

필라멘트 간격을 상태 변수로 간주하는 임계값 효과[75,76,77]와 같은 새로운 기능을 도입한 모델[78,79,80,81]이 검토되었습니다. 단극 장치와 온도 효과[82,83,84]를 설명하는 일부 모델이 자세히 검토됩니다. 장치 성장 역학에 기반한 물리적 모델[85, 86]도 고려됩니다. 이와 함께 바이폴라 소자만을 고려한 모델[87,88,89], CF 크기의 변화[90, 91], 기타 많은 요인[92, 93]이 고려된다. 논의된 모든 모델에 대한 간결한 분석이 표 1에 나와 있습니다.

Joglekar [94], Biolek [95], Benderli-Wey [96], Shin [97], Prodromakis [98, 99] 등과 같은 창 기능 구현을 기반으로 하는 다양한 모델 또한 다양한 모델과 이를 극복하기 위해 후속 모델에서 사용하는 방법이 포괄적으로 제시되었습니다. RRAM 모델링을 개선하기 위해 Wang과 Roychowdhury[100]가 수행한 중요한 작업도 전체 RRAM 모델링 커뮤니티를 위한 올바른 방향으로 상당한 추진력이 있기 때문에 심도 있게 검토되었습니다. 이러한 예와 함께 다양한 플랫폼의 장치에 대한 시뮬레이션 및 검증 연구를 다룹니다. 현 단계에서 RRAM과 멤리스터 모델에 대한 가장 종합적인 리뷰입니다. 모델에 대한 설명은 바이폴라 소자와 유니폴라 소자를 설명하는 것으로 나누어져 있다. 창 함수 구현 모델은 별도의 섹션에서 설명합니다.

이전에 RRAM 장치 메커니즘[46, 101,102,103,104,105], 제조 기술[106,107,108,109], 재료 스택[110,111,112,113]에 대한 여러 검토가 있었으며 당시의 일부 모델에 대한 간략한 논의가 있었습니다[114]. 아주 최근에 Villena et al. [115]는 모든 RRAM 모델링의 이론을 결합하고 최적화 모델을 제안했습니다. 본 연구에서는 다양한 결점에 대한 솔루션과 함께 다양한 모델링 기법에 더 초점을 맞추었다. 유사 컴팩트 모델로 분류될 수 있는 경계 조건 모델에 대한 포괄적인 논의도 논의되었습니다. 이 작업에서는 모델 개발자에게 상당한 도움이 될 수 있는 몇 가지 중요한 모델링 기술이 조사되었습니다. 또한 SPICE[116, 117]와 같은 RRAM 모델에 대한 다양한 시뮬레이션 기술 및 플랫폼에 대한 논의가 포함되어 있어 매우 중요합니다. 우리의 작업은 RRAM 모델링 커뮤니티의 상당한 격차를 메우는 것을 목표로 합니다.

바이폴라 장치용 RRAM 모델

추아 모델

1971년 Leon O. Chua는 멤리스터[1]가 저항, 커패시터, 인덕터에 이어 네 번째 기본 요소라는 아이디어를 제시했습니다. 멤리스터의 기본 특성은 플럭스 제어(φ ) 또는 충전 제어(q ) 및 유형 g(φ,q)의 관계에 의해 정의됩니다. ) =0.

Chua는 멤리스터의 전압을 다음과 같이 정의했습니다. [1]:

$$ v(t)=M\left(q(t)i(t)\right) $$ (1)

어디에

$$ M(q)=d\varphi (q)/ dq $$ (2)

플럭스 제어 멤리스터를 통해 흐르는 전류는 다음과 같이 공식화되었습니다. 1 :

$$ i(t)=W\left(\varphi (t)v(t)\right) $$ (3)

어디에

$$ W\left(\varphi \right)=dq\left(\varphi \right)/ d\varphi $$ (4)

여기서 매개변수 M (q ) 및 W (φ )은 저항 및 컨덕턴스와 유사한 단위를 가지므로 각각 증분 멤리스턴스와 증분 멤덕턴스로 정의됩니다. φ-q 3개의 멤리스터 장치에 대한 곡선이 그림 1에 나와 있습니다. 이 곡선은 3가지 유형의 멤리스터를 생성하는 기본 멤리스터-저항(M-R) 회로에 의해 생성됩니다. φ-q 이러한 장치에 대한 분산은 각각 그림 1a-e에 나와 있습니다. 그림 1b–f는 해당 I-V를 나타냅니다. 동일한 3개의 멤리스터의 관계.

<그림>

f 플럭스 전하(ϕ -q ) 세 가지 다른 멤리스터에서 얻은 곡선 [1]

위에 제시된 방정식은 다음과 같이 단순화될 수 있습니다[1]:

$$ v=R(w)\times i $$ (5) $$ \frac{dw}{dt}=i $$ (6)

여기서 w 장치의 상태 변수이고 R 장치의 내부 상태에 따라 달라지는 일반화된 저항입니다.

t 시점에서의 증분 멤리스턴스(memductance) 값 0 t에서 완전한 멤리스터 전류(전압)의 시간 적분에 따라 달라집니다. =− t t = 0 . 따라서 이것은 memristor가 t 모든 순간에 일반 저항으로 작동한다는 사실로 해석됩니다. 0 , 그러나 저항(컨덕턴스) 값은 장치 전류(전압)의 완전한 과거 이력에 따라 달라지므로 메모리 저항이라는 이름이 정당화됩니다.

흥미롭게도 지정된 멤리스터 전압 v 시점에서 ( ) 또는 현재 i ( ), 멤리스터는 선형 시변 저항으로 동작합니다. 그러나 φ-q의 경우 곡선은 직선입니다. 즉, M (q ) =R 또는 W (φ ) =G , 멤리스터는 선형 시불변 저항기처럼 작동합니다. 따라서 멤리스터 장치는 선형 네트워크 이론에서 사용할 수 없지만 매개변수의 현재 상태가 과거 상태에 종속되는 회로를 정의하는 데 사용할 수 있습니다.

이후 1976년 Chua와 Kang[2]은 많은 비선형 동적 시스템을 포함하는 멤리스터 시스템을 포함하는 멤리스터 개념을 일반화했습니다. 방정식 [2]로 설명되었습니다.

$$ v=R\left(w,i\right)\times i $$ (7) $$ \frac{dw}{dt}=f\left(w,i\right) $$ (8)

여기서 w 상태 변수 세트로 정의됩니다. R 그리고 f 시간의 명시적 기능입니다. 멤리스터와 멤리스틱 시스템의 기본적인 차이점은 나중에 플럭스가 전하에 의해 더 이상 고유하게 정의되지 않는다는 것입니다. Memristive 시스템은 전압 강하가 0일 때 장치에 전류가 흐르지 않는다는 점에서 일반적인 동적 시스템과 구별될 수 있습니다.

멤리스터 방정식은 장치 모델링에서 멤리스터를 사용한 첫 번째 사례인 Chua[1]에 의해 임계값 스위치의 가변 상태를 정의하는 데 합리적으로 사용되었습니다. Chua에 의한 멤리스터의 공식화는 데이터를 저장하기 위해 기본 회로 요소를 사용하는 새로운 종류의 장치와 다양한 애플리케이션의 토대를 마련했습니다. 멤리스터의 이러한 기본 개념은 RRAM이 유망한 후보인 미래의 비휘발성 메모리 애플리케이션을 위한 새로운 아키텍처의 설계로 이어졌습니다. 기본적으로 멤리스터 모델을 기반으로 하는 RRAM 장치 및 이를 정의하는 모델의 작동을 설명하는 이론이 상당히 많습니다.

자속 전하 모델의 매우 흥미로운 적용은 유니폴라 RRAM을 정의하고 SPICE에서 구현하기 위해 [118] 사용하는 것입니다. 자속-전하 방정식의 단순성으로 인해 수정이 거의 없이 회로 시뮬레이터에 쉽게 통합될 수 있습니다. SPICE 모델은 HfO2의 실험 데이터에 대해 테스트되었습니다. 기반 단극 RRAM 장치. 실험적으로 얻은 정규화된 q에 맞추기 위해 제안된 비선형 관계 -φ 값은 [118]로 제공됩니다.

$$ q\left(\varphi \right)={q}_r\times \min \left(1,{\left(\frac{\varphi }{\varphi_r}\right)}^n\right) $$ (9)

여기서 φ r RESET 지점에서의 자속입니다. 이 값이 q일 때 (φ ) =q r 가 획득되면 CF가 사라지고 CF와 관련된 전류가 다시 0으로 설정됩니다. 이는 장치가 HRS에 있음을 의미합니다. 장치의 단극 스위칭 특성을 재현하는 모델의 능력을 조사하기 위해 표준 바이어스 스윕 작업이 수행됩니다. 리셋 상태에서 디바이스에 인가된 전압은 제로 바이어스에서 LRS에 도달할 때까지 점진적으로 증가하고 바이어스는 다시 0볼트로 스위프됩니다. LRS 전류는 [118]로 주어진 Chua 모델 [1]의 전류 관계의 수정된 형식을 사용하여 모델링됩니다.

$$ 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)

HRS 전류는 열전자 방출에 의해 제어되는 것으로 가정하므로 해당 상태의 전류는 다음과 같이 모델링됩니다.

$$ i(v)={I}_A\left({e}^{\frac{v}{v_A}}-1\right) $$ (11)

임계값 효과도 모델에서 고려됩니다. 문턱 전압 효과는 접촉 효과로 인해 발생한다고 가정했습니다. SET 및 RESET 프로세스 모두에서 자속 계산을 위한 전압 임계값을 포함하여 고려할 수 있습니다. 수정된 전류는 [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{배열 }\오른쪽. $$ (12)

여기, ϕ r ϕ s RESET 및 SET 플럭스는 각각입니다. 이러한 방정식은 커패시터 네트워크로 구성된 SPICE 호환 회로로 구현될 수 있습니다. SPICE 구현 결과는 모델이 거의 동일한 멤리스터 특성을 재현할 수 있는 실험 결과를 밀접하게 따르는 것으로 나타났습니다. 단극 장치를 모델링하는 데에도 사용되는 Chua 자속 전하 모델[1]의 사용을 검증합니다.

선형 이온 드리프트 모델

Chua가 멤리스터를 공식화한 후 수십 년 동안 상당한 격차가 있는 가운데 2008년 HP Labs[3]의 연구원들은 멤리스터 장치에 관한 흥미로운 발견을 했습니다. Chua가 멤리스터와 같은 요소의 존재를 공식화했지만 21세기 초반에 RRAM 장치를 제조하기 위한 여러 노력이 보고되었지만 이후에 개발된 실현 가능한 회로나 모델은 없었습니다. Strukov 등이 이끄는 HP Labs 팀. [3]은 고체 상태의 전자 및 이온 수송이 외부 전압 바이어스 하에서 함께 결합되는 멤리스턴스가 자연적으로 발생하는 기능적 나노 스케일 멤리스틱 시스템을 실현했습니다. 이러한 시스템은 다른 나노스케일 전자 장치와 유사한 전류 및 전압 특성 사이의 히스테리시스 관계를 보여 멤리스티브 시스템 및 유사한 시스템 설계에 대한 근본적인 이해로 이어집니다.

산화물(TiO2 ) 두께 D 두 개의 Pt 전극 사이에 끼워져 있습니다. 히스테리시스 -V 스위칭 곡선은 시뮬레이션된 곡선과 비교되었습니다. 그 당시에는 이러한 장치의 정확한 메커니즘이 완전히 이해되지 않았지만 저항성 스위칭 메모리가 멤리스틱 시스템으로 분류된 최초의 사례 중 하나였습니다.

TiO2의 개략적인 장치 구조 - 기반 멤리스터는 그림 2a[3]에 나와 있으며, 여기서 R이라고 하는 두 개의 가변 저항이 직렬로 연결되어 있습니다. 켜기 이는 도펀트 농도가 더 높은 반도체 영역의 낮은 저항입니다. 도펀트 농도가 낮을수록 R이라고 하는 다른 부분의 저항이 높아집니다. 꺼짐 . 인가 전압 v 사이의 관계 ( ) 및 시스템 i를 통한 전류 ( ) 옴 전자 전도도 및 평균 이온 이동도를 갖는 균일한 장에서의 선형 이온 드리프트로 인해 [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)

멤리스터에 대한 결합된 가변 저항 모델이 제시됩니다. (V) 전압계와 (A) 전류계로 구성된 단순화된 등가 회로. , 시간 t에 따른 인가 전압(파란색) 및 결과 전류(녹색) 전형적인 멤리스터에 대해서도 제시되어 있다. b에서 적용된 전압은 v입니다. 0 죄(v 0 ) 저항비는 R입니다. 꺼짐 /R 켜기 =160, c 적용된 전압은 ±v입니다. 0 20 ) 및 R 꺼짐 /R 켜기 =380, 여기서 ω0 는 주파수이고 v 0 인가 전압의 크기입니다. 숫자 1–6은 적용된 전압의 연속적인 파동과 i–v의 해당 루프에 대해 레이블이 지정됩니다. 곡선. 각 플롯에서 축은 무차원이며 전압, 전류, 시간, 플럭스 및 전하가 v 단위로 표시됩니다. 0 =1V, i 0 ≡ v0 /R 켜기 =10mA, t 0 ≡ 2π0 2v v 0 =10m/s, v 0 0 그리고 0 0 , 각각. 용어 i 0 장치를 통해 가능한 최대 전류를 나타내며 t 0 균일한 필드 v에서 전체 장치 길이에 걸친 도펀트의 선형 드리프트에 필요한 최단 시간입니다. 0 / , 예를 들어 D =10 nm 및 μV =10 −10 cm 2 s −1 V −1 . 선택된 매개변수에 대해 적용된 바이어스는 두 저항 영역 중 어느 하나를 강제로 무너뜨리지 않는다는 점에 유의해야 합니다. 예:w / 0 또는 1에 접근하지 않습니다(b의 중간 플롯에 점선으로 표시됨). 및 c ). 또한 파선 i–v b의 플롯 는 스위프 주파수가 10배 증가할 때 관찰되는 히스테리시스 붕괴를 보여줍니다. i–v 삽입 b의 플롯 및 c 이 예에서 전하는 멤리스터에 있어야 하므로 플럭스의 단일 값 함수임을 보여줍니다[3]

위의 방정식 자체는 비선형이지만 장치의 저항은인가 전압 v에 따라 선형적으로 변합니다. ( ), 따라서 모델에 대한 선형성의 귀속. Strukov et al.에 의해 정의된 장치. [3] 상태 변수 w의 특정 범위에 대해서만 완벽한 멤리스터 역할을 합니다. . 상태 변수는 [3]으로 정의됩니다.

$$ \frac{dw(t)}{dt}={\mu}_v\frac{R_{\mathrm{ON}}}{D}i(t) $$ (14)

Eq.에서 Chua[1]가 제안한 시스템의 Memristance (1)은 위의 두 식을 이용하여 정의된다. (13) 및 (14) [3]:

$$ M(q)={R}_{\mathrm{OFF}}\left(1-\frac{\mu_v{R}_{\mathrm{ON}}}{D^2}q(t)\ 오른쪽) $$ (15)

위 식에서 (15), q -의존 용어는 멤리스턴스에 대한 주요 기여입니다. 이 특정 현상이 왜 그렇게 오랫동안 숨겨져 있었는지에 대한 흥미로운 분석은 자기장이 메커니즘에서 명확한 역할을 하지 않았기 때문입니다. 멤리스터가 간단히 구현되기 위해서는 전압과 전류의 적분 사이에 비선형 관계가 존재해야 합니다.

에쿠스. (13)-(15)는 또한 양극 스위칭의 기본 사항을 통합합니다. 즉, 장치는 두 극성의 전압을 인가하여 한 상태에서 다른 상태로 전환합니다. 결과적으로 양극성 히스테리시스를 나타내는 장치 I -V 관계는 이러한 방정식에 의해 모델링될 수 있으며 따라서 이러한 장치를 멤리스티브 시스템으로 분류할 수 있습니다. 이러한 거동은 유기막[119,120,121,122,123], 칼코게나이드[124,125,126], 금속 산화물[127,128,129], 유전체 산화물[130,131,132], teamO 자체, 페로브스카이트[119,120,121,122,123], 2 [3] 시스템과 유사한 바이폴라 스위칭 특성을 관찰했으며, 활성 영역을 통한 불순물 또는 불순물 움직임이 저항의 이러한 극적인 변화의 원인으로 나타났습니다. 이것은 그림 2b, c에 나와 있으며, 전류는 전압의 변화에 ​​따라 급격한 감소와 급격한 상승을 보입니다.

물리적으로 이 두 단말 장치의 활성 영역은 0에서 D 범위 내에서 작동합니다. , 산화물 층의 두께, 그래서 상태 변수 w 또한 두께 사이에 경계가 있습니다. 그림 3은 w/D의 변화를 나타냅니다. 매개변수에 대한 시간이 0 및 D 범위를 벗어나지 않음 [삼]. 저항의 급격한 변화 또는 스위칭은 이러한 경계에 도달하는 장치로 인해 발생합니다. 이 조건을 모델링하기 위해 적절한 경계 조건이 사용됩니다. 특히 경계에서 장치에서 특정 이상 현상이 관찰됩니다. 사용 가능한 변경에 대해 동적 상태 변수의 비율에 일정하지 않은 변경이 있습니다. 또한 이온 이동도는 중간보다 경계에서 훨씬 적습니다. 이것은 경계에서 비선형 도펀트 드리프트 효과에 기인합니다. 따라서 이러한 효과를 적절하게 설명하기 위해 특정 창 기능의 변형을 사용하여 장치의 경계를 정의합니다. HP 팀은 상태 변수 Eq에 곱한 창 함수를 제안했습니다. (9) [3]으로 주어짐:

$$ f(x)=\raisebox{1ex}{$w\left(1-w\right)$}\!\left/ \!\raisebox{-1ex}{${D}^2$}\right . $$ (16) <그림>

시뮬레이션된 전압 구동 멤리스티브 장치. 동적 음의 차동 저항을 사용한 시뮬레이션. 동적 음의 차동 저항이 없는 시뮬레이션. 비선형 이온 드리프트에 의해 지배되는 시뮬레이션. a의 상단 플롯에서 , b , 및 c , 전압 자극(파란색) 및 정규화된 상태 변수의 해당 변화 w / (빨간색)은 시간에 따라 표시됩니다. 모든 경우에 하드 스위칭은 w일 때 발생합니다. / 0과 1(점선)의 경계에 가깝게 접근하고 질적으로 다른 i -v 히스테리시스 모양은 w의 특정 의존성으로 인한 것입니다. / 경계 근처의 전기장에서. d 비교를 위해 실험 i–v Pt–TiO2 − x의 플롯 -Pt 장치가 제시됨 [3]

이 모델은 미래의 RRAM 모델을 위한 토대를 마련한 데 기인할 수 있습니다. 양극성 히스테리시스 I를 갖는 2단자 반도체 장치에도 사용할 수 있습니다. -V 관계. 멤리스터의 메커니즘을 참조하여 RRAM 장치에 대한 수많은 미래 모델이 개발되었습니다.

비선형 이온 드리프트 모델

HP[3]에 의해 개발된 선형 이온 드리프트 모델은 주로 멤리스터 장치의 벌크 영역에서 선형 드리프트 효과를 보여주었습니다. 그들은 경계에서 일부 비선형 효과를 관찰했지만 포괄적으로 정의하지는 않았습니다. 인가 전압에 대한 도펀트 드리프트의 비선형 의존성은 Yang et al.에 의해 관찰되고 공식화되었습니다. [46] 2008. 그들은 비선형 효과를 정확하게 설명하는 전류-전압 관계를 제안했습니다. 나중에 Eero Lehtonen과 Mika Laiho가 개선하고 추가했습니다[68].

멤리스티브 장치의 전도는 공간적으로 이질적인 금속/산화물 전자 장벽에 의해 제어된다는 것이 Yang et al.에 의해 보고되었습니다. [46]. 스위칭은 이 전자 장벽을 통해 전도성 채널을 형성하거나 용해시키기 위해 네이티브 도펀트로 작용하는 양으로 하전된 산소 결손의 드리프트에 의해 발생합니다. 공석의 농도는 경계 또는 금속/산화물 계면에서 더 높습니다. ON 및 OFF 스위칭은 상단 인터페이스에서만 발생했으며 이는 상단 전극이 활성 전극으로 작용함을 나타냅니다.

산화티타늄 기반 멤리스터의 스위칭 특성에 대한 산소 결손의 영향은 그림 4[46]에 나와 있습니다. 극성에 의해 정의된 반대 스위칭을 보여줍니다. 또한 그림 4c와 같이 상단 인터페이스에 추가 공석을 추가하면 스위칭 곡선이 변경되어 멤리스틱 장치에서 비저항 인터페이스의 지배적인 역할이 확인됩니다. 이것은 인터페이스에서 발생하고 장치 스위칭을 제어하는 ​​비선형성 효과의 기초를 형성합니다.

<그림>

박막 TiO2 − x 제어된 산소 결핍 프로파일이 있는 장치는 스위칭 메커니즘을 확인하는 데 사용됩니다. 샘플 I 및 II는 15nm TiO2의 역층 시퀀스를 포함합니다. 및 15nm TiO2 − x (더 많은 공석) 레이어. 이것은 I-V의 반대 극성을 나타냅니다. 처녀 상태의 곡선. . 이 두 샘플의 스위칭 극성도 서로 반대입니다. . 이 두 샘플의 상단 인터페이스에 5nm Ti 레이어를 추가하여 더 많은 공석이 도입되면 I-V 곡선은 완전히 다른 방식으로 변화하여 박막 장치에서 비옴 인터페이스의 지배적인 역할을 확인합니다[46]

Yang et al. [46]은 멤리스티브 장치가 인가된 전류 또는 전압의 시간 적분에 따라 상태를 변경하는 동적 저항기 역할을 한다는 사실을 설명했습니다. 그들은 동적 상태 변수를 설명하는 관계를 제공하지 못했습니다. 제안된 전류-전압 관계는 [46]으로 설명될 수 있습니다.

$$ I={w}^n\베타 \sinh \left(\alpha v\right)+\chi \left({e}^{\gamma v}-1\right) $$ (17)

여기서 β, γ, n , 및 χ는 피팅 상수입니다. 위의 방정식에서 첫 번째 항 β sinh(αv )은 [1] 전자가 얇은 잔류 전자 장벽을 통해 터널링되는 멤리스터의 ON 상태에 가깝습니다. 0(OFF) ~ 1(ON) 범위에서 디바이스의 상태 변수로 정의됩니다. 방정식의 두 번째 부분은 피팅 상수로 작용하는 다른 매개변수를 사용하여 장치의 OFF 상태를 근사화합니다. 매개변수 n 여기에서 상태 간 전환을 수정하는 데 사용되는 자유 매개변수 역할을 합니다. n을 조정하는 동안 비선형 효과가 나타납니다. -V 제작된 장치의 곡선은 식을 사용하여 모델링됩니다. (16). 14 ≤ n에서 최적의 피팅을 얻습니다. ≤ 22. 이는 효과적인 공석 드리프트 속도가 장치에 인가된 전압에 따라 매우 비선형적으로 의존한다는 증거로 해석될 수 있습니다. 따라서 경계/인터페이스에서 대부분의 도펀트 드리프트 효과는 본질적으로 비선형으로 이해될 수 있습니다.

상태 변수 w의 역학을 설명하는 관계 SPICE[116, 117]를 사용하는 이 모델은 Lehtonen과 Laiho[68]에 의해 제안되었습니다. w의 시간 도함수 [68]로 모델링되었습니다:

$$ \frac{dw}{dt}=a\times f(w)\times g(v) $$ (18)

여기, a 상수, f :[0, 1] → R 는 제안된 창 함수이고 g:R → R은 선형 드리프트 모델에서 이전에 제안된 선형 함수로 간주됩니다(여기서 R 실수)를 나타냅니다. 저자는 Yang et al.이 제안한 멤리스터의 작동을 모방하기 위해 솔루션에서 시연했습니다. [46], (v ) must be a non-linear, odd, and monotonically increasing function. A non-linear function which was proposed was [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 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 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. , d The numerical derivative w ˙ of the data in ac 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 그리고 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 그리고 , 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)

어디에

$$ {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, 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 ( ), 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 ( )) and f (x ( )), are responsible for the change in the state variable. 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 ( )) 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 그리고 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 ( ( )) 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 p ) is a window function which limits the value of f ( ) to 0 when x ( ) = 1 and v ( ) > 0. f n ( ,w n ) is a similar window function which does not allow the value of w ( ) 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. -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]. R의 값 ONR OFF are set as 50 Ω and 1 kΩ, and an ideal rectangular window function is applied in Eqs. (38) and (37). -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 켜기 그리고 꺼짐 act as the current thresholds. Functions f 켜기 and f 꺼짐 are window functions which bound the internal state variable x ( ) within [w 켜기 , w 꺼짐 ]. 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.

-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.

-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 켜기 and f 꺼짐 act as window functions which bound the internal state variable w within [w 켜기 , w 꺼짐 ]. 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]. Yakopcic model [73]. BCM model [99]. 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 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)

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. -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].

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. 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 O5 )-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)

where 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. V density (n D ) map. Calculated profiles of b n D , , 및 d for states A (1.0 V), B (1.7 V), and C (3.0 V). The position of z  = 15 nm indicates the Ta2 O5 /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 ). 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−V . 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].

-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)

, 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 α 그리고 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 그리고 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 . 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)

여기, 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 h are the activation energies for oxygen species generation and hopping respectively. Similarly, a 그리고 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)

여기, d mn is the distance between m 그리고 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 (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 ). 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 그리고 α 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.

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; -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.

Applied triangular voltage signal (blue) and corresponding device current (red) versus time. The RRAM current against the applied voltage. 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 영화. Schematic of C-AFM measurements (b ) Typical I -V curves of TiO2 films using C-AFM tip as a top electrode. At V tip = 8 V, mapping of the current flow through the surface just after the forming operationshows locally distributed conducting regions. d TiO2 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 is the activation energy, k 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)

여기서 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); 0.3 V, f 0.375 V, g 0.45 V, and h 0.475 V for sample Q (400 °C annealed); 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 = 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 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.

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 작업. 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 . 여기, 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 ( ) = 2v0 sin(ω 0 /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 2D where 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 ( ) did not go beyond its limits. The comparators clipped w ( ) 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 ( ) reached the upper boundary of the device, the device is connected to a voltage source of value D , which effectively clips w ( ) at D . This operation occurs when the voltage bias is positive.

At negative voltage bias when w ( ) approaches the lower boundary, the circuit is connected to ground, thus clipping w ( ) 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 MINR MAX . Thus, in mathematical terms it means that R ∈ [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 최대 (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 최대  > 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. 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 ( ) ∈ (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)

and

$$ 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 ( ) there exists exactly one solution. Values of the non-negative parameters v th,1v th,0 determine the occurrence of such transitions. The conditions in the first Eq. (125) are established when x ( ) obeys the boundary conditions x  = 1 and x  = 0. But the conditions in the second Eq. (126) is established when the function x ( ) 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. -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 ( ) and current i ( ). A state variable s ( ) 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 ( ) at a particular time instant t is governed by the history or the state of v ( ). 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 ( ) 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 clipminF 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 w1F 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.


나노물질

  1. Elatec:USB 포트가 없는 장치에 대한 사용자 인증 및 액세스 제어 활성화
  2. SARS-CoV-2에 대한 새로운 통찰력을 제안하는 모델링 및 시뮬레이션
  3. 원자층 증착에 의해 제조된 Pt 및 TiN 코팅 기판 상의 HfO2/TiO2/HfO2 삼중층 구조 RRAM 장치의 양극성 저항 스위칭 특성
  4. 3상 전기 변색 장치를 위한 침지 코팅 공정 엔지니어링 및 성능 최적화
  5. 세륨 다이아몬드 절단의 분자 역학 모델링 및 시뮬레이션
  6. 보로펜의 안정성 및 STM 이미지에 대한 제1원칙 연구
  7. HT29 및 SPEV 세포주에 대한 Au 나노입자의 영향에 대한 체외 연구
  8. Abaqus의 금속 재료 모델링
  9. 안전 장치 및 고려 사항
  10. 3D 모델링 및 시뮬레이션의 힘으로 제조 공정 혁신 촉진