Теплофизика высоких температур, 2021, T. 59, № 6, стр. 826-835

ВВЕДЕНИЕ

Однокомпонентная плазма (one-component plasma – OCP) состоит из одного вида заряженных частиц, погруженных в жесткий нейтрализующий фон. OCP, состоящую из вырожденных электронов (однородный электронный газ, uniform electron gas − UEG), можно использовать как фундаментальную модель простых металлов. Хотя UEG не соответствует реальной физической системе, тем не менее точное описание его свойств имеет решающее значение для теории функционала плотности (density functional theory – DFT) [1, 2]. В рамках DFT сильно взаимодействующая многоэлектронная система преобразуется в эффективную невзаимодействующую систему частиц посредством введения эффективного внешнeго потенциала, который содержит все обменные и корреляционные эффекты и который можно аппроксимировать обменно-корреляционной энергией системы [3, 4].

Надежные данные для UEG при нулевой температуре были получены с помощью квантового метода Монте-Карло [5, 6] и широко использовались в расчетах методом DFT. Улучшенные данные для основного состояния представлены в [7]. Однако все эти результаты ограничены нулевой температурой и неприменимы к сильно неидеальным плазменным средам.

Надежные расчеты свойств плазменных сред при конечных температурах можно выполнить в рамках фейнмановской формулировки квантовой механики с помощью метода Монте-Карло, который позволяет рассчитать термодинамические величины, представленные в виде интегралов по траекториям (path integral Monte Carlo – PIMC [8–12]. За последние два десятилетия для фермионных систем были получены многочисленные результаты, включая, например, свойства коррелированных электронов в квантовых точках [13, 14] или в плотной плазме [15–17]. Эти результаты использовались в DFT в качестве реперных для построения приближенных зависящих от температуры выражений для локальной (спиновой) плотности (the local (spin-)density approximation – L(S)DA) и для более сложных градиентных приближений [18–24]. Энергия и другие результаты расчетов методом Монте-Карло для UEG важны для понимания формализма отклика плотности [25–32].

Основная трудность квантового моделирования системы фермионов с помощью PIMC – “проблема знаков”, связанная с требованием антисимметризации матрицы плотности. В результате все термодинамические величины представляются как суммы положительных и отрицательных слагаемых, относящиеся к четным и нечетным перестановкам. Результат определяется малой разностью этих двух больших вкладов. Численные расчеты в этом случае сильно затруднены. Для преодоления этой трудности были разработаны различные подходы, но “проблема знаков” для сильно коррелированных фермионов так и не была полностью решена за последние 50 лет.

В [33, 34] была использована вигнеровская формулировка квантовой механики, чтобы избежать антисимметризации матричных элементов и тем самым избежать “проблемы знаков” и реализовать принцип запрета Паули. Однако этот подход неприменим при высоком вырождении фермионов.

В [35], чтобы избежать “проблемы знаков”, был развит метод Монте-Карло с ограничениями (restricted fixed-node path-integral Monte Carlo – RPIMC) для моделирования UEG при конечной температуре. В RPIMC учитываются только положительные перестановки, поэтому точность результатов трудно количественно оценить [36].

Более интересными подходами являются методы permutation blocking path integral Monte Carlo (PB-PIMC) и configuration path integral Monte Carlo (CPIMC) [22]. Основная идея CPIMC состоит в представлении матрицы плотности в виде интеграла по траекториям в пространстве чисел заполнения. Однако оказывается, что метод CPIMC также подвержен “проблемe знаков”. В PB-PIMC сумма по перестановкам представлена в виде обменного детерминанта, который может быть вычислен прямыми методами линейной алгебры, что позволит частично преодолеть “проблему знаков”. Однако недостатком этого метода является то, что обменные детерминанты представляют собой знакопеременныe функции, что ухудшает точность расчетов методом PIMC и порождает новую “проблему знаков детерминанта”.

В данной работе, чтобы избежать “проблемы знака детерминанта”, для сильно вырожденной фермионной системы разработано приближение, в котором обменный детерминант преобразован в детерминант Грама, не принимающий отрицательных значений. Чтобы повысить точность этого приближения и учесть интерференционные эффекты кулоновского и обменного взаимодействия электронов, разработано новое представление мaтрицы плотности в виде интегралов по траекториям, в котором взаимодействие включено в обменный детерминант. Для выполнения расчетов разработан новый метод (modified fermionic path integral Monte Carlo approach – MFPIMC), который позволяет в значительной степени избежать “проблемы знаков”. В работе исследуется квантовое упорядочение электронов в кулоновских системах в каноническом ансамбле при конечных температурах. Наблюдаемое упорядочение электронов вызвано эффектом исключенного объема обменных “дырок”, которые возникают за счет ферми-отталкивания и кулоновского взаимодействия электронов с этими положительно заряженными “дырками”. В результате формируются нейтральные обменно-корреляционные экситоны в “море ферми-электронов”, возможность образования которых рассматривалась в работах [37–39].

Концепция обменно-корреляционных “дырок” позволяет как построить интуитивно понятную физическую картину явления, так и получить количественные результаты для обменно-корреляционной энергии. Более того, обменно-корреляционная “дырка” является одним из фундаментальных понятий теории функционала плотности [37–39].

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

Понятие “обменная дырка” может быть уточнено путем добавления других корреляций, таких как экранирование. Экранировка обычно слабее обменного взаимодействия, тем не менее можно говорить об обменно-корреляционных “дырках” с эффективным зарядом. Обменная “дырка” может быть положительно заряжена, что нейтрализует заряд электрона на нескольких атомных расстояниях. Дополнительный электрон может присоединяться к обменной дырке, образуя нейтральный обменный экситон. Экранированное кулоновское взаимодействие между электронами и нейтральные обменные экситоны позволяют использовать модели слабо взаимодействующей ферми-жидкости в металлах, несмотря на высокую плотность электронов. Учет обменно-корреляционных экситонов важен для построения аналитических и численных моделей теории функционала плотности [37–39]. Данный эффект не возникает для электронов с противоположным спином и симметричной волновой функцией.

В выполненных MFPIMC-расчетах появление обменно-корреляционных экситонов подтверждается парными функциями распределения для двухкомпонентной плазмы и однородного электронного газа. В дальнейшем для анализа трехчастичных корреляций, связанных с экранировкой и обменными экситонами, рассматриваются парные функции распределения для электронов с параллельными и антипараллельными спинами, а также парные функции распределения для положительных зарядов и электронов. Анализ парных корреляций позволяет получить оценки плотностей заряда, средних расстояний и плотностей обменной “дырки” и электронов. Такой эффект никогда ранее не наблюдался при моделировании с использованием стандартного PIMC.

МНОГОЧАСТИЧНАЯ МАТРИЦА ПЛОТНОСТИ

Рассмотрим нейтральную систему квантовых электронов и положительных классических (для простоты) зарядов с гамильтонианом $\hat {H} = \hat {K} + {{\hat {U}}^{c}}$, содержащим кинетическую $\hat {K}$ и кулоновскую ${{\hat {U}}^{c}} = \hat {U}_{{pp}}^{c} + \hat {U}_{{ee}}^{c} + U_{{ep}}^{c}$ энергии. Термодинамические свойства системы частиц в каноническом ансамбле с заданной температурой $T$ и фиксированным объемом $V$ полностью описываются оператором плотности $\hat {\rho } = {{{\text{e}}}^{{ - \beta \hat {H}}}}$ со статистической суммой $Z({{N}_{e}},{{N}_{p}},V;\beta )$ = $\frac{1}{{{{N}_{e}}!{{N}_{p}}!}}\sum\nolimits_\sigma \int_V {{\text{d}}x\rho (x,\sigma ;\beta )} ,$ где β = = 1/k_BT и $\rho (x,\sigma ;\beta ) = \left\langle {x\left| {{{{\text{e}}}^{{ - \epsilon \hat {H}}}}} \right|x} \right\rangle $ обозначает диагональные элементы матрицы плотности в координатном представлении.

Здесь $x = \{ {{x}_{e}},{{x}_{p}}\} $ (отнесенные к тепловой длине волны ${{\tilde {\lambda }}_{a}} = \sqrt {\frac{{2\pi {{\hbar }^{2}}\beta }}{{{{m}_{a}}M}}} $) – пространственные координаты электронов и положительных зарядов, т.е. ${{x}_{a}} = \{ {{x}_{{1,a}}}, \ldots ,{{x}_{{{{N}_{a}},a}}}\} $ и $\sigma = \{ {{\sigma }_{{1,e}}}, \ldots ,{{\sigma }_{{{{N}_{e}},e}}}\} $ – спиновые степени свободы электронов.

Как правило, точные матричные элементы матрицы плотности неидеальных квантовых систем частиц неизвестны, но могут быть получены с помощью интегралов по траекториям, которые основаны на операторном тождестве ${{{\text{e}}}^{{ - \beta \hat {H}}}} = {{{\text{e}}}^{{ - \epsilon \hat {H}}}} \cdot {{{\text{e}}}^{{ - \epsilon \hat {H}}}} \ldots {{{\text{e}}}^{{ - \epsilon \hat {H}}}}$ $(\epsilon = {\beta \mathord{\left/ {\vphantom {\beta M}} \right. \kern-0em} M})$ [8]. Тождество включает $M$ одинаковых высокотемпературных матриц плотности с температурой $MT$:

где индекс $m = 0, \ldots ,M$ маркирует недиагональные матричные элементы высокотемпературных матриц плотности. Здесь каждый высокотемпературный сомножитель представлен в виде $\left\langle {{{x}^{{(m)}}}\left| {{{{\text{e}}}^{{ - \epsilon \hat {H}}}}} \right|{{x}^{{(m + 1)}}}} \right\rangle \approx $ $ \approx \left\langle {{{x}^{{(m)}}}\left| {{{{\text{e}}}^{{ - \epsilon {{{\hat {U}}}_{m}}}}}} \right|{{x}^{{(m + 1)}}}} \right\rangle \rho _{0}^{{(m)}}$, ($\rho _{0}^{{(m)}} = \left\langle {{{x}^{{(m)}}}\left| {{{{\text{e}}}^{{ - \epsilon \hat {K}}}}} \right|{{x}^{{(m + 1)}}}} \right\rangle = $ $ = {{{\text{e}}}^{{ - \pi {{{\left| {{{x}^{{(m)}}} - {{x}^{{(m + 1)}}}} \right|}}^{2}}}}}$) – матрица плотности идеaльной системы частиц) с ошибкой порядка ${1 \mathord{\left/ {\vphantom {1 {{{M}^{2}}}}} \right. \kern-0em} {{{M}^{2}}}}$, возникающей из-за пренебрежения коммутатором ${{{{\epsilon }^{2}}} \mathord{\left/ {\vphantom {{{{\epsilon }^{2}}} 2}} \right. \kern-0em} 2}\left[ {K,{{U}^{c}}} \right]$. В пределе $M \to \infty $ погрешность всего произведения стремится к нулю $({{ \propto 1} \mathord{\left/ {\vphantom {{ \propto 1} M}} \right. \kern-0em} M})$, и возникает точное представление статистической суммы через интеграл по траекториям. Теперь каждая частица представлена траекторией, состоящей из набора $M$ координат ${{x}^{{(m)}}}$. Спин порождает спиновую часть матрицы плотности ($\mathcal{S}$) с обменными эффектами, которые учитываются благодаря операторам перестановки $\hat {P}$, действующим на координаты электрона в ${{x}^{{(M)}}}$ и проекции спина $\sigma $. Суммирование ведется по всем перестановкам с четностью ${{\kappa }_{P}}$. Предполагается, что в термодинамическом пределе основной вклад в сумму по спиновым переменным дает член с равным количеством (${{{{N}_{e}}} \mathord{\left/ {\vphantom {{{{N}_{e}}} 2}} \right. \kern-0em} 2}$) электронов с одинаковой проекцией спина [40, 41], поэтому ${\text{det}}\left\| {\Psi (x)} \right\| = $ ${\text{det}}\left\| {{{{\text{e}}}^{{ - \pi {{{\left| {x_{{k,e}}^{{(M)}} - x_{{t,e}}^{{(0)}}} \right|}}^{2}}}}}} \right\|_{1}^{{{{{{N}_{e}}} \mathord{\left/ {\vphantom {{{{N}_{e}}} 2}} \right. \kern-0em} 2}}} \times $ $ \times \,\,{\text{det}}\left\| {{{{\text{e}}}^{{ - \pi {{{\left| {x_{{k,e}}^{{(M)}} - x_{{t,e}}^{{(0)}}} \right|}}^{2}}}}}} \right\|_{{{{{{N}_{e}}} \mathord{\left/ {\vphantom {{{{N}_{e}}} 2}} \right. \kern-0em} 2}}}^{{{{N}_{e}}}}$.

