2. Experimental set-up

Our experimental set-up is sketched in figure 2(a). An X-ray point source emits radiation vertically from above through the container with the ferrofluid, which is placed midway between a Helmholtz pair of coils. Underneath the container, an X-ray camera records the radiation passing through the layer of ferrofluid. The intensity at each pixel of the detector is directly related to the height of the fluid above that pixel. Therefore, the full surface topography can be reconstructed after calibration (Richter & Bläsing Reference Richter and Bläsing2001; Gollwitzer et al. Reference Gollwitzer, Matthies, Richter, Rehberg and Tobiska2007).

Figure 2. Set-up of the apparatus for measuring the Rosensweig instability. Scheme of the assembled set-up (a), and exploded view of the container with the coils (b).

The container, which holds the ferrofluid sample, is depicted in figure 2(b). It is a regular octagon machined from aluminium with a side length of 77 mm and two concentric inner bores with a diameter of 140 mm. These circular holes are carved from above and below, leaving only a thin base in the middle of the vessel with a thickness of 2 mm. On top of the octagon, a circular lid of aluminium is placed, which closes the hole from above (see figure 2 b). Each side of the octagon is equipped with a thermoelectric element QC-127-1.4-8.5MS from Quick-Ohm. These are powered by a 1.2 kW Kepco KLP-20-120 power supply. The hot side of the thermoelectric elements is connected to water-cooled heat exchangers. The temperature is measured at the bottom of the aluminium container with a Pt100 resistor. The temperature difference between the centre and the edge of the bottom plate does not exceed 0.1 K at temperature $10.0\,^{\circ }\text{C}$ measured at the edge of the vessel. A closed loop control, realised using a computer and programmable interface devices, holds the temperature ${\it\theta}$ of the vessel constant with an accuracy of 10 mK.

The container is surrounded by a Helmholtz pair of coils, thermally isolated from the vessel with a ring made from a flame-resistant composite material (FR-2). As shown in figure 2, the diameter of the coils is intentionally smaller than the diameter of the vessel in order to introduce a magnetic ramp. With this arrangement, the magnetic field strength falls off towards the border of the vessel, where it reaches 80 % of its value in the centre. The ramp is introduced in order to minimise the differences between the amplitudes obtained in a finite and an infinite geometry. For details of the magnetic ramp, see Knieling, Rehberg & Richter (Reference Knieling, Rehberg and Richter2010). Filling the container to a height of 5 mm with ferrofluid enhances the magnetic induction in comparison with the empty coils for the same current  $I$ . Therefore $B(I)$ is measured immediately beneath the bottom of the container, at the central position, and serves as the control parameter in the following.

3. Properties of the magnetic fluid

For the experiments, we use the ferrofluid APG E32 from Ferrotec Co. Table 1 lists the important material parameters. The density ${\it\rho}$ was measured utilising an electronic density meter (DMA 4100) from Anton Paar. The surface tension was measured by means of a ring tensiometer (Lauda TE 1) and corroborated by the pendant drop mehode (Dataphysics OCA 20).

The viscosity ${\it\eta}$ has been measured in the temperature range of $-5\,^{\circ }\text{C}\leqslant {\it\theta}\leqslant 20\,^{\circ }\text{C}$ , using a commercial rheometer (MCR 301, Anton Paar) with a shear cell featuring a cone–plate geometry. At room temperature, the magnetic fluid with a viscosity of ${\it\eta}=2~\text{Pa}~\text{s}$ is 2000 times more viscous than water. The value of ${\it\eta}$ can be increased by a factor of 9 when the liquid is cooled to $-5\,^{\circ }\text{C}$ . The temperature-dependent viscosity data can be well described by the Vogel–Fulcher law (Rault Reference Rault2000)

with ${\it\eta}_{0}=0.48~\text{mPa}~\text{s}$ , ${\it\psi}=1074~\text{K}$ and ${\it\theta}_{0}=-107.5\,^{\circ }\text{C}$ , as reported in detail by Gollwitzer, Rehberg & Richter (Reference Gollwitzer, Rehberg and Richter2010). For the present measurements, we chose a temperature of ${\it\theta}=10\,^{\circ }\text{C}$ , where the viscosity amounts to ${\it\eta}=4.48~\text{Pa}~\text{s}$ according to (3.1).

