Статистико-механічне моделювання властивостей полярних і магнітних рідин

У вступі обґрунтовано актуальність теми дисертації, сформульовано мету і задачі роботи, відзначено її наукову новизну і практичне значення, визначено особистий внесок здобувача та наведено інформацію щодо апробації одержаних результатів.

Основна частина дисертації розпочинається з огляду літератури у першому розділі, де висвітлено важливі попередні досягнення та описано сучасний стан досліджень динамічних і фазових властивостей простих, полярних і магнітних рідин. Обговорено базисні теоретичні методи та комп'ютерно-адаптовані підходи і наголошено на області їх застосування. Звернуто увагу на недоліки, труднощі і переваги існуючих теорій. Таким чином окреслюється коло проблем, що вимагали детального теоретичного вивчення і які стали предметом розгляду в цій дисертаційній роботі.

У другому розділі досліджуються найпростіші моделі полярних рідин, де електрична взаємодія між частинками задається потенціалом точкових диполів. Слідуючи ідеї скороченого опису Боголюбова, методом НСО Зубарєва [2] отримано узагальнені рівняння переносу

1)

для середніх значень динамічних змінних, що визначають стан системи в час t. У набір включено оператори густини числа частинок , трансляційного і кутового імпульсів, енергії та дипольного моменту . Вектор br, t) визначається локальними значеннями оберненої температури , поступальної vr t) і кутової wr t) гідродинамічних швидкостей, хімічного потенціалу r t) та електричного поля Er t), а r r t t) позначає матрицю функцій пам'яті. Усереднення у правій частині 1) здійснюється за квазірівноважною функцією розподілу , де m і J - це маса і момент інерції частинок, відповідно. Рівняння узагальненої гідродинаміки 1) є формально точними і можуть описувати як слабо-, так і сильнонерівноважні процеси.

У випадку малих відхилень від рівноваги рівняння 1) спрощуються та записуються в просторі хвильового вектора k і частоти для рівноважних часових кореляційних функцій за відсутності зовнішнього електричного поля як

iIk) k Fk Fk 2)

де - статичні функції, , а L - оператор Ліувілля системи. Для розв'язання 2) було застосовано метод УКМ, згідно якого у набір динамічних змінних, окрім базисних величин , включаються також їхні кінетичні складові вищих порядків Lak), , …, . На розширеному наборі можна застосувати марківське наближення k k 0 і звести задачу до знаходження власних значень колективних мод) узагальненого гідродинамічного оператора Gk)k) k 0. Тоді кожний елемент матриці Fk t виражається через суму експоненційних доданків , що пов'язані з існуванням певної моди . Більше того, гідродинамічні та дипольні флуктуації для слабонерівноважних станів стають відносно незалежними і можуть досліджуватись окремо. Зокрема, для дипольних флуктуацій, коли ak)pk), у частотному представленні отримано аналітичний вираз для динамічного фактора Кірквуда у формі ланцюгового дробу

3)

де - елементи частотної матриці k), які залежать лише від статичних кореляційних функцій різних порядків s12S, а - елемент матриці функцій пам'яті, що визначається кореляційним часом . Останній також можна виразити через статичні функції, використовуючи співвідношення при достатньо великому числі мод S. Знаючи f k , повздовжню і поперечну компоненти діелектричної проникності можна обчислити як . Подібним чином було отримано аналітичні вирази для узагальнених коефіцієнтів переносу, таких як теплопровідність та зсувна і об'ємна в'язкості.

Рис. 1. Спектр колективних дипольних збуджень рідини Штокмайєра для поперечних a) і повздовжніх b) флуктуацій у чотиримодовому S=4) описі. Дифузійні моди, дійсна і уявна частини пропагаторних мод позначені літерами D, P і W, відповідно. Результати для S=2 показано тонкими кривими.

Як ілюстрацію дієвості запропонованого підходу було розглянуто процес діелектричної релаксації для моделі Штокмайєра. Числові розрахунки виконувалися при густині і температурі, що відповідають термодинамічному стану води при нормальних умовах. Спектр колективних дипольних збуджень, отриманий у дво- і чотиримодовому S2 і 4) описах, зображено на Рис. 1. У випадку поперечних флуктуацій можна чітко ідентифікувати дебаївську моду , яка добре відокремлена від усіх інших збуджень при малих значеннях хвильового вектора k. Мода зумовлена дифузійним механізмом релаксації і в границі k0 може бути записана через коефіцієнти обертової і трансляційної дифузій як , де NV - густина числа частинок, а ck) - двочастинкова пряма кореляційна функція. Для повздовжніх флуктуацій можна легко розрізнити дві пропагаторні моди , які проявляють квазічастинкову диполяронну) поведінку при малих k, оскільки , де і визначають частоту і згасання колективних збуджень, відповідно. Диполяронні моди якісно відтворюються вже в рамках формалізму двох змінних S2), який окрім дифузійних процесів враховує і ефекти вільного руху. Якщо дипольні взаємодії вважати слабкими, то у гідродинамічній границі знайдемо , де d - дипольний момент частинки, - її розмір, , а - інтенсивність леннард-джонсівського потенціалу.

