Nonlocal effects in the method of integral equations of the theory of liquids

Abstract

The Ornstein-Zernike integral equations for thermodynamically equilibrium liquids are considered, considering irreducible diagrams. An algorithm for calculating the first irreducible diagram is formulated for a spatially homogeneous fluid. For liquids bordering a solid surface, an equation is obtained that implicitly considers all irreducible diagrams.

Full Text

Введение

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

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

Постановка задачи

Статистическая физика термодинамически равновесных жидкостей базируется на цепочке уравнений Боголюбова—Борна—Грина—Кирквуда—Ивона для частичных функций распределения [1—5], либо эквивалентной ей цепочке уравнений Орнштейна—Цернике для одночастичной F1 и двухчастичной F12 функций распределения:

  $\begin{array}{c}{\omega }_1=n{\displaystyle \int F_2{C}_{12}^{\left(1\right)}d\left(2\right)+\mu h_{12}}=\\ =C_{12}^{\left(2\right)}+n{\displaystyle \int F_3{C}_{13}^{\left(2\right)}{h}_{23}d\left(3\right)}.\end{array}$  (1)

Заданными параметрами в (1) являются потенциальная энергия частицы во внешнем поле Φ(ri), энергия межмолекулярного взаимодействия $\Phi \left(\left|r_i-r_j\right|\right)$, температура T и численная плотность n, где σ — характерный размер молекулы. Численная плотность изменяется в пределах $\text{0}\le n=\frac{N}{V}{\sigma }^3\le \approx \mathrm{1,}$ что соответствует изменению концентрации от идеального газа до плотной жидкости. Для систем, межмолекулярное взаимодействие в которых задано парными потенциалами, наиболее важными являются одночастичная $F_1=F_1\left(r_1\right)=\exp (-{\Phi }_1/kT+{\omega }_1)$ и двухчастичная $F_2=F_{12}\left({\overrightarrow{r}}_1,{\overrightarrow{r}}_2\right)=\exp (-{\Phi }_{12}/kT+{\omega }_{12})$ функции распределения. Знание этих функций дает возможность вычислить структурные и термодинамические характеристики вещества. В свою очередь, в (1) входят парная $h_{ij}\left({\overrightarrow{r}}_i,{\overrightarrow{r}}_j\right)=F_{ij}\left({\overrightarrow{r}}_i,{\overrightarrow{r}}_j\right)-1$ и прямые $C_{ij}^{\left(k\right)}\left({\overrightarrow{r}}_i,{\overrightarrow{r}}_j\right)$ корреляционные функции. Интегрирование ведется по координатам i-й частицы $d\left(i\right)=d{\overrightarrow{r}}_i$; где n — плотность; Gi = exp(–Φi / kT + ωi) — одночастичная функция распределения, описывающая положение частицы в лабораторной системе координат; Φi –потенциальная энергия частицы во внешнем поле; ωi — одночастичный термический потенциал; µ — химический потенциал, определяемый из условия перехода к пространственно-однородной системе; hij = [exp(–Φij / kT + ωij) — 1] — парная корреляционная функция, связанная с двухчастичной функцией распределения соотношением Gij = GiGj,(1 + hij); ωij — двухчастичный термический потенциал, учитывающий опосредованное взаимодействие двух частиц через их окружение; $C_{ij}^{\left(k\right)}$ — прямые корреляционные функции:

 $C_{ij}^{\left(1\right)}=h_{ij}-{\omega }_{ij}-\frac{1}{2}{h}_{ij}\left({\omega }_{ij}+B_{ij}^{\left(1\right)}\right)$,            (2)

 $C_{ij}^{\left(2\right)}=h_{ij}-{\omega }_{ij}+B_{ij}^{\left(2\right)}$.    (3)

Здесь $B_{ij}^{\left(1\right)}\left(r_{ij}\right)$, $B_{ij}^{\left(2\right)}\left(r_{ij}\right)$ — бридж-функционалы (бесконечные функциональные ряды неприводимых диаграмм):

    (4)

 

   (5)

 

