К теории динамической восприимчивости мягких магнитных коллоидов

Full Text

ВВЕДЕНИЕ

Коллоидно-дисперсные магнитополимерные композиты (феррогели, магнитные эластомеры) состоят из полимерной среды и внедренных в нее нано- или микроразмерных магнитных частиц. Эти многофункциональные интеллектуальные материалы привлекают интерес исследователей и практиков благодаря богатому набору уникальных физических свойств, важных для многих существующих и перспективных высоких технологий. Обзоры работ по этим направлениям можно найти, например, в [1–14]. Эти системы используются как материалы для изготовления магнито-полимерных актуаторов и стабилизаторов механических напряжений; нано- и микрокапсул для транспорта и высвобождения лекарств в организме; для создания магнитоуправляемых матриц роста, инженерии и регенерации биологических тканей. Внедрение магнитных наночастиц в больные клетки используется в магнито-гипертермическом методе терапии раковых и других опухолевых заболеваний.

Динамическая реакция магнитных гелей и подобных им мягких композитных сред на меняющееся магнитное поле является одной из наиболее интересных, с научной точки зрения, и важных в прикладном смысле характеристик этих систем. Для частиц в жидких магнитных коллоидах (феррожидкостях), а также частиц, иммобилизованных в достаточно жестком материале (твердые магнитные коллоиды), задача о динамическом перемагничивании полностью иммобилизованной частицы рассматривалась в работах [15–18] на основе различных приближенных методов решения уравнения Фоккера–Планка для плотности вероятности различных ориентаций магнитного момента частицы. Теоретические модели динамической восприимчивости композита с неподвижными взаимодействующими феррочастицами были предложены в [19, 20].

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

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

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

2. ОСНОВНЫЕ ПРИБЛИЖЕНИЯ И МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Рассмотрим ферромагнитную однодоменную частицу, помещенную в магнитное поле H. Схематическая ее иллюстрация приведена на рис. 1. Обозначим µ – единичный вектор, направленный вдоль магнитного момента частицы; ν – единичный вектор вдоль оси легкого намагничивания частицы; ν₀ – начальное, до приложения поля, значение вектора ν.

Мы предполагаем, что частица находится в немагнитной упруго-вязкой среде. Под упруго-вязкой мы будем понимать среду, которая при малых статических нагрузках ведет себя как упругое тело (эффекты ползучести на больших временах не учитываем), а при динамических нагрузках механическая энергия в ней диссипирует. Необходимо отметить, что мягкие материалы чрезвычайно разнообразны по своим физическим и реологическим свойствам и сколько-нибудь общее уравнение их реологического состояния в литературе отсутствует. Поэтому, чтобы сконцентрировать внимание на принципиальной стороне проявления упругих и диссипативных эффектов при динамической деформации среды, мы будем описывать ее в рамках простейшего уравнения Кельвина-Фойгхта $\mathrm{\Sigma }=G\gamma +\eta \dot{\gamma },$ где ∑ – напряжение, G – статический модуль упругости, η – вязкость среды, γ – величина ее деформации.

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

Учитывая это, введем нормализованную на единицу функцию распределения (плотность вероятности) w(µ,ν|β) по ориентациям µ и ν при заданном ν₀, где β – угол между векторами H и ν₀ (рис. 1). В рамках приближения Кельвина–Фойгхта реологического поведения несущей среды функция w может быть найдена как решение уравнения Фоккера-Планка:

$\frac{\partial w}{\partial t }=-\frac{1}{2{\tau }_{\text{D}}}{\text{div}}_{\mu }{j}_{\mu }-\frac{1}{2{\tau }_{\eta }}{\text{div}}_{\nu }{j}_{\nu }.$ (1)

 

Рис. 1. Иллюстрация феррочастицы и используемой системы координат.

 

Здесь и далее τ_D – характерное время вращательной диффузии магнитного момента частицы относительно ее тела, определяемое внутренней кристаллической структурой частицы; ${\mathrm{\tau }}_{\mathrm{\eta }}=3\mathrm{\eta }{v}_{\text{p}}\text{/}kT$ – характерное время вращательной диффузии тела частицы [22] в несущей среде; v_p – объем частицы; Ψ – полярный угол между векторами ν и ν₀, показанный на рис. 1; θ – полярный угол между вектором момента µ и вектором оси легкого намагничивания ν; нижние индексы µ и ν отмечают операторы и векторы, относящиеся к соответствующим векторам; j_µ и j_ν – безразмерные плотности потоков вероятности изменения векторов µ и ν соответственно. Они могут быть записаны в следующем виде (см. аналогичные формы в [16–18]:

$\begin{array}{c}{j}_{\mu }=-\left(w{\text{grad}}_{\mu }u+{\text{grad}}_{\mu }w\right),\\ j_{\nu }=-\left(w{\text{grad}}_{\nu }u+{\text{grad}}_{\nu }w\right),\end{array}$ (2)

где u – безразмерная (отнесенная к тепловой энергии kT) энергия частицы:

$\mathrm{u}\left(\mathrm{\theta },\mathrm{\psi },\mathrm{\phi }\right)={\mathrm{u}}_{\text{h}}\left(\mathrm{\theta },\mathrm{\phi }\right)+{\mathrm{u}}_{\mathrm{\sigma }}\left(\mathrm{\theta },\mathrm{\psi }\right)+{\mathrm{u}}_{\text{e}}\left(\mathrm{\psi }\right).$ (3)

Первый член в правой части (3) – безразмерная энергия Зеемана взаимодействия частицы с полем Н, равная:

$\begin{array}{c}{u}_{\text{h}}=-\left(\mathbf{\mu }\cdot \mathit{h}\right)=-h{\mathrm{\mu }}_{\text{h}} ,\\ \mathit{h}={\mathrm{\mu }}_0\frac{m}{kT}\mathbf{H}.\end{array}$ (4)

Здесь µ₀ – магнитная проницаемость вакуума; µ_h – проекция вектора µ на поле H. Параметр h может рассматриваться как безразмерное напряжение поля H.

Безразмерная энергия магнитной анизотропии частицы в известном приближении Стонера-Вольфара может быть записана в виде:

$\begin{array}{c}{u}_{\sigma }=-\sigma {\left(\mathbf{\mu }\cdot \mathbf{\nu }\right)}^2,\\ \sigma =\frac{K{v}_{\text{p}}}{kT},\end{array}$ (5)

где K – параметр этой анизотропии.

Потенциал u_e – безразмерная энергия упругой деформации несущей среды при повороте частицы на угол y. Мы оценим эту энергию в простейшем квадратичном приближении:

$\begin{array}{c}{u}_{\text{e}}=\frac{1}{2}g{\psi }^2,\\ g~\frac{G{v}_{\text{p}}}{kT},\end{array}$ (6)

где g – безразмерный модуль упругости несущей среды.

Отметим, что если вектор ν направления оси анизотропии частицы фиксирован (т. е. τ_h→∞ или g→∞), то уравнение (1) совпадает с хорошо известным уравнением Фоккера-Планка для ориентационного вектора µ момента полностью неподвижной ферромагнитной частицы (см., например, [15, 17, 18]). Ниже будет показано, что, когда неелевская переориентация вектора µ невозможна (σ→∞) и упругие силы в среде отсутствуют (g = 0), уравнение (1) совпадает с уравнением [16] вращения частицы в жидкости. Форма (1) является естественным обобщением этих уравнений на рассматриваемый “двухвекторный” процесс.

Мы ограничимся исследованием частиц с высокой магнитной анизотропией, для которых выполняются неравенства $\sigma \gg \mathrm{1,} \sigma \gg h.$ К примеру, для частиц окислов железа, часто используемых в различных приложениях, параметр магнитной анизотропии K ≈ 14 кДж/м³ [22]. Следовательно, при комнатных температурах неравенство σ > 10 выполняется, когда диаметр частицы d_p > 17–18 нм. Такой размер частиц часто встречается в экспериментах и практических применениях магнитных нанокомпозитов, во-первых, в силу практически неизбежной и часто широкой полидисперсности синтезируемых систем магнитных наночастиц (см., например, [22]); во-вторых, композиты с диаметром частиц, превышающим 17–18 нм, часто специально изготавливаются, чтобы получить сильный отклик композита на магнитное поле (см., например, [23, 24]). Отметим, что очень часто структурные особенности поверхности ферромагнитных наночастиц вносят дополнительный и очень существенный вклад в общую энергию их магнитной анизотропии и условие σ > 10 выполняется для частиц, размер которых меньше 10 нм (см., например, [25, 26]).

Учитывая, что намагниченность насыщения окислов железа М_p ≈ 400–500 кА/м (таблицы физических характеристик магнитных наночастиц можно найти в [22]), при учете только объемной магнитной анизотропии получаем, что неравенство σ > h выполняется при напряженности магнитного поля в диапазоне H < 22–27 кA/м. Отметим, что многие эксперименты (например, по биомедицинским применениям феррогелей) выполнены в этих рамках напряженности магнитного поля. Учет поверхностной магнитной анизотропии частиц может существенно расширить область полей, для которых выполняется условие σ > h .

Для дальнейшего удобно ввести систему координат с осью Ozʹ, направленной вдоль вектора ν, и осью Oxʹ в плоскости, образуемой векторами ν и H. Обозначим θ и φ полярный и азимутальный углы вектора µ в этой системе координат (рис. 2).

 

Рис. 2. Система координат с полярной осью, направленной вдоль вектора ν оси легкого намагничивания частицы. Азимутальный угол не показан для простоты. Знаки 1 и 2 обозначают полусферические потенциальные области, обсуждаемые в тексте.

 

Потенциал u_µ = u_h + u_σ в этой системе координат может быть представлен в виде:

$\begin{array}{c}{u}_{\mathrm{\mu }}=-h\left({\mathrm{\zeta }}_{}\text{cos}\mathrm{\theta }+{\mathrm{\zeta }}_{}\text{sin}\mathrm{\theta }\text{cos}\mathrm{\varphi }\right)-\mathrm{\sigma }{\text{cos}}^2\mathrm{\theta },\\ {\mathrm{\mu }}_{\text{h}}={\mathrm{\zeta }}_{}\text{cos}\mathrm{\theta }+{\mathrm{\zeta }}_{}\text{sin}\mathrm{\theta }\text{cos}\mathrm{\varphi },\end{array}$ (7)

где ζ = h/h – единичный вектор вдоль поля h. Поскольку выполняется сильное неравенство h << σ, минимальные значения потенциала u_µ приблизительно соответствуют θ = 0 и равны

$u_{\mu 1}=-h{\zeta }_{}-\sigma , u_{\mu 2}=h{\zeta }_{}-\sigma$ (8)

соответственно.

При заданном азимутальном угле f максимальное значение потенциала u_µ соответствует примерно полярному углу $\mathrm{\theta }=\mathrm{\pi }/2$ и равно

$u_{\mu \text{max}}=-h{\zeta }_{}\text{cos}\varphi .$ (9)

Член ${\text{div}}_{\mu }{j}_{\mu }$ уравнения (1) в системе координат (θ, φ) может быть представлен в виде:

${\text{div}}_{\mu }{j}_{\mu }=\frac{1}{\text{sin}\theta }\frac{\partial }{\partial \theta }{j}_{\mu \theta }+\frac{1}{\text{sin}\theta }\frac{\partial }{\partial \varphi }{j}_{\mu \varphi },$ (10)

где функции

$\begin{array}{c}{j}_{\mu \theta }\left(\theta ,\varphi \right)=-\text{sin}\theta \left(w\frac{\partial u}{\partial \theta }+\frac{\partial w}{\partial \theta }\right),\\ j_{\mu \varphi }\left(\theta ,\varphi \right)=-\frac{1}{\text{sin}\theta }\left(w\frac{\partial u}{\partial \varphi }+\frac{\partial w}{\partial \varphi }\right)\end{array}$

являются полярной и азимутальной компонентами вектора j_µ соответственно.

Мы будем отмечать все физические величины, относящиеся к областям “1” $\left(0<\mathrm{\theta }<\frac{\mathrm{\pi }}{2}\right)$ и “2” $\left(\frac{\mathrm{\pi }}{2}<\mathrm{\theta }<\mathrm{\pi }\right)$ на рис. 2 нижними индексами 1 и 2. Пусть P₁ и P₂ вероятности того, что момент частицы принадлежит областям “1” и “2” соответственно:

$\begin{array}{c}{P}_1\left(\mathrm{\psi }\right)={\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }\text{/}2}\mathrm{w}(\mathrm{\theta },\mathrm{\varphi }|\mathrm{\psi })\text{sin}\left(\mathrm{\theta }\right)\mathrm{d\theta }d\mathrm{\varphi },\\ P_2\left(\mathrm{\psi }\right)={\int }_0^{2\mathrm{\pi }}{\int }_{\mathrm{\pi }\text{/}2}^{\mathrm{\pi }}\mathrm{w}(\mathrm{\theta },\mathrm{\varphi }|\mathrm{\psi })\text{sin}\left(\mathrm{\theta }\right)\mathrm{d\theta }d\mathrm{\varphi }.\end{array}$ (11)