Рис. 2. Часові кореляційні функції моделі Штокмайєра для двох різних значень хвильового вектора ).

Дані МД для повздовжніх і поперечних компонент зображені кружечками і квадратиками, відповідно. Результати дво-, три- і п'ятимодового наближень показано довго-, короткоштриховими і суцільними кривими.

Інші моди можна визначити при S?2, коли додатково враховуються взаємодія між диполями та кінетичні процеси. Вони є важливими для кількісного опису поперечних і повздовжніх флуктуацій, особливо у проміжній області k. Для прикладу на Рис. 2 приведено нормовані часові кореляційні функції k t f k t f k , отримані при S2,3 і 5, що порівнюються з даними МД розділ 4). Як бачимо, монотонне згасання поперечних флуктуацій задовільно відтворюється вже при S2. Повздовжні флуктуації, де спостерігаються сильні диполяронні осциляції з повільним згасанням, описуються у цьому випадку лише якісно. Зі збільшенням числа мод до S5 наближені i МД дані практично співпадають.

Схожу швидку збіжність результатів з ростом числа мод встановлено і для повздовжньої та поперечної компонент узагальненої діелектричної проникності. Зокрема, дійсну і уявну частини функції зображено на Рис. 3. Видно, що при S5 можна говорити про кількісне відтворення діелектричної проникності в усій області спектру.

Рис. 3. Дійсна a) і c)) та уявна b) і d)) частини повздовжньої компоненти діелектричної проникності рідини Штокмайєра в залежності від частоти для двох різних значень хвильового вектора. Результати дво-, три- і п'ятимодового наближень показано довго-, короткоштриховими і суцільними кривими, відповідно.

Предметом розгляду третього розділу є нерівноважні і діелектричні властивості більш реалістичних моделей полярних рідин типу взаємодіючих силових центрів, що враховують атомну структуру молекул. На відміну від дипольних систем таким рідинам притаманні складніші динамічні процеси, які відбуваються на суттєво різних часових масштабах і зумовлені появою кулонівських взаємодій. Як наслідок, застосування стандартного підходу УКМ у рамках звичайного марківського наближенням призводить до значних труднощів, пов'язаних зі збіжністю результатів [11,15]. Щоб обійти ці труднощі, запропоновано модифікацію методу УКМ, яка базується на математично строгому формалізмі і дає змогу послідовно врахувати немарківські ефекти при моделюванні кінетичних ядер пам'яті.

Основна ідея модифікованого підходу полягає у включенні до набору динамічних величин ak), таких як густини числа частинок, імпульсу, енергії, поляризації і їх вищі кінетичні складові, також додаткових, так званих ультраповільних квазізмінних , що являють собою інтеграли по часу від базисних величин. Наприклад, у випадку поляризаційних флуктуацій, коли ak)Pk), для розширеного набору отримаємо

4)

де - вектор поляризації. Тут позначає положення атома a, що несе заряд і належить молекулі i з центром мас в , а N і M - загальне число молекул і силових центрів у кожній молекулі, відповідно. Вектор Pk) відрізняється від густини дипольного моменту і прямує до останнього лише в довгохвильовій границі, коли , де - дипольний момент i -ї молекули. Складнішими, ніж у випадку дипольних моделей, будуть тут і вирази для гідродинамічних величин.

Використання модифікованого набору 4) веде до перенормування системи рівнянь 2) для часових кореляційних функцій і власних значень узагальненого гідродинамічного оператора. У результаті знайдемо, що часову кореляційну функцію поляризаційних флуктуацій можна записати як суму S S експонент з певними ваговими множниками, які асоціюються з вкладом відповідної моди , де 12S S . У частотному зображенні f k t приймає форму ланцюгового дробу

5)

з ядром пам'яті у виді Паде-апроксимації

6)

