Analytical formula and numerical calculation of the second harmonic of dynamic susceptibility in concentrated ferrofluids

Full Text

Введение

Суспензия магнитных наночастиц, диспергированных в немагнитной жидкости-носителе, называется магнитной жидкостью, или феррожидкостью [1]. Уникальное свойство такого коллоидного раствора заключается в сочетании текучести, характерной для жидкостей, со способностью к намагничиванию. Средний размер наночастиц составляет около 10 нм. Такую систему можно рассматривать как ансамбль одинаковых сферических диполей с диаметром магнитного ядра $d$ и магнитным моментом $m = v M_s$, где $M_s$ — намагниченность насыщения, а $v = \pi d^3 / 6$ — объем магнитного ядра. Под действием внешнего магнитного поля магнитные моменты, противодействуя тепловому движению, преимущественно ориентируются по направлению поля, что приводит к возникновению намагниченности системы. 

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

В переменном магнитном поле с амплитудой $h$ и частотой $\omega$ намагниченность можно представить в виде ряда по гармоникам:
\[\begin{equation}
    \frac{M(t)}{\rho m} = \sum_{p=1}^{\infty} \chi_p (h, \omega) \, e^{i p \omega t}, 
\tag{1}
\end{equation}\]
где $\rho$ — концентрация частиц. Коэффициенты ряда (1), зависящие от амплитуды и частоты переменного поля, определяют $p$-ю гармонику динамической восприимчивости. 

Динамическая восприимчивость представляет собой ключевой параметр, характеризующий отклик ансамбля магнитных частиц на внешнее магнитное поле. Исследование магнитного отклика имеет важное значение для многочисленных прикладных областей. В частности, мнимая часть динамической восприимчивости позволяет оценить интенсивность диссипации энергии, что рассматривается в медицинских приложениях как потенциальный метод терапии онкологических заболеваний [2–4]. Для клинического применения этого метода требуется повышение его эффективности и безопасности. Установлено, что наличие постоянного поля способствует более эффективной генерации тепла [5]. Однако в этом случае в спектре динамической восприимчивости возникают новые эффекты, обусловленные комбинацией переменного и постоянного полей [6]. 

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

Динамическая восприимчивость представляет собой активную область теоретических исследований. В классической работе [9] получено выражение для первой гармоники (линейной восприимчивости), описывающее отклик ансамбля невзаимодействующих дипольных частиц на слабое переменное поле:
\[\begin{equation}
    \chi_D = \frac{\chi_L}{1+i\omega\tau_B},
\tag{2}
\end{equation}\]
где $\tau_B = \pi d^3 \eta / (2 k_B T)$ — характерное время броуновской релаксации; $\eta$ — динамическая вязкость окружающей жидкости; $k_B T$ — тепловая энергия системы; $\chi_L = 8\lambda\varphi$ — статическая восприимчивость Ланжевена, пропорциональная константе дипольной связи $\lambda = \mu_0 m^2 / (4 \pi  d^3 k_B T)$ и объемной доле частиц $\varphi = \rho \pi d^3 / 6$; $\rho$ — объемная концентрация магнитных частиц.

Выражение для третьей гармоники (нелинейного отклика) в случае малых амплитуд переменного поля без учета межчастичных взаимодействий приведено в работе [10]. Особого внимания заслуживают исследования [11, 12], посвященные влиянию диполь-дипольных взаимодействий на первую и третью компоненты динамического отклика. В работах [12, 13] предложены упрощенные аппроксимационные выражения для первой и третьей гармоник, учитывающие как межчастичные взаимодействия, так и произвольные амплитуды переменного поля. Описание первой и третьей гармоник для произвольных амплитуд поля (без учета межчастичных взаимодействий) представлено в исследовании [14]. Влияние постоянного подмагничивающего поля на первую компоненту динамического отклика рассмотрено в работе [15]. Выражения для второй компоненты динамического отклика получены в исследованиях [16, 17], однако в этих работах взаимодействие между частицами не учитывалось.

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

1. Математическая модель

Рассматривается ансамбль из $N$ идентичных сферических магнитных диполей, диспергированных в жидкости-носителе. Вектор магнитного момента $k$-й частицы описывается выражением
\[
\mathbf{m}_k = m_k \hat{\mathbf{m}}_k,
\] 
где $\hat{\mathbf{m}}_k = (\sin \theta_k \cos \varphi_k, \sin \theta_k \sin \varphi_k, \cos \theta_k)$ — единичный вектор ориентации, а $m_k$ — модуль магнитного момента. Радиус-вектор частицы задается как 
\[
\mathbf{r}_k = r_k \hat{\mathbf{r}}_k,
\]
где $\hat{\mathbf{r}}_k = (\sin \zeta_k \cos \psi_k, \sin \zeta_k \sin \psi_k, \cos \zeta_k)$ — единичный вектор направления, а $r_k$ — длина радиус-вектора. Здесь $\theta_k$, $\varphi_k$, $\zeta_k$, $\psi_k$ — углы в сферической системе координат. 

