Magnetic micro-droplet in rotating field: numerical simulation and comparison with experiment

TL;DR

This study develops a 3D boundary-integral algorithm to simulate magnetic droplet shapes in a rotating magnetic field, validating results with experiments and revealing key shape transitions and behaviors.

Contribution

The paper introduces a novel numerical method for predicting 3D equilibrium shapes of magnetic droplets under rotating fields, validated against experimental data.

Findings

Good agreement between simulations and experiments

Identification of oblate-prolate shape transition

Observation of star-fish like equilibrium shapes

Abstract

Magnetic droplets obtained by induced phase separation in a magnetic colloid show a large variety of shapes when exposed to an external field. However, the description of shapes is often limited. Here we formulate an algorithm based on three dimensional boundary-integral equations for strongly magnetic droplets in a high-frequency rotating magnetic field, allowing us to find their figures of equilibrium in three dimensions. The algorithm is justified by a series of comparisons with known analytical results. We compare the calculated equilibrium shapes with experimental observations and find a good agreement. The main features of these observations are the oblate-prolate transition, the flattening of prolate shapes with the increase of magnetic field strength and the formation of star-fish like equilibrium shapes. We show both numerically and in experiments that the magnetic droplet behaviour may be described with a tri-axial ellipsoid approximation. Directions for further research are mentioned, including the dipolar interaction contribution to the surface tension of the magnetic droplets, account for the large viscosity contrast between the magnetic droplet and the surrounding fluid.

Equations

$$-\nabla p+\frac{\partial {\mathbf{\sigma }}_k^v}{\partial x_k}+\frac{\partial {\mathbf{T}}_k}{\partial x_k}=0;\nabla \cdot \mathbf{v}=0,$$

$$E=-\frac{1}{2}\int \mathbf{M}\cdot {\mathbf{H}}_{0}dV+\gamma S,$$

$${\delta }_{L}f=f^{\prime }({\mathbf{x}}^{\prime })-f(\mathbf{x});{\mathbf{x}}^{\prime }=\mathbf{x}+\mathbf{\xi }(\mathbf{x})$$

$${\delta }_L\frac{\partial f}{\partial x_i}=\frac{\partial {\delta }_{L}f}{\partial x_i}-\frac{\partial {\xi }_k}{\partial x_i}\frac{\partial f}{\partial x_k},$$

$${\delta }_L{H}_i=\frac{\partial {\delta }_L\psi }{\partial x_i}-H_k\frac{\partial {\xi }_k}{\partial x_i}.$$

$$\frac{1}{4\pi }{B}_i{\delta }_L{H}_i=\frac{1}{4\pi }{B}_i\frac{\partial {\delta }_L\psi }{\partial x_i}-\frac{1}{4\pi }{H}_i{B}_k\frac{\partial {\xi }_i}{\partial x_k}.$$

$$\frac{1}{4\pi}\int B_{i}\delta_{L}H_{i}dV=\int\xi_{i}\frac{\partial}{\partial x_{k}}\Bigl{(}\frac{H_{i}B_{k}}{4\pi}\Bigr{)}dV-\int\Bigl{(}\frac{H_{i}^{i}B_{k}^{i}}{4\pi}-\frac{H_{i}^{e}B_{k}^{e}}{4\pi}\Bigr{)}\xi_{i}n_{k}dS~{}.$$

$$-\int\xi_{i}\frac{\partial}{\partial x_{k}}\Bigl{(}\frac{H^{2}}{8\pi}\delta_{ik}\Bigr{)}dV+\int\xi_{i}n_{k}\Bigl{(}\frac{(H^{i})^{2}}{8\pi}-\frac{(H^{e})^{2}}{8\pi}\Bigr{)}\delta_{ik}dS=0$$

$$\frac{1}{4\pi }\int B_i{\delta }_L{H}_{i}dV=\int {\xi }_i\frac{\partial T_{ik}}{\partial x_k}dV-\int {\xi }_i{n}_k(T_{ik}^i-T_{ik}^e)dS.$$

$$\int v_i({\sigma }_{ik}^i+T_{ik}^i-{\sigma }_{ik}^e-T_{ik}^e)n_{k}dS+\frac{1}{4\pi }\int \mathbf{B}\cdot \frac{{\delta }_L\mathbf{H}}{\delta t}dV-\int \frac{\partial v_i}{\partial x_k}{\sigma }_{ik}^{v}dV=0.$$

$$(\sigma^{i}_{ik}+T^{i}_{ik}-\sigma^{e}_{ik}-T^{e}_{ik})n_{k}=-\gamma\Bigl{(}\frac{1}{R_{1}}+\frac{1}{R_{2}}\Bigr{)}n_{i}~{},$$

$$\frac{d}{dt}(\tilde{F}+\gamma S)=-\int \frac{\partial v_i}{\partial x_k}{\sigma }_{ik}^{v}dV.$$

$$\tilde{F}+\frac{1}{8\pi }\int H_0^{2}dV=-\frac{1}{2}\int \mathbf{M}\cdot {\mathbf{H}}_{0}dV,$$

$$\frac{d}{dt}\Bigl{(}\gamma S-\frac{1}{2}\int\boldsymbol{M}\boldsymbol{\cdot}\boldsymbol{H}_{0}dV\Bigr{)}=-\int\frac{\partial v_{i}}{\partial x_{k}}\sigma^{v}_{ik}dV~{}.$$

$$-\nabla \tilde{p}+\eta \Delta \mathbf{v}=0,$$

$$(-\tilde{p}\delta_{ik}+\sigma^{v(i)}_{ik})n_{k}-(-p\delta_{ik}+\sigma^{v(e)}_{ik})n_{k}=\frac{\mu-1}{8\pi}\boldsymbol{H}^{(i)2}n_{i}+2\pi M^{2}_{n}n_{i}-\gamma\Bigl{(}\frac{1}{R_{1}}+\frac{1}{R_{2}}\Bigr{)}n_{i}~{}.$$

$$\boldsymbol{f}=\frac{\mu(\mu-1)}{8\pi}H^{(i)2}_{n}\boldsymbol{n}+\frac{\mu-1}{8\pi}H^{(i)2}_{t}\boldsymbol{n}-\gamma\Bigl{(}\frac{1}{R_{1}}+\frac{1}{R_{2}}\Bigr{)}\boldsymbol{n}~{}.$$

$$v_{i}(\boldsymbol{x})=\frac{1}{8\pi\eta}\int\Bigl{(}\frac{\delta_{ik}}{|\boldsymbol{x}-\boldsymbol{y}|}+\frac{(x_{i}-y_{i})(x_{k}-y_{k})}{|\boldsymbol{x}-\boldsymbol{y}|^{3}}\Bigr{)}f_{k}(\boldsymbol{y})dS_{\boldsymbol{y}}~{}.$$

$$\frac{E}{2\pi\gamma R^{2}_{0}}=K^{2/3}\Bigl{(}1+\frac{1}{2K^{2}\sqrt{1-1/K^{2}}}\log{\frac{1+\sqrt{1-1/K^{2}}}{1-\sqrt{1-1/K^{2}}}}\Bigr{)}-\frac{H^{2}_{0}R_{0}}{12\pi\gamma}\frac{1}{1/(\mu-1)+n_{1}}~{},$$

$$n_{3}=\frac{K^{2}}{(K^{2}-1)^{3/2}}\Bigl{(}(K^{2}-1)^{1/2}-\arctan{\sqrt{K^{2}-1}}\Bigr{)}~{}.$$