Каждая диаграмма является многомерным интегралом. Линии соответствуют парным корреляционным функциям hij(rij). Черные точки означают умножение на nGi с интегрированием по соответствующей координате $d{\overrightarrow{r}}_i$. Светлые точки соответствуют аргументам rij. Таким образом, значения бридж-функционалов зависят не только от аргумента rij, но и от распределения плотности в окрестности точек ${\overrightarrow{r}}_i,{\overrightarrow{r}}_j$. При решении конкретных задач нелокальные бридж-функционалы заменяют локальной функцией $B_{ij}^{\left(1\right)}\left(r_{12}\right)=B_{ij}^{\left(1\right)}\left(h\left(r_{12}\right)\right)$, $B_{ij}^{\left(2\right)}\left(r_{12}\right)=B_{ij}^{\left(2\right)}\left(h\left(r_{12}\right)\right)$. Логическое обоснование такой замены отсутствует, тем не менее такой подход до сих пор активно используется и приводит к приближенным интегральным уравнения для парной корреляционной функции h(r12).

Интегральные уравнения для пространственно-однородных систем. Локальное приближение

Важное значение имеют пространственно-однородные изотропные среды (объемные жидкости в отсутствии внешних полей и вдали от ограничивающих поверхностей). В этом случае F1 = 1, и уравнение для одночастичной функции распределения сводится к определению избыточного химического потенциала µ. Уравнение (3) определяет парную корреляционную функцию h12(r12), зависящую от расстояния между центрами двух частиц $r_{12}=\left|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2\right|$. При вычислении бридж-функций в этих рядах оставляют слагаемые, которые сводятся к интегралам типа свертки.

В результате получается связь между прямой и парной корреляционной функцией (замыкание). Это сводит соотношение ОЦ к приближенным нелинейным интегральным уравнениям, которые решаются численно. Наиболее известными из них являются гиперцепное, Перкус—Йевика, Роджерса—Янга, Мартынова—Саркисова [1—5]. Аналитическое решение нелинейных уравнений получено только для системы твердых сфер в аппроксимации Перкус—Йевика [13]. Как правило, для решения нелинейных уравнений приходится прибегать к численным методам. Основной вывод по результатам численных расчетов сводится к тому, что любое из перечисленных замыканий дает лучшие результаты по сравнению с другими, в зависимости от вида потенциала межмолекулярного взаимодействия, температуры и плотности [1—5]. Таким образом, на данный момент нет однозначных физических критериев, позволяющих отдать предпочтение какому-либо приближению (замыканию). Оценка погрешности решения приближенных нелинейных уравнений проводится сравнением с данными численного эксперимента, являющимися эталоном точности.

Попытки получить все более точные интегральные уравнения [1—5] в рамках замены бридж-функционалами бридж-функциями, так и не привели к каким-либо значимым результатам.

Пространственно-однородные системы. Нелокальные эффекты

Рассмотрим первую неприводимую диаграмму в (5):

      (6)

 

Выражая парную корреляционную функцию h(r34) через ее фурье-компоненту hχ с помощью обратного фурье-преобразования получим

 $\begin{array}{c}{I}_{12}\left(r_{12}\right)={\displaystyle \iint d{\overrightarrow{r}}_{13}d{\overrightarrow{r}}_{14}h\left(r_{13}\right)h\left(r_{23}\right)h\left(r_{14}\right)h\left(r_{24}\right)4\pi }\times \\ \times {{\displaystyle \int }}_0^{\infty }\frac{{\chi }^{2}d\chi }{{(2\pi )}^3}{h}_{\chi }\frac{\sin \chi r_{34}}{\chi r_{34}}.\end{array}$            (7)

При вычислении внутреннего интеграла применим теорему сложения для цилиндрических функций. Ось z декартовой системы координат направляем вдоль вектора ${\overrightarrow{r}}_{12}$, проходящего через центры частиц 1 и 2:

 $\frac{\sin \chi r_{34}}{\chi r_{34}}={{\displaystyle \sum }}_{\mathcal{l}=0}^{\infty }\left(2\mathcal{l}+1\right)j_{\mathcal{l}}\left(\chi r_{13}\right)\left(\chi r_{14}\right)P_{\mathcal{l}}\left(\cos \cos \gamma \right)$, (8)

 $\begin{array}{c}{P}_{\mathcal{l}}\left(\cos \gamma \right)=P_{\mathcal{l}}\left(\cos {\theta }_{13}\right)P_{\mathcal{l}}\left(\cos {\theta }_{14}\right)+\\ +2{{\displaystyle \sum }}_{m=1}^{\mathcal{l}}\frac{\left(\mathcal{l}-m\right)!}{\left(\mathcal{l}+m\right)!}{P}_{\mathcal{l}}^m\left(\cos {\theta }_{13}\right)P_{\mathcal{l}}^m\left(\cos {\theta }_{14}\right)\times \\ \times \cos m\left({\phi }_{13}-{\phi }_{14}\right),\end{array}$ (9)