Частицы находятся в бесконечно протяженном цилиндре, ориентированном вдоль оси $Oz$. Такая геометрия позволяет пренебречь краевыми эффектами и считать, что внешнее магнитное поле совпадает с полем внутри образца. Магнитное поле имеет вид
\[
\mathbf{H} = (h_{0} + h e^{i \omega t})\hat{\mathbf{H}},
\]
где $h$ — амплитуда переменного магнитного поля, $\omega$ — его частота, $h_0$ — напряженность постоянного поля, а $\hat{\mathbf{H}} = (0, 0, 1)$ — единичный вектор направления. Введем обозначение $x_k = \cos \theta_k$, где $\theta_k$ — угол между векторами $\hat{\mathbf{m}}_k$ и $\hat{\mathbf{H}}$.

Энергия зеемановского взаимодействия $k$-го магнитного момента с внешним полем определяется выражением
\[\begin{equation}
U_H = -\mu_0 (\mathbf{m}_k \cdot \mathbf{H}) = -\mu_0 m (h_{0} + h e^{i \omega t}) x_k.
\tag{3}
\end{equation}\]

Потенциал диполь-дипольного взаимодействия между $k$-й и $j$-й частицами имеет вид
\[\begin{equation}
U_d = -\frac{\mu_0 m^2}{4 \pi r_{k,j}^3} \left[3 (\hat{\mathbf{m}}_k \cdot \hat{\mathbf{r}}_{k,j}) (\hat{\mathbf{m}}_j \cdot \hat{\mathbf{r}}_{k,j}) - (\hat{\mathbf{m}}_k \cdot \hat{\mathbf{m}}_j)\right],
\tag{4}
\end{equation}\]
где $r_{k,j}$ — расстояние между центрами частиц.

Поскольку система обладает осевой симметрией, ориентация $k$-го магнитного момента зависит только от угла $\theta_k$. Эволюцию ориентации магнитного момента описывает уравнение Фоккера—Планка:
\[\begin{equation}
2 \tau_B \frac{\partial W}{\partial t} = \frac{\partial}{\partial x_k} 
\Bigl[ (1-x_k^2) \Bigl( \frac{\partial W}{\partial x_k} + W \frac{\partial U}{\partial x_k} \Bigr) \Bigr].
\tag{5}
\end{equation}\]

Решение уравнения (5) представляет собой функцию плотности вероятности ориентации магнитного момента $W = W(t, x_k)$, определяющую наиболее вероятное направление магнитного момента относительно внешнего поля в момент времени $t$. Данное решение удовлетворяет условию нормировки:
\[\begin{equation}
    \int_{-1}^1 W(t, x_k) dx_k = 1.
\tag{6}
\end{equation}\]

Для идеальной системы невзаимодействующих частиц нормированная на тепловую энергию потенциальная энергия $U=U(t, x_k)$ сводится к зеемановской энергии (3):
\[\begin{equation}
    U = U^{id} = U_{H} / k_B T = -(\xi_{0} + \xi e^{i \omega t})x_k,
\tag{7}
\end{equation}\]
где $\xi_0 = \mu_0 m h_0 / (k_B T)$ и $\xi = \mu_0 m h / (k_B T)$ — параметры Ланжевена для постоянного и переменного полей соответственно. 