Чтобы получить уравнение для P₁(Ψ), мы проинтегрируем уравнение (1) так же, как w в первом соотношении (11). Учитывая периодическую зависимость компоненты j_µφ от угла φ, после вычислений получаем:

$\begin{array}{c}\frac{\partial {\mathrm{P}}_1}{\partial \mathrm{t} }=-\frac{1}{2{\mathrm{\tau }}_D}{\mathrm{J}}_{\mathrm{\mu }}-\frac{1}{2{\mathrm{\tau }}_{\mathrm{\eta }}}{\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }\text{/}2}{\text{div}}_{\mathrm{\nu }}{\mathrm{j}}_{\mathrm{\nu }}\text{sin}\left(\mathrm{\theta }\right)d\mathrm{\theta }d\mathrm{\varphi },\\ {\mathrm{J}}_{\mathrm{\mu }}={\int }_0^{2\mathrm{\pi }}{\mathrm{j}}_{\mathrm{\mu \theta }}\left(\mathrm{\pi }\text{/}\mathrm{2,}\mathrm{\varphi }\right)d\mathrm{\varphi }.\end{array}$ (12)

Физический смысл величины J_µ – поток вероятности перехода вектора µ из области “1” в область “2” через плоскость θ = π/2.

Чтобы оценить J_µ, мы используем идеи классической теории Крамерса диффузии броуновской частицы через потенциальный барьер. А именно, представим:

$j_{\mathrm{\mu \theta }}\left(\mathrm{\theta },\mathrm{\varphi }\right)=-\text{sin}\left(\mathrm{\theta }\right){\mathrm{e}}^{-\mathrm{u}\left(\mathrm{\theta },\mathrm{\varphi }\right)}\frac{\partial }{\partial \mathrm{\theta }}\left(\mathrm{w}\left(\mathrm{\theta },\mathrm{\varphi }\right){\mathrm{e}}^{\mathrm{u}\left(\mathrm{\theta },\mathrm{\varphi }\right)}\right).$ (13)

Далее, после простых преобразований уравнения (13), приходим к соотношению:

${\int }_0^{\mathrm{\pi }/2}\frac{{\mathrm{j}}_{\mathrm{\mu \theta }}\left(\mathrm{\theta },\mathrm{\varphi }\right)}{\text{sin}\left(\mathrm{\theta }\right)}{\mathrm{e}}^{\mathrm{u}\left(\mathrm{\theta },\mathrm{\varphi }\right)}\mathrm{d\theta }=-\left({\mathrm{w}}_2\left(\mathrm{\varphi }\right){\mathrm{e}}^{{\mathrm{u}}_2}-{\mathrm{w}}_1\left(\mathrm{\varphi }\right){\mathrm{e}}^{{\mathrm{u}}_1}\right),$ (14)

где

$\begin{array}{c}{u}_1=u\left(\mathrm{0,}\mathrm{\varphi }\right)=u_{\mathrm{\mu }1}+\frac{1}{2}{\mathrm{g\psi }}^2=-h{\mathrm{\zeta }}_{\mathrm{z}\text{'}}-\mathrm{\sigma }+\frac{1}{2}{\mathrm{g\psi }}^2,\\ u_2=u\left(\mathrm{\pi },\mathrm{\varphi }\right)=u_{\mathrm{\mu }2}+\frac{1}{2}{\mathrm{g\psi }}^2=h{\mathrm{\zeta }}_{\mathrm{z}\text{'}}-\mathrm{\sigma }+\frac{1}{2}{\mathrm{g\psi }}^2,\\ {\mathrm{w}}_1\left(\mathrm{\varphi }\right)=\mathrm{w}\left(\mathrm{0,}\mathrm{\varphi }\right); {\mathrm{w}}_2\left(\mathrm{\varphi }\right)=\mathrm{w}\left(\mathrm{\pi },\mathrm{\varphi }\right),\\ u_{\mathrm{\mu }1}=-h{\mathrm{\zeta }}_{\mathrm{z}\text{'}}-\mathrm{\sigma },u_{\mathrm{\mu }2}=h{\mathrm{\zeta }}_{\mathrm{z}\text{'}}-\mathrm{\sigma }.\end{array}$

Поскольку магнитная анизотропия частицы предполагается сильной (σ >> 1), функция exp(u(θ)) имеет острый максимум в точке θ_µax ≈ π/2. Используя стандартные соображения метода перевала, получаем:

$\begin{array}{c}{\int }_0^{\pi /2}\frac{j_{\mu \theta }\left(\theta \right)}{\text{sin}\left(\theta \right)}{e}^{u\left(\theta ,\varphi \right)}d\theta \approx \frac{j_{\mu \theta }\left(\frac{\pi }{2},\varphi \right)}{\text{sin}\left({\theta }_{\text{max}},\varphi \right)}{e}^{u_{\text{max}}}\sqrt{\frac{\pi }{\Omega }} \approx \\ \approx j_{\mu \theta }\left(\frac{\pi }{2},\varphi \right)e^{u_{\text{max}}}\sqrt{\frac{\pi }{\Omega }},\\ u_{\text{max}}=u_{\mu \text{max}}+\frac{1}{2}g{\psi }^2=-h{\zeta }_{}\text{cos}\varphi +\frac{1}{2}g{\psi }^2,\\ \Omega =-\frac{1}{2}\frac{{\partial }^{2}u}{\partial {\theta }^2}{|}_{\theta ={\theta }_{\text{max}}}\approx \sigma .\end{array}$ (15)

Комбинируя уравнения (12)–(15), приходим к следующему соотношению:

$\begin{array}{c}\frac{\partial P_1}{\partial t }=\frac{1}{2{\tau }_{\text{I}}}{\int }_0^{2\pi }\left(\begin{array}{c}{w}_2\left(\varphi \right)e^{h{\zeta }_{}+h{\zeta }_{}\text{cos}\varphi }-\\ -w_1\left(\varphi \right)e^{-h{\zeta }_{}+h{\zeta }_{}\text{cos}\varphi }\end{array}\right)d\varphi -\\ -\frac{1}{2{\tau }_{\eta }}{\int }_0^{2\pi }{\int }_0^{\pi /2}{\text{div}}_{\nu }{j}_{\nu }\text{sin}\left(\theta \right)d\theta d\varphi ,\\ {\tau }_{\text{I}}={\tau }_{\text{D}}\sqrt{\frac{\pi }{\Omega }}{e}^{\sigma } .\end{array}$ (16)

По порядку величины параметр τ_I – характерное время перехода момента частицы через потенциальный барьер магнитной анизотропии из области “1” в область “2”. Поскольку σ >> 1, то выполняется сильное неравенство τ_I >> τ_D. Заметим, что по порядку величины τ_D равно времени релаксации функции распределения w к своей равновесной форме в каждой из областей 1 и 2. Так как τ_I >> τ_D, следуя методу Крамерса, можно принять, что за время перехода момента частицы через потенциальный барьер в каждой из этих областей успевает установиться практически равновесное распределение по ориентациям момента, т. е. равновесное распределение по углам θ, φ. Заметим, что это не означает равновесного распределения между областями “1” и “2”. Таким образом, в областях “1” и “2” соответственно получаем:

$\begin{array}{c}{w}_1\left(\varphi \right)\approx {w_1}^0{e}^{u_1-u\left(\mathrm{0,}\varphi \right)}={w_1}^0,\\ w_2\left(\varphi \right)\approx {w_1}^0{e}^{u_2-u\left(\pi ,\varphi \right)}={w_2}^0,\\ w_{\mathrm{1,2}}^0=w_{\mathrm{1,2}}\left(\varphi =0\right).\end{array}$

Следовательно, уравнение (16) может быть переписано в виде:

$\begin{array}{c}\frac{\partial P_1}{\partial t }=\frac{1}{2{\tau }_{\text{I}}}Q\left(h,{\zeta }_{}\right)[{w_2}^0{e}^{h{\zeta }_{}}-{w_1}^0{e}^{-h{\zeta }_{}}] -\\ -\frac{1}{2{\tau }_{\eta }}{\int }_0^{2\pi }{\int }_0^{\pi \text{/}2}{\text{div}}_{\nu }{j}_{\nu }\text{sin}\left(\theta \right)d\theta d\varphi ,\\ Q\left(h,{\gamma }_{}\right)={\int }_0^{2\pi }{e}^{h{\zeta }_{}\text{cos}\varphi }d\varphi .\end{array}$ (17)

В рамках квазиравновесного приближения Крамерса в каждой потенциальной области выполняются следующие соотношения:

$\begin{array}{c}w\left(\theta ,\varphi \right)\approx w_1^0{e}^{u_1-u\left(\theta ,\varphi \right)}\approx \\ \approx w_1^0{e}^{-\sigma \left(1-{\text{cos}}^2\left(\theta \right)\right)-h\left({\zeta }_{}\left(1-\cos \left(\theta \right)\right)-{\zeta }_{}\sin \left(\theta \right)\text{cos}\varphi \right) },\\ 0<\theta <\frac{\pi }{2}\end{array}$ (18)

и

$\begin{array}{c}w\left(\theta ,\varphi \right)\approx w_2^0{e}^{u_2-u\left(\theta ,\varphi \right)}\approx \\ \approx w_2^0{e}^{-\sigma \left(1-{\text{cos}}^2\left(\theta \right)\right)+h\left({\zeta }_{}\left(1+\cos \left(\theta \right)\right)+{\zeta }_{}\sin \left(\theta \right)\text{cos}\varphi \right) },\\ \frac{\pi }{2}<\theta <\pi .\end{array}$ (19)

Используя эти соотношения, можно следующим образом оценить вероятности Р₁ и Р₂:

$\begin{array}{c}{P}_1={w_1}^0{q}_1, P_2={w_2}^0{q}_2,\\ q_1={\int }_0^{2\pi }{\int }_0^{\pi /2}{e}^{-\sigma \left(1-{\text{cos}}^2\left(\theta \right)\right)}\times \\ \times e^{-h\left({\zeta }_{}\left(1-\cos \left(\theta \right)\right)-{\zeta }_{}\sin \left(\theta \right)\text{cos}\varphi \right) }\text{sin}\left(\theta \right)d\theta d\varphi ,\\ q_2={\int }_0^{2\pi }{\int }_{\pi /2}^{\pi }{e}^{-\sigma \left(1-{\text{\hspace{0.17em}cos}}^2\left(\theta \right)\right)}\times \\ \times e^{h\left({\zeta }_{}\left(1+\cos \left(\theta \right)\right)+{\zeta }_{}\sin \left(\theta \right)\text{cos}\varphi \right)}\text{sin}\left(\theta \right)d\theta d\varphi .\end{array}$ (20)

