Нелинейные моды и транспорт энергии в полимерных цепях Текст научной статьи по специальности «Физика»

===============^=^===^= ТЕОРИЯ

УДК 541.64:539.2

НЕЛИНЕЙНЫЕ МОДЫ И ТРАНСПОРТ ЭНЕРГИИ В ПОЛИМЕРНЫХ ЦЕПЯХ1 © 2005 г. Л. И. Маневич, А. В. Савин

Институт химической физики им. H.H. Семенова Российской академии наук 119991 Москва, ул. Косыгина, 4 Поступила в редакцию 09.04.2004 г. Принята в печать 15.11.2004 г.

Представлены результаты аналитического и численного исследования динамики нелинейных возбуждений в молекулярных цепях. Показано, что существование и распространение случайных возбуждений в рассмотренных моделях существенно зависит от свойств нелинейных нормальных мод. Так, нелинейные колебания в зоне поглощения линеаризованной системы, инициированные термостатом, соответствуют локализованным нелинейным нормальным модам - бризерам. В области локализации таких мод происходит периодическое сжатие (растяжение) межатомных связей, которое сопровождается уменьшением (увеличением) валентных углов. Установлено, что бризеры существуют в термализованной цепи, и рассчитана температурная зависимость их вклада в теплоемкость. Изучен процесс переноса тепла в одномерной решетке связанных ротаторов, в которой ориентаци-онное взаимодействие соседних звеньев описывается периодическим потенциалом. Показано, что конечная теплопроводность данной решетки обусловлена существованием локализованных стационарных нелинейных возбуждений (ротобризеров), препятствующих движению фононов. На примерах моделей Френкеля-Конторовой, ф-4 и sinh-Гордон продемонстрировано, что характер теплопроводности гармонической одномерной решетки с непараболическим потенциалом подложки определяется спектром нелинейных возбуждений в данной цепи. Нелинейные возбуждения определяют конкретный механизм рассеяния энергии, который обеспечивает сходимость теплопроводности. Для моделей sin-Гордон и ф-4 (при низких температурах) - это рассеяние фононов на термализованной решетке топологических солитонов, а для sinh-Гордон и ф-4 (при высоких температурах) - рассеяние на локализованных высокочастотных дискретных бризерах (для ф-4 механизм рассеяния энергии изменяется с ростом температуры).

ВВЕДЕНИЕ

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

1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (коды проектов 04-02-17306 и 04-03-32119).

E-mail: [email protected] (Маневич Леонид Исаакович); [email protected] (Савин Александр Васильевич).

ципиально новые особенности. Указанные особенности могут быть поняты, если процесс распространения рассматривается в терминах неравновесной статистической механики. С такой альтернативной точки зрения случайные колебания в молекулярной цепи соответствуют глобальному, а локализованные случайные возбуждения - локальному тепловому воздействию. Распространение случайных возмущений (перенос механической энергии) вдоль цепи ассоциируется с процессом теплопроводности.

Исторически эта проблема восходит к знаменитому численному эксперименту Ферми-Паста-Улама (1955 г.) [1], который был предпринят для обоснования модели теплопроводности в диэлектриках, предложенной Ре1ег1з (1929 г.) [2]. Эта модель связывает теплопроводность с так называемыми процессами переброса, в которых слабые

х

Рис. 1. Схематическое представление моноклинной структуры кристаллического ПЭ. Показан центральный транс-зигзаг (0) и шесть ближайших к нему цепей (1-6). Для центральной цепи даны локальные системы координат.

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

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

Численные эксперименты показывают, что в нелинейной осцилляторной цепи возможны два

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

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

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

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

ЛОКАЛИЗОВАННЫЕ НЕЛИНЕЙНЫЕ КОЛЕБАНИЯ (БРИЗЕРЫ) В ЗИГЗАГООБРАЗНОЙ ПОЛИМЕРНОЙ ЦЕПИ

Как известно, кристалл ПЭ состоит из параллельно упакованных зигзагообразных цепей молекулы ПЭ (СН^. Схематически структура кристалла представлена на рис. 1. Каждая отдельная макромолекула ПЭ в кристалле находится в кон-формации транс-зигзаг, т.е. ее скелет имеет плоскую зигзагообразную структуру, которая характеризуется шагом р0 = 1.53 А (равновесной длиной валентной связи Н2С-СН2) и углом зигзага 60= 113° (равновесным валентным углом СН2-СН2-СН2). Все макромолекулы в моноклинной решетке находятся в параллельных плоскостях.

Рассмотрим отдельную молекулу ПЭ. Так как при исследовании низкоэнергетических нелинейных динамических процессов движение атомов водорода относительно основной цепи не являет-

ся существенным, в дальнейшем воспользуемся а косинус п-го торсионного угла приближением объединенных атомов, т.е. будем

считать каждую группу СН2 одной частицей мае- cos(5„) =

сыМ= 14тр (гар - масса протона). = (~bnAbn+hl - bn<2bn + h2 +Ьп,гЬп + 1Ж^ + и

(4)

Рассмотрим динамику только центральной це- где пи (к = 0). Выберем систему координат так, чтобы в положении равновесия узлы транс-зигзага имели координаты

<и , 1чп, ,г, и) „ (и) ,

= (-1) У2, Уп =0, г„ = п1г,

где п = 0, ±1, ±2, ... - номер узла цепи, 1Х = = росо8(0о/2) - поперечный шаг зигзага, /г = с/2 = = р0 зш(0о/2) - продольный шаг зигзага. Перейдем от абсолютных координат п-го звена транс-зигзага хп, уп, 1п к относительным координатам

, , . л + 1 , (Ок , - чп+ 1

ип = (-1) (х„-хп ), = (-1) >>„,

W =7

ry п

.(0)

п = 0, ±1, ±2, ...

Схематически локальные системы координат даны на рис. 1.

ЬпА = ап-1,2ап,Ъ + ап,2ап-1,г

Ьп,2 ~ ап- 1,1ап, 3 + ап, \ап-\, 3

Ьп, 3 ~ ап-1,2ап, 1 ~ ап,2ап-1, 1 О п.1 IЛ 1_ 2 .1/2

Следуя работам [3,4], примем потенциалы валентной связи, валентного угла и торсионного угла в виде

У(р„) = £>0{1 -ехр[-а(р-р0)]}2

U( 6Я) = |y(cos0 - cos0o)2

W(5„) = Cj + C2cos6n + C3cos35„

Функцию Гамильтона цепи можно записать Здесь D0 = 334.72 кДж/моль, а = 1.91 Á"1, у =

= 130.122 кДж/моль, С, = 8.37 кДж/моль, С2 = = 1.675 кДж/моль, С3 = 6.695 кДж/моль.

следующим образом:

Н = Щ±М(й2п+<г2п + м,2п) +

п ^

+ v(p„) + [/(e„) + w(5n)}

(i)

Здесь первое слагаемое характеризует кинетиче-

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

скую энергию п-го узла, второе - энергию дефор- возбуждений в квазиодномерных цепях осцилля-мации п-й валентной связи, третье - энергию де- тоРов' основанный на комплексной формули-формации п-го валентного угла, а последнее - Ровке Уравнений движения с последующим ис-

пользованием многомасштабных асимптотических разложений. Длинноволновые солитонные возбуждения в кристаллах ПЭ и ГГГФЭ изучены в работах [6-9]. Влияние геометрии цепи на динамику нелинейных возбуждений детально анализируется в работе [10]. Коротковолновые солито-

энергию деформации п-го торсионного угла. Длина п-й валентной связи

, 2 2 2 ч1/2 Р п = («п.1 +«п,2 + «п,з)

(2)

, , , , , ,, ны (бризеры) в одномерных цепях обсуждаются в

(an,l = un + un+l-lx,ant2 = vn + vn+l,an<3 = wni.l-wn + lz), 4 ^ \ ' * _

и т ,„,„„ статьях [11-18], а их моделирование в кристалле

ПЭ - в работе [19]. Обзорная статья [20] содер-

косинус п-го валентного угла

cos(6„) =

- (ап- 1,1ап, 1 + ап-1,2ап,2~ап-1,Зап,з)^Рп-1рп' ВЫСОКОМОЛЕКУЛЯРНЫЕ СОЕДИНЕНИЯ Серия А том 47 № 5 2005

жит анализ результатов теоретического, числен-(3) ного и экспериментального исследования физических закономерностей, обусловленных наличи-

л

к/2

0.01

0.01

о 0.02

0.01

О

/ 'Ч (а) 1 1

; (б) 1_)

(в)

| и (Г) ь , J -

О

(д)

400

800

1200 ш, см

г-1

Рис. 2. Дисперсионные кривые (О = 0)а(#) (7), (о = (2), <в = «<,(<?) (3) для транс-зигзага ПЭ с подложкой, образованной шестью неподвижными соседними цепями (а), а также плотность распределения тепловых колебаний р по частотному спектру ш при температуре Т = 1 (б), 100 (в), 200 (г) и 300 К (д). Штриховкой отмечена область частот, соответствующая дискретным бризерам.

ем локализованных нелинейных возбуждений в полимерных кристаллах.

Малоамплитудные колебания изолированного транс-зигзага были рассмотрены в работах [6, 21,22], а с учетом взаимодействия с неподвижной подложкой - в работе [7]. В этом случае колебания можно разделить на плоские (колебания в плоскости зигзага) и поперечные (с выходом из плоскости). Плоские колебания в свою очередь разделяются на низкочастотные акустические и высокочастотные оптические колебания. Соответственно дисперсионная кривая зигзага имеет

три ветви: ветвь плоских акустических фононов (О = (йа(д), ветвь плоских оптических фононов со = = «><,(<?) и ветвь поперечных колебаний ю = сог(д). Вид дисперсионных кривых показан на рис. 2а.

Эти кривые получены на основе решения линейных уравнений, соответствующих функции Гамильтона

Н = \%{Щй1+*1 + *2Л) + К1(рЛ-р0)г +

(5)

+ к2(вп - 0О)2 + ад<р„+1 - Ф») - (ф„-1 - Фл-г)]2}

при линеаризованных относительно смещений деформациях [6]

Ар « = Р«-Ро = = (н'„+1-н'и)8т(0о/2)-(мп + 1 + мл)со8(0о/2) Д6„ = 0„-0о = [(>Ул+1-Н'п_1)СО8(0о/2) + + (ип + у + 2ип + ип _ I) вт (0о/2) ]/р0

(6)

Рассмотрим вначале динамику в плоскости транс-зигзага (плоская модель макромолекулы с "объединенными" атомами представлена на рис. 3) при частотах, лежащих в окрестности нижней границы первой оптической ветви ИК-спектра

К,

М

соь {%!!) +4 К

2 8Ш (0о/2)]

(7)

(8)

н'у±1 = и/(^х)±е8т(9о/2)^ +

1 2 . 2,л ,лчЭ2м;

+ -£ вт (0о/2)—+ ...,

2 ос,

п-2

Из анализа линеаризованных уравнений движения следует, что при частотах, близких к Юц трансверсальные и продольные смещения существенно отличаются по порядку величины, так что ~ еи„, где е = к17 (е 1), к - волновое число. Таким образом, введенный малый параметр 8 характеризует отношение характерной длины волны в рассматриваемой области спектра к расстоянию между соседними группами СН2. Малость этого отношения позволяет перейти к континуальному приближению при моделировании как линейной, так и слабо нелинейной динамики.

Введем безразмерные пространственную, а также временную координаты

£ = ег/ро, т = «V,

где г = Д = 7ро81п(0о/2) в точках, в которых расположены объединенные атомы, и степенные разложения

Э и

ит = т)±еяп(0о/2)^| +

1 2 . 2.л .„.Э2и + -е яш (0о/2)—+...

1 6с,

Рис. 3. Структура плоской макромолекулы полиэтилена (СН^ (а), выбор системы локальных координат (б) и плоская механическая модель системы (в).

где безразмерные смещения групп СН2, интерполированные в одномерном континууме; х), х) определяются соотношениями