Для учета межчастичных взаимодействий применяется метод среднего магнитного поля первого порядка (МСП-1), предложенный в работе [18]. Метод основан на предположении, что на $k$-ю частицу помимо внешнего магнитного поля действует дополнительное поле, создаваемое $j$-й частицей во всем объеме системы. В рамках данного подхода диполь-дипольные взаимодействия учитываются с точностью до вклада парных корреляций. Таким образом, потенциальная энергия системы взаимодействующих частиц включает два вклада: энергию взаимодействия магнитных моментов с внешними полями и энергию диполь-дипольных взаимодействий:
\[\begin{equation}
U = U^{int} = \biggl(U_{H} + \rho \int\!\!\!\int U_d(k,j) \Theta(k, j) W^{id} \, d\hat{\mathbf{m}}_j \, d\hat{\mathbf{r}}_j \biggr) / (k_B T),
\tag{8}
\end{equation}\]
где $W^{id}$ — решение уравнения (5) с потенциальной энергией (7). Интегралы по $d\hat{\mathbf{m}}_j$ и $d\hat{\mathbf{r}}_j$ представляют собой усреднение по всем возможным ориентациям магнитного момента и положениям $j$-й частицы соответственно:
\[\begin{equation}
    \begin{gathered}    
\int d\hat{\mathbf{r}}_j = \lim_{R\to\infty} \int_0^{2\pi} d\psi_j \int_0^{\pi} \sin\zeta_j d\zeta_j \int_0^{R/\sin\zeta_j} r_j^2 dr_j,
    \\
    \int d\hat{\mathbf{m}}_j = \frac{1}{4\pi} \int_{-1}^1 dx_j \int_0^{2\pi} d\varphi_j.
    \end{gathered}
\tag{9}
\end{equation}\]
Функция Хевисайда $\Theta(k,j) = \Theta(|\mathbf{r}_j - \mathbf{r}_k|)$ исключает возможность перекрытия $j$-й и $k$-й частиц.

Подстановка выражений (4) и (9) в  (8) с последующим интегрированием дает
\[\begin{equation}
    U^{int} = -\biggl[\xi_0 + \xi e^{i \omega t} + \chi_L \int_{-1}^{1} W^{id} x_j dx_j\biggr]x_k.
\tag{10}
\end{equation}\]

Намагниченность системы определяется как усредненная по всем возможным ориентациям проекция магнитного момента $k$-й частицы на направление магнитного поля с весовой функцией $W$:
\[\begin{equation}
 M(t)= \rho m \int (\hat{\mathbf{m}}_k \cdot \hat{\mathbf{H}}) W d\hat{\mathbf{m}}_k = \rho m \int_{-1}^1 W x_k dx_k.
\tag{11}
\end{equation}\]

Действительная и мнимая части второй гармоники находятся через коэффициенты ряда Фурье:
\[\begin{equation}
\mathcal{R}[\chi_2] = \frac{\omega}{\pi \xi^2}\int_{0}^{2\pi/\omega} \!\!\! M(t) \cos(2\omega t) dt, \quad 
\mathcal{I}[\chi_2] = \frac{\omega}{\pi \xi^2}\int_{0}^{2\pi/\omega} \!\!\! M(t) \sin(2\omega t) dt,
\tag{12}
\end{equation}\]
где $\mathcal{R}[x]$ и $\mathcal{I}[x]$ обозначают действительную и мнимую части величины $[x]$ соответственно.

2. Аналитическая формула второй гармоники в случае малых амплитуд переменного поля

Для аналитического решения уравнения (5) применяется метод разделения переменных. Функция плотности вероятности $W(t, x_k)$ представляется в виде бесконечного ряда по системе ортогональных функций (далее индекс $k$ опущен для упрощения записи):
\[\begin{equation}
    W(t, x) = W_0(\xi_0, x) + \sum_{l=1}^{\infty} T_l(t) P_l(x),
\tag{13}
\end{equation}\]
\[\begin{equation}
    W_0(\xi_0, x) = \frac{\xi_0}{2 sh\xi_0} e^{\xi_0 x}.
\tag{14}
\end{equation}\]

В качестве базиса разложения использованы полиномы Лежандра $P_l(x)$. Выражение (14) описывает равновесное распределение ориентаций магнитных моментов. Подстановка (13) в уравнение (5) с учетом свойств полиномов Лежандра приводит к бесконечной системе рекуррентных уравнений для коэффициентов $T_l$:
\[\begin{equation}
2\tau_B T'_l = -T_{l+1} \frac{\partial U}{\partial x} \frac{l(l+1)}{2l+3} + T_{l-1} \frac{\partial U}{\partial x} \frac{l(l+1)}{2l-1} - l(l+1)T_l - \xi e^{i\omega t} f_l(\xi_0),
\tag{15}
\end{equation}\]
где функция $f_l(\xi_0)$ определяется выражением
\[\begin{equation}
f_l(\xi_0) = \int_{-1}^1 P_l(x) \frac{\partial}{\partial x} \bigl[(1-x^2) W_0(\xi_0, x)\bigr] dx.
\tag{16}
\end{equation}\]

Решение для функций $T_l(t)$ ищется в виде разложения по гармоникам:
\[\begin{equation}
T_l(t) = \sum^{\infty}_{p=1} \chi_{l,p} e^{i p \omega t},
\tag{17}
\end{equation}\]
где коэффициенты разложения $\chi_{l,p}$ зависят от параметров $\xi_0$, $\xi$ и $\omega$.

