2.1 Surface potentials

For a ferrofluid column of height $h(x,y)$ along the $z$ direction and positioned in $(x,y)$ , the surface potentials due to gravity and surface tension, and the volume potential due to the magnetic interactions are respectively written (Rosensweig Reference Rosensweig1985, Reference Rosensweig1987; Bacri & Elias Reference Bacri, Elias, Bourgine and Lesne2011):

with $\unicode[STIX]{x1D70C}$ the fluid density, $g$ the gravitational acceleration, $\unicode[STIX]{x1D70E}$ the fluid surface tension coefficient, $\boldsymbol{H}_{\boldsymbol{a}}$ the applied field, $\boldsymbol{M}$ the ferrofluid magnetization and $\boldsymbol{H}$ the magnetic field modified by the presence of the fluid. In this theoretical study, we neglect the effect of the lying surface, assuming it either superhydrophobic or absent (if for instance the droplet is immersed in another fluid). Note that even if the surface is not superhydrophobic, its contribution to the total energy may be negligible.

Figure 3.

Diagram of the studied system.

For the geometries under consideration (e.g. flat layer or ellipsoidal droplet(s)), the induced field $\boldsymbol{H}$ and the applied field $\boldsymbol{H}_{\boldsymbol{a}}$ can be related through the following self-consistent equation (Rosensweig Reference Rosensweig1985):

where $N$ is the demagnetization factor, which encodes the interactions between the magnetic particles inside the fluid. When dealing with a non-uniform field, relation (2.2) is not well defined. In particular, it is not obvious that a demagnetization factor can be defined. Appendix A shows that in an inhomogeneous field it is not possible to obtain a local relationship of type (2.2), but it is possible to define a ‘demagnetizing series’ $N_{n}(a;\boldsymbol{r})$ , as:

with $a$ any height inside the fluid distribution. In this approach, the demagnetization factor $N$ in relation (2.2) is the first term of the series, $N_{0}(a;\boldsymbol{r})$ . Considering averages, one can show that the second term of the series vanishes, so that relation (2.2) is also true at first order.

Note that a thorough study of the effect of a non-uniform magnetic field on the magnetization is beyond the scope of our paper. We focus on the effect of the magnetic forces (due to the magnetic field derivatives), and not on the effect of the non-uniform magnetization. This would be an interesting and, as far as we know, new study. However, in most cases, the magnetization curve can be approximated by a linear function (using the differential susceptibility), so that terms beyond the first order in (2.3) are not relevant.

We assume the spatial extension of the magnet in the horizontal plane $(x,y)$ to be large compared to that of the ferrofluid. We can thus suppose the magnetic field to be uniform inside the droplets along the horizontal plane, and the magnetization aligned with the external field along the vertical axis. In particular, this amounts to neglecting the horizontal gradient.

For simplicity, we will consider in this section a linear relation for the magnetization:

and we will reduce the non-uniformity of the magnetic field to its gradient:

A development taking into account the full Taylor expansion in a generalized regime for the magnetization is proposed in appendix B. Note that in order to compare with experiments we will need to use the generalized regime for the magnetization. The final expression given by equation (B 18) is also suitable for numerical simulations, if one wants for instance to include higher-order derivatives of the magnetic field.

In the simplified approach, the magnetic volume potential reads:

and the magnetic surface potential is:

Integrating by parts the second term leads to:

where

The first term of (2.9) is the usual magnetic potential considered in several papers on the Rosensweig instability, e.g. Friedrichs & Engel (Reference Friedrichs and Engel2001), Lange et al. (Reference Lange, Richter and Tobiska2007). It cannot be directly integrated because $H$ depends on the interface of the fluid distribution through its boundary conditions. The second term is the gradient contribution: it comes from the homogeneous volume force applied in the whole ferrofluid distribution by the magnet. The third term is a squared gradient contribution, which is opposed to the gradient contribution. The last term is the contribution from the inhomogeneity of the magnetization. Appendix A shows that in most cases this inhomogeneity is not expected to come from the inhomogeneity of the applied field: its main origin should be the geometry of the sample, as for uniformly magnetized samples.

The total surface potential reads:

We can already predict from the above expression three different regimes for the external parameters, leading to different behaviours of the ferrofluid:

First, in the limit case of a weak gradient and a high demagnetization, the physics of the ferrofluid is driven by its demagnetization, which corresponds to the study of Rosensweig, and spawns a normal-field instability. We can then keep the leading-order terms only, i.e. the first 3 terms, and extract a critical magnetic field at which the instability appears, characterized by a modified field-dependent wavelength $\unicode[STIX]{x1D706}$ . We investigate this in § 2.2.

Second, in the high gradient and low demagnetization regime, the gradient plays the main role, essentially through the term $F_{M}$ . In this regime, a single droplet can be observed to fission. We focus on this phenomenon in § 2.5.

Third, in the ‘in-between’ regime, there exist various instabilities leading to first- or second-order phase transitions in the evolution of a droplet’s aspect ratio. This is described in §§ 2.3 and 2.4.

2.2 Rosensweig instability in a non-uniform magnetic field

The first observed behaviour of the ferrofluid distribution when slowly increasing the magnetic field is the so-called Rosensweig instability. From the more general expression of the surface potential given by (B 18), taking the demagnetization to be homogeneous inside the fluid and equal to the one in the usual case of an uniformly magnetized layer, as prescribed by our generalization of the demagnetizing field given in appendix A, and keeping only the first-order terms in $h$ , we can write the surface potential as:

where

and $\tilde{\unicode[STIX]{x1D712}}$ is the differential susceptibility:

Equation (2.12) is mathematically equivalent to the surface potential usually written when describing the Rosensweig instability of a ferrofluid submitted to an uniform field $H_{a}(0)$ within an energetic formalism first used by Gailitis (Reference Gailitis1977). The non-uniformity of the field reduces, under our assumptions, to a modified contribution of the gravity within  $\unicode[STIX]{x1D6FC}$ .

From here, writing the height $h$ as a superposition of plane waves, it can be shown that the amplitude $A$ of a perturbation at the surface of the fluid grows as:

with the frequency $\unicode[STIX]{x1D714}$ related to the wavevector $k$ through the following dispersion relation (e.g. Cowley & Rosensweig Reference Cowley and Rosensweig1967; Rosensweig Reference Rosensweig1985; Salin Reference Salin1993; Abou Reference Abou1998; Andelman & Rosensweig Reference Andelman and Rosensweig2009):

with $M\equiv M(H(H_{a}(0)))$ the magnetization, $r$ the dimensionless permeability ratio, $\unicode[STIX]{x1D707}_{c}$ the chord permeability, $\unicode[STIX]{x1D707}_{t}$ the tangent permeability:

This dispersion relation is valid without assuming  a  linear relation for the magnetization (Cowley & Rosensweig Reference Cowley and Rosensweig1967; Andelman & Rosensweig Reference Andelman and Rosensweig2009). Even though the magnetization is often considered linear, the case of a nonlinear magnetization has also been studied by e.g. Knieling et al. (Reference Knieling, Richter, Rehberg, Matthies and Lange2007), Lange et al. (Reference Lange, Gollwitzer, Maretzki, Rehberg and Richter2016).

If $\unicode[STIX]{x1D714}^{2}>0$ , any disturbance will develop at the interface of the fluid, which signals an instability resulting in the creation of a pattern of spikes. The neutral stability is therefore defined by $\unicode[STIX]{x1D714}=0$ , and the flat layer is unstable if there exist a wavevector fulfilling this condition, i.e. the discriminant of (2.16) is positive:

This defines the critical magnetization (reached at a critical field $H_{a}^{c}$ ), above which the normal-field instability is observed:

For $M>M^{c}$ , a pattern of spikes develops at the surface of the fluid. The distance between two spikes is related to the wavelength of the most unstable mode $\unicode[STIX]{x1D714}_{m}$ , i.e. the one with the maximum growth rate (e.g. Bashtovoi, Krakov & Reks Reference Bashtovoi, Krakov and Reks1985):

Together with the dispersion relation (2.16), equation (2.20) gives the characteristic wavevector $k_{0}$ and the corresponding characteristic wavelength $\unicode[STIX]{x1D706}$ of the pattern. Placing ourselves at the point of appearance of the instability, we get:

where we introduced the capillary wavelength $\unicode[STIX]{x1D706}_{c}$ (which is the characteristic wavelength of the normal-field instability in uniform field (Rosensweig Reference Rosensweig1985)) and the characteristic wavelength from the magnetic field gradient $\unicode[STIX]{x1D706}_{M}$ , which are defined by:

These expressions reflect the competition between gravity and surface tension, which tend to flatten the distribution, and the magnetic interactions, which tend to increase the interface area. From (2.22) we can distinguish two regimes: a low gradient regime, providing us with $\unicode[STIX]{x1D706}\sim \unicode[STIX]{x1D706}_{c}$ , in which the uniform contribution is prevalent and the gradient only slightly increases the value of the critical magnetization, and a high gradient regime, in which $\unicode[STIX]{x1D706}\sim \unicode[STIX]{x1D706}_{M}$ , dominated by the gradient. If the gradient is too high the instability might never appear (in particular if the critical magnetization is above the saturated magnetization). In this second regime, the critical wavelength has been identified in Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013) as an important criterion in the description of the fission of ferrofluid droplets.

The dispersion relation (2.16) has been obtained by mapping the problem of a non-uniformly magnetized ferrofluid layer to the problem of an uniformly magnetized layer, and applying the linear stability analysis known in the latter case. Interestingly, the full linear stability analysis in the case of an inhomogeneous applied field has been carried out by Zelazo & Melcher (Reference Zelazo and Melcher1969), and leads to the same dispersion relation. This theoretical result has later been experimentally corroborated by Catherall et al. (Reference Catherall, Benedict, King and Eaves2003) using liquid oxygen in a magnetic field gradient. They measured the critical magnetic field and the peak separation as function of an effective gravity $\tilde{g}$ defined, in our notations, as:

The relations $M^{c}\propto \tilde{g}^{1/4}$ and $\unicode[STIX]{x1D706}\propto \tilde{g}^{-1/2}$ were in very good agreement with the experiments, which shows that the main impact of a magnetic field gradient on a paramagnetic fluid is to redefine the gravitational acceleration. This is also summarized by Berkovsky et al. (Reference Berkovsky, Medvedev and Krakov1993), who substitute the gravitational force by an effective value:

Interestingly, such mathematical analogy is used in Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013) to describe the fission of a ferrofluid droplet in a magnetic field gradient with the magnetic wavelength $\unicode[STIX]{x1D706}_{M}$ , although they give no further justification to it. We will investigate this in § 2.5.

In the case of a layer of fluid, we see that the gradient introduces two important corrections. First the critical magnetization $M^{c}$ is field dependent and increases with the gradient: the gradient stabilizes the distribution. Stabilizing effects of non-uniform magnetic fields were indeed observed by Zelazo & Melcher (Reference Zelazo and Melcher1969) and are known in the case of a ferrofluid column surrounding a current-carrying wire (e.g. Rannacher & Engel Reference Rannacher and Engel2006). Second, the characteristic wavelength $\unicode[STIX]{x1D706}$ is also field dependent, so that the distance between the spikes at onset of the instability is impacted by the gradient, in particular not expected to be equal to the capillary wavelength as shown by Catherall et al. (Reference Catherall, Benedict, King and Eaves2003).

Moreover, pursuing the linear analysis beyond the threshold in the thick-film regime allows one to find the characteristic wavelength of the evolved pattern, corresponding to the wavevector of the most unstable mode, i.e. the solution of (2.20):

which implies that the peaks should get closer and closer as the magnetic field is increased. This result was already obtained by Abou, de Surgy & Wesfreid (Reference Abou, de Surgy and Wesfreid1997), Lange, Reimann & Richter (Reference Lange, Reimann and Richter2001) in the case of an uniform magnetic field. Equation (2.26) shows that the increasing gradient of a non-uniform magnetic field slows down this effect. Note that the evolved pattern, made of sharped spikes, is nonlinear so that a quantitative agreement between (2.26) and experiments is not expected. In order to make relevant quantitative predictions, one has to carry a nonlinear analysis, as done in e.g. Gailitis (Reference Gailitis1977), Friedrichs & Engel (Reference Friedrichs and Engel2001), Lange et al. (Reference Lange, Richter and Tobiska2007). This is left for future work. We nevertheless investigated this behaviour qualitatively by performing an elementary experiment, considering a thin layer of ferrofluid and slowly increasing the magnetic force by approaching a magnet from below. We observed, beyond a threshold, the appearance of spikes in a hexagonal pattern similar to the one typically observed beyond the normal-field instability threshold, as shown in figure 4. Then, when the magnet is approached beyond this threshold, the peaks are observed to get closer and closer, as predicted by (2.26).

Figure 4.

Rosensweig instability experiment in the case of a non-uniform magnetic field. A large neodymium permanent magnet slowly approaches a thin layer of ferrofluid (with depth of a few tenths of a millimetre) from below. We see in picture (a) the apparition of the spikes characteristic of the normal-field instability, located above the centre of the magnet. The number of spikes increases as the magnet gets closer, as shown in picture (b). In (c) we observe that, because of the growing magnetic field gradient, the spikes are getting closer to one another. The spikes are then separated in individual droplets since the original ferrofluid layer is significantly thin.

Figure 5.

(a) Evolution of the field-dependent wavelength as function of the distance from a typical neodymium magnet of intrinsic magnetization $M_{0}=950~\text{kA}~\text{m}^{-1}$ , height 4 cm, radius 1 cm and various saturated magnetizations. Equation (2.27) is used to obtain the field strength and field gradient as function of the distance. The other parameters are $\unicode[STIX]{x1D70E}=50~\text{mN}~\text{m}^{-1}$ , $\unicode[STIX]{x1D70C}=2000~\text{kg}~\text{m}^{-3}$ , $\unicode[STIX]{x1D712}=1$ . (b) Dispersion relation corresponding to $M_{s}=20~\text{kA}~\text{m}^{-1}$ (green curve of a), at different distances from the magnet.

Importantly, since the critical magnetization increases with the gradient, there is no guarantee that the condition (2.18) is always fulfilled after the instability threshold if one keep increasing the field beyond the critical value. A first-order phase transition may occur at some point, making the whole distribution suddenly collapse to a flat interface. This is shown in figure 5 where is plotted the value of the wavelength as a function of the distance to a typical neodymium magnet (detailed input values are given in the caption), where $1/r$ was put to 1 for simplicity and the magnetic field is computed using the expression of the field created by a finite solenoid (e.g. Camacho & Sosa Reference Camacho and Sosa2013):