где θ13, θ14 и φ13, φ14 — полярные и азимутальные углы, задающие ориентацию векторов ${\overrightarrow{r}}_{13}$ и ${\overrightarrow{r}}_{14}$ соответственно; $P_{\mathcal{l}}\left(\cos \gamma \right)$ — полиномы Лежандра, $P_{\mathcal{l}}^m\left(\cos {\theta }_{ij}\right)$ — присоединенные полиномы Лежандра.

При интегрировании по угловым переменным φ13, φ14, слагаемые, содержащие присоединенные полиномы Лагранжа, обращаются в нуль, что приводит к следующему выражению:

 $I_{12}\left(r_{12}\right)=4\pi {{\displaystyle \int }}_0^{\infty }\frac{{\chi }^{2}d\chi }{{(2\pi )}^3}{h}_{\chi }{{\displaystyle \sum }}_{\mathcal{l}=0}^{\infty }\left(2\mathcal{l}+1\right){\left[I_{\mathcal{l}}\left(r_{12},\chi \right)\right]}^2,$ (10)

  $I_{\mathcal{l}}\left(r_{12},\chi \right)=\int d{\overrightarrow{r}}_{13}h\left(r_{13}\right)h\left(r_{23}\right)j_{\mathcal{l}}\left(\chi r_{13}\right)P_{\mathcal{l}}\left(\cos {\theta }_{13}\right)$. (11)

При вычислении внутреннего интеграла парную корреляционную функцию h(r23) выразим через ее фурье-образ, тем самым:

 $\begin{array}{c}{I}_{\mathcal{l}}\left(r_{12},\chi \right)=4\pi {{\displaystyle \int }}_0^{\infty }\frac{k^{2}dk}{{(2\pi )}^3}\int d{\overrightarrow{r}}_{13}h\left(r_{13}\right)\times \\ \times \frac{\sin k{r}_{23}}{k{r}_{23}}{j}_{\mathcal{l}}\left(\chi r_{13}\right)P_{\mathcal{l}}\left(\cos {\theta }_{13}\right).\end{array}$ (12)

Применяя теорему сложения для цилиндрических функций, аналогично (8) и выполняя интегрирование по угловым переменным, получим

  $\begin{array}{l}{I}_{\mathcal{l}}\left(r_{12},\chi \right)=\frac{2}{\pi }{{\displaystyle \int }}_0^{\infty }{k}^{2}dk{h}_k{j}_{\mathcal{l}}\left(k{r}_{12}\right)\times \\ \times {{\displaystyle \int }}_0^{\infty }{r}_{13}^{2}h\left(r_{13}\right)j_{\mathcal{l}}\left(\chi r_{13}\right)j_{\mathcal{l}}\left(k{r}_{13}\right).\end{array}$           (13)

Меняя местами порядок интегрирования и выражая фурье-образ hχ с помощью прямого преобразования через парную корреляционную функцию, придадим (13) следующий вид:

 $I_{\mathcal{l}}\left(r_{12},\chi \right)=2\pi {{\displaystyle \int }}_0^{\infty }{r}_{13}^{2}d{r}_{13}h\left(r_{13}\right)j_{\mathcal{l}}\left(\chi r_{13}\right)K_{\mathcal{l}}\left(r_{12},r_{13}\right)$, (14) 

    $K_{\mathcal{l}}\left(r_{12},r_{13}\right)=\frac{1}{r_{12}{r}_{13}}{{\displaystyle \int }}_{\left|r_{12}-r_{13}\right|}^{r_{12}+r_{13}}rdrh\left(r\right)P_{\mathcal{l}}\left(\frac{r_{12}^2+r_{13}^2-r^2}{2r_{12}{r}_{13}}\right)$. (15)

