The normal field instability under side-wall effects: comparison of experiments and computations

We consider a single spike of ferrofluid, arising in a small cylindrical container, when a vertically oriented magnetic field is applied. The height of the spike as well as the surface topography is measured experimentally by two different technologies and calculated numerically using the finite element method. As a consequence of the finite size of the container, the numerics uncovers an imperfect bifurcation to a single spike solution, which is forward. This is in contrast to the standard transcritical bifurcation to hexagons, common for rotational symmetric systems with broken up-down symmetry. The numerical findings are corroborated in the experiments. The small hysteresis observed is explained in terms of a hysteretic wetting of the side wall.


Introduction
Rotational symmetric systems with broken up-down symmetry become first unstable due to a transcritical bifurcation to hexagons, which is hysteretic (Cross & Hohenberg 1993). Examples are non-Boussinesq Rayleigh-Bènard convection (Busse 1962) and chemical reactions of the Turing type (Turing 1952). The same is true, when a layer of magnetic fluid is exposed to a normal magnetic field. Above a certain critical induction B c , a hexagonal pattern of liquid spikes appears on the surface of the fluid. This striking phenomenon was first reported by Cowley & Rosensweig (1967) and described in terms of a linear stability analysis. This observation 40 years ago triggered numerous efforts to describe also the nonlinear aspects of the phenomenon theoretically. Gailitis (1977), later Friedrichs & Engel (2001) and Friedrichs (2002) and most recently Bohlius et al. (2006) used the principle of free energy minimisation to predict the pattern ordering, wavelength and final amplitude of the peaks on an infinitely extended surface. The numerical computations by Boudouvis (1987) and Boudouvis et al. (1987) specify quantitatively the hysteresis in spike height in unbounded ferrofluid pools. Matthies & Tobiska (2005) calculate also the dynamics of an infinite periodic lattice of peaks.
The experiments, however, are performed with limited amounts of fluid. A finite layer depth in the vertical dimension has been incorporated into the theory by Friedrichs & Engel (2001). According to Lange et al. (2000), the infinitely deep limit is well approximated if the depth exceeds at least the wavelength of the pattern. However, in the horizontal dimension the finite container size has not been considered. Therefore, experimental realizations approximated this limit of an infinitely extended layer by several different approaches. Abou et al. (2001) used a very large aspect ratio, whereas Gollwitzer et al. (2007) employed an inclined container edge. Richter & Barashenkov (2005) used a magnetic ramp to minimize the influence of the border for the Rosensweig instability, whereas Embs et al. (2007) independently applied it to the Faraday instability in ferrofluid.
The question arises, what happens if the container size is intentionally reduced until only a single spike is left. In this case, all symmetries are kept, nonetheless the character of the bifurcation may change. Indeed, in experiments before the seminal work of Bacri & Salin (1984), it was difficult to uncover a hysteresis due to the small container size.
Although there have been numerous experiments in small containers merely because they are simple and cheap, a systematic study of the influence of the constrained geometry on the bifurcation is missing. One reason is, that model descriptions which deal with the finite container size are rare. So far, we know only the work by Friedrichs & Engel (2000) where the free surface is modeled by a four parameter function to fit the measurements by Mahr & Rehberg (1998). In this case, a highly susceptible fluid still showed a hysteretic transition.
In the present article, we demonstrate, that a transcritical bifurcation to hexagons, as found by Gollwitzer et al. (2007), e.g., becomes an imperfect supercritical bifurcation to a single spike if the container size is reduced sufficiently. This is the outcome of a
Material properties of the ferrofluid APG 512a (Lot F083094CX) from Ferrotec Co. † The absolute error of the measurement is unknown. The error given here is taken from the analysis by Harkins & Jordan (1930) numerical model which is able to calculate the stable and unstable solutions for given container size and fluid parameters. It also takes into account the side-wall effects, namely the wetting and the fringing field. We compare the numerical results with our measurements of the surface topography. For the first time, we apply two different techniques which are capable of recording the amplitude (Megalios et al. 2005) and also the full topography of the fluid surface (Richter & Bläsing 2001), to the same experiment.
In the following two sections, we give an overview over the experimental methods. Subsequently, we describe the numerical computations. Finally, we compare all three results.