As the magnet approaches (from right to left in the plot), we can first see the appearance of the instability, i.e. the formation of spikes, followed by a smooth but quick decrease of the wavelength, i.e. the decrease of the distance between the spikes. Then, as the contribution of the gradient gets higher, the magnetization may return below the threshold so there is no more instability (the frequency $\unicode[STIX]{x1D714}$ is imaginary). The perturbed interface therefore collapses to the initial flat distribution. Note that in this theoretical example, for $M_{s}\lesssim 16.4~\text{kA}~\text{m}^{-1}$ , the normal-field instability never appears.

Such transition was not observable in our qualitative experiment shown in figure 4 because the magnetization of our ferrofluid was too high, and neither was the morphological transition from the hexagonal pattern to a square pattern.

The kind of phase transition evidenced here, due to the competition between the uniform magnetic field contribution and the magnetic field gradient, will have its equivalent in the case of a single ferrofluid droplet studied in § 2.3. The uniform contribution, which elongates the spikes, will equivalently stretch the droplet, while the gradient contribution, which flattens the spikes, will equivalently flatten the droplet until inducing a first-order phase transition from a prolate ellipsoid to an oblate ellipsoid. The fundamental mechanism of such transition is exactly the same as for the transition evidenced in figure 5.

Just as this last remark shows a bridge between the peaks creation and the evolution of the shape of a droplet, a similar connection can be made with the droplet’s fission phenomenon. Indeed, our qualitative experiment displayed the fission of the peaks once the magnet was close enough to the ferrofluid, as if they were droplets lying on the ferrofluid film they originated from. This phenomenon is similar to the one evidenced by Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013) for a single droplet on a hyperhydrophobic surface, and on which we focus in § 2.5. We can also see on picture (c) of figure 4 that, starting with a very thin layer, we end up at high field with separated spikes which form individual droplets. This change of topology, recently studied by Bushueva et al. (Reference Bushueva, Kostarev and Lebedev2011), is due to the fact that the layer is not thick enough for peaks to appear, hence several separated droplets are created instead of several peaks. Bushueva et al. (Reference Bushueva, Kostarev and Lebedev2011) indeed observed that the critical magnetic field imposing such breakup is proportional to the thickness of the layer.

2.3 Evolution of the shape of a single ellipsoidal droplet

We now turn to the main study of this paper and investigate the behaviour of a ferrofluid droplet in a non-uniform magnetic field, again created by approaching a magnet. We do not consider a layer of fluid anymore but an ellipsoidal droplet. The theoretical modelling is in the same spirit as the energetic study done by Bacri & Salin (Reference Bacri and Salin1982).

Figure 6.

Typical distributions adopted by the ferrofluid (various initial volumes and external fields), showing the different phases leading to the fissioning of a droplet into several sub-droplets as the applied magnetic field and magnetic field gradient are increased. (a) Vertical stretching. (b) Splitting phase. (c) Distribution of fluid after fission. (d) Creation of several droplets. The amount of fluid is huge in the second picture in order to capture the transition, but it is not representative of our experiments.

The experimental observation we want to explain is the following: when the magnet is approached, moving along the vertical direction, the droplet lengthens and slims until it suddenly splits into smaller droplets. If we bring the magnet even closer, the droplets split again, and so on, creating henceforth numerous droplets, whose volume shrinks with each division. As observed by Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013), the splitting is irreversible since the droplets’ inertia spawns a distance of separation which can no longer be overcome. Figure 6 illustrates the steps leading to the fission in several droplets. Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013) provide high-quality videos of the phenomenon under interest (see in particular movies S1 and S3).

The physical reason underlying this phenomenon is similar to the one responsible for the appearance of the peaks of the normal-field instability: beyond the threshold, the cost of energy of the increased surface of the droplets is smaller than the loss of magnetic energy due to their increased interface.

In this section we focus on the evolution of the aspect ratio, i.e. before the fission occurs. In the next section we investigate the fission phenomenon.

We start the modelling by computing the previously derived potential for an ellipsoidal geometry. In most of the previous studies, the linear regime $M=\unicode[STIX]{x1D712}H$ was assumed for the magnetization. In our case (and also for the set-up of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013)), we are working with high magnetic fields and the magnetization is in an intermediary regime, closer to saturation. We should therefore use the generalized potential, given in appendix B by (B 18). We keep the first and second magnetic field derivatives, as well as the gradient squared term, but not the product of the gradient and the second derivative nor higher-order terms. As justified in appendix A, we take the demagnetization factor equal to the demagnetization factor of an uniformly magnetized ellipsoid, which is (Bacri & Salin Reference Bacri and Salin1982; Brancher & Zouaoui Reference Brancher and Zouaoui1987; Tyler Reference Tyler2010):

where $e\equiv \sqrt{1-K^{2}}$ is the eccentricity, $K\equiv L/A$ is the inverse aspect ratio, $A$ and $L$ are respectively the major and minor axes.

Under those approximations, the general potential of the ellipsoidal fluid distribution reduces to:

where $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D70C}g+\tilde{F_{M}}$ and:

$V=\unicode[STIX]{x03C0}AL^{2}/6$ is the ellipsoid volume, ${\mathcal{A}}(z)=L^{2}\unicode[STIX]{x03C0}(A-z)z/A^{2}$ is the area of the section of the ellipsoid at height  $z$ , $S=(AL\unicode[STIX]{x03C0}/2)(K+(1/e)\arcsin e)$ is the ellipsoid surface with $e$ its eccentricity.

$A$ and $L$ are related to the inverse aspect ratio thanks to the following:

where $R_{0}$ defined through $V\equiv 4/3\unicode[STIX]{x03C0}R_{0}^{3}$ is the radius of the equivalent sphere.

Finally the potential $\tilde{U} \equiv U/(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E})$ is:

where $\unicode[STIX]{x1D706}=2\unicode[STIX]{x03C0}\sqrt{\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D6FC}}$ , given by (2.22), and

In order to compute the magnetic integral in (2.32), we use the same approximation as in appendix B, generalizing the linear regime for the magnetization as follows:

which is valid either for a small demagnetization factor (which is in general the case for an ellipsoid, $N\lesssim 0.3$ ) or when the second derivative of the magnetization is negligible (it is in particular exact in the linear regime).

We separate the magnetic integral into two linear parts, defining $H_{\ast }$ such that:

where $\unicode[STIX]{x1D712}$ is the initial susceptibility and $\tilde{\unicode[STIX]{x1D712}}\equiv \unicode[STIX]{x2202}M/\unicode[STIX]{x2202}H|_{H_{a}(0)}$ the differential susceptibility.

The magnetization curve is therefore approximated by two connected straight lines:

The integral becomes:

where

Finally the potential is:

where

Note that $B_{\unicode[STIX]{x1D706}}$ is shape dependent because of the presence of $\tilde{\unicode[STIX]{x1D712}}N$ in $\tilde{F_{M}}$ . However, this product is very small since we are usually close to the saturated regime. We therefore neglect this shape dependency in the analytical study.