Элементы матрицы плотности включают эффективное парное взаимодействие $U = \frac{1}{2}\sum\nolimits_{k,t} {{{\Phi }_{{kt}}}} $, которое аппроксимируется потенциалом Кельбга, задаваемым выражением: ${{\Phi }_{{ab}}}({{x}_{{ab}}};\epsilon ) = \frac{{{{e}_{a}}{{e}_{b}}}}{{{{{\tilde {\lambda }}}_{{ab}}}{{x}_{{ab}}}}} \times $ $ \times \,\,\left[ {1 - {{{\text{e}}}^{{ - x_{{ab}}^{2}}}} + \sqrt \pi {{x}_{{ab}}}\left\{ {1 - {\text{erf}}({{x}_{{ab}}})} \right\}} \right].$

Здесь ${{\tilde {\lambda }}_{{ab}}}{{x}_{{ab}}} = \left| {{{q}_{{k,a}}} - {{q}_{{t,b}}}} \right|{{\tilde {\lambda }}_{e}}$, $\tilde {\lambda }_{{ab}}^{2} = {{2\pi {{\hbar }^{2}}\epsilon } \mathord{\left/ {\vphantom {{2\pi {{\hbar }^{2}}\epsilon } {{{m}_{{ab}}}}}} \right. \kern-0em} {{{m}_{{ab}}}}}$, $\tilde {\lambda }_{{ab}}^{2} = {{2\pi {{\hbar }^{2}}\epsilon } \mathord{\left/ {\vphantom {{2\pi {{\hbar }^{2}}\epsilon } {{{m}_{{ab}}}}}} \right. \kern-0em} {{{m}_{{ab}}}}}$, ${1 \mathord{\left/ {\vphantom {1 {{{m}_{{ab}}}}}} \right. \kern-0em} {{{m}_{{ab}}}}} = {1 \mathord{\left/ {\vphantom {1 {{{m}_{a}}}}} \right. \kern-0em} {{{m}_{a}}}} + {1 \mathord{\left/ {\vphantom {1 {{{m}_{b}}}}} \right. \kern-0em} {{{m}_{b}}}}$, ${1 \mathord{\left/ {\vphantom {1 {{{m}_{{ab}}}}}} \right. \kern-0em} {{{m}_{{ab}}}}} = $ $ = {1 \mathord{\left/ {\vphantom {1 {{{m}_{a}}}}} \right. \kern-0em} {{{m}_{a}}}} + {1 \mathord{\left/ {\vphantom {1 {{{m}_{b}}}}} \right. \kern-0em} {{{m}_{b}}}}$, $\tilde {\lambda }_{e}^{2} = {{2\pi {{\hbar }^{2}}\epsilon } \mathord{\left/ {\vphantom {{2\pi {{\hbar }^{2}}\epsilon } {{{m}_{e}}}}} \right. \kern-0em} {{{m}_{e}}}}$, ${\text{erf(}}x)$ – стандартная функция ошибок. При ${{x}_{{ab}}} \gtrsim {{\lambda }_{e}}$ потенциал Кельбга совпадает с кулоновским, а при ${{x}_{{ab}}} \to 0$ конечен в силу своей квантовой природы [40].