Measurements of the material properties
We used the ferrofluid APG 512a (Lot F083094CX) from Ferrotec for all experiments. It is based on an ester with a very low vapor pressure, suitable for vacuum pumps. It has an excellent long-term stability. Over one year, the critical induction has not changed by more than 3 %. In contrast to less stable magnetic liquids, the formation of agglomerates in the tips of the Rosensweig spikes was not observed. After applying magnetic fields for an hour, the field was switched off. Neither the visual inspection nor the X-ray images unveiled any agglomerates at the site of the spikes.
The Rosensweig instability is a counterplay between gravitational and surface terms on the one hand and magnetic forces on the other hand (Cowley & Rosensweig 1967). Therefore, a set of basic material properties of the fluid is necessary for a comparison with the theory, namely the surface tension σ, the density ρ and the magnetization curve M(H). These quantities are summarized in table 1.
The surface tension has been measured using a commercial ring tensiometer (LAUDA TE1). This device wets a ring made from platinum wire, pulls it off the fluid surface and determines the maximum force acting on the ring, from which the surface tension can be computed following du Noüy (1919). According to an analyis by  & Jordan (1930), the error for this method is smaller than 0.25 %, given that the density of the fluid and the geometry of the ring are known with sufficient accuracy. The density ρ has been measured using a commercial vibrating-tube densimeter (DMA 4100 by Anton Paar). This device enables us to determine the density with an error of 0.01 %.

Harkins
The contact angle θ c was determined with the contact angle system OCA 20 (Dataphysics) by optical means. Three measurements were performed at the inner side wall of the container, which was tilted by 90 . The difference between advancing and receding angle could not be measured in this way.
The magnetization curve M(H) of the ferrofluid has the biggest influence on the surface pattern. It has been meticulously measured using a fluxmetric magnetometer consisting of a Helmholtz pair of sensing coils with 6800 windings and a commercial integrator (Lakeshore Fluxmeter 480). The sample has been held in a spherical cavity with a diameter of 12.4 mm. This spherical shape ensures a homogeneous magnetic field inside the sample and an exact homogeneous demagnetization factor of 1/3, which makes it possible to get accurate results over the whole range of H. Figure 1 shows the magnetization curve of our ferrofluid. The solid and the dashed line provide two analytic approximations to M(H). The dashed line is a fit with the Langevin approximation for monodisperse colloidal suspensions and provides the constitutive equation (Rosensweig 1985) Only the data points within a range of H ∈ [0 . . . 10 kA/m] have been taken into account for the estimation of the adjustable parameters p = 14.6 kA/m and τ = 0.24 m/kA. A satisfying fit of the whole curve with this equation is not possible, because real ferrofluids consist of magnetic particles with a broad size distribution (Popplewell & Sakhnini 1995). We therefore make use of a model for dense polydisperse magnetic fluids put forward by Ivanov & Kuznetsova (2001). The solid line in figure 1 displays the best fit with that model. The saturation magnetization given in table 1 is extrapolated from there. This extrapolation indicates together with the manufacturer information M S ≈ 26 kA/m±10 %, that this model is very well fitted to our dense ferrofluid (cf. Ivanov et al. 2007).

Measurements of the surface pattern
We fill a cylindrical container, machined from aluminum, with the ferrofluid and expose it to a magnetic induction ranging from B = 7.6 mT to B = 37.7 mT. The depth of the container amounts to 20 mm and the diameter is 29.7 mm. This diameter is chosen by trial such that only one single spike emerges in the centre of the vessel for all magnetic inductions we apply. From the weight of the filled container and the density, we calculate the amount of ferrofluid filled into the container to V 0 = 6.387 ml, which is equivalent to a filling height of D = 9.22 mm. Two complementary experimental methods were used to determine the height of the emerging spike in the centre of the vessel: the X-ray method by Richter & Bläsing (2001) and the laser method by Megalios et al. (2005), which are described in the following.