точного кінетичного ядра пам'яті у наближенні S -го порядку. Бачимо, що при S 0 функція явно залежить від частоти, а отже може описувати суттєво немарківські ефекти. Для S 0 модифікований підхід зводиться до стандартного формулювання методу УКМ див. розділ 2). Коефіцієнти і , які з'являються у рівнянні 6), пов'язані з кореляційними часами різних порядків 012S . Таким чином, правила сум для кореляційних функцій будуть задовольнятися точно до порядку 2S у часовому просторі і до порядку 2S , включно, у частотному представленні. Кореляційні часи можна обчислити, застосовуючи рекурентні співвідношення при достатньо великих S.

Рис. 4. Основні поляризаційні моди, отримані для моделі води TIP4P.

Числові розрахунки проводились для моделі води TIP4P M4) при густині mNV1 г/см і температурі T293 K. Спектр поляризаційних мод, отриманий в рамках модифікованого підходу при S 1 і S3 для повздовжніх флуктуацій, подано на Рис. 4. Видно, що існують дві пропагаторні моди і дві чисто дифузійні моди . Пропагаторні моди описують лібраційні збудження в флуктуаціях поляризаційної густини , де і - частота коливань і коефіцієнт згасання, відповідно. Лібраційні збудження зумовлені швидким обертовим рухом молекул в електричному полі, створюваним сусідніми молекулами. Вони зникають при k3Е, коли . Релаксація міжмолекулярних моментів сил визначається модою дебаївського типу . Мода відноситься до дебаївської релаксації на великих часах. Величини і тісно пов'язані з коефіцієнтами трансляційної і обертової дифузії. При малих і проміжних k спостерігається така ієрархія частот:

.

Це свідчить про співіснування динамічних процесів, що відбуваються на суттєво різних масштабах часу.

Рис. 5. Часові кореляційні функції поляризаційної густини для моделі води TIP4P, отримані в рамках звичайного S 0) і немарківського S 1) підходів. Результати, що відповідають S3,5 і 7, зображено довго-, короткоштриховими і неперервними кривими, відповідно. Дані МД показано кружечками.

Нормовану часову кореляційну функцію k t f k t f k , а також узагальнену діелектричну проникність , що отримані в рамках звичайного S 0) і немарківського S 1) підходів для S3,5 і 7, зображено на Рис. 5 і 6, відповідно, разом із даними комп'ютерного експерименту див. розділ 4). Для функції k t можна спостерігати згасаючі дуже швидкі лібраційні коливання на коротких часах, які далі переходять у майже чисте спадання на проміжних інтервалах і нарешті в надзвичайно повільну дебаївську релаксацію на великих часах. Бачимо, що звичайна S 0) теорія УКМ нездатна задовільно описати часову залежність k tпри S3. Тут навіть високомодові наближення S5 і S7 не дають змоги суттєво покращити результат. З іншого боку, застосування модифікованого підходу вже з S 1 веде до якісного відтворення k t і при S3. Більше того, врахування немарківських ефектів з S 1 значно прискорює збіжність теорії, так що отримані значення для k t і при S5 і S7 стають майже нерозрізненними від даних МД. Подібні висновки було зроблено і для узагальненої зсувної в'язкості k , результати розрахунків якої подано на Рис. 7.

Рис. 6. Дійсна a) і c)) та уявна b) і d)) частини діелектричної проникності для моделі води TIP4P. Результати у марківському S 0 і S3) та немарківському S 1 при S3 і 5) наближеннях зображено штрих-пунктирними, штриховими і неперервними кривими, відповідно. Дані МД подано кружечками.

Рис. 7. Дійсна a) і c)) та уявна b) і d)) частини узагальненої зсувної в'язкості в одиницях Pa s) для моделі води TIP4P, отримана в рамках модифікованого S 3) підходу УКМ при S2 неперервні криві), у порівнянні з даними МД кружечки).

Четвертий розділ присвячено розвитку комп'ютерно-адаптованих теорій як для дипольних, так і атомних моделей полярних рідин. Оскільки комп'ютерний експеримент має справу з скінченними системами, то на великих відстанях взаємодія між частинками не може бути врахована явно. Для мінімізації скінченно-розмірних ефектів методом поля відгуку всі частинки поза сферою взаємодії замінюються неперервним середовищем з діелектричної постійної . У якості середовища можна вибрати провідник, де , бо полярні рідини характеризуються великим значенням . Тоді загальне макроскопічне поле, яке створюється частинками системи і середовищем, буде рівне з rrR, де позначає модифікований тензор дипольної взаємодії, , а - радіус сфери обрізання. Застосовуючи метод збурень і враховуючи наявність середовища, знайдено флуктуаційні формули для розрахунку узагальненої діелектричної проникності . Зокрема, для повздовжніх флуктуацій у канонічному ансамблі отримано

