Журнал вычислительной математики и математической физики, 2021, T. 61, № 5, стр. 885-894

1. ВВЕДЕНИЕ

Одним из центральных вопросов, возникающих при разработке аналоговых и цифровых сверхпроводниковых микросхем, является создание численных алгоритмов, решающих задачу вычисления распределения сверхпроводящих токов в таких структурах, а также значений собственных и взаимных индуктивностей их отдельных частей. Это связано с тем, что работоспособность и производительность устройств на основе сверхпроводниковых квантовых интерферометров (СКВИДов) находятся в критической зависимости от индуктивностей отдельных компонентов таких устройств.

Расчет реальных индуктивностей особенно важен для устройств, преобразующих магнитный поток в ток или напряжение. Примером таких технических решений являются сверхпроводниковые антенны на основе би-СКВИДов и Д-СКВИДов. В этом случае преобразование потока в напряжение должно оставаться линейным с высокой точностью без использования цепей обратной связи [1]–[3].

В настоящее время особый интерес представляют нелинейные элементы сверхпроводниковых искусственных нейронных сетей с магнитным представлением информации. В этом случае преобразование магнитного потока в ток (функция активации) для нейрона должна с хорошей точностью совпадать с сигмоидальной функцией, либо с гиперболическим тангенсом [4].

Целью данной работы является разработка усовершенствованного вычислительного алгоритма и программы экстракции (расчета) рабочих и паразитных индуктивных коэффициентов сверхпроводниковых микросхем, в частности методики расчета индуктивностей для СКВИДов и устройств на основе СКВИДов. Эффективность разработанного подхода показана на примере одного из таких устройств, искусственного нейрона с сигмовидной активационной функцией [5]–[7]. В отличие от полупроводниковых аналогов сверхпроводниковые нейроны оказываются существенно менее сложными, и могут быть относительно просто изготовлены с использованием существующей микроэлектронной технологии. Однако их работоспособность требует высокой точности реализации нужного соотношения между значениями индуктивностей элементов схемы. В данной работе мы демонстрируем возможности и пределы применимости нашего вычислительного алгоритма и программы для вычисления индуктивностей [8] и [9] на примере решения актуальной проблемы — проектирования упомянутого сверхпроводникового нейрона.

Математическая модель, лежащая в основе расчетов пространственного распределения токов в сверхпроводниковых микроэлектронных структурах, основана на совместном решении уравнений Лондонов и Максвелла [10]. В стационарном случае такой подход приводит к интегральному уравнению относительно векторного потенциала [11], содержащему трехмерные интегралы.

Наш подход содержит ряд модификаций этого интегрального уравнения. Во-первых, используя то обстоятельство, что в рассматриваемых структурах толщина слоев невелика по сравнению с лондоновской глубиной проникновения магнитного поля в сверхпроводник, мы преобразуем объемные токи в листовые токи. Во-вторых, мы учитываем конечную толщину слоев в функции Грина, которая является ядром интегрального уравнения. В-третьих, для точного и простого учета межслойных соединений и джозефсоновских переходов мы разделяем ток на потенциальную и вихревую составляющие и представляем вихревую часть тока через функцию тока. Это сводит проблему к решению гиперсингулярного интегрального уравнения с хорошими математическими свойствами, что значительно упрощает численное решение. При таком подходе точно учитываются как магнитная, так и кинетическая части индуктивности фрагментов микросхем.

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

Вычислительный алгоритм реализован в нашей программе 3D-MLSI [9].

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

2. МОДЕЛЬ НЕЙРОНА

Искусственные свехпроводниковые нейроны [5]–[7] могут быть спроектированы на базе структур с двумя сверхпроводящими слоями над сверхпроводящим же экраном. Результат проектирования сверхпроводящего нейрона в виде упрощенной схемы показан на фиг. 1. Он включает в себя три разделенных диэлектриками металлических слоя: плоскость заземления M0, слой разводки M1 и слой с верхним покрытием M2. Верхний слой М2 в разных точках характеризуется разными расстояниями до подложки (разными высотами) и частично замыкается на М1, если М1 и М2 перекрываются. В наших расчетах мы используем упрощенную модель, показанную на фиг. 2. В этой модели слой M2 имеет постоянную высоту. Замыкания слоев M1–M2 представлены в виде пары токовых контактов прямоугольной формы.