$$\displaystyle\frac{H^{2}_{0}R_{0}}{\gamma}=\pi\Bigl{[}\sqrt{K^{2}-1}(\Bigl{(}-1+\frac{2(K^{2}-1)}{\mu-1}\Bigr{)}+K^{2}\arctan{\sqrt{K^{2}-1}}\Bigr{]}^{2}$$

$$\displaystyle\times\frac{\Bigl{(}2\sqrt{K^{2}-1}K(1+2K^{2})+(1-4K^{2})\log{\frac{K+\sqrt{K^{2}-1}}{K-\sqrt{K^{2}-1}}}\Bigr{)}}{K^{7/3}(K^{2}-1)^{2}(-3\sqrt{K^{2}-1}+(2+K^{2})\arctan{\sqrt{K^{2}-1}})}~{}.$$

$$\frac{E}{2\pi\gamma R^{2}_{0}}=\Bigl{(}\frac{c^{2}}{ab}\Bigr{)}^{2/3}\Bigl{(}1+\frac{b}{\sqrt{a^{2}-c^{2}}}\Bigl{(}\rm{E}(k,m)\Bigl{(}\frac{a^{2}}{c^{2}}-1\Bigr{)}+\rm{F}(k,m)\Bigr{)}\Bigr{)}-\frac{1}{6}Bm\chi\Bigl{(}\frac{1}{1+4\pi\chi n_{x}}+\frac{1}{1+4\pi\chi n_{y}}\Bigr{)}~{},$$

$$n_z(c/a,c/b)=\frac{1}{2}{\int }_0^{\infty }\frac{du}{(u+1{)}^{3/2}((c^{2}u/a^2+1)(c^{2}u/b^2+1){)}^{1/2}};$$

$$n_x(c/a,c/b)=n_z(a/c,a/b);n_y(c/a,c/b)=n_z(b/a,b/c).$$

$$\int_{\partial\Omega}\frac{\partial\psi}{\partial n_{\boldsymbol{x}}}\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}dS_{\boldsymbol{x}}-\int_{\partial\Omega}\psi(\boldsymbol{x})\frac{\partial}{\partial n_{\boldsymbol{x}}}\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}dS_{\boldsymbol{x}}=\left\{\begin{array}[]{ll}4\pi\psi(\boldsymbol{y}){}&if~{}\boldsymbol{y}\in\Omega\\ 2\pi\psi(\boldsymbol{y}){}&if~{}\boldsymbol{y}\in\partial\Omega\\ 0{}&if~{}\boldsymbol{y}\ni\Omega\end{array}\right.$$

$$\psi (\mathbf{y})=\frac{2{\mathbf{H}}_0\cdot \mathbf{y}}{\mu +1}-\frac{1}{2\pi }\frac{\mu -1}{\mu +1}\text{vbox}-{\int }_{\partial \Omega }\psi (\mathbf{x})\frac{\partial }{\partial n_{\mathbf{x}}}\frac{1}{|\mathbf{x}-\mathbf{y}|}d{S}_{\mathbf{x}},$$

$$\mathbf{B}\cdot {\mathbf{n}}_{\mathbf{y}}={\mathbf{B}}_0\cdot {\mathbf{n}}_{\mathbf{y}}+{\mathbf{n}}_{\mathbf{y}}\cdot {\int }_{\partial \Omega }\mathbf{K}\times {\nabla }_{\mathbf{x}}\frac{1}{|\mathbf{x}-\mathbf{y}|}d{S}_{\mathbf{x}}={\mathbf{B}}_0\cdot {\mathbf{n}}_{\mathbf{y}}$$

$$-\frac{\mu -1}{4\pi }{\mathbf{n}}_{\mathbf{y}}\cdot {\int }_{\partial \Omega }({\mathbf{n}}_{\mathbf{x}}\times \mathbf{H})\times {\nabla }_{\mathbf{x}}\frac{1}{|\mathbf{x}-\mathbf{y}|}d{S}_{\mathbf{x}},$$

$$\psi (\mathbf{y})\left(1-\frac{\mu -1}{\mu +1}\right)=\frac{2{\mathbf{H}}_0\cdot \mathbf{y}}{\mu +1}-\frac{1}{2\pi }\frac{\mu -1}{\mu +1}{\int }_{\partial \Omega }\left[\psi (\mathbf{x})-\psi (\mathbf{y})\right]\frac{\partial }{\partial n_{\mathbf{x}}}\frac{1}{|\mathbf{x}-\mathbf{y}|}d{S}_{\mathbf{x}},$$

$${\psi }_j\left(1-\frac{\mu -1}{\mu +1}\right)=\frac{2{\mathbf{H}}_0\cdot {\mathbf{x}}_0}{\mu +1}-\frac{1}{2\pi }\frac{\mu -1}{\mu +1}\sum _{i\ne j}[{\psi }_i-{\psi }_j]\frac{({\mathbf{x}}_j-{\mathbf{x}}_i)\cdot {\mathbf{n}}_j}{|{\mathbf{x}}_i-{\mathbf{x}}_j{|}^3}{S}_j,$$

Full text

1 Introduction

The behaviour of droplets under the action of magnetic fields of different configurations is an important issue in many fields, including microfluidics (Seeman et al., 2012), mechanics of tissues (Douezan et al., 2011; Frasca et al., 2014), studies of dynamic self-assembly (Timonen et al., 2013) and many others. After the successful synthesis of magnetic fluids in the late sixties by (Rosensweig, 1985) an exciting story about droplets of magnetic fluid began. At first the elongation of the droplets in an external field was observed and explained by Arhipenko et al. (1978). The characteristic hysteresis at a sufficiently high magnetic susceptibility of the magnetic fluid was described by Bacri & Salin (1982). This problem is quite similar to the deformation of a conducting liquid droplet studied by Taylor (1964), which also described the so-called Taylor cone and the jet formation from the tip of droplet. Recently, progress has been made in studying the dynamics of magnetic drop deformation (Rowghanian et al., 2016).

Phase separation of magnetic colloids enabled researchers to obtain magnetic micro-droplets with a very high magnetic susceptibility allowing the observation of different new phenomena. One particularly interesting case is the behaviour of a magnetic micro-droplet under the action of a rotating magnetic field (Bacri et al., 1994). At low and intermediate rotating field frequencies synchronous and asynchronous rotation of the droplet with respect to the rotating field was observed by Sandre et al. (1999). Depending on the magnetic field strength and frequency, various shapes were seen (Janiaud et al., 2000). A family of stationary droplet shapes is observed in the present work as shown in figure 1. They are formed under the action of a high frequency rotating field. A detailed description of the observed sequence of bifurcations is given below.

In a high frequency field the characteristic time of the shape deformation is much larger than the field period, therefore it is possible to average the magnetic field energy with respect to the field period. As a result the droplet shape is a figure of equilibrium and is determined by the balance of the ponderomotive forces of the self-magnetic field and the capillary forces (Bacri et al., 1994). On this basis in (Bacri et al., 1994; Cebers & Lacis, 1995) the re-entrant transition for the sequence of shapes oblate-prolate-oblate was described by applying the virial method and assuming that prolate ellipsoid of revolution arises at the oblate-prolate-oblate bifurcations. Later these results were extended in (Morozov & Lebedev, 2000; Morozov et al., 2002) considering the general ellipsoid with three non-equal axes at the oblate-prolate-oblate bifurcations. An essential assumption of these models is the ellipsoidal shape of the droplet.

Our aim is to answer the following question. Is it possible to describe the set of arising bifurcations by a simple model of magnetic fluid (Rosensweig, 1985)? For that we develop a numerical algorithm for a droplet in a rotating field which is free from the limitations of the ellipsoidal approximation. We compare the numerical results with experimental ones (including those shown in figure 1) and also with analytical results, where it is possible. We organize the paper as follows. In § 2 we formulate the model and give the available theoretical results. The numerical algorithm is based on a boundary integral equation technique. It is described with its tests in § 3. The implementation written in the Julia programming language can be found in Erdmanis (2016). In § 4 the numerical analysis of the magnetic droplet shape bifurcations is carried out. The presentation of the experimental results and their comparison with numerical and theoretical results are carried out in § 5, before a short discussion and conclusion.

2 Model and analytical results

2.1 Equation of motion and energy functional

Equations of motion of the magnetic fluid reads (Rosensweig, 1985; Blums et al., 1997)

[TABLE]

where $\sigma^{v}_{ik}=\eta\left(\frac{\partial v_{i}}{\partial x_{k}}+\frac{\partial v_{k}}{\partial x_{i}}\right)$ is the viscous stress tensor, $\eta$ is the viscosity, $p$ is the pressure, \boldsymbol{e}_{i}\boldsymbol{\cdot}\boldsymbol{T}_{k}=\frac{1}{4\pi}\Bigl{(}H_{i}B_{k}-\frac{1}{2}H^{2}\delta_{ik}\Bigr{)} is the Maxwell stress tensor and $\frac{\partial\boldsymbol{T}_{k}}{\partial x_{k}}=(\boldsymbol{M}\boldsymbol{\cdot}\nabla)\boldsymbol{H}$. According to the equations of magnetostatics $\nabla\boldsymbol{\cdot}\boldsymbol{B}=0;~{}\boldsymbol{B}=\boldsymbol{H}+4\pi\boldsymbol{M}$ and $\nabla\times\boldsymbol{H}=4\pi\boldsymbol{J}/c$ ($\boldsymbol{J}$ is the current density in external coils creating the magnetic field and equal to zero throughout the volume of the magnetic fluid.) In the framework of the linear response, magnetization $\boldsymbol{M}$ expresses as $\boldsymbol{M}=\chi\boldsymbol{H}$, with magnetic susceptibility $\chi$ such as $\mu=1+4\pi\chi=const$, $\mu$ being the magnetic permeability. The equilibrium shapes of the magnetic fluid droplet are usually analysed considering the energy functional (Bacri & Salin, 1982) (see (Landau & Lifshitz, 1960) for the derivation of the magnetic field contribution)

[TABLE]

where $\boldsymbol{H}_{0}$ is the field strength of external homogeneous magnetic field, $S$ is the surface area of a droplet and $\gamma$ is the surface tension (supposed here independent on $H$). Although this approach, as it should be expected, is reasonable, for the completeness of further analysis we illustrate below that it directly follows from the dynamic equations (1) and the boundary conditions describing the evolution of the shapes of the droplet.

Considering the Lagrange variation of some variable $f$ ($\boldsymbol{\xi}(\boldsymbol{x})$ is the displacement of material element with a radius vector $\boldsymbol{x}$)

[TABLE]

and its derivative

[TABLE]

while the current of external source is fixed according to (4) we have ($\psi$ being the potential of the magnetostatic field defined by $\boldsymbol{H}=\boldsymbol{\nabla}\psi$)

[TABLE]

For simplicity we further assume that the fluid is incompressible ($\nabla\boldsymbol{\cdot}\boldsymbol{\xi}=0$). According to Eq.(5) we have

[TABLE]

Integrating the relation Eq.(6) over the entire volume and accounting for the boundary conditions on the surface of a droplet $\delta_{L}\psi^{i}=\delta_{L}\psi^{e}$ and $B_{n}^{i}=B_{n}^{e}$ we obtain (by superscripts i and e we denote the values on the inside and outside boundaries of the droplet respectively)

[TABLE]

Taking into account the identity valid for an incompressible fluid

[TABLE]

the relation (7) may be rewritten as follows

[TABLE]

Using $\xi_{i}=v_{i}\delta t$, integrating the dot product of Eq.(1) with the fluid velocity field and accounting for the incompressibility we obtain ($\sigma_{ik}=-p\delta_{ik}+\sigma^{v}_{ik}$)

[TABLE]

The balance of forces on the boundary of the droplet gives

[TABLE]

where $1/R_{1,2}$ are principal curvatures. As a result we conclude that the variation of the sum of the surface energy $\gamma S$ and the thermodynamic potential at fixed current in the external coils $\tilde{F}=-\frac{1}{8\pi}\int\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{H}dV$ equals the amount of the dissipated energy per unit time due to the viscosity of the fluid

[TABLE]

Instead of the thermodynamic potential $\tilde{F}$, the potential obtained by adding a constant value $\frac{1}{8\pi}\int H_{0}^{2}dV$ is usually considered (Landau & Lifshitz, 1960), thus leading to

[TABLE]

and

[TABLE]

The energy $\gamma S-\frac{1}{2}\int\boldsymbol{M}\boldsymbol{\cdot}\boldsymbol{H}_{0}dV$, which is usually minimized to describe the equilibrium shapes of the magnetic droplets (Bacri & Salin, 1982; Blums et al., 1997; Afkami et al., 2010), is in fact decreasing with time due to the viscous dissipation and reaching the minimal value at equilibrium (Bacri & Salin, 1983).

The arguments given above show that the figures of equilibrium of the magnetic fluid droplets may be obtained by numerical simulation of the evolution of shapes of a viscous droplet suspended in a surrounding viscous fluid, as described by Eq. (1) and the boundary conditions Eq.(11). For simplicity we will here consider the case of both the magnetic and the surrounding fluids having equal viscosities.

2.2 Dynamics of droplet shape in an applied field

We now formulate the set of equations which are solved numerically for the simulation of the droplet shape dynamics in a rotating field. The equation of motion (1) in the case of a constant magnetic susceptibility may be rewritten in the form of the Stokes equation

[TABLE]

where $\tilde{p}=p-\frac{(\mu-1)\boldsymbol{H}^{2}}{8\pi}$. The boundary condition (11) accounting for $T^{(i)}_{nn}-T^{(e)}_{nn}=-2\pi M_{n}^{2}$ and $T^{(i)}_{tn}=T^{(e)}_{tn}$ may be rewritten as follows

[TABLE]

From the relation (16) we see that the fluid motion is determined by the action of the normal force on the boundary of the droplet

[TABLE]

Its velocity for equal drop and surrounding fluid viscosities is given as an integral over surface forces as a special case of the Stokes BIE Pozikridis (2003)

[TABLE]

The normal component of the liquid velocity (18) on the interface gives its position in the next time step. We see that the dynamics is known if the normal $H^{(i)}_{n}$ and tangential $H^{(i)}_{t}$ components of the magnetic field strength on the internal side of the droplet boundary are found. These components are given by the solution of equations of magnetostatics $\nabla\boldsymbol{\cdot}\boldsymbol{H}=0;~{}\nabla\times\boldsymbol{H}=0$ with boundary conditions $H^{(i)}_{t}=H^{(e)}_{t}$, $\mu H^{(i)}_{n}=H^{(e)}_{n}$ and $\boldsymbol{H}(\boldsymbol{x})\rightarrow\boldsymbol{H}_{0}$ at $|\boldsymbol{x}|\rightarrow\infty$. This problem of magnetostatics may be efficiently solved by applying a boundary integral equation technique which we describe in the next section.

2.3 Ellipsoidal figures of equilibrium

The result (14) allows us to find some of the simplest figures of equilibrium of the magnetic fluid droplet in a high frequency rotating magnetic field. If the period of the rotating field is much smaller than the characteristic time $\gamma/\eta R_{0}$ of the droplet shape relaxation, where $R_{0}$ is the radius of initial spherical droplet, we may average the magnetic energy with respect to the period of the rotating field and consider the minimum of the energy functional $\gamma S-\frac{1}{2}\int\boldsymbol{M}\boldsymbol{\cdot}\boldsymbol{H}_{0}dV$. As an approximation we assume that the droplet has the shape of a general ellipsoid with the semi-axes $a,b,c$ and its interface is determined by the equation $x^{2}/a^{2}+y^{2}/b^{2}+z^{2}/c^{2}=1$.

As the droplet always has an axial symmetry at a small field strength we consider it here as an oblate ellipsoid of revolution with the semi-axes $a=b>c$. The demagnetizing field coefficients are $n_{1}=n_{2}$ in the plane $(x,y)$ of the rotating field and $n_{3}$ in the perpendicular direction. As a result the energy of the droplet reads ($K=a/c>1$)

[TABLE]

where $n_{1}=(1-n_{3})/2$ and

[TABLE]

We obtain the axis ratio $K$ for the equilibrium shape by minimizing $E$ with respect to $K$. After a long but straightforward calculation the following equation is obtained

[TABLE]

Relation (2.21) was derived in (Morozov & Lebedev, 2000) in a more compact form using the demagnetizing field coefficients of an ellipsoid and its adjoint.

Figure 2 shows for several values of the magnetic permeability the analytic dependence of the ratio $c/R_{0}=K^{-2/3}$111We have chosen to plot $c/R_{0}$ instead of $K$ since it remains a well defined quantity when the droplet takes a tri-axial shape. on the magnetic Bond number $Bm=H^{2}_{0}R_{0}/\gamma$ ($R^{3}_{0}=abc$). It shows that if the magnetic permeability is large enough, there is a hysteresis of $c/R_{0}$ as a function of the field. At infinite magnetic permeability an equilibrium shape, in the form of an oblate ellipsoid of revolution, only exists for magnetic fields below a critical field value.

In this respect the behavior is similar to that of a conducting droplet in an electrical field. An equilibrium shape, elongated along the field, only exists if the electrical field strength is below a critical value (Taylor, 1964). Above it there are jets appearing from the tips of the droplet.

Since in our case a fingering is observed experimentally (see figure 1 and Bacri et al. (1994)) on the perimeter of the disc-like droplet then it may be related to the threshold of existence of oblate equilibrium shape for large values of the magnetic permeability.

By Bacri et al. (1994) (see Cebers & Lacis (1995) for more details) it was found during the analysis of the experimental results that the oblate shape loses its stability with respect to the symmetry destroying perturbations, the long axis of the ellipsoid of revolution remaining in the plane of the rotating field. The main conclusions included a critical value of the magnetic permeability for this instability and the re-entrant character of this oblate-prolate transition. This approach was further developed by Morozov & Lebedev (2000) where a general ellipsoid with the semi-axes ($a>b>c$) for symmetry destroying perturbations was considered. The main features of the oblate-prolate transition, apart from the prediction of hysteresis in the oblate-prolate and the prolate-oblate transitions, qualitatively remained unchanged. Here we are using these results (see figure 3) for the comparison with the numerical simulation of the three-dimensional droplet and experiments (see ahead in figure 10) . For completeness we give the main points of the analysis following to Morozov & Lebedev (2000) approach. The energy of the droplet reads

[TABLE]

where $k=\arcsin{\sqrt{1-c^{2}/a^{2}}};m=(1-(c/b)^{2})/(1-(c/a)^{2})$ and $\rm{E},\rm{F}$ are the elliptic integrals. We are using relations (Beleggia et al., 2006) that express the demagnetizing field coefficients $n_{x},n_{y}$ in terms of the demagnetizing field coefficient $n_{z}$:

[TABLE]

[TABLE]

The equilibrium shapes are found by minimizing the energy of the droplet with respect to $c/a$ and $c/b$. The bifurcation diagram in coordinates $H^{2}_{0}R_{0}/\gamma$ and $b/a$ for several $\mu$ values ($\mu=6,10,17$) is given in figure 3. A hysteresis in the oblate-prolate-oblate transition at large values of the magnetic field strength emerges for large values of $\mu$ which is in accord with the results obtained by Morozov & Lebedev (2000). For small field values the hysteresis is small. The calculations are in agreement with (Morozov & Lebedev, 2000) and show that the transition from oblate to a triaxial ellipsoid exists above some critical value of the magnetic susceptibility (in the interval of $\mu$ $[5;6]$).

3 Numerical algorithm

3.1 Boundary integral equations (BIE) for magnetic field

Since $\nabla\times\boldsymbol{H}=0$, we recall that the magnetic field is potential $\boldsymbol{H}=\nabla\psi$ and the potential due to $\nabla\boldsymbol{\cdot}\boldsymbol{H}=0$ satisfies the Laplace equation $\Delta\psi=0$. This equation can also be written in a boundary integral form (Pozikridis, 2003)

[TABLE]

where $\Omega$ is an arbitrary region where $\Delta\psi=0$ and $\partial\Omega$ denotes the boundary of the fluid drop. If $\boldsymbol{y}\in\partial\Omega$ the second integral from the left side must be interpreted in the Cauchy principal value sense which we denote here as \mathchoice{{\vbox{\hbox{\textstyle\boldsymbol{-}}}\kern-3.75pt}}{{\vbox{\hbox{\scriptstyle\boldsymbol{-}}}\kern-3.25pt}}{{\vbox{\hbox{\scriptscriptstyle\boldsymbol{-}}}\kern-2.29166pt}}{{\vbox{\hbox{\scriptscriptstyle\boldsymbol{-}}}\kern-1.875pt}}\!\int.

Applying the boundary conditions $\psi^{(i)}=\psi^{(e)}$ and $\mu\frac{\partial\psi^{(i)}}{\partial n}=\frac{\partial\psi^{(e)}}{\partial n}$ for the boundary of fluid drop $\partial\Omega$ we obtain the following boundary integral equation (BIE) for the magnetic potential (see appendix A)

[TABLE]

where $\boldsymbol{y}\in\partial\Omega$. After solving this equation for $\psi$ we can use numerical differentiation to obtain the tangential components of the magnetic field $(\boldsymbol{I}-\boldsymbol{n}\otimes\boldsymbol{n})\boldsymbol{H}$.

Numerous boundary integral equations exist for the normal component of the magnetic field. One possibility is to use (25) for the normal component of the field with a known $\psi$. However this gives an integral equation of the first kind and numerical experiments show a high sensitivity to the precision of the $\psi$ values. Another way is to formulate the boundary integral equation for the normal field as previously done for fluid droplets suspended in an electric field by Baygents et al. (1998) (see also appendix § A for derivation). However, it requires the integration of a strong singularity which is hard to do numerically in the general 3D case. Eventually we settled on interpreting a perturbation of the magnetic field from the droplet as coming from the surface currents on its interface $4\pi\boldsymbol{K}/c=\boldsymbol{n}\times(\boldsymbol{B}^{(e)}-\boldsymbol{B}^{(i)})=-(\mu-1)\boldsymbol{n}\times\boldsymbol{H}$.

The magnetic field outside the fluid droplet is given by the Biot-Savart integral $\boldsymbol{B}=\boldsymbol{B}_{0}+c^{-1}\int_{\partial\Omega}\boldsymbol{K}\times\nabla_{\boldsymbol{x}}\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}dS_{\boldsymbol{x}}$, which for the normal component of the magnetic field in the presence of an external field $\boldsymbol{H}_{0}$ gives

