Российские нанотехнологии, 2024, T. 19, № 2, стр. 178-182

ВВЕДЕНИЕ

Физический механизм резистивного переключения в нитридных мемристорах (НМ) – эволюция проводящих филаментов в их рабочей зоне – представляет собой процесс диффузии вакансий азота в AlN под действием электрического поля и градиента температуры в системе [1]. Согласно [2] оптимальные условия функционирования НM достигаются в результате сложного компромисса между составом, структурой и геометрией гетероструктуры и ограничениями на диапазоны внешних условий. Исчерпывающие реальные исследования подобных многопараметрических задач по понятным причинам не только весьма затратны, но и невыполнимы практически. Поэтому единственным продуктивным сценарием является их сочетание с численным экспериментом, основанным на понимании физических процессов в системе, формализованных в математическую модель, которая алгоритмизована и воплощена в соответствующей программе. В рассматриваемом случае такой инструмент должен обеспечить самосогласованное описание процессов эволюции распределения в системе электростатического потенциала, температуры и плотности вакансий под действием переменного внешнего смещения и при определенных граничных условиях на все эти поля. Для процессов переноса, как известно, оптимальным по соотношению качества результатов, скорости счета и требований к вычислительным ресурсам, является гидродинамическое приближение. Но для того чтобы применить его к описанию исследуемой системы, необходимо знать тензор диффузии вакансий азота в AlN.

В [3] из первых принципов рассчитаны высоты активационных барьеров при перескоках вакансий азота на соседние узлы кристаллической решетки AlN в зависимости от их заряда и направления движения в рабочей области рассматриваемых здесь НМ – слое кристаллического нитрида алюминия со структурой вюрцита. В табл. 1 приведены величины активационных барьеров при диффузии вакансий азота в кристаллическом AlN со структурой вюрцита в зависимости от их заряда и направления перескока.

В настоящей работе на основе полученных данных рассчитан тензор диффузии $\hat {D}$ вакансий азота в этой системе в отсутствие внешнего электрического поля. Для анализа зависимости его компонент от локальной температуры среды расчеты проведены для трех температур – 300, 500 и 1000 К. Как будет видно из результатов, этого массива данных вполне достаточно для ее характеризации. При этом рассмотрен лишь случай вакансий $V_{{\text{N}}}^{q}$ с зарядом $q = - 1$e, которые, согласно табл. 1, в основном определяют кинетику филаментов в НМ.

МЕТОД ИССЛЕДОВАНИЯ

Рассмотрим диффузию заряженных вакансий азота в кристаллическом AlN со структурой вюрцита как термически активированное скачкообразное движение по узлам N кристаллической решетки. Если все ближайшие узлы решетки заполнены, возможны 12 направлений прыжков заряженных вакансий азота: шесть из них расположены в плоскости (0001) кристалла, еще по три направления ориентированы вверх и вниз относительности плоскости (0001) попарно симметрично относительно нее. Хотя рассматриваемая система имеет объемную симметрию $C_{{6{v}}}^{4}$, возможные пути диффузии заряженных вакансий $V_{{\text{N}}}^{q}$ (рис. 1) характеризуются тригональной симметрией ${{C}_{{3{v}}}}$.

В предположении, что перенос вакансий в системе обусловлен только градиентом их концентрации n(r, t), компоненты соответствующей плотности потока частиц j(r, t) подчиняются закону Фика:

где $\alpha ,\beta = (x,y,z)$, ${{D}_{{\alpha \beta }}} \equiv {{\left( {\hat {D}} \right)}_{{\alpha \beta }}}$, r – радиус-вектор физически бесконечно малой области системы, в которой концентрация вакансий азота в момент времени $t$ есть n(r, t). Выберем ортогональную систему координат, в которой оси $Ox$ и $Oy$ лежат в плоскости (0001), а ось $Oz$ ориентирована вдоль направления [0001] (рис. 1). В силу симметрии системы тензор диффузии в этой системе координат имеет диагональный вид, ${{D}_{{\alpha \beta }}} = {{\delta }_{{\alpha \beta }}}{{D}_{{\alpha \alpha }}}$, причем ${{D}_{{xx}}} = {{D}_{{yy}}} \ne {{D}_{{zz}}}$. Предполагая, что число ионов азота (а значит, и количество вакансий $V_{{\text{N}}}^{q}$) в системе сохраняется, подставим вектор плотности потока вакансий из (1) в уравнение непрерывности ${{\partial n} \mathord{\left/ {\vphantom {{\partial n} {\partial t}}} \right. \kern-0em} {\partial t}} + \operatorname{div} {\mathbf{j}} = 0$, и получим для $n({\mathbf{r}},t)$ уравнение диффузии в главных осях тензора диффузии: где подразумевается суммирование по индексу $\alpha $.