Подстановка разложения (17) и выражения для потенциальной энергии (7) в систему (15) приводит к бесконечной системе уравнений относительно коэффициентов $\chi_{l,p}$ для случая невзаимодействующих частиц:
\[\begin{multline}
\bigl(2 i k \omega \tau_B + l(l+1)\bigr) \chi_{l,p} = -\xi \frac{l(l+1)}{2l+3} \chi_{l+1,p-1} + \xi \frac{l(l+1)}{2l-1} \chi_{l-1,p-1} -{}
\\
{}-\xi_0 \frac{l(l+1)}{2l+3} \chi_{l+1,p} + \xi_0 \frac{l(l+1)}{2l-1} \chi_{l-1,p} - \delta_{1,p}\,\xi f_l(\xi_0),
\tag{18}
\end{multline}\]
где $\delta_{1,p}$ — символ Кронекера.

Выражение для намагниченности (11) после подстановки разложения (13) с учетом свойства ортогональности полиномов Лежандра принимает вид
\[\begin{equation}
    M(t) = \rho m \int_{-1}^1 x W dx = \rho m \biggl(L(\xi_0) + \frac{2}{3} \sum^{\infty}_{p=1} \chi_{1,p} e^{i p\omega t} \biggr),
\tag{19}
\end{equation}\]
где $L(x) = cth x - 1/x$ — функция Ланжевена.

Действительная и мнимая части второй гармоники, вычисляемые согласно (12), после подстановки (19) выражаются через коэффициент $\chi_{1,2}$:
\[\begin{equation}
\begin{gathered}
\mathcal{R}[\chi_2] = \frac{\omega}{\pi \xi^2}\int_{0}^{2\pi/\omega} M(t) \cos(2\omega t) dt = \frac{k\chi^2_L}{\xi^2}\mathcal{R} [\chi_{1,2}], \\
\mathcal{I}[\chi_2] = \frac{\omega}{\pi \xi^2}\int_{0}^{2\pi/\omega} M(t) \sin(2\omega t) dt = \frac{k\chi^2_L}{\xi^2}\mathcal{I} [\chi_{1,2}],
\end{gathered}
\end{equation}\]
где $k = 12m/\rho$. Таким образом, вторая гармоника динамической восприимчивости полностью определяется коэффициентом $\chi_{1,2}$.

Система уравнений (18) содержит бесконечное число неизвестных. Для получения аналитического решения необходимо ограничить систему конечным числом уравнений. Учитывая условие малости амплитуды переменного поля ($\xi \ll 1$), можно показать, что вклад гармоники с номером $p$ пропорционален $\xi^p$. 

Анализ показывает, что для определения гармоник выше первой необходимо рассматривать как минимум два уравнения системы. Введем приближение, полагая $\chi_{l,p} = 0$ при $l > 2$ и $p > 2$. Для случая невзаимодействующих диполей решение системы (18) дает следующие выражения для коэффициентов:
\[\begin{equation}
\chi^{id}_{1,1} = \frac{\xi }{2} \frac{\frac{1}{5}\xi_0 f_2(\xi_0) - f_1(\xi_0)(3 + i\omega {\tau_B })}{(1+i\omega {\tau_B})(3+i\omega {\tau_B} ) + \frac{1}{5}\xi_0^2},
\tag{20}
\end{equation}\]
\[\begin{equation}
\chi^{id}_{1,2} = -\frac{\xi^2}{10}\frac{f_2(\xi_0)\left(\frac{1}{5}\xi_0^2-(1+i\omega\tau_B)(3+2i\omega\tau_B)\right)-3f_1(\xi_0)\xi_0(2+i\omega\tau_B)}{\left[(1+2i\omega\tau_B)(3+2i\omega\tau_B)+\frac{1}{5}\xi_0^2\right]\left[(1+i\omega\tau_B)(3+i\omega\tau_B)+\frac{1}{5}\xi_0^2\right]},
\tag{21}
\end{equation}\]
где функции $f_1(x)$ и $f_2(x)$, полученные из (16), имеют вид
\[
f_1(x) = -\frac{3L(x)}{x}, \quad f_2(x) = -\frac{15L_3(x)}{x}, \quad L_3(x) = 1 - \frac{3L(x)}{x}.
\]