This potential, function of four dimensionless control parameters $R_{0}^{2}B_{\unicode[STIX]{x1D706}}$ , $R_{0}^{3}B_{\unicode[STIX]{x1D735}}$ , $R_{0}B_{m,1}$ , $R_{0}B_{m,2}$ contains a rich structure depending on the relative strength of the parameters. The evolution of the droplet’s shape is driven by the competition between the gradient, the gravity and the surface tension contributions, which tend to flatten the droplet, and the squared gradient, the second derivative of the magnetic field and the demagnetization, which tend to stretch the droplet. For simplicity, we will now restrict the analysis to the ‘weak gradient regime’ defined by $R_{0}B_{\unicode[STIX]{x1D735}}\ll B_{\unicode[STIX]{x1D706}}$ . We therefore fix $B_{\unicode[STIX]{x1D735}}$ to zero. Qualitatively, $B_{\unicode[STIX]{x1D735}}$ enhances the effect of the demagnetization: $K$ decreases when $B_{\unicode[STIX]{x1D735}}$ is increased, hence the droplet stretches more. Importantly, we notice that the potential (2.40) is not bounded from below in the high gradient limit $R_{0}B_{\unicode[STIX]{x1D735}}\gg B_{\unicode[STIX]{x1D706}}$ . According to (2.40), there exists a threshold beyond which the aspect ratio diverges to infinity, which is of course non-physical. In fact, the limitation to the gradient and second derivative breaks down when the height of the droplet is very high, since the Taylor expansion of the magnetic field is done in terms of $z$ . Higher-order terms in $h^{n},n>3$ should be taken into account in order to counter the squared gradient contribution and suppress the divergence. Note however that gradients created by usual magnets are much too weak to encounter such issue.

Table 1.

Experimental data are estimated from the work of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013). Theoretical values are taken from solutions of (2.42) in which the term $B_{\unicode[STIX]{x1D735}}$ was not neglected (numerical resolution without any assumption, also taking into account the varying demagnetization in $\tilde{F_{M}}$ and the variation of $\tilde{\unicode[STIX]{x1D712}}$ although they have no significant impact). The susceptibility of the ferrofluid is $\unicode[STIX]{x1D712}=3.3$ . The initial volume is $V_{0}=20~\unicode[STIX]{x03BC}$ l and the corresponding spherical radius is $R_{0}=1.7$  mm.

Table 1 provides orders of magnitudes for the four dimensionless numbers, from the data provided by the work of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013). We notice that $B_{\unicode[STIX]{x1D735}}$ can indeed be neglected, as well as $\tilde{\unicode[STIX]{x1D712}}$ if the field is not very small. Note however that both terms will be taken into account as additional corrections in the numerical resolution of § 2.6.

Neglecting $B_{\unicode[STIX]{x1D735}}$ , the minimization of the potential with respect to $K$ gives:

For $\tilde{\unicode[STIX]{x1D712}}\ll \unicode[STIX]{x1D712}$ , this can be simplified to: (This is also true in the linear regime $\tilde{\unicode[STIX]{x1D712}}=\unicode[STIX]{x1D712}$ but then the product $N(K)\tilde{\unicode[STIX]{x1D712}}$ cannot be neglected anymore.)

where

For $B_{\unicode[STIX]{x1D706}}=0$ and $\unicode[STIX]{x1D712}\gg 1$ we retrieve the diagram showing the bistability observed by Bacri & Salin (Reference Bacri and Salin1982), which is shown in figure 7. The aspect ratio increases with the magnetic field strength. For small susceptibilities, this increase is continuous, as in the experimental and numerical study done by Zhu et al. (Reference Zhu, Nguyen, Ramanujan and Huang2011). For strong susceptibilities, there is a first-order phase transition at some point: the aspect ratio jumps from a low value to a higher value. This was observed by Bacri & Salin (Reference Bacri and Salin1982) for a droplet in an immiscible fluid, where the ferrofluid was characterized by an effective susceptibility $\unicode[STIX]{x1D712}=(\unicode[STIX]{x1D707}_{2}-\unicode[STIX]{x1D707}_{1})/\unicode[STIX]{x1D707}_{1}$ where $\unicode[STIX]{x1D707}_{2}$ and $\unicode[STIX]{x1D707}_{1}$ are respectively the ferrofluid permeability and the surrounding fluid permeability. Such very strong susceptibilities are not unrealistic in the case of very concentrated ferrofluids as the one used by Bacri & Salin (Reference Bacri and Salin1984).

Figure 7.

Aspect ratio instability for $B_{\unicode[STIX]{x1D706}}=0$ . Solid lines show stable minima, dotted lines show unstable maxima. $R_{0}B_{m}$ is used as control parameter.

Figure 8.

(a) Modification of the aspect ratio instability for $B_{\unicode[STIX]{x1D706}}>0$ , with the apparition of a threshold and the suppression of the stable branch (solid line) at $K\sim 1$ for high values of $R_{0}^{2}B_{\unicode[STIX]{x1D706}}$ . Here, $\unicode[STIX]{x1D712}=50$ . The yellow curve is the same as the yellow curve of figure 7 ( $\unicode[STIX]{x1D712}=50$ , $R0^{2}B_{\unicode[STIX]{x1D706}}=0$ ). (b) First-order phase transition for $\unicode[STIX]{x1D712}=50$ , $R_{0}^{2}B_{\unicode[STIX]{x1D706}}=1.5$ (corresponding to the green curve). From highest to lowest curve: $R_{0}B_{m}=15,18,20,21,22,24$ . A shift has been applied to the curves in order to have the same value at $K=1$ .

For $R_{0}^{2}B_{\unicode[STIX]{x1D706}}\neq 0$ , equation (2.43) provides modified diagrams as plotted in figure 8. We observe the appearance of a threshold for the control parameter $R_{0}B_{m}$ : for a given $R_{0}^{2}B_{\unicode[STIX]{x1D706}}\neq 0$ , the stable branch does not start at $R_{0}B_{m}=0$ , hence nothing happens ( $K$  is equal to 1) until the magnetic field strength reaches the threshold. Beyond this threshold, there is either the same type of first-order phase transition as before, for small values of $R_{0}^{2}B_{\unicode[STIX]{x1D706}}\neq 0$ , or a new transition. In the first case, the first part of the branch is stable, hence the aspect ratio increases continuously until it reaches the unstable branch and jumps to a much higher value. In the second case, once the threshold is reached, the system enters right away into an unstable state, and thence jumps immediately to a much higher aspect ratio. In both cases, before the critical threshold, the second stable equilibrium is separated from the initial one by a barrier of potential, as shown in figure 8(b). This barrier decreases as $B_{m}$ increases until it vanishes when $B_{m}$ reaches the threshold.

The minimization condition can be rewritten to link the aspect ratio to the parameter  $B_{\unicode[STIX]{x1D706}}$ :

Figure 9(a) shows the corresponding diagram, i.e. the evolution of the inverse aspect ratio $K$ as function of the control parameter $R_{0}^{2}B_{\unicode[STIX]{x1D706}}$ . The droplet flattens as $R_{0}^{2}B_{\unicode[STIX]{x1D706}}$ increases, and may reach an unstable branch: we therefore identify a third possibility of first-order phase transition, now with $R_{0}^{2}B_{\unicode[STIX]{x1D706}}$ as control parameter. Although the droplet starts to flatten continuously, it suddenly collapses to a shape closer to that of a sphere ( $K=1$ ) at some point, and then evolves continuously towards a spherical shape as shown by figure 9. For high enough values of $R_{0}B_{m}$ the droplet directly collapses to $K\geqslant 1$ once the instability is reached. Since the potential (2.40) was derived assuming a prolate ellipsoidal shape, it is not possible for now to investigate the behaviour for $K>1$ . It will be the purpose of the next section, where we will see that the transition in fact occurs from a prolate ( $K<1$ ) to an oblate ( $K>1$ ) ellipsoid, with nothing special about the spherical shape ( $K=1$ ).