Предположим, что тензор диффузии в (2) зависит явно только от концентрации диффузанта: ${{D}_{{\alpha \alpha }}} = {{D}_{{\alpha \alpha }}}\left( {n({\mathbf{r}},t)} \right)$, а профиль концентрации диффундирующих частиц в начальный момент времени зависит лишь от координаты вдоль одного из его главных направлений – пусть для определенности ${{\left. {n({\mathbf{r}},t)} \right|}_{{t = o}}} = n(x,0)$. Тогда уравнение (2) становится одномерным, а подстановка Больцмана ${{\left( {x - {{X}_{M}}} \right)} \mathord{\left/ {\vphantom {{\left( {x - {{X}_{M}}} \right)} {2\sqrt t }}} \right. \kern-0em} {2\sqrt t }} = \xi $, где ${{X}_{M}}$ – координата некой фиксированной точки оси $Ox$, позволяет свести его к обыкновенному нелинейному дифференциальному уравнению вида

На этом свойстве уравнения (2) основан классический метод Матано–Больцмана [4], который позволяет решить обратную задачу и определить неизвестную компоненту тензора диффузии по известным из эксперимента значениям $n(x,t)$ для так называемой диффузионной пары – пространственно-однородной системы толщиной ${{L}_{ \prec }} + {{L}_{ \succ }}$ с одномерным начальным распределением концентрации диффундирующих частиц вида $n(x,0)$ = = $n{}_{ \prec }\theta \left( {{{L}_{ \prec }} - x} \right)$ + $n{}_{ \succ }\theta \left( {x - {{L}_{ \succ }}} \right)$ и при условии, что в ходе эксперимента они не успевают достичь ее внешних границ, т.е. что ${{\left. {{{\partial n(x,t)} \mathord{\left/ {\vphantom {{\partial n(x,t)} {\partial x}}} \right. \kern-0em} {\partial x}}} \right|}_{{{{L}_{ \pm }}}}} = 0.$ В этом случае любая плоскость заданной концентрации перемещается в эксперименте со скоростью ~$\sqrt t $. При этом параметр ${{X}_{M}}$ выбирается из условия, чтобы поток диффундирующих частиц через плоскость $x = {{X}_{M}}$ (называемую плоскостью Матано) был равен нулю: ${{\left. {{\mathbf{j}}({\mathbf{r}},{{t}_{0}})} \right|}_{{x = {{X}_{M}}}}}$ ≡ ${{j}_{x}}(x,{{t}_{0}}){{{\mathbf{e}}}_{x}}$ = 0, где ${{t}_{0}}$ – момент измерения, ${{{\mathbf{e}}}_{x}}$ – единичный вектор оси $Ox$.

Соответствующее дискретизированное решение обратной задачи диффузии имеет вид, например, [5]:

где ${{n}_{{\max }}} = \max \left( {n({{x}_{1}}),...,n({{x}_{N}})} \right)$, n_min = $\min \left( {n({{x}_{1}})} \right.$, ..., $\left. {n({{x}_{N}})} \right)$, $N$ – количество точек расчетной сетки на оси $Ox$ и предполагается, что и положение плоскости Матано, и все $n({{x}_{i}})$, $i = 1,...,N$ отвечают одному и тому же моменту времени $t = {{t}_{0}}$.