[TABLE]

which we also derive in appendix § A without the use of the surface current interpretation. Since the integral on the right side of relation (3.3) only has a cross product of the field and the normal, it allows us to calculate the normal field component from the known tangential components. The numerical algorithm for this calculation is described in § 3.2.

3.2 Numerical implementation

Solving the BIE for the magnetic potential in the form (26) or integrating (3.3) would require singular integral quadrature which makes these equations hard to implement numerically. We overcome this difficulty using a singularity subtraction technique, which is a common practise in the application of BIE (Keaveny & Shelley, 2011). Specifically, for the magnetic potential we have an identity \mathchoice{{\vbox{\hbox{\textstyle\boldsymbol{-}}}\kern-3.75pt}}{{\vbox{\hbox{\scriptstyle\boldsymbol{-}}}\kern-3.25pt}}{{\vbox{\hbox{\scriptscriptstyle\boldsymbol{-}}}\kern-2.29166pt}}{{\vbox{\hbox{\scriptscriptstyle\boldsymbol{-}}}\kern-1.875pt}}\!\int\frac{\partial}{\partial n_{\boldsymbol{x}}}\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}dS_{\boldsymbol{x}}=-2\pi (see appendix § B) which we use to rewrite (26) in the following form

[TABLE]

where \mathchoice{{\vbox{\hbox{\textstyle\boldsymbol{-}}}\kern-3.75pt}}{{\vbox{\hbox{\scriptstyle\boldsymbol{-}}}\kern-3.25pt}}{{\vbox{\hbox{\scriptscriptstyle\boldsymbol{-}}}\kern-2.29166pt}}{{\vbox{\hbox{\scriptscriptstyle\boldsymbol{-}}}\kern-1.875pt}}\!\int\rightarrow\int since we reduced the singularity by one order. The discretization of the drop surface into plane triangles enables us to solve it for arbitrary geometries.