This third phase transition corresponds to the one described in the previous section in the case of a thin layer of ferrofluid, and shown by figure 5. It is interesting to note the similarity of those two first-order phase transitions despite the discrepancies in the considered system and the employed formalism.

Figure 9.

(a) Instability for $B_{\unicode[STIX]{x1D706}}>0$ , with $\unicode[STIX]{x1D712}=20$ . Solid line parts are stable, dotted line parts are unstable. (b) First-order phase transition for $\unicode[STIX]{x1D712}=20$ , $R_{0}B_{m}=5$ (corresponding to the blue curve). From lowest to highest minimum: $R_{0}^{2}B_{\unicode[STIX]{x1D706}}=0.32,0.33,0.34,0.36,0.38$ . A shift has been applied to the curves in order to have the same value at $K=1$ .

Figure 10 shows the phase diagram $(R_{0}B_{m},R_{0}^{2}B_{\unicode[STIX]{x1D706}})$ , i.e. the value of the aspect ratio $1/K$ corresponding to the deepest minimum of the potential for each couple of parameters $(R_{0}B_{m},R_{0}^{2}B_{\unicode[STIX]{x1D706}})$ . A single transition line is evidenced. For small susceptibilities, the transition is second order (e.g. blue and red curves of figure 7), i.e. the aspect ratio starts increasing or decreasing continuously as the transition line is crossed. For high susceptibilities, the transition is first order and corresponds to the previous diagrams $(1/K,R_{0}B_{m})$ and $(K,R_{0}^{2}B_{\unicode[STIX]{x1D706}})$ (figures 8 and 9). We see here that the previously evidenced transitions correspond in fact to the same transition, namely the crossing of the line in the phase diagram. It can be either crossed from low aspect ratio to high aspect ratio, increasing $B_{m}$ , or from high aspect ratio to low aspect ratio, increasing $B_{\unicode[STIX]{x1D706}}$ : the gradient contribution enhances the effect of gravity and flattens the droplet, while the demagnetization contribution tends to increase the interface area.

Figure 10.

Phase diagram for $\unicode[STIX]{x1D712}=5$ (a) with second-order phase transition, and for $\unicode[STIX]{x1D712}=50$ , (b) with first-order phase transition. The colours correspond to the value of $\log 1/K$ computed at each couple of parameters $(R_{0}B_{m},R_{0}^{2}B_{\unicode[STIX]{x1D706}})$ by numerically inverting equation (2.43). The deep blue area corresponds more generally to an oblate ellipsoid ( $K\geqslant 1$ ).

Note that, in general, the behaviour of the droplet as function of an approaching magnet will be a combination of the previous conclusions, following a trajectory $(R_{0}B_{m},R_{0}^{2}B_{\unicode[STIX]{x1D706}})$ in the phase diagram, which may cross the phase transition line several times. Numerical examples of such realistic processes will be shown in § 2.6. It would be very interesting to design future experiments to investigate the predicted transition. One could for instance use a combination of two fields: one uniform and the other non-uniform. A very precise control of the increase of the magnetic field at low field would probably be required, as well as a very good hyperhydrophobic surface (or, better, an immiscible fluid). We can predict an experimental observation in 4 steps as the magnet approaches. First, the droplet does not react because the field is below the initial threshold. Second, the droplet stretches to an aspect ratio which can be quite high depending on the field configuration. Third, the droplet shape suddenly collapses to an oblate ellipsoid. Fourth, the aspect ratio slowly decreases until fissions occur. It should be possible to observe more complicated scenarios with several successive transitions, as we will see in § 2.6 through a numerical resolution.

To conclude this section, we show a quantitative comparison of our predicted aspect ratios and the data obtained by Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013). Figure S4 of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013) provides pictures of droplets submitted to a given field configuration. From these pictures we estimated the aspect ratio of the droplets and compared it to our theoretical prediction (taking into account the gradient squared term and the second derivative of the magnetic field in the term $B_{\unicode[STIX]{x1D735}}$ ). For the few cases where the number of created droplets is bigger than 1, we use the generalized formula which will be derived in the next section. The results, shown in table 1, prove a good agreement, keeping in mind that the experimental aspect ratios are estimates. It is interesting to notice the discrepancy for high value of $R_{0}^{2}B_{\unicode[STIX]{x1D706}}$ : the theory predicts an oblate ellipsoid (the phase transition has happened), although the aspect ratio remains bigger than 1 in the observations. This may be due to neglected contributions, like the inhomogeneity of the demagnetization which becomes non-negligible for high gradients. A thorough experimental study of the previously evidenced instabilities remains necessary to properly test our theoretical predictions, although those first results are encouraging.

2.4 Fate of the spherical droplet

We have seen that when the magnetic gradient increases faster than the uniform contribution, the aspect ratio eventually reaches 1. The remaining question is the following: what happens next? One may guess that the gradient will continue to flatten the droplet, transiting from a prolate to an oblate ellipsoid ( $K>1$ ). It is not possible to properly answer this question using the potential (2.40) because it has been derived assuming a prolate ellipsoidal shape. However, it can be easily adapted using the analytic continuation of the previous geometrical formula. In particular, the eccentricity, surface and demagnetization factor are given by (for the demagnetization, see Beleggia, De Graef & Millev Reference Beleggia, De Graef and Millev2006):

Discarding again the $B_{\unicode[STIX]{x1D735}}$ term and assuming $\tilde{\unicode[STIX]{x1D712}}\ll \unicode[STIX]{x1D712}$ , the potential has the same shape as before:

where

Figure 11 shows the continuation of the plots of figures 8 and 9. It reveals that the value $K=1$ has in fact nothing special and the aspect ratio keeps evolving after the transition: the gradient indeed keeps flattening the droplet. Moreover, in general, the previously evidenced first-order phase transition happens between a prolate to an oblate ellipsoid, hence jumps over the spherical shape.

Basically, we therefore identified a transition from a prolate shaped to an oblate shaped droplet, which should be observable. Note however that we have not yet taken into account the possibility of the droplet’s fission, which can occur before the prolate–oblate transition. The study of the fission phenomenon is the purpose of the next section.

2.5 Fission instability

In the previous section we have written the potential of a single ferrofluid droplet (2.40). It can be generalized to an ensemble of $n$ droplets of volume $V/n$ and spherical radius $n^{-1/3}R_{0}$ :

Starting from a single droplet, it is energetically favourable to fission the initial droplet if the following ‘rupture condition’ holds:

where $\unicode[STIX]{x1D716}$ is a phenomenological energy per unit surface taking into account the irreversibility of the phenomenon. It can be put to zero in first approximation.

Figure 11.

(a) Evolution of the potential for $\unicode[STIX]{x1D712}=20$ , $R_{0}B_{m}=5$ and, from bottom curve to top curve, $R_{0}^{2}B_{\unicode[STIX]{x1D706}}=0.32,0.33,0.34,0.36,0.38,0.42,0.5,0.6$ (continuation of the plot of figure 9). This shows a prolate–prolate first-order phase transition followed by a smooth evolution from a prolate to an oblate shape of the droplet. (b) Evolution of the potential for $\unicode[STIX]{x1D712}=50$ , $R_{0}^{2}B_{\unicode[STIX]{x1D706}}=1.5$ and, from top curve to bottom curve, $R_{0}B_{m}=15,18,20,21,22,24$ (continuation of the plot of figure 8). This shows a prolate–oblate first-order phase transition which jumps over the spherical shape. To improve the readability, a shift has been applied to both plots in order to have the same value at $K=1.5$ (a) and $K=2.5$  (b).