7)

де . Схожі формули одержано і у випадку мікроканонічного ансамблю та поперечних флуктуацій. Вони узагальнюють знайдені раніше [14] часткові розв'язки, коли k0 або 0, і можуть використовуватися для довільних значень k та .

Суттєвою відмінністю комп'ютерно-адаптованих формул 7) від відомих із звичайної теорії макроскопічних систем, коли R і , полягає у тому, що тут повністю беруться до уваги усі деталі, при яких виконуються симуляції, а саме, розмір зразка, граничні умови та спосіб обрізання потенціалу взаємодії. Для дипольних моделей модифікація тензора взаємодії приводить до зміни міжчастинкового потенціалу у вигляді

8)

де , тоді як при . Для моделей типу взаємодіючих силових центрів густину дипольного моменту p необхідно замінити на вектор поляризації P див. розділ 3). Як наслідок, кулонівська взаємодія між атомами a і b молекул i та j модифікується у геометрії поля відгуку для до

 9)

при , якщо . Рівняння 9) точно враховує просторовий розподіл атомів у молекулі для поля відгуку середовища. У цьому полягає основна відмінність від попередніх наближених підходів [15], де істинне поле середовища для атомних моделей підмінялось полем відгуку 8) системи точкових диполів.

Для атомних моделей розвинуто також метод підсумовування за Евальдом. Тут для нівелювання ефектів скінченного розміру вводяться періодичні копії симуляційного зразка, який має форму куба з ребром , і розглядаються взаємодії як між реальними атомами, так і їх зображеннями. Це можна інтерпретувати як ефективну взаємодію лише між молекулами базисної комірки, але з модифікованим кулонівським оператором

,

де при rR i для rR, а сума проводиться за векторами з цілими компонентами при . У результаті отримано флуктуаційну формулу типу 7), де необхідно замінити на з . При кожному фіксованому оптимальні значення параметра вибирались з умови мінімізації . Тоді буде якнайменше відрізнятися від часової кореляційної функції f k t нескінченної системи, а відносна похибка обчислень складатиме . Зокрема для і Rl/2 отримано 5.929/l з 0.132, що є цілком достатньо для більшості обчислень.

Запропонований підхід застосовано до кількох моделей води в рамках методів МД і МК. Досліджено самоузгодженість комп'ютерно-адаптованих формул 7) і показано, що система вже з декількох сотень молекул N200) може забезпечити об'ємні тобто ті, що відповідають безмежній системі) значення діелектричної проникності k на основі дипольних чи поляризаційних флуктуацій , що отримуються в комп'ютерних симуляціях для скінченного зразка. Вперше розраховано f k і в усій області зміни k і . Частину цих результатів вже було приведено на Рис. 2, 3, 5 і 6 у контексті їх порівняння з передбаченнями теорії УКМ розділи 2 і 3).

Рис. 8. Часові кореляційні функції моделі Штокмайєра для повздовжніх кружечки) та поперечних квадратики) компонент дипольного моменту скінченної системи. Результати для безмежної системи подано неперервними і штриховими кривими.

Рис. 9. Фактор Кірквуда води TIP4P, отриманий в рамках різних підходів.

Часові кореляційні функції скінченної системи, які обчислені в МД симуляціях при N 256 для моделі Штокмайєра, і функції k t f k t f k , що відповідають безмежній системі NR), подано на Рис. 8. Значення k t були відтворені, використовуючи співвідношення 7) і самоузгодженість теорії, тобто незалежність від розмірів системи при N 200. Видно, що різниця між k t і є суттєвою, особливо для малих k. Отже, ефекти скінченного розміру дійсно можуть значно спотворювати діелектричні величини і тільки застосування комп'ютерно-адаптованого підходу дає змогу отримати для них коректні значення.

Фактор Кірквуда скінченної N 256, R 9.9 Е) системи, одержаний в МД симуляціях для моделі TIP4P із застосуванням поля відгуку силових центрів, і значення gk) для безмежної системи зображено на Рис. 9 штриховою і суцільною кривими, відповідно. Результат для gk), отриманий в рамках підсумовування за Евальдом і поля відгуку точкових диполів, показано кружечками і точковою кривою. Як бачимо, поле відгуку взаємодіючих силових центрів веде практично до тих самих величин, що і громіздка техніка Евальда. З іншого боку, відхилення між цими результатами і значеннями, що відповідають полю відгуку точкових диполів, сягають порядку 20. Це засвідчує важливість врахування просторового розподілу заряду у молекулі при конструюванні поля відгуку для атомних моделей полярних рідин.