Параметры q_1,2 без труда могут быть вычислены численно. Однако полезно также получить их аналитические оценки. Вследствие сильного неравенства s >> 1, экспонента ${\mathrm{e}}^{-\mathrm{\sigma }(1-{\text{cos}}^2(\mathrm{\theta }\left)\right) }$ имеет резкий максимум при θ = 0 и p. Поскольку s >> h, экспоненты ${\mathrm{e}}^{-\mathrm{h}\left({\mathrm{\zeta }}_{}\right(1-\cos \left(\mathrm{\theta }\right))-{\mathrm{\zeta }}_{}\sin (\mathrm{\theta }\left)\text{cos}\mathrm{\varphi }\right)}$ и ${\mathrm{e}}^{\mathrm{h}\left({\mathrm{\zeta }}_{}\right(1+\cos \left(\mathrm{\theta }\right))+{\mathrm{\zeta }}_{}\sin (\mathrm{\theta }\left)\text{cos}\mathrm{\varphi }\right)}$ меняются с углом θ намного медленнее, чем $e^{-\sigma (1-{\text{cos}}^2(\theta )}$. Учитывая это, используя стандартные соображения метода перевала, представив ${\cos }^2\mathrm{\theta }\approx 1-{\mathrm{\theta }}^2$ и $\mathrm{sin\theta }\approx \mathrm{\theta }$ вблизи θ = 0, можно получить оценку:

$\begin{array}{c}{q}_1={\int }_0^{2\pi }{\int }_0^{\pi /2}{e}^{-\sigma (1-{\text{cos}}^2(\theta \left)\right)}\times \\ \times e^{-h\left({\zeta }_{}\right(1-\cos \left(\theta \right))-{\zeta }_{}\sin \left(\theta \right)\text{cos}\varphi ) }\text{sin}\left(\theta \right)d\theta d\varphi \approx \\ \approx {\int }_0^{2\pi }{\int }_0^{\pi /2}{e}^{-\sigma {\theta }^2 }\theta d\theta d\varphi =\pi (1-e^{-\pi \sigma /2})\approx \frac{\pi }{\sigma }.\end{array}$ (21)

Аналогично

$q_2\approx \frac{\pi }{\sigma }.$ (22)

Используя уравнения (20)–(22) в уравнении (17), приходим к следующему уравнению:

$\begin{array}{c}\frac{\partial {w_1}^0}{\partial t }=\frac{1}{2{\tau }_{\text{II}}}\left[\right({w_2}^0{e}^{h{\zeta }_{}}-{w_1}^0{e}^{-h{\zeta }_{}}] -\\ -\frac{\sigma }{2\pi {\tau }_{\eta }}{\int }_0^{2\pi }{\int }_0^{\pi /2}{\text{div}}_{\nu }{j}_{\nu }\text{sin}\left(\theta \right)d\theta d\varphi ,\\ {\tau }_{\text{II}}=\frac{\pi {\tau }_{\text{I}}}{\sigma Q}\approx {\tau }_{\text{D}}\sqrt{\frac{1}{\pi {\sigma }^3}}\frac{e^{\sigma }}{Q}.\end{array}$ (23)

Используя аналогичные соображения, получаем:

$\begin{array}{c}\frac{\partial {w_2}^0}{\partial t}=-\frac{1}{2{\tau }_{\text{II}}}[{w_2}^0{e}^{h{\zeta }_{}}-{w_1}^0{e}^{-h{\zeta }_{}}] -\\ -\frac{\sigma }{2\pi {\tau }_{\eta }}{\int }_0^{2\pi }{\int }_{\pi \text{/}2}^{\pi }{\text{div}}_{\nu }{j}_{\nu }\text{sin}\left(\theta \right)d\theta d\varphi .\end{array}$ (24)

Рассмотрим сейчас интегральный член в (23). В системе координат, показанной на рис. 1, имеем:

$\begin{array}{c}{\int }_0^{2\pi }{\int }_0^{\pi /2}{\text{div}}_{\nu }{j}_{\nu }\text{sin}\left(\theta \right)d\theta d\varphi =\\ ={\int }_0^{2\pi }{\int }_0^{\pi /2}\frac{1}{\text{sin}\psi }\left(\frac{\partial }{\partial \psi }{j}_{\nu \psi }\right)\text{sin}\left(\theta \right)d\theta d\varphi =\\ ={\int }_0^{2\pi }\frac{1}{\text{sin}\psi }\frac{\partial }{\partial \psi }{\int }_0^{\pi /2}{j}_{\nu \psi }\text{sin}\left(\theta \right)d\theta .\end{array}$ (25)

Здесь и ниже Ψ и φ – полярный и азимутальный углы вектора ν в системе координат с полярной осью, направленной вдоль вектора ν₀.

Принимая во внимание уравнение (2), получаем:

$\begin{array}{c}{\int }_0^{\pi /2}{j}_{\nu \psi }\text{sin}\left(\theta \right)d\theta =\\ =-\text{sin}\psi {\int }_0^{\pi \text{/}2}{e}^{-u}\frac{\partial }{\partial \psi }\left(w{e}^u\right)\text{sin}\left(\theta \right)d\theta .\end{array}$ (26)

В системе координат (Ψ, φ) имеем:

$\begin{array}{c}{j}_{\nu \psi }\left(\theta ,\varphi \right)=-\text{sin}\psi \left(w\frac{\partial u}{\partial \psi }+\frac{\partial w}{\partial \psi }\right),\\ j_{\mu \phi }\left(\psi ,\phi \right)=-\frac{1}{\text{sin}\psi }\left(w\frac{\partial u}{\partial \phi }+\frac{\partial w}{\partial \phi }\right).\end{array}$ (27)

Комбинируя уравнения (26) и (14), приходим к соотношению:

$\begin{array}{c}{\int }_0^{\pi \text{/}2}{e}^{-u}\frac{\partial }{\partial \psi }\left(w{e}^u\right)\text{sin}\left(\theta \right)d\theta \approx \\ \approx {\int }_0^{\pi \text{/}2}{e}^{-u}\frac{\partial }{\partial \psi }\left({w_1}^0{e}^{u_1}\right)\text{sin}\left(\theta \right)d\theta =\\ ={\int }_0^{\pi \text{/}2}{e}^{-u}\text{sin}\left(\theta \right)d\theta \cdot \frac{\partial }{\partial \psi }\left({w_1}^0{e}^{u_1}\right) .\end{array}$ (28)

Далее:

$\begin{array}{c}{\int }_0^{\pi \text{/}2}{e}^{-u}\text{sin}\left(\theta \right)d\theta =e^{-\frac{1}{2}g{\psi }^2}\times \\ \times {\int }_0^{\pi \text{/}2}{e}^{\sigma {\text{cos}}^2\left(\theta \right)}{e}^{h\left({\zeta }_{}\cos \left(\theta \right)+{\zeta }_{}\sin \left(\theta \right)\text{cos}\varphi \right) }\text{sin}\left(\theta \right)d\theta .\end{array}$ (29)

Снова, как и в (21), (22), используя соображения метода перевала, получаем

$\begin{array}{c}{\int }_0^{\pi \text{/}2}{e}^{\sigma {\text{cos}}^2\left(\theta \right)}{e}^{h\left({\zeta }_{}\cos \right(\theta )+{\zeta }_{}\sin (\theta \left)\text{cos}\varphi \right) }\text{sin}\theta d\theta \approx \\ \approx e^{\sigma +h{\zeta }_{}}{\int }_0^{\pi \text{/}2}{e}^{-\sigma {\theta }^2}\theta d\theta =\\ \frac{1}{2\sigma }{e}^{\sigma +h{\zeta }_{}}(1-e^{-\sigma {\left(\pi \text{/}2\right)}^2})\approx \frac{1}{2\sigma }{e}^{\sigma +h{\zeta }_{}}=\frac{1}{2\sigma }{e}^{-u_{\mu 1}}.\end{array}$ (30)

Следовательно, интеграл в (25) может быть оценен следующим образом:

$\begin{array}{c}{\int }_0^{2\pi }{\int }_0^{\pi /2}{\text{div}}_{\nu }{j}_{\nu }\text{sin}\left(\theta \right)d\theta d\varphi \approx \\ \approx \frac{\pi }{\sigma } \frac{1}{\text{sin}\psi }\frac{\partial }{\partial \psi }\left(\text{sin}\psi e^{-u_{\mu 1}}\frac{\partial }{\partial \psi }\left({w_1}^0{e}^{u_{\mu 1}}\right)\right).\end{array}$ (31)

Уравнение (23) сейчас может быть представлено в форме:

$\begin{array}{c}\frac{\partial {w_1}^0}{\partial t }=\frac{1}{2{\tau }_{\text{II}}}\left[\right({w_2}^0{e}^{h{\zeta }_{}}-{w_1}^0{e}^{-h{\zeta }_{}}] -\\ -\frac{1}{2{\tau }_{\eta }}\frac{1}{\text{sin}\psi }\frac{\partial }{\partial \psi }\left(\text{sin}\psi e^{-u_{\mu 1}}\frac{\partial }{\partial \psi }\left({w_1}^0{e}^{u_{\mu 1}}\right)\right) . \end{array}$ (32)

Аналогично, уравнение (24) может быть переписано в виде:

$\begin{array}{c}\frac{\partial {w_2}^0}{\partial t }=-\frac{1}{2{\tau }_{\text{II}}}\left[\right({w_2}^0{e}^{h{\zeta }_{}}-{w_1}^0{e}^{-h{\zeta }_{}}] -\\ -\frac{1}{2{\tau }_{\eta }}\frac{1}{\text{sin}\psi }\frac{\partial }{\partial \psi }\left(\text{sin}\psi e^{-u_{\mu 2}}\frac{\partial }{\partial \psi }\left({w_2}^0{e}^{u_{\mu 2}}\right)\right).\end{array}$ (33)

В рамах используемого приближения σ >> 1, σ >> h, в областях 1 и 2, проиллюстрированных на рис. 2, вектор µ магнитного момента частиц, в основном, направлен вдоль и противоположно вектору ν оси легкого намагничивания частицы. Следовательно, при заданном угле β между H и ν₀ (рис. 1), статистически средняя проекция вектора µ на направление поля H может быть представлена в виде:

$\begin{array}{c}〈{\mathrm{\mu }}_{\mathrm{h}}〉\left(\mathrm{\beta }\right)=〈{\mathrm{\mu }}_1〉\left(\mathrm{\beta }\right)+〈{\mathrm{\mu }}_2〉\left(\mathrm{\beta }\right),\\ 〈{\mathrm{\mu }}_1〉\left(\mathrm{\beta }\right)={\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }}{\mathrm{\mu }}_{\text{h}}{\mathrm{w}}_1^0\text{sin}\mathrm{\psi }d\mathrm{\psi }d\mathrm{\phi },\\ 〈{\mathrm{\mu }}_2〉\left(\mathrm{\beta }\right)=-{\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }}{\mathrm{\mu }}_{\text{h}}{\mathrm{w}}_2^0\text{sin}\mathrm{\psi }d\mathrm{\psi }d\mathrm{\phi }.\end{array}$ (34)

В рамках обсужденного приближения выполняется приближенное равенство $\left(\mathbf{\mu }\mathbf{\cdot }\mathbf{H}\right)\approx \left(\mathbf{\nu }\mathbf{\cdot }\mathbf{H}\right).$ Следовательно, можно записать:

${\mathrm{\mu }}_{\text{h}}={\mathrm{\zeta }}_{}=\text{cos}\mathrm{\beta }\text{cos}\mathrm{\psi }+\text{sin}\mathrm{\beta }\text{sin}\mathrm{\psi }\text{cos}\mathrm{\phi }.$ (35)

Наша цель сейчас определить функции w⁰_1,2. В общем случае решение кинетических уравнений (32), (33) неизвестно. Чтобы найти приближенные решения, введем вероятности

${\mathrm{W}}_{\mathrm{1,2}}\left(\mathrm{t}\text{|}\mathrm{\beta }\right)={\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }}{\mathrm{w}}_{\mathrm{1,2}}^0\text{sin}\mathrm{\psi d\psi d\phi }$ (36)

того, что магнитный момент частицы находится соответственно в областях “1” и “2” на рис. 2. Учитывая, что магнитная анизотропия частицы предполагается сильной, т. е. то, что неравенства σ >> 1, σ >> h выполняются, справедливо сильное неравенство τ_II >> τ_h. Это означает, что релаксация ориентации вектора ν, связанная с поворотом частицы в окружающей среде, происходит намного быстрее, чем неелевский переход момента частицы через потенциальный барьер ее внутренней магнитной анизотропии из области “1” в область “2”.

В рамках этого приближения можно записать

$\begin{array}{c}{\mathrm{w}}_{\mathrm{1,2}}^0\left(\mathrm{\psi },\mathrm{t}\text{|}\mathrm{\beta }\right)\approx {\mathrm{W}}_{\mathrm{1,2}}\left(\mathrm{t}\text{|}\mathrm{\beta }\right){\mathrm{f}}_{\mathrm{1,2}}\left(\mathrm{\psi },\mathrm{t}\text{|}\mathrm{\beta }\right),\\ {\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }}{\mathrm{f}}_{\mathrm{1,2}}\text{sin}\mathrm{\psi d\psi d\phi }=\mathrm{1,}\end{array}$ (37)