Equations (2.51) and (2.52) gather two mechanisms capable of splitting a droplet of magnetic fluid. When the field is uniform, the last term of (2.51) can trigger the instability through the variation of the demagnetization factor, providing that the initial volume is sufficiently high. This topological instability in uniform fields has been evidenced in the experiments by Barkov & Berkovsky (Reference Barkov and Berkovsky1980) and theorized by Berkovsky & Kalikmanov (Reference Berkovsky and Kalikmanov1985) using an energetic approach similar to ours. An important conclusion of these works is that there exists a minimum initial volume allowing the fission, typically of order 0.1 ml for usual magnetic intensities and fluid properties (Berkovsky et al. Reference Berkovsky, Bashtovoi, Mikhalev and Rex1987, Reference Berkovsky, Medvedev and Krakov1993; Blums, Cebers & Maiorov Reference Blums, Cebers and Maiorov1997). This is why the splitting effect was not observed in the experiments of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013), who stated that it is generically not observable with uniform magnetic fields, which is only true for small droplets. In the case of droplets of typical volume 10  $\unicode[STIX]{x03BC}$ l, the variation of the last term of (2.51) is indeed too small to trigger the fission of the initial droplet, and the key mechanism becomes the magnetic force created by the non-uniformity of the field, i.e. the first term of (2.51).

In the following, we therefore assume for simplicity that $K_{1}=K_{2}$ , i.e. the aspect ratio evolves continuously while the fission occurs. We computed that the discontinuity, if it exists, is less than 0.1, hence does not introduce any significant correction. This implies that the last term of (2.51) is negligible and the fission phenomenon is indeed observable only when the applied field is non-uniform, as experimentally concluded in Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013).

Equation (2.52) reduces to the following condition:

where

For $B_{\unicode[STIX]{x1D735}}=0$ , $\unicode[STIX]{x1D716}=0$ , $K=1$ , ${\mathcal{C}}(K)\simeq 8$ and (2.53) is the fission criterion intuited by Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013): the droplet fissions if the characteristic wavelength is smaller than the initial diameter.

Starting with a distribution of $n$ droplets of spherical radius $n^{-1/3}R_{0}$ , the rupture condition becomes:

Theoretically, if this condition holds, we should obtain a distribution of $2n$ droplets (In practice there is always one droplet bigger than the others, hence with a bigger $R_{0}$ , hence which splits first.). Hence the number of droplets for a given $\unicode[STIX]{x1D706}$ is:

where $V_{0}$ is the initial volume of fluid, $K_{c}$ is the inverse aspect ratio which minimizes the potential at the fission.

The theoretical number of droplets predicted by (2.56) does not depend on the uniform contributions $B_{m,1}$ and $B_{m,2}$ : the creation of droplets is essentially driven by the competition between surface tension, which tends to keep a small number of droplets, and the gradient, which tends to create as much as small droplets as possible. The squared gradient tends to elongate the droplets: it introduces a correction which slows down the droplets fission.

The relation $n\propto V_{0}\unicode[STIX]{x1D706}^{-3}$ is in agreement with the experimental results of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013). They indeed identified this proportionality up to a large number of created droplets (see figure S2D of the supplementary materials of Timonen et al. Reference Timonen, Latikka, Leibler, Ras and Ikkala2013). This suggests first that the coefficient ${\mathcal{C}}(K_{c})$ is only weakly dependent of $n$ , i.e. $K_{c}$ is an almost scale-invariant criterion, and second that the contribution from $B_{\unicode[STIX]{x1D735}}$ is small. In first approximation, the slope of the law $n\propto \unicode[STIX]{x1D706}^{-3}$ is indeed constant and moreover universal in the sense that it does not depend on the field configuration or fluid properties. This explains why the slopes are roughly the same in the vicinity of the origin in figure S2D of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013). We can for instance adjust the experimental data provided by figure 2C of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013), with $\unicode[STIX]{x1D706}$ given by (2.22) and a free coefficient as fitting parameter. Figure 12 shows the comparison between these data (blue points) and the proportionality law (red curve). The global factor is 0.38. This factor can be theoretically obtained taking $\unicode[STIX]{x1D716}\simeq 0.5\unicode[STIX]{x1D70E}$ in (2.52), which means that the energetic cost of the irreversibility should be about five times less than the total surface tension energy of the initial droplet, which is big. We will however show in the next section that the squared gradient contribution, taken into account in the green line fit, brings the global factor closer to 1.

Figure 12.

Comparison between our theoretical prediction for the evolution of the number of created droplets and the experimental data (blue points) taken from Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013). The red line shows the agreement obtained with the proportionality law $n\propto \unicode[STIX]{x1D706}^{-3}$ parametrized using the corresponding parameters (magnet M3) given in the supplementary materials of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013), and with the slope as fitting parameter. The green dashed line shows the theoretical prediction obtained by taking into account $B_{\unicode[STIX]{x1D735}}$ in (2.56) and a supplementary global factor of 0.75.

We observe in Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013), that the slope of (2.56) is in fact not universal: it weakly depends on the magnet. This is due to the squared gradient contribution which is not negligible anymore and slows down the fission process. One can check that the hierarchy of the slopes in figure S2D of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013) is indeed in agreement with the hierarchy of the gradients. For a larger number of droplets, the proportionality is not even verified anymore: nonlinear deviations are observed. The numerical resolution of the next section will show that those are not due to the squared gradient contribution. We suspect that they are due to the radial variation of the field in the horizontal plane since for a large number of droplets the droplets are less confined, hence not aligned with the centre of the magnet anymore.

Figure 13.

Numerical resolution for the evolution of the aspect ratio, obtained by solving (2.42), taking into account both $B_{m,1}$ and $B_{m,2}$ and also taking into account the squared gradient and the second magnetic field derivative in $B_{\unicode[STIX]{x1D735}}$ . (a) $\unicode[STIX]{x1D712}=10$ , $M_{0}=320~\text{kA}~\text{m}^{-1}$ . (b) $M_{s}=200~\text{kA}~\text{m}^{-1}$ , $M_{0}=320~\text{kA}~\text{m}^{-1}$ . (c) $\unicode[STIX]{x1D712}=40$ , $M_{s}=200~\text{kA}~\text{m}^{-1}$ . Between (a) and (b) the red curve is the same. Between (b) and (c) the green curve is the same.

2.6 Numerical resolution

In a first part, we solve (2.40) numerically in order to obtain the evolution of the aspect ratio $1/K$ as function of the field configuration, which evolves when decreasing the distance between the fluid and the magnet. The squared gradient and second magnetic field derivative are taken into account in $B_{\unicode[STIX]{x1D735}}$ . $B_{m,1}$ and $B_{m,2}$ are both taken into account without further assumption. In figure 13(a), we use as input a ferrofluid of susceptibility $\unicode[STIX]{x1D712}=10$ , surface tension $\unicode[STIX]{x1D70E}=50~\text{mN}~\text{m}^{-1}$ , density $\unicode[STIX]{x1D70C}=2000~\text{kg}~\text{m}^{-3}$ and various saturated magnetizations. The magnetic field is computed using equation (2.27) with a radius $R=2$  cm, a height $h=4$ cm and an intrinsic magnetization $M_{0}=320~\text{kA}~\text{m}^{-1}$ (e.g. a ferrite magnet). The droplet has a volume of 10  $\unicode[STIX]{x03BC}$ l. As the distance $d$ to the magnet is decreased, for $M_{s}\lesssim 140~\text{kA}~\text{m}^{-1}$ , the aspect ratio is observed to evolve continuously, first increasing as the uniform contribution increases and then decreasing once the gradient contribution becomes predominant. But for higher values of the saturated magnetization we see the apparition of a succession of various phase transitions: while the magnet is approached (from right to left on the plot), $1/K$ first becomes bigger than one (prolate ellipsoid) in a second-order phase transition and then continuously increases, until $R_{0}^{2}B_{\unicode[STIX]{x1D706}}$ becomes high enough to have a first-order phase transition evidenced in figure 9, i.e. $1/K$ collapses to a value smaller than one (oblate ellipsoid). This transition is caused by the magnetization becoming saturated, hence drastically reducing the increase of the uniform contribution. Then when the contributions of $R_{0}B_{m,1}$ and $R_{0}B_{m,2}$ become high enough again, a second first-order phase transition occurs, of the kind evidenced in figure 8, i.e. $1/K$ suddenly jumps to a higher value, and then decreases again because of the gradient domination, until becoming again smaller than one.