НОВОЕ ПРЕДСТАВЛЕНИЕ МАТРИЦЫ ПЛОТНОСТИ

Недостатком уравнения (1) является знакопеременный детерминант ${\text{det}}\left\| {\Psi (x)} \right\|$, который порождает “проблему знаков”, экспоненциально ухудшающую точность расчетов с помощью PIMC. Далее для простоты рассмотрим бесспиновые фермионы. Обобщение для фермионов со спиновой степенью свободы тривиально. Чтобы избежать “проблемы знаков”, заменим переменные интегрирования ${{x}^{{(m)}}}$ на ${{q}^{{(m)}}}$: ${{x}^{{(m)}}} = (Px - x)\frac{m}{M} + x + {{q}^{{(m)}}}$ для любой конкретной перестановки $P$ в стандартном представлении матрицы плотности через интеграл по траекторям (1) [8, 11, 42]:

где ${{\eta }^{{(m)}}} \equiv {{q}^{{(m)}}} - {{q}^{{(m + 1)}}}$.

Из выражения для потенциальной энергии $\epsilon U\left( {(Px - x)\frac{m}{M} + x + {{q}^{{(m)}}}} \right)$ в уравнении (2) следует, что к каждой координате вектора ${{x}_{k}}$ может быть добавлено слагаемое $({{x}_{{Pk}}} - {{x}_{k}})\frac{{2m}}{M}$, если $Pk \ne k$. Парные потенциалы в $U$ для $Pk = k$ и $Pt = t$ обозначим далее символами ${{\Phi }_{{kt}}} = \Phi (\left| {{{r}_{{kt}}} + q_{{kt}}^{{(m)}}} \right|)$, в противном случае используется обозначение ${{\tilde {\Phi }}_{{kt}}} = \Phi \left( {\left| {{{r}_{{tk}}}\frac{{2m}}{M} + {{r}_{{kt}}} + q_{{kt}}^{{(m)}}} \right|} \right)$. Разность матричных элементов $\Delta {{\tilde {\Phi }}_{{kt}}} = {{\tilde {\Phi }}_{{kt}}} - {{\Phi }_{{kt}}}$ не равна нулю, только если перестановка изменяет число $k$ или/и $t$ ($Pk \ne k$ или/и $Pt \ne t$). Таким образом, ненулевые элементы $\Delta {{\tilde {\Phi }}_{{kt}}}$ любой данной перестановки $P$ образуют строки и столбцы с соответствующими номерами $k$ и $t$.