Иу = £ро"(^ X), = £р0и>(§, X)

(9)

При этом деформации транс-зигзага в его плоскости принимают вид

Ар = ер0Др„ = = £ро] -СО8(0о/2)

Эм

2 и + £8Ш(0о/2)ггё +

1 2 . 2,л + -£ вш (0о/2)—+ ...

2 Эн* 2 2 + £вт (0о/2)^ + 2£81п (0о/2)и + ...

(10)

А0 = еА0„ = £

г)

48ш(0о/2 )и + г2$т(в0/2)—2 +

ЭГ

2£8Ш0О^ + 4£8Ш0ОМ2 + ...

где не выписанные члены имеют более высокий (начиная с третьего) порядок малости по параметру £.

В силу соотношений (9) и (10) параметр £ определяет также относительную малость интенсивности возбуждения.

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

Подставляя выражение (14) в первое из уравнений (12), получим

д2и . 2 од 2 3 2 2Э2М „ /1СЧ

+ и - 4а ем + 8ре и -г с —5 + ... =0 (15)

Эх

ЭГ

(с2 = с 1 + Т2).

дт) 1эх;

Мсо,

^др' + ^дё2)

(П)

Согласно работам [5,20], введем комплексные функции

х) = у + ш, х) = у-ш, (16)

После подстановки выражений (10) в формулу где ^ = Эи/Эт. Уравнение в частных производных

(11) и варьирования функции Лагранжа, соответ- второго порядка (15) можно записать как систему

ствующей гамильтониану (11), приходим к двум дВух уравнений первого порядка по времени связанным нелинейным уравнениям с частными производными относительно функций м(£,, т), X):

2 2 д и . 2оГ>23 2 2Э и Эм>

—- + и -Лаеи + 8ре и -е с1

Эх

Э2И>

ди ,.

= 0 (12)

Здесь

а - 3 2 Мсо2

8т608ш(60/2)^1 -

8К2 К

М со, V

_ ^зте^т^во^)/ _ АК.

1 - о *

Мсо2 V

^хРо

Эу 2 гЭ'ы . 2 о0 2 3 -г— = е с —г-м + 4аем -8ре и Эх э^

Эм Эх

(17)

= V

После подстановки функции (16) в систему (17) и ряда преобразований получим уравнение