Таким образом, получаем окончательное выражение для первой неприводимой диаграммы:

 $\begin{array}{c}{I}_{12}\left(r_{12}\right)=2{\pi }^2{{\displaystyle \sum }}_{\mathcal{l}=0}^{\infty }\left(2\mathcal{l}+1\right){{\displaystyle \int }}_0^{\infty }{r}^{2}drh\left(r\right)K_{\mathcal{l}}\left(r_{12},r\right)\times \\ \times {{\displaystyle \int }}_0^{\infty }{R}^{2}dRh\left(R\right)K_{\mathcal{l}}\left(r_{12},R\right)K_{\mathcal{l}}\left(r,R\right).\end{array}$ (16)

Заметим, что первый член ряда (I = 0) является уточнением аппроксимации Мартынова—Саркисова. Действительно, в этом случае полином Лежандра равен единице и первая неприводимая диаграмма принимает вид

 $\begin{array}{l}{I}_{12}\left(r_{12}\right)={{\displaystyle \int }}_0^{\infty }{r}_{13}^{2}d{r}_{13}h\left(r_{13}\right)K_0\left(r_{12},r_{13}\right)\times \\ \times {{\displaystyle \int }}_0^{\infty }{R}^{2}dRh\left(R\right)K_0\left(r_{12},R\right)K_0\left(r_{13},R\right).\end{array}$ (17)

Делая предположение, что основной вклад в неприводимую диаграмму (6) вносит значение h(r34) = –1 и соответственно в (17) полагая K0(r13, R) = –1, получаем аппроксимацию Мартынова—Саркисова [1]:

 $I_{12}\left(r_{12}\right)=-\frac{1}{2}{\omega }_{12}^2$. (18)

Таким образом, (17) является обобщением аппроксимации Мартынова—Саркисова. Естественно, для полного вычисления неприводимой диаграммы I12(r12) в (16) необходимо учитывать все члены ряда по полиномам Лежандра.

Пространственно-неоднородные системы. Синглетное приближение

В случае пространственно-неоднородных систем (молекулярная система вблизи твердой поверхности) для вычисления микроструктуры вещества и его термодинамических параметров необходимо знать обе функции распределения — одночастичную F1 и двухчастичную F12. Функция F1 зависит только от одной переменной z1 — удаления частицы от поверхности. Однако функция F12 теперь зависит от трех переменных: расстояния между центрами частиц r12 и удаления каждой из них от поверхности — z1 и z2. В синглетном приближении [11] полагают, что F12 зависит только от переменной r12 и система уравнений (2) и (3) распадаются на два независимых уравнения. Второе уравнение системы решается численно для пространственно-однородной системы (граничное условие для F12) и определяет парную корреляционную функцию h12(r12). Первое уравнение системы в зависимости от замыканий между прямой и парной корреляционной функции, приводит к нелинейным интегральным уравнениям для одночастичной функции распределения F1(z1). Обзор численных решений различных синглетных уравнений приведен в работе [11]. Отметим, что аналитическое решение уравнения для одночастичной функции распределения возможно лишь для частных случаев, например для одномерной и двумерной задачи [12].

Модифицированное синглетное уравнение. Суммирование неприводимых диаграмм

В трехмерном случае для молекулярной системы, граничащей с твердой поверхностью, аналитическое решение получить не удается. В подобной ситуации возникло предположение, что можно учитывать бесконечный ряд неприводимых диаграмм в уравнении ОЦ так, чтобы они взаимно компенсировали друг друга и в результате получалось линейное уравнение. Такая идея была реализована в наших работах [14—16], в соответствии с которой мы полагаем:

 $G_1\left(z_2\right)C_{12}^{\left(1\right)}\left(r_{12},z_1,z_{12}\right)=C_{12}^{\left(\mathrm{1,0}\right)}\left(r_{12}\right)+S_{12}^{\left(1\right)}\left(r_{12},z_1,z_{12}\right)$, (19)

