CC BY
25
Вычислительные технологии
Том 15, № 6, 2010
Конечно-разностный метод вычисления основного состояния наноструктур в шестизонном
k • р-приближении*
В. П. Жуков, М.П. Федорук Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: zukov@ict.nsc.ru, mife@ict.nsc.ru
Предложен компактный, дивергентный и эрмитовый конечно-разностный аналог оператора Шредингера для шестизонной к • рмодели. Описан неявный итерационный метод поиска основного и первого возбужденного состояний этого конечно-разностного оператора. Показано, что обобщение шестизонной к • рмодели на случай зависящих от пространственных координат эффективных масс имеет различный вид. Спектр некоторых вариантов такого обобщения нефизичен.
Ключевые слова: шестизонная к • рмодель, квантовая точка, оператор Шредингера, эффективная масса, неявная конечно-разностная схема, собственная функция.
Введение
В настоящее время успешно выращиваются полупроводниковые наноструктуры, представляющие собой включение германия в виде правильной четырехгранной пирамиды в матрицу кремния (квантовая точка) [1]. Эти структуры имеют различные приложения. В частности, обсуждается возможность их применения в качестве элементной базы компьютеров, в том числе квантовых компьютеров, открывающих принципиально новые возможности в вычислениях за счет использования квантово-механичееких эффектов [2]. В связи с этим встает вопрос о свойствах наноструктур, при изучении которых возникает потребность в вычислении энергетического спектра и волновых функций, соответствующих определенному энергетическому состоянию. При этом с практической точки зрения наибольший интерес представляют основное и первое возбужденные состояния. В настоящей работе описывается численный метод, позволяющий проводить вычисления этих состояний в приближении хорошо зарекомендовавшей себя шестизонной к • рмодели [3, 4].
Особенностью квантовых точек пирамидальной формы является низкая (по сравнению со сферически или цилиндрически симметричными структурами) симметрия. При решении квантово-механичееких задач с низкой симметрией применение разложений уравнения Шредингера по каким-либо в достаточной мере простым функциям [5, 6] нерационально. В этом случае более удобен конечно-разностный метод. Отметим, что использованный в работах [3, 5] простейший итерационный метод для поиска собственных функций и собственных значений конечно-разностного аналога оператора Шредин-
*Работа выполнена при поддержке Интеграционного проекта СО РАН № 43 и РФФИ (грант № 08-01-00264-а).
гера требует неприемлемо большого расчетного времени, а для алгоритма, предложенного в [7, 8], в котором применяется разложение по большому числу построенных определенным образом сеточных функций, необходима значительная память компьютера, В настоящей работе предложены конечно-разностный аналог оператора Шрединге-ра и итерационный метод, позволяющий эффективно вычислять энергию и волновую функцию основного и первого возбужденного состояний наноструктуры, состоящей из кремниевой матрицы с включением германия, описываемых шеетизонной к • р-моделью в трехмерном случае, на персональных компьютерах за разумное время. Представленный конечно-разностный оператор Шредингера является эрмитовым, дивергентным (что важно для задач с разрывными коэффициентами) и компактным (для аппроксимации используются только ближайшие узлы).
Исходно к • р-модель выведена для структур с постоянными эффективными массами [3, 4], Значения этих масс для германия и кремния различаются в несколько раз, в связи с чем возникает необходимость обобщить данную модель на случай изменяющихся в пространстве эффективных масс. Это можно сделать различными способами [3, 7], В настоящей работе вопросу о виде оператора Шредингера в шеетизонной к • р-модели в случае непостоянных масс уделено значительное внимание.
Результаты расчетов и аналитические исследования, представленные в работе, привели к выводу о необходимости введения определенных ограничений при использовании шеетизонной к • р-модели для структур, состоящих из веществ с изменяющимися в пространстве эффективными массами, в частности, в случае широко распространенной пары германий/кремний. При численном моделировании учет этих ограничений требуется только при расчетах наноструктур, но строго говоря ограничения имеют место при любых (в том числе и очень больших) характерных размерах полупроводниковой структуры при наличии резкой границы между составляющими ее веществами.
1. Постановка задачи
Одноэлектронная волновая функция Ф, описывающая распределение электронов в системе, состоящей из кремниевой матрицы с германиевыми квантовыми точками, представляет собой вектор из шести компонент {ф\,'ф2,'ф3,'ф4,'фъ,'ф6)- Стационарное уравнение Шредингера для нее имеет вид [3, 4]
Е Ф = ЯФ =
Н, 0 0 Н,
Ф
£ 0 0 £
Ф
Ао
3
0 —г 0 0 0 1
г 0 0 0 0 —г
0 0 0 —1 г 0
0 0 —1 0 г 0
0 0 —г —г 0 0
1 г 0 0 0 0 )
Ф + и Ф. (1)
Здесь в первых двух членах справа под "О" понимается нулевая матрица 3 х 3,
где столбцы hft, h^, матрицы Hß имеют вид
( д_
дх
К1
L
д
дх
д
д
д
дх
+ ду у ду д
ду,
+
д
dz
M
д_\ \
dz
/'— | V— ) +(1
д
ду \ дх
V (
д д д д ^\NTz К*
h<2' =
д
дх
д
ду
\
д ( д \ д (м д \ д ( д ду \ ду) дх I дх) dz 1 dz
V /
h<3)
дд ß— I N-
^dz I dx
дд дд
\
^ dz ( dy) ^ ^ ^ ду ^ 'dz
KnI"
V
9 I L 9 I 9 (M 9 I 9 (M 9 dz \ dz J dx \ dx J dy l ду _
/
Коэффициенты Ь, М N £гг, и, До являются вещественными функциями координат, при этом Ь, М, N Д0 ^ 0 Коэффициенты Ь, М, N связаны с массой частицы в различных направлениях, е связана с деформационными эффектами, Д0 — со спин-орбитальным взаимодействием, и — добавка к потенциальной энергии, обеспечивающая разрыв зон на гетерогранице германий—кремний. Роль коэффициента ц (0 < ^ < 1) будет объяснена в следующем разделе.
Задача состоит в отыскании собственных функций и соответствующих им уровней энергии Е в предположении, что на бесконечности Ф = 0, Поскольку представляющие интерес локализованные волновые функции достаточно быстро стремятся к нулю с
Ф=0
области 0 < х < х0, 0 < у < у0, 0 < z < z0.
Проведем нормировку задачи. Выберем в качестве масштаба длины величину I в нм, времени — г, = 2ш!2/П = 1.7275555 ■ 10-14(ш/ше)/2гат] в с. Здесь П = 1.0546 ■ 10-34 Дж ■ с — постоянная Планка, ш — масса носителя, ше = 9.1094 ■ 10-31
кг
масса
LM
N необходимо вычислять в единицах h2/(2me), а под До, <% и U надо понимать эти величины, выраженные в электронвольтах (1 эВ = 1.6022 • 10-19 Дж) и умноженные на
u0 = (2m/2/h2)1.6022 • 10-19 = 26.24587(m/me)/fnm],
m/me
m/me = 1
2. Проблемы k • p-модели и выбор параметра ц
Наличие коэффициента ц в (2) связано с попыткой обобщения k • p-модели на случай непостоянных в пространстве коэффициентов L, M и N. Как легко видеть, при N = const параметр ц в (2) исчезает, Строгого вывода уравнений в случае N = const не существует. Из требования естественной симметрии уравнений относительно поворотов
N = N
взять линейную комбинацию выражений
которые при N = ^^^ те отличаются друг от друга. Под —& понимаются векторы (—4,—5,—б) в зависимости от того, па какие компоненты Ф данная часть оператора действует. Оба оператора (3), (4) эрмитовы. Выражение (3) соответствует ц = 1, а (4) — ц = 0,
Исследования выявили следующий недостаток шеетизонной к • р-модели в случае непостоянных коэффициентов. Существует диапазон значений параметра ц, определяющийся величинами Ь, М и N входящих в гетероетруктуру материалов, в котором спектр оператора Гамильтона Н является нефизическим, А именно, качественная за-Е
собственной функции оператора Н, от характерного пространственного масштаба собственной функции Я имеет вид, изображенный на рис, 1, При малых Я-1 энергия растет с увеличением Я-1, однако начиная с некоторого значения она уменынаетея. При Я — 0 существует асимптотика
Волновая функция при Я — 0 оказывается локализованной около границы контакта двух материалов с разными параметрами (в нашем случае на гранях пирамиды) в слое толщиной ~ Я и имеет характерный масштаб осцнлляцпй вдоль контактной грани-Я
являетея нефизической. Характерная величина пространственного масштаба Я*, при Е
наибольшим пространственным масштабом, имеет порядок нанометра.
Для случая контакта германия и кремния подобное явление наблюдается при ц* <
Ц < 1, Ц* ~ 0.4, а при ц, близком к пулю, оно не существует. Однако это не озна-
ц
ли взять пару гипотетических материалов с Ь1 = М1 = (Ьое + Мое)/2, N1 = N06
и Ь2 = М2 = (Ь$1 + Мз!)/2 N2 = N3;, то оператор (1), (2) будет иметь указанный де-цц
эффект отсутствует. Если взять два материала с параметрами германия и кремния,
но в качестве использовать величину, значительно меньшую ее реального значе-
ц
К изложенным выводам приводит анализ результатов аналитических решений за-
НН
(3)
(4)
1
Е ~ —, Е < 0.
(германия), помещенного в другой материал (кремний), и задачи о контакте полупространств, занимаемых разными материалами. Решение данных задач будет приведено в отдельной работе. Полученные результаты полностью подтверждаются приведенным ниже численным решением задачи о квантовой точке германия в матрице кремния.
Наличие пефизического решения в какой-то степени можно попять из следующих соображений. Легко показать, что собственные значения оператора — Нм при постоянных Ь, М и N для собственных фу нкций Ф ~ егк(х+у) имеют в ид Е = (Ь+М—N )к2, Пр и N < Ь + М это соответствует случаю свободной частицы с массой порядка (Ь + М — N)-1 > 0, при N > Ь + М — случаю отрицательной массы и минимальной энергии, равной —то при к = то. При контакте таких материалов как германий и кремний на контактной границе величина ~ 34.14 существенно превышает Ьд; + Мд; « 9.44. Таким образом, имеет место контакт двух материалов с N > Ь + М, Заметим, что для типичных полупроводников (в том числе германия и кремния) 0 < Ь + М — N < Ь + М, т. е. ситуация близка к предельно допустимому случаю Ь + М = N.
Рассмотрим проявление пефизических эффектов при решении задач па собственные значения численными методами. Следует отметить, что использование разностных методов ограничивает наименьший масштаб задачи шагом разностной сетки к, а при использовании разложений по каким-либо функциям — характерным масштабом наиболее сильно меняющейся в пространстве функции. Поэтому при поиске основного состояния оператора Н со значениями ц, дающими нефизический эффект в случае структур с масштабом выше нескольких десятков нанометров, в качестве основной будет выделяться гармоника с наименьшим пространственным масштабом. Уменьшение шага сетки к будет создавать иллюзию сходимости разностного решения к этой гармонике. Однако, как
< Я*
границе мелкомасштабные гармоники и численное решение будет носить явно пефи-зический характер. Поскольку многомерные задачи требуют больших вычислительных ресурсов, то возможна ситуация, когда вычисления производятся но па достаточно мелкой сетке, дающей правильный с точки зрения соответствия дифференциальной задаче (по но физике!) результат. Аналогичное может произойти и при решении рассматривае-
Ф
Рис. 1. Качественная зависимость энергии Е связанного состояния от характерного пространственного масштаба К волновой функции этого состояния при неудачном выборе параметра у
мало, то нефизическое решение может себя не проявить. Указанные особенности численного решения подтверждаются решением нижеприведенных задач.
При больших характерных размерах задачи ограниченность мелкого масштаба шагом разностной сетки является своеобразным фильтром, регуляризатором для некорректной исходной дифференциальной задачи. Однако для систем с размерами порядка нескольких нанометров шаг сетки сравним с критическим масштабом и требуется выбор конкретного значения
Заметим также, что уравнения более сложной восьмизонной модели [4] имеют структуру, подобную структуре уравнений рассматриваемой шестизонной модели, поэтому при описании структур с меняющимися в пространстве эффективными массами в данной модели могут возникать аналогичные проблемы.
3. Конечно-разностный аналог оператора Шредингера
Простейший конечно-разностный аналог оператора Шредингера (1), в котором все компоненты волновой функции вычисляются в одних и тех же точках при зависящем от координат коэффициенте М, оказывается не эрмитовым. Для исправления этой ситуации будем вычислять различные компоненты волновой функции и задавать различные коэффициенты в разных узлах: -01, —4, ё11 в точках г + 1/2, ], к, —2, —5, е22 в точках г,] + 1/2, к, —з, —6, е33 в точках г,], к + 1/2, Коэффициенты Ь, и, До и недиагональные элементы ё (е^ при г = ]) будем задавать в точках г,], к, а М — в точках г + 1/2,] + 1/2, к + 1/2, При вычислении членов с сомножителем ц в (2) функция N вычисляется в точках г + 1/2,] + 1/2, к + 1/2, а при вычислении членов с сомножителем 1 — ц — в точках г,], к. Поскольку N является заданной функцией, вычисление этого параметра в разных точках не представляет трудности.
Вычисление различных величин на сдвинутых сетках позволяет не только естественным и компактным образом аппроксимировать оператор И^, но и обеспечивает эрмито-вость конечно-разностного аналога этого оператора. Например, для первой компоненты волновой функции будем иметь (для краткости рассмотрим случай ^ = 1)
д^ д^ ду у ду дх у дг у дг дх ^ ^ ^ _ Ухг+1 ,],к — Чхг,],к % ¿+1/2,7+1/2,к ~ % ¿+1/2,7 — 1/2,к Ог+1/2,7,£+1/2 ~ & ¿+1/2,7,*:-1/2
кх ку кг
г "1
Qxí,j,k к 7 )
кх
_ Mi+1/2,j+1/2,k+1/2 + Mi+1/2,j+1/2,k-1/2 "1 i+1/2,j+1,k — "1 i+1/2,j,k Чу ¿+1/2,7+1/2,А: - -^--^--Ь
Д+1/2,7+1/2, к+1/2 + Д+1/2,7+1/2, Лг—1/2 —2 + 1,7+1/2,& — —2г,]+1/2,к
2 кх
_ Мг+1/2,^+1/2,к+1/2 + Мг+1/2.-1/2,к+1/2 —1 г+1/2,.+1,к — —1 г+1/2,^,к
^¿+1/2,^,^+1/2 — -^--Л--
Дн-1/2,;7+1/2,А+1/2 + ^¿+1/2,^-1/2,к+1/2 4,Зг+1^,к+1/2 ~ Фз^,к+1/2
Проверим эрмитовоеть конечно-разностного аналога оператора Н^ т. е. выполнение равенства (Ф|ЯМ|Ф) = (Ф|Я*|Ф), Под (Ф|ЯМ|Ф) понимается
е ф е ф
г=о,го-1 г=1,го-1
^=1,^0-1 ^'=0,^0-1 к=1,ко-1 к=1,ко-1
,.,к+1/2 "Г" •••
г=1,го — 1
.7 = 1 ,.70-1 к=0,к0-1
Здесь верхние индексы (цифры в скобках) относятся к номеру компонент. Будем учитывать граничные условия
—1 г=-1/2,.,к = ——1 г^А^к, —1 г=г0+1/2,.,к = ——1 г=го-l/2j,k,
—1 г+1/2,.=0,к = —1 г+1/2,.=.0,к = —1 г+1/2,.,к=0 = г+1/2,.,к=к0 = 0 (аналогично для других компонент).
Для краткости докажем эрмитовоеть только для членов, связанных с N при ц = 1,
Ф
ф1 г+1/2,. {N¡+1/2,7+1/2 (—2 г+1,.+1/2 — —2 г.+ 1/2) —
1/2 2
г+1,7-
1/2 — —2 г,. -1/2 )}+
г=0,г0- 1 . =1,. 0-1
+ £ ф2 г,.+ 1/2^г+1/2.+1/2(—1 г+1/2.+1 — —1 г+1/2,.)
г=1,г0-1 .=0,70-1
—^г-1/2,.+1/2 (—1 г-1/2,.+1 — —1 г-1/2,.)}. (5)
Проведя замену индексов, первую сумму можно записать в виде
ф1 г+1/2,. N+1/2,7+1/2 (—2 г+1,.+1/2 — —2 г,.+1/2) —
г=0,г0 - 1 7=1,.0-1
ф1 г+1/2,.+1^г+1/2,.+1/2(—2 г+1,.+1/2 — —2 г,.+1/2).
г=0,г0-1 .=0,.0-2
С учетом того что ф1 г+1/2,.=0 = 0 первую го этих сумм можно продлить до ] = 0, с учетом ф1 г+1/2,.=.0 = 0 вторую сумму можно продлить до ] = — 1, В итоге получим выражение
N+1/2,7+1/2 (ф1 г+1/2,. — ф1 г+1/2,.+1)(—2 г+1,.+1/2 — —2 г,.+1/2).
г=0,г0-1 7=1,.0-1
Проделав аналогичные преобразования для второй суммы в (5), вместо этого уравнения получим
У^ Ní+1/2j+1/2{(01 i+1/2,j - ф1 i+1/2,j+l)(^2i+1,j+1/2 - ^2i,j+1/2) +
i=0,io — 1 j = 1,jo —1
+ (ф2íj+1/2 - ф2í+1,j+1/2)(^1 i+1/2,j+1 - i+1/2,j
Это выражение инвариантно относительно замены Ф на Ф, что и означает эрмитовость. Аппроксимация членов ЛФ в (1), связаныых с U и с диагональными компонентами ё, не вызывает затруднений. Но вследствие вычисления компонент Ф па сдвинутых сетках для построения эрмитовою конечно-разностного аналога членов с Д0 и недиагональными компонентами ё, зависящими от координат, понадобятся дополнительные усилия. Поступим следующим образом.
Для аппроксимации выражения вида Д0^т в некоторой точке pqr, не совпадающей с точками задания функций А0 и фт, вычислим в этой точке величину [^f7\i]pqr путем усреднения величин у/А0 ¿jfc в точку pqr. С помощью аналогичного усреднения y/A0¿¿,fc вычислим величину [л/Ао]еднфт в точках egh, в которых вычисляется га-компонента волновой функции. Затем усредним [л/Ао]еднфт в точку pgr. Обозначим результат через
[[л/А^]едкФт]рдг- В КНЧОПИО А0фт В ТОЧК Ü ЩГ ПрИМвМ ВбЛИЧИНУ [y/K^}pqr[[y/K¿}egh^m}pqr.
Продемонстрируем это на примере сопряженных членов —Д0^2, который необходимо вычислить в точке i + 1/2, j (индекс k для краткости опуекаем), и Д0^1; который необходимо вычислить в точке i, j + 1/2. В этом случае имеем
/ А , ч + л/^oí+lj /2ÍJ+I/2 + /2ÍJ-I/2 + /2Í+IJ+I/2 + /2Í+IJ-I/2
\ 0 т2 jí+l/2j' =-----^-,
, V Aoi,j + л/ Aoíj+I где /2 ¿¿+1/2 =-2-4>2i,j+l/2,
0ÍJ + \J AoiJ+l fl г+1/2 j + fl i—1/2, j + fl í+1/2 j+1 + /lí-l/2j+l
От 1 г j+1/2 - ----^-,
\/Д0 + \/Д0 i+1,j , ВДе f 1 i+1/2, j = ---V'l г+1/2 j-
Проводя исследования аналогичные сделанным для членов, связанных с N, можно
Д0
Поскольку функции ёу могут менять знак, то во избежание вычисления корня из отрицательного числа и появления комплексных величин в исходно вещественном операторе представим недиагональные компоненты ёу в виде
+—
ёу ёу — ёу,
где ё+ = ёу, если ёу > 0, и ё+ = 0 если ёу < 0 ё— = — ёу, если ёу < 0, и ё— = 0, если ёу > 0, Аппроксимация членов с ё+ и ё— производится так же, как и членов с Д0,
4. Итерационный метод
Для поиска основного энергетического уровня (наименьшего собственного значения Е) и соответствующей ему волновой функции Ф используем итерационный процесс, который по существу представляет собой решение уравнения теплопроводности (уравнения Шредингера с мнимым временем) на установление. Поскольку конечно-разностный аналог оператора Н эрмитов, то при достаточном числе итераций будет выделяться гармоника с наименьшим Е. Для увеличения параметра итераций т (шага по "времени") и устойчивости итераций применим следующую неявную схему:
—1 - —9 д2—1 д2—1 д ( д—Л
^ =ьо-ш+м+ш (м°ж I + + ^
—2 - — 9 д2—2 д2т—2 д ( д—2 ,
v1 = +ш I+ «>+
ф3- —9 д2—3 д2—3 д / д—3 ,
= + + - ] + Мфг^Ж,-Ж) + ЕЖ
ф6 - № д2— д2ф д / дфЛ ~
= + + _ ^.j + /,«,„.... С) +
L0(z) = max L(x,y,z), M0(z) = max M(x,y,z), (6)
x,y x,y
Ф"+1 = «М+г' <**> = ///***«ьфь. (7)
Здесь под дифференциальными операторами подразумеваются их конечно-разностные аналоги, n — номер итерации. Под /т(Ф) понимаетея функции — НФ за
вычетом членов с M0, L0, Уравнения (6) решались путем разложения в ряды Фурье по x y z
Величина Ер связана с энергией
Еп = (Фп(ЯФп))/(ФпФп) (8)
следующим образом, В начале итераций ЕП=0 = 0, При |En — En-1| < SE |Еп| полагается ЕП+1 = Еп, Здесь SE — малое число, В приведенных ниже расчетах оптимальным оказалось SE ~ 10-3, Итерации прекращаются при выполнении условия
тах(|(Фп)*(ЯФп — ЕпФп)|)1/2 <S max |Фп| (9)
x,y,z x,y,z
во всех точках расчетной области. Здесь S — малое число.
Член с Еп в итерационном процессе (6) необходим в силу того что используемая схема является схемой расщепления, которую можно представить в виде
Ф = (1 + тР + тЕпР )Фп, (10)
где оператор Р зависит от т и Ер и совпадает с конечно-разностным аналогом оператора Шредингера Hh только в пределе
lim Р = Hh. (11)
т
Соответственно собственные функции и собственные значения операторов P и Hh от-
т
ных значений и функций оператора Шредингера может быть велика, поскольку для более быстрой сходимости требуются достаточно большие т. Однако, если Ер равно собственному числу оператора Hh, то это число и собственная функция оператора Hh, соответствующая этому числу, являются собственным числом и собственной функцией оператора Р. Но другие собственные числа и функци и операторов P и Hh не совпадают при конечном т. Изменение Ер в процессе итераций (6) приводит к тому, что Ер стремится к собственному числу, а Фр — к постоянной функции, определяемой уравнением (H — Ер)Ф = 0, т.е. к требуемому решению.
Действительно, правая часть равенств (6) при Ф = Фр переходит в выражение (Hh + Ер)Фр. Поэтому еели Ер равно собственному числу, а Фр — собственной функции оператора то (1 + тР + тЕр)Фр = Фр независимо от т.
Интересно заметить, что, по-видимому, оператор Р является эрмитовым. При Д0 = 0 это предположение не очевидно и следует из результатов расчетов. Исследовать подроб-
PH
5. Ускорение и вычисление первого возбужденного состояния
Для данной задачи оказалась эффективной стандартная процедура ускорения. Она заключается в том, что после большого числа итераций и, когда остаются только две наиболее плохо затухающие гармоники, в силу нормировки (7) функцию Фп можно представить в виде
Фп+т - Фо + штФь (12)
где Ф0 — интересующая пас функция, Ф1 — погрешность, от которой необходимо избавиться, т — декремент затухания. Предполагая ортогональность Ф0 и Ф1; из (12) имеем
[[[ |ФП+3 - Фп+2|2<1х<1у<1г
ь? « Ш--, (13)
/ / / |ФП+2 - Фга+1| 2dxdydz
(Фр+3 — ^Фга+2)/(1 — w).
Функция
Ф1 - (Фп+2 - Фп+3)/(т2(1 - т)) (14)
в (12) представляет собой приближенное значение первого возбужденного состояния. Ф1
даться следующая картина. Первые итерации дадут с некоторой точностью энергию
и
Ф0
Ф1
Фо > Ф1.
Возможен алгоритм вычисления первого возбужденного состояния, состоящий в проведении итераций (6) до тех пор, пока не затухнут все гармоники, кроме Ф0 и Ф^ При этом время от времени проводится вычитание Фо с помощью формул, аналогичных формулам (12)—(14). После вычитания остается Ф1, Более подробно метод выглядит следующим образом.
1. Проводим итерации (6) до тех пор, пока не выделится основная гармоника Ф0, Это означает, что выполняется условие (9). Типичное 6 в приведенных ниже расчетах составляло 10-(4^7),
2. Вычисляем Ф1 по формулам (13), (14) и приближенное значение энергии первого возбужденного состояния Е1 по формуле (8), в которой в качестве Фп нспользуетея Ф1.
3. Проводим итерации (6) при Ерп = E1 = const. Такой выбор при вычислении первого возбужденного состояния необходим по причинам, указанным в разделе 4 настоящей работы. Важным фактором (см. ниже) является также то, что эти итерации проводятся уже без нормировки (7) на каждом шаге. Получающуюся последовательность функций обозначим через Фп, Итерации продолжаются до тех пор, пока Ф0 опять не станет преобладающей.
Признак того что Ф0 преобладает при ЕП = const, не равной энергии основного уровня, имеет вид не (9), а
шах(|(Фn+1 - Фп)/т - ефn| < 6 max |Фn|, (15)
где
/ (Фп+1фп) е = т —=—=--1
у (фпфп)
Заметим, что условия (9) и (15) будут выполняться одновременно только в том случае, если величина Ер с высокой точностью равна энергии, а Фп — волновой функции первого возбужденного состояния. Это имеет место при начальных после этапа 2 итерациях (затем будет расти гармоника, соответствующая основному состоянию) и служит основанием для прекращения итераций и объявления Фп волновой функцией первого
6
4. Вычисляем Ф1; но уже не по формулам (13), (14), а с помощью следующего алгоритма. Полагаем, что после большого числа итераций n функцию Фп можно представить в виде
Фп+т « <™Фо + (16)
В силу того что при вычислении первого возбужденного состояния Еп с возрастающей точностью стремится к точному значению энергии этого состояния, w1 стремится к единице. Поэтому положим в (16) w1 = 1:
Фп+т « w^0 + Ф1,
iff |Фп+3 - Ф^pdxdydz
соответственно
w
JJJ 1Фп+2 - Ф^'dxdydz
Ф1 « (Фп+3 - w0^+2)/(1 - W0).
Далее вычисляем приближенное значение энергии первого возбужденного состояния Е1 то функции Ф1 по формуле (8) и возвращаемся на этап 3,
В результате при первых после этапа 2 итерациях ЕП будет стремиться к точному значению энергии, а Фп — к волновой функции первого возбужденного состояния.
Необходимость отказа от нормировки (7) связана со следующим. Формула (12) в обозначениях (16) имеет вид
~ (а^о'Па^)1^ ~
где
|Ф0|2dxdydz, а^ =11 / |Фi|2dxdydz.
При поиске волновой функции первого возбужденного состояния Ф1 в отличие от алгоритма поиска основного состояния отношение Ф1/Ф0 на этапе 4 остается хотя и малым, но конечным. Поэтому наличие второго члена в (17) будет давать конечную ошибку и описываемый алгоритм при использовании формул (12)—(14) сходиться не будет. Данный алгоритм позволяет вычислить первое возбужденное состояние с точностью 8 = 10-(2^3) за число итераций, сравнимое с числом итераций, необходимым для вычисления Ф0 с точностью 8 = 10-10.
Возможен другой способ поиска первого возбужденного состояния. Если бы точное значение энергия Е1 было известно, то волновую функцию основного состояния можно вычислить с высокой точностью, а затем провести итерации при Ер = E1 = const, вычитая на каждом шаге компоненту волновой функции, соответствующую основному состоянию, В результате будет выделяться волновая функция первого возбужденного состояния. Но так как значение Е1 не известно, то необходимо производить итерации по Е1; каждый раз вычисляя с высокой точностью волновую функцию основного состояния, которая в случае неявной схемы зависит от Ер и т. Кроме того, в интересующей нас конфигурации волновая функция основного состояния многократно вырождена (см. Приложение), Соответственно необходимо вычислять разнообразные волновые функции основного состояния, отвечающие различным симметриям. Добавим, что схема расщепления (6) вносит некоторое неравноправие в вычисление компонент вектора Ф при значении Ер, не равном энергии основного состояния. Это может привести к тому, что связь между волновыми функциями, соответствующими симметрии относительно замены жу ^ уж, будет иметь намного более сложный вид, чем в случае оператора
6. Результаты расчетов при постоянных Ь, М, N и А0
Рассмотрим случай постоянных Ь = 31.5, М = 5.75, N = 34.14 в безразмерных единицах, До = 0.289 эВ, Эти значения соответствуют параметрам германия. Тензор напряжений е = 0, Потенциал и = — 0.54 эВ в некоторой области П и и = 0 вне ее. Область П состоит из смачивающей поверхноети г0/2 < г < г0/2 + = 0.54 нм, и квантовой
точки — области, имеющей форму правильной пирамиды высотой Л,р = 1 нм. Основание пирамиды находится в плоскости г = г0/2 + и представляет собой квадрат со стороной 10 нм. Грани основания ориентированы вдоль осей координат. Вершина пирамиды расположена па прямой х = я0/2, у = у0/2.
Результаты расчетов демонстрируют, что в этом случае волновая функция основного состояния локализована в центральной части пирамиды. На рис, 2 приведено распределение плотности вероятности |Ф|2 = ^ ф^фа, Штрихом показаны контуры области П,
а=1,6
Из расчетов следует, что |Ф|2 обладает симметрией
|Ф(х,у)| = |Ф(-х,у)| = |Ф(х, у) | = |Ф(-х, у) | =
= |Ф(у,х)| = |Ф( у,х)| = |Ф(у, -х)| = |Ф(-у, х)|.
При этом итерации начинались с заведомо несимметричной функции.
Поскольку реальные размеры матрицы кремния велики, а расчетная область конечна, то первая серия расчетов имела цель найти размер области, при котором влияние границ будет несущественным. Оказалось, что при х0 х у0 х г0 = 40 х 40 х 20 нм дальнейшее увеличение расчетной области приводит к изменению значений энергии Е и функции |Ф|2 в пятом знаке. Ниже данные будут приводиться для области 40 х 40 х 20 нм. Рассмотрим, как результаты расчетов зависят от числа узлов копечпо-разпоетпой сетки г0, ]0, к0 по х, У г соответственно. Расчеты на сетках 32 х 32 х 80, 64 х 64 х 160, 128 х 128 х 320 дают значения энергии основного уровня Е = -0.35414, -0.33287, -0.32160, что отвечает первому порядку сходимости по пространственным шагам. Полагая, что вычисленная при шаге Н энергия равна сумме точного значения энергии и величины, пропорциональной Н, из этих данных получаем точное значение Е « -0.309 эВ, Распределение |Ф|2 вдоль вертикальной прямой, проходящей через вершину пирамиды, и прямой х = х0/2, г = г0/2 для этих сеток представлено на рис, 3, Штрихпунк-тирная, пунктирная и сплошная линии — сетки 32 х 32 х 80, 64 х 64 х 160 и 128 х 128 х 320, Видно, что имеет место сходимость и погрешность расчетов на сетке 64 х 64 х 160 составляет по более 7%.
Следует отметить, что погрешность по шагам пространственной сетки имеет первый порядок. Отличие от второго порядка сходимости связано с тем, что функция и разрывна, а узлы пространственной сетки пе совпадают с поверхностью разрыва.
Рис. 2. Распределение |Ф|2 в плоскостях у = У0/2 (проходит через высоту пир амиды) иг = 20/2 + (основание пирамиды); штрих — контур области, занимаемой германием
Рис. 3. Распределение |Ф|2 вдоль высоты пирамиды (а) и вдоль прямой х = Х0/2, г = 20/2 (б)', штрихпупктир, пунктир и сплошная линии соответствуют сеткам 32 х 32 х 80, 64 х 64 х 160 и 128 х 128 х 320
Рис. 4. Распределение |Ф|2 первого возбужденного состояния в плоскостях у = У0/2 (проходит через высоту пирамиды) иг = 20/2 + Нт (основание пирамиды)
В приведенных расчетах итерации прекращались при шах(|(Фп)*(НФп—ЕпФп)|)1/2 <
6 = 10-10, ^то соответствовало ошибкам округления. При этом оптимальное т ~ 0.06. Требовалось несколько сотен итераций, что для сетки 64 х 64 х 160 составляет 1.5 ч на РС аеег5710 (1 Ггц), Заметим, что дня практических пужд итерации можно прекращать уже при 6 ~ 10-5, ^^этение Е вычисляется с высокой точностью при еще больших 6.
На рис. 4 показано распределение |Ф|2 первого возбужденного состояния в рассматриваемом случае. В точке, расположенной примерно в центре расчетной области, |Ф|2 = 0, При удалении от этой точки |Ф|2 увеличивается, достигая максимума, а затем уменьшается. При этом максимумы вдоль оси г существенно меньше максимумов вдоль х у 0.25
7. Результаты расчетов при Ь, М, N А0 = const и ц = 1
Рассмотрим случай, когда в области О (см, предыдущий раздел) находится германий с параметрами £Се = 31.5, МСе = 5.75 (безразмерные значения), Дове = 0.289, ио,^ = — 0.54 эВ, окруженный кремнием с параметрами = 6.22, Мд; = 3.22 N1 = 8.28 (безразмерные значения), Д0^ = 0.044 ^ = 0 эВ, Значение параметра = 34.14, Расчет с такими параметрами сталкивается с трудностями, евязапыми с эффектами, описанными в разделе 2, Поэтому исследуем варианты, в которых параметр N в области О равен некоторому значению Расчетная области, как и в предыдущем случае, имела размер 40 х 40 х 20 нм.
При N < 17 распределение |Ф|2 и значение Е принципиально не отличаются от описанных в разделе 6, Например, на сетке 64 х 64 х 80 Е с увеличением N от 0 до 17 плавно уменьшается от —0.351 до — 0.374 эВ, Однако при N « 18 решение качественно меняется. Максимумы распределения |Ф|2 достигаются на границе германий—кремний, на распределении |Ф|2 появляются мелкомасштабные осцилляции (рис, 5), Энергия Е достигает значения — 0.47 эВ, При дальнейшем увеличении N Е резко уменьшается (например, при N = 20 Е = — 0.765 эВ), а осцилляции |Ф|2 на границе германий—кремний увеличиваются. При N > 18 сходимость существенно ухудшается. Эти результаты находятся в полном соответствии с аргументами, высказанными в разделе 2,
При разрывных коэффициентах делать какие-либо выводы о сходимости численного решения к точному затруднительно. Поэтому были проведены расчеты с плавным переходом значений параметров, соответствующих кремнию, к значениям параметров, соответствующих германию. Кроме того, исследована аналогичная двумерная задача, позволяющая использовать более мелкую пространственную сетку. Эти расчеты также указывают па то, что появление пефизичеекого решения является свойством исходной дифференциальной задачи, а пе дефектом разностного алгоритма. Сравнение результатов расчетов с коэффициентами с различной степенью сглаживания и без сглаживания показывают, что при N < 17 предложенная в работе схема имеет в случае разрывных коэффициентов первый порядок аппроксимации по шагу пространственной сетки К,
Из расчетов с квантовой точкой различных размеров следует, что при равном реальному значению дня германия, пефизичеекое решение возникает только при достаточно мелком шаге К. Если шаг К большой, то в качестве гармоники с наименьшей
:=г„/2+й
20 40
30 X
40 40
Рис. 5. Распределение |Ф|2 в плоскостях у = уо/2 (проходит через высоту пир амиды) иг = го/2 + (основание пирамиды) при разрывном N (Ж* = 18) и ц = 1
энергией выделяется гармоника с характерным размером порядка размеров пирамиды. Это также полностью соответствует результатам, изложенным в разделе 2,
8. Результаты расчетов при Ь, М, N А0 = const и ^ = 0
Расчеты при ц = 0 и параметрах, соответствующих германию и кремнию, показывают наличие физического решения с характерным размером порядка размеров включения германия, подобного решению, представленному на рис, 2, Сходимость имеет первый
Н
значение энергии основного состояния, вычисленное по расчетам на последовательности -0.356
9. Расчеты при ё = 0
Методические расчеты, проведенные с ё = 0 показали, что наличие членов с ё не приводит к каким-либо вычислительным трудностям или принципиальным эффектам, подобным эффекту, описанному в разделе 2, Исследование физических особенностей решения, связанных с ё = 0, не входит в цели настоящей работы.
Заключение
В работе построен конечно-разностный аналог оператора Шредингера (1) для шеети-зонной к • р-модели, Этот конечно-разностный оператор имеет дивергентную форму, компактен (для вычисления Н/Ф в узле конечно-разностной сетки используется Ф только в ближайших узлах) и эрмитов.
Если входящие в (1) коэффициенты зависят от пространственных координат гладким образом, то предложенный оператор аппроксимирует (1) со вторым порядком точности, Если коэффициенты разрывны, то порядок аппроксимации падает до первого.
Предложен неявный итерационный метод, позволяющий эффективно находить волновую функцию и энергию основного и первого возбужденного состояний оператора (1), Показано, что при описании полупроводниковых структур с изменяющимся в пространстве параметром N вид шеетизонной к •р-модели неоднозначен, В случае полупроводниковых структур, в которых имеет место контакт таких материалов, в которых N одного из них оказывается существенно больше Ь + М другого (к таким структурам относится, в частности, чрезвычайно распространенная структура кремний—германий), выбор значения ц ограничен. При неправильном выборе значения ^ оператор Шредингера имеет нефизический спектр. При численном решении задач в системах с большими характерными размерами (сотни нанометров) нефизическая часть спектра может не проявляться при достаточно больших значениях шагов конечно-разностной сетки. Однако при расчете наноструктур правильный выбор параметра ^ имеет принципиальное значение. Заметим также, что уравнения воеьмизонной к • р-модели имеют структуру, аналогичную структуре уравнений шестизонной модели. Поэтому в воеьмизонной модели также могут проявляться подобные эффекты.
Выбор значения ^ для конкретной гетероетруктуры должен основываться либо на более глубоком физическом анализе "из первых принципов", который, по-видимому, невозможен из-за своей сложности, либо на экспериментальных данных.
Авторы выражают благодарность С, Б, Медведеву и Ю.Н, Моро ко г,у за полезные обсуждения материала представленной работы.
Приложение. Симметрия задачи
Пусть коэффициенты До, Ь, М, N и симметричны относительно отражения координат ж, у (в том числе одной из координат) и замены ж ^ у у ^ ж, а тензор ё имеет симметрию
ё12(ж, у) = ё12( Ж, у) = ё12(ж, -у), ё1з(ж, у) = ё13( ж, у) = ё1з(ж, -у), ё23(ж, у) = ё1з(у,ж),
ё2з(ж,у) = ё2з( ж, у) = -ё2з(ж, -у).
Под точками ж = у = 0 подразумеваются точки вертикальной оси, проходящей через центр основания пирамиды. Для краткости зависимость от г здесь опущена.
Пусть вектор-функция Ф с компонентами фт является решением (1) с энергией состояния (не обязательно основного) Е. Тогда вектор-функция Ф(1) с компонентами Ф^, соответству-
ж -ж
Ф(11)(ж,у) = -ф4(-ж,у), Ф21)(ж,у)= фб(-ж,у), Фз1)(ж,у)= фб(-ж,у),
ф41)(ж,у) = -ф1(-ж,у), ф51)(ж,у)= ф2(-ж,у), ф61)(ж,у)= фз(-ж,у),
Е
симметриям выражения в более компактном виде:
Ф(1)(ж,у) = (-^4, фб, фб, -^1, ф2, фз)(-ж,у), Ф(1) = Рх^-хФ, Ф(2)(ж,у) = (ф4, -фб, фб, -ф1, ф2, -фз)(ж, -у), Ф(2) = Ру^-уФ, Ф(з) (ж, у) = (ф1, ф2, -фз,, -ф4, -фб, фб)( ж, -у), Ф(з) = -Ру^-уРх^-хФ,
Ф(4)(ж, у) = (фб, ф4, фб, -гф2, -гфь -гфз)(у, ж), Ф(4) = Рху^ухФ Ф(б)(ж, у) = (ф2, -ф1, -фз, ¿фб, -#4, -¿фб)(у, -ж), Ф(б) = -¿Рх^-хРху^ух Ф(б)(ж, у) = (ф2, -ф1, фз, -гфб, гф4, -гфб)(-у,ж), Ф(б) = ¿Ру^-уРху^ухФ
Ф(7)(ж,у) = (фб, ф4, -фб, гф2, гф1, -гфз)(-у, -ж), Ф(7) = -Ру^-уРх^-хРху^ухФ.
Еще одной симметрией является симметрия относительно направления спина. Она состоит
ФЕ
Ф(ж,у,г) = (-ф6, ф5, -ф6, ф6, ф2, фз)(ж,у,г) также является решением с той же энергией состояния.
Список литературы
fl] Voigtlander В. Fundamental processes in Si/Si and Ge/Si epitaxy studied by scanning tunneling microscopy during growth // Surface Sci. Rep. 2001. Vol. 43. P. 127-254.
[2] Берман Г.П., Дулен Г.Д., Майньери P., Цифринович В.И. Введение в квантовые компьютеры. Москва; Ижевск: Ин-т компьютерных исследований, 2004. 188 с.
[3] Ненашев А.В. Моделирование электронной структуры квантовых точек Ge в Si. Дисс. ... канд. физ.-мат. наук. Новосибирск, ИФП СО РАН, 2004. 242 с.
[4] Вир Г.Л., Пикус Е.Г. Симметрия и деформационные эффекты в полупроводниках. М.: Наука, 1972. 584 с.
[5] Кунин С. Вычислительная физика. М.: Наука, 1982. 518 с.
[6] Цуканов А.В. Зарядовый кубит па основе одноэлектронной донорской пары с оптическим и электростатическим управлением // Тр. физико-технолог, ин-та. 2008. Т. 19. С. 7-22.
f7] Stier О., Bimberg D. Modeling of strained quantum wires using eight-band k • p theory // Phvs. Rev. B. 1997. Vol. 55, No. 12. P. 7726-7732.
[81 Stier O., Grundmann M., Bimberg D. Electronic and optical properties of strained quantum dots modeled bv 8-band k • p theory // Ibid. 1999. Vol. 59, No. 8. P. 5688-5701.
Поступила в редакцию 28 мая 2010 г., с доработки — 11 августа 2010 г.