The simplest quadrature for triangles $\int_{\Delta}fdS=(f_{1}+f_{2}+f_{3})S_{\Delta}/3$ (trapezoidal quadrature) is used and it allows us to rearrange the integration efficiently. Denoting by $S_{j}$ the $1/3$ of the sum of the area of all triangles with a common vertex $j$ we can rewrite (28) in a discrete form

[TABLE]

which we solve numerically with the LAPACK library. Using a triangulation of the ellipsoid (for which the normal vectors and the solution of the magnetostatic problem are known) we see that the solution is only moderately sensitive to the precision of the normal vectors. Therefore, we use a simple and robust method for calculating the vertex normal vectors (Jin et al., 2005) which is

[TABLE]

where the summation is over the triangles which have a vertex $j$. In relation (3.6) $\boldsymbol{n}_{i}$ is the normal to the plane triangle and $\alpha_{i}$ is the angle of the triangle at the vertex $j$.

In order to obtain the tangential components of the magnetic field we use the numerical differentiation of the obtained potential. For each vertex $j$ we write a system of equations

[TABLE]

where $i$ are the neighbour vertices. This system is overdetermined for the three components of $\nabla\psi$ therefore the linear least squares method is applied. Tangential field components are then projected out of the estimated gradient $\rm{P}\boldsymbol{H}_{j}=(\boldsymbol{I}-\boldsymbol{n}_{j}\otimes\boldsymbol{n}_{j})(\nabla\psi)_{j}$, where $P=(\boldsymbol{I}-\boldsymbol{n}_{j}\otimes\boldsymbol{n}_{j})$ is the projection operator. The known tangential field components according to the relation (3.3) give us the normal component of the field as a Biot-Savart quadrature.

Since the quadrature in (3.3) is strongly singular it is inefficient to do it numerically in this form. The general procedure for subtracting the singularity in the Biot-Savart integral was given by Pozrikidis (2000), but it requires calculating the curvature by adding another issue of complexity. Instead of that we have developed a regularisation technique which does not require a curvature calculation.

In the appendix § B we prove the identity

[TABLE]

