CC BY
234
УДК 537.533:004.94:519.63 Вестник СПбГУ. Сер. 10, 2013, вып. 2
К. А. Никифоров, Н. В. Егоров
МОДЕЛИРОВАНИЕ ЭМИССИОННЫХ ПРОЦЕССОВ В СРЕДЕ MATLAB
1. Введение. Физическая модель. В работе рассматриваются многоэмиттерные полевые эмиссионные катоды различного типа и оценивается их эффективность с позиций математического и компьютерного моделирования в среде МАТЬАВ. Полевой эмиссионный катод является одним из перспективных источников электронов для многих практически важных приложений [1, 2].
Изучаются два типа многоэмиттерных систем: с острийными коническими и лезвийными прямолинейными эмиттерами. Методами математического и компьютерного моделирования исследуется взаимодействие между эмиттерами в массиве и анализируется влияние геометрических параметров на величину тока полевой электронной эмиссии. Характеризуются также такие геометрические параметры как расстояние между эмиттером и коллектором (¿), радиус кривизны эмиссионной поверхности (г), расстояние между эмиттерами (Ь) (рис. 1). Зависимость плотности тока полевой электронной эмиссии от напряженности электрического поля, вызывающего ее, рассчитывается в рамках теории Фаулера-Нордгейма [1]. Многоэмиттерная система рассматривается как периодически структурированный массив эмиттеров. Вершина конического эмиттера аппроксимируется сферической поверхностью, а край лезвийного - цилиндрической (рис. 1). Неоднородность физических условий на эмиссионной поверхности (шероховатость на атомарном масштабе, изменения работы выхода и пр.) учитывается в коэффициенте усиления локального поля [1]. Предполагается, что наноструктуриро-ванные детали поверхности одинаковы по всей обследуемой области, а коэффициент усиления постоянен. Влиянием пространственного заряда пренебрегается и используется электростатическое приближение.
2. Математическая модель. Постановка задачи. Основная цель - моделирование электрических полей и токов в ячейке эмиттера, рассматриваемой в двумерной геометрии.
Никифоров Константин Аркадьевич — кандидат физико-математических наук, доцент, Санкт-Петербургский государственный университет; e-mail: nikiforov_k@mail.ru.
Егоров Николай Васильевич — доктор физико-математических наук, профессор, Санкт-Петербургский государственный университет; e-mail: robby7@mail.ru.
© К. А. Никифоров, Н. В. Егоров, 2013
7
Рис. 1. Схематическое изображение конических и лезвийных эмиттеров
Показаны ячейки для моделирования с условиями симметрии (периодичности) на боковых поверхностях.
В соответствии с принятой физической моделью процессы эмиссии описываются следующими уравнениями:
^Р = 0, (1)
E = -V?, (2)
(3)
S
В (1)—(3) p - электростатический потенциал, E - вектор напряженности электрического поля, I - эмиссионный ток, S - эмиссионная поверхность, A и B - постоянные, связанные с фундаментальными физическими величинами согласно соотношениям A = e3/(8n h), В = 8тгл/2^/(3eh), h - постоянная Планка, e - заряд электрона, m -масса покоя свободного электрона, t(y), v(y) - эллиптические функции Нордгейма аргумента у = \/е3Е/(47Г£оФ2), для которых используются аппроксимации t2(y) « 1.11619, v(y) « 0.95 — y2, £q - электрическая постоянная, E - величина напряженности внешнего электрического поля, Ф - работа выхода материала, являющаяся мерой энергии связи электронов с твердым телом.
Уравнение (1) решается в замкнутой области G, соответствующей геометрии ячейки, на границе Г которой заданы краевые условия для p: 1) на границах проводников (электродов)
(p\ri = (i = const,
2) на границе ячейки (т. е. на боковой поверхности, на линиях или плоскостях симметрии) - условие равенства нулю нормальной производной от потенциала (нормальной составляющей электрического поля)
-О
3) на границе диэлектриков - условие непрерывности нормальной производной потенциала
где п - нормаль к поверхности границы; ^ - потенциал г-го проводника с границей Г®; Га - граница сред с различными диэлектрическими проницаемостями е\ и £2 (по разные стороны от границы); Ге - граница ячейки вакуумного микроприбора.
Ввиду осевой симметрии для ячейки с коническим эмиттером применяется цилиндрическая система координат, в которой уравнение Лапласа (1) переписывается в следующем виде (с учетом осевой симметрии задачи д/дв = 0, где 6 - угловая координата):
1 дер д2ср д2ср р др др2 дх2
В случае прямолинейных лезвийных эмиттеров в уравнении (1) также используется двумерная постановка (д/дх = 0, ось г направлена вдоль лезвия эмиттера), координаты декартовы:
дх2 ду2
3. Решение задачи средствами MATLAB и MATLAB PDE Toolbox. Численное решение задачи основано на методе конечных элементов с реализацией в пакете PDE Toolbox в среде MATLAB [3-5].
Рис. 2. Графический интерфейс пользователя: основное рабочее окно с редактором параметров (слева) и окно визуализации адаптивной сетки (справа) На вставке внизу: детализация сетки вблизи эмиссионной поверхности.
Программная реализация в виде отдельно стоящего (stand alone) приложения с графическим интерфейсом пользователя (рис. 2) учитывает специфику эмиссионных систем:
• вычислительная область сложной формы включает границу эмиттера с большой кривизной поверхности и малыми размерами, что приводит к значительному разбросу характерных размеров в одной геометрической конфигурации;
• экспоненциальная зависимость плотности тока от напряженности поля требует повышенной точности при учете граничных условий на эмиттере.
Решение имеет быстро изменяющийся градиент в области эмиссии (на вершине катода), поэтому конечноэлементная сетка должна сгущаться в окрестности вершины эмиттера, чтобы скорость сходимости решения к точному существенно не снизилась и не произошло увеличения числа неизвестных - размерности конечноэлементной системы.
При адаптивном построении сетки используется индикатор ошибки, включающий норму невязки уравнения и скачки градиента конечноэлементного решения, поскольку они связаны с одной из основных в данной задаче физических величин - напряженностью электрического поля. Этот индикатор реализован в функции pdejmps, которая, в свою очередь, вызывается в функции adaptmesh, предназначенной для адаптивного решения задачи. При измельчении сетки выбран способ разбиения конечных элементов
по наибольшей стороне (параметр rmethod=longest), при котором взамен подлежащего разбиению треугольника (на основе значения индикатора ошибки) образуются два новых треугольника делением наибольшей его стороны пополам.
Другие параметры adaptmesh, управляющие ходом адаптивных вычислений, могут изменяться, в том числе и для учебных целей, и задаются пользователем с помощью графического интерфейса: число получающихся треугольников (конечных элементов) -параметр тах^ которое может быть и неограниченным (тахЬ=1п£); выбор треугольников для разбиения на очередном шаге измельчения сетки (параметр ^1р1ск) может производиться либо по относительной погрешности (реализуется функцией pdeadgsc), при котором разбиению подлежат конечные элементы со значением индикатора ошибки, большем заданной погрешности с масштабирующим множителем, либо по условному критерию качества (функция pdeadworst), когда разбиваются треугольники со значением индикатора ошибки, превышающим заданную часть этого показателя для худшего из всех треугольников (задается коэффициентом отношения между 0 и 1); максимальное количество итераций измельчения сетки (параметр ngen), которое может быть и неограниченным (ngen=inf).
В общем случае, если начальная сетка не задана, adaptmesh вызывает функцию initmesh для инициализации сетки, при этом используются значения по умолчанию для параметров initmesh. Важнейшие из параметров initmesh - скорость роста размера треугольников (параметр hgrad со значением между 1 и 2), характеризующая степень неравномерности начальной сетки, и максимальный размер стороны треугольников (параметр hmax). Чтобы иметь возможность изменять эти параметры, при создании сетки функция initmesh вызывается явно, и получившаяся триангуляция передается для уточнения в функцию adaptmesh. Параметры задаются также через графический интерфейс.
Функция adaptmesh возвращает массивы p, e, ^ содержащие построенную триангуляцию ^ - координаты вершин, e - информацию о граничных ребрах триангуляции, t - инцидентную матрицу), а также массив и с двумерным конечноэлементным решением краевой задачи для дифференциального уравнения в частных производных эллиптического типа:
—V • (е^и) + аи = /,
здесь входные параметры - скалярные величины е, а, / задаются в соответствии с поставленной задачей (уравнение Лапласа): а = / = 0. В цилиндрических координатах (уравнение (4)) е = £х, в декартовых (уравнение (5)) е = £, где £ - относительная диэлектрическая проницаемость среды, х - первая координата двумерного пространства.
Чтобы построенное методом конечных элементов решение уравнений (4), (5) соответствовало указанным в п. 2 граничным условиям задачи, в функцию adaptmesh передаются массив dl, содержащий матрицу декомпозиционной геометрии с информацией о граничных участках вычислительной области, и массив Ь, включающий матрицу граничных условий.
Моделирование вольт-амперных характеристик эмиссионных структур требует многократного запуска алгоритма в цикле по изменяющимся значениям электродных напряжений. При этом изменения вносятся только в массив Ь, в соответствии с новыми величинами напряжений, а конечноэлементная сетка-триангуляция заново не перестраивается. Функция adaptmesh уже не вызывается каждый раз в цикле, а производится сборка матрицы и вектора правой части конечноэлементной системы линейных алгебраических уравнений при помощи функции assempde, которая и решает данную
систему. Такая организация вычислений значительно снижает время счета без существенных потерь точности, при условии, что разброс значений напряжений вдоль вольт-амперной характеристики невелик.
На основе рассчитанных узловых значений электростатического потенциала определяется распределение напряженности поля при помощи функций pdegrad (вычисляет градиент конечноэлементного решения в средних точках треугольников) и pdeprtni (интерполирует значения в узлы сетки). Картина поля показана на рис. 3.
а б
2, МКМ
р, мкм
Рис. 3. Картина эквипотенциальных (а) и силовых (б) линий электрического поля (в цилиндрических координатах) в ячейке конического эмиттера (напряжение - 30 В, шаг эквипотенциалей - 2 В)
В соответствии с блок-схемой алгоритма после определения потенциала и напряженности поля производится расчет значений эмиссионного тока I в граничных узлах вычислительной области по формуле (3). Интеграл находится по квадратурным формулам Гаусса-Лобатто с использованием функции quadl.
4. Результаты моделирования. Приведем примеры расчетов, выполненные с помощью разработанного программного комплекса. В них использовались следующие параметры: а = 54.7° (это значение угла конуса получается в результате анизотропного травления кремния [6]), высота эмиттеров Н = 4 мкм, работа выхода Ф = 4.01 эВ (сильно легированный кремний п-типа), коэффициент усиления локального поля 7 = 7.5.
4-1- Распределение потенциала и напряженности поля. Расстояние между соседними эмиттерами - это один из параметров, определяющих взаимодействие между эмиттерами и их электрическое поле. Как показано на рис. 4, а, напряженность электрического поля на вершинах (лезвиях) уменьшается с приближением эмиттеров друг к другу, так как взаимодействие между эмиттерами усиливается.
Влияние расстояния с! на электрическое поле также проиллюстрировано на рис. 4, а. С уменьшением расстояния между эмиттером и коллектором напряженность поля на вершинах эмиттеров все меньше зависит от расстояния между ними. Следовательно, ослабевает эффект взаимной экранировки и эмиттеры начинают работать как индивидуальные острия (лезвия). Изменение с! имеет больший эффект у конических эмиттеров, нежели у лезвийных.
Е, ГВ/м
4.5 4
3.5
3
2.5 2
1.5 1
0.5
—*—4 ♦ -ж
Л" -г-' *
А
» ♦ * ♦ ♦ ♦ ♦ -Н т
Е, В/м
II
4 8 12 16 20 24 Ъ, мкм
Рис. 4- Зависимость электрического поля на вершине эмиттера от межэмиттерного расстояния (6) и расстояния эмиттер-коллектор (й) при постоянном радиусе г = 10 нм (а) и величина напряженности электрического поля на поверхности от вершины эмиттера до середины между ними при напряжении 5 В (б)
Семейство кривых I соответствует коническим эмиттерам, II — лезвийным. о: напряжение 30 В, на вставке указаны значения d в мкм; б — по оси абсцисс отложена горизонтальная координата точек поверхности.
На рис. 4, б показано электрическое поле эмиттера, построенное в точках поверхности начиная от вершины до середины между соседними остриями (лезвиями). На графике видны резкое увеличение напряженности вблизи вершины эмиттера, благодаря так называемому «эффекту острия», и минимум у его основания из-за «впадины», образованной боковой поверхностью эмиттера и плоскостью подложки.
4.2. Эмиссионный ток. На рис. 5, а показано распределение нормализованной плотности тока вдоль поверхности эмиттера. Как видно, при удалении от вершины острия (края лезвия) плотность тока спадает очень быстро и основная часть эмиссионного тока производится очень малой областью вблизи вершины эмиттера. Таким образом, эмиссионный ток сильно локализован вблизи острия.
Другая важная деталь показана на рис. 5, б: увеличение тока эмиттера в связи с уменьшением расстояния между эмиттером и коллектором (с!). Когда с! уменьшается до величины порядка радиуса кривизны острия, ток лезвийного эмиттера начинает превышать ток с отрийного (рис. 5, б). Такое поведение может быть обьяснено тем, что при уменьшении с! плотность тока принимает сопоставимые по порядку величины для обоих случаев и количество тока (3) становится больше в случае лезвийного эмиттера из-за большей площади эмиссии. Важная особенность, видная на рис. 5, б, состоит
а
б
log10(y/y0),*1000
/, А
-0.16
-0.08
-0.04
-0.12
0
0 0.05 0.1 0.15 0.2 0
0.05 0.1 0.15 0.2 0.25
d, мкм
х, мкм
Рис. 5. Плотность эмиссионного тока на поверхности эмиттера от вершины к основанию (а) и зависимость тока эмиссии от расстояния эмиттер-коллектор при r = 0.01 мкм, напряжение 5 В (б)
jo — плотность тока на вершине; линия I соответствует острийному эмиттеру, II — лезвийному.
в том, что скорость изменения эмиссионного тока в зависимости от d для лезвийных эмиттеров выше, чем для острийных, что опять связано с большей площадью эмиссии. Отсюда следует, что лезвийные эмиттеры имеют более высокую чувствительность к межэлектродному расстоянию, чем острийные.
Это свойство, очевидно, справедливо для любой формы лезвий (не только прямолинейных [7]), что в практическом применении имеет большое значение для некоторых сенсорных устройств, например датчиков давления с коллектором в виде мембраны [5].
5. Заключение. Представлены физическая и математическая модели, результаты моделирования потенциала, напряженности электрического поля и эмиссионного тока многоэмиттерных систем в программном комплексе, разработанном в среде MATLAB и MATLAB PDE Toolbox. Описаны алгоритмы и функции метода конечных элементов на неравномерной сетке из PDE Toolbox, примененные в программной реализации.
Работа выполнена с использованием программного обеспечения и вычислительного оборудования Ресурсного центра «Вычислительный центр» СПбГУ.
Литература
1. Егоров Н. В., Шешин Е. П. Автоэлектронная эмиссия. Принципы и приборы. М.: Издат. Дом «Интеллект», 2011. 704 с.
2. Денисов В. П., Егоров Н. В., Карпов А. Г. Модели транспорта электронов в твердых телах. СПб.: Изд-во С.-Петерб. ун-та, 2011. 204 с.
3. Ануфриев И. Е. Применение PDE TOOLBOX при изучении некоторых разделов вычислительной математики // Труды III Всерос. науч. конференции «Проектирование инженерных и научных приложений в среде MATLAB». СПб., 2007. С. 43-54.
4. Никифоров К. А., Егоров Н. В. Научно-учебный программный комплекс для конечноэлементно-го моделирования диодных и триодных структур вакуумной микро/наноэлектроники // Труды V Меж-дунар. науч. конференции «Проектирование инженерных и научных приложений в среде MATLAB». Харьков, 2011. С. 659-679.
5. Nikiforov K. A., Egorov N. V. Program Complex for Vacuum Nanoelectronics Finite Element Simulations // Proc. of XXIII Russian Particle Accelerator Conference. St. Petersburg, 2012. P. 409-411.
6. Lee H. C., Huang R. S. A Theoretical Study of Field Emission Array for Microsensors // IEEE Transactions on electron devices. 1992. Vol. 39. P. 313—324.
7. Nikiforov K. A., Antonova L. I., Egorov N. V. e. a. Non-gated Field Emission Array as Low-Energy Electron Source: Experiment and Simulation // Proc. of XXIII Russian Particle Accelerator Conference. St. Petersburg, 2012. P. 218-220.
Статья поступила в редакцию 20 декабря 2012 г.