В предельном случае $\xi_0 \to 0$ коэффициент (20) воспроизводит формулу Дебая (2), тогда как коэффициент (21) стремится к нулю. Таким образом, аналитическое решение уравнения (5) для невзаимодействующих диполей в случае малых амплитуд переменного поля принимает форму
\[\begin{equation}
W^{id}(t,x) = W_0 + (\chi^{id}_{1,1}e^{i\omega t} + \chi^{id}_{1,2}e^{i2\omega t})P_1(x) + \cdots
\tag{22}
\end{equation}\]

Рассмотрим вклад диполь-дипольных взаимодействий в потенциальную энергию. Подстановка решения (22) в выражение (10) с последующим интегрированием дает модифицированное выражение для потенциальной энергии, учитывающее межчастичные взаимодействия:
\[\begin{equation}
U = U^{int} = -\bigl[\widetilde{\xi}_0 + A e^{i\omega t} + B e^{i2\omega t}\bigr]x,
\tag{23}
\end{equation}\]
где введены следующие обозначения:
\[
\widetilde{\xi}_0 = \xi_0 + \chi_L L(\xi_0),
\]
\[
A = \xi \Bigl(1 + \frac{\chi_L}{3} \frac{\frac{1}{5}\xi_0 f_2(\xi_0) - f_1(\xi_0)(3 + i\omega {\tau_B})}{(1+i\omega {\tau_B})(3+i\omega{\tau_B}) + \frac{1}{5}\xi_0^2}\Bigr), 
\]
\[
B = -\xi^2\frac{\chi_L}{15}\frac{f_2(\xi_0)\left(\frac{1}{5}\xi_0^2-(1+i\omega\tau_B)(3+2i\omega\tau_B)\right)-3f_1(\xi_0)\xi_0(2+i\omega\tau_B)}{\left[(1+2i\omega\tau_B)(3+2i\omega\tau_B)+\frac{1}{5}\xi_0^2\right]\left[(1+i\omega\tau_B)(3+i\omega\tau_B)+\frac{1}{5}\xi_0^2\right]}.
\]

Решение снова ищем в виде разложения (13), используя модифицированное начальное распределение $W_0(\widetilde{\xi}_0, x)$. Подстановка потенциальной энергии (23) в систему (15) приводит к новой системе рекуррентных соотношений для коэффициентов $\chi_{l,p}$, учитывающей межчастичные взаимодействия:
\[\begin{multline}
\bigl(2ik\omega\tau_B + l(l+1)\bigr)\chi_{l,p} = -A\frac{l(l+1)}{2l+3}\chi_{l+1,p-1} + A\frac{l(l+1)}{2l-1}\chi_{l-1,p-1} - {} 
\\
{}-\widetilde{\xi}_0\frac{l(l+1)}{2l+3}\chi_{l+1,p} + \widetilde{\xi}_0\frac{l(l+1)}{2l-1}\chi_{l-1,p} - {}
\\
{}-B\frac{l(l+1)}{2l+3}\chi_{l+1,p-2} + B\frac{l(l+1)}{2l-1}\chi_{l-1,p-2} - \delta_{1,p}Af_l(\widetilde{\xi}_0) - \delta_{2,p}Bf_l(\widetilde{\xi}_0).
\tag{24}
\end{multline}\]

Решение системы (24) находится аналогичным методом. Вторая гармоника динамической восприимчивости с учетом диполь-дипольных взаимодействий принимает вид
\[\begin{multline}
\chi_{1,2} = -\frac{A^2}{10}\frac{f_2(\widetilde{\xi}_0)\bigl(\frac{1}{5}\widetilde{\xi}_0^2-(1+i\omega\tau_B)(3+2i\omega\tau_B)\bigr)-3\widetilde{\xi}_0f_1(\widetilde{\xi}_0)(2+i\omega\tau_B)}{\bigl[(1+2i\omega\tau_B)(3+2i\omega\tau_B)+\frac{1}{5}\widetilde{\xi}_0^2\bigr]\bigl[(1+i\omega\tau_B)(3+i\omega\tau_B)+\frac{1}{5}\widetilde{\xi}_0^2\bigr]} + {}
\\
{}+\frac{B}{2}\frac{\frac{1}{5}\widetilde{\xi}_0f_2(\widetilde{\xi}_0) - f_1(\widetilde{\xi}_0)(3+2i\omega\tau_B)}{(1+2i\omega\tau_B)(3+2i\omega\tau_B) + \frac{1}{5}\widetilde{\xi}_0^2}.
\tag{25}
\end{multline}\]

3. Численное решение уравнения Фоккера—Планка