where $\rm{P}H_{\boldsymbol{y}}$ is a constant vector. This identity enables us to evaluate the strongly singular quadrature on the left side of relation (3.8) with weakly singular quadratures on its right side. Using this formula for (3.3) gives us weakly singular form of Biot-Savart quadrature

[TABLE]

We integrate both integrals of (33) with the trapezoidal rule for non-singular integrals and use a polar transformation for weakly singular integrals (Pozikridis, 2003). Specifically, when the singularity is on $\boldsymbol{x}_{1}$ we have the quadrature formula

[TABLE]

where we use a Gauss-Legendre quadrature with $10$ points for the $\rho$ and $\chi$ variables.

Besides the magnetic field another important quantity for the surface force calculation is the curvature as can be seen in (17). However, it turns out that the explicit mean curvature calculation can be avoided since we are interested in an integral containing a curvature over a closed surface. In (Pozrikidis, 2001), (Zinchenko et al., 1997) and appendix § C the following identity is derived for the value of the integral at a surface point with the radius vector $\boldsymbol{y}$ ($\boldsymbol{r}=\boldsymbol{x}-\boldsymbol{y}$)

[TABLE]

Putting this formula in (17), (18) and using the identity $\int\left(\frac{\boldsymbol{I}}{r}+\frac{\boldsymbol{r}\otimes\boldsymbol{r}}{r^{3}}\right)\boldsymbol{\boldsymbol{\cdot}}\boldsymbol{n}_{\boldsymbol{x}}dS_{\boldsymbol{x}}=0$ for the singularity subtraction we get the normal component of the velocity

[TABLE]

[TABLE]

Using the trapezoidal integration formula for all triangular elements for both integrals we obtain the normal component of the velocity which we use to proceed to the next time step.

3.3 Averaged equations for high frequency

As noted in section § 1 we are interested in the behaviour of the droplets in a rotating field with a high frequency. In this limit we can neglect the droplet shape variation during the field period. At equal droplet and surrounding fluid viscosities the numerical simulations have shown that the high frequency behaviour establishes at $\omega\gg 10R_{0}\eta/\gamma$.

To derive the corresponding high frequency equations let us set, as previously, the magnetic field rotation plane as the $xy$ plane

[TABLE]

Because of the superposition principle we can write a solution of the magnetostatic equations

[TABLE]

where we set $\boldsymbol{H}_{x}=\boldsymbol{H}(\boldsymbol{x},\boldsymbol{e}_{x})$ and $\boldsymbol{H}_{y}=\boldsymbol{H}(\boldsymbol{x},\boldsymbol{e}_{y})$. With this solution we can now write the normal component of the instantaneous magnetic surface force (36) as

[TABLE]

Averaging this force over a single period gives us

[TABLE]

which can be inserted in (3.13) to obtain the average normal velocity in high frequency rotating field.

3.4 Mesh maintenance and generation

Within a small number of time steps, at which the surface markers were advanced, an uneven triangulation is produced. It makes the representation of the geometry unsatisfactory due to degenerate triangular elements affecting the precision of the force and velocity calculation. To overcome this difficulty mesh relaxation and stabilization methods have been used for viscous drop simulations in Cristini et al. (2001) and Zinchenko et al. (1999). Here we use a general purpose mesh stabilization/relaxation algorithm by Brochu & Bridson (2009), which is also available as an open source library named ElTopo.

This algorithm collapses edges which are smaller than $l_{min}$, splits those which are larger than $l_{max}$ and flips edges if after an operation they become shorter. The algorithm also moves each vertex to its average of neighbour vertices while the movement of vertex is restricted to its estimated normal plane.

To ensure that these operations do preserve sharp features in places with a high curvature they are allowed only if the volume changes introduced by any operation are smaller than $\gamma_{vol}$. Also, to avoid generation of new degenerate triangles, an operation is not allowed if it introduces triangles with an angle either smaller than $\alpha_{min}$ or larger than $\alpha_{max}$. In our studies we used the following set of mesh configuration maintenance parameters

[TABLE]

where we set $\xi$ as the average edge length of the initial mesh.

We used a simple Euler method for advancing markers to new positions and used a sufficiently small time step so that the volume changes are less than $1\%$ when the equilibrium is obtained. We observed that large time steps could also induce a random noise and that the main limitation for the step-size comes from the higher chance of getting a degenerate mesh in the middle of a simulation.

For starting the simulation we need an initial triangulation. We used the Distmesh library developed by Persson & Strang (2004) for generating the initial triangulation from the analytical signed distance functions taking only the average edge length $\xi$ as an argument. With this algorithm we successfully generated high quality meshes for spheres, ellipsoids and star-like shapes.

4 Numerical simulation results

4.1 Worm-like equilibrium figures

The compression of the droplet perpendicularly to the plane of the rotating field was used as a test for the developed algorithm along with smaller tests for droplet extension in a constant field, relaxation dynamics and magnetic field calculation for the ellipsoid.

The procedure for finding the figures of equilibrium is as follows. Starting with initial triangulation we let the simulation algorithm relax the shape until the following conditions are met

[TABLE]

where $\Delta t$ is the time-step, $v_{i}$ is the maximum velocity at the step $i$ and $v_{start}$ is the velocity when the relaxation of triangulation begins. First condition ensures that the velocities start to decay exponentially which is expected when the surface is close to an equilibrium. The second condition assuming that equilibrium is approached exponentially requires that the maximum deviation from true equilibrium (when $t\to\infty$), for any point on the surface, is smaller than $1/100$ of $R_{0}$.

For obtaining figures of equilibrium we used two methods. In the simplest one (ellipsoid relaxation) we started from an ellipsoid with semi-axes calculated by the energy (22) minimization. In the second more complex one (quasi-static simulation) we start from triangulation of a sphere and increase the magnetic field quasi-statically ($\Delta Bm=1$) letting the triangulation relax at each step until condition (42) is met.

We induced a spontaneous symmetry breaking of an axial figure, in a quasi-static simulation, by stretching the triangulation in $x$ axis direction and compressing in the $y$ and $z$ axes directions to conserve the volume at the beginning of each quasi-step. When the difference between the lengths in the $x$ and $y$ directions became larger than $0.01R_{0}$ we disabled this part of the algorithm. By disabling we kept triangulation more regular which, as we observed, slightly affects simulated equilibrium figures. We observed that more time steps are necessary for a quasi-step when the transition to a tri-axial figure is approached. The figures of equilibrium obtained with this algorithm at $\mu=10$ for different magnetic Bond numbers are shown in figure 4.

However, some differences are observed depending on how the axial symmetry destroying perturbations are imposed - ellipsoid relaxation or quasi-static simulation. Equilibrium shapes calculated by ellipsoid relaxation method are in excellent agreement with analytical results and also with experiments (as we will see in §  5). The numerical experiments show that the difference between the figures of equilibrium obtained by ellipsoidal relaxation and by quasi-static simulation is due to an overly coarse triangulation and the quality of the mesh.

4.2 Star-like equilibrium figures

Even if we do not manage to induce the star-like instability with numerical simulations directly from a worm-like figure, we are able to make a simulation of a star-like equilibrium figure, if we initially start with a mesh which has a broken symmetry at the starting point. As a seed we use the following surface

[TABLE]

which we meshed with the Distmesh library (Persson & Strang, 2004). We instantly set a magnetic field and depending on its strength we obtain the star-like equilibrium figures (see figure 5).

For small magnetic Bond numbers $Bm<19$ the perturbation decays and we obtain an axisymmetric oblate shape. At $Bm=30$ the 6-star perturbation still decays, but due to some imperfections in the initial mesh we obtain a triangular shape. Starting from magnetic Bond number $40$ we see that the 6-star perturbation grows for larger magnetic Bond numbers. Qualitatively the simulated star-like figures for different magnetic Bond numbers differ by their sizes in the $(x,y)$ plane but keep the length of the peaks approximately constant.

With higher magnetic Bond numbers up to $100$ we were not able to show that the initial higher order perturbation decays. Instead some higher order perturbations determine the shape of the droplet as we expected from the experiment. However, we did calculate energies for the stars. Within error bars they agree with the quasi-static simulations for $Bm\leq 50$. However at $Bm=70$ we get an energy $E/4\pi R_{0}^{2}\gamma=-3.15$ for the star, where for the worm-like shape we obtain $E/4\pi R_{0}^{2}\gamma=-2.98$.222 For numerical energy calculation we used formula (4.4). We tested this formula with ellipsoidal meshes. Comparing with analytical formula (2.19) the error for the magnetic part was found to be about $2\%$ for a mesh with $3000$ triangles. Thus, the star-like shape has a smaller energy than the worm-like shape for magnetic Bond numbers $Bm\geq 70$. We can also conclude that by increasing the magnetic Bond number, the number of stable shapes increases while in the experiment we observe only two, when a hysteresis at the re-entrant transition happens.