We did the same kind of simulation but fixing the saturated magnetization at $M_{s}=200~\text{kA}~\text{m}^{-1}$ , and increasing the susceptibility, while taking the other parameters unchanged. The results are shown on figure 13(b). We still observe the phase transitions described above, but the first one becomes first order: for high enough values of $\unicode[STIX]{x1D712}$ , the aspect ratio almost immediately jumps to a higher value (see figure 8).

Finally, after fixing $\unicode[STIX]{x1D712}=40$ , we increase the intrinsic magnetization of the magnet in figure 13(c). We can see that the two middle phase transitions disappear while the first first-order phase transition remains (the other small discontinuities are numerical errors).

These numerical results show that taking into account all contributions does not prevent the occurrence of the previously evidenced phase transitions with realistic ferrofluid and magnet parameters, and should therefore be observable. Furthermore, it shows that in realistic situations, all the instabilities and phase transitions can appear as coupled or superposed. Still, it is necessary to check their existence through experimental studies.

Table 2.

Characteristics of the magnets. $R$ : radius, $h$ : width, $M_{0}$ : intrinsic magnetization.

In a second part, we focus on a numerical study of the fission behaviour. For a given field configuration, the minimization of the potential of (2.51) with respect to $K$ gives the droplets aspect ratio for a fixed number of droplets and (2.56) gives the number of droplets for a fixed aspect ratio. We simulated the behaviour of a ferrofluid droplet submitted to an increased field by numerically solving this system. As inputs, we used the characteristics of our magnets and magnetization curve (see table 2, figures 15, 16), as well as the characteristics of the magnets and fluid used by Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013).

First we need to calibrate (2.56). In the case of the magnets of Timonen et al., this is done by correcting the numerical solution obtained in the case of their magnet M3, comparing it to the experimental data of Timonen et al. and the simple proportionality law plotted in figure 12. The numerical resolution corresponds to the green dotted line. The fit is obtained by introducing a factor equal to 0.75 multiplying ${\mathcal{C}}(K_{c})$ in (2.56). This is equivalent to taking $\unicode[STIX]{x1D716}\simeq 0.1\unicode[STIX]{x1D70E}$ , hence an irreversibility about 20 times smaller than the surface tension energy, which is realistic and shows that the loss of energy due to the irreversibility may satisfactorily explain the small discrepancy between the experimental and predicted factors. In the case of our magnets, the factor is determined by comparing the numerical resolution to the experiments which will be performed in § 4.2. We need to multiply ${\mathcal{C}}(K_{c})$ by 1.35 (the same factor for our four magnets). We attribute the difference with the factor 0.75 obtained for the magnets of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013) to the difference of geometry. Since we are working on a weakly hydrophobic surface, our droplets do not have a proper ellipsoidal shape, hence global factors are likely to be introduced in all geometrical relations used to derive the potential (2.29), e.g. equations (2.28), (2.31) and (2.33).

Figure 14.

Numerical resolution of the evolution of the aspect ratio (solid lines) and the number of droplets (dotted lines) as functions of $\unicode[STIX]{x1D706}^{-3}$ , for an initial droplet of volume $10~\unicode[STIX]{x03BC}$ l. This simulates an approach from $d\sim 8$  cm to $d\sim 0.5$  cm. (a) The input parameters correspond to the experiments from the work of Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013), in the cases of their four magnets described in the supplementary materials. The correction factor is 0.75. (b) The input parameters correspond to our non-diluted ferrofluid and our four magnets described in § 3.2. The correction factor is 1.35.

Figure 15.

Crosses: experimental magnetization curve. Dotted lines: fit obtained with (3.2), and $M_{s}$ and $\unicode[STIX]{x1D712}$ as fitting parameters.

Figure 16.

(a) Magnetic field strength of our four magnets as a function of the distance from their surface on the vertical axis. (b) Ratio of the second derivative of the magnetic field over the gradient, for a typical droplet height $A=5$  mm.

Taking into account both factors in the resolution, we obtain the results in figure 14(a) for the magnets used by Timonen et al. (Reference Timonen, Latikka, Leibler, Ras and Ikkala2013) and on the right for our magnets, whose characteristics are detailed in § 3.2. The gradient squared is taken into account in $B_{\unicode[STIX]{x1D735}}$ , and the variation of the demagnetization factor with the droplet shape is taken into account in $B_{\unicode[STIX]{x1D706}}$ . The neglected quantities are therefore the horizontal gradient, the higher-order derivatives and the inhomogeneity of the demagnetization.

The eight behaviours from the plots are similar, and allow us to observe four steps. First: the aspect ratio increases. Second: after reaching its maximum, it starts decreasing. Third: the fission occurs. Fourth: other fissions occur while the aspect ratio converges towards a constant value around one. This process perfectly describes the experimental observation (see movie S1 of Timonen et al. Reference Timonen, Latikka, Leibler, Ras and Ikkala2013). It can be understood in a qualitative way: first the gradient is negligible, hence the shape is driven by the demagnetization which increases the aspect ratio, exactly as observed by Bacri & Salin (Reference Bacri and Salin1982), Afkhami et al. (Reference Afkhami, Tyler, Renardy, Renardy, Pierre, Woodward and Riffle2010), Zhu et al. (Reference Zhu, Nguyen, Ramanujan and Huang2011). The relevant diagram is the one plotted in figure 7, with a curve of the blue one type ( $\unicode[STIX]{x1D712}\sim 5.4$ ). Then the gradient becomes non-negligible and its contribution opposes the demagnetization: the aspect ratio goes backwards and the droplet flattens. The relevant diagram showing this decrease is the one plotted in figure 9. Then the gradient becomes so high that it is energetically favourable to create several droplets, resulting in an increased surface but smaller height. The aspect ratio then quickly reaches 1 after the first fissions occurred. Eventually, the aspect ratio is fixed to 1 and the number of droplets converges to the law $n\propto \unicode[STIX]{x1D706}^{-3}$ , with, as explained above, the same slope in our case (the squared gradient correction is negligible) but different slopes in the Timonen et al. case (the squared gradient correction is not negligible and induces a hierarchy). Note that the ‘Flat neodymium’ magnet is not strong enough to break the droplet after the first fission.

3.1 Characteristics of the ferrofluid