X-ray method
The X-ray apparatus comprises a stable X-ray point source, that emits radiation vertically from above through the fluid layer. The container is placed midway between a water cooled Helmholtz pair of coils, which generate a DC magnetic field of up to 40 mT. Directly below the container an X-ray camera with 512 × 512 pixels is located, which measures the transmitted intensity at every pixel in one plane underneath the fluid (Richter & Bläsing 2001). This setup is depicted in figure 2. The transmitted intensity of the X-rays is directly related to the height of the fluid above every corresponding pixel. To calibrate this relation, we use a wedge of known size, fill it with ferrofluid and place it in the empty container. In this calibration image, we therefore know the height of the fluid. Figure 3 shows the calibration data from the wedge. These are then fitted with an overlay of three exponential functions  After applying the inverse of (2) to an arbitrary image from the detector, we finally end up with a complete three-dimensional surface topography of the filled container. A reconstruction for one specific field is displayed in figure 4. A survey of the surface topography for different fields is provided by the related movie.   Megalios et al. (2005). The solid lines denote the path of the light reflected on the surface. The path sketched by the dashed lines is used to verify that the laser operates.

Laser method
The laser method developed by Megalios et al. (2005) enables precise, relative measurements of the extrema of the surface topography. Figure 5 depicts the principle of operation. The container with the ferrofluid is situated in a long solenoid, which generates a vertical magnetic field. The solenoid is 33 cm long with an internal diameter of 13 cm and an external diameter of 14 cm. It has 1124 windings and produces up to 21 mT at its centre, with a variation of less than 1 % at the experimental region.
A laser beam is directed at the fluid surface through a semitransparent mirror, which splits the beam into two parts. One part is deflected sideways onto a test point and serves as an indicator, whether the laser operates correctly (the dashed path in figure 5). The other beam is focused on the fluid surface and the reflected light is deflected by the semitransparent mirror onto a photodiode detector (the solid path in figure 5).
The position of the focal spot can be adjusted by means of a micrometre screw. The maximum of the reflected intensity is reached, when the direction of the beam coincides with the normal vector of the surface at the focus spot and the distance of the lens from the surface is equal to the focal length. In normal operation, the beam is oriented vertically -thus the intensity of the signal reaches its maximum when the focus spot hits an extreme point of the surface, namely the top of a spike or the minimum in the centre of the meniscus. By tracing the maximum intensity of the reflected beam and recording the position of the laser optics, we get the position of the spike with micrometre resolution relative to some reference point. Also the absolute height of the spike above the bottom of the container can be determined by setting the reference point at the top of the container edge.