При последующих расчетах мы будем считать, что лондоновская глубина проникновения $\lambda $ [10] для всех слоев одинакова и составляет $0.075$ мкм. Нижний слой M0 является базовой заземляющей плоскостью. Он собирает все возвратные токи из более высоких слоев. Толщину этого слоя будем считать равной $0.3$ мкм ($4\lambda $). Следующий слой – M1. Он не является сплошным и содержит несколько участков, по которым протекает ток. Полагаем, что его высота относительно верхней части нижнего слоя составляет $0.3$ мкм, а толщина – $0.11$ мкм (менее $2\lambda $). Верхний слой M2 находится на высоте $0.61$ мкм и имеет толщину $0.45$ мкм. Между верхним и нижним слоями имеется три джозефсоновских контакта и множество межслойных контактов. Оставшиеся после изготовления нейрона участки слоев представлены на фиг. 2. Средний и верхний слои отдельно представлены на фиг. 3. Там же приведена и конечно-элементная сетка. Ширина полос в М2 составляет $10$ мкм. Верхний и нижний слои довольно толстые по сравнению с $\lambda $. Это ограничивает точность расчетов в нашей математической модели на уровне 5% [8].

Физическая постановка задачи наряду с описанием проводящей структуры должна содержать определение токов, для которых будут вычислены индуктивности. Эти токи, как правило, образуют токовые петли, замыкающиеся через заземляющий слой и проходящие через определенные контакты, в том числе через джозефсоновские переходы. Последовательность таких контактов определяет ток, для которого должны быть вычислены индуктивные коэффициенты. Два таких токовых контура представлены на фиг. 2. В свою очередь заданные процессы в контактах позволяют формулировать краевые и иные условия для постановки математической задачи. В результате каждый токовый контур дает одну строку (столбец) в матрице индуктивности. Диагональные элементы этой матрицы являются собственными индуктивностями, а внедиагональные элементы являются взаимными индуктивностями.

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

3. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ДЛЯ ТОКОВ И ИНДУКТИВНОСТЕЙ

Вычислительный алгоритм этой работы является улучшенной и усовершенствованной версией алгоритмов, предложенных в [8] и [9].

Введем некоторые обозначения.

Пусть ${{N}_{c}}$ – число односвязных проводников. Проводник расположен в некотором слое металлизации. Каждый проводник занимает трехмерную область:

здесь ${{S}_{m}}$ есть 2D-проекция ${{V}_{m}}$ на плоскость.

Пусть $\partial {{S}_{{ext,m}}}$ – граница проводника $m$, исключая границы возможно существующих отверстий. Толщина проводника $m$ составляет ${{t}_{m}} = h_{m}^{1} - h_{m}^{0}$. Пусть ${{N}_{h}}$ – общее количество всех отверстий во всех проводниках. Будем называть $\partial {{S}_{{h,k}}}$ границей отверстия $k$.

Всего имеется ${{N}_{T}}$ токовых контактов (терминалов) во всех проводниках, не менее двух терминалов для каждого проводника. Джозефсоновские переходы образуют пару терминалов на разных проводниках. Терминалы на $\partial {{S}_{{ext,m}}}$ являются внешними терминалами. Терминалы внутри ${{S}_{m}}$, например, в форме прямоугольников, являются внутренними терминалами. Джозефсоновский переход образует два внутренних терминала.

Обозначим через ${{N}_{t}}$ число замкнутых или открытых путей для интересующих нас полных токов. Каждый путь содержит некоторую цепочку терминалов. Ток может переходить от проводника к проводнику, а также между слоями через терминалы. Всего имеем $N = {{N}_{t}} + {{N}_{h}}$ независимых полных токов, так как существуют замкнутые токи вокруг отверстий.

В качестве результата нас интересует вычисление $N \times N$ матрицы индуктивности $L$ [8].

Математическая модель основана на стационарных уравнениях Лондонов и формулы Био–Савара для магнитного поля ${\mathbf{H}}$, а также выражения для полной энергии $E$:

где $\lambda $ является лондоновской глубиной проникновения магнитного поля в сверхпроводник [10], [11]. Из этих выражений следует объемное интегральное уравнение для фазы ${{\varphi }_{n}}(r)$ сверхпроводящего параметра порядка и плотности тока ${{j}_{n}}(r)$: Неизвестными являются фаза и плотность тока. Объемное интегральное уравнение, дополненное условием $\nabla {{{\mathbf{j}}}_{n}} = 0$, является основой численных методов для расчета токов ${{{\mathbf{j}}}_{n}}(r)$, а затем индуктивностей [12]–[15]. В нашем подходе мы используем некоторые преобразования (2), позволяющие ввести токи вокруг отверстий и существенно упростить численное решение.

Существует ряд специфических трудностей в моделировании сверхпроводящих токов. Первой трудностью является сложная топология распределения тока в объемной структуре. Известно [10], [11], что уравнение Лондонов приводит к концентрации тока вблизи поверхности сверхпроводников. Точный учет этого эффекта довольно сложен алгоритмически и приводит к продолжительным вычислениям, однако при некоторых условиях (тонкий проводник) его общее влияние на значения индуктивности является умеренным.

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

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

В нашем подходе мы используем возможность упростить задачу для случая сверхпроводящих слоев относительно небольшой, но конечной толщины. Обычно толщина сверхпроводящих слоев составляет порядка 1–3 лондоновской глубины проникновения. В этом случае для проводника $m$ мы можем использовать модель листового тока ${{{\mathbf{J}}}_{m}}(x,y) = ({{J}_{{m,x}}}(x,y),{{J}_{{m,y}}}(x,y))$. Полагая ${{j}_{z}} \approx 0$, получаем