Аналитическое представление функции $W$ с использованием конечного числа уравнений в системе (17) ограничивает область параметров, где выражение (25) адекватно описывает вторую гармонику. Для определения границ применимости формулы (25) было выполнено численное решение уравнения Фоккера—Планка с последующим расчетом второй гармоники.

Численное решение уравнения (5) проводилось с использованием конечно-разностной схемы, предложенной в работе [19]. Изначально разработанная для уравнений конвекции-диффузии, данная схема успешно адаптируется для уравнения Фоккера—Планка, где первый член интерпретируется как диффузионный, а второй — как конвекционный. Сходимость и устойчивость схемы доказаны в [19], а ее применение продемонстрировано в работах [20, 21] для систем с неелевской релаксацией.

Численное решение получено для диапазона частот $10^{-2} \leqslant \omega\tau_B \leqslant 10^{2}$ с логарифмическим шагом 0.1 по десятичным порядкам. Пространственно-временная сетка задавалась следующим образом:
\[\begin{equation}
    \begin{gathered}
        t_n = t_{n-1} + h_t, \quad t_0 = 0, \quad h_t = \frac{2\pi}{\omega} 10^{-4}; \\
        x_m = x_{m-1} + h_x, \quad x_0 = -1 + \frac{h_x}{2}, \quad h_x = 0.001,
    \end{gathered}
\end{equation}\]
где $h_x$ и $h_t$ — шаги по пространственной и временной координатам соответственно. Индексы $n$ и $m$ изменяются от 0 до $N_n = 10^4$ и $N_m = 2/h_x - 1$ соответственно. Значения функции в узлах сетки обозначены как $W_{n,m} = W(t_n,x_m)$.

Дискретный аналог уравнения (5) согласно [19] имеет вид
\[\begin{equation}
\frac{e^{-\delta h_t}W_{n,m}-W_{n-1,m}}{h_t} + (D + C + \delta)(e^{-\delta h_t}W_{n,m}) = 0,
\tag{26}
\end{equation}\]
где операторы диффузии $D$ и конвекции $C$ определяются следующими выражениями:
\[
DW_{n,m} = -f\Bigl(x_m + \frac{h_x}{2}\Bigr)\frac{W_{n,m+1}-W_{n,m}}{h_x^2} +  f\Bigl(x_m - \frac{h_x}{2}\Bigr)\frac{W_{n,m}-W_{n,m-1}}{h_x},
\]
\[
CW_{n,m} = v\Bigl(t^*,x_m + \frac{h_x}{2}\Bigr)\frac{W_{n,m+1}+W_{n,m}}{2h_x} - v\Bigl(t^*,x_m - \frac{h_x}{2}\Bigr)\frac{W_{n,m}+W_{n,m-1}}{2h_x},
\]
\[
f(x) = 1-x^2, \quad t^* = \frac{t_{n+1}+t_n}{2}, \quad v(t,  x) = \frac{\partial U}{\partial x}(1-x^2).
\]

Численное решение включает следующие этапы:

1. решение уравнения (5) с потенциальной энергией (7) для получения $W^{id}$ — распределения без учета межчастичных взаимодействий;
2. вычисление интегрального вклада диполь-дипольных взаимодействий:
      \[\begin{equation}
       I = \chi_L \int_{-1}^{1} W^{id} x dx;
       \end{equation}\]
3. повторное решение уравнения (5) с модифицированной потенциальной энергией (10), учитывающей взаимодействия;
4. численное интегрирование для определения намагниченности (11) и компонент восприимчивости (12) методом трапеций.

Параметр стабилизации $\delta$ в схеме (26) выбирается согласно [19]:
\[
\delta = \frac{1}{2} \max |v'_x|.
\]
В расчетах использовались следующие значения:

• $\delta = \xi + \xi_0/2$ — для невзаимодействующих частиц;
• $\delta = \xi + \xi_0/2 + I/2$ — с учетом взаимодействий частиц. 

Начальное распределение задается выражением (14). Для достижения термодинамического равновесия, исключающего зависимость от начальных условий, моделирование проводится до $t_f = 10T$, где $T = 2\pi/\omega$ — период переменного поля. Установившееся решение определяется усреднением $W$ на интервале $[9T,10T]$ с нормировкой, обеспечивающей выполнение (6).

4. Сравнение аналитической формулы и численного решения

На рис. 1 представлены зависимости действительной (a) и мнимой (b) частей динамической восприимчивости (21) для трех значений напряженности постоянного поля: $\xi_0 = 0.1$, $0.5$ и $1$. Численное решение уравнения (5), показанное точками (символами), получено при амплитуде переменного поля $\xi = 0.1$. Наблюдается хорошее согласие аналитического и численного решений при $\xi_0 \lesssim 0.5$. Расхождение при больших значениях $\xi_0 > 0.5$ обусловлено усечением системы (18) до двух уравнений.