где f_1,2(Ψ,t |φ) – плотность вероятности отклонения вектора ν от ν₀ на угол Ψ (рис. 1), если магнитный момент частицы находится в области 1 или 2 соответственно.

Подставляя (37) в (32) и (33), после преобразований получаем:

$\begin{array}{c}{W}_{\mathrm{1,2}}\frac{\partial f_{\mathrm{1,2}}}{\partial t}+f_{\mathrm{1,2}}\frac{\partial W_{\mathrm{1,2}}}{\partial t}=\\ =\frac{1}{2{\tau }_{\text{II}}}\left(W_{\mathrm{2,1}}\right(t\text{|}\phi )f_{\mathrm{2,1}}{e}^{h{\mu }_{\text{h}}}-W_{\mathrm{1,2}}(t\text{|}\phi \left)f_{\mathrm{1,2}}{e}^{-h{\mu }_{\text{h}}}\right)+\\ +\frac{W_{\mathrm{1,2}}}{2{\tau }_{\eta }}\left[\begin{array}{c}\frac{1}{\text{sin}\psi }\frac{\partial }{\partial \psi }\left(\text{sin}\psi e^{-u_{\mu \mathrm{1,2}}}\frac{\partial }{\partial \psi }\left(f_{\mathrm{2,1}}{e}^{u_{\mu \mathrm{1,2}}}\right)\right)+\\ +\frac{1}{{\text{sin}}^2\psi }\frac{\partial }{\partial \beta }\left(e^{-u_{\mu \mathrm{1,2}}}\frac{\partial }{\partial \beta }\left(f_{\mathrm{2,1}}{e}^{u_{\mu \mathrm{1,2}}}\right)\right)\end{array}\right].\end{array}$ (38)

Поскольку переход момента частицы между потенциальными областями “1” и “2” является относительно редким событием (τ_II >> τ_h), мы можем искать функции f_1,2 в квазистационаром приближении, предполагая, что за характерное время изменения f_1,2 вероятности W_1,2 остаются практически неизменными. В этом приближении из уравнения (38) приходим к уравнению

$\begin{array}{c}\frac{\partial f_{\mathrm{1,2}}}{\partial t}=\frac{1}{2{\tau }_{\eta }}\times \\ \times \left[\begin{array}{c}\frac{1}{\text{sin}\psi }\frac{\partial }{\partial \psi }\left(\text{sin}\psi e^{-u_{\mathrm{1,2}}}\frac{\partial }{\partial \psi }\left(f_{\mathrm{1,2}}{e}^{u_{\mathrm{1,2}}}\right)\right)+\\ +\frac{1}{{\text{sin}}^2\psi }\frac{\partial }{\partial \beta }\left( e^{-u_{\mathrm{1,2}}}\frac{\partial }{\partial \beta }\left(f_{\mathrm{1,2}}{e}^{u_{\mathrm{1,2}}}\right)\right)\end{array}\right].\end{array}$ (39)

Отметим, что в случае отсутствия упругих сил в среде (g = 0), уравнение (39) совпадает с уравнением Фоккера–Планка для ферромагнитной частицы в жидкой среде.

В течение относительно длительного времени между двумя переходами момента частицы из области “1” в область “2” угол Ψ успевает принять практически все значения в промежутке от 0 до p, а азимутальный угол j – все значения от 0 до 2p. Учитывая это, усредним (38) по всем Ψ и j, соответствующих областям 1 и 2. Учитывая, что функция f_1,2 нормализована на единицу, приходим к следующему соотношению

$\begin{array}{c}\frac{\partial W_1}{\partial t}=\left(W_2\right(t\text{|}\beta )G_2-W_1(t\text{|}\beta \left)G_1\right),\\ \frac{\partial W_2}{\partial t}=-\left(W_2\right(t\text{|}\beta )G_2-W_1(t\text{|}\beta \left)G_1\right),\\ G_1=\frac{1}{2}{\int }_0^{2\pi }{\int }_0^{\pi }\frac{f_1}{{\tau }_{II}}{e}^{-h{\mu }_{\text{h}}}\text{sin}\psi d\psi d\phi ,\\ G_2=\frac{1}{2}{\int }_0^{2\pi }{\int }_0^{\pi }\frac{f_2}{{\tau }_{II}}{e}^{h{\mu }_{\text{h}}}\text{sin}\psi d\psi d\phi .\end{array}$ (40)

Напомним, что величина τ_II определена в (23).

По основному свойству вероятности должно выполняться равенство W₁ + W₂ = 1. Следовательно, уравнение (40) может быть представлено в виде

$\frac{\partial {\mathrm{W}}_1\left(\mathrm{t}\text{|}\mathrm{\beta }\right)}{\partial \mathrm{t}}=[{\mathrm{G}}_2-{\mathrm{W}}_1(\mathrm{t}\text{|}\mathrm{\beta }\left)\right({\mathrm{G}}_2+{\mathrm{G}}_1\left)\right].$ (41)

Дальнейшая программа действий должна быть такой. Во-первых, необходимо решить уравнение (39) для функций f_1,2. Затем – решить релаксационное уравнение (41) для W₁. После этого, используя соотношения (34)–(37), можно найти статистически среднюю проекцию магнитного момента ⟨µ_h⟩.

Строгое аналитическое решение уравнения (39) неизвестно. Чтобы найти приближенное решение, мы воспользуемся методом “эфективного поля”, предложенного в [16]. Этот метод, являющийся вариантом общего метода пробных функций, был многократно успешно применен для описания кинетики перемагничивания и магнитовязких эффектов в жидких магнитных коллоидах (см. обсуждение, например, в [27–30]). Ключевая идея этого метода состоит в том, чтобы представить неравновесные функции f_1,2 в виде равновесных функций распределения, однако не в реальном магнитном поле H, а в некотором эффективном H_ef, компоненты которого необходимо определить. Эти компоненты определяются из уравнений для первых статистичских моментов функции f_1,2.

Равновесная функция f₁ в реальном поле H может быть представлена в больцмановской форме

$\begin{array}{c}{f_1}^{\text{eq}}\left(\psi \right|\phi )=\frac{1}{Z}\text{exp(}-u_{\mu 1}\text{)}=\frac{1}{Z}\text{exp}\left(h{\mu }_{\text{h}}-\frac{1}{2}g{\psi }^2\right),\\ Z={\int }_0^{\pi }\text{exp}\left(h{\mu }_{\text{h}}-\frac{1}{2}g{\psi }^2\right)\text{sin}\psi d\psi .\end{array}$

Следуя идеям [16], мы представим искомую неравновесную функцию f₁ в виде

$\begin{array}{c}{f}_1\left(\psi \right|\beta )=\frac{1}{Z_1^{\text{ef}}}\exp \left(h_1^{\text{ef}}{\mu }_{\text{h}}-\frac{1}{2}g{\psi }^2\right),\\ Z_1^{\text{ef}}={\int }_0^{\pi }\exp \left(h_1^{\text{ef}}{\mu }_{\text{h}}-\frac{1}{2}g{\psi }^2\right)\text{sin}\psi d\psi .\end{array}$ (42)

Аналогично для f₂ используем пробную функцию

$\begin{array}{c}{f}_2\left(\psi \right|\beta )=\frac{1}{Z_2^{\text{ef}}}\exp \left(-h_2^{\text{ef}}{\mu }_{\text{h}}-\frac{1}{2}g{\psi }^2\right),\\ Z_2^{\text{ef}}={\int }_0^{\pi }\exp \left(h_2^{\text{ef}}{\mu }_{\text{h}}-\frac{1}{2}g{\psi }^2\right)\text{sin}\psi d\psi .\end{array}$ (43)