The magnetisation, $\boldsymbol{M}(\boldsymbol{H}_{F})$ , where $\boldsymbol{H}_{F}$ is the magnetic field in the fluid, has been measured by means of a fluxmetric magnetometer (Lakeshore Model 480). A spherical cavity with a volume of 1 ml was constructed from Perspex and filled up to the brim with the magnetic fluid. By means of X-rays we checked for a bubble-free filling. Figure 3(a) presents the measured data together with a fit (solid line) by the equation

utilised by Fröhlich (Reference Fröhlich1881) and Vislovich (Reference Vislovich1990), where $M=|\boldsymbol{M}|$ , $M_{\mathit{S}}$ is the saturation magnetisation, $H_{\mathit{F}}$ is the magnitude of the applied magnetic field in the fluid and ${\it\mu}_{\mathit{r},i}$ is the initial magnetic permeability.

Fitting the data with (3.2) we obtain for the initial magnetic permeability ${\it\mu}_{\mathit{r},i}=4.57$ , and for the saturation magnetisation $M_{\mathit{S}}=23.15~\text{kA}~\text{m}^{-1}$ . Equation (3.2) is used to calculate the chord permeability

and the tangent permeability

which are displayed in figure 3(b) by a dashed red (blue) line, respectively. The effective permeability

is the geometric mean of the two permeabilities (Cowley & Rosensweig Reference Cowley and Rosensweig1967) and is marked in figure 3(b) by the black solid line.

Figure 3. (a) Measured magnetisation, $M$ , curve versus the field in the fluid, $H_{F}$ , and (b) related chord (red dashed), tangent (blue dashed) and effective (black solid line) permeabilities. The dotted line marks ${\it\mu}_{\mathit{eff}}$ at $H_{c}$ .

By means of ${\it\mu}_{\mathit{eff}}(M)$ , the material parameters ${\it\sigma}$ and ${\it\rho}$ , and the gravitational acceleration $g=9.81~\text{m}~\text{s}^{-2}$ , the critical magnetisation for a semi-infinite layer of ferrofluid is approximated by the implicit equation (Rosensweig Reference Rosensweig1985)

By numerically solving the implicit equation, we obtain $M_{\text{c}}=6.2303~\text{kA}~\text{m}^{-1}$ , $H_{\text{c}}=2.3872~\text{kA}~\text{m}^{-1}$ and $B_{\text{c}}={\it\mu}_{0}(H_{\text{c}}+M_{\text{c}})=10.83~\text{mT}$ . Note that this value is just an approximation for a semi-infinite layer of ferrofluid. Revisiting figure 3(b), we see that ${\it\mu}_{\mathit{eff}}$ drops considerably for increasing field. At the critical field, we find ${\it\mu}_{\mathit{eff}}(H_{c})=3.2$ (cf. the dotted line). The critical induction found by fitting amplitude equations to the experimental data (Gollwitzer et al. Reference Gollwitzer, Rehberg and Richter2010) is about 4 % higher (cf. table 1).

Figure 4. (a) The bifurcation diagram (black line) of this pattern is obtained by fitting equation (11) in Gollwitzer et al. (Reference Gollwitzer, Rehberg and Richter2010) to about 170 000 measured data. Via the path indicated by the arrows, one obtains localised patches, as shown in (c). The magnetic induction $B$ is kept fixed during (1) for a time delay of ${\it\tau}_{1}$ , then $B$ is switched to (2) for a time delay of ${\it\tau}_{2}$ before finally being switched to (3) until equilibrium is reached at point (3b). (b) A reconstruction of the height field of the domain-covering hexagons at point (1b). (c) Gives a reconstruction of the height field of a three-spike patch at point (3b) in the bifurcation diagram.

We stress that the mapping between a linear and a nonlinear $M(H)$ by (3.5) is only valid in the linear regime of pattern formation, i.e. for the onset of the instability, but not for the fully developed nonlinear patterns.

4. Experimental results

In the investigated regime, the plain surface and domain-covering regular up-hexagonal (i.e. hexagons whose maximum amplitude is positive) patterns are bistable ( $10.95~\text{mT}\leqslant B\leqslant 11.23~\text{mT}$ ), as shown by the bifurcation diagram in figure 4(a), which stems from a fit to about 170 000 measured data points (Gollwitzer et al. Reference Gollwitzer, Rehberg and Richter2010). Figure 4(b) displays an experimental realisation of the domain-covering hexagonal pattern observed at the upper branch, the infinite extension of which has the symmetry of the dihedral group $\mathbb{D}_{6}$ .