Для верификации результатов на графиках приведены данные из работы [16] (крестики), полностью совпадающие с выражением (21), что подтверждает корректность полученных результатов.

Рис. 1. Действительная (a) и мнимая (b) части динамической восприимчивости (21) для значений напряженности постоянного поля $\xi_0 = 0.1$, $0.5$ и $1$ и амплитуде переменного поля $\xi = 0.1$: линии — аналитическое решение по (21); символы — численные данные; крестики — результаты [16]
[Figure 1. Real (a) and imaginary (b) parts of the dynamic susceptibility (21) for DC field strengths $\xi_0 = 0.1$, $0.5$, and $1$ and AC field amplitude $\xi = 0.1$: lines — analytical solution by (21); symbols — numerical data; crosses — results from [16]]

На рис. 2 показаны зависимости действительной (a) и мнимой (b) частей восприимчивости (25), учитывающей межчастичные взаимодействия, при $\xi = \xi_0 = 0.1$ для различных значений восприимчивости Ланжевена: $\chi_L = 0.5$, $1$ и $3$. Численное решение уравнения (5) представлено точками. В случае $\chi_L=0.5$, соответствующем разбавленной феррожидкости со слабыми корреляциями, результаты (25) (сплошные линии) близки к решению (21) для невзаимодействующей системы (пунктирные линии). Наблюдается хорошее соответствие аналитического и численного решений во всем исследованном диапазоне $\chi_L \lesssim 3$.

Рис. 2. Действительная (a) и мнимая (b) части восприимчивости (25) для $\chi_L = 0.5$, $2$ и $3$: сплошные линии — аналитическое решение по (25); символы — численные данные; пунктирные линии — аналитическое решение по (21). Значения амплитуды переменного поля и напряженности постоянного поля равны: $\xi = \xi_0 = 0.1$
[Figure 2. Real (a) and imaginary (b) parts of the susceptibility (25) for $\chi_L = 0.5$, $2$ and $3$: solid lines — analytical solution by (25); symbols — numerical data; dashed lines — analytical solution by (21). The AC field amplitude and DC field strength are equal: $\xi = \xi_0 = 0.1$]

Рассмотрим поведение статической восприимчивости и сдвиг положения максимума мнимой части в зависимости от напряженности постоянного магнитного поля $\xi_0$ при различных значениях восприимчивости Ланжевена $\chi_L$. 

На рис. 3, a представлены следующие данные:  сплошными линиями показана зависимость действительной части аналитического решения (25) в статическом пределе ($\omega\tau_B \to 0$) для значений $\chi_L = 1$ и $2$; точками обозначены результаты численного решения уравнения (5) при малой частоте $\omega\tau_B = 0.01$ для тех же значений $\chi_L$; для сравнения пунктирной линией приведена зависимость, полученная по формуле (21), не учитывающая межчастичные взаимодействия, и сопутствующие численные данные.

Анализ полученных зависимостей позволяет выявить существование выраженного локального максимума вблизи значения $\xi_0 \approx 1$. В области относительно слабых полей ($\xi_0 \lesssim 1$) учет диполь-дипольных взаимодействий приводит к существенному возрастанию статической восприимчивости второй гармоники, причем наиболее значительный эффект наблюдается в окрестности указанного максимума. 
При превышении напряженности поля значения $\xi_0 > 1$ влияние межчастичных корреляций постепенно ослабевает. При дальнейшем увеличении $\xi_0$ основную роль начинает играть взаимодействие магнитных моментов с внешним постоянным полем, что приводит к практическому совпадению статических восприимчивостей взаимодействующей и невзаимодействующей систем.

Рис. 3. Зависимость статической восприимчивости (a) и положения максимума мнимой части (b) от величины напряженности постоянного поля $\xi_0$. Восприимчивость Ланжевена принимает значения $\chi_L = 1$ и 2. Сплошные линии обозначают максимум, найденный по формуле (25), которая учитывает межчастичные взаимодействия. Пунктирные линии — максимум, определенный по формуле (21) для системы без учета взаимодействий. Символы — результат численного решения уравнения (5)
[Figure 3. Dependence of static susceptibility (a) and position of the imaginary part maximum (b) on DC field strength $\xi_0$. The Langevin susceptibility takes values $\chi_L = 1$ and $2$. Solid lines show maxima calculated by (25) accounting for interparticle interactions. Dashed lines — maxima from (21) for non-interacting systems. Symbols — numerical solution of equation (5)]