3. УРАВНЕНИЯ КИНЕТИКИ ПЕРЕМАГНИЧИВАНИЯ ЧАСТИЦЫ

Умножая обе части уравнения (39) на µ_h и интегрируя их по правилам сферической системы координат, после несложных преобразований приходим к следующим уравнениям:

$\begin{array}{c}2{\tau }_{\eta }\frac{d}{dt}{\left\{{\mu }_{\mathrm{h}}\right\}}_1=-{\left\{{\mu }_{\mathrm{h}}\right\}}_{1}1+\\ +\cdots {\int }_0^{2\pi }{\int }_0^{\pi }\frac{\partial {\mu }_{\text{h}}}{\partial \psi }\cos \psi f_{1}d\psi d\phi +\\ +h{\left\{{\left(\frac{\partial {\mu }_{\text{h}}}{\partial \psi }\right)}^2\right\}}_1-g{\left\{\psi \frac{\partial {\mu }_{\text{h}}}{\partial \psi }\right\}}_1\end{array}$ (44)

и

$\begin{array}{c}2{\tau }_{\eta }\frac{d}{dt}{\left\{{\mu }_{\mathrm{h}}\right\}}_2=-{\left\{{\mu }_{\mathrm{h}}\right\}}_2+\\ +{\int }_0^{2\pi }{\int }_0^{\pi }\frac{\partial {\mu }_{\text{h}}}{\partial \psi }\cos \psi f_{2}d\psi d\phi -\\ -h{\left\{{\left(\frac{\partial {\mu }_{\text{h}}}{\partial \psi }\right)}^2\right\}}_2-g{\left\{\psi \frac{\partial {\mu }_{\text{h}}}{\partial \psi }\right\}}_2,\end{array}$ (45)

где

${\left\{\dots \right\}}_{\mathrm{1,2}}={\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }}\dots {\mathrm{f}}_{\mathrm{1,2}}\mathrm{sin\psi }d\mathrm{\psi }d\mathrm{\phi }.$ (46)

Учитывая (35), получаем $\frac{\partial {\mathrm{\mu }}_{\text{h}}}{\partial \mathrm{\psi }}=-\text{cos}\mathrm{\beta }\text{sin}\mathrm{\psi }+\text{sin}\mathrm{\beta }\text{cos}\mathrm{\psi }\text{cos}\mathrm{\phi }.$

Используя здесь функции f_1,2 в виде (42) и (43), а также соотношение (35), мы выражаем {µ_h}_1,2 через $h_{\mathrm{1,2}}^{\text{ef}}$ и приходим к замкнутым уравнениям относительно $h_{\mathrm{1,2}}^{\text{ef}}$. Решив эти уравнения, находим функции распределения f₁ и f₂ в форме (42) и (43). Далее, используя уравнения (40), (41), определяем плотности вероятностей W₁ и W₂ = 1 – W₁. Подставляя полученные результаты в (34), (37), определяем искомую среднюю проекцию < µ_h > (β) вектора момента частицы µ на приложенное поле H.

Отметим, что уравнения (44), (45), вообще говоря, не могут быть записаны в форме классических уравнений Дебая. Это значит, что при изменении поля релаксация намагниченности частиц не происходит по простому экспоненциальному закону с определенным временем релаксации, как это имеет место, например, в случае иммобилизованных частиц [15, 17, 18], а также частиц, вращающихся вместе со своим магнитным моментов в вязкой среде [16, 27].

В общем случае уравнения для $h_1^{\text{ef}}$ и $h_2^{\text{ef}}$ нелинейны и не могут быть решены аналитически. Однако они могут быть решены численно. Далее мы рассмотрим подробно случай слабых внешних полей, когда выполняется неравенство h < 1.

4. ДИНАМИЧЕСКАЯ ВОСПРИИМЧИВОСТЬ В СЛАБЫХ ПОЛЯХ

В этой части мы рассмотрим слабые магнитные поля, когда выполняются сильные неравенства h_1,2 << 1 и, следовательно, $h_{\mathrm{1,2}}^{\text{ef}}$ << 1.

В линейном приближении по отношению h_1,2 и $h_{\mathrm{1,2}}^{\text{ef}}$ соотношения (38), (39) могут быть переписаны в виде

$\begin{array}{c}{f}_1\left(\psi \right|\beta )=\frac{1}{Z_0}\text{exp}\left(-\frac{1}{2}g{\psi }^2\right)\times \\ \times (1+h_1^{\text{ef}}({\mu }_{\text{h}}-\text{cos}\beta {〈\text{cos}\psi 〉}^0\left)\right)\end{array}$ (47)

и

$\begin{array}{c}{f}_2\left(\psi \right|\beta )=\frac{1}{Z_0}\text{exp}\left(-\frac{1}{2}g{\psi }^2\right)\times \\ \times (1-h_2^{\text{ef}}({\mu }_{\text{h}}-\text{cos}\beta {〈\text{cos}\psi 〉}^0\left)\right)\end{array}$ (48)

соответственно, где снова ${\mathrm{\mu }}_{\mathrm{h}}=\text{cos}\mathrm{\beta }\text{cos}\mathrm{\psi }+\text{sin}\mathrm{\beta }\text{sin}\mathrm{\psi }\text{cos}\mathrm{\phi }$ и

$\begin{array}{c}{〈\dots 〉}^0=\frac{1}{Z^0}{\int }_0^{2\pi }{\int }_0^{\pi }\dots \exp \left(-\frac{1}{2}g{\psi }^2\right)\sin \psi d\psi d\phi ,\\ Z^0={\int }_0^{2\pi }{\int }_0^{\pi }\exp \left(-\frac{1}{2}g{\psi }^2\right)\sin \psi d\psi d\phi .\end{array}$ (49)

Используя соотношения (48) и (49) в (44) и (45), после простых, но громоздких преобразований, приходим к одному и тому же линейному дифференциальному уравнению относительно $h_1^{\text{ef}}$ и $h_2^{\text{ef}}$

$\begin{array}{c}2{\tau }_{\eta }A\left(\beta ,g\right)\frac{d}{dt}{h}^{\text{ef}}=h-h^{\text{ef}},\\ A=\frac{2{\text{cos}}^2\beta \left({〈{\text{cos}}^2\psi 〉}^0-{\left({〈\text{cos}\psi 〉}^0\right)}^2\right)+{\text{sin}}^2\beta {〈{\text{sin}}^2\psi 〉}^0}{2{\text{cos}}^2\beta {〈si{n}^2\psi 〉}^0+{\text{sin}}^2\beta \left(1+{〈{\text{cos}}^2\psi 〉}^0\right)}.\end{array}$ (50)

Здесь и далее

$h_1^{\text{ef}}=h_2^{\text{ef}}=h^{\text{ef}}.$

Пусть поле Н по закону $H=H_0\cos \omega t.$ Безразмерное магнитное поле запишем в виде

$h\left(t\right)=h_0\text{cos}\omega t,$ (51)

где $h_0={\mu }_0\frac{\mu }{kT}{H}_0$ – безразмерная амплитуда поля.

Решение уравнения (50) может быть представлено в виде:

$\begin{array}{c}{h}^{\text{ef}}=({\kappa }^{\text{'}}\text{cos}\omega t+{\kappa }^{\text{'}\text{'}}sin\omega t)h_0,\\ {\kappa }^{\text{'}}=\frac{1}{1+{\left(2\omega A{\tau }_{\eta }\right)}^2}; {\kappa }^{\text{'}\text{'}}=\frac{2\omega A{\tau }_{\eta }}{1+{\left(2\omega A{\tau }_{\eta }\right)}^2}.\end{array}$ (52)

Здесь 2Aτ_h является временем вязкой релаксации магнитного момента частицы при ее вращении в окружающей среде. Отметим, что в чисто жидкой среде это время равно τ_h. Поэтому коэффициент А, зависящий от модуля g, характеризует влияние упругости несущей среды на этот релаксационный процесс. Некоторые результаты расчета А приведены на рис. 3.

 

Рис. 3. Зависимость множителя A в (48.50) от угла b, показанного на рис. 1 для трех безразмерных модулей упругости несущей среды. Кривые 1, 2, 3: g = 0, 0.5 и 1.5 соответственно.

 

Отметим, что в отсутствие упругих эффектов (g = 0) параметр A равен ½ и не зависит от угла β между H и ν₀, как это и должно быть из физических соображений. Кривые на рис. 3 симметричны относительно β = π/2.

На рис. 3 показано, что увеличение безразмерного модуля упругости g уменьшает время 2Aτ_h вязкой релаксации. Отметим, что для частицы с диаметром d_p ~ 20 нм при комнатных температурах значения g = 0.5 и 1.5 соответствуют значениям модуля упругости несущей среды G ~ 0.5 и 1.5 кПа соответственно. Это вполне соответствует мягким гелям.

Подставляя $h_1^{\text{ef}}=h_2^{\text{ef}}=h^{\text{ef}}$ из (52) в (47) и (48), получаем явные выражения для неравновесных функций распределения f₁ и f₂.

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

$\begin{array}{c}\frac{\partial W_1}{\partial t}=\frac{1}{2{\tau }_{\text{N}}}\left[\begin{array}{c}1+h\text{cos}\left(\beta \right)\left({〈\cos \left(\psi \right)〉}^0+{〈\text{sin}\left(\psi \right)〉}^0\right)-\\ -2\left(1+h\text{cos}\left(\beta \right){〈\sin \left(\psi \right)〉}^0\right)W_1\end{array}\right],\\ {\tau }_{\text{N}}={\tau }_{\text{D}}\sqrt{\frac{\pi }{{\sigma }^3}}\exp \left(\sigma \right).\end{array}$(53)

Здесь τ_ν – характерное время неелевского перемагничивания неподвижной частицы при исчезающе малых магнитных полях [15].

Уравнение (53) может быть решено аналитически. В линейном приближении по h = h₀cos(ωt) это решение имеет вид:

$\begin{array}{c}{W}_1=\frac{1}{2}\left(1+h_0{\chi }_{\text{N}}^0\left(\begin{array}{c}\frac{1}{1+{\left(\omega {\tau }_{\text{N}}\right)}^2}\cos \omega t+\\ +\omega {\tau }_{\text{N}}\frac{\omega {\tau }_{\text{N}}}{1+{\left(\omega {\tau }_{\text{N}}\right)}^2}\sin \omega t\end{array}\right)\right),\\ {\chi }_{\text{N}}^0=\left({〈\cos \left(\psi \right)〉}^0+2{〈sin\left(\psi \right)〉}^0\right){\cos }^2\left(\phi \right){〈\cos \left(\psi \right)〉}^0.\end{array}$ (54)

Статистически средняя компонента <µ_h> единичного вектора µ (см. уравнение (34)) может быть представлена в виде:

$\begin{array}{c}〈{\mathrm{\mu }}_{\text{h}}〉\left(\mathrm{\beta }\right)=<{\mathrm{\mu }}_{\text{h1}}>\left(\mathrm{\beta }\right)+<{\mathrm{\mu }}_{\text{h2}}>\left(\mathrm{\beta }\right),\\ 〈{\mathrm{\mu }}_{\text{h1}}〉={\mathrm{W}}_1\left(\mathrm{t}\right){\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }}{\mathrm{\mu }}_{\text{h}}{\mathrm{f}}_1\left(\mathrm{\psi },\mathrm{\phi }\text{|}\mathrm{\beta }\right)\text{sin}\mathrm{\psi d\psi d\phi },\\ 〈{\mathrm{\mu }}_{\text{h2}}〉=-{\mathrm{W}}_2\left(\mathrm{t}\right){\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }}{\mathrm{\mu }}_{\text{h}}{\mathrm{f}}_2\left(\mathrm{\psi },\mathrm{\phi }\text{|}\mathrm{\beta }\right)\text{sin}\mathrm{\psi d\psi d\phi }=\\ =-\left(1-{\mathrm{W}}_1\left(\mathrm{t}\right)\right){\int }_0^{\mathrm{\pi }}{\mathrm{\mu }}_{\text{h}}{\mathrm{f}}_2\left(\mathrm{\psi }\text{|}\mathrm{\phi }\right)\text{sin}\mathrm{\psi d\psi }=\\ =-\left(1-{\mathrm{W}}_1\left(\mathrm{t}\right)\right){\int }_0^{2\mathrm{\pi }}{\int }_0^{\mathrm{\pi }}{\mathrm{\mu }}_{\text{h}}{\mathrm{f}}_2\left(\mathrm{\psi },\mathrm{\phi }\text{|}\mathrm{\beta }\right)\text{sin}\mathrm{\psi d\psi d\phi }.\end{array}$ (55)

Используя уравнение (43), получаем:

${\mu }_1=W_1\left(t\right)\left[\begin{array}{c}\cos \left(\beta \right){〈\cos \left(\psi \right)〉}^0+h^{\text{ef}}\times \\ \times \left(\begin{array}{c}{\cos }^2\left(\beta \right){〈co{s}^2\left(\psi \right)〉}^0+\\ +\frac{1}{2}{\sin }^2\left(\beta \right){〈si{n}^2\left(\psi \right)〉}^0-\\ -{\left(\cos \left(\beta \right){〈cos\left(\psi \right)〉}^0\right)}^2\end{array}\right)\end{array}\right].$ (56)

Аналогично

$\begin{array}{c}{\mu }_2=-(1-W_1\left(t\right))\times \\ \times \left[\begin{array}{c}\cos \left(\beta \right){〈\cos \left(\psi \right)〉}^0-\\ -h^{\text{ef}}\left(\begin{array}{c}{\cos }^2\left(\beta \right){〈{\text{cos}}^2\left(\psi \right)〉}^0+\\ +\frac{1}{2}{\sin }^2\left(\beta \right){〈{\text{sin}}^2\left(\psi \right)〉}^0{〈{\text{sin}}^2\left(\psi \right)〉}^0-\\ \text{-}{\left(\cos \left(\beta \right){〈cos\left(\psi \right)〉}^0\right)}^2-\end{array}\right)\end{array}\right].\end{array}$ (57)

Используя уравнения (53)–(57), приходим к соотношению