4.3 Numerical tests and results

For the ratio $c/R_{0}$, which characterizes the compression of the droplet in a rotating field and is calculated according to relation (2.21), we see a qualitatively good agreement between the numerical results and the analytic ellipsoidal approximation in a wide range of $\mu$ values as shown in figure 2. Specifically when $\mu=4$ the axially symmetric state of the droplet remains stable - the result which was predicted analytically for the droplet with an oblate ellipsoidal shape.

On the other hand while the droplet keeps an ellipsoidal shape we may use the formula (2.22) for calculating the axis ratio in the plane of rotating field as can be seen in figure 6. We observe that the axis ratio $b/a$ starts to deviate from $1$ exactly where the ellipsoidal approximation predicts a transition from an oblate ellipsoid to a tri-axial one. We also see that both simulations very well agree with the analytic ellipsoidal approximation as expected since droplet in simulation keeps a quasi-ellipsoidal shape. For large magnetic Bond numbers the quasi-static simulation becomes unreliable because of mesh quality loss and in figure 6 we can see its growing difference with ellipsoid relaxation simulation.

Figure 6 also presents two calculated (ellipsoid relaxation) equilibrium figures with ellipses drawn on the top of them to show better their deviations. We can see that figures of equilibrium starts to deviate significantly from the ellipsoidal shape for large magnetic fields, forming sharp tips as it is also observed in the experiments (figure 1). Unexpectedly the simulation closely follows the analytic ellipsoidal approximation even for large magnetic fields where the droplet presents tips. Thus we can use the ellipsoidal approximation quantitatively for adjusting the experimental data.

We also calculate energies of the equilibrium figures and compare them with the energies found analytically in the ellipsoidal approximation (22). This latter should be larger or equal due to constraints. The total instantaneous energy of the drop is the sum of the surface energy $E_{S}=\gamma S$ and the magnetization energy $E_{M}$ (13)

[TABLE]

Since we consider the high frequency limit of a rotating field, then according to the superposition principle (39), we write the average energy per period in the following way

[TABLE]

where $\psi_{x}(\boldsymbol{x})=\psi(\boldsymbol{e}_{x},\boldsymbol{x})$ is the potential for a unit field directed along the $x$ axis. Numerical integration with the trapezoidal rule (45) for each time-step showed that the energy decreases monotonously until it reaches its equilibrium value.

We plot equilibrium energies in figure 7 where we see that the energies of equilibrium shapes obtained in the quasi-static simulation and in the ellipsoidal approximation (22) are within errorbars equal to the analytically calculated energy. The main contribution to the error comes from its magnetic part. For mesh with $3000$ triangles it equals to $2~{}\%$.

As a last numerical result, we show here that the critical slowing down is a useful tool for the calculation of the critical magnetic field associated to the threshold of tri-axial ellipsoidal shape. By analogy with (Cebers & Zemitis, 1983) we can estimate the critical Bond number by calculating the decrement of the symmetry-breaking perturbations going to zero due to the critical slowing down at the bifurcation value of the magnetic Bond number. Using a simulation with $\mu=10$, which is close to the value for droplets used in experiments, we estimate the inverse of the logarithmic decrement $\tau$ just before we jump to the next quasi-step. This decrement is drawn in figure 8 in a normalized form. We clearly see that the drop shape with an axial symmetry becomes unstable around $Bm=19..20$ agreeing with result from ellipsoidal approximation and simulation (figure 6).

5 Experiments and their comparison with numerical simulation

5.1 Preparation of magnetic droplets

Magnetic droplets can be prepared in various ways. The most common method is to disperse a magnetic fluid in an immiscible liquid (Zakinyan & Dikansky, 2011). Recently an interesting system was used where drops in glycerin were formed by magnetic nanoparticles with a biocompatible polydimethylsiloxane (PDMS) coating without a carrying liquid (Afkami et al., 2010). However, here we use magnetic drops that are obtained by inducing a phase separation in a magnetic fluid via an increase of ionic strength. The water based magnetic fluid is prepared in the PHENIX laboratory, using the precipitation method (Massart, 1981), and has maghemite nanoparticles with a size $d=7.1$ nm and polydispersity $\sigma_{\rm PDI}=0.32$, volume fraction $\Phi=2.8$ %, volume susceptibility $\chi=0.016$ (CGS units) and saturation magnetization $M_{\rm sat}=8.4$ G at $H_{sat}=10$ kOe, as given by magnetization measurements. The colloidal particles are stabilized with citrate ions from trisodium citrate, which is added in excess, creating a large initial ionic strength $I_{\rm citr}=0.18$ M present in the magnetic fluid. By adding $0.05$ M of sodium chloride a phase separation is induced and droplets of concentrated phase are formed in the dilute phase. As the droplets have a larger concentration of magnetic nanoparticles than the surrounding liquid in coexistence and a low surface tension, droplets respond to an external magnetic field of low amplitude. A volume of $\approx 25$ $\mu$l of the phase separated fluid is put in quartz cuvettes of $0.1$ mm thickness (Helma analytics). The sample is then investigated using the experimental setup.

5.2 Experimental setup

The experimental setup is based on a ZEISS Axio Observer.D1 microscope equipped with a custom-built coil system. Coils are powered by two KEPCO BOP 20-10ML power supplies and controlled by NI DAQ system. The system has been calibrated with a gauss-meter, to link the magnetic field with current measurements. Images are recorded with AVT Guppy F-046B (8-bit grayscale, resolution $640\times 480$ px, framerate $3.75$ Hz). All devices are governed by a custom-built program in NI LabVIEW. The system enables us to apply constant and alternating fields in the image plane. A more detailed description of the experimental system and methods can be found in Kitenbergs (2015).

Recorded images are processed with MATLAB®. They are first thresholded to binary images for detecting individual objects. Then, to remove noise, dirt and other small drops, the largest object is selected and its perimeter is traced. The longest distance between two points of the perimeter defines $2a$, while $2b$ is measured perpendicular to the midpoint of $2a$.

5.3 Properties of magnetic droplets in constant field

One way to experimentally obtain the physical properties of a magnetic drop, namely surface tension $\gamma$, magnetic permeability $\mu$ and viscosity $\eta$, is to make an elongation-relaxation experiment with subsequent analysis, as shown in figure 9.

For the elongation part, the magnetic drop is gradually deformed under a stepwise increase of a homogeneous and static magnetic field, with its direction in the image plane (figure 9 (a), parallel to the long side of the image). After each step an image and current measurements are recorded once the steady state has been reached. Image and data processing allows to form a plot of the ratio $a/b$ as a function of magnetic field, as shown in figure 9 (b). The data are then adjusted with the equation from Bacri & Salin (1982) and Cebers (1985), assuming an ellipsoidal shape of revolution and a constant volume for the droplet:

[TABLE]

where $e$ is the ellipse eccentricity $(e=\sqrt{1-\frac{b^{2}}{a^{2}}})$ and $N$ its demagnetizing coefficient, which can be expressed for a prolate ellipsoid of revolution in the following way

[TABLE]

Equation (46) allows us to obtain the $\gamma$ and $\mu$ values.

For the relaxation part, the magnetic field is cut to zero after a sufficient elongation and the droplet is allowed to relax back to a spherical shape. Its relaxation is recorded over time (figure 9 (c)). Close to a spherical shape, the parameter $(a-b)/a$ follows an exponential decay with a characteristic time $\tau_{\eta}$ as illustrated by figure 9 (d). In these experiments $(a-b)/a<0.02$ is the resolution limit, $\tau_{\eta}$ is described with an equation from Dikansky et al. (1990):

[TABLE]

where $R_{0}$ is the initial drop radius, $\eta_{c}$ is the viscosity of the drop (concentrated phase) and $\eta_{d}$ is the viscosity of the surrounding liquid (diluted phase). $\eta_{c}$ can be found by taking $\eta_{d}$ as the viscosity of water ($\eta_{d}=0.01$ P) and using the surface tension $\gamma$ from the elongation part. The phase-separated magnetic fluid droplets used in these experiments have a radius $(R_{0}\pm\Delta R_{0})=(10.0\pm 0.5)$ $\mu$m and properties $(\gamma\pm\Delta\gamma)=(6.0\pm 0.5)\boldsymbol{\cdot}10^{-3}~{}\rm{erg}/\rm{cm}^{2}$, $(\mu\pm\Delta\mu)=(10.5\pm 0.5)$ and $(\eta_{c}\pm\Delta\eta_{c})=(10\pm 1)$ P. However, let us note that according to (Bashtovoi et al., 1987) equation (5.1) is only valid up to $a/b\approx 7$, thus hereafter we use calculated $\gamma$ deduced from the compression of the oblate shape.

