В силу несжимаемости жидкости ее объем жидкости не, поэтому $\langle \Theta \rangle = 0$

Индуктивность магнитного поля перпендикулярна поверхности парамагнетика, поэтому её величина одинакова во всём пространстве.
\[H_0(+0) = \frac{B_0}{\mu_0}, \quad H_0(-0) = \frac{B_0}{\mu \mu_0} \quad \Rightarrow \quad \Delta H = \frac{B_0}{\mu_0} \frac{\mu - 1}{\mu}\]

$\vec{b} = \mu \mu_0 \vec{h}$. Если мы берем дивергенцию от правой части в области, которая не включает границу, то $\operatorname{div} \vec{b} = \mu \mu_0 \operatorname{div} \nabla \phi = \mu \mu_0 \left(\partial_x^2 + \partial_y^2 + \partial_z^2\right) \phi = \mu \mu_0 \Delta \phi = 0$. Получаем уравнение Пуассона в обоих областях.

Нормаль к поверхности $\vec{n}$ имеет координаты $(- \partial_x \zeta, -\partial_y \zeta,1)$. Поэтому граничное условие $\vec{B}_\uparrow \cdot \vec{n}= \vec{B}_\downarrow \cdot \vec{n}$ сводится к $b_{z\uparrow} = b_{z \downarrow}$.

Нормаль к поверхности $\vec{n}$ имеет координаты $(- \partial_x \zeta, -\partial_y \zeta,1)$. Поэтому граничное условие $\vec{H}_\uparrow \times \vec{n}= \vec{H}_\downarrow \times \vec{n}$ сводится к системе уравнений
\[\begin{cases}
h_{x\uparrow} + H_{0\uparrow} \partial_x \zeta = h_{x\downarrow} + H_{0\downarrow} \partial_x \zeta \\
h_{y\uparrow} + H_{0\uparrow} \partial_y \zeta = h_{y\downarrow} + H_{0\downarrow} \partial_y \zeta
\end{cases}.\]

Напомним, что $\phi(x,y,z)$ создает малые возмущения магнитного поля во всем пространстве, в т.ч. при $z = \pm \infty$.

Функциональный вид $u_i(z)$ следует из уравнения Пуассона для $\phi$. $\Delta (u_i(z) \cos (k_{ix} x + k_{iy} y)) = ( \partial_z^2 - k_i^2) u_i(z) \cos (k_{ix} x + k_{iy} y)) $, то есть $u_{i\uparrow}(z) = u_{i\uparrow} e^{-kz}$ и $u_{i\downarrow}(z) = u_{i\downarrow} e^{kz}$ (в каждом случае остаётся только одна экспонента из условия малости возмущения при удалении на бесконечность).

Запишем граничные условия из \textbf{B3}:
\[ b_{z\uparrow}(\zeta) = \mu_0 \partial_z \phi_\uparrow = -\mu_0 \sum u_{i\uparrow} k_i e^{-k_i \zeta} \cos (k_{ix} x + k_{iy} y)\]\[ b_{z\downarrow}(\zeta) = \mu \mu_0 \partial_z \phi_\downarrow = \mu \mu_0 \sum u_{i\downarrow} k_i e^{k_i \zeta} \cos (k_{ix} x + k_{iy} y)\]В ведущем приближении экспоненты $e^{\pm k\zeta}\simeq 1$ и тогда
\[-u_{i\uparrow} = \mu u_{i\downarrow}.\]Граничные условия из \textbf{B4} выглядят абсолютно одинаково:
\[
\begin{split}\sum u_{i\uparrow} k_{ix} e^{-k \zeta} \sin (k_{ix} x + k_{iy} y) + \Delta H \sum a_i k_{ix} \sin (k_{ix} x + k_{iy} y) =\\= \sum u_{i\downarrow} k_{ix} e^{k \zeta} \sin (k_{ix} x + k_{iy} y)
\end{split}.\]

Снова в ведущем приближение $e^{\pm k \zeta} \simeq 1$ и
\[u_{i\uparrow} + \Delta H a_i = u_{i\downarrow}.\]Решим систему уравнений для $u_{i\uparrow}$ и $u_{i\downarrow}$ и получим
\[ u_{i\uparrow} = -a_i \frac{\mu \Delta H}{\mu + 1}, \quad u_{i\downarrow} = a_i \frac{\Delta H}{\mu + 1}.\]

Запишем разность энергий:
\[
\Delta U_B = - \frac{\mu-1}{2\mu \mu_0\sqrt{\sigma \rho g}}\left(\left \langle \int\limits_{-\infty}^{\zeta} b_z B_0 \, dz \right\rangle + \left \langle \int\limits_{-\infty}^{\zeta} B_0^2 \, dz \right\rangle - \left \langle \int\limits_{-\infty}^{0} B_0^2 \, dz \right\rangle\right).
\]
\[
\Delta U_B = - \frac{\mu-1}{2\mu \mu_0\sqrt{\sigma \rho g}}\left(\left \langle \int\limits_{-\infty}^{\zeta} b_z B_0 \, dz \right\rangle + \left \langle \int\limits_{0}^{\zeta} B_0^2 \, dz \right\rangle\right)
\]