Отметим, что последнее условие ставит под сомнение саму возможность осуществления реального эксперимента в такой постановке, но не представляет трудности для его численного моделирования, пусть и сопряженного с определенными методическими ограничениями. Действительно, нахождение массива концентраций $\left\{ {n({{x}_{i}},{{t}_{0}})} \right\}$ заряженных вакансий методом молекулярной динамики затруднено, так как временной шаг молекулярной динамики ~10^–15 c, тогда как в исследуемой системе имеют место высокобарьерные события с характерными временами ~10^–9 c и более (табл. 1). Поэтому для ее моделирования существенно более эффективен кинетический метод Монте-Карло [6], который численно описывает процесс хаотического движения системы через последовательность квазиравновесных состояний, отвечающих реализации наименьших энергетических барьеров вдоль траектории конфигурационного пространства, ведущей к состоянию термодинамического равновесия. Здесь рассмотрены процессы диффузии вакансий азота при фиксированных температурах в сверхъячейке кристаллического нитрида алюминия со структурой вюрцита, состоящей из 50 × 50 × 200 элементарных ячеек со структурой вюрцита, и реализован алгоритм моделирования канонического ансамбля атомов, подробно описанный, к примеру, в [6].

РЕЗУЛЬТАТЫ И ИХ ОБСУЖДЕНИЕ

Некоторые результаты моделирования процесса диффузии вакансий в главных направлениях тензора диффузии при температуре 1000 К представлены на рис. 2. Начальная высота концентрационной ступеньки в этих расчетах составляла от 1 до 100% вакансий азота на элементарную ячейку нитрида алюминия с шагом 9%. При этом 100% начальной высоты концентрационной ступеньки вакансий соответствует структуре, состоящей из 175 × 10⁴ атомов и 25 × 10⁴ вакансий с зарядом –1е.

Результаты расчета концентрационной зависимости главных компонент тензора диффузии на основе формул (4), (5) приведены на рис. 3, где сплошной линией изображен результат подгонки зависимостью

а $\alpha = x,z$. Величины параметров подгонки даны в табл. 2, 3.

Как видно, параметры подгонки существенным образом зависят от температуры, и дальнейшее применение рассмотренного метода расчета компонент тензора диффузии к решению конкретных задач потребует более детального анализа этих зависимостей.

Согласно результатам выполненного расчета диффузия вакансий азота в кристаллическом нитриде алюминия со структурой вюрцита имеет существенно анизотропный характер, так что в отсутствие внешнего электрического поля перенос вакансий в исследованном здесь интервале температур 300–1000 К происходит преимущественно в плоскости (0001). Нетрудно видеть, что в данном интервале температур вклад в диффузию перескоков, не лежащих в этой плоскости, ничтожен, так что ${{D}_{{xx}}} = {{D}_{{yy}}}$ ~ $\exp ({{ - \hbar {{\omega }_{{{\text{jamp}}}}}} \mathord{\left/ {\vphantom {{ - \hbar {{\omega }_{{{\text{jamp}}}}}} {{{k}_{{\text{B}}}}T}}} \right. \kern-0em} {{{k}_{{\text{B}}}}T}})$, где частота перескоков ${{\omega }_{{{\text{jamp}}}}} \approx 1.43$ эВ (табл. 1).

ЗАКЛЮЧЕНИЕ

В контексте мемристивных применений указанный выше факт означает, что в рабочем режиме НМ перенос вакансий через границы пленки AlN, выращенной вдоль плоскости (0001), практически отсутствует, но может иметь место под действием импульса формовки. Кинетика происходящих при этом процессов требует специального рассмотрения, однако априори ясно, что оптимальная степень нестехиометричности рабочей области мемристора по азоту является существенным параметром, выбор которого зависит, в частности, и от меры участия электродов в ее вакансионном балансе.

Отметим, что полученные результаты позволяют восстановить и зависимости компонент тензора диффузии вакансий от градиента потенциала, что дает возможность применить их для моделирования процессов формирования и эволюции вакансионных филаментов в исследуемой системе в рамках самосогласованного гидродинамического приближения.

Работа выполнена при финансовой поддержке Российского научного фонда, проект № 22-29-00535.

Численное моделирование выполнено на ресурсах Федерального центра коллективного пользования научным оборудованием “Комплекс моделирования и обработки данных исследовательских установок мега-класса” НИЦ “Курчатовский институт”.