В сумме по перестановкам в уравнении (2) каждая перестановка $P$ состоит из композиции циклических перестановок без общих компонентов [43], поэтому без ограничения общности можно рассматривать только циклические перестановки. В свою очередь, циклическую перестановку можно представить как композицию парных циклических перестановок $(k,t)$ [43].

Рассмотрим приближение, позволяющее ослабить “проблему знаков” и использовать преимущества прямых методов линейной алгебры для вычисления детерминантов, описывающих обменное взаимодействие. Поясним с помощью таблицы парных потенциалов взаимодействия используемое приближение на примере простой перестановки, состоящей из двух нетривиальных циклических перестановок и четырех тождественных циклов ${{P}_{c}} = (1,3,5),(4,7),$ $(2,2),(6,6),(8,8),(9,9)$ для девяти электронов. В табл. 1 заменим парные потенциалы ${{\tilde {\Phi }}_{{kt}}}$ на ${{\Phi }_{{kt}}}$ в строке $k$ и столбце $t$, кроме общего элемента на их пересечении, который пометим звездочкой $ \star {{\tilde {\Phi }}_{{kt}}}$. Тогда в табл. 2 в строке $k$ и столбце $t$ все матричные элементы $\Delta {{\tilde {\Phi }}_{{kt}}} = {{\tilde {\Phi }}_{{kt}}} - {{\Phi }_{{kt}}}$ равны нулю, кроме $ \star \Delta {{\tilde {\Phi }}_{{kt}}}$.

В результате этого приближения в уравнении (2) для ${{P}_{c}}$ имеем

где $\phi _{{kt}}^{*} = \exp \{ {{ - \pi {{{\left| {{{r}_{{kt}}}} \right|}}^{2}}} \mathord{\left/ {\vphantom {{ - \pi {{{\left| {{{r}_{{kt}}}} \right|}}^{2}}} M}} \right. \kern-0em} M}\} $$\exp \left\{ { - \frac{1}{2}\sum\nolimits_{m = 0}^{M - 1} {(\epsilon \star {\kern 1pt} \Delta {{{\tilde {\Phi }}}_{{kt}}})} } \right\}$ и $\phi _{{kk}}^{*} = 1$. (Произведение $\phi _{{31}}^{*}\phi _{{15}}^{*}\phi _{{53}}^{*}$ будет учтено симметричным циклом, как в матрице $3 \times 3$.)

Благодаря этому приближению частицы любого заданного цикла взаимодействуют с остальными частицами с помощью парного потенциала ${{\Phi }_{{kt}}}$ вместо ${{\tilde {\Phi }}_{{kt}}}$, а частицы, входящие в циклическую перестановку, последовательно взаимодействуют с помощью потенциала ${{\tilde {\Phi }}_{{kt}}}$. В табл. 3 приведены парные межчастичные потенциалы для полученного приближения.

Сумма по всем перестановкам сворачивается в детерминант, и в результате проведенных преобразований диагональные матричные элементы матрицы плотности можно переписать в виде интеграла по “замкнутым” траекториям $\{ {{q}^{{(0)}}}, \ldots ,{{q}^{{(M)}}}\} $:

где ${\text{det}}\left\| {\psi (x)} \right\| = {\text{det}}\left\| {{{\phi }_{{kt}}}} \right\|_{1}^{{{{N}_{e}}/2}}\left\| {{\text{det}}\left| {{{\phi }_{{kt}}}} \right.} \right\|_{{{{N}_{e}}/2}}^{{{{N}_{e}}}}$, q⁽⁰⁾ = q^(M)== 0,

(4)

$\begin{gathered} {{\phi }_{{kt}}} = \exp \left\{ {{{ - \pi {{{\left| {{{r}_{{kt}}}} \right|}}^{2}}} \mathord{\left/ {\vphantom {{ - \pi {{{\left| {{{r}_{{kt}}}} \right|}}^{2}}} M}} \right. \kern-0em} M}} \right\} \times \\ \times \,\,\exp \left\{ { - \frac{1}{2}} \right.\sum\limits_{m = 0}^{M - 1} {\left( {\epsilon \Phi \left( {\left| {{{r}_{{tk}}}\frac{{2m}}{M} + {{r}_{{kt}}} + q_{{kt}}^{{(m)}}} \right|} \right)} \right.} - \\ \left. {\left. { - \,\,\epsilon \Phi \left( {\left| {{{r}_{{kt}}} + q_{{kt}}^{{(m)}}} \right|} \right)} \right)} \right\}, \\ \end{gathered} $

а ${{r}_{{kt}}} \equiv ({{x}_{k}} - {{x}_{t}})(k,t = 1, \ldots ,{{N}_{e}})$.