The experiments are carried out with an aqueous ionic ferrofluid, with cobalt ferrite particles ( $\text{CoFe}_{2}\text{O}_{4}$ ) with a mean size of about 10 nm (see Neveu-Prin et al. Reference Neveu-Prin, Tourinho, Bacri and Perzynski1993), of density $\unicode[STIX]{x1D70C}=2150\pm 10~\text{kg}~\text{m}^{-3}$ and surface tension $\unicode[STIX]{x1D70E}=50\pm 2~\text{mN}~\text{m}^{-1}$ , provided by Sophie Neveu (Université Pierre et Marie Curie, CNRS) and Florence Elias (Laboratoire Matière et Systèmes Complexes, Université Paris Diderot, CNRS). We repeated the experiments with its 2 times diluted version, of density $\unicode[STIX]{x1D70C}=1570\pm 10~\text{kg}~\text{m}^{-3}$ and surface tension $\unicode[STIX]{x1D70E}=65\pm 2~\text{mN}~\text{m}^{-1}$ .

The magnetization curve has been measured by the ferrofluid manufacturer using a vibrating sample magnetometer. It is given in figure 15. The magnetic susceptibility and saturated magnetization for the non-diluted ferrofluid are $\unicode[STIX]{x1D712}=5.4$ and $M_{s}=55~\text{kA}~\text{m}^{-1}$ .

The Langevin equation (e.g. Behrens et al. Reference Behrens, Bönnemann, Modrow, Kempter, Riehemann, Wiedenmann, Odenbach, Will, Thrams, Hergt and Odenbach2009):

with $M_{s}$ the saturated magnetization, $H$ the magnetic field, $m$ the magnetic moment of a particle, $k_{B}T$ the thermal energy of a particle and $\unicode[STIX]{x1D712}$ the initial susceptibility, describes only approximatively the magnetization curves of our original ferrofluid and its diluted version. This well-known discrepancy is mostly attributed to the polydispersity of the fluid, and also to the interactions between the particles in concentrated fluids (Shliomis Reference Shliomis1974; Berkovsky et al. Reference Berkovsky, Medvedev and Krakov1993), even if the latter can be taken into account in more sophisticated models (Ivanov et al. Reference Ivanov, Kantorovich, N Reznikov, Holm, Pshenichnikov, V Lebedev, Chremos and J Camp2007). We chose consequently to fit our data, as plotted in figure 15, using the Vislovich approximation (Vislovich Reference Vislovich1990):

where $H_{T}$ is the magnetic field strength such that $M(H_{T})=M_{s}/2$ . The magnetization curve of the diluted ferrofluid is obtained by multiplying the one of the non-diluted ferrofluid by the ratio of the densities (i.e. a factor 0.73). Relation (3.2) is used to determine the magnetization $M(H_{a}(0))$ and the differential susceptibility $\tilde{\unicode[STIX]{x1D712}}$ in the theoretical expressions.

3.2 Magnets

Our experiments were carried out using four different types of axisymmetric permanent magnets: two piles of neodymium magnets, labelled ‘Neodymium 1’ and ‘Neodymium 2’, one pile of ferrite magnets labelled ‘Ferrite’ and one large flat neodymium magnet labelled ‘Flat neodymium’. Their respective dimensions are listed in table 2. Their calibrations are shown in figure 16(a). The values of the magnetic field are measured using a Hirst model GM07 gaussmeter, for the field on the $z$ axis as a function of the distance $d$ to the magnet. The fits are obtained using the analytic expression of the magnetic field created by a cylindrical magnet (e.g. Camacho & Sosa Reference Camacho and Sosa2013), which is analogous to the magnetic field created by a finite solenoid:

with $R$ the magnet radius, $h$ the magnet height and $M_{0}$ the intrinsic magnetization to be adjusted. The order of magnitudes for the latter are tabulated: $M_{0}\sim 950~\text{kA}~\text{m}^{-1}$ for the neodymium magnets and $M_{0}\sim 200~\text{kA}~\text{m}^{-1}$ for the ferrite magnet. The values used to obtain the fits are indicated in table 2.

We can infer from this plot a clear difference between the behaviours of the ‘Flat neodymium’ magnet and the other magnets. In the case of the ‘Flat neodymium’, the gradient is not a monotonous function of $d$ : its maximum value is reached at $d=1.15~\text{cm}\neq 0$ , and it is almost vanishing at $d=0$ . This difference in the generated magnetic field predicts, starting from the same theoretical expressions, the observation of clear discrepancies in the behaviour for small values of $d$ between the different magnets. These discrepancies should play the role of indicators to check the robustness of our theory in the experiments.

For our magnets, we typically have:

where $A\sim 5$  mm is the typical height of the droplets. Figure 16(b) shows the ratio of the second derivative and the gradient for a typical height of 5 mm. The second derivative is in general small but non-negligible: it will introduce an observable correction. It becomes negative for the ‘Flat neodymium’ magnet for $d\lesssim 1.2$  cm (this is why, in absolute value, the derivative is discontinuous at $d\simeq 1.2$  cm).

3.3 Overall procedure

We determined the droplet’s volume thanks to a Denver Instruments Pinnacle scale, model PI-403, sensitive to 1 mg. We estimated the uncertainty resulting of manipulation errors to be approximately 3 mg. We then placed the droplet on a horizontal support, vertically aligned with one of our four magnets placed underneath. The support was maintained using a metal arm and its horizontality was checked using a spirit level. We slowly raised the magnet using a horizontal lifting plate, being careful to stay as much as possible in the quasistatic regime, until we observed the fission of the droplet, i.e. two separated droplets like in figure 6(c). Note that in figure 6(b), two non-separated peaks are observed but this is only due to the large amount of fluid required to get a picture of the transition. In practice, the splitting occurs quasi instantaneously: when reaching the splitting instability, an infinitesimal displacement of the approaching magnet is sufficient to spawn a clearly visible gap between the two daughter droplets (as in figure 6 c). Note that in our experiment this ‘infinitesimal’ displacement is always smaller than the uncertainty on the measure of the distance to the magnet (1 mm), so that there is no ambiguity on the fission criterion.

We can then obtain the distance between the magnet and the fluid at first fission for a given droplet mass, with a precision of 1 mm. We carried out this experiment for each magnet with both ferrofluids (diluted and non-diluted), with the exception of the ‘Flat neodymium’ magnet in the case of the diluted ferrofluid. In this case, the magnet needed to be approached very close, and the vertical magnetic force in the droplet was too weak compared to the horizontal magnetic field gradient. The droplet was observed to move toward the edge of the magnet and split earlier than expected, leading to a huge discrepancy in comparison with all the other cases, rendering its study irrelevant.

Our surface, made out of a thin polystyrene plastic cup, is not superhydrophobic. Yet, it is hydrophobic enough for the splitting of the droplets to be clearly observed, and the created droplets are completely separated, as seen in figure 6 (except figure 6(b), which is not representative of our experiments). The additional surface tension term in the potential is negligible and we will therefore neglect this issue while carrying out the experiments.

In the following, we focus the study on the first fission, measuring the distance between the magnet and the fluid at first splitting for different initial droplets’ volumes.

Figure 17.

Initial droplet’s volume as a function of the distance from the magnet at first fission, compared between the magnets and between the two concentrations of ferrofluid. Each colour corresponds to a specific magnet: blue: neodymium 1, red: neodymium 2, green: ferrite, black: flat neodymium. Circles: original non-diluted ferrofluid, triangles: diluted ferrofluid. The theoretical curves are obtained by numerically inverting equation (2.56) and introducing a global factor of 1.35, while taking into account the squared gradient and second magnetic field derivative into the term $B_{\unicode[STIX]{x1D735}}$ .