\[
\Delta U_B = - \frac{\mu-1}{2\mu \mu_0\sqrt{\sigma \rho g}}\left(\left \langle \int\limits_{-\infty}^{\zeta} b_z B_0 \, dz \right\rangle + \left \langle \zeta\right\rangle B_0^2\right)
\]

С учётом $\left<\zeta\right>$ = 0 получаем итоговый ответ.

\[
\begin{split}
\Delta U_B = - \frac{\mu-1}{2\mu \mu_0\sqrt{\sigma \rho g}} \left\langle \int\limits_{-\infty}^{\zeta} b_z B_0 \, dz \right\rangle = - \frac{(\mu-1)B_0}{2\sqrt{\sigma \rho g}} \left\langle \int\limits_{-\infty}^{\zeta} \, d\phi_{\downarrow}\right\rangle = - \frac{(\mu-1)B_0}{2\sqrt{\sigma \rho g}} \left\langle \phi_\downarrow(\zeta) \right\rangle = \\ =
- \frac{(\mu-1)^2 B_0^2}{2 \sqrt{\sigma \rho g} \mu (\mu+1) \mu_0} \left\langle \zeta e^{k \zeta} \right\rangle \simeq - \frac{(\mu-1)^2 B_0^2}{2 \sqrt{\sigma \rho g} \mu (\mu+1) \mu_0}k \left\langle \zeta^2 \right\rangle
\end{split}\]

Рассмотрим $\zeta$ как сумму:
\[\zeta = \frac{A}{6} \left( e^{i\vec{k}_1 \vec{r}} + e^{-i\vec{k}_1 \vec{r}} +e^{i\vec{k}_2 \vec{r}} + e^{-i\vec{k}_2 \vec{r}} + e^{i\vec{k}_3 \vec{r}} + e^{-i\vec{k}_3 \vec{r}} \right),\]тогда $\zeta^2$ состоит из 36-ти слагаемых вида $e^{i(\pm \vec{k}_i \pm \vec{k}_j) \vec{r}}$. При усреднении по площади обнулются все слагаемые, где $\pm \vec{k}_i \pm \vec{k}_j \neq 0$, а остальные шесть равны единице. Поэтому
\[\langle \zeta^2 \rangle = \frac{A^2}{6}.\]Теперь $\nabla \zeta$:
\[ \nabla \zeta = \frac{iA}{6}\left[ \vec{k}_1 \left( e^{i\vec{k}_1 \vec{r}} - e^{-i\vec{k}_1 \vec{r}} \right) + \vec{k}_2 \left( e^{i\vec{k}_2 \vec{r}} - e^{-i\vec{k}_2 \vec{r}} \right) + \vec{k}_3 \left( e^{i\vec{k}_3 \vec{r}} - e^{-i\vec{k}_3 \vec{r}} \right)\right].\]Тогда
\[
\begin{split}
(\nabla \zeta)^2 = -\frac{A^2 k^2}{36}\left[ \left( e^{i\vec{k}_1 \vec{r}} - e^{-i\vec{k}_1 \vec{r}} \right)^2 + \left( e^{i\vec{k}_2 \vec{r}} - e^{-i\vec{k}_2 \vec{r}} \right)^2 + \left( e^{i\vec{k}_3 \vec{r}} - e^{-i\vec{k}_3 \vec{r}} \right)^2 \right] -\\- \frac{A^2}{36} \left[ \vec{k}_1 \cdot \vec{k}_2 \left( e^{i\vec{k}_1 \vec{r}} - e^{-i\vec{k}_1 \vec{r}} \right) \left( e^{i\vec{k}_2 \vec{r}} - e^{-i\vec{k}_2 \vec{r}} \right) + \dots \right].
\end{split}\]
Легко видеть, что после усреднения по площади обнулятся все слагаемые со скалярными произведениями $\vec{k}_i \cdot \vec{k}_j, i \neq j$. Среди оставшихся будет шесть равных $-1$, поэтому
\[ \langle (\nabla \zeta)^2 \rangle = \frac{k^2 A^2}{6}.\]

Обозначим
\[ \frac{(\mu-1)^2 B_0^2}{\sqrt{\sigma \rho g} \mu (\mu+1) \mu_0} = 2k_0\]. Среди всех $k$ выберем такой, при котором выражение в скобках принимает минимальное значение:
$k_\mathrm{min} = k_0$, при этом
\[U = 1 + \frac{a^2}{12} \left( 1 - k_0^2 \right),\] то есть при $k_0 > 1$ наличие гексагональной решетки с $k=k_0$ начинает быть энергетически выгодным.