где $S_{12}^{\left(1\right)}\left(r_{12},z_1,z_{12}\right)$ учитывает вклад всех неприводимых диаграмм. Однако в работах [14—16] не было учтено, что в химическом потенциале µ, входящем в уравнения ОЦ, также необходимо выделять неприводимые диаграммы. С учетом этого обстоятельства подстановка (19) в уравнения ОЦ с последующим интегрированием в цилиндрической системе приводит к уравнениям

 $\begin{array}{c}\left(e^{{\omega }_1\left(z_1\right)}-1\right)-2\pi n{{\displaystyle \int }}_0^{\infty }d{z}_2\left(e^{{\omega }_1\left(z_2\right)}-1\right)\times \\ \times {{\displaystyle \int }}_{\left|z_{12}\right|}^{\infty }{r}_{12}d{r}_{12}{C}_{12}^{\left(\mathrm{1,0}\right)}\left(r_{12}\right)=\\ =-2\pi n{{\displaystyle \int }}_{z_1}^{\infty }d{z}_{12}{{\displaystyle \int }}_{z_{12}}^{\infty }{r}_{12}d{r}_{12}{C}_{12}^{\left(\mathrm{1,0}\right)}\left(r_{12}\right),\end{array}$ (20)

$\begin{array}{l}2\pi n{{\displaystyle \int }}_0^{\infty }d{z}_2\left(e^{{\omega }_1\left(z_2\right)}-1\right){{\displaystyle \int }}_{\left|z_{12}\right|}^{\infty }{r}_{12}d{r}_{12}{S}_{12}^{\left(1\right)}\left(r_{12},z_1,z_2\right)-\\ -2\pi n{{\displaystyle \int }}_{z_1}^{\infty }d{z}_{12}{{\displaystyle \int }}_{z_{12}}^{\infty }{r}_{12}d{r}_{12}{S}_{12}^{\left(1\right)}\left(r_{12}\right)={\omega }_1\left(z_1\right)-\left(e^{{\omega }_1\left(z_1\right)}-1\right)\end{array}$. (21)

Уравнение (20) является интегральным уравнением Фредгольма второго рода. Его ядро и правая часть зависят только от прямой корреляционной функции пространственно-однородной жидкости. Алгоритм решения уравнения (20) приведен нами в [14—16]. Рассмотрим уравнение (21), которое также должно решаться в синглетном приближении, то есть $S_{12}^{\left(1\right)}\left(r_{12},z_1,z_{12}\right)\to S_{12}^{\left(1\right)}\left(r_{12}\right)$. В результате получается уравнение

$\begin{array}{c}2\pi n{{\displaystyle \int }}_0^{\infty }d{z}_2\left(e^{{\omega }_1\left(z_2\right)}-1\right){{\displaystyle \int }}_{\left|z_{12}\right|}^{\infty }{r}_{12}d{r}_{12}{S}_{12}^{\left(1\right)}\left(r_{12}\right)-\\ -2\pi n{{\displaystyle \int }}_{z_1}^{\infty }d{z}_{12}{{\displaystyle \int }}_{z_{12}}^{\infty }{r}_{12}d{r}_{12}{S}_{12}^{\left(1\right)}\left(r_{12}\right)=\\ ={\omega }_1\left(z_1\right)-\left(e^{{\omega }_1\left(z_1\right)}-1\right).\end{array}$ (22)

Неизвестной функцией является сумма всех неприводимых диаграмм — $S_{12}^{\left(1\right)}\left(r_{12}\right)$. Дифференцированием по переменной z1, уравнение (22) сводится к уравнению Вольтерра первого рода:

$\begin{array}{l}2\pi n{{\displaystyle \int }}_{-z_1}^{\infty }xdx\left(e^{{\omega }_1\left(z_1+x\right)}-1\right)S_{12}^{\left(1\right)}\left(\left|x\right|\right)+\\ +2\pi n{{\displaystyle \int }}_{z_1}^{\infty }xdx{S}_{12}^{\left(1\right)}\left(x\right)={\omega }_1^{\left(1\right)}\left(z_1\right)\left(1-e^{{\omega }_1\left(z_1\right)}\right)\end{array}$, (23)

которое может быть решено стандартными методами.

Заключение

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

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

Авторы признательны Ю. Д. Фомину и В. Н. Рыжову за беседы, способствовавшие написанию данной публикации.

Исследование проведено при частичной финансовой поддержке Российского фонда фундаментальных исследований (проект № 18-02-00523а).