При возмущении потенциальной энергии (${{\Phi }_{{kt}}}$ вместо ${{\tilde {\Phi }}_{{kt}}}$) изменение соответствующей свободной энергии имеет порядок ${{{{n}_{e}}\lambda _{e}^{3}} \mathord{\left/ {\vphantom {{{{n}_{e}}\lambda _{e}^{3}} 2}} \right. \kern-0em} 2}\int_V d{{r}_{{12}}}\int dq_{{12}}^{{(1)}}, \ldots ,$ $dq_{{12}}^{{(M - 1)}}$ ${{g}_{{ee}}}({{r}_{{12}}},q_{{12}}^{{(1)}} \ldots q_{{12}}^{{(M - 1)}})$$\left( {{\text{exp(}} - {\kern 1pt} \sum\nolimits_{m = 0}^{M - 1} {\epsilon \Delta {{{\tilde {\Phi }}}_{{12}}}} ) - 1} \right) \sim $ $ \sim {{{{n}_{e}}\lambda _{e}^{3}} \mathord{\left/ {\vphantom {{{{n}_{e}}\lambda _{e}^{3}} {2{{g}_{{ee}}}}}} \right. \kern-0em} {2{{g}_{{ee}}}}}(({{n}_{e}}\lambda _{e}^{3}{{)}^{{ - 1/3}}})$$( - \beta \Delta \tilde {\Phi }(({{n}_{e}}\lambda _{e}^{3}{{)}^{{ - 1/3}}})) \sim {{{{n}_{e}}\lambda _{e}^{3}} \mathord{\left/ {\vphantom {{{{n}_{e}}\lambda _{e}^{3}} {2{{g}_{{ee}}}}}} \right. \kern-0em} {2{{g}_{{ee}}}}}$ × $ \times \,(({{n}_{e}}\lambda _{e}^{3}{{)}^{{ - 1/3}}})({{\beta {{e}^{2}}} \mathord{\left/ {\vphantom {{\beta {{e}^{2}}} {{{\lambda }_{e}}}}} \right. \kern-0em} {{{\lambda }_{e}}}})$ и им можно пренебречь для небольших $\Delta \tilde {\Phi }$ и ${{n}_{e}}\lambda _{e}^{3} \lesssim 1$. Здесь ${{\beta {{e}^{2}}} \mathord{\left/ {\vphantom {{\beta {{e}^{2}}} {{{\lambda }_{e}}}}} \right. \kern-0em} {{{\lambda }_{e}}}} = \sqrt {{{\beta {\text{Ry}}} \mathord{\left/ {\vphantom {{\beta {\text{Ry}}} \pi }} \right. \kern-0em} \pi }} $ и ${{g}_{{ee}}}$ — парная функция распределения электронов. Эта оценка может быть получена в рамках обобщения техники разложения Майера для неидеальных систем частиц с помощью алгебраического подхода, развитого в [44, 45].

Приближение (3) воспроизводит оба предела сильно вырожденной и невырожденной систем фермионов. В классическом пределе из-за множителя $\exp [{{ - \pi {{{\left| {Px - x} \right|}}^{2}}} \mathord{\left/ {\vphantom {{ - \pi {{{\left| {Px - x} \right|}}^{2}}} M}} \right. \kern-0em} M}]$ основной вклад вносит тождественная перестановка, и разности потенциальных энергий в экспоненте уравнения (3) равны нулю. В то же время для сильно вырожденной плазмы, где тепловая длина волны частицы больше среднего межчастичного расстояния и траектории сильно запутаны, потенциальная энергия в уравнении (3) слабо зависит от перестановки, что позволяет заменить все перестановки $P$ на тождественную $E$, как в [33, 34, 40, 41].

При промежуточных вырождениях детерминант в уравнении (3) учитывает интерференционные эффекты кулоновского и обменного взаимодействий электронов, а также позволяет снизить влияние “проблемы знаков” за счет использования прямых методов вычисления определителей. Заметим, что при умеренном вырождении (${{n}_{e}}\lambda _{e}^{3} \lesssim 1$) из-за ферми-отталкивания основной вклад в обменный детерминант дают парные перестановки [33, 34]. Следует подчеркнуть, что элементы обменной матрицы ${{\phi }_{{kt}}}$ и приближенная матрица плотности (3) совпадают с точными аналогичными функциями для двух фермионов:

(5)