5.4 Properties of magnetic droplets in rotating field

Magnetic drop deformation under a rotating field is also studied here, using a frequency $500$ Hz. It is notably larger than the frequency of the droplet breathing mode $\gamma/R_{0}\eta_{c}\approx 0.5$ Hz fulfilling the high frequency criterion. Magnetic field strength is increased and decreased stepwise, allowing to reach equilibrium after each step. Equilibrium shape is recorded as a function of the field, together with the characteristic time to reach equilibrium.

Both in numerical simulations and experiments we see that the oblate magnetic droplet becomes unstable with a respect to axial symmetry breaking perturbations at the critical magnetic Bond number. Above the threshold value of the magnetic Bond number it assumes a worm-like shape for small fields and a star-like shape for larger magnetic Bond numbers, possibly by a bifurcation of a different kind (see figure 1). The droplets axes ratio $b/a$ obtained experimentally in the range of moderate field strength is shown in figure 10. Separate experiments, indicated with thin dashed lines, show a good agreement up to a field $H^{2}\approx 550$ Oe2. Above that, the droplet shape can remain elongated, but relaxes to the energetically more preferable star shape (as shown in section 4.2), if the deformation via the increase of the field is large enough.

The magnetic field is increased by small steps. Nevertheless, as shown in figure 11(a) the shape perturbations and their dynamics (black dots) can be visually detected, which allows to trace their characteristic decay time $\tau$ (dashed line). As in the numerical simulation (see figure 8) the critical slowing down for the decay of perturbations in the case of an oblate shape is shown in figure 11(b). It leads to a critical field is $H_{c}^{2}=125~{}\rm{Oe}^{2}$. Experimentally it is possible to observe the critical slowing down of the prolate shape as shown in figure 11(c), when approaching prolate-oblate transition, where the critical field is $H_{c,2}^{2}=765~{}\rm{Oe}^{2}$.

The compression of the oblate shape is determined from the semi-axes of the visible projection of the droplet approximated by a 3D ellipsoid according to $c/R_{0}=R^{2}_{0}/ab$. By numerical experimentation with formula (2.21) we have found that semi-axes ratio $1/K$ is a linear function of $H_{0}^{2}$ in a broader field range than $c/R_{0}=K^{-2/3}$. Thus from experiment we deduce $1/K=(ab/R_{0}^{2})^{-3/2}$. In figure 12 $1/K$ is plotted as a function of $H^{2}_{0}$. It is adjusted with the relation

[TABLE]

and $H^{2}_{*}=303~{}\rm{Oe}^{2}$. While according to (2.21) $H^{2}_{*}$ reads

[TABLE]

In section 4 it is shown that the ellipsoidal approximation for the calculation of the axis ratio $b/a$ is valid even when the droplet becomes non-ellipsoidal. Thus we use the ellipsoid approximation to fit the data, adjusting $\mu$ with constrained $R_{0}/\gamma$ by relation (50). In figure 10 we plot the ratio $b/a$ obtained analytically from the ellipsoidal approximation for $\mu=10,11,12$. We see that theory gives a very good agreement when $\mu=11$. This value also agrees with the elongation experiment we did in § 5.3.

When we take $\mu=11$ and $R_{0}=10~{}\mu\rm{m}$, we get the surface tension from (50) $\gamma=8.0\cdot 10^{-3}~{}\rm{erg}/\rm{cm}^{2}$ ($R_{0}/\gamma=0.125~{}\rm{Oe}^{-2}$ in reduced numerical unit). This value is slightly larger than the one calculated from the elongation experiment. But as we see in figure 10 that it fits better the axes ratio data. We take this value to evaluate magnetic Bond number in the rotating field experiments.

Another peculiar aspect, for which we can use the results obtained thus far, is for the calculation of the critical magnetic Bond number $Bm_{c}$. According to figure 11 we have $H_{c}^{2}=125~{}\rm{Oe}^{2}$. It results in $Bm_{c}=15.6$ ($R_{0}/\gamma=0.125~{}\rm{Oe}^{-2}$), which is close to the value obtained numerically $Bm_{c}=16.2$. Since energy differences between axial shape and tri-axial are very small near the critical magnetic Bond number then the deviation can either come from numerical errors or some small physical effects in experiments.

As the typical properties of the droplets are now determined we use in the following the magnetic Bond number $Bm$ to characterize the state of the droplet. The rich family of characteristic drop shapes observed during these experimental measurements is displayed in the introduction of this work (see figure 1).

5.5 Comparison of experimental and numerical results

We first compare qualitatively the equilibrium figures obtained experimentally with the numerical simulation results from section 4 in figure 13. Numerous effects can be seen in the experiments and in the simulations. (i) For magnetic Bond numbers under $Bm<20$ droplet becomes flatter with increasing field strength as shown in figure 2 which is very well described with analytical results. Figure  12 shows that this droplet flattening $1/K=c/a=f(H_{0})$ is almost a linear function of field strength square. (ii) At $Bm_{c}=20$, there is a spontaneous symmetry breaking . (iii) For $20<Bm<50$, the worm-like shape elongation increases with the field strength. (iiii)For $50<Bm<80$ a non-ellipsoidal regime where the droplet becomes wider, flatter and forms sharp tips is observed. (iiiii) At $Bm=85$ a re-entrant transition is observed experimentally. We were not able to simulate it numerically since the droplet becomes so flat that it would need more than $10000$ triangles to have a width of at least three elements. It would require a very time consuming simulation.

To explore the star-like shapes we started with an already broken symmetry mesh and put corresponding magnetic field on. We obtained shapes qualitatively similar to the experimental star-like shapes (see figure 5). A more detailed investigation of hysteresis at the re-entrant transition is pending for a future publication. Here we only point out that its observation in numerical simulations is not easy due to the strong compression of the droplet at large magnetic Bond number values.

As an important result we present a quantitative comparison of equilibrium figures obtained experimentally and numerically by drawing calculated equilibrium shapes on top of experimental images as shown in figure 14. Overall we see a good agreement except that simulation tips are a bit shorter either because of numerical errors or because $\mu$ is a bit larger in experiment.

It is also worth noting the differences between experimental data and elliposidal approximation at larger fields, as visible in figure 10 for $H^{2}>550$. For this we see several possible explanations. First, the drop shape at larger fields is no more ellipsoidal, but takes a rather flattened shape with sharper tips. Second, the surface tension of the concentrated magnetic particle system, might become anisotropic under large fields. Third, owing to different field increase steps, as visible in Fig.10, perturbations of different magnitude are imposed on the droplet changing the value of field strength at which the re-entrant transition is observed.

6 Conclusions and discussion

As a general comment we would like to emphasize that the description of the highly magnetic droplets under study here, which are issued from a phase separation, is not obvious. They are in fact based on a new kind of soft magnetic matter with strong interparticle interaction. We show here that many features of their under-field behaviour may be however described by a simple model of magnetic fluid. Moreover it is shown here that the boundary integral equation technique allows one to develop efficient numerical algorithms for the simulation of the equilibrium shapes of these magnetic droplets submitted to an applied field either static or rotating. Numerical and experimental results show a good agreement at equilibrium and as well during shape relaxation processes. Under rotating field, many similar features are observed both in numerical simulations and experiments, namely - instability of the axially symmetric droplet shape with respect to the symmetry destroying perturbations (leading to the oblate-prolate transition), flattening of the prolate shape as the rotating magnetic field strength is increased, formation of the star-fish like equilibrium shapes.

The experimental and numerical results show that the behaviour of the magnetic droplets may be well described by approximating their shape by a tri-axial ellipsoid. However the observed deviations from the simple ellipsoidal approximation, such as the formation of sharp tips, challenge for their theoretical description. Besides that, further development of this work is connected with accounting for the large viscosity contrast between the magnetic droplet and the surrounding fluid. It is possible by applying the boundary integral equation technique to the solution of Stokes equation. Another important issue which should be mentioned with a respect to these strongly magnetic droplets, issued from a phase separation, is the contribution of the magnetic dipolar interactions to their surface tension. An estimate shows that this contribution is just of the same order of magnitude as the surface tension determined in experiments and it may depend on the orientation of the interface with respect to the applied field.