Эм/ . 2 гЭ2(х1/ — \1/*)

Эх г Э£2

+ ае(у - у* )2 - г р£2(у - у* )3 = 0

Замена переменных у = еп(р(С1, х) позволяет привести это уравнение к виду

Эф . 2 гЭ (<р-<р*е <г) , ¿х л -¿хч2 —|х т-2- + :е с———-- + аг(<ре ~(р*е )е -

Эх Э?2 (18)

-фе (<реп-(р*е ") е п = 0

с, = —

К^ш(60/2)

А/со,

? 8 К-у ? сое (60/2) +-^¡п (в0/2)

Введем наряду с "быстрым временем" х = х0 "медленные времена"

С другой стороны, в принятых безразмерных переменных

Х[ — ех0, х2 — £ х0, и степенное разложение

(19)

Э2н>

= - IV ■

Эх

Следовательно,

Эм

н' = +еу^г + ...

(13)

ф = ф0 + £ф1 + £ ф2+ •••

(20)

После подстановки разложения (20) в формулу (18) с учетом выражения (19), выделения членов одного и того же порядка относительно мало-го параметра е и приравнивания их нулю прихо-

дим к выражениям ВЫСОКОМОЛЕКУЛЯРНЫЕ СОЕДИНЕНИЯ Серия А том 47 № 5 2005

Эфо

Эт0

= О

-3i'T„

Эф! Эф0 2 "о i i ""О - --"'о л ^ + ^ + -2а|ф0|е + а(ф£)<? =0

2

Эф2 ^ Эф, Эфо 2э Фо /х0

-31Тп

-2афоФГ^ -2аф0ф*е -2афоф,е +

•о г 3 2,t0 + ф[ф0« '

3|фо|2фо + 3|фо|2ф^

-2¿t„

+ (Фо*)^"4,Ч°] = О

Интегрирование первого уравнения и условие отсутствия секулярных членов в остальных уравнениях этой системы приводят к соотношениям Фо = ФоСЧ, х2' ■••). Эфо/Эт! = 0, следовательно Фо = Фо(х2» тз> •••)• Тогда решение второго уравнения может быть записано следующим образом:

2 <*о , . , ,2 ->Ч I, * .2 -3(Х0

ф! = Шф0е +2т|ф0| е -^(Фо) «

Условие отсутствия секулярных членов в третьем уравнении может быть представлено в виде

Эфо + ic2<fjPo + /16 Эт, эг;2 ^3

а2-Зр)|ФоГ

Фо = 0 (21)

Итак, приходим к известному точно интегрируемому нелинейному уравнению Шредингера. Если

(22)

то это уравнение имеет локализованное солито-ноподобное решение - бризер:

Фо(С. ъ) = (у) ехр('27 + ' х sechj^51/2^- + vx2j

где

со

= (тЧ

(24)

Возбуждения такого типа характерны и для более простых квазиодномерных цепей [11-17]. Чтобы проверить допущения, принятые при выводе уравнения (21), а также устойчивость бризе-ров при столкновении с локализованными возбуждениями и при термических возмущениях, было выполнено численное исследование.

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ

ЛОКАЛИЗОВАННЫХ НЕЛИНЕЙНЫХ КОЛЕБАНИЙ

Гамильтониану цепи в трехмерном случае соответствует система уравнений движения

э н ,... э н .... э н

Ми„ = -т—, Мч. = —, Мм>„ = -тг—,

Эм„ " дн>п (25)

п = 0, ±1, ±2, ..., где гамильтониан Я задан формулой (1).

Для численного моделирования динамики рассмотрим конечную цепь из N = 200 звеньев. На концах цепи введем вязкое трение, обеспечивающее поглощение излучаемых фононов. Система уравнений движения (25) с п= 1,2, ..., N и бризе-роподобным начальным условием интегрировалась численно. Численное интегрирование показало, что после излучения фононов, связанного с наличием небризерной части в начальном условии, в цепи формируется дискретный бризер (устойчивое нелинейное локализованное периодическое колебание). В области локализации этого нелинейного возбуждения наблюдается периодическое изменение валентных связей, сопряженное с изменением валентных углов - рис. 46, 4в. Колебание происходит в плоскости зигзага со смещением узлов цепи в поперечном направлении - рис. 4а.

Такое локализованное колебание является устойчивым нелинейным возбуждением, которое характеризуется частотой со, энергией Е и безразмерной шириной

(23)

Г N

L = 2

Ln= 1

1/2

где точка пс = ^ппр„ задает положение центра

возбуждения, а последовательность рп = EJE -распределение энергии вдоль цепи. Нелиней-

п

Рис. 4. Локализованное плоское колебание транс-зигзага цепи ПЭ. а - колебания показаны схематически, толщина линии соответствует амплитуде колебания; б, в - изменения длины валентных связей р„ (б)

и валентных углов 9„ (в) для десяти разных значений времени. Частота бризера со = 820.5 см-1, энергия Е = 26.4 кДж/моль, ширина £ = 4.28.

ность проявляется в зависимости частоты от амплитуды колебания (с ростом амплитуды частота монотонно снижается).

Зависимость энергии Е и ширины бризера Ь от его частоты со представлена на рис. 5. Частотный спектр бризера расположен левее нижней границы спектра оптических фононов.

ТЕПЛОВЫЕ КОЛЕБАНИЯ трянс-ЗИГЗАГА КАК ИСТОЧНИК ДИСКРЕТНЫХ БРИЗЕРОВ

Рассмотрим тепловые колебания транс-зигзага. Для этого возьмем конечную цепь из N звеньев. Концевые /У0 звеньев с каждого конца цепи поместим в термостат с температурой Т. Динамика системы задается системой уравнений Ланжевена

Miin = - — + Г пМйп

Mvn = -b£ + 4n-TnMun (2б)

Mwn = + \Мйп, п = 1,2, ...,N,

где гамильтониан системы определяется формулой (1), случайные силы rj„, С,п и коэффициенты трения Г„ равны нулю для N0<n<N- N0. При n<N0aN-N0<n<N коэффициенты Гл = Г, а t,n, Т|„, - случайные нормально распределенные силы, описывающие взаимодействие п-й молекулы с термостатом. Коэффициент трения Г = 1 /tr, где tr - время релаксации скорости молекулы.

Случайные силы имеют корреляционные функции

<$»(il)$m(f2)> = <Л„('1)Пт('2)> = (Uh)Uh)) =

= 2мгквтьптНь-ь)

(UtMh)) = = <4n(tl)Uh))=o

l<n, m<N0, N-N0<n, m<N

Здесь кв - постоянная Больцмана, Т- температура термостата (для внутреннего участка цепи N0<n<N-N0 случайные силы отсутствуют).

Систему уравнений (26) интегрировали численно стандартным методом Рунге-Кутта четвертого порядка точности с постоянным шагом интегрирования At. При численной реализации дель-та-функция 5(г) = 0 при |;| > At/2 и 5(f) = 1/Дt при 111 < At/2, т.е. шаг численного интегрирования соответствует времени корреляции случайной силы. Поэтому для корректного использования уравнения Ланжевена необходимо, чтобы выполнялось условие At < tr. Примем tr = 0.1 пс, a At = 0.001 пс.

Рассмотрим распределение кинетической энергии тепловых колебаний молекул цепи по частотам. Для этого проинтегрируем систему уравнений (26) при N = 500, N0 = N/2 с начальными условиями, соответствующими основному состоянию цепи, в течение времени t = 10ir. Термализацию цепи удобно контролировать, следя за изменением среднего значения кинетической энергии системы

N

М , .2 .2 . 2Ч

я*(0 = 2лД(м"+ vn + wJ

я = I

В начальный момент времени энергия Ек(0) = 0. Затем она начинает увеличиваться, что сопровождается флуктуационными изменениями, амплитуда которых пропорциональна обратной длине цепи 1/N. Рост практически прекращается при t = 5tr. Кинетическая энергия достигает нужного среднего значения Ък^Т/2 и в дальнейшем лишь флуктуирует около данного значения. Таким образом, в момент времени t = 10tr координаты и

скорости системы {и„, йп, v„, vn, wn, wn , будут заведомо задавать одно из множества возможных состояний термализованной цепи. В

о), см-1

Рис. 5. Зависимость энергии Е (а) и ширины L (б)

дискретного бризера от его частоты со.

дальнейшем связь с термостатом можно отключить (для этого достаточно в уравнениях движения (26) положить Г„, t,n, Т|я, = 0) и рассмотреть динамику изолированной цепи, используя полученное термализованное состояние в качестве начального. Найденные динамические параметры цепи будут зависеть от используемого термализо-ванного состояния цепи. Чтобы найти их точные значения, нужно провести усреднение по всем возможным термализованным состояниям цепи. Для получения нового, независимого от предыдущего, состояния термализованной цепи достаточно повторить в течение времени t = 10/г описанное выше численное интегрирование системы (26), используя при этом уже другую независимую реализацию случайных сил §„(i), Пп(/)> Так, при вычислении плотности распределения кинетической энергии молекул цепи по частотам р(со) было взято 1000 независимых реализаций термализа-ции цепи. Вид плотности распределения при разных значениях температуры дан на рис. 2 (плотность распределения нормирована условием

jV(co)dco = 3).

Рис. 6. Выделение дискретных бризеров из тепловых колебаний зигзага (N ~ 500, Т = 200 К). Условие поглощающих концов (N0 = 50). Дана зависимость распределения энергии по цепи Еп от времени г.

т„ к \

200ч--"" » ""

Рис. 7. Разрушение бризера в термализованной цепи (N=200, = 10, Т= 10 К, частота бризера

ю = 820.5 см-1). Дана зависимость от времени < текущих локальных значений температуры (кинетической энергии звеньев цепи) Тп.

При температуре Т = 1 К плотность распределения практически совпадает с плотностью распределения линеаризованной системы. Здесь колебания остаются линейными, термализуются лишь фононы. С повышением температуры растет амплитуда тепловых колебаний и начинает проявляться их ангармонизм. При Т = 100 К уже заметен сдвиг плотности распределения за нижнюю границу спектра высокочастотных оптических фононов, который увеличивается с ростом температуры. Появляются высокочастотные колебания, не соответствующие оптическим фоно-нам. По частотам колебания попадают в частотный интервал дискретного бризера [соь, со0(0)], поэтому их можно идентифицировать как бризеры. Доля энергии бризеров определяется как интеграл

<йо(0)

Рь = | р(ю)Ж»

<°ь

Вклад бризеров в тепловые колебания цепи растет с повышением температуры (при Т = 1 К рь = 0.002, при Т= 100 Крь = 0.106), достигает максимального значения рь = 0.115 при Т = 200 К, а затем начинает уменьшаться (при Т = 300 К вклад бризеров р = 0.083).

Выделим бризеры из тепловых колебаний. Для этого возьмем цепь из N = 500 звеньев с кон-

цевыми звеньями (Л/0 = 50), погруженными в термостат с температурой Т. После полной термали-зации положим температуру термостата Т = 0 и рассмотрим процесс стока тепловой энергии из внутреннего участка цепи (А^ < п < N - А^). Релаксация энергии при Т = 200 К показана на рис. 6. Хорошо видно образование нескольких подвижных локализованных возбуждений. Их анализ показал, что это дискретные бризеры с частотами со ~ Таким образом, дискретные бризеры присутствуют в тепловых колебаниях.

Для изучения взаимодействия дискретного бризера с тепловыми фононами в центре конечной цепи (№ = 200) возбуждался стационарный бризер с частотой со = 820.5 см-1, а концы цепи помещались в термостат температуры Т = 10 К. Динамика бризера показана на рис. 7. Видно, что бризер проявляет тенденцию к разрушению, как только центр цепи начинает термализоваться. За 20 пс бризер теряет примерно половину своей энергии.

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

ной - рь растет при Т > 200 К и уменьшается при Т> 200 К. Наибольшее значение достигается при Т= 200 К. Численное моделирование показывает, что при этой температуре бризеры наиболее отчетливо выделяются из тепловых колебаний.

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

ОДНОМЕРНАЯ РЕШЕТКА

С ПЕРИОДИЧЕСКИМ ПОТЕНЦИАЛОМ МЕЖДУУЗЕЛЬНОГО ВЗАИМОДЕЙСТВИЯ

Примером одномерной решетки одинаковых ротаторов может служить линейная макромолекула, в которой кроме продольных деформаций возможны еще крутильные (внутренние вращения вокруг жестких валентных связей). Здесь потенциал междуузельного взаимодействия оказывается периодической функцией (поворот одного звена на 360° переводит его в исходное состояние).

В качестве модели такой решетки целесообразно рассмотреть цепь атомов, расположенных друг от друга на фиксированном расстоянии /. Следуя работам [23, 24], предположим, что молекулы могут совершать только повороты вокруг оси цепи (рис. 8). Пусть переменная ф„(г) задает поворот п-го мономера цепи в неподвижной системе координат. Тогда функция Гамильтона системы в безразмерной форме будет иметь вид

Я = Х|^ + г/(фп + 1-ф„)1, (27)

п '

где точка обозначает дифференцирование по времени I, а потенциал междуузельного взаимодействия (потенциал вращения) £/(ф) - неотрицательная периодическая функция с периодом 2к,

Рис. 8. Механическая модель одномерной решетки ротаторов с периодическим потенциалом междуузельного взаимодействия.

удовлетворяющая условиям Щ0) = 0, и\0) = 0, £/"(0) = 1. Для определенности примем потенциал взаимодействия в простейшем виде

£/(ф) = 1 - совф (28)

При численном моделировании теплопроводности цепи обычно используется термостат Нозе-Хувера [25,26]. Однако в работе [27] показано, что его использование не обеспечивает корректную термализацию цепи. Поэтому воспользуемся классическим термостатом Ланжевена.

Рассмотрим конечную цепь из N звеньев с концами, погруженными в термостат Ланжевена, с безразмерными температурами Г+ и Т_. Соответствующая система уравнений движения цепи имеет вид

п = 1, ...,ЛГ0 фп = ^(ф„ + 1-фл)-^(ф„-фп_1),

п = + 1, ...,

фп = ^(Фп+1-Фл)-^(Ф„-ф„-1)-уфП + П„,

п = М-Ы0+1,

где функция Дф) = ¿/[/(ф)М|), /У0 - длина погруженных в термостат концевых участков цепи, у- коэффициент релаксации, Т|„ - моделирующий взаимодействие с термостатом белый гауссов-ский шум «£„(*)> = о, <%(/)> = о, Шн)цп({2)) = 0, (ЫПШЬ)) = 2уТ+ЬпДч - Н), <Г|„(?1)Л*('2)> =

= 2уГ_5„,5(*2-г1)).

Систему уравнений движения (29) интегрировали численно. После наступления теплового

0.004

0.002-

Рис. 9. Распределение локального теплового потока Зп (а) и локальной температуры Тп (б) в цепи с периодическим потенциалом взаимодействия. N = 500, ЛГ0 = 50, Т+ = 0.11, Т_ = 0.09, время усреднения ; = 5 х 106.

Данный способ термализации цепи позволяет решить проблему граничных условий. Распределение в цепи теплового потока У„ и температурного профиля Тп показывает (рис. 9), что на внутреннем участке цепи < п < N - Ы0 наблюдается свободный тепловой поток. Температурный градиент имеет линейный профиль, а локальный тепловой поток не зависит от номера звена цепи: ./„ = У. Это дает возможность определить коэффициент теплопроводности, используя только внутренний участок цепи:

к (W,) = JNi/(Tn - Тк_и),

(32)

где А^ = N - 2N0 - длина участка. Предельное значение

к = lim к(АМ

(33)

соответствует коэффициенту теплопроводности цепи при температуре Т = (Т+ + Т_)/2.

Коэффициент теплопроводности можно также найти по формуле Грина-Кубо [28]

i

кс = lim lim fC(x)dx (34)

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

Т„ = <<j>Ü(i)>« =

и локальный тепловой поток

(30)

Jn= UM, = limU jn(x)dx, (31)

/ —, GO Г J

гдejn = -[F(<t>„ +1 - ф„) + ЯФ„ - ф„_ i)]фл /2. Принимали значения у = 0.1, N0 = 50, N = 150, 200, 300,500, 700,900, 1300, 1700, 2100, 2400 и начальные условия, соответствующие основному состоянию цепи {ф„ = 0, ф„ = 0 i. За время t = 105 решетка приходила в тепловое равновесие с термостатом. Средние значения (30), (31) находили из последующей динамики системы в течение времени t = = 106-107.

Здесь функция корреляции С(?) = (У/хЩх - 0)т, а = ^„./л(0 _ общий поток тепла в цепи.

Для нахождения функции корреляции с(1) рассмотрим конечную циклическую цепь изЛ/ = 4000 звеньев, полностью погруженную в термостат Ланжевена. После наступления теплового равновесия цепи с термостатом рассматривали уже динамику изолированной цепи. Для увеличения точности определения функции корреляции использовали ее среднее значение по 500 различным реализациям начальной термализации цепи.

Численное моделирование динамики показало, что при температуре Т > 0.2 цепь ротаторов обладает конечной теплопроводностью. При этих значениях температуры удается численно проверить сходимость последовательности к(Ы1) при ЛГ, —- °о. Так, при Т = 0.3 с ростом /V, теплопроводность к(Л^) стремится к конечному значению к = 31.8 (рис. 10), т.е. при указанном значении температуры цепь обладает конечной теплопроводностью. При меньших температурах сходи-

Рис. 10. Зависимость коэффициента теплопроводности к от длины внутреннего участка цепи Л^. Точки а, б - полученные численно значения при температуре Т = 0.2 (Г+ = 0.21, Т_ = 0.19) и 0.3 (Т+ = 0.33, Т_ = 0.27). Аппроксимирующие прямые имеют наклон 8 = 0 (7) и 0.26 (2). Прямая 1 соответствует к = 31.8.

мость обеспечить не удается, но утверждать, что здесь наблюдается бесконечная теплопроводность, все же нельзя. Возможно, имеется более медленная сходимость, и для больших значений длины цепи N^ последовательность к(Л/,) уже имеет конечный предел. Тем не менее, при рассматриваемой температуре наблюдается резкое изменение поведения коэффициента теплопроводности.

Изучение поведения функции корреляции с(/) при / —- °° подтвердило утверждение о конечности теплопроводности цепи. При температуре Т >0.3 функция корреляции стремится к нулю с экспоненциальной скоростью, откуда следует сходимость интеграла в формуле (34) и, следовательно, конечность значения кс. Два использованных способа определения коэффициента теплопроводности дают очень близкие значения.

Вычисление теплопроводности по формуле Грина-Кубо (34) показало, что величина к монотонно уменьшается с увеличением температуры. Вид температурной зависимости приведен на рис. 11; при Т —► оо коэффициент теплопроводности к экспоненциально убывает.

.л_I_I_

1 2 3

Т

Рис. 11. Зависимость коэффициента теплопроводности к от температуры цепи Т. На вставке показано экспоненциальное убывание к при Т—►оо.

Рассмотрим теперь зависимость теплового потока I от разности температур АТ на примере конечной цепи из N = 300 звеньев с нулевой температурой на правом конце (Т_ = 0). Тогда разность температур АТ = Т+. Величину теплового потока здесь можно определить как работу силы трения на правом конце цепи:

N

Jo = l X (35)

n = N-N¡ + l

Численное интегрирование системы уравнений движения (29) показывает, что эта величина совпадает с локальным тепловым потоком во внутреннем участке цепи: /„ = J0 при п = = + 1, ..., N — /V,. Такое совпадение подтверждает корректность определения локального теплового потока (31).

Зависимость теплового потока У0 от разности температуры ДТ представлена на рис. 12. Как видно, тепловой поток монотонно растет с увеличением температуры. При температуре Т = Тг= 1.3 он достигает максимального значения (происходит насыщение). Дальнейшее повышение температуры приводит к постепенному уменьшению потока. Такое необычное поведение системы связано со специфическими свойствами фононов и с наличием в цепи при больших значениях темпера-

АТ

Рис. 12. Зависимость теплового потока /0 от разности температур АТ. N = 300, = 50, Т+ = АТ, 71 = 0.

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

Численное исследование системы уравнений движения показало, что в рассматриваемой цепи периодические волны (нелинейные фононы: фи(/) = ф(г), где волновая переменная г - п - 5 -скорость волны, а ф(г) - периодическая функция с периодом Ь) ограничены по амплитуде (для каждой длины волны с периодом Ь > 2 существует максимальное значение амплитуды А(Е)). Амплитуда периодической волны не может превышать его (А < А(Ь)). Удобно также определить амплитуду относительных поворотов молекул цепи Аф = = тах |ф(г)|, где функция ф(г) = ф(г + 1) - ф(г) задает профиль относительных поворотов. Вблизи максимального значения А(Ь) амплитуда относительных поворотов Ау достигает своего естественного максимального значения равного л.

С увеличением амплитуды энергия волны монотонной растет. Вблизи максимального значения амплитуды энергия волны также достигает наибольшего значения. Таким образом, в цепи связанных ротаторов периодические волны (колебательные моды) обладают свойством насыщения. Существует максимальное значение энергии Е(Ц, соответствующее температуре Т(Ь) = 1.8 [23,24]. При температуре Т0~ 1.8 в цепи термали-зованы практически все нелокализованные колебательные моды. Происходит насыщение колебательных мод, они становятся неустойчивыми и начинают сбрасывать лишнюю энергию. При

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

Ограниченность энергии междуузельного взаимодействия приводит к возможности существования в цепи локализованных ротационных мод, детально изученных для более общей системы в работе [29]. Ротационная мода соответствует вращению одной молекулы при почти неподвижных соседях. Динамика моды в первом приближении аппроксимируется вращением одной молекулы при неподвижных соседях. В этом приближении движение молекулы описывается уравнением маятника. При энергии Е > тах£/(ф) = 2 движение молекулы соответствует равномерному вращению. С ростом энергии Е частота вращения со будет монотонно увеличиваться от 0 до -н». Движение соседних молекул слабо влияет на вращение, которое однако проявляется со временем - при частотах меньше некоторой предельной величины. Численное моделирование вращения показало, что такое локализованное возбуждение будет устойчивым только при частоте со > со0 = 2.173. Зависимость энергии ротационной моды от частоты хорошо аппроксимируется параболой со2/2, соответствующей свободному вращению одной молекулы цепи. С ростом энергии теплоемкость моды монотонно уменьшается и при Е —► °° стремится к 1/2 (к теплоемкости изолированного ротатора).

Таким образом, в рассматриваемой системе кроме фононов (нелинейных кооперативных колебательных мод), имеющих частотный спектр 0 < со < 2, существуют локализованные ротационные моды с частотным спектром со > со0 > 2. Появление в цепи ротационных мод должно приводить к уменьшению теплоемкости системы.

Ротационные моды (локализованные вращения) являются стационарными возбуждениями цепи, они не могут участвовать в переносе тепловой энергии и должны приводить к уменьшению теплопроводности цепи. Изучим их взаимодействие с тепловыми фононами. Для этого рассмотрим цепочку из N = 300 молекул с температурой на левом конце Т+ = 0.02 > 0, а на правом Т_ = 0. Возьмем в начальный момент времени цепочку в

основном состоянии ф„ = 0, ф„ = 0 и сформируем в центре цепи на узле п — N/2 ротационную моду, имеющую энергию Е = 16. Рассмотрим теплопе-ренос в цепи. На рис. 13 показаны зависимость

? х 10"

Рис. 13. Взаимодействие ротационной моды с тепловыми фононами. а - зависимость распределения энергии Еп в цепи, б - зависимость энергии Е моды от времени

распределения энергии в цепи Еп и зависимость энергии моды (локализованного вращения) Е от времени Как видно, локализованная мода в центре цепи препятствует прохождению всех тепловых фононов. В результате происходит термали-зация только левой половины цепи - правая половина остается не термализованной (рис. 13а). Давление фононов на локализованную моду в центре цепи приводит к ее постепенному разрушению. Энергия моды монотонно уменьшается (рис. 136). В момент времени г = 28000 локализованное возбуждение разрушается, и в цепи возникает тепловой поток. Проведенное моделирование позволяет сделать вывод о невозможности формирования теплового потока при наличии локализованных возбуждений и о конечности времени их жизни в термализованной цепи.

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

10000

Рис. 14. Релаксация энергии термализованной цепи при температуре Т= 1 (а), 2 (б) и 3 (в). Дана зависимость локальной энергии Е от номера звена цепи п и времени I. В начальный момент времени система имеет температуру Г, далее рассматривается динамика цепи с поглощающими концами. N = 300, А/о = 50, Т+ = Т_ = 0.

глощение энергии (для этого достаточно в системе уравнений движения (29) положить Т+ = Т_ = 0). Численное моделирование динамики показало наличие в цепи локализованных ротационных мод. При температуре Т = 1 они обладают малым временем жизни, поскольку в результате взаимодействия с фононами быстро разрушаются, не препятствуя откачке энергии из цепи (рис. 14а). При температуре Т= 2 локализованные вращения характеризуются большим временем жизни и уже существенно влияют на механизм откачки энергии. Откачка происходит только после разрушения концевых локализованных мод. В результате этого процесса в цепи остается только одно локализованное состояние (рис. 146). Осо-

бенно хорошо заметна скачкообразность откачки энергии при температуре Т = 3 (рис. 14в).

Таким образом, теплоперенос в цепи можно качественно описать как последовательность случайных локальных перебросов энергии, происходящих при каждом разрушении стационарной локализованной моды. Этот механизм объясняет уменьшение теплопроводности цепи при увеличении температуры. Действительно, с повышением температуры растет лишь энергия локализованных возбуждений, увеличивается время их жизни и, следовательно, уменьшается частота локальных перебросов энергии. Этот механизм вместе с эффектом насыщения фононов объясняет также уменьшение потока энергии 30 с увеличением температуры при Т > Т0 (рис. 12). Теплоемкость цепи монотонно уменьшается с ростом температуры и стремится к значению 1/2, соответствующему теплоемкости системы несвязанных частиц. В пределе высоких температур рассматриваемая система превращается в цепь несвязанных ротаторов с нулевой теплопроводностью.

ТЕПЛОПРОВОДНОСТЬ ОДНОМЕРНОЙ РЕШЕТКИ С ПОТЕНЦИАЛОМ ПОДЛОЖКИ

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

Ж = ^\\мх1 + \к(хп + 1-х^ +и(хп)\, (36)

п '

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

Для удобства численного моделирования перейдем к безразмерным смещениям ип = 2шп/а

(а - шаг цепи), времени х = гЩм и энергии

Н = 4к2Ж/Ка2. Тогда функция Гамильтона (36) преобразуется к виду

Н = ц\ип+\(»п + 1-"п)2+У(ип)\ (37)

и '

Здесь штрих обозначает дифференцирование по безразмерному времени х, а безразмерный потенциал подложки V(un) = 4n2U(auJ2)/Ka2. Безразмерная температура T=4n2kBQ/Ka2, где кн - постоянная Больцмана, 6 - температура в размерных единицах.

Рассмотрим четыре наиболее характерные типа потенциала подложки:

гармонический потенциал

V(u) = (38)

потенциал sin-Gorgon

V(u) = е[1 + cos(m)], (39)

потенциал ф-4

V(u) = 2е[(м/я)2 - l]2 (40)

и потенциал sinh-Гордон

V(u) = (Do[ch(M)-l] (41)

Потенциалы sin-Гордон и ф-4 широко используются в нелинейной динамике. Параметр е > 0 определяет высоту потенциального барьера между основными состояниями, а его обратное значение g = 1/е характеризует кооперативность системы. Потенциалы (39) и (40) имеют одинаковое расстояние между соседними минимумами 2к и одинаковую высоту барьера 2е. Параметр (Do у потенциалов (38) и (41) соответствует минимальной частоте линейных колебаний решетки.

Теплопроводность цепи с периодическим потенциалом подложки (39) изучена в работах [30, 31], цепи с двухъямным потенциалом (40) - в работах [31,32], а с потенциалом (41) - в статье [33].

1 <п<^

Методы нахождения коэффициента теплопроводности

Нашей задачей является анализ теплопроводности на конечном участке молекулярной цепи из N звеньев. Для этого нужно рассмотреть бесконечную цепочку, левая часть которой (п < 0) погружена в термостат с температурой Т+, а правая часть (п > АО - в термостат с температурой 71 (Т+ > 71). Чтобы промоделировать процесс переноса тепла на конечном среднем участке цепи, возьмем конечную цепь из N. + N + звеньев, у которой первые звеньев помещены в термостат с температурой Т+, а последние А1 звеньев - в термостат с температурой 71 (рис. 15). Так как в рассматриваемых моделях используется гармонический потенциал междуузельного взаимодействия, характер граничных условий существенно не влияет на транспорт энергии по внутреннему участку цепи < п < + N. В силу симметричности потенциала междуузельного взаимодействия здесь не будет происходить тепловое расширение цепи, поэтому можно использовать как условие свободных концов (рис. 15а), так и условие закрепленных концов (рис. 156). Численное моделирование показало, что при = 40 и любом потенциале подложки граничные условия не оказывают заметного влияния на процесс переноса тепла по внутреннему участку цепи. Поэтому при моделировании теплопереноса в конечной цепи использовалось условие свободных концов с А/± = 40.

В большинстве работ для численного моделирования теплопроводности [30,32,34, 35] вводится детерминированный термостат Нозе-Хувера [25, 26] с = Л/_ = 1. Однако он предназначен для описания равновесного состояния термализован-ной системы и не всегда пригоден для описания неравновесного состояния. Поэтому будем использовать более адекватный физической постановке задачи стохастический термостат Ланжеве-на (как и в случае цепи ротаторов).

Рассмотрим конечную цепь со свободными концами (1 < п < N + + АО, у которой концевые А/± звеньев погружены в термостат Ланжевена с

(а) |

^ ......^ 1ЦЦ.~ ш».^ Шиь-уШШ^

Ои^ишищ^

Т —Т ' : т = т

т=т

У(и)

хмллллллллллллл/

Рис. 15. Модель конечной молекулярной цепи из + N + Ы_ звеньев с левыми звеньями, погруженными в термостат с температурой Т = Т+ и правыми N. звеньями в термостате с температурой Т= 71. Показана модель со свободными (а) и с закрепленными (б) концами. Потенциал подложки У(и) соответствует дискретной моделе Френкеля-Конторовой.

температурой Т±. Динамика цепи описывается системой уравнений движения

и'п = ип+1 -ип-Р(ип)-уип + С при п = 1 и"п = ип + 1-2ип + ип_1-Г(и„)-уип + Сп

при п = 2,..., Л^

"л = "л+1-2мн + м„-1-£■(«„)

при п = АГ+ + 1, + (42)

и» = ип+1-2ип + ип_1-Р(ип)-уип + ^'п

при п = + + Ы +

при п = А/+ + М + АГ_,

где функция Дм) = IШ(и)Ш, коэффициент трения у= 1/хг, хг- время релаксации скорости узлов цепи

погруженных в термостат, а ^ - моделирующий взаимодействие с термостатом гауссовский белый шум с условиями нормировки

<&!)> = <^(Х1)§1(Х2)> = 0

<^(х,)^(х2)> = 2уТ±дпк5(х2- х,)

Мгновенное значение безразмерной темпера-

,2

туры п-го узла цепи /и(х) = м„ (х). Чтобы найти ло-

0.010 -

0.005 -

Поэтому из уравнения непрерывности Кп = 7Я -]п_у следует, что мгновенный локальный поток энергии )п = -ип(ип+1-ип_ ,)/2.

Систему уравнений движения (42) интегрировали численно. Задавали значения у = 0.1, = 40, N = 10, 20, 40, 80, 160, 320, 640 и начальные условия, соответствующие основному состоянию цепи. За время х = 105 концевые звенья цепи приходили в тепловое равновесие с термостатом и на внутреннем участке цепи формировался стационарный тепловой поток. Дальнейшую динамику системы (42) рассматривали в течение времени х = 107. На соответствующем временном интервале рассчитывали средние значения локальных температур

т

Тп = </„(т))т = (45)

о

и средние значения теплового потока

Рис. 16. Распределение локального теплового потока /„ (а) и локальной температуры Тп (б) в цепи с периодическим потенциалом подложки (39). Параметр е = 1, N = 320, Ы± = 40, Т+ = 2.1,

Т_ = 1.9, время усреднения х = 107.

кальные значения потока энергии }п, рассмотрим распределение энергии по звеньям цепи

К = ^¿"п + и'п+1) + У(ип) + У(ми+1) +

+ ("л+1-"п)2]

(43)

Продифференцировав уравнение (43) по времени х, получим соотношение

I

К = ип(1))т = 11т \ (46)

Если разность температур АТ=Т+- Т_ мала, то данный способ термализации концов цепи позволяет избежать температурных скачков на границах термостатов [36]. Характерные распределения теплового потока /„ и локальной температуры Тп показаны на рис. 16. На внутреннем участке цепи < п < М+ + N локальный тепловой поток постоянен (7„ = У), а температурный профиль линеен. Это позволяет определить коэффициент теплопроводности, используя лишь внутренний участок цепи

к(Л0 = /(N-1)/(7у + 1-Гдг+дг)

(47)

К = ^{ип[и"п + Р(ип)] + м), + 1["» + 1 + Р(ип + 1)]} + + (ип+1-ип)(ип + 1-ип),

которое, учитывая систему уравнений движения (42), можно преобразовать к виду

К = + - ип) - ип(ип + 1 - ип_:)] (44)

Пусть прямая Т=ап + Ь наилучшим образом аппроксимирует зависимость Т = Тп на внутреннем участке цепи < п < + М, тогда точность вычисления коэффициента к(Л0 можно немного улучшить, положив к(Л0 = /а.

Предельное значение

к = Нт к(Л0

(48)

соответствует коэффициенту теплопроводности при температуре Т= (Т+ + Т_)/2. Здесь вопрос о конечности теплопроводности сводится к вопросу о существовании конечного предела (48). Если последовательность к(N) при N —► ос стремится к бесконечности, то при данном значении температуры цепочка обладает бесконечной теплопроводностью.

Другой способ нахождения коэффициента теплопроводности к связан с приведенной выше формулой Грина-Кубо, которую удобно переписать в виде

т

к = lim f lim —JJs(s)Js(0))ds, (49)

О

где N - число звеньев циклической цепи, полный тепловой поток

N п= 1

усреднение {•) выполняется по всем возможным термализованным состояниям цепи. Здесь вопрос о конечности теплопроводности сводится к вопросу о сходимости интеграла

Jc(x)dx, (50)

о

т.е. к оценке скорости убывания корреляционной функции

С(т) = lim —Ц(/5(х)/Д0))

Численно мы можем найти только корреляционную функцию для конечной цепи

Cw(x) = -2(J,(s)Us-x))s (51)

NT

Для достаточно больших значений N функция корреляции Cw(x) с высокой точностью аппроксимирует функцию C(t). Для численного анализа поведения корреляционной функции в большинстве случаев достаточно принять N = 4000.

Оба метода нахождения коэффициента теплопроводности (48), (49) дополняют друг друга. При

(б) 1.

" 1 ' I

I

0 80 160 240

п

Рис. 17. Распределение локального теплового потока Jn (а) и локальной температуры Тп (б) в цепи с гармоническим потенциалом подложки (38). Параметр «о = 1, N = 160, АГ± = 40, Г+ = 2.1, Т_ = 1.9, хг = 10, время усреднения х = 107.

корректном проведении вычислений мы должны получить одно и то же значение.

Теплопроводность цепи с гармоническим потенциалом подложки

Из рассматриваемых моделей цепочка с гармоническим потенциалом подложки (38) является самой простой. Соответствующая система уравнений движения линейна и поэтому вполне интегрируема. Транспорт энергии в данном случае происходит за счет подвижности не взаимодействующих между собой фононов, энергия которых не рассеивается. Тепловой поток ] не зависит от длины цепи N, а зависит лишь от разности температур А Т. Здесь не образуется наклонный температурный профиль, т. е. на внутреннем участке цепи локальная температура постоянна: Тп = (Г+ + Т_)/2 (рис. 17). Поэтому в силу определения (47), коэффициент теплопроводности к(А0 = т.е. система обладает бесконечной теплопроводностью, соответствующая функция корреляции С(х) не зави-

Рис. 18. Зависимость логарифма коэффициента теплопроводности 1п к(/У) (1,2) и к(Аг) (3,4) от логарифма длины внутреннего участка конечной цепи 1п/У (Л/+ = Ы_ = 40) с периодическим потенциалом подложки (39). Параметр е = 1, Т= 0.2 (1), 200 (2), 3 (3) и 20 (4). Точки - полученные численно значения, прямые - их наилучшие линейные аппроксимации.

сит от времени и, следовательно, поэтому интеграл (50) заведомо расходится.

Периодический потенциал подложки

Характер динамики цепи с периодическим потенциалом подложки (39) зависит от значения температуры. При низких температурах Т < е потенциал подложки хорошо аппроксимируется од-ноямным гармоническим потенциалом (38) с

соо = л/ё. Здесь транспорт энергии осуществляют слабо взаимодействующие между собой фононы. При температурах Т ~ е в цепи образуется хаотическая решетка топологических солитонов, и динамика приобретает совершенно иной характер. При высоких температурах Т > е происходит отрыв цепи от подложки, и динамика снова описывается лишь слабо взаимодействующими фонона-ми. Поэтому и теплопроводность цепи будет су-

1пС(т)

X

Рис. 19. Экспоненциальное убывание функции корреляции С(т) в цепи с периодическим потенциалом подложки (39). Параметр е=1,Г=3(7) и 20 (2).

щественно зависеть от значения приведенной температуры Т = Т/е.

Динамика цепи также сильно зависит от значения параметра кооперативности g = 1/е. Чем больше кооперативность, тем меньше будет плотность возникающей решетки топологических солитонов и тем слабее эффекты рассеяния на ней фононов. В пределе g —► (е —► 0) дискретная система уравнений движения преобразуется в континуальное уравнение вт-Гордон.

Таким образом, для трех предельных значений параметров дискретная модель Френкеля-Конто-ровой становится вполне интегрируемой системой: при Т —- 0 - цепью с гармоническим потенциалом подложки; при Т —► оо - изолированной гармонической цепью; при g —► оо - континуальной системой вш-Гордон. Во всех пределах система имеет бесконечную теплопроводность. Возникает естественный вопрос: что происходит вблизи этих пределов?

Положим вначале для определенности # = 1

(е = 1, Т = 7) и рассмотрим, как ведет себя последовательность к(А0 при увеличении N (Ы = 10, 20, 40,80,160,320,640) для различных значений Т. Из рис. 18 видно, что при малых (Г= 0.2) и больших (Т = 200) значениях температуры коэффициент теплопроводности к(А0 растет как №■, для Т-20-как \пЫ, а для Т=3 стремится к постоянному значению к = 18.5. Отсюда можно заключить только, что при Т = 3 цепочка имеет конечную тепло-

С(т)

XX 1СГ3

Рис. 20. Функции корреляции С(т) в цепи с периодическим потенциалом подложки (39). Параметр € = 1, Т = 0.2 (/) и 200 (2).

проводность. Для остальных трех значениях температуры остается вопрос, что будет при больших величинах N (вообще говоря, последовательность к(А0 может иметь конечный предел). Ответить на него однозначно путем прямого численного моделирования невозможно. Но здесь полезно рассмотреть поведение функции корреляции С(х) при х —- °о.

Численное моделирование показало, что при температуре Т= 3 функция корреляции стремится к нулю с экспоненциальной скоростью (рис. 19, кривая 1). Интеграл (22) сходится, и значение коэффициента теплопроводности, полученное по формуле Грина-Кубо (49), к = 17.5 согласуется с величиной к = 18.5, найденной прямым вычислением. При Т= 20 функция корреляции на интервале времени 0 < х < 800 также экспоненциально стремится к нулю (рис. 19, кривая 2). Если предположить, что эта зависимость сохранится и при х > 800, то по формуле Грина-Кубо мы получим значение коэффициента теплопроводности к = 77.4. Сравним эту величину с зависимостью к(Л0, представленной на рис. 18 (кривая 4). Здесь максимальное значение к(640) = 75.9, но, как видно из рисунка, тенденция к замедлению скорости роста к(Л0 не наблюдается, и вид зависимости свидетельствует скорее о ее бесконечном увеличении. Но для однозначного ответа на вопрос о конечности или бесконечности теплопроводности нужно еще найти к(М) при трех-четырех значениях N = 1280, 2560, 5120, 10240, для чего пока не хватает вычислительных мощностей современных компьютеров.

8

Рис. 21. Области в пространстве параметров (#, Т), где для конечных цепочек длины N < 640 с периодическим потенциалом подложки (39) (модель Френкеля-Конторовой) выполняется (заштрихованная область) и не выполняется (незаштрихованная) нормальная теплопроводность. Кривая 1 - граница раздела этих областей. Интервал 2 соответствует значениям параметров, рассмотренных в работе [30]. Для конечных цепочек (М < 640) с потенциалом подложки ф-4 (40) (модель ф-4) конечная теплопроводность наблюдается только в области, находящейся над прямой 3.

Еще труднее сделать однозначное утверждение о бесконечности теплопроводности при температурах Т=0.2 и 200. Вид функции корреляции для этих значений Г дан на рис. 20. Функция корреляции убывает очень медленно, что не позволяет однозначно определить характер убывания. Если полученную зависимость с(х) экстраполировать на времена х > 8000 экспонентой, то формула Грина-Кубо даст значения теплопроводности к = 1016 при Т = 0.2 и к = 2252 при Т = 200. Здесь тоже для однозначного ответа о бесконечности теплопроводности нужно, как минимум, найти к(А0 для четырех отмеченных значений длины цепи N. Но, с другой сторойы, при Т = 200 и N = 640 логарифм теплопроводности 1пк(Л/) = = 7.9 > 1п(2252) = 7.7, а вид зависимости 1пк(Л0 (рис. 18, кривая 2) не свидетельствует о замедлении скорости роста, т.е. эта зависимость тоже указывает скорее на бесконечное значение коэффициента теплопроводности.

к

Т

Рис. 22. Зависимость коэффициента теплопроводности к от приведенной температуры Т = Г/е для цепи с периодическим потенциалом подложки (39) при е = 3 (7), 5 (2) и 10 (3).

Рассмотрим теперь поведение конечной последовательности к(Л0 (М = 10,20,40,80,160,320,64) при других значениях параметра кооперативнос-ти. Результаты численного анализа даны на рис. 21.

Пространство параметров (#, Т) разделяется на две области. В первой (заштрихованной) области последовательность к(Л0 сходится (к(160) = = к(320) = к(640)), а во второй (незаштрихован-ной) области - монотонно растет. Таким образом, в первой области параметров дискретная цепочка вт-Гордон имеет конечную теплопроводность, а во второй - либо бесконечную, либо конечную, но достаточно высокую теплопроводность (здесь для конечных цепочек длины N < 640 не выполняется закон Фурье).

Первая область ограничена по оси g. Существует конечное значение параметра кооператив-ности g0 > 1 такое, что при g > g0 конечная последовательность к(ЛГ) монотонно растет при всех значениях температуры. Это связано с тем, что при сильной кооперативности дискретная модель яш-Гордон становится близкой к континуальной, являющейся вполне интегрируемой системой и обладающей бесконечной теплопроводностью. При каждом фиксированном значении коопера-

Т

Рис. 23. Зависимость безразмерной теплоемкости с от приведенной температуры Т - Т/е для цепи с периодическим потенциалом подложки (39) (1-5) и для цепи с потенциалом ф-4 (40) (6-11). е = 10 (1, 6), 5 (2, 7), 3 (3,8), 1 (4, 9) и 0.5 (5, 10). Штриховая кривая 11 - зависимость для цепи с потенциалом подложки втЬ-Гордон (41), частота Сйц = 1.

тивности g < g0 для конечных цепей длины N < 640 закон Фурье выполняется только в конечном интервале температур 0 < Ть < Т < Ти < При уменьшении кооперативности (g —0) верхняя граница этого интервала стремится к бесконечности (Ти —- а нижняя - к нулю (Ть —- 0) пропорционально g.

Вид зависимости коэффициента теплопроводности к от приведенной температуры Т показан на рис. 22. Внутри интервала [Ть, Тн] существует

критическая температура Тт, при которой коэффициент теплопроводности имеет минимальное

значение. При Т < Тт теплопроводность к моно-

тонно уменьшается, а при Т > Тт - монотонно увеличивается.

Для выяснения механизма реализации" нормальной теплопроводности проанализируем зависимость теплоемкости системы с = (Н)/ЫТ, где (Н) - среднее значение энергии циклической цепи из N молекул при температуре Т, от приведенной

температуры Т (рис. 23). Гармоническая цепочка всегда имеет теплоемкость, равную единице. Следовательно, отличие указанной величины от единицы может служить мерой нелинейности динамики при данном значении температуры. Для рассматриваемой модели потенциал подложки имеет отрицательный энгармонизм, поэтому теплоемкость системы должна быть больше единицы для всех значений температуры. Теплоемкость стремится к единице при Т —► 0 и Г —► <» и имеет только один максимум при температуре Тс. Данная температура совпадает с температурой Тт, при которой система имеет наименьшую теплопроводность.

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

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

Совпадение значений Тт и Тс позволяет заключить, что здесь нормальная теплопроводность обусловлена рассеиванием фононов на со-литонной решетке и зависит как от плотности этой решетки, так и от способности одного соли-

£

Рис. 24. Области в пространстве параметров (§, Т), где для конечных цепочек длины N < 640 с периодическим потенциалом подложки (39) и

периодом / = 2 л л/2 выполняется (заштрихованная область) и не выполняется (незаштрихо-ванная) нормальная теплопроводность. Для сравнения штриховой линией дана граница этих областей для модели Френкеля-Конторо-вой с соразмерными периодами (I = 2к).

тона к рассеиванию фононов. При сильной коо-перативности g > g0 взаимодействие между соли-тонами и фононами близко к упругому, вследствие чего теплопроводность перестает сходиться. С уменьшением кооперативности ширина соли-тонов уменьшается, плотность их решетки возрастает, взаимодействие фононов и солитонов становится все более неупругим. Появляется конечный интервал температур [Ть, Ти], в котором наблюдается сходимость теплопроводности. Для

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

Рассмотрим теперь несоизмеримую модель Френкеля-Конторовой, в которой период решетки не совпадает с периодом потенциала подложки. Пусть безразмерный потенциал подложки является периодической функцией (39) с периодом 2я, а цепь имеет безразмерный период I = 2щ.

1пС(х)

X

Рис. 25. Экспоненциальное убывание функции корреляции С(х) в цепи с потенциалом подложки ф-4 (40). Параметр е = 1, Т= 20 (/), 10 (2), 5 (5) и 3(4).

Тогда в системе уравнений движения (42) функция Ди„) будет иметь вид

Р(ип) = £и(ип + п1)

Для определенности примем, что отношение

периодов <7 = //2л = 72, когда период решетки / несоизмерим с периодом потенциала подложки 2л. Как хорошо известно [37, 38], в такой системе решетка солитонов существует при любом значении температуры, что способствует образованию нормальной теплопроводности.

На рис. 24 в пространстве параметров (#, Т) показана область значений (заштрихована), в которой наблюдается сходимость конечной последовательности к(А0, N = 10, 20, 40, 80, 160, 320, 640. Для сравнения штриховой линией выделена эта же область при совпадении периодов (/ = 2л). Рисунок демонстрирует, что несоизмеримость периодов не приводит к качественному изменению ситуации с теплопроводностью системы. Наблюдается лишь сдвиг вниз по температурной шкале

Т области с нормальной теплопроводностью. Обнаружение нормальной теплопроводности при более низких температурах связано с тем, что здесь уже нет необходимости в появлении решетки солитонов (она существует при всех температурах). Но наличия одной только солитонной ре-

1пС(х)

1пт

Рис. 26. Степенное убывание функции корреляции С(х) в цепи с потенциалом подложки ф-4 (40). Параметр е = 1, Т= 1 (1) и 0.5 (2). Наклон асимптот (прямых линий) а соответствует скорости степенного убывания х-01. Степень а = 1.2 (/) и 1.02(2).

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

Теплопроводность цепи с двухъямным потенциалом подложки

Рассмотрим теплопроводность цепи с потенциалом подложки ф-4 (40). Здесь анализ поведения конечной последовательности к(ЛГ), N = 10,20,40, 80, 160, 320, 640 показывает, что цепь обладает нормальной теплопроводностью при температуре Т > Т0 = 3#/2 (Т> 1.5) (рис. 21).

Для выяснения характера теплопроводности при температурах Т < То рассмотрим зависимость поведения функции корреляции С(х) от значения температуры. Примем для определенности параметр кооперативности # = 1 (е= 1). Поведение функции корреляции при высоких температурах представлено на рис. 25. При х —- <» функция корреляции с экспоненциальной скоростью стремится к нулю. Скорость убывания монотонно растет с повышением температуры. Это

подтверждает вывод о конечности теплопроводности при Т > То. При низких температурах функция корреляции уже убывает как показательная функция Т"® (рис. 26). С уменьшением температуры показатель а монотонно уменьшается (чем меньше температура, тем медленнее убывает функция корреляции). Так, при Т = 1 показатель а = 1.2 > 1, поэтому интеграл (50) сходится, и цепь обладает конечной теплопроводностью. При температуре Т = 0.5 показатель а= 1.02, здесь цепь также должна обладать конечной теплопроводностью. Но при Т < 0.5 показатель а < 1, интеграл (50) расходится, и цепь имеет бесконечную теплопроводность. Таким образом, для дискретной модели ф-4 можно утверждать,

что при низких температурах ( Т < 0.5#) цепь имеет бесконечную теплопроводность, а для высоких ( Т > 0.5g) - конечную.

Вид зависимости коэффициента теплопроводности к от приведенной температуры Т показан на рис. 27. Легко заметить, что при сильной кооперативное™ (при малых значениях параметра е) температурная зависимость для дискретной модели ф-4 существенно отличается от таковой для дискретной модели вт-Гордон. Совпадение наблюдается только при слабой кооперативное™ (^ < 0.5) и низкой температуре. Здесь с ростом температуры теплопроводность монотонно уменьшается и затем начинает увеличиваться. Но для модели ф-4 теплопроводность достигает локального максимума (Т = Т\), затем монотонно

убывает и при Т —► °о стремится к нулю. С ростом кооперативности относительная величина локального максимума теплопроводности уменьшается, и при сильной кооперативности максимум исчезает (точка локального максимума становится точкой перегиба).

Для выяснения механизма реализации нормальной теплопроводности полезно рассмотреть зависимость безразмерной теплоемкости с от

приведенной температуры цепи Т (рис. 236). При

Т —► 0 теплоемкость с —► 1. С повышением температуры теплоемкость монотонно растет, достигает своего максимального значения при

температуре Тс, затем монотонно уменьшается и становится меньшей единицы. Температура мак-

к

Рис. 27. Зависимость коэффициента теплопроводности к от приведенной температуры Т для цепи с двухъямным потенциалом подложки (40) при е = 4 (/), 2 (2) и 1 (5) и для цепи с потенциалом подложки ятЬ-Гордон (41) при частоте

(00 = 1 (4). Т = Г/« (1-3) и Т (4).

симума теплоемкости Тс оказывается близкой к

температуре минимума теплопроводности Ту. Такое поведение теплоемкости связано с особенностью потенциала ф-4. При низких температурах в первую очередь проявляется отрицательный ангармонизм, обусловленный двухъямнос-тью потенциала, в результате чего теплоемкость системы становится больше единицы, а при высоких (Т > 1) - положительный ангармонизм, вследствие которого теплоемкость системы становится меньше единицы.

Рассмотрим распределение по частотам тепловых колебаний цепи. Примем для определенности е = 4 (^ = 1/4) и найдем распределение для трех наиболее характерных значений температуры: Т = = 0.4,10,100. Для цепи с гармоническим потенциалом подложки (38) вид распределения не зависит от температуры и имеет форму

£(<о) = 2О)/71,/((Й2-(ОоК«1-«2), (52)

где максимальная частота со,2 = 4 + «о. При е- 4

значения частоты следующие: со0 = 4/ял/е = 2.546, со, = 3.238. Как видно из рис. 28а, при температуре

Е

ш

Рис. 28. Распределение по частотам ш энергии тепловых колебаний Е в цепи с двухъямным потенциалом подложки (40) при температуре Т= 0.4 (а), 10(б) и 100 (в). Параметре = 4. Штриховой линией показано распределение энергии в цепи с гармоническим потенциалом (38) с

(Оо = 4/яТЁ.

7 = 0.4 распределение тепловых колебаний по частотам в цепи ф-4 практически совпадает с распределением в гармонической цепи (52). Это свидетельствует о том, что при данном значении температуры колебания цепи остаются практически гармоническими, термализуются только фоно-ны. При Т= 10 > е распределение сдвигается влево за наименьшую частоту линейных колебаний (йо (рис. 286). Эта низкочастотная составляющая может быть связана с колебаниями солитонной решетки. Для более высокой температуры Т = 100 §> е распределение уже существенно сдвигается вправо за максимальную частоту линейных колебаний (Й! (рис. 28в). Сдвиг может быть объяснен только термализацией высокочастотных дискретных бризеров.

Таким образом, для низких температур Т < То = 0.5g динамика системы близка к линейной. Транспорт энергии реализуется за счет подвижности слабо взаимодействующих фононов, что приводит к расходимости теплопроводности. При

более высоких температурах Т > То теплопро-

,САВИН 1пС(х)

X

Рис. 29. Экспоненциальное убывание функции корреляции С(х) в цепи с потенциалом подложки втЬ-Гордон (41), частота соц = 1, при Т-10 (7), 7 (2) и 5 (5).

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

высоких температур Т >Т\ -с наличием высокочастотных дискретных бризеров. Следовательно, в цепи ф-4 имеется два механизма реализации нормальной теплопроводности: солитонный (рассеивание фононов на термализованной решетке топологических солитонов) и бризерный (рассеивание фононов на локализованных дискретных бризерах).

Теплопроводность цепи с потенциалом подложки зтИ-Гордон

Теплопроводность данной цепи была детально исследована в работе [33]. В этой модели потенциал подложки (41) является одноямной симметрической функцией с положительным энгармонизмом. Здесь анализ поведения конечной последовательности к(А0, N = 10, 20, 40, 80 160, 320, 640 тоже показывает, что цепь имеет конечную теплопроводность при температуре, превышающей некоторое пороговое значение (Т > Т0> 0). Это подтверждается тем, что при высоких температурах функция корреляции С(т) стремится к нулю с экспоненциальной скоростью при х —► °° (рис. 29).

Коэффициент теплопроводности монотонно уменьшается с ростом температуры и при Т—► о« с экспоненциальной скоростью стремится к нулю (рис. 27, кривая 4). Положительный энгармонизм потенциала приводит к тому, что с повышением температуры теплоемкость системы монотонно уменьшается (рис. 23, кривая 11). Исследование частотного спектра тепловых колебаний показывает, что с увеличением температуры спектр монотонно сдвигается в область высоких частот. Все это позволяет говорить о том, что для данной цепи имеется только один механизм реализации нормальной теплопроводности - рассеяние энергии на дискретных высокочастотных бризерах. Бризеры - локализованные высокочастотные нелинейные колебания цепи. Они малоподвижны и поэтому сами не участвуют в переносе тепловой энергии, но являются эффективными рассеивате-лями энергии фононов. Сдвиг частотного спектра тепловых колебаний выше максимальной частоты линейных колебаний (0] свидетельствует о том, что бризеры эффективно аккумулируют энергию тепловых колебаний. С повышением температуры растет концентрация бризеров, что и приводит к монотонному убыванию коэффициента теплопроводности.

Цепочка с потенциалом подложки

У(и) = ри4/4 (53)

(модель ф4) также обладает нормальной теплопроводностью [32, 39]. Потенциал (53), как и потенциал втЬ-Гордон (41), является одноямной симметричной функцией с положительным энгармонизмом. Поэтому здесь нормальная теплопроводность также связана с существованием дискретных бризеров, и коэффициент теплопроводности к(7) —»- 0 при Т —► °о. При р = 2 теплопроводность к(7) ~ Г"135 [39].

ЗАКЛЮЧЕНИЕ

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

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

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

В модели Френкеля-Конторовой (дискретный вт-Гордон) при исследовании конечных цепей ограниченной длины сходимость теплопроводности наблюдается только для определенной области значений кооперативное™ и температуры системы. Для высоких и низких температур или сильной кооперативное™ теплопроводность либо очень медленно сходится, либо расходится (стандартные численные методы не позволяют здесь дать однозначный ответ). Но есть основания полагать, что теплопроводность расходится, хотя область в пространстве параметров, при которых цепь имеет конечную теплопроводность, существенно больше области, приведенной на рис. 21.

В отличие от модели Френкеля-Конторовой для моделей ф-4 и втЬ-Гордон удается показать, что при низких температурах цепь характеризуется бесконечной теплопроводностью, а при высоких - конечной. Здесь существует критическая температура, при которой наблюдается переход цепи из состояния с бесконечной теплопроводностью в состояние с конечным ее значением.

Результаты проведенного анализа позволяют высказать предположение, что любая одномерная решетка (с потенциалом подложки или без

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Fermi Е., Pasta J., Ulam S. Los Alamos Report No. LA-1940. USA, 1955.

2. Peierls R.E. Quantum Theory of Solids. London: Oxford Univ. Press, 1955.

3. NoidD.W., Sumpter B.G., Wunderlich B. // Macromole-cules. 1991. V. 24. P. 4148.

4. Sumpter B.G., Noid D.W., Liang G.L., Wunderlich B. // Adv. Polym. Sci. 1994. B. 116. S. 29.

5. Manevitch LI. // Nonlinear Dynamics. 2001. V. 25. P. 95.

6. Manevitch L.I., Savin A.V. I I Phys. Rev. E. 1997. V. 55. №4. P. 4713.

7. Savin A.V., Manevitch LI. // Phys. Rev. B. 1998. V. 58. №17. P. 11386.

8. Savin A.V., Manevitch L.I. // Phys. Rev. E. 2000. V. 61. № 6. P. 7065.

9. Savin A.V., Manevitch LI. // Phys. Rev. B. 2001. V. 63. № 224303.

10. Савин A.B., Маневич Л.И., Христиансен ПЛ., Зо-лотарюк А.В. // Успехи физ. наук. 1999. Т. 169. № 3. Р. 255.

11. Sievers A J., Takeno S. // Phys. Rev. Lett. 1988. V. 61. № 8. P. 970.

12. MacKay R.S., Aubry S. // Nonlinearity. 1994. V. 7. P. 1623.

13. Aubry S. // Physica D. 1997. V. 103. P. 201.

14. Flach S., Willis СЯ. // Phys. Rep. 1998. V. 295. P. 181.

15. Kopidakis G., Aubry S. // Phys. Rev. Lett. 2000. V. 84. № 15. P. 3236.

16. Kopidakis G., Aubry S., Tsironis G.P. // Phys. Rev. Lett. 2001. V. 87. № 165501.

17. LeitnerDM.//Phys. Rev. Lett. 2001. V. 87. № 188102.

18. Kopidakis G., Aubry S. // Physica B. 2001. V. 296. P. 237.

19. Savin A.V., Manevitch LI. // Phys. Rev. B. 2003. V. 67. № 144302.

20. Маневич Л.И. // Высокомолек. соед. С. 2001. V. 43. № 12. С. 2215.

21. KirkwoodJ.G. //J. Chem. Phys. 1939. V. 7. P. 506.

22. PitzerK. //J. Chem. Phys. 1940. V. 8. P. 711.

23. Савин A.B., Гендельман О.В. // Физика твердого тела. 2001. Т. 43. № 2. С. 341.

24. Gendelman O.V., Savin A.V. // Phys. Rev. Lett. 2000. V. 84. №11. P. 2381.

25. Nose S. //J. Chem. Phys. 1984. V. 81. P. 511.

26. Hoover W.G. // Phys. Rev. A. 1985. V. 31. P. 1695.

27. FillipovA.,Hu В., Li В., ZeltserA. //J. Phys., Math. Gen. 1998. V. 31. P. 7719.

28. Kubo R., Toda M., Hashitsume N. I I Statistical Physics П, Springer Ser. Solid State Sci. 1991. V. 31.

29. TakenoS.,PeyrardM.HPhysicaD. 1996.V.92.P. 140.

30. Ни В., Li В., Zhao H. // Phys. Rev. E. 1998. V. 57. № 3. P. 2992.

31. Savin A.V., Gendelman O.V. // Phys. Rev. E. 2003. V. 67. № 041205.

32. Ни В., Li В., Zhao H. // Phys. Rev. E. 2000. V. 61. № 4. P. 3828.

33. Tsironis G.P., Bishop A.R., Savin A.V., Zolota-ryukA.V. Ц Phys. Rev. E. 1999. V. 60. № 6. P. 6610.

34. LepriS.,LiviR., PolitiA. //Phys. Rev. Lett. 1997. V. 78. № 10. P. 1896.

35. Lepri S., Livi R., Politi A. // Physica D. 1988. V. 119. P. 140.

36. Aoki K, Kusnezov D. I I Phys. Rev. Lett. 2001. V. 86. № 18. P. 4029.

37. Покровский ВЛ., Талапов АЛ. // Журн. эксперим. и теорет. физики. 1978. Т. 75. № 3. С. 1151.

38. Theodorou G., Rice Т.М. // Phys. Rev. В. 1978. V. 18. № 6. P. 2840.

39. Aoki К., Kuznezov D. // Phys. Lett. A. 2000. V. 265. P. 250.

Nonlinear Modes and Energy Transfer in Polymer Chains L. I. Manevich and A. V. Savin

Semenov Institute of Chemical Physics, Russian Academy of Sciences, ul. Kosygina 4, Moscow, 119991 Russia

Abstract—An analytical and numerical study on the dynamics of nonlinear excitations in molecular chains was performed. It was shown that the existence and propagation of random excitations in terms of the model used strongly depend on the properties of nonlinear normal modes. For example, heat bath-initiated nonlinear oscillations in the attenuation zone correspond to localized nonlinear normal modes, the breathers. In the breather localization region, periodic compression (tension) of interatomic bonds takes place, which is accompanied by a decrease (increase) in bond angles. It was found that breathers could exist in a thermalized chain, and the temperature dependence of their contribution to heat capacity was calculated. The heat transfer process in a one-dimensional lattice of coupled rotators, in which orientation interactions of neighboring units is described by a periodic potential, was studied. It was shown that a finite thermal conductivity in this lattice is due to the existence of localized stationary nonlinear excitations (rotobreathers), which scatter phonons. Using the Frenkel-Kontorova, ф-4, and sinh-Gordon models as examples, it was shown that the behavior of thermal conductivity of a harmonic one-dimensional lattice model with a nonparabolic form of substrate potential is determined by the spectrum of nonlinear excitations in the chain. Nonlinear excitations determine a particular mechanism of energy dissipation, which provide convergence of thermal conductivity. This is phonon scattering by the theimalized lattice of topological solitons for the sin-Gordon and ф-4 models at low temperatures and by localized high-frequency breathers for the sinh-Gordon and ф-4 models at high temperatures (in the case of the phi-4 model, the mechanism of energy scattering changes with a rise in temperature).