Побудова ефективних МД алгоритмів для знаходження числових розв'язків рівнянь руху в багаточастинкових системах становить зміст п'ятого розділу. Тут проведено узагальнення декомпозиційного підходу на атомні і молекулярні рідини за одночасної присутні як трансляційних, так і орієнтаційних ступенів вільності. Детально розглянуто випадок полярних рідин, для яких оператор Ліувілля можна записати як

 10)

де W - лінійна антисиметрична матриця, і задають положення і орієнтацію молекули i з масою m і тензором інерції J, що рухається з поступальною швидкістю і обертається з кутовою швидкістю при наявності кутового моменту кількості руху , а і позначають сили і моменти сил, які діють на молекули i12N) згідно атомних взаємодій. Рівняння руху частинок d/dtLt) має формальний розв'язок , де - набір фазових змінних, tt/l - величина часового кроку, а l - загальне число кроків.

Ідея декомпозиційного підходу полягає у розщепленні LAB оператора Ліувілля на його кінетичну A і потенціальну B частини. Тоді оператор еволюції можна представити як добуток експонент

11)

де коефіцієнти , і при фіксованому K2 підбирають у такий спосіб, щоб забезпечити найвищий можливий порядок Q точності, позначає комутатор двох операторів, а - це похибка обчислення. Основна перевага факторизації 11) полягає в тому, що, на відміну від повного оператора еволюції , дію його кінетичної і потенціальної складових на довільну точку фазового простору можна виконати аналітично:

12)

Тут і - модифіковані сили і їх моменти з градієнтними вкладами , де і при . Співвідношення 12) задають зміщення і -обертання координат вільний поступальний і обертовий рух) та швидкостей рівноприскорений рух) молекул, які не змінюють об'єму у фазовому просторі, а отже точно забезпечують такі важливі риси гамільтонових систем як їх симплектичність і самоспряженість, незважаючи на наближений характер факторизації 11).

Таблиця 1: Декомпозиційні алгоритми

Було проведено повну класифікацію і детальне виведення всіх декомпозиційних алгоритмів аж до шостого порядку включно. Сімейство алгоритмів другого і четвертого порядків подано в Табл. 1, де літери A та B ) або C ) позначають відповідні експоненти, які входять в 11). Серед них тільки два алгоритми, а саме, Верле №1) і Фореста-Рутха №11), були відомі раніше [5]. Бачимо, що нові схеми під №5, 8, 14 і 15 є безумовно кращі, оскільки забезпечують найвищу ефективність при заданому числі обчислень сил на один крок. Зокрема, алгоритм №5 є найкращим у випадку другого порядку точності і перевищує ефективність схеми Верле у три рази. Якщо ж порівнювати схеми четвертого порядку, то видно, що запропонований алгоритм №8 дає змогу при тих самих затратах ) комп'ютерного часу покращити точність інтегрування рівнянь руху майже на два порядки по відношенню до схеми Фореста-Рутха.

У цьому розділі запропоновано також модифікацію методу декомпозицій, що базується на канонічному перетворенні фазового простору , де і при , а і . Записуючи рівняння руху у нових фазових змінних і застосовуючи факторизацію 11) при , отримано узагальнений алгоритм Верле

13)

Тут дія експоненційних операторів виконується у трансформованому просторі F і визначається рівняннями типу 12) з формальною заміною на F, а перехід до істинних змінних здійснюється оберненим перетворенням . Часові похідні, які виникають в і , обчислювалися за допомогою скінченнорізницевої схеми

 .

Це не вимагає додаткових затрат, оскільки використовуються ті ж сили, що і при пропагуванні Ft). Параметри вибирались так, щоб порядок точності зріс з Q2 до 4. У результаті знайдено і , з . Другий варіант є особливо ефективний, оскільки долає бар'єр з до для алгоритмів четвертого порядку, зменшуючи, таким чином, затрати комп'ютерного часу втричі.

Рис. 10. Затрати часу для різних алгоритмів в залежності від їх похибки.