Acknowledgements

Authors are grateful to Delphine Talbot for providing us with the magnetic fluid sample and Rudolfs Livanovics for fruitful discussions. The research was funded by a project of the University of Latvia ”Nano and quantum technologies, innovative materials - Soft materials in electromagnetic field”.

Appendix A Boundary integral equations

A.1 Boundary integral equation for potential

For the potential boundary integral equation derivation we are going to use formula (25) which follows from the Green formula. With subscripts $\psi^{(i)}$ and $\psi^{(e)}$ we are denoting the limiting values of the potential approaching to the surface from both sides of the droplet boundary. For the limiting value from the inside we use relation (25) when $\boldsymbol{y}\in\partial\Omega$

[TABLE]

and for the external region $\Omega_{e}=\infty-\Omega$ when $\boldsymbol{y}\in\partial\Omega$

[TABLE]

where for the first two integrals we changed the direction of the normal vector so it is directed outwards of $\Omega$. The last two integrals of (52) according to (25) correspond to the external field contribution to the potential. Since here we are interested in case of homogenous field whose potential is $\boldsymbol{H}_{0}\boldsymbol{\cdot}\boldsymbol{y}$ then we rewrite (52)

[TABLE]

Multiplying equation (51) by $\mu$ and adding it to (53) and using the boundary conditions $\psi^{(i)}=\psi^{(e)}$ and $\frac{\partial\psi^{(e)}}{\partial n}=\mu\frac{\partial\psi^{(i)}}{\partial n}$ we obtain the boundary integral equation for the potential

[TABLE]

where we omit the suffixes since the potential is continuous.

A.2 Boundary integral formula for normal component

Looking for the solution of the Laplace equation in the form of a single layer potential it is possible to derive the boundary integral equation for the normal component of the magnetic field strength. Numerical calculations showed that the solution of this equation is ill-behaved and so we have chosen an alternative approach deriving an integral relation which allows us to calculate the normal component of the magnetic field strength from the known tangential component.

Let us consider the identity

[TABLE]

where $\boldsymbol{B}$ is the magnetic induction, $\boldsymbol{Q}=\nabla_{\boldsymbol{x}}\times\frac{\boldsymbol{a}}{|\boldsymbol{x}-\boldsymbol{y}|}$ ($\boldsymbol{a}$ is an arbitrary constant vector. It is easy to check that $\nabla_{\boldsymbol{x}}\times\boldsymbol{Q}=(\boldsymbol{a}\boldsymbol{\cdot}\nabla_{\boldsymbol{x}})\nabla_{\boldsymbol{x}}\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}$. Integrating relation (55) throughout the arbitrary region $\Omega$ in dependence on the position of the point with the radius vector $\boldsymbol{y}$ we have

[TABLE]

where $\boldsymbol{n}$ is the external normal to the boundary of the region and the integrals in the case $\boldsymbol{y}\in\partial\Omega$ are taken as the Cauchy principal value integrals.

Applying relation (56) for the region $\Omega_{e}=\infty-\Omega$, where $\Omega$ is the region occupied by the magnetic droplet at $\boldsymbol{y}\in\partial\Omega$ we have (at $\partial\Omega$ we take $\boldsymbol{n}$ as the external normal to the surface of the droplet)

[TABLE]

According to relation (55) last two integrals correspond to contribution of the external field. Since we have a homogeneous field $\boldsymbol{H}_{0}$ (without droplet) everywhere then we rewrite relation (57) as

[TABLE]

where for simplicity we can say that $\boldsymbol{H}_{0}$ corresponds also to magnetic field far away of droplet.

Taking the projection of relation (56) at $\boldsymbol{y}\in\partial\Omega$ and (58) on the normal to the interface and accounting for the boundary conditions $\boldsymbol{B}^{(e)}\boldsymbol{\cdot}\boldsymbol{n}=\boldsymbol{B}^{(i)}\boldsymbol{\cdot}\boldsymbol{n}$ and $\boldsymbol{H}^{(i)}\times\boldsymbol{n}=\boldsymbol{H}^{(e)}\times\boldsymbol{n}$ we obtain the following relation for the calculation of the normal component of the magnetic field

[TABLE]

A.3 Boundary integral equation for normal component of field

While we had a lack of success with direct calculation of the normal field with the boundary integral equation in our studies we have derived a corresponding equation for 3D. Here we briefly give the derivation as it is analogous and uses previous results of derivations (54) and (59).

Again taking the normal projection of relation (56) at $\boldsymbol{y}\in\partial\Omega$ and (58) and using the boundary conditions $\boldsymbol{B}^{(e)}\boldsymbol{\cdot}\boldsymbol{n}=\boldsymbol{B}^{(i)}\boldsymbol{\cdot}\boldsymbol{n}$ and $\boldsymbol{H}^{(i)}\times\boldsymbol{n}=\boldsymbol{H}^{(e)}\times\boldsymbol{n}$ we subtract $\int_{\partial\Omega}(\boldsymbol{B}\times\boldsymbol{n}_{\boldsymbol{x}})\times\nabla_{\boldsymbol{x}}\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}dS_{\boldsymbol{x}}$ and obtain

[TABLE]

where $B_{n}=\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{n}$. This equation is remarkably similar to the boundary integral equation for potential (54) except for normal derivative inside the integral. Using the same numerical procedure as outlined in § 3 (plane triangles, trapezoidal integration for all triangles) we report that this integral equation gives large errors. This is clearly because we don’t integrate the singular quadrature properly which is cumbersome for 3D.

Appendix B Singularity subtraction

Consider identity with $\boldsymbol{y}\ni\Omega$

[TABLE]

which can easily be proven with the divergence theorem and using $\Delta_{\boldsymbol{x}}\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}=0$. For taking the limit $\boldsymbol{y}\to\partial\Omega$ we will split the integration domain into two parts. One integral is for all surface except for the excluded region with a radius $\varepsilon$ centred at $\boldsymbol{y}$ which in the limit $\varepsilon\to 0$ gives us a principal value integral. The other one is around the half sphere whose value can be calculated analytically

[TABLE]

Taking the limit $\varepsilon\to 0$ we conclude that limiting value of the integral (61) is

[TABLE]

which we use for the singularity subtraction in (28).

When $\boldsymbol{y}\ni\Omega$ we have another identity

[TABLE]

since the curl of a potential field is [math]. For taking the limit $\boldsymbol{y}\to\boldsymbol{x}$ we again split the integration domain into two parts. Integral around the half sphere is [math] as the normal vector $\boldsymbol{n}_{\boldsymbol{x}}$ is collinear with $\nabla_{\boldsymbol{x}}\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}$ when $\boldsymbol{y}$ is at the center of half sphere. So we conclude that the principal value integral when $\boldsymbol{y}\in\partial\Omega$ is

[TABLE]

For deriving the formula used in (33) we start with the identity

[TABLE]

Putting in $\boldsymbol{A}\to\nabla_{\boldsymbol{x}}\frac{1}{|\boldsymbol{x}-\boldsymbol{y}|}$, $\boldsymbol{B}\to\boldsymbol{n}_{\boldsymbol{x}}$, $\boldsymbol{C}\to\rm{P}\boldsymbol{H}_{\boldsymbol{x}}$ gives us the following integral identity

[TABLE]

Using the singularity subtraction identities (65) and (63) for both integrals on the right side we obtain

[TABLE]

The right side is now weakly singular and thus can be calculated by means of trapezoidal quadrature for a nonsingular integral and using a transformation of variables for weakly singular elements.

Appendix C Identity

To transform the integral

[TABLE]

we consider the two terms in (69) separately.

[TABLE]

Multiplying (70) by an arbitrary constant vector $\boldsymbol{a}$ and using the relation

[TABLE]

we have

[TABLE]

Last term on the right-hand side of (71) is zero since

[TABLE]

The first term on the right-hand side is transformed as follows

[TABLE]

The first term on the right-hand side of relation (72) is zero since the integration is carried out over a closed surface. Since the vector $\boldsymbol{a}$ is arbitrary we obtain

[TABLE]

The second term in (69) is transformed similarly

[TABLE]

[TABLE]

Using the relation

[TABLE]

we obtain

[TABLE]

The first term on the right-hand side of (76) may be transformed as follows

[TABLE]

The first term on the right-hand side of relation (77) is zero since the integration is carried out over the closed surface. The second term on the right-hand side of relation (77) equals

[TABLE]

As a result for $\boldsymbol{I}$ we have

[TABLE]

[TABLE]

Since the integral of right-hand side equals zero we arrive at the identity used in the text

[TABLE]