$〈{\mathrm{\mu }}_{\text{h}}〉\left(\mathrm{\beta }\right)={\mathrm{h}}_0\left({\mathrm{\chi }}^{\text{'}}\right(\mathrm{\beta })\mathrm{cos\omega t}+{\mathrm{\chi }}^{\text{'}\text{'}}(\mathrm{\beta }\left)\mathrm{sin\omega t}\right),$ (58)

где

$\begin{array}{c}{\chi }^{\text{'}}\left(\beta \right)={\chi }_{\text{N}}^0\frac{1}{1+{\left(\omega {\tau }_{\text{N}}\right)}^2}+{\chi }_{\eta }^0\frac{1}{1+{\left(2\omega A{\tau }_{\eta }\right)}^2},\\ {\chi }^{\text{'}\text{'}}\left(\beta \right)=\omega \left(\begin{array}{c}{\chi }_{\text{N}}^0\frac{1}{1+{\left(\omega {\tau }_{\text{N}}\right)}^2}{\tau }_{\text{N}}+\\ +{\chi }_{\eta }^0\frac{1}{1+{\left(2\omega A{\tau }_{\eta }\right)}^2}2A{\tau }_{\eta }\end{array}\right),\\ {\chi }_{\text{N}}^0=\left({〈\cos \left(\psi \right)〉}^0+2{〈\text{sin}\psi 〉}^0\right){\cos }^2\left(\beta \right){〈\cos \left(\psi \right)〉}^0,\\ {\chi }_{\eta }^0={\cos }^2\left(\beta \right)\left({〈{\cos }^2\left(\psi \right)〉}^0-{\left({〈\cos \left(\psi \right)〉}^0\right)}^2\right).\end{array}$ (59)

Параметры χʹ(β) и χʺ(β) могут рассматриваться как реальная и мнимая части комплексной восприимчивости, связывающей ⟨µ_h⟩ и h при заданном угле β между полем H и вектором ν₀ (рис. 1). Отметим, что в предельных случаях, когда первые или вторые слагаемые доминируют в правых частях (59), эти соотношения соответствуют дебаевскому закону релаксации магнитного момента по неелевскому или вязкому механизму соответственно. В общем случае релаксация магнитного момента имеет существенно недебаевский характер.

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

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

$\overline{{\mathrm{\chi }}^{\text{'}}}=\frac{1}{2}{\int }_0^{\mathrm{\pi }}{\mathrm{\chi }}^{\text{'}}\left(\mathrm{\beta }\right)\mathrm{sin\beta d\beta },\overline{{\mathrm{\chi }}^{\text{'}\text{'}}}=\frac{1}{2}{\int }_0^{\mathrm{\pi }}{\mathrm{\chi }}^{\text{'}\text{'}}\left(\mathrm{\beta }\right)\mathrm{sin\beta d\beta }.$ (60)

Некоторые результаты расчетов χʹ(β) и χʺ(β) для β = 0 и π/2, а также их усредненные по всем β значения (60) приведены на рис. 4.

 

Рис. 4. Реальная и мнимая части комплексной восприимчивости при t_N = 10 τ_h (a) и (б): угол β равен 0 и π/2 соответственно; (в) – усредненные по β компоненты $\overline{\mathbf{\chi }}$, определенные в (60). Кривые 1, 2, 3 – безразмерный коэффициент упругости g = 0.1, 0.5, 1.5 соответственно.

 

При β = 0 (параллельная ориентация ν₀ и H, рис. 4a), а также при усреднении по всем возможным β (рис. 4в), увеличение безразмерного модуля g ведет к увеличению реальных восприимчивостей χʹ(0) и $\overline{\mathrm{\chi }}$ʹ при малых частотах ω и к их уменьшению при достаточно больших ω. Для этих ориентаций ν₀ относительно H при малых g мнимые восприимчивости χʺ(0) и $\overline{\mathrm{\chi }}$ʺ имеют единственный максимум, соответствующий вязкой релаксации частицы. При достаточно больших g каждая из этих мнимых восприимчивостей имеет один максимум, соответствующий неелевской релаксации магнитного момента. При промежуточных значениях g мнимые части восприимчивости имеют по два максимума. Таким образом, релаксация намагниченности, даже при слабых полях, имеет существенно недебаевский характер.

Напротив, если векторы ν₀ и H взаимно перпендикулярны (рис. 4б), увеличение g уменьшает реальную часть ${\mathrm{\chi }}^{\text{'}}\left(\frac{\mathrm{\pi }}{2}\right)$ восприимчивости для малых ω и увеличивает для больших частот. Мнимая часть ${\mathrm{\chi }}^{\text{'}\text{'}}\left(\frac{\mathrm{\pi }}{2}\right)$ имеет только один максимум, соответствующий вязкой релаксации для всех исследованных значений g. Отсутствие второго максимума, соответствующего неелевской релаксации, объясняется выбранным приближением сильной магнитной анизотропией частицы.

Оценим теперь параметры системы, соответствующие выбранному при расчетах соотношению τ_N = 10 τ_η. Для магнетитовых частиц, наиболее часто используемых в приложениях, параметр объемной магнитной анизотропии $\mathrm{K}\approx \frac{14\text{ кДж}}{{\text{м}}^3},$ а время τ_D ~ 10^–9 с (см., например, [22]). Простые расчеты показывают, что для частиц диаметром 20 нм при учете только ее объемной магнитной анизотропии безразмерный параметр магнитной анизотропии σ ≈ 14 (т. е. условие сильной анизотропии σ >> 1 выполняется) и соотношение τ_N = 10 τ_h соответствуют вполне реалистическому значению вязкости несущей среды η ≈ 3⋅10^–3 Па⋅с. Если учитывать вклад поверхностных эффектов в магнитную анизотропию частицы, то значение полного параметра σ может существенно увеличиться, а время τ_N возрасти на порядки величины. Поэтому выбранное условие τ_N = 10 τ_h будет выполняться и для систем с намного более высокой вязкостью.

МАГНИТНАЯ ГИПЕРТЕРМИЯ

Переменное магнитное поле вызывает нагрев магнитных частиц вследствие диссипации магнитной энергии при их перемагничивании – так называемый эффект магнитной гипертермии [27]. Этот эффект вызывает значительный интерес исследователей, биоинженеров и врачей клиницистов, поскольку он является перспективным для лечения опухолевых, включая онкологические, заболеваний [7–10, 23, 24, 27, 31–36]. Насколько нам известно, в литературе теоретические модели магнитной гипертермии развиты или применительно к жидким суспензиям магнитных наночастиц, или для жестких сред, в которых магнитные частицы полностью обездвижены. Однако, как уже обсуждалось, в мягких биологических средах перемагничивание частиц может происходить комбинированно как их частичное вращение в упруго-вязкой среде и по сценарию неелевского перехода через барьер магнитной анизотропии.

Используя общий термодинамический подход [27], в случае слабых полей h << 1 интенсивность генерации тепла (диссипации магнитной энергии) в частице можно представить в виде

$p=kT\frac{h_0^2}{{\tau }_{\eta }} \xi \left(\beta \right),\xi \left(\beta \right)=\omega {\tau }_{\eta }{\chi }^{\text{'}\text{'}}\left(\beta \right)$ (61)

при фиксированном угле β, или

$p=kT\frac{h_0^2}{{\tau }_{\eta }}\overline{\xi } ,\overline{\xi }=\omega {\tau }_{\eta }\overline{{\chi }^{\text{'}\text{'}}}$ (62)

для среды с хаотической начальной ориентацией осей легкого намагничивания частиц. Отметим, что для случаев предельно мягких (g = 0) и жестких (g → ∞) систем соотношения (61), (62) для интенсивности тепловыделения совпадают с таковыми в [27] для частиц, вращающихся в вязких жидкостях, и иммобилизованных частиц в жестких средах соответственно.

Некоторые результаты расчетов безразмерных параметров ξ(β) и $\overline{\mathrm{\xi }}$ показаны на рис. 5. Для β = 0 и для хаотической ориентации осей частиц увеличение коэффициента жесткости g несущей среды ведет к увеличению безразмерной генерации тепла ξ для малых частот w и к уменьшению ξ при достаточно высоких частотах. При β = π/2 увеличение g уменьшает значение ξ.

 

Рис. 5. Безразмерные параметры ξ(β) и $\overline{\mathbf{\xi }}$ генерации тепла частицей как функции от частоты поля. Значки а, б, в и 1, 2, 3 означают то же, что на рис. 4.

 

Форма зависимости параметров ξ(β) и $\overline{\mathrm{\xi }}$, следовательно, интенсивности тепловыделения р от частоты поля w определяется тем, что при низких значениях g (очень мягкий несущий гель) доминирующий вклад в генерирование тепла оказывает механизм вращения частицы в окружающей упруго-вязкой среде. При увеличении g (жесткость геля растет) относительный вклад неелевского механизма увеличивается. Вначале он доминирует только при малых частотах; при достаточно больших g неелевский механизм доминирует при всех исследованных нами частотах. Во всех случаях при увеличении w параметр ξ стремится к насыщению. Некоторые результаты расчетов предельных (ω → ∞) значений ξ(β) и $\overline{\mathrm{\xi }}$ показаны на рис. 6.

 

Рис. 6. Предельные (ω → ∞) значения ξ(β) и $\overline{\mathbf{\xi }}$ как функции от безразмерного модуля жесткости композита g. Кривые 1 и 2 – значения ξ(β) при β = 0 и π/2 соответственно; кривая 3 – усредненное значение $\overline{\mathbf{\xi }}$.

 

ЗАКЛЮЧЕНИЕ

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

Предполагается сильная магнитная анизотропия частицы (σ >> 1). Неелевская релаксация магнитного момента частицы исследована на основе математически регулярного метода Крамерса-Брауна диффузии магнитного момента частицы через потенциальный барьер магнитной анизотропии. Вращательная вязкая релаксация момента частицы – при помощи метода эффективного поля, хорошо описывающего релаксационные процессы в коллоидных магнитных жидкостях. Модель развита для широкого диапазона напряженности приложенного магнитного поля Н, ограниченного только условием, что энергия взаимодействия частицы с полем меньше энергии ее магнитной анизотропии (в принятых обозначениях h < σ).

В результате получены уравнения, позволяющие определить статистически среднюю проекцию ⟨µ_h⟩ единичного вектора µ магнитного момента частицы на поле Н. Это дает возможность рассчитывать макроскопическую намагниченность композита, пропорциональную ⟨µ_h⟩.

Более подробно рассмотрен случай слабых полей, когда справедлива линейная зависимость ⟨µ_h⟩ от Н. Редуцированную комплексную магнитную восприимчивость χ мы определяем как коэффициент пропорциональности между ⟨µ_h⟩ и безразмерным магнитным полем h. В рамках принятых приближений найдены зависимости реальной χʹ и мнимой χʺ частей χ от частоты поля, угла β между полем H и начальной ориентацией оси легкого намагничивания частицы. Для β = 0 (параллельная ориентация этих векторов) и при хаотической начальной ориентации этих осей частиц в очень мягких материалах мнимая восприимчивость χʺ, как функция от частоты поля ω, имеет один максимум, соответствующий вязкой релаксации магнитного момента частицы. Для достаточно жестких композитов мнимая компонента χʺ(ω) имеет один максимум, соответствующий неелевскому перемагничиванию частицы. При промежуточных значениях жесткости материала χʺ(ω) имеет два максимума. При β = π/2 (поле и ось магнитной анизотропии частицы перпендикулярны) магнитная часть восприимчивости имеет один максимум, соответствующий вязкой релаксации намагниченности.

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа поддержана грантом Министерства науки и образования РФ, проект FEUZ 2023-0020.

СОБЛЮДЕНИЕ ЭТИЧЕСКИХ СТАНДАРТОВ

В данной работе отсутствуют исследования человека или животных.

КОНФЛИКТ ИНТЕРЕСОВ

Авторы не имеют конфликта интересов.

About the authors

А. Ю. Зубарев

Уральский федеральный университет им. Б. Н. Ельцина

Email: Antoniusmagna@yandex.ru
Russian Federation, 620002, Екатеринбург, просп. Ленина, 51

Л. Ю. Искакова

Уральский федеральный университет им. Б. Н. Ельцина

Email: Antoniusmagna@yandex.ru
Russian Federation, 620002, Екатеринбург, просп. Ленина, 51

А. Ю. Мусихин

Уральский федеральный университет им. Б. Н. Ельцина

Author for correspondence.
Email: Antoniusmagna@yandex.ru
Russian Federation, 620002, Екатеринбург, просп. Ленина, 51

References

Supplementary files