Залежність обчислювальних затрат від похибки R симуляцій, що проводилися різними алгоритмами для моделі води TIP4P, зображено на Рис. 10, де кружечки відповідають часовим крокам зліва направо) t1234 і 5 fs. Бачимо, що навіть для найменших значень 200, модифікована схема 13) з PS1) забезпечує прийнятну точність 1, а алгоритм Верле VT) веде при цьому до відносно великої похибки R7. При збільшенні переваги модифікованого підходу як при , так і з PS2), стають ще більш очевидними. Наприклад, алгоритми №5 VO) і №8 GL), які є кращі за стандартні схеми VT і Фореста-Рутха FR), відчутно поступаються інтегратору PS1. Зокрема, при 500 для схем FR, VT, VO, PS2 і GL маємо R4, 1, 0.4 0.2 і 0.1%, відповідно, тоді як алгоритм PS1 може покращити точність обчислень більш ніж на два порядки і звести похибку майже до нульового рівня R0.02 без жодних додаткових затрат.

Детальне дослідження фазових властивостей магнітних рідин із спіновою взаємодією проведено у шостому розділі. Для моделей Ізинга і1), XY і2) та Гайзенберга і3) запропоновано модифікацію теорії середнього поля, яка враховує м'якість потенціалу відштовхування r) між частинками на малих відстанях. Феромагнітна взаємодія між спінами задається потенціалом Юкави , де Jr)exp[r)/]r, а s - це вектор спіну вимірності і. Враховуючи, що у наближенні середнього поля радіальна функція розподілу , знайдено рівняння стану

,

де позначає тиск системи відліку, - повне магнітне поле, яке включає зовнішнє B, а - намагніченість системи, що залежить від h. З нерівності 0 при B=0 отримано температуру фазового переходу пара-феромагнетик, а з умов механічної і хімічної рівноваг визначено критичну температуру фазового переходу газ-рідина та густини співіснування при .

Залежності від B, отримані в рамках звичайної HSMF) і модифікованої SСMF) версій середнього поля, зображено на Рис. 11 у порівнянні з даними МК кружечки). Як бачимо, обидва підходи якісно відтворюють функції , хоча відхилення від значень МК у випадку теорії SСMF є меншими порядку 25%). Розглянуто також і версію ASСMF), коли при обчисленні діаметр частинок модифікується так, щоб результати теорії збігалися з даними МК у границі B, де значення стає незалежним від і. Версія ASСMF приводить до майже точних результатів для в усій області зміни B для всіх трьох моделей. Видно, що характер залежності суттєвим чином визначається типом спінових взаємодій. Функція монотонно спадає з ростом поля B для моделі Ізинга, тоді як для взаємодій XY і Гайзенберга маємо немонотонну поведінку. Така складна залежність є проявом кореляцій між трансляційними та спіновими ступенями вільності у магнітних рідинах. Розраховувались також фазові діаграми , де, однак, спостерігалось лише якісне узгодження з результатами МК навіть у рамках підходу ASСMF.

Рис. 11. Залежність критичної температури від зовнішнього поля і ), отримана для зліва направо) моделей Ізинга, XY та Гайзенберга на основі різних версій теорії середнього поля, у порівнянні з результатами МК.

Щоб прийти до кількісного опису фазових діаграм, розвинуто метод інтегральних рівнянь для моделі Ізинга і вперше запропоновано його узагальнення для магнітної XY рідини. Таке узагальнення базується на анізотропному рівнянні типу Орнштейна-Церніке

14)

для прямої c і повної h двочастинкових кореляційних функцій, а також на першому рівнянні з ланцюжка Боголюбова-Борна-Гріна-Кірквуда-Івона

15)

для одночастинкової функції ), де кут задає орієнтацію спіна щодо поля. Для замикання використано середньо-сферичне наближення expuhcb)h1 з bln1) і hc[Ir)Jr), де - міжчастинковий потенціал взаємодії, який включає немагнітне притягання Ir)Jr)/Y. Кореляційні функції було розкладено за гармоніками і як і . Тоді рівняння 14) і 15) спрощуються і в k-просторі набирають алгебраїчного вигляду відносно коефіцієнтів розкладу

16)

де , а з . Розв'язуючи рівняння 16) і використовуючи віріальний тиск густини співіснування газ-рідина визначались з конструкції Максвелла

, де .

Трьох гармонік 0nm3) було достатньо для того, щоб забезпечити збіжність результатів.

Рис. 12. Фазові діаграми газ-рідина ) для моделі XY при декількох значеннях Y і . Перехід пара-феромагнетик при ) показано точковою кривою.