На рис. 3, b представлена зависимость характерной частоты $\omega^{*}\tau_B$, соответствующей максимуму мнимой части динамической восприимчивости, от безразмерной напряженности постоянного поля $\xi_0$ для систем с восприимчивостью Ланжевена $\chi_L = 1$ и $2$. Сплошные линии отражают аналитическую зависимость, полученную по формуле (25), учитывающей межчастичные взаимодействия, тогда как точки соответствуют результатам численного моделирования на основе уравнения (5). Для сравнения приведены аналогичные зависимости, рассчитанные по формуле (21) без учета диполь-дипольных взаимодействий (пунктирные линии), и соответствующие им численные данные (точки).

Анализ представленных зависимостей демонстрирует следующие особенности поведения системы. В области слабых полей ($\xi_0 \lesssim 1.5$) увеличение восприимчивости Ланжевена вызывает сдвиг положения максимума в область более низких частот, что обусловлено усилением роли межчастичных взаимодействий. В точке $\xi_0 \approx 1.5$ наблюдается баланс между взаимодействием магнитных моментов с внешним полем и диполь-дипольными взаимодействиями. При дальнейшем росте напряженности поля ($\xi_0 > 1.5$) преобладающее влияние внешнего поля приводит к смещению максимума мнимой части в область высоких частот.

5. Заключение

В данной работе получена аналитическая формула для второй гармоники динамической восприимчивости ансамбля взаимодействующих магнитных частиц в слабом переменном поле. Сравнение теоретических результатов с численным решением уравнения Фоккера—Планка показало, что  при значениях параметра Ланжевена $\xi_0 \lesssim 0.5$, характеризующего отношение магнитной энергии к тепловой, между ними наблюдается хорошее соответствие. При $\xi_0 > 0.5$ наблюдается снижение точности аналитического описания. 

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

Анализ выявил существенное влияние межчастичных взаимодействий на динамический отклик системы. В зависимости от напряженности постоянного поля $\xi_0$ диполь-дипольные взаимодействия могут как усиливать, так и ослаблять отклик. Направление сдвига частоты $\omega^{*}\tau_B$, соответствующей максимуму мнимой части восприимчивости, также определяется величиной приложенного постоянного поля.

Конкурирующие интересы. Конкурирующих интересов нет.
Авторский вклад и ответственность. Работа полностью выполнена автором. Автор несет полную ответственность за предоставление законченной версии рукописи в печать. Окончательная версия рукописи одобрена автором.
Финансирование. Исследование выполнено за счет гранта Российского научного фонда № 23–12–00039, https://rscf.ru/project/23-12-00039.
Благодарность. Выражаю глубокую благодарность моему научному руководителю д.ф.-м.н. Е. А. Елфимовой.

About the authors

Mikhail S. Rusanov

Author for correspondence.
Email: mikhail.rusanov@urfu.ru
ORCID iD: 0000-0001-7439-8179
SPIN-code: 9081-7507
Scopus Author ID: 57306389700
ResearcherId: JNS-8758-2023
https://www.mathnet.ru/rus/person229346

Research Engineer; Lab. of Mathematical Modeling of Physicochemical Processes in Multiphase Media; Dept. of Mathematics, Mechanics and Computer Science; Institute of Natural Sciences and Mathematics

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Figure 1. Real (a) and imaginary (b) parts of the dynamic susceptibility (21) for DC field strengths $\xi_0 = 0.1$, $0.5$, and $1$ and AC field amplitude $\xi = 0.1$: lines — analytical solution by (21); symbols — numerical data; crosses — results from [16]

Download (281KB)

Indexing metadata

3. Figure 2. Real (a) and imaginary (b) parts of the susceptibility (25) for $\chi_L = 0.5$, $2$ and $3$: solid lines — analytical solution by (25); symbols — numerical data; dashed lines — analytical solution by (21). The AC field amplitude and DC field strength are equal: $\xi = \xi_0 = 0.1$

Download (268KB)

Indexing metadata

4. Figure 3. Dependence of static susceptibility (a) and position of the imaginary part maximum (b) on DC field strength $\xi_0$. The Langevin susceptibility takes values $\chi_L = 1$ and $2$. Solid lines show maxima calculated by (25) accounting for interparticle interactions. Dashed lines — maxima from (21) for non-interacting systems. Symbols — numerical solution of equation (5)

Download (271KB)

Indexing metadata