Конечная толщина проводящих слоев эффективно учитывается путем усреднения уравнения (2) по толщине слоя и введения некоторой простой функции Грина вместо потенциала простого слоя с ядром ${1 \mathord{\left/ {\vphantom {1 {\left| {r - r{\text{'}}{\kern 1pt} } \right|}}} \right. \kern-0em} {\left| {r - r{\text{'}}{\kern 1pt} } \right|}}$, $r = (x,y)$:

В этом выражении $\lambda _{m}^{s} = {{\lambda }^{2}}{\text{/}}{{t}_{m}}$ – лондоновская глубина проникновения для тонких сверхпроводящих пленок.

Затем с помощью оператора ротора ${{\nabla }_{ \bot }} = ({{\partial }_{y}}, - {{\partial }_{x}})$ интегральное уравнение для листовых токов преобразуется в гиперсингулярное интегральное уравнение. Гиперсингулярная форма уравнения предпочтительнее для численного решения, так как его интегральный оператор является симметричным и положительно-определенным.

Более того, в (2) мы исключаем член $\nabla {{\varphi }_{n}}(r)$. Для этого плотность тока сверхпроводимости представляется в виде суммы двух компонент:

Полагаем Будем полагать ${\mathbf{J}}_{m}^{{{\text{pot}}}}(r)$ терминальным током возбуждения, ${\mathbf{J}}_{m}^{{{\text{rot}}}}(r)$ представляет собой ток экранирования магнитного поля и ток вокруг отверстий.

Окончательно приходим к следующему интегральному уравнению для функции тока ${{\psi }_{n}}(r)$, $n = 1,\; \ldots ,\;{{N}_{c}}$:

Пусть ${{I}_{{h,k}}}$ – полный ток, циркулирующий вокруг отверстия $k$. Тогда граничными условиями для функции тока являются [8]: Здесь ${{C}_{{h,k}}}$ – двумерная граница дырки $k$. В настоящее время отверстия в физическом дизайне рассматриваемых моделей нейронов отсутствуют. Далее, ${{C}_{{ext,m}}}$ – внешняя граница проводника $m$. Подробно краевые условия для функции тока представлены в [8]. Если в схеме нет отверстий, то на всех границах ${{\psi }_{m}}(r) = 0$.

Для потенциала ${{\varphi }_{n}}(r)$ в (8) ставится вторая краевая задача

где ${\mathbf{n}}$ – нормаль. Функция ${{f}_{n}}(r)$ определяет модель терминала как распределенный источник или сток. Пусть ${{I}_{k}}$ – полный ток через терминал $k$ с областью ${{T}_{k}} \in {{S}_{m}}$ и площадью $\left| {{{T}_{k}}} \right|$, тогда: Форма терминалов может быть произвольной. Для контактов большой площади можно моделировать лондоновскую неоднородность тока с помощью кольцевых терминалов [9].

Таким образом, для достижения общей цели необходимо последовательно решить ряд задач. Сначала решаются ${{N}_{c}}$ задач (10), (11) с заданным ${{I}_{k}}$ для ${\mathbf{J}}_{m}^{{{\text{pot}}}}(r)$. После этого решаются задачи (7)–(9). Наконец, полная энергия рассчитывается с использованием выражения:

Матрица индуктивности содержит все индуктивные коэффициенты и рассчитывается с помощью выражения для полной энергии [8], [9].

Для однослойной структуры ядра ${1 \mathord{\left/ {\vphantom {1 {\left| {r - r{\text{'}}{\kern 1pt} } \right|}}} \right. \kern-0em} {\left| {r - r{\text{'}}{\kern 1pt} } \right|}}$ и одного проводника уравнение для функции тока имеет вид

Для оператора в смысле главного значения имеем В [16] доказано, что для функций, принимающих нулевое значение на границе, оператор $L$ положителен.

Формула (15) важна для доказательства свойств интегрального оператора и корректности метода конечных элементов. Кроме того, представление интегрального оператора в гиперсингулярной форме может быть использовано для быстрой оценки и вычисления матричных элементов в методе конечных элементов для удаленных друг от друга точек сетки (см. [17]).

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

4. ВЫЧИСЛИТЕЛЬНЫЙ АЛГОРИТМ И ПРОГРАММА

Для решения обеих краевых задач (10), (11) и (7), (8) используется метод конечных элементов. Мы используем треугольные сетки и линейные лагранжевы конечные элементы [18]. Такой выбор приводит к простому и быстрому алгоритму сборки общей матрицы метода конечных элементов.

Для уравнения Пуассона (10), (11) используется аналогичный метод конечных элементов. Этот алгоритм описан в литературе.

Для краткости рассмотрим задачу с одним проводником. Тогда билинейная форма задачи для функции тока имеет вид

Сеточная аппроксимация функции тока с главными краевыми условиями имеет вид

где суммирование проводится по всем точкам сетки $J$, включая граничные точки, $\psi _{j}^{h}$ – приближенные значения функции тока в точках сетки, $u_{j}^{h}(r)$ – базисные функции интерполяции метода конечных элементов.

В нашем вычислительном алгоритме используются линейные лагранжевы конечные элементы (P1). Эти элементы обладают достаточной гладкостью и удобны для аппроксимации гиперсингулярного оператора (15).

Такой метод конечных элементов приводит к системе линейных уравнений для $\psi _{i}^{h}$, $i \in I$, $I$ – внутренние точки сетки:

Матрица в (18) плотная, симметричная и при достаточно малом шаге сетки положительно-определенная. Размерность этой матрицы равна числу внутренних точек сетки ${{N}_{I}}$. Граничные значения $\psi (r)$ формируют правую часть как главное краевое условие.

Таким образом, матричные элементы конечно-элементной матрицы $K$ имеют вид

где ${{S}_{i}}$, ${{S}_{j}}$ – носители базисных функций конечно-элементной аппроксимации $u_{i}^{h}$, $u_{j}^{h}$. Эффективный метод вычисления этих элементов описан в [8] и [9], [17].

Обоснование описанного метода конечных элементов для гиперсингулярного уравнения с ядром ${1 \mathord{\left/ {\vphantom {1 {{{{\left| {r - r{\text{'}}{\kern 1pt} } \right|}}^{3}}}}} \right. \kern-0em} {{{{\left| {r - r{\text{'}}{\kern 1pt} } \right|}}^{3}}}}$ было дано в [16].

Для решения линейных уравнений описанного метода конечных элементов мы используем разложение Холецкого для плотных матриц. Метод конечных элементов для (10), (11) основан на применении разреженных матриц.

Для построения треугольных сеток наша программа содержит собственный генератор. Также можно использовать известный генератор треугольных сеток Triangle [19].

В алгоритме есть несколько вычислительных процедур, занимающих много времени процессора. Первой процедурой является вычисление матрицы $K$ (19) и правой части $F(r)$ в (8). Второй процедурой является расчет полной энергии, определяемой выражением (13). Чтобы ускорить эти вычисления, мы вводим матрицу Галеркина $T$ взаимодействий между треугольниками на сетке метода конечных элементов, которая является общей частью обоих вычислений. Пусть ${{\Delta }_{i}}$ и ${{\Delta }_{j}}$ – треугольные ячейки, тогда элементы матрицы взаимодействия $T$ являются четырехкратными интегралами:

Таким образом, вычислительный алгоритм основан на двух матрицах, $T$ и $K$. Матрица $T$ имеет размерность, равную числу ячеек в сетке ${{N}_{t}}$. Вторая матрица, $K$, является матрицей системы линейных уравнений метода конечных элементов, размерность этой матрицы – число внутренних узлов сетки ${{N}_{I}}$. ${{N}_{I}}$ и ${{N}_{t}}$ в практически интересных случаях могут быть большими (несколько тысяч) или очень большими (десятки тысяч). В настоящее время мы вычисляем эти плотно заполненные матрицы и храним в явной форме. В данной работе мы демонстрируем решение одной практической задачи с помощью этого подхода. В то же время мы рассматриваем разработанный вычислительный подход как основу для разработки и реализации алгоритмов скелетонного типа и малоранговой аппроксимации для работы с плотно заполненными матрицами большой размерности.

Программная реализация численного алгоритма содержит вычислительные модули, написанные на C++ и работающие в Linux или Windows. Также нами разработан графический Windows интерфейс для визуализации входных данных и распределения тока [9].

5. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Рассмотрим модель нейрона для многослойного персептрона, см. фиг. 2. Исходные данные для расчетов первоначально были подготовлены в САПР в формате DXC (DXF), а затем в полуавтоматическом режиме преобразованы во входные данные для программы 3D-MLSI.

Схема содержит множество возможных токовых структур, которые можно использовать для прямого или косвенного вычисления схемных и паразитных индуктивностей. Для краткости рассмотрим пример с двумя замкнутыми токами. Во-первых, это горизонтальная замкнутая токовая петля с собственной индуктивностью ${{L}_{{11}}}$. Эта петля моделирует измерительный СКВИД с двумя джозефсоновскими контактами и возвратным током через заземленную поверхность. Второй контур с собственной индуктивностью ${{L}_{{22}}}$ – это вертикальная (на фигуре) петля с одним джозефсоновским переходом и обратным током через плоскость заземления. Обе петли показаны на фиг. 3.

Для расчета описанных индуктивностей с точностью 5–6% решалась задача с размерностями матриц ${{N}_{t}} = 12\,722$ и ${{N}_{I}} = 5770$. Всего система уравнений с матрицей $K$ решалась 3 раза. Были получены значения ${{L}_{{11}}} = 9.05$, ${{L}_{{12}}} = 1.8$ и ${{L}_{{22}}} = 26.3$ пикогенри. Линии тока для расчета ${{L}_{{11}}}$ представлены на фиг. 4. Ток через первый контур создает близкие к прямым линии тока. Токи экранирования на других частях конструкции образуют вихревые структуры. Полученные результаты были далее использованы при моделировании процессов в таком нейроне, а также для дальнейшей оптимизации его структуры.

Описанные вычисления оказались достаточно быстрыми и заняли несколько минут на процессоре Intel i7 в одноядерном режиме. Те же вычисления с более плотной сеткой и ${{N}_{I}} = 18\,694$ и ${{N}_{t}} = 42\,842$ заняли всю 16-гигабайтную память и потребовали несколько часов счета.

5. ВЫВОДЫ

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

Вычислительный метод настоящей работы, программа и методология расчетов могут быть довольно просто использованы для вычисления индуктивностей схем небольших и умеренных размеров. Для моделирования больших схем требуется дополнить программу средствами автоматической подготовки входных данных, использующими данные систем автоматизированного проектирования и разработки (САПР) микросхем.

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

Авторы благодарят В. Больгинова за плодотворные обсуждения возможных реализаций моделей сверхпроводниковых нейронов.