Owing to the high viscosity of the magnetic fluid, the formation of Rosensweig patterns takes 60 s. This allows us to vary the magnetic induction $B$ (the control parameter) and the pattern amplitude $A$ almost independently. In this way we can enter the region around the unstable branch of the instability, where homoclinic snaking is suspected, from the side, as sketched in figure 4.

Starting from a flat layer, $B$ is always switched to an overcritical value $B_{1}=11.45~\text{mT}$ , cf. position (1a) in figure 4. With time, the amplitude of the domain-covering hexagon pattern increases along path (1). After ${\it\tau}_{1}=60.0~\text{s}$ this pattern is stabilised (1b). Then $B$ is switched to a value $B_{2}=10.74~\text{mT}$ before the fold of the hexagons and the magnetic fluid is allowed to relax from (2a) towards the flat state. Before the magnetic fluid flattens completely, $B$ is switched at (2b) to a value $B_{3}$ in the neighbourhood (3a) of the unstable branch of the imperfect transcritical bifurcation (cf. equations (1) and (2) in Gollwitzer et al. Reference Gollwitzer, Rehberg and Richter2010), marked by a dashed line and shaded region in figure 4. Now the pattern relaxes on path (3a) to (3b), where the height of the patch finally settles down (equilibrates) to a value approximately 1.5 mm higher (depending on the size of the patch, i.e. for larger patches, the height tends to that for the domain-covering hexagons) than the domain-covering hexagons. A movie (available at http://dx.doi.org/10.1017/jfm.2015.565) shows the birth of a three-spike patch. Each patch was at least stable for 120 s, as recorded by the measuring programme. Moreover, we have checked the long-term stability of exemplary patches after generation and observed no decay for more than 2 h.

Starting at the edge of the bistable region defined by the fold of the domain-covering cellular hexagon spikes and increasing $B_{3}$ , we obtain a sequence of localised states ranging from solitary spikes to patches containing multiple spikes/hexagon cells on an otherwise flat fluid surface with two to 24 cells, as shown in figure 5 for ${\it\tau}_{2}=10.0~\text{s}$ . The height of the spikes/cells of a patch are not uniform; we observe that spikes near the centre of the patch are higher than those on the edge of the patch. The distance between spikes within a patch is approximately the same as the spatial period of the domain-covering hexagon pattern. Remarkably, the spikes do not spread out in space as time increases despite their magnetic dipole–dipole repulsion, originating from the fact that the spikes are magnetised in parallel. In this way the patches are not an ensemble of weakly interacting solitary spikes in the sense of Richter & Barashenkov (Reference Richter and Barashenkov2005).

Figure 5. Localised states for ${\it\tau}_{2}=10.0~\text{s}$ at (a)  $B=11.0725~\text{mT}$ , (b)  $B=11.0857~\text{mT}$ , (c)  $B=11.1120~\text{mT}$ , (d)  $B=11.1648~\text{mT}$ , (e)  $B=11.1780~\text{mT}$ , (f)  $B=11.1911~\text{mT}$ , (g)  $B=11.2043~\text{mT}$ , (h)  $B=11.2175~\text{mT}$ and (i)  $B=11.2307~\text{mT}$ .

Figure 6. The number of spikes versus $B_{3}$ for different delays ${\it\tau}_{2}=1.0~\text{s}$ (♢), 7.5 s (▵), 8.5 s (○), 9.5 s (▫) and 10.0 s (▿) exhibits plateaux. The plateaux at 3, 10, 12, 14 and 21 spikes, where patches are forming complete patches, are marked by horizontal lines.

Moreover, it is interesting to note the mechanism by which hexagon cells are added to large patches as $B_{3}$ is increased. Going from figure 5(f) to (i), we see that spikes are added in the middle of the lower front in order to form a completed patch with dihedral symmetry $\mathbb{D}_{2}$ . It is also possible to see completed patches with $\mathbb{D}_{6}$ symmetry, as shown in figure 1.

Similar patches have been uncovered in the 2D Swift–Hohenberg equation as a signature of homoclinic snaking (cf. Lloyd et al. Reference Lloyd, Sandstede, Avitabile and Champneys2008; Escaff & Descalzi Reference Escaff and Descalzi2009; Kozyreff & Chapman Reference Kozyreff and Chapman2013). Also the mechanisms governing the growth of a patch, as observed above, resemble the ones described by these authors.

As shown in figure 6, the number of spikes in a patch increases monotonically with $B_{3}$ , with clearly defined plateaux. The data related with figure 5 are marked by ▿ and have been recorded with a delay of ${\it\tau}_{2}=10.0~\text{s}$ . Reducing ${\it\tau}_{2}$ to 9.5 s (▫), 8.5 s (○) and 7.5 s (▵) shifts the curves to lower magnetic inductions: because the starting amplitudes are higher, a lower $B_{3}$ is sufficient to generate a localised patch. Most importantly the prominent plateaux at $3,10,12,14,\ldots$ cells exist for all these values of ${\it\tau}_{2}$ . For the very short delay of ${\it\tau}_{2}=1.0~\text{s}$ , the curve (♢) hardly exhibits plateaux (apart from the plateau at 27 spikes, which is the maximal size of a hexagonal pattern in the set-up). That is because, for ${\it\tau}_{2}=1.0~\text{s}$ , the regular hexagonal pattern (see 1b in figure 4) has hardly decayed, and we cannot sensitively explore the phase space where homoclinic snaking occurs.

5. Governing equations

In order to elucidate the snaking process, we numerically solve the full three-dimensional (3D) dynamical equations. The field in the fluid is modelled by Maxwell’s equations with no free currents and a Young–Laplace equation for the free-surface. Since there are no free currents, the magnetic field in the fluid, $\boldsymbol{H}_{F}$ , can be written in terms of a gradient field $\boldsymbol{H}_{F}=h\hat{\boldsymbol{z}}-\boldsymbol{{\rm\nabla}}{\it\phi}$ , where $h$ is the magnitude of the applied uniform magnetic field strength in the fluid in the vertical direction $\hat{\boldsymbol{z}}$ and the non-uniform part of the field is described by $-\boldsymbol{{\rm\nabla}}{\it\phi}$ . The magnetisation is given by the linear equation $M(H_{F})=({\it\mu}_{r}-1)H_{\mathit{F}}$ , where $H_{\mathit{F}}=|\boldsymbol{H}_{\mathit{F}}|$ and ${\it\mu}_{r}=\text{const.}$ is the magnetic permeability of the magnetic fluid relative to that of free space, ${\it\mu}_{0}$ . Thus the magnetic induction follows from the linear relation $\boldsymbol{B}={\it\mu}_{r}{\it\mu}_{0}\boldsymbol{H}_{F}$ .

Of course, a linear magnetisation law can only be a rough approximation of the real $M(H_{\mathit{F}})$ measured in the experiment (§ 3). Despite that, previous analytical studies, which are based on a linear $M(H_{\mathit{F}})$ , can successfully predict all regular planforms on the ferrofluids and their underlying bifurcations (see Gailitis Reference Gailitis1977; Friedrichs & Engel Reference Friedrichs and Engel2001; Friedrichs Reference Friedrichs2002; Bohlius, Brand & Pleiner Reference Bohlius, Brand and Pleiner2011). Likewise in our numerical study we aim to catch the decisive qualitative features seen in the experiment, and to elucidate whether they are connected to homoclinic snaking.

The magnetic field in the air above the magnetic fluid is given by $\boldsymbol{H}_{A}={\it\mu}_{r}h\hat{\boldsymbol{z}}-\boldsymbol{{\rm\nabla}}{\it\chi}$ and the steady-state free-surface is given by $z={\it\zeta}(x,y)$ , with ${\it\zeta}=0$ corresponding to the flat state. The equations are non-dimensionalised using the critical wavenumber $k_{c}=\sqrt{{\it\rho}g/{\it\sigma}}$ for the length scale and $\sqrt{{\it\rho}g{\it\sigma}}$ for the unit of pressure, where ${\it\rho}$ is the fluid density and ${\it\sigma}$ is the surface tension. Furthermore, the amplitude $h$ of the uniform field strength inside the magnetic fluid is rescaled using ${\it\mu}_{r}$ ,

such that the instability occurs at $\mathscr{H}_{\mathit{c}}=1$ . We note that, when comparing (5.1) with the related equation (5.5) of Silber & Knobloch (Reference Silber and Knobloch1988), we have rescaled with respect to ${\it\mu}_{0}$ . The bifurcation parameter is then

We also note that ${\it\epsilon}$ depends on ${\it\mu}_{r}$ . For a static magnetic fluid in a vertical field, the dimensionless equations are

where $D$ is the depth of the fluid and $p^{\ast }$ is an effective pressure in the magnetic fluid. Equation (5.3a ) is the stationary Navier–Stokes equation for the magnetic fluid, and (5.3b ) and (5.3c ) are derived from Maxwell’s equations assuming a linear magnetisation law (see Cowley & Rosensweig Reference Cowley and Rosensweig1967). We consider a domain above the fluid to a height $D$ for computational ease. At the interface $z={\it\zeta}$ , the boundary conditions are given by

where $[\cdot ]_{-}^{+}$ is the value of the jump across the interface of the quantity in the brackets, $\hat{\boldsymbol{n}}$ is the unit normal to the interface and ${\it\kappa}$ is the mean curvature of the surface (see Silber & Knobloch Reference Silber and Knobloch1988):

No-flux boundary conditions are imposed on the top and bottom of the domain, i.e. ${\it\phi}_{z}(\boldsymbol{x},-D)={\it\chi}_{z}(\boldsymbol{x},+D)=0$ where $\boldsymbol{x}=(x,y)$ . Solving (5.3a ) for the pressure, within a constant, $C$ , at the free-surface yields

The pressure can then be eliminated and the boundary conditions written in terms of the fields ${\it\zeta}(x,y),{\it\phi}(x,y,z)$ and ${\it\chi}(x,y,z)$ , as well as two parameters, ${\it\mu}_{r}$ and ${\it\epsilon}$ .

The system of equations possesses a free energy given by

made up of magnetic, hydrostatic and surface energies, respectively, with conservation of mass, with boundary conditions at the top and bottom of the domain, and with the boundary condition at the interface ${\it\phi}(\boldsymbol{x},{\it\zeta})-h{\it\zeta}(1-{\it\mu}_{r})={\it\chi}(\boldsymbol{x},{\it\zeta})$ as constraints, where ${\it\Omega}\subset \mathbb{R}^{2}$ is the space domain of ${\it\zeta}(\boldsymbol{x})$ and $\boldsymbol{x}=(x,y)$ . We note that Gailitis (Reference Gailitis1977) and Bohlius, Pleiner & Brand (Reference Bohlius, Pleiner and Brand2006) use other energy formulations without the boundary conditions and conservation of mass. The constant $C$ is a Lagrange multiplier for the mass constraint: $\int _{{\it\Omega}}{\it\zeta}\,\text{d}\boldsymbol{x}$ equal to a constant. Taking variations of ${\mathcal{E}}$ with respect to the free-surface ${\it\zeta}$ yields the interface equation (5.4a ), and taking variations with respect to ${\it\phi}$ and ${\it\chi}$ , using the constraint, yields the Laplace equations (5.3b ) and (5.3c ) and the no-flux boundary conditions ${\it\phi}_{z}(\boldsymbol{x},-D)={\it\chi}_{z}(\boldsymbol{x},+D)=0$ . The advantage of this free energy is that it represents the whole system and we can now numerically compute energies taking into account the magnetic fields.

In order to fix the pressure constant, $C$ , we assume that the base state $({\it\chi},{\it\phi},{\it\zeta})=(0,0,0)$ solves the equations corresponding to the mass constraint $\int _{{\it\Omega}}{\it\zeta}\,\text{d}\boldsymbol{x}=0$ . We choose to study the case of large depth of magnetic fluid and we set the depth, $D$ , to be 10. For ${\it\mu}_{r}=2$ , this yields a depth of 16.42 mm, i.e. about three times the depth of the fluid in the experiment, which is sufficiently large that the domain-covering hexagons do not change as the depth is increased. We have ignored the effect of the shallow layer in our numerics since one would have to solve for an additional magnetic field below the magnetic fluid and apply suitable compatibility conditions.

Carrying out a coordinate transformation to flatten out the free-surface (Twombly & Thomas Reference Twombly and Thomas1983) given by

and

we can perform a Legendre transformation in order to define a Hamiltonian. Arbitrarily taking the $x$ direction, we define the momenta coordinates ${\it\alpha}=\partial {\mathcal{E}}/\partial {\it\phi}_{x},{\it\beta}=\partial {\mathcal{E}}/\partial {\it\chi}_{x}$ and ${\it\eta}=\partial {\mathcal{E}}/\partial {\it\zeta}_{x}$ . Then the Hamiltonian for the system is given by

on the constraint manifold ${\mathcal{M}}=\{{\it\phi}(\boldsymbol{x},0)-h{\it\zeta}(1-{\it\mu}_{r})-{\it\chi}(\boldsymbol{x},0)=0,{\it\phi}_{x}(\boldsymbol{x},0)-h{\it\zeta}_{x}(1-{\it\mu}_{r})-{\it\chi}_{x}(\boldsymbol{x},0)=0\}$ ; see Groves & Toland (Reference Groves and Toland1997) and Groves, Lloyd & Stylianou (Reference Groves, Lloyd and Stylianou2015) for the similar case of water waves.

The Hamiltonian (5.10) is a constant along solutions in the $x$ direction and provides a wavelength selection principle (see Lloyd et al. (Reference Lloyd, Sandstede, Avitabile and Champneys2008) for the same principle applied to the Swift–Hohenberg equation). In particular, if we evaluate ${\mathcal{H}}$ along any localised planar (hexagon) front connecting to the trivial state ${\it\zeta}=0$ and the domain-covering hexagon pattern, we find that ${\mathcal{H}}={\mathcal{H}}(\mathbf{0})$ along the limiting domain-covering hexagon pattern so that it belongs to the same level set. Generically, ${\mathcal{H}}$ will only equal ${\mathcal{H}}(\mathbf{0})$ at finitely many spatially periodic orbits in the family of spatially periodic patterns and will therefore serve as a selection principle that involves only the spatially periodic patterns; see Beck et al. (Reference Beck, Knobloch, Lloyd, Sandstede and Wagenknecht2009) for more details.

6. Numerical results

The 3D equations (5.3b,c ) and (5.4a ) are solved by first flattening out the free-surface using (5.8) and (5.9) so that the computational domains are fixed cuboids. We apply Neumann boundary conditions on the four sizes of ${\it\Omega}$ in order to remove the translational invariance in $x$ and $y$ that would cause problems with Newton solvers. This choice of boundary conditions on ${\it\Omega}$ limits us to cellular hexagons, symmetric planar fronts and $\mathbb{D}_{2}$ symmetric patches involving cellular hexagons computed on the positive quadrant; see Lloyd et al. (Reference Lloyd, Sandstede, Avitabile and Champneys2008) for more details. To remove the continuous symmetry due to the magnetic potentials being defined only up to a constant, we add an extra parameter ${\it\lambda}$ to the Laplace equations (5.3b,c ) to be solved along with ${\it\phi},{\it\chi}$ and ${\it\zeta}$ and we append an extra equation (see Twombly & Thomas Reference Twombly and Thomas1983)

We now have a well-posed boundary value problem with all the continuous symmetries removed, allowing us to use Newton methods. We discretise in MATLAB using fourth-order central finite differences (for planar fronts) or a pseudo-spectral Fourier method (for hexagon spikes and patches) in $x$ and $y$ directions and Chebyshev pseudo-spectral method in $z$ with a Newton–GMRES method (Kelley Reference Kelley2003) and secant pseudo-arclength continuation (Krauskopf, Osinga & Galan-Vioque Reference Krauskopf, Osinga and Galan-Vioque2007). For the Newton–GMRES scheme, one needs a preconditioner. For this we precompute a sparse incomplete LU-decomposition (with the Crout algorithm Saad (Reference Saad1996)) of the linear system with a flat interface at onset and this is fixed for the continuation. The incomplete LU-decomposition is computed with a drop tolerance of $1\times 10^{-6}$ and we set a maximum of 20 iterations for GMRES with 40 Arnoldi restarts and a tolerance of $10^{-4}$ . Typical spatial step sizes in $x$ and $y$ directions are 0.25–0.5 for the Fourier pseudo-spectral method and 20 Chebyshev collocation points in $z$ and a fixed arclength step size of 0.01. The numerics are validated by comparing with normal-form results by Silber & Knobloch (Reference Silber and Knobloch1988). In particular, we check that the change of sub/supercriticality of rolls and squares occurs at the theoretically predicted values of ${\it\mu}_{r}$ ; see the appendix for more details. Fold loci are found by carrying out parameter slices and then doing a linear interpolation about the fold values. The Maxwell point is computed by appending the Hamiltonian and energy conditions to the system.

In the following we present the results obtained with the procedures sketched above. We start with the hexagon Maxwell curve (§ 6.1), which is succeeded by hexagon planar fronts (§ 6.2), hexagon patches (§ 6.3) and rhomboid patches (§ 6.4).

7. Comparison between experimental and numerical results

We now turn to comparing our experimental and numerical results. As stated in § 5, our model is based on a linear magnetisation law, whereas we presented in § 3 the nonlinear $M(H_{\mathit{F}})$ of the real ferrofluid, with an effective permeability ranging from ${\it\mu}_{\mathit{eff}}=4.57$ at $B=0~\text{mT}$ down to 3.2 at $B_{c}=10.8~\text{mT}$ . Therefore, a direct matching of all experimental and numerical results cannot be expected. For ${\it\mu}_{r}\rightarrow 1$ , the linear magnetisation law becomes a more accurate approximation of the nonlinear magnetisation law and the model predicts that homoclinic snaking should persist, albeit with an exponentially small width in which the snaking exists.

In order to provide a comparison with the experimental results in § 4, we instead chose a magnetic permeability ${\it\mu}_{r}=\text{const.}$ for our model, which captures the qualitative features of the experiment. In a first attempt, we chose ${\it\mu}_{r}$ to be the effective permeability at the critical point ${\it\mu}_{\mathit{eff}}(H_{c})=3.2$ (cf. figure 3 b) for the model: ${\it\mu}_{r}=3.2$ . In figure 12(a), we path-follow the two-spike patch, as depicted in figure 12(b), bifurcating from the flat state to the upper branch. As seen from figure 12(a), for ${\it\mu}_{r}=3.2$ , localised structures exist from $B=10.3$ to 11.5 mT, which is a significantly larger range (about 200 %) than that observed in figure 6. Moreover, the height of the spikes is about 8 % larger than that observed experimentally.

Figure 12. (a) Numerical bifurcation diagram for ${\it\mu}_{r}=3.2$ for a two-spike configuration, depicted in (b), that has bifurcated off the flat state. The state in (b) was computed for $B=11.3~\text{mT}$ . We terminated the branch early.

In a second attempt, we select a smaller value of magnetic permeability, namely ${\it\mu}_{r}=2$ . Here we find that, while the critical value of the induction is far from that experimentally observed, the width of the bistability region between the flat state and the domain-covering hexagons is very close (approximately 0.25 mT in both experiments and numerics); see figure 13. Furthermore, we find from the numerics at the fold of the domain-covering hexagons that the predicted height of the hexagons is 1.2 mm (from figure 4 we find the measured height of 1.3 mm) and at onset the height of the stable hexagons is predicted to be 2.4 mm (and from figure 4 we find the measured height of 2.4 mm), i.e. in the bistability region the numerics are within 5 % error of that experimentally observed. Hence, it is possible that phenomenologically ${\it\mu}_{r}=2$ is the appropriate value to take in the linear magnetisation law for the numerics for localised hexagon patches.

Figure 13. Maximum height of the domain-covering pattern (black line) computed numerically for ${\it\mu}_{r}=2$ at the critical wavenumber $k_{\mathit{c}}\approx 608.94~\text{m}^{-1}$ versus the magnetic induction. The maximum height of the rhomboid patches is marked by a brown line.

The heights of the spikes in the patches are found numerically to be highest at the core of the patch and decrease as the spikes get closer to the edge of the patch. As shown in figure 13(a), the maximum height of the rhomboid patches is slightly larger than the domain-covering hexagons, similar to that seen in the experiments.

The numerical observation of various coexisting hexagon patches in figure 9, in the bistability region of the flat state and the domain-covering hexagons, agrees well with the experimental results shown in figure 6, where for the same $B_{3}$ values several different patches were observed. In particular, we are able to find numerically patches (figure 10 a–d) that look to be identical to those seen in figure 5(b,d,f,i). Our numerical methods imposed that the patches should have $\mathbb{D}_{2}$ symmetry; hence we were unable to locate any of the other patches seen in figure 5. However, owing to the existence of a Maxwell point nearby, we believe hexagon spikes/cells can be added to the patches to break the $\mathbb{D}_{2}$ symmetry provided you are in the $\langle 11\rangle$ hexagon front snaking region. Hence, we have a possible explanation for the existence of all the various types of patches observed experimentally.

Perhaps the largest difference between the numerical and experimental results is that the sequence of plateaux for fixed ${\it\tau}_{2}$ suggests that one can climb up the snakes for increasing $B_{3}$ , i.e. the snakes are not vertical as in figures 9 or 10 but slanted, as found by Firth, Columbo & Scroggie (Reference Firth, Columbo and Scroggie2007) and Dawes (Reference Dawes2008) for a conserved quantity. Slanted snaking may emerge in the presence of a nonlinear magnetisation law and finite-domain effects, and could provide an explanation for the emergence of larger and larger patches seen in figure 10.

8. Conclusion and outlook

Utilising a magnetic pulse sequence, we prepare a ferrofluidic layer to be situated next to the unstable branch of the Rosensweig instability. There, localised hexagon patches emerge spontaneously. The measured sequence of patches indicates homoclinic snaking. This is corroborated by numerics, where we investigate the Young–Laplace equation coupled to the Maxwell equations. Homoclinic snaking is demonstrated for two types of planar hexagon fronts as well as hexagon and rhomboid patches. The experimentally observed sequence of patterns can be interpreted as an imperfect mixture of those prototype scenarios. The numerically unveiled Maxwell point of the hexagons provides an energy argument for those various regular and irregular patches observed, since near the Maxwell point the energy of a single hexagon spike is the same as the flat state and there is no energy penalty in adding/removing hexagon spikes. The width of the experimentally observed plateaux is about 0.5 % of the control parameter, which is comparable to the snaking interval of the numerical prototypes. Moreover, the height of the patches in experiment and numerics agrees within 10 %. Therefore, we believe that our model captures all the most important phenomena of the system. In particular, as the relative magnetic permeability tends to unity, the validity of our model increases. It thus unveils an organising centre for the emergence of localised hexagon patches seen in the experiment.

Figure 14. (a) Bifurcation diagram of the domain-covering hexagons for ${\it\mu}_{r}=2$ on the domain $(x,y)\in [-2{\rm\pi},2{\rm\pi}]\times [-2{\rm\pi}/\sqrt{3},2{\rm\pi}/\sqrt{3}]$ with $N_{x}=N_{y}=20$ , $D=10$ and varying $N_{z}$ . (b) Bifurcation diagram of the domain-covering hexagons for ${\it\mu}_{r}=2$ on the domain $(x,y)\in [-2{\rm\pi},2{\rm\pi}]\times [-2{\rm\pi}/\sqrt{3},2{\rm\pi}/\sqrt{3}]$ with $N_{z}=20$ and varying $N_{x}$ ( $N_{y}=N_{x}$ ). (c) Relative error on a semi-log scale of the L2 norm of the interface $\Vert {\it\zeta}\Vert _{2}^{2}$ for the domain-covering hexagons on the upper branch at $({\it\epsilon},{\it\mu}_{r})=(0,2)$ with $N_{z}=20,D=10$ and varying $N_{x}=N_{y}$ with respect to the hexagon computed for $N_{z}=22,N_{x}=N_{y}=22$ . (d) Relative error on a semi-log scale of the L2 norm of the interface $\Vert {\it\zeta}\Vert _{2}^{2}$ for the domain-covering hexagons on the upper branch at $({\it\epsilon},{\it\mu}_{r})=(0,2)$ with $N_{x}=N_{y}=20,D=10$ and varying $N_{z}$ with respect to the hexagon computed for $N_{z}=22,N_{x}=N_{y}=22$ . From both (c) and (d) we observe rapid (geometric) convergence due to the spectral accuracy of the numerical pseudo-spectral methods.

It still remains to explore how various stable configurations of patches connect via unstable branches both experimentally (by quasi-statically varying the magnetic induction $B$ once a patch has been generated) and theoretically. The discrepancy between the experimental indications for slanted snaking and the vertical snaking obtained numerically may be solved by incorporating a nonlinear magnetisation law, finite-size and shallow-layer effects or by carrying out more experiments to determine the precise extent of the plateaux.

Figure 15. Normal-form results of Silber & Knobloch (Reference Silber and Knobloch1988) showing the changes in sub- and supercriticality of the bifurcations for rolls, squares and up-hexagons, i.e. hexagons whose maximum amplitude is positive.

Theoretically, there are now some very interesting directions one could take. With the existence of an energy functional for the full system that incorporates the boundary conditions, one could look at applying variational and spatial-dynamical techniques that have been developed in the water–wave literature (see e.g. Groves & Sun Reference Groves and Sun2008; Buffoni et al. Reference Buffoni, Groves, Sun and Wahlén2013).