$\begin{gathered} \rho ({{x}_{1}},{{x}_{2}}) \equiv \int d{{q}^{{(1)}}} \ldots d{{q}^{{(M - 1)}}} \times \\ \times \,\,\left\{ {\exp \left[ { - \sum\limits_{m = 0}^{M - 1} \left[ {\pi {{{\left| {{{\eta }^{{(m)}}}} \right|}}^{2}} + \epsilon \Phi \left( {\left| {{{r}_{{12}}} + q_{{12}}^{{(m)}}} \right|} \right)} \right]} \right]} \right. - \\ - \,\,\exp \left[ { - \pi \frac{{{{{\left| {{{r}_{{12}}}} \right|}}^{2}}}}{M} - \pi \frac{{{{{\left| {{{r}_{{21}}}} \right|}}^{2}}}}{M} - \sum\limits_{m = 0}^{M - 1} } \right. \times \\ \times \,\,\left[ {\pi {{{\left| {{{\eta }^{{(m)}}}} \right|}}^{2}} + \frac{\epsilon }{2}\Phi \left( {\left| {{{r}_{{21}}}\frac{{2m}}{M} + {{r}_{{12}}} + q_{{12}}^{{(m)}}} \right|} \right) + } \right. \\ \left. {\left. { + \,\,\left[ {\frac{\epsilon }{2}\Phi \left( {\left| {{{r}_{{12}}}\frac{{2{{m}^{{^{{}}}}}}}{{{{M}_{{}}}}} + {{r}_{{21}}} + q_{{21}}^{{(m)}}} \right|} \right)} \right]} \right]} \right\} \equiv \\ \equiv \int d{{q}^{{(1)}}} \ldots d{{q}^{{(M - 1)}}}\exp \left\{ { - \sum\limits_{m = 0}^{M - 1} \pi {{{\left| {{{\eta }^{{(m)}}}} \right|}}^{2}}} \right\} \times \\ \times \,\,\left\{ {\exp \left\{ { - \sum\limits_{m = 0}^{M - 1} \varepsilon \Phi \left( {{{r}_{{12}}} + q_{{12}}^{{(m)}}} \right)} \right\} - \exp \left\{ { - \frac{\pi }{M}2r_{{12}}^{2}} \right\} \times } \right. \\ \left. { \times \,\,\exp \left\{ { - \sum\limits_{m = 0}^{M - 1} \varepsilon \Phi \left[ {{{r}_{{12}}}\left( {1 - {{2m} \mathord{\left/ {\vphantom {{2m} M}} \right. \kern-0em} M}} \right) + q_{{12}}^{{(m)}}} \right]} \right\}} \right\}, \\ \end{gathered} $

где $q_{{12}}^{{(m)}} = q_{1}^{{(m)}} - q_{2}^{{(m)}}$, ${\text{det}}\left\| {\psi (x)} \right\| = 1 - {{\phi }_{{12}}}{{\phi }_{{21}}}$.

Более того, для произвольного числа фермионов в условиях преобладания парного обменного взаимодействия (${{n}_{e}}\lambda _{e}^{3} \lesssim 1$) разложение Лапласа обменного детерминанта по минорам и парным обменным блокам показывает, что обменный детерминант в уравнении (3) стремится к произведению соответствующих точных двухфермионных детерминантов $\prod\nolimits_ (1 - {{\phi }_{{kt}}}{{\phi }_{{tk}}})$, см. уравнение (5). Это подтверждает надежность приближения матрицы плотности (3) и позволяет проводить уверенные расчеты термодинамических величин.

Рассмотрим особенности поведения обменного детерминанта для ${{N}_{e}} = 2$ и $M = 2$ на малых межчастичных расстояниях. Основной вклад в разность потенциальной энергии двух одинаковых частиц определяется первым дифференциалом в ${{\phi }_{{kt}}}$. После сокращения подобных членов имеем

$\begin{gathered} \Phi \left( {\left| {\left( {{{{\tilde {r}}}_{{12}}} + \left( {{{r}_{{12}}} + q_{{12}}^{{(1)}}} \right)} \right)} \right. - \Phi \left( {\left| {\left( {{{r}_{{12}}} + q_{{12}}^{{(1)}}} \right)} \right|} \right)} \right. \approx \\ \approx \left\langle {\nabla \Phi \left( {\left| {\left( {{{r}_{{12}}} + q_{{12}}^{{(1)}}} \right)} \right|({{{\tilde {r}}}_{{12}}})} \right.} \right\rangle = \\ = \left| {\nabla \Phi \left( {\left| {\left( {{{r}_{{12}}} + q_{{12}}^{{(1)}}} \right)} \right|} \right) \cdot \left| {{{{\tilde {r}}}_{{12}}}} \right|\cos {\kern 1pt} {{\theta }_{1}}} \right., \\ \Phi \left( {\left| {\left( {{{{\tilde {r}}}_{{21}}} + \left( {{{r}_{{21}}} + q_{{21}}^{{(1)}}} \right)} \right)} \right.} \right. - \Phi \left( {\left| {\left( {{{r}_{{21}}} + q_{{21}}^{{(1)}}} \right)} \right|} \right) \approx \\ \approx \left\langle {\nabla \Phi \left( {\left| {\left( {{{r}_{{21}}} + q_{{21}}^{{(1)}}} \right)} \right|} \right.({{{\tilde {r}}}_{{21}}})} \right\rangle = \\ = \left| {\nabla \Phi \left( {\left| {\left( {{{r}_{{21}}} + q_{{21}}^{{(1)}}} \right)} \right|} \right) \cdot \left| {{{{\tilde {r}}}_{{21}}}} \right|\cos {\kern 1pt} {{\theta }_{2}}} \right.. \\ \end{gathered} $

Здесь угловые скобки $\left\langle {\,|\,} \right\rangle $ обозначают скалярное произведение двух векторов. Из-за косинуса угла между векторами ${{\tilde {r}}_{{kt}}} = {{r}_{{kt}}}$ и $\nabla \Phi \left( {\left| {({{r}_{{kt}}} + q_{{kt}}^{{(1)}}} \right|} \right)$ определитель является осциллирующей функцией и можно ожидать осцилляций в поведении парных функций распределения на малых межчастичных расстояниях. Для больших межчастичных расстояний $|{{r}_{{kt}}}|$ потенциал взаимодействия приближается к нулю и детерминант является монотонной функцией.

РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

В двухкомпонентной плазме (two component plasma – TCP) все заряды из-за взаимодействия коррелированы, в то время как в модели UEG положительные частицы должны быть некоррелированными, чтобы моделировать нейтрализующий фон [42, 46]. Далее плотность электронов характеризуется параметром Бракнера ${{r}_{s}} = {a \mathord{\left/ {\vphantom {a {{{a}_{{\text{B}}}}}}} \right. \kern-0em} {{{a}_{{\text{B}}}}}}$, где $a = [{3 \mathord{\left/ {\vphantom {3 {(4\pi {{n}_{e}})}}} \right. \kern-0em} {(4\pi {{n}_{e}})}}{{]}^{{1/3}}}$, ${{n}_{e}}$ – плотность электронов, а ${{a}_{{\text{B}}}}$ – радиус Бора. В проведенных расчетах методом MFPIMC использовался алгоритм Метрополиса [11, 40]. Для проверки сходимости расчетов варьировалось максимально доступное число частиц (в диапазоне ${{N}_{e}} = 30, \ldots ,50$) и количество высокотемпературных сомножителей (в диапазоне $M = 20{\text{--}}30$). Для анализа влияния межчастичного взаимодействия в обменной матрице проведено сравнение внутренней энергии и парных функций распределения (pair distribution functions – PDF) для матрицы $\left\| {{{\phi }_{{kt}}}} \right\|$ (4) (далее используется обозначение A1) и приближения ${{\tilde {\phi }}_{{kt}}} \approx {{{\text{e}}}^{{ - \pi {{{\left| {{{x}_{k}} - {{x}_{t}}} \right|}}^{2}}}}}$ (далее называется A2), которое часто используется в литературе.

ПАРНЫЕ ФУНКЦИИ РАСПРЕДЕЛЕНИЯ

Начнем с физического анализа пространственного расположения электронов и положительных частиц. Парные функции распределения ${{g}_{{ab}}}(R)$ определяются как

где $a$ и $b$ – типы кулоновских частиц ($a = e,p$; $b = e,p$). PDF зависит только от разности координат частиц $r = \left| {{{{\mathbf{R}}}_{1}} - {{{\mathbf{R}}}_{2}}} \right|$ из-за трансляционной инвариантности системы.

На рис. 1 приведены PDF для UEG и TCP при ${{r}_{s}} = 2$ и $T = 12{\text{Ry}}$. В невзаимодействующей классической системе зарядов ${{g}_{{ep}}}(R) = {{g}_{{pp}}} \equiv 1$, тогда как межчастичные взаимодействия и квантовая статистика приводят к перераспределению частиц (см. линии 1–4 на рис. 1). Для TCP кулоновское притяжение приводит к возрастанию ${{g}_{{ep}}}(r)$ на малых расстояниях, в то время как кулоновское отталкивание, наоборот, уменьшает ${{g}_{{pp}}}(r)$. Благодаря экранировке обе функции на больших расстояниях стремятся к единице. В UEG жесткий нейтрализующий фон моделируется как идеальный газ некоррелированных классических положительных зарядов, равномерно распределенных в пространстве, поэтому ${{g}_{{pp}}}$ и ${{g}_{{ep}}}$ тождественно равны единице, что демонстрируют линии 3 и 4 на рис. 1a. Для TCP линии 1 и 2 демонстрируют ${{g}_{{ee}}}$, рассчитанные согласно (3) для матрицы $\left\| {{{\phi }_{{kt}}}} \right\|$ (A1), в то время как линии 3 и 4 представляют ${{g}_{{ee}}}$, соответствующие приближению $\left\| {{{{\tilde {\phi }}}_{{kt}}}} \right\|$ (A2).

Отметим различие в поведении ${{g}_{{ee}}}$, полученных для расчетов в случае A1 и его приближения A2. В отличие от приближения A2, выражение A1 более точно учитывает кулоновское и ферми-отталкивания электронов, что и приводит к появлению пика на PDF ${{g}_{{ee}}}$, который отсутствует в расчетах A2 (ср. линии 1 и 3 на рис. 1). PDF ${{g}_{{ee}}}$ в случае A2 (линии 3) ведет себя стандартным известным в литературе образом.

Для электронов с антипараллельными спинами для расчетов в случаях A1 и A2 этот эффект отсутствует (линии 2 и 4 на рис. 1б) и поведение PDF определяется эффективным потенциалом Кельбга, не имеющим кулоновской особенности в нуле за счет учета квантовых эффектов, в частности учета туннелирования на расстояниях порядка тепловой длины волны электрона. Естественно, что линии 2 и 4 практически совпадают. Небольшие колебания PDF вызваны статистической ошибкой метода Монте-Карло.

Рассмотрим физические причины и особенности поведения ${{g}_{{ee}}}$ на расстояниях порядка тепловой длины волны электрона (рис. 1б), которые возникают за счет интерференционных эффектов кулоновского и обменного взаимодействий электронов, учтенных в детерминанте в уравнении (3). Наблюдаемое в рачетах MFPIMC упорядочение электронов в кулоновских системах вызвано эффектом исключенного объема обменных “дырок”, которые возникают за счет ферми-отталкивания и кулоновского взаимодействия электронов. Сопоставление PDF ${{g}_{{ee}}}$ и ${{g}_{{ep}}}$ позволяет сделать вывод о формировании на расстояниях порядка длины волны электрона компактных трехчастичных комплексов, состоящих из двух электронов и положительно заряженной частицы. В результате в “дырке” в среднем на малых расстояниях от электрона может преобладать положительный заряд, который за счет присутствия второго электрона будет скомпенсирован при увеличении этого расстояния до длины волны электрона. Для проверки этого утверждения рассмотрим относительные вероятности найти две частицы сорта $a$ и $b$ на расстоянии $r$ друг от друга. Эти вероятности, пропорциональные ${{r}^{2}}{{g}_{{ab}}}(r)$, показаны на рис. 2 в логарифмическом масштабе для $T = 12{\text{Ry}}$ и $T = 6{\text{Ry}}$ (${{r}_{s}} = 2$). Оценить заряд такого комплекса можно с помощью интегралов от относительных вероятностей нахождения двух частиц сорта $a$ и $b$ на расстоянии $r$ друг от друга $\left( { \sim {{r}^{2}}{{g}_{{ab}}}(r)} \right)$ по формуле ${{Z}_{{{\text{tot}}}}}(\tilde {r}) \sim \int_0^{\tilde {r}} {4\pi {{r}^{2}}} ({{g}_{{ep}}}(r) - {{g}_{{ee}}}(r))dr$ при $\tilde {r} \sim {{\lambda }_{e}}$. На рис. 2 линиями 1–4 показаны относительные вероятности для одноименно и разноименно заряженных частиц. Поведение линии 1 на рис. 3 позволяет понять распределение суммарного заряда в трехчастичном комплексе и подтверждает его электронейтральность, так как полный осредненный заряд “дырки”, создаваемой парой электронов и присутствием в ней положительно заряженной частицы, стремится к нулю при $\tilde {r} \to {{\lambda }_{e}}$. Это согласуется с теоретическим результатом, согласно которому заряд “обменной дыры” вокруг электрона равен по модулю и противоположен по знаку заряду электрона [47]. Таким образом, наблюдаемое в расчетах MFPIMC упорядочение электронов связано с формированием квазинейтральных обменно-корреляционных экситонов в “море ферми-электронов”. Возникновение нейтральных обменно-корреляционных экситонов обсуждалось в [37–39], где в качестве причины этого явления также рассматривался эффект исключенного объема [48] и взаимодействие электронов с положительно заряженной обменно-корреляционной дыркой, появляющейся благодаря ферми-отталкиванию электронов на расстояниях порядка длины волны.

Сопоставление рассчитанных по MFPIMC относительных вероятностей нахождения двух частиц сорта $a$ и $b$ на расстоянии $r$ друг от друга для TCP и UEG при температуре $T = 6{\text{Ry}}$ проведено на рис. 4. Необходимо отметить очень хорошее совпадение рассчитанных вероятностей для системы электронов с одинаковым и противоположным спинaми, несмотря на абсолютно различный характер положительных зарядов, нейтрализующих кулоновскую систему частиц, и соответствующих парных функций распределения.

ВНУТРЕННЯЯ ЭНЕРГИЯ

Расчет внутренней энергии производился с помощью выведенного ранее выражения [42, 46]. При ${{r}_{s}} = 2$ и ${T \mathord{\left/ {\vphantom {T {{\text{Ry}}}}} \right. \kern-0em} {{\text{Ry}}}} = 6$ имеем следующие результаты для внутренней энергии согласно формуле (3) с выражением A1 для $\psi $:

${\text{энергия}}\,\,{\text{TCP}}{{{{E}_{{{\text{tcp}}}}}} \mathord{\left/ {\vphantom {{{{E}_{{{\text{tcp}}}}}} {({{k}_{{\text{B}}}}T)}}} \right. \kern-0em} {({{k}_{{\text{B}}}}T)}} = 1.419,$

${{{\text{UEG}}\,\,{{E}_{{{\text{ueg}}}}}} \mathord{\left/ {\vphantom {{{\text{UEG}}\,\,{{E}_{{{\text{ueg}}}}}} {({{k}_{{\text{B}}}}T)}}} \right. \kern-0em} {({{k}_{{\text{B}}}}T)}} = 1.356$.

Для внутренней энергии, рассчитанной по формуле (3) с приближением A2 для $\psi $, получено:

Внутренняя энергия для TCP и UEG оказывается ниже для обменной матрицы $\psi $ в приближении A2 по сравнению с A1. Это связано с тем, что PDF для выражения A1 выше, чем в приближении A2.

Разница заметная, но не очень большая для интегральных термодинамических характеристик, поскольку произведение ${{r}^{2}}{{g}_{{ab}}}(r)$ уменьшает удельные вклады в интегралы на малых межчастичных расстояниях.

ЗАКЛЮЧЕНИЕ

Приближение A2 справедливо для сильно вырожденной системы. В этом случае определитель $\psi (x)$ эквивалентен определителю Грама [43] ${\text{det}}\left\| {\left\langle {\left. {\chi ({{x}_{t}})} \right|\chi ({{x}_{k}})} \right\rangle } \right\|$ для линейно независимой системы векторов ${{N}_{e}}\left| {{{\chi }_{p}}({{x}_{k}})} \right\rangle = \left| {{{{\text{e}}}^{{i\langle p|{{x}_{k}}\rangle /\hbar }}}} \right\rangle $ со скалярным произведением 〈χ(x_t)|χ(x_k)〉 = $\int dp{{e}^{{{\text{i}}\langle p|{{x}_{k}} - {{x}_{t}}\rangle /\hbar }}} \times $ $ \times \,\,{{{{{\text{e}}}^{{ - \beta {{p}^{2}}/2m}}}} \mathord{\left/ {\vphantom {{{{{\text{e}}}^{{ - \beta {{p}^{2}}/2m}}}} {\int dp{{{\text{e}}}^{{ - \beta {{p}^{2}}/2m}}}}}} \right. \kern-0em} {\int dp{{{\text{e}}}^{{ - \beta {{p}^{2}}/2m}}}}} = {{{\text{e}}}^{{ - \pi {{{\left| {{{r}_{{kt}}}} \right|}}^{2}}/\lambda _{e}^{2}}}}$. Определитель Грама всегда неотрицателен [43] и позволяет решить “проблему знака” в PIMC-моделировании фермионов в данном приближении.

Для сравнения рассмотрим аналогичный матричный элемент ${{{\text{e}}}^{{ - \pi {{{\left| {x_{k}^{{(M)}} - x_{t}^{{(0)}}} \right|}}^{2}}}}}$ в обменном детерминанте $\det \Psi (x)$ уравнения (1), который обычно используется в стандартном PIMC. Отметим, что матричные элементы в $\psi (x)$ имеют другую структуру по сравнению с матричными элементами в $\Psi (x)$, которые определяются аналогичным скалярным произведением, но для другого набора векторов $2{{N}_{e}}\left| {{{{\text{e}}}^{{{\text{i}}\langle p|x_{k}^{{(M)}}\rangle /\hbar }}}} \right\rangle $ и $\left| {{{{\text{e}}}^{{{\text{i}}\langle p|x_{k}^{{(0)}}\rangle /\hbar }}}} \right\rangle $. Обменный детерминант $\det \Psi (x)$ – знакопеременная функция, порождающая “проблему знака” в уравнении (1), что приводит к снижению точности расчетов PIMC.

Таким образом, модифицированный метод MFPIMC демонстрирует квантовое упорядочение электронов на малых расстояниях, возникающее от кулоновского взаимодействия электронов с обменными “дырками”. Это явление известно в литературе [37, 38] и связано с образованием обменных экситонов, модифицирующих парные функции распределения.

Авторы признательны профессору М. Боницу за плодотворные обсуждения. Теоретический подход, вывод основных уравнений и алгоритмическая реализация метода MFPIMC поддержаны Российским научным фондом (грант № 20-42-04421). Выражаем благодарность ЦКП “Суперкомпьютерный центр ОИВТ РАН”, МСЦ РАН и ЦКП “Дальневосточный вычислительный ресурс” за предоставленное компьютерное время.