Governing equations and computational analysis
A scheme of a small cylindrical ferrofluid pool in a vertical magnetic field is shown in figure 6. The surrounding air and the embedded ferrofluid are denoted by (a) and (b), respectively. The applied field can be produced either by a pair of Helmholtz coils or a solenoid of suitable dimensions. It is uniform i. e., of constant strength and vertical orientation, in a region far away from the pool. The field uniformity, however, is disturbed in the neighborhood of the pool, due to the demagnetizing field of the pool itself. Therefore, the applied magnetic field could be taken uniform only away from the pool. In the following, the magnetic field distribution and the free surface deformation are taken as axially symmetric about the r = 0 axis (cf. figure 6).
The field distribution in regions (a) and (b) is governed by the equations of magnetostatics. The Gauss law for the magnetization reads where B is the magnetic induction. Since the magnetic field H is irrotational it can be derived from a magnetostatic potential H ≡ ∇u both inside and outside the ferrofluid and, provided that the materials are isotropic, it is parallel to B and so is the magnetization The magnetic permeability µ is constant in non-magnetic media, namely µ a = µ 0 = 4π × 10 −7 H/m; inside the ferrofluid, it depends on the field. Two different constitutive equations are used to account for the field dependence on the magnetization. The first one comes from Langevin's theory for monodisperse colloidal suspensions (equation 1). The second one comes from a polydisperse model by Ivanov & Kuznetsova (2001) that is based on the assumption of a gamma distribution of particle diameters. Writing equation (3) in terms of the magnetostatic potential, u, and taking into account the equation (4) yields inside the non-magnetic phase (a) and inside the magnetic phase (b), respectively. Equlibrium is governed by force balance along the ferrofluid free surface which is stated by the magnetically augmented Young-Laplace equation of capillarity where g is the gravitational acceleration, ρ is the density, σ is the surface tension and ζ is the vertical displacement of the free surface parametrized by the radial coordinate r, i.e., ζ = ζ(r). The upper limit H bs of the integral in the magnetization term is the field strength in the ferrofluid, evaluated at the free surface, i.e. at z = ζ(r).
The reference pressure K is constant at the free surface. The unit normal to the free surface n and the local mean curvature of the free surface 2ℵ are n = −ζ r e r + e z 1 + ζ 2 r , 2ℵ = 1 r d dr where e r and e z are mutually orthogonal unit vectors along the r-and z-axis, respectively, and ζ r ≡ dζ/dr. The reference pressure K in equation (6) is determined by the constraint, that the ferrofluid volume is of fixed amount i.e. we assume an incompressible liquid. The coordinate system, i.e. the location of the z = 0 line, is chosen such that C = 0. The set of the governing equations (4), (6) and (8), needs to be solved for the magnetostatic potential u a (r, z) and u b (r, z), the free surface shape ζ(r) and the reference pressure K, taking into account the following boundary conditions (see also figure 6): u a = u b , µn · ∇u b = µ 0 n · ∇u a at z = ζ(r) and 0 ≤ r ≤ R 0 (10a,b) The subscripts r and z denote differentiation with respect to r and z, respectively. Equations (9) are the conditions that the shape of the free surface and the magnetostatic potential are axially symmetric. Equations (10) are statements of the continuity of the potential and of the normal component of the magnetic induction across interfaces between two media with different magnetic permeabilities. A datum for the potential is fixed by (11). Equations (12) are the conditions that the magnetic field be uniform far away from the pool. A contact angle θ c is prescribed by equation (13) and reflects the wetting properties of the magnetic liquid in contact with the solid wall of the container. The governing equations give rise to a nonlinear, free boundary problem, owing to the presence of the free surface, the location of which enters the equations nonlinearly and is unknown a priori. An additional nonlinearity comes from the constitutive equation for the magnetization of the fluid. Such a problem is only amenable to Here we will only outline the application of the method. Details can be found in previous works by Boudouvis et al. (1988) and Papathanasiou & Boudouvis (1998). The domain is tessellated into nine-node quadrilateral elements between vertical spines and transverse curves whose intersections with each spine are located at distances that are proportional to the displacement of the interface along that spine. The tessellation creates a mesh of nodes and at each node a finite element basis function is assigned that is unity at that node and zero at all other nodes. As the basis functions, we choose a quadratic polynomials of the independent variables r and z. A sample of the computational mesh is shown in Figure 7. The dependent variables u a (r, z), u b (r, z) and ζ(r) are approximated in terms of a truncated set of the finite element basis functions. The governing equations are reduced with Galerkin's method to a set of nonlinear algebraic equations for the values of the unknowns u a , u b and ζ at the nodes and for the value of K.
At fixed values of the physical parameters, the nonlinear algebraic equation set is solved by Newton iteration. The parameter of interest here is the applied magnetic induction B 0 , which appears in the boundary conditions (cf. equation 12a). Solution families, i.e., solutions at sequences of different values of B 0 are systematically traced with first-order continuation. The computational results reported are obtained with a mesh of 24000 nodes. The sensitivity of the spike height to further mesh refinement is practically negligible; namely, less than 0.2 % when doubling the mesh density. Three to five Newton steps are needed for the convergence at each value of the continuation parameter.

Comparison of experimental and numerical results
We record the surface profile for 200 different magnetic inductions B 0 ∈ [7.6 mT . . . 37.7 mT] with the X-ray method described above. Unevitably, the magnetic induction in the neighborhood of the container is distorted and emphasized in comparison to the induction for empty coils. The given values B 0 denote the spatially averaged magnetic induction below the container at a vertical distance of 23.8 mm from the bottom of the fluid layer. This corresponds to the boundary of the computational domain in the calculations. The height and position of the extreme point of the surface topography in the centre has been determined by fitting a paraboloid to a small circular region of the surface with a diameter of 1.5 mm. Figure 8 shows the resulting central heightĥ(B 0 ). The red solid line marks the data for increasing induction. The magnetic fluid first rises at the edge of the container, thus the level of fluid in the centre of the vessel drops. The central height then corresponds to the minimum level in the centre. At a magnetic induction of around 16 mT, a single spike emerges in the centre that continues to rise for increasing induction. The central height then corresponds to the height of this spike. A small hysteresis is found when decreasing the field again. See the blue line in figure 8.
Using the laser method, we performed measurements of the spike height for 29 different magnetic inductions from 0 to 25.65 mT, which are also plotted in figure 8. By focusing the laser beam on top of the container edge, a reference point was taken to get absolute values for the central heightĥ above the bottom of the container, denoted by the open squares. They differ from the X-ray measurement by a shift of 0.7 mm. However, if the reference point is adjusted by this shift, we find a nice coincidence of the data points from both methods, as shown by the open circles in figure 8.
The inaccuracy of the reference point of the laser method can well be explained from the fact, that the laser beam cannot be focused precisely onto the top edge of the container. The vessel is machined from aluminum with a quite rough surface and diffuses the incident laser beam, which leads to the observed shift. This has been experimentally verified by comparing the height measurements of the bare aluminium and a ferrofluid surface at the same level. The difference in the reading is large enough to explain the shift between the laser data and the X-ray data. On the other hand, the accuracy of the reference point of the X-ray measurement depends on the positioning of the calibration wedge. The resolution of this position is limited by the lateral resolution of the detector, which leads to an estimate of the absolute error of 0.2 mm. Due to the roughness of the aluminium vessel, the X-ray data seem to be more precise than the laser data concerning the absolute height in the present experiment. Relatively, both yield practically the same result.
Further deviations may stem from the different ambient temperature at which the data were taken. Whereas in Bayreuth, the lab temperature was stabilized at 21 ± 1 C, the temperature in Athens was 30 C. This leads to a reduced magnetisation and may be the origin of a reduced spike height for higher fields (cf. figure 8).
After successfully comparing the results of the two measurement techniques, we now present the numerical predictions. The results of the computational analysis are depicted in figure 9 together with the X-ray data. The value corresponding to the central heightĥ of the measurements is the height at the axis of symmetry h| r=0 = ζ| r=0 + D, where D denotes the filling level. Two computational equilibrium paths are shown for two different models for the magnetziation M(H). The green dashed line displays the numerical result using Langevin's equation (1), while the black line makes use of the model by Ivanov & Kuznetsova (2001). For magnetic inductions up to 17 mT there is no noticable difference between both results. This is explained by the coincidence of both magnetization laws up to an internal field of H ≈ 10 kA/m, as shown in figure 1. For higher fields, however, Langevin's law is no longer valid. This leads to a rather big deviation between both theoretical curves. Regardless of the underlying magnetisation curve, the numerical solutions show a continuous behaviour ofĥ in the full range of B 0 . In particular, no turning points are traced on the curve. This is mentioned, since a pair of turning points, if existed, should imply a hysteresis in surface deformation observed when increasing and then decreasing B 0 . Thus, from the theoretical calculations we do not expect any hysteresis. We stress this fact, because in the case of an infinitely extended container, a hysteresis is both expected theoretically (cf. Friedrichs & Engel 2001, e.g.) and found in the experiment and numerical calculations ( Gollwitzer et al. 2007, e.g.). Moreover, we tested the influence of the contact angle. A computation for θ c = 20 does not deviate more than 60 µm from the above calculation over the full range of B. Thus, the spike height does only weakly depend on the contact angle.
Next, we compare the measurements with the numerical results. In the full range, the experimental data coincide very well with the more advanced computations taking into account the magnetization law by Ivanov & Kuznetsova (2001). The difference is within only 1 % of the absolute height, except near the threshold, where it amounts to 6 %. This is natural for a sigmoidal function, where close to the steep increase the error can become arbitrarily large, when there are unertainties of the control parameter. The difference between the thresholds in theory and experiment amounts to at most 3 % as can be seen from Figure 9 (b). It shows an enlarged view ofĥ in the immediate vicinity of the threshold. In opposite to the numerical results, we observe in the experiment a small hysteresis between the data for increasing B 0 , as marked by upwards triangles, and the one for decreasing B 0 , denoted by downward triangles. The hysteresis is in the range of ∆B 0 = 0.2 mT. The origin for this hysteresis is a priori not clear.
Note, that Gollwitzer et al. (2007), e.g., measure a hysteresis of ∆B ∞ = 0.17 mT for the same ferrofluid as in our case in a container with a diameter of ≈ 10 × λ c . Remarkably, this value is in the same range as the one observed above. If the hysteresis in our experiment would be of the same nature, it should be much smaller due to the imperfection caused by the container edge (Cross & Hohenberg 1993). Therefore we suspect another mechanism. One candidate is a hysteretical wetting of the cylindrical wall. The difference between the advancing and the receding contact angle can be as large as 10 for a surface that has not been specially prepared (de Gennes 1985, Dussan 1979. Moving the contact line always costs energy. This may explain the small hysteresis of the spike height for increasing and decreasing induction. In our experiment, this effect is important, because the interfacial area between the fluid and the vertical wall is comparable to the free surface of the fluid. The availability of the complete surface topography from the X-ray method allows us not only to compare the central height, but also the shape of the spike or meniscus, respectively. Examples for the free surface shapes at selected values of the magnetic induction are shown in figure 10. The experimental data have been averaged angularly around the centre of the observed spike, which is not exactly in the centre of the container in the experiment. This off-centre distance is rather small (0.05 mm), however it must be taken into account, otherwise the averaging would disturb the shape of the spike.
Similiarly to the comparison of the height alone, there are only slight differences between the computations and the experimentally observed shape. Most notably, we discern a drop at the edge of the container. The reason for this difference is twofold: first, the angular average does not work well near the container border, because the centre of the spike is off-axis, as explained before. Second, the X-ray method has problems to accurately detect the height near the border, where the container wall shadows the X-rays. Further, for magnetic inductions near the threshold (cf. figure 10 (b) and (c)), the height of the tip differs by ≈ 1 mm. The reason is probably a slight shift of the critical induction (cf. figure 9), where the height is very sensitive to small changes of the induction B.
The hysteresis, already observed from the central height alone, manifests itself by a difference of the surface profiles for increasing and decreasing induction. Far away from the threshold, both profiles match nearly perfectly (figure 10 (a) and (d)), while there is a clear difference near the threshold ( figure 10 (b) and (c)). The tip of the spike is considerably smaller for an increasing magnetic induction, while the level of the fluid near the container edge is higher. This corroborates a hysteretical wetting to be responsible for the hysteresis.
Apart from these differences, the deviation between the computed and measured profiles is around 1 %.

Discussion and Conclusion
For a rotational symmetric system with broken up-down symmetry we have reduced the container size until only a single entity remains. In the case of the Rosensweig instability this is a single spike of ferrofluid. Whereas for our fluid, the extended system exhibits a transcritical bifurcation to hexagons, here an imperfect bifurcation sets in, and the axisymmetric free surface deformation evolves supercritically and monotonically. This is at least the outcome of the monolithic finite element approach, which takes into account the side-wall effects, namely the wetting and the fringing field, as well as the polydispersity of the fluid. We find a convincing agreement between theory and two independent measurement techniques, the errors being within 3 % without any adjustable parameter. The nonetheless observed hysteresis is due to the wetting. Our findings immediately raise the issue of what is "in between", regarding the structure of the solution space -that is surface deformation versus applied field and other key parameters -as the size of the pool grows laterally. This will be tackled in a forthcoming publication by Spyropoulos et al. (2009).
A transition from a transcritical backward to an imperfect forward bifurcation under spatial constraints has also been reported by Peter et al. (2005). Their spatial stripe forcing simultaneously breaks the rotational and translational symmetry. In our case it is sufficient to break the translational symmetry.
To conclude, we have quantitatively compared numerics and experiments of the Rosensweig instability in a system of finite size. This is a specific example, how external constraints may change a perfect transcritical bifurcation to an imperfect supercritical one.