Localised Radial Patterns on the Free Surface of a Ferrofluid

This paper investigates the existence of localised axisymmetric (radial) patterns on the surface of a ferrofluid in the presence of a uniform vertical magnetic field. We formally investigate all possible small-amplitude solutions which remain bounded close to the pattern's centre (the core region) and decay exponentially away from the pattern's centre (the far-field region). The results are presented for a finite-depth, infinite expanse of ferrofluid equipped with a linear magnetisation law. These patterns bifurcate at the Rosensweig instability, where the applied magnetic field strength reaches a critical threshold. Techniques for finding localised solutions to a non-autonomous PDE system are established; solutions are decomposed onto a basis which is independent of the radius, reducing the problem to an infinite set of nonlinear, non-autonomous ODEs. Using radial centre manifold theory, local manifolds of small-amplitude solutions are constructed in the core and far-field regions, respectively. Finally, using geometric blow-up coordinates, we match the core and far-field manifolds; any solution that lies on this intersection is a localised radial pattern. Three distinct classes of stationary radial solutions are found: spot A and spot B solutions, which are equipped with two different amplitude scaling laws and achieve their maximum amplitudes at the core, and ring solutions, which achieve their maximum amplitudes away from the core. These solutions correspond exactly to the classes of localised radial solutions found for the Swift-Hohenberg equation. Different values of the linear magnetisation and depth of the ferrofluid are investigated and parameter regions in which the various localised radial solutions emerge are identified. The approach taken in this paper outlines a route to rigorously establishing the existence of axisymmetric localised patterns in the future.


Introduction
In an experiment, a quiescent layer of ferrofluid in a flat, horizontal container is subjected to a uniform vertical magnetic field of magnitude h, see Figure 1a). As the applied magnetic induction exceeds a critical value, h c , the homogeneous state becomes unstable to small perturbations and regular cellular spikes of patterns appear. The pattern of spikes most commonly forms in a hexagonal lattice. This phenomenon is known as the 'Rosensweig' instability, and has been a subject of interest since the 1960s (see Cowley & Rosensweig [12], Rosensweig [35]). In 2005, Richter & Barashenkov [34] were the first to show that they could induce solitary axisymmetric spikes, as seen in Figure 1b), using a local perturbation of the applied magnetic field. These radial 'spots' persisted once the local magnetic perturbation was removed and exhibited some notable properties. The spikes emerged away from the centre of the container, were able to be moved around the container, and did not appear to feel the effects of the container boundary. This suggests the spots are well localised, possibly decaying exponentially fast to the flat state. This experimental evidence was further supported by a collection of numerical results (see Lloyd et al. [27] for pseudo-spectral numerical methods and experimental results, and Lavrova et al. [26], Cao & Ding [8] for finite-element simulations). Despite this, there is no analytical theory for proving the existence of exponentially localised axisymmetric solutions to the ferrohydrostatic problem. Figure 1: a) Static localised patterns appear on the free surface between a non-magnetisable fluid (transparent) and a ferrofluid (shaded) as the strength of a vertical applied magnetic field is varied. b) Localised axisymmetric peaks, termed 'spot A' solutions, have been observed experimentally (See Richter [33], reproduced with permission).
Theoretical Approaches to Ferrofluids: The first attempt at an analytic study of domain-covering cellular patterns of the Rosensweig instability was by Gailitis [15], who substituted a cellular free-surface ansatz into a hypothetical free energy for the system (an infinite-depth ferrofluid with a linear magnetisation law). This resulted in finding regions of existence for various lattice patterns (squares, stripes and hexagons), and was later extended to finite-depth ferrofluids by Friedrichs & Engel [14].
Twombly & Thomas [40] considered static, doubly-periodic solutions of the ferrohydrostatic equations near the onset of instability, extending the previous work of Zaitsev & Shliomis [45] to a finite-depth ferrofluid. These results, close to the bifurcation point, were supplemented by other studies of two-dimensional periodic free-surfaces; see Silber & Knobloch [39] for normal-form analysis, Bohlius et al. [5,6] for deriving amplitude equations near onset and Groves & Horn [18,23] for a Dirichlet-Neumann formulation and local bifurcation theory. With the exception of the Groves & Horn work, all these studies included a linear magnetisation law.
The study of localised solutions to the ferrohydrostatic equations has had very little attention to date. Currently, the only rigorous result is a proof of existence for one-dimensional localised planar waves by Groves et al. [19]. This work examined a finite-depth ferrofluid with a fully nonlinear magnetisation law, using variational methods to formulate a spatial Hamiltonian system. Using a combination of centre-manifold reduction techniques and normal-form theory, they were then able to prove the existence of spatially localised one-dimensional solutions. Notably, the nonlinearity of the magnetisation law does not qualitatively affect any of the results, and so it was concluded that a linear magnetisation law is a reasonable assumption for analytical studies.
Localised Radial Patterns: The ferrofluid problem strongly resembles that of gravity-capillary waves [7,24]; as a result, water wave methodologies provide a robust framework to study ferrofluids, notably in the case of ferrofluid jets [4,20]. However, localised solutions in water waves are created via the introduction of a uniform flow, creating an inherent symmetry-breaking that makes localised radial solutions impossible. Hence, new techniques must be developed in order to find localised radial solutions in the ferrofluid problem.
Recently, there has been some progress in the study of localised radial patterns in reaction diffusion systems; for example, in 2003 Scheel [38] introduced radial centre manifold and normal-form theory for stationary and time-periodic axisymmetric patterns in general reaction diffusion systems. Prototypical pattern forming systems have been analysed to prove the existence of localised radial solutions to the Swift-Hohenberg equation [28][29][30], the complex Ginzburg-Landau equation [31], and others [13]. However, the non-autonomous nature Figure 2: a) The free surface of a localised radial pattern lies on the intersection of solutions that are bounded as r → 0, called the 'core manifold' W cu − , and solutions that decay exponentially as r → ∞, called the 'far-field manifold' W s + . b) The far-field manifold is constructed as perturbations from the decaying solutions that lie on the autonomous centre manifold W c + . The mappings from W c + to W s + are derived from foliation theory, as seen in [31].
of radial patterns complicates any theoretical understanding one may have for two-or three-dimensional patterns (as discussed in [25]). In particular, the standard centre manifold reduction procedure [21,32,42] for small amplitude patterns does not apply. One instead has to construct two sets of small amplitude solutions (one set where solutions remain bounded as r → 0 which we call the core manifold, W cu − , and the other set where solutions decay exponentially to the flat state as r → ∞, which we call the far-field manifold, W s + ) and look for intersections of these two sets, see Figure 2(a), which we describe in more detail below. For each problem, these sets must be constructed on a case-by-case basis.
The present contribution is strongly motivated by the results found for the Swift-Hohenberg equation by Lloyd & Sandstede [28], who proved the existence of localised rings and elevated spots using radial centremanifold and normal-form theory of Scheel [38]. Following this, McCalla & Sandstede [30] employed geometric blow-up methods to prove the existence of a class of localised depressed spots, called spot B.
Overview: In the present paper we formally identify all possible small-amplitude localised radial patterns that bifurcate from the flat state in the ferrohydrostatic problem, providing a framework that guides future rigorous results on the existence of localised axisymmetric patterns. The expected properties of a localised radial solution u(x), where x = (r, θ, z) are cylindrical polar coordinates, are the following: 1. Axisymmetric: u(x) is independent of θ, i.e. u(x) = u(r, z).
3. Localised: u(r, z) → u 0 exponentially as r → ∞, where u 0 is the flat state of the system. Property 1 will inform our formulation of the radial ferrohydrostatic problem; we construct a free-surface problem with two immiscible fluids of equal arbitrary depth D separated by an interface at z = η(r), see Figure 1, with a constant magnetic field of strength h applied vertically. There exists a critical magnetic field strength, h c , such that for values of h > h c the flat quiescent state destabilises and patterns form, known as the Rosensweig instability. Hence, we define ε := h − h c to be the bifurcation parameter of the system; the bifurcation point should occur at ε = 0 and localised solutions should emerge for 0 < | ε| ≪ 1. The ferrofluid, taken to be the lower fluid, is assumed to have a linear magnetisation law, such that M = (µ − 1)H, where µ > 1 is the magnetic permeability of the ferrofluid. Using variational methods, we have formulated the radial ferrohydrostatic problem, for finite depth and a linear magnetisation law, as a quasilinear non-autonomous PDE of the form d dr u = L (r) u + F (u, ε, r), (1.1) Figure 3: A sketch of the matching process in the geometric blow-up coordinates. The centre coordinates of the farfield manifold are written in 'rescaling' coordinates, known as the 'rescaling chart', in order to identify solutions that decay exponentially as r → ∞. These solutions are then tracked back to the point r = δ0ε − 1 2 , where δ0 > 0 is fixed. In order to retain control of the point r = δ0ε − 1 2 as ε → 0, solutions are tracked in 'transition' coordinates, known as the 'transition chart', to the point r = r0, where they are matched with the core manifold.
where L(r) is a non-autonomous linear differential operator containing ∂ zz terms, ε is a rescaled bifurcation parameter such that ε ∝ ε, and F contains all the nonlinear terms of the system. The problem is nonautonomous due to the presence of 1/r terms that also create an apparent singularity at r = 0. However, in the limit as r → ∞, the system (1.1) reduces to the problem studied by Groves et al. [19]. Properties 2 & 3 will be satisfied for local regions of r, specifically r ∈ [0, r 0 ] and r ∈ [r ∞ , ∞), called the 'core' and 'far-field', respectively. Here r 0 , r ∞ > 0 are fixed constants, and r 0 is chosen to be larger than r ∞ such that these regions overlap, see Figure 2a). In order to analyse (1.1) in these regions, we decompose our solution onto a z-dependent eigenbasis, found from the spectrum of the linear operator L(r) in the limit r → ∞. The eigenvalues of the matrix L ∞ := lim r→∞ L(r) are the roots of the linear dispersion relation ∆(λ) := µ(µ − 1) 2 (µ + 1) λ sin 2 (λD) − λ 2 − ρ 0 gσ µ 2 h 4 sin(λD) cos(λD), (1.2) where ρ 0 , g, σ are the relative density of the two fluids, gravitational constant, and surface tension, respectively. The autonomous linear problem u r = L ∞ u, related to (1.1) as r → ∞, exhibits a Hamiltonian-Hopf bifurcation at h = h c . Then, at the bifurcation point ε = 0, a complex number λ ∈ C is an eigenvalue of L ∞ if and only if λ ∈ {±ik} ∪ {±λ n } n∈N , where k, λ n ∈ R, for n ∈ N. The pair of imaginary eigenvalues {±ik} have double algebraic multiplicity, and we will treat the wavenumber k as a parameter. By projecting onto the respective eigenvectors for each eigenvalue, the PDE system reduces to an infinite set of non-autonomous ordinary differential equations. Writing the full solution in this 'spectral' decomposition, we can then solve the system in both the 'core' and 'far-field' regions, subject to Properties 2 and 3, respectively.
We first construct the set of all small-amplitude solutions that remain bounded as r → 0, which we call the 'core' manifold and denote as W cu − , using notation from dynamical systems theory for a centre-unstable manifold. We note that this is a submanifold of the standard local centre-unstable manifold W cu − , containing the set of all solutions that grow at most algebraically as r → 0. Then, we construct the set of all solutions that exponentially decay to the flat state as r → ∞, called the 'far-field' manifold and denoted as W s + , using the notation for a standard stable manifold. In order to perform our intended analysis, we parametrise the 'far-field' manifold with respect to the decaying elements on an autonomous centre manifold W c + , see Figure  2b), using the theory of foliations seen in [31] for the Ginzburg-Landau equation.
Next we need to isolate the decaying solutions of W c + , and so we introduce blow-up coordinates seen in [30]  Solutions have the following behaviour: they behave like the Bessel function J0(kr) in the core region (blue), maintain algebraic decay in the transition region (red), and exhibit exponential decay away from the core (yellow).
in order to restrict W c + to just its elements that decay as r → ∞. More specifically, we will track solutions through two coordinate charts, called the 'rescaling chart' and the 'transition chart'. By rescaling r = s/ √ ε, where ε is the bifurcation parameter of the system, we construct the 'rescaling chart' which allows us to establish a parametrisation for decaying solutions that is valid for all s ∈ [δ 0 , ∞) for fixed δ 0 > 0. At the bifurcation point, i.e. for ε = 0, we lose control of the point r = δ 0 / √ ε; we construct the 'transition chart' which allows us to maintain control of the solution as ε → 0, see Figure 3.
Once the 'core' and 'far-field' manifolds have both been constructed, we evaluate both at the large, but finite, value r = r 0 and apply asymptotic matching to identify any intersections between the two. Any solution that lives on the intersection of the 'core' and 'far-field' manifolds satisfies properties 1, 2 & 3 and so is a small-amplitude localised radial pattern.
Through this procedure, we identify four types of localised radial patterns; up-spots, down-spots, up-rings, and down-rings. Here, the prefix 'up-' or 'down-' indicate whether the maximum point on the free surface is an elevation or a depression, and the class 'spot' or 'ring' indicates whether that maximum is at, or away from the core respectively.
Main Results: The primary results of this work are related as follows; firstly, we formally show the existence of up-spot solutions, termed 'spot A' in [30], that were identified experimentally in [34]: Existence of Spot A: Fix µ > 1, δ 0 > 0, and (D, k) such that then for 0 < ε ≪ 1, the ferrohydrostatic equation (1.1) has a stationary localised axisymmetric solution u(r, z). These solutions remain close to the trivial state u ≡ 0, and, for fixed r 0 > 0, the profile of the height Figure 5: a) The sign of c3 is plotted for the scaled wavenumber kD := kD, and magnetic permeability M0 := µ−1 µ+1 , where c3 < 0 only in the shaded regions. The lower shaded region, corresponding to shallow-depth solutions, is unphysical in our current model, and so we disregard this. The heteroclinic orbit q+(s) exists in the shaded regions, i.e. when c3 < 0. b) A sketch of the envelope function q(s), which exhibits algebraic growth for small s and exponential decay for large s. Any solution u(s) that lies close to the connecting orbit q+(s) is bounded by q(s).
We note that spot A does not require stripes to bifurcate subcritically, only that ν > 0. The restriction ν > 0 can be solved numerically; we find a range of values for k given a fixed D, see Figure 4a) for a plot of the parameter space where spot A solutions emerge. The solution has its maximum at r = 0 and behaves like the Bessel function J 0 (kr) near to the core, as depicted in Figure 4b). While we anticipate the spot A to be initially unstable as it bifurcates subcritically, we expect it to undergo a fold where it restabilises and it is this spot that we believe is observed experimentally.
We also show the existence of a class of up-ring and down-ring solutions, where the pattern's maximum exists away from their centre. These solutions, termed 'rings', have been found in the Swift-Hohenberg equation [28] with a scaling of ε 3 4 . To establish the existence of ring solutions, we first introduce re-scaled parameters k D := kD and M 0 := µ−1 µ+1 ; we find a restriction on k D > 0 in terms of a fixed 0 ≤ M 0 < 1 such that ring solutions emerge from the flat state.
As before, we can numerically plot the parameter space in which ring solutions emerge from the flat state, as shown in Figure 5a). The restriction c 3 < 0 is identical to the one-dimensional problem, where c 3 < 0 is the condition for homoclinic solutions to emerge. These solutions are restricted by a homoclinic envelope function q( √ εr), which exponentially decays as r → ∞; an illustration of this function is depicted in Figure   5b), while the free-surface profile of a ring solution is plotted in Figure 6a). The ring solutions identified in this paper are closely related to the one-dimensional homoclinic solutions in [19], corresponding to planar waves in the three-dimensional ferrohydrostatic problem.
The two previous classes of solutions, spot A and rings, were found for the prototypical Swift-Hohenberg by Lloyd & Sandstede in 2009 [28], and have 'standard' scaling laws ε  [30] also showed the existence of a down-spot solution, with a depression at its centre. This solution, termed 'spot B', was found to have an unexpected scaling of ε  axisymmetric solution u(r, z). These solutions remain close to the trivial state u ≡ 0, exhibit exponential decay as r → ∞, and, for fixed r 0 > 0, the profile of the height of the interface is as ε → 0 uniformly, where m, b D > 0, ν, and q(ε 1 2 r) are defined in (8.4), (4.7), and (6.13), respectively, and δ 1 , δ 2 > 0 are appropriate small constants.
The free-surface profile of this spot B solution is plotted in Figure 6b). By taking a fixed depth D > 0, we can plot the restriction of the rescaled wavenumber k D , as seen in Figure 5a), as a restriction of the wavenumber k in terms of the rescaled magnetic permeability M 0 . Then, we can combine Figures 4a) & 5a) into one plot for the existence regions of all localised radial solutions in (k, M 0 )-parameter space. This is shown in Figure 7a) for a finite depth D > 0 and in Figure 7b) for the limit as D → ∞. For sufficiently large fixed D, spot B and ring solutions can only occur for M 0 ≥ M c ≈ 0.56, or equivalently µ ≥ µ c ≈ 3.54. However, spot A and spot B solutions require that k > k c = 1/6, meaning that spot A solutions emerge for a larger range of parameters than other localised radial patterns; see Figure 7b).
Outline: The outline of the rest of the paper is as follows.
Section 2 (Formulation) We formulate the radial ferrohydrostatic problem with two fluids of finite depth; the ferrofluid lies below a non-magnetisable fluid, where the ferrofluid obeys a linear magnetisation law.
The key variables are axisymmetric magnetic potentials φ ± (r, z) of the upper and lower fluid respectively, as well as the height of the free surface η(r), such that the two fluids form an interface at z = η(r). Here, r and z are the radius and the height variables in cylindrical polar coordinates, respectively. The physical system is defined by the following equations: Maxwell's equations for both fluids with continuity conditions on the interface; and the steady Euler equations with a magnetic forcing term, solved on the interface in conjunction with a normal stress continuity condition. These equations can be found from variations of a 'free energy', or Lagrangian, L of the form where variations are taken with respect to φ − , φ + , and η, and subject to the constraint φ − = φ + at z = η(r).
Here, µ 0 is the magnetic permeability of free space, ∇ is the axisymmetric gradient operator, and the other parameters were defined previously for (1.2).
There are various problems making analytical methods difficult to apply: • The domains of the system, i.e. the upper and lower fluids, are non-uniform in r due to the interface at z = η(r).
• The magnetic potentials have translational symmetry: if (φ ± ) are solutions, then (φ ± + c) are also solutions, for any arbitrary constant c.
• The system has nonlinear boundary conditions, so standard invariant manifold theory cannot be applied in its current state.
To deal with these problems, multiple coordinate transformations are applied. We begin by non-dimensionalising the problem and applying a 'flattening' transformation often seen in water wave literature. We rescale the height variable z → z(r) with respect to the free surface η(r), such that each domain defines a uniform rectangle in z coordinates, see Figure 8, at the expense of increasing the system's nonlinearity. We also isolate the constant applied magnetic field by defining φ + = ψ + + µhz, φ − = ψ − + hz and perform a Legendre transformation to turn our second-order system for three variables into a first-order system for six variables. We remove the translational symmetry by defining the 'average' of ψ ± and removing it from the system; this imposes an extra constraint on the system and adds to the nonlinearity. We deal with the nonlinear boundary conditions by employing a coordinate transformation seen in [19,Section 4]. This transformation leaves the linear problem invariant, and the new coordinates still satisfy the constraints of the system, now equipped with linear boundary conditions. Then, we define u = ( ψ − , ψ + , η, α − , α + , γ) ⊺ , where ψ ± , η are the transformed variables related to the respective original variables ψ ± , η and the transformed variables α ± , γ are related to the Legendre conjugates of ψ ± , η, respectively. We can then write the full system in the form (1.1), with linear boundary conditions B 1 u = 0.
Section 3 (Decomposition) We want to be able to construct local solutions in both the 'core' and 'far-field' regions and match them at a fixed r = r 0 . Therefore, we would like to decompose solutions of (1.1) into a basis of z-dependent vectors with r-dependent amplitudes. This means that, in both the 'core' and 'far-field' regions, we can match the amplitudes without having to worry about our basis changing across regions. Hence, in Section 3, we construct a basis of eigenmodes for the matrix L ∞ := lim r→∞ L(r), where L(r) is the non-autonomous linear operator for (1.1). The eigenvalues of L ∞ are found from the dispersion relation (1.2), with corresponding eigenmodes that are found to form a complete z-dependent basis. Hence we can write, where · indicates complex conjugation, {e, f } correspond to the double mulitplicity eigenvalue ik, and {e ±n } correspond to the respective real eigenvalues ±λ n , for n ∈ N. Then, projecting onto each of these eigenmodes, we reduce our PDE system to infinitely-many nonlinear coupled ODEs, where the linear problem decouples into disjoint pairs, , d d r a n = λ n a n − 1 2r (a n − a −n ) , (1.6) Section 4 (Core Parametrisation) In this section, we discuss solutions near the core and construct a parametrisation for the core manifold, W cu − , which contains all small-amplitude solutions that are bounded as r → 0. To do this, we first find explicit solutions to the linear problem (1.6), noting that half the solutions are bounded as r → 0 and half are not. Taking a general solution to (1.6) with coefficients d 1 , d 2 , d 3 , d 4 , c 1,n , and c 2,n , we apply the method of variation of constants [43] to find a general solution for the nonlinear problem; here, d 1 , d 2 , and c 1,n are coefficients of the asymptotically stable solutions as r → ∞ and d 3 , d 4 , and c 2,n are coefficients of the unstable solutions as r → ∞, for all n ∈ N. Then, for d 3 = d 4 = c 2,n = 0 and sufficiently small constants d 1 , d 2 , and c 1,n for all n ∈ N, we define a parametrisation that captures all small-amplitude solutions of (1.1) with decomposition (1.5) that are bounded as r → 0.
We explicitly find the coefficient ν in front of d 2 2 in the parametrisation of the amplitude b(r 0 ); we will require this to be positive for spot solutions to emerge from the flat state, and we find sgn(ν) = sgn( ν), where ν is defined in (1.3).
Section 5 (Far-field Decomposition) In this section, we discuss solutions in the far-field, constructing a parametrisation for the far-field manifold W s + , which contains all small-amplitude solutions that decay exponentially as r → ∞.
We first extend the system by defining σ := 1 r , where 0 < σ ≪ 1 in the far-field region. This makes the system (1.6) autonomous, and so standard invariant manifold theory can be applied; we note, for the linear problem, σ = σ * forms an invariant subspace, where σ * is an arbitrary constant. We initially use the same variation of constants approach as in the core problem to construct a centre-stable manifold W cs + | σ=σ * for a fixed σ = σ * > 0, containing all small-amplitude solutions that grow at most algebraically as r → ∞. For sufficiently small solutions, the centre and stable submanifolds W c + , W s + ⊂ W cs + form a transverse intersection, i.e. they cover the entire manifold W cs + . Here, W c + is an autonomous centre manifold, containing all smallamplitude solutions that grow at most algebraically for all r, and W s + is the stable manifold containing all small-amplitude solutions that exponentially decay to zero as r → ∞.
In order to be able to match the far-field manifold W s + to the core manifold W cu − , we parametrise q ∈ W s + as a perturbation from a base point p ∈ W s c , where W s c ⊂ W c + is the set of all small-amplitude solutions on the centre manifold that exhibit exponential decay as r → ∞, such that q → p as r → ∞. Then, to fully parametrise the far-field manifold, we need to identify initial conditions for a point p ∈ W c + on the centre manifold such that p ∈ W s c , i.e. they decay exponentially to zero as r → ∞. We restrict the far-field system onto the centre manifold, applying a normal-form transformation to make the next analytical step tractable. Section 6 (Blow-up Coordinates) We formulate the geometric blow-up coordinates introduced in [30]. These coordinates are required to parametrise decaying solutions on the centre manifold. However, we also need to control the far-field parametrisation as r approaches r 0 at the bifurcation point ε = 0, reconciling the algebraic behaviour in the core with the exponential decay in the far-field.
To maintain control of the solutions as E → 0, we introduce a second transformation called the 'transition' chart [30,Section 2.4]. By scaling our variables by σ and introducing 'exponential time' ρ = ln r, the system can be seen to have two equilibria P ± . Finding explicit solutions to (6.5), we establish heteroclinic orbits q ± (s) that connect Q − in the rescaling chart to P ± in the transition chart, respectively. In the transition chart, the dynamics of the system is well-defined close to these equilibria, and so we will use P + and P − as a guide to track our solutions back to the point ρ = ln r 0 . Then, substituting the parametrisation of the centre coordinates back into the foliation parametrisation, we will be able to match the far-field manifold with the core manifold at the point, as seen in the next section. Section 7 (Matching) In the final section, we follow the matching procedure for each localised radial pattern: spot A, spot B and rings. For spot A solutions, we track solutions from Q − in the rescaling chart, staying close to the orbit q − (s), to some fixed point E = δ 0 > 0. Then, in the transition chart, solutions are tracked backwards in ρ, close to the equilibrium P − up to the point ρ = ln r 0 . For spot B solutions, we track solutions from Q − in the rescaling chart, staying close to the orbit q + (s), to E = δ 0 . Then, in the transition chart, solutions are tracked close to the equilibrium P + up to some point ρ = ρ 1 > ln r 0 . Following this, solutions transfer from the neighbourhood of P + to the neighbourhood of P − , and then solutions are tracked close to P − up to the point ρ = ln r 0 . For ring solutions, we track solutions from Q − in the rescaling chart, staying close to the orbit q + (s), to E = δ 0 . Then, in the transition chart, solutions are tracked close to the equilibrium P + up to the point ρ = ln r 0 . In each case, the full far-field parametrisation is then matched with the core manifold at the point r = r 0 , to find the full profile of each localised radial pattern.
Novelty: The novelty of this work can be identified as the following: 1. We establish the formal existence of localised radial patterns for a free-boundary problem.
2. We provide an approach to reduce the ferrofluid problem to an infinite system of ODEs.
3. Using methods from ODE theory, we find four classes of localised radial patterns; up-spots (called spot A), down-spots (called spot B), up-rings, and down-rings.
Notably, this analytically verifies the existence of up-spot (or spot A) solutions which have been observed experimentally [34], but also predicts the existence of down-spot (or spot B), up-ring, and down-ring solutions. These later solutions have not been observed in the Rosensweig instability experiment, to the author's knowledge, but have been proven to exist in the Swift-Hohenberg equation [28,30]. Future Work: This work serves as a framework for a rigorous proof of existence of each localised solution, which is in progress. Many parts of this work have been proven rigorously, however there are still some areas that prove to be more difficult. For example, the existence of infinitely many foliations and the regularity of our solutions are both non-trivial exercises to prove, in part due to the quasilinearity of the system, and the non-autonomous nature of the problem. There has recently been some progress made for quasilinear systems (see Chen et al. [10]) that may prove fruitful in this problem.
As mentioned in the previous subsection, we find an appropriate basis decomposition that allows us to reduce the problem to an ODE system. However, this is an approach that is tailored for our specific problem, whereas we would prefer to find a generalisation of these methods for a general quasilinear PDE of this form. The recent work of Beck et al. [3] for semilinear systems in radial domains may be useful for future progress in this area of study.
Some other areas of interest include: applying recent studies of radial invasion fronts (see Castillo-Pinto [9]) to the ferrofluid problem, where there is an opportunity for experimental validation; looking at the question of homoclinic snaking (see [27]) for localised radial patterns on a ferrofluid; and studying stability of our solutions, which is non-trivial due to the form of the Euler equations in our formulation of the problem.

The Radial Ferrohydrostatic Problem
In order to formulate the radial ferrohydrostatic problem, we will assume that all of the relevant variables are axisymmetric. Therefore, for any function f (r, θ, z), d dθ f (r, θ, z) = 0, where (r, θ, z) denote cylindrical polar coordinates. Then, the full 3-dimensional, finite-depth cylinder illustrated in Figure 1 can be reduced to a 2-dimensional infinite strip in (r, z) coordinates. Hence, we consider a domain consisting of two quiescent immiscible inviscid incompressible fluids separated by a free surface z = η (r), as shown in Figure 9. The upper fluid has density ρ + and unit relative permeability, while the lower fluid is a ferrofluid with density ρ − . Denoting H + , H − and B + , B − as the respective magnetising and induction fields of the fluids, we assume that these fields satisfy the relation where µ 0 is the magnetic permeability of free space and the ferrofluid is assumed to have a linear magneti- For a static fluid with no electric field, Maxwell's equations state that the magnetising and induction fields are irrotational and solenoidal, respectively. Defining magnetic potential functions φ − , φ + , we can write are the axisymmetric gradient and divergence operators for axisymmetric cylindrical polar coordinates (r, z) with related unit vectors ( r, z). Substituting in the linear relation (2.1), (2.2) reduces to Laplace's equation in both domains with the following magnetic continuity equations at z = η(r) Here n = ∇(z−η(r)) |∇(z−η(r))| is the unit normal vector to the interface, and it follows that at z = η(r), where subscripts denote partial derivatives. The ferrohydrostatic Euler equations for an incompressible fluid are given by , where g is the acceleration due to gravity, and p + , p − are the respective pressures of the upper and lower fluid. Integrating over their respective domains and equating on the interface, these equations are equivalent to , and b ± are the respective Bernoulli constants for each domain D ± . The ferrohydrostatic boundary condition at z = η(r) is the continuity of the normal stress is the mean curvature of the free surface. By substituting (2.6) into (2.5) to eliminate the relative pressure p 0 , we find That is, φ − , φ + , and η must satisfy 3), (2.4) and (2.7), and corresponds to a flat free surface and a uniform magnetic field with applied field strength h. In the standard Rosensweig instability experiment, h has a critical value at which domain-covering peaks begin to form; we call this value h c and define the bifurcation parameter of the problem to be ε := h − h c such that the Rosensweig bifurcation occurs at ε = 0. Also, choosing inhomogeneous Neumann boundary conditions where variations are taken with respect to φ − , φ + , and η, and subject to the constraint φ − = φ + at z = η(r).
Hamiltonian formulation: One of the first difficulties encountered in this problem is the r-dependent nature of the domains D ± , due to the boundary z = η(r). In order to apply analytical techniques to solve (2.9), we need the domains to be fixed for all r. This problem is often overcome in the water waves literature via a 'flattening' transformation of the variable z, as seen for the ferrofluid problem in Twombly & Thomas [40], to which maps the respective domains D + , D − to the infinite strips R × (0, D) and R × (−D, 0), and the free surface z = η(r) to z = 0. The 'flattened' variables ψ ± (r, z) = ψ ± (r, z) satisfy the equations with boundary conditions and and the tildes have been dropped for notational simplicity. Note that equations (2.13)-(2.17) can be obtained from the new variational principle δ L = 0, where and the variations are taken with respect to ψ − , ψ + , η satisfying the constraint (2.15) at z = 0. We consider L as an action functional where L is the integrand on the right-hand side of equation (2.18). We carry out a Legendre transformation and introduce conjugate variables α − , α + , γ defined by and thus the Hamiltonian function H (ψ − , ψ + , η, α − , α + , γ) is given by Hamilton's equations are given explicitly by 26) and are complemented by the constraint (2.15) and boundary conditions The boundary conditions (2.27) and (2.28) are found as a result of taking variations of (2.20) with respect to ψ − , ψ + and η, subject to (2.15). The condition (2.29) is equivalent to saying [ψ − r − ψ + r ] z=0 − (µ− 1)η r = 0 in the original, non-'flattened' coordinates. We note that our equations are invariant under the transformation ψ − → ψ − + c, ψ + → ψ + + c for any constant c. To eliminate the translational symmetry of the problem, and Then, we find that Note that A α is a conserved quantity, and so we set A α = 0, without loss of generality. Dropping the hats for notational simplicity, we observe that (2.21) and (2.22) from Hamilton's equations become Note that equations (2.30),(2.31) and (2.23) appear to be singular as r → 0, while equations (2.24)-(2.26) appear unbounded as r → ∞. In order to resolve the singular behaviour as r → ∞, we introduce the new variables This transformation again appears to be singular as r → 0, however we recall that the condition (2.8) implies that 1 r f r remains bounded as r → 0 for any arbitrary function f (r, z), and α ± , γ are Legendre conjugates associated with ψ ± r , η r , respectively. In particular, from (2.19), we see that which remain bounded as r → 0. Applying the transformation (2.32), our equations become where we have removed the tildes for notational simplicity, with constraints (2.15) and and boundary conditions (2.27) and Here, we have defined ν 0 := d dr A ψ , so In the form of Υ defined following (2.12), we consider h = h c + ε as a perturbation from the critical field strength for the Rosensweig instability. Thus, The ferrohydrostatic problem can then be written as where u = (ψ − , ψ + , η, α − , α + , γ) ⊺ and g(u, ε, r) is a nonlinear function containing (up to and including) second-order derivatives in z. The function g is given by the right-hand side of (2.33)-(2.38), and can be written as a non-autonomous linear part, L(r)u, and a collection of nonlinear terms F (u, ε, r). In addition, we have the linear constraints (2.15) and (2.39) and nonlinear boundary conditions which we write as B(u) = 0, where B is given by the left-hand sides of (2.27), (2.40), and (2.41). In order to apply centremanifold theory to our problem, and to keep consistent with the one dimensional problem in [19], we require the problem to have linear boundary conditions. We can overcome this problem via the transformation u = Gu, where with η and γ left unchanged, and τ and ν 1 are defined so that the constraint (2.39) is satisfied (see [19] for an almost identical transformation). Then, u := Gu satisfies the equations with boundary conditions B 1 u = 0, where B 1 u is the linearisation of B(u) and we have dropped the tildes for notational simplicity. The transformation G is a near-identity operator, leaving the linear part of g invariant; therefore the related linearised system is with linear constraints (2.15) and (2.39). We define the autonomous linear operator L ∞ := lim r→∞ L(r), which is equivalent to the linear operator of the one-dimensional ferrohydrostatic problem found in [19].

'Spectral' Decomposition
We wish to decompose u onto an r-independent basis. We recall that L(r) → L ∞ as r → ∞, and we will show that the eigenmodes of L ∞ , which are independent of r, are a good choice of basis. Here L ∞ is a linear operator acting on a subset of the space X , where Through explicit calculations, one can show that a non-zero complex number λ is an eigenvalue of the linear problem u r = L ∞ u if and only if λ = ς D , and ς satisfies the dispersion relation There are two eigenvalues with double algebraic multiplicity at λ = ±ik, and infinitely-many real eigenvalues λ = ±λn, with π 2D < λn < λn+1 for all n ∈ N. b) Purely imaginary eigenvalues plotted as functions of the parameters M and Υ0. A Hamiltonian Hopf bifurcation occurs at points on the curve C, and localised patterns emerge in the shaded region.
furthermore, even though ς = 0 satisfies (3.1), λ = 0 is not an eigenvalue of L ∞ when restricted to X . A purely imaginary number λ = ik is an eigenvalue of (3.

1) if and only if
This equation has either zero, one or two pairs ±k D of solutions, as seen in Figure 10b). A Hamiltonian Hopf bifurcation takes place for values for any k D ∈ (0, ∞); at this point, two pairs of simple, purely imaginary eigenvalues become complex by colliding on the imaginary axis at ±ik, and forming two Jordan chains of length 2. We therefore set The spectrum of L ∞ now consists of two eigenvalues ±ik with double algebraic multiplicity and a countably infinite family {±λ n } n∈N , see Figure 10a), where {λ n } n∈N is the set of all positive real solutions of ∆ 0 (λ n D) = 0. Furthermore, {λ n } n∈N is an increasing sequence, and λ 1 > π 2D (see Appendix 8.1). The limiting linear system u r = L ∞ u has a Hamiltonian structure, equipped with a symplectic two-form Ω, defined as where Ω(f , e) = Ω(e −n , e n ) = 1, Ω(e, f ) = Ω(e n , e −n ) = −1, and the 'symplectic products' of all other distinct combinations are zero. Here, 1 denotes the identity matrix and overbars denote the complex conjugate. The set {e, e, f , f } ∪ {e n } n∈Z\{0} can be proven to be a Riesz basis in X [17], such that the set is a complete basis that converges independently of its ordering (called unconditional convergence). This means that we can decompose our solution u in terms of a 'spectral decomposition' where a, b ∈ C, a ±n ∈ R for all n ∈ N, and the dynamics of u is tied directly to the dynamics of the amplitudes a, b, and a ±n . We introduce the projections with respective complements Q j := 1 − P j . By applying P 0 to the full system (2.42) and using the fact that {e, e, f , f } ∪ {e n } n∈Z\{0} is a Riesz basis, we find that a(r), b(r) satisfy the complex amplitude equations where F (ε, r, a, b, a 1 , a −1 , a 2 , . . . ) = F (u, ε, r) for u defined in (3.4). Similarly, by applying P n for each n ∈ N respectively, we find that d d r a n = λ n a n − 1 2r (3.8)

Core Solutions
We first wish to construct the 'core' manifold, i.e. the set of all small-amplitude solutions to (2.42) that remain bounded as r → 0. We begin by investigating the linear problem of (3.5) and (3.6), i.e. with F ≡ 0, thus By performing a change of coordinates into the real and imaginary parts of both a and b respectively, a simple calculation provides the general linear solution where J ν , Y ν are ν-th order Bessel functions of the first and second kind, respectively. Similarly, for the linear problem related to (3.7) and (3.8) d d r a n = λ n a n − 1 2r (a n − a −n ) , d d r a −n = −λ n a −n + 1 2r (a n − a −n ) , we find, by defining b n := an+a−n 2 and b −n := an−a−n 2 and solving the resulting equations, the general linear solution a n = a n := 2 j=1 c j,n W j,n (r), (4.2) for a n := (a n , a −n ) ∈ R 2 , where W 1,n (r) = λ n 2 I 0 (λ n r) + I 1 (λ n r), I 0 (λ n r) − I 1 (λ n r) , for n ∈ N. Here, I ν and K ν are ν-th order modified Bessel functions of the first and second kind, respectively. We note that, in both cases, the respective linear solutions V j and W j,n split into two types; solutions (V 1 , V 2 and W 1,n ) that are bounded as r → 0, and solutions (V 3 , V 4 and W 2,n ) that are unbounded as r → 0; see Table 1. By denoting ·, · as the complex dot product {x n y n + x n y n } , then the respective eigenmodes V i (r) and W i,n (r) satisfy where δ ij is the Kronecker delta, and V * j (r), W * i,n (r) are the respective adjoint vectors of V j (r), W i,n (r), defined by We apply the method of variation of constants [43, §2.II] in conjunction with the general linear solution (4.1), to find a general solution to the full nonlinear system (3.5)-(3.6) for the complex amplitudes a, b. This solution has the form where F := −Ω F , f , Ω (F , e) and e(z), f (z) are the eigenmodes defined in (3.3) for the imaginary eigenvalues ik and a large fixed radial value r 0 > 0. Similarly, for the pairs of real-valued amplitudes a n = (a n , a −n ), we take the linear solution (4.2) and apply the method of variation of constants to derive a general solution for the full nonlinear system (3.7)-(3.8), which has the form a n (r) = 2 j=1 W * j,n , a n W j,n (r) + W 1,n (r) r r0 W * 1,n , F n ds + W 2,n (r) r 0 W * 2,n , F n ds, (4.4) , Ω (F , e n )) and e n (z), e −n (z) are the eigenmodes defined in (3.3) for the respective real eigenvalues λ n , −λ n . Using the r → 0 asymptotic forms for the Bessel functions detailed in Table 1, we note that equipped with the standard ℓ 1 -norm Table 1: Expansions of the Bessel functions Jν (r), Yν(r), Iν(r), Kν(r) for r → 0 and r → ∞ (see [1], (9.1.10-11), (9.6.10-11), §9.2 and §9.7).

Far-field Solutions
In this section, we now turn our attention to the far-field region, where r is larger than some fixed r 1 , and construct the 'far-field' manifold; the set of all small-amplitude solutions that decay to zero exponentially Figure 11: The centre-stable manifold W cs + is parametrised by the foliation F s (p, σ) over base points p in the centre manifold W c + . Stable fibres are displayed in respective Poincaré sections σ = σ1, σ2, with base points p1, p2 ∈ W c + . Solutions related to initial data q1 ∈ F s (p1, σ1) approach the solution related to the initial condition p1 on the centre manifold exponentially as r → ∞.
fast as r → ∞. We augment (2.42) by introducing the variable σ := 1 r so that the extended system becomes By introducing σ and treating it as a unconstrained function of r, (5.1) is now considered to be an autonomous system. This is because the 1/r terms in L(r) in (2.42) are now of 'quadratic' order, and hence they are absorbed into the nonlinearity, so the remaining linear operator L ∞ is independent of r. By applying the same 'spectral decomposition' as detailed in Section 3, we arrive at the extended versions of (3.5)-(3.8), namely ) − Ω(F , f ), d dr a n = λ n a n − σ 2 (a n − a −n ) − Ω(F , e −n ), for n ∈ N, with the extra condition (5.1). Linearising about the trivial state a = b = a n = a −n = σ = 0 for all n ∈ N, the system has the form d dr a = ika + b, d dr a n = λ n a n , d dr We note that the equation for σ decouples from the rest of the system and so we can solve this to find σ = σ * , where σ * is an arbitrary constant. Therefore, for sufficiently small solutions, we can look at the spatial dynamics of the linear problem (5.2) on the Poincaré section σ = σ * , where σ * is constant, see Figure  11. For a fixed σ, the rest of the linear problem (5.2) has four centre directions given by a ±n = 0 for all n ∈ N, a countably infinite set of stable directions given by a = b = a n = 0 for all n ∈ N, and a countably infinite set of unstable directions given by a = b = a −n = 0 for all n ∈ N; see Figure 12a). Standard parametrisations would define the stable manifold by displacements in each of the stable directions, and similarly for the centre and unstable manifolds with the centre and unstable directions, respectively.
We also note that the four-dimensional centre manifold can be decomposed further: two directions define solutions that remain bounded for all r, while the other two directions define solutions that grow algebraically as r → ∞. Then, the expected four-dimensional centre manifold W c + , containing all solutions with at most algebraic growth for all r ∈ [r ∞ , ∞), is also expected to have a two-dimensional submanifold W c + ⊂ W c + , which contains all solutions on the centre manifold that are bounded as r → ∞.
The following theory is well-defined in infinite dimensions; however, in order to discuss the dimensions of each manifold, we will refer to a countably infinite set of dimensions as being n-dimensional. Conceptually, we think of this as an order n-truncation as n becomes very large. Formally, we would expect to find a (4 + n)-dimensional centre-stable manifold W cs + , containing all solutions that have at most algebraic growth as r → ∞; an n-dimensional stable manifold W s + , containing all solutions that exponentially decay as r → ∞; a four-dimensional centre manifold W c + , containing all solutions with at most algebraic growth for all r ∈ [r ∞ , ∞); and a two-dimensional submanifold W c + , that contains solutions on W c + that remain bounded as r → ∞; see Figure 12b) for a sketch of this idea.
We want to parametrise the far-field stable manifold W s + , containing all small-amplitude solutions that exponentially decay to zero as r → ∞, which we would expect to be n-dimensional. However, there are a few problems: • The core manifold W cu − parametrised in §4 is (2 + n)-dimensional, so there are not currently enough parameters to match with.
• In order to perform standard analytical techniques, we wish to apply normal-form transformations to make our system more tractable. However, these are hard to formulate unless working on a centre manifold [38].
Hence, we need to find a way to parametrise the far-field manifold in (2 + n) dimensions, such that any required analysis can be applied on the centre manifold. In [31], McQuighan & Sandstede faced a similar problem for the case of a finite-dimensional ODE system, which they overcame by utilising the theory of stable foliations; see [2,11].
The motivating theory behind this can be seen as follows; we already established that we expect to find a (4 + n)-dimensional centre-stable manifold. By definition, the centre and stable manifolds should both be contained within the centre-stable manifold, such that we could parametrise the stable manifold as an n-dimensional manifold contained in the larger centre-stable manifold. This parametrisation would isolate all of the solutions that decay exponentially as r → ∞, but there would only be n parameters, causing the matching problem.
We could parametrise the stable manifold with respect to the coordinates of the centre-stable manifold, that is, displacements in all of the centre and stable directions. By performing analysis on each parameter to ensure they decay to zero as r → ∞, we could reduce this manifold to (2 + n)-dimensions, and so the matching problem would be solved. However, formulating the required analysis remains a difficult task.
We will instead choose to parametrise the stable manifold in relation to the four-dimensional centre manifold. By taking a point in the centre manifold, parametrised by the four centre manifold parameters, and then taking displacements from that point in each of the stable directions, we can construct a (4 + n)-dimensional parametrisation that covers the entire centre-stable manifold close to the trivial state u ≡ 0. As before, in order to isolate the far-field stable manifold and reduce the parametrisation to (2 + n) dimensions, we would need to ensure that every parameter decays to zero as r → ∞. We know that the stable displacements tend to zero, and so the final analysis can be conducted on the remaining centre manifold parameters. This allows us to perform normal-form transformations on the four-dimensional centre manifold parameters in order to find the set of solutions that lie on the centre manifold and decay exponentially as r → ∞, which we call the stable part of the centre manifold and denote as W s c ⊂ W c + . Since any solution in W s c also remains bounded as r → ∞, we note that W s c ⊂ W c + , and so W s c can be parametrised as a two-dimensional manifold. In particular, as illustrated in Figure 11, the entire centre-stable manifold W cs + | σ=σ * can be written as a foliation, i.e. a union of stable fibres, p {F s (p, σ * )} over base points p in the full centre manifold W c + | σ=σ * ⊂ W cs + | σ=σ * , for σ = σ * fixed. Then, a trajectory beginning at (q, σ * ) will converge to zero as r → ∞ if and only if its related base point (p, σ * ) on the centre manifold does. The stable manifold, defined as the set of all solutions that exponentially tend to zero as r → ∞, is hence given by the union of all stable fibres associated with solutions on the centre manifold that exponentially decay as r → ∞, We begin by finding a parametrisation of all stable fibres F s (p, σ * ) with respect to their base points p on the full centre manifold W c + | σ=σ * . Following this, we analyse the dynamics of the system restricted to W c + | σ=σ * and determine initial conditions required for exponentially decaying solutions. This parametrises the twodimensional manifold W s c ⊂ W c + and reduces the parametrisation of p {F s (p, σ * )} to (2 + n)-dimensions. Finally, we substitute our initial conditions into the full infinite-dimensional parametrisation of the stable foliations, now with exponentially decaying base points, which exactly defines the stable far-field manifold W s + , and we match these solutions to the core solutions derived in Section 4.

Parametrisation of the far-field manifolds
By the same methods as outlined in Section 4, we can construct a local 'centre-stable' manifold which contains all solutions to the full far-field system which have mild exponential growth 1 as r → ∞. with respective adjoint solutions We note that { V 3 , V 4 } are bounded as r → ∞, { V 1 , V 2 } grow algebraically as r → ∞, W 2,n (r) decays exponentially as r → ∞, and W 1,n (r) grows exponentially as r → ∞, for any n ∈ N. Applying the method of variation of parameters and solving the respective equations similar to (4.3), (4.4) for small enough (d, c 2 ), the parametrisation of the centre-stable manifold takes a form similar to (4.6), where O σ (·) denotes the standard Landau symbol where the bounding constants may depend on the value of σ. We introduce the centre manifold W c + , which contains all solutions that are bounded for all r, and evaluate it on the Poincaré section σ = σ * , We note that any solution that remains bounded for all r, i.e. is on the centre manifold W c + , must remain bounded as r → ∞, by definition. Hence, the set of all solutions bounded for all r is contained within the set of all solutions bounded as r → ∞, and so W c + ⊂ W cs + for any σ = σ * . Thus, the centre and centre-stable parametrisations should coincide when c 2,n is slaved to the centre coordinates d for all n ∈ N; this is a key assumption which we do not prove in this work. That is, when c 2,n = O σ |d| |ε| + |d| , we find We now wish to parametrise W cs + | σ=σ * as the union of strong stable fibres p F s (p, σ * ) over base points p ∈ W c + | σ=σ * . We refer to [11] for background on stable foliations, [31] for a finite dimensional example, and Figure 13 for an illustration. We introduce a base point p ∈ W c + | σ=σ * with the decomposition such that, since p ∈ W c + | σ=σ * , we can use the parametrisation (5.3) to write for all n ∈ N. Then, we can write the stable fibre of trajectories related to the point p as  Figure 13: Parametrisation of a strong stable fibre F s (p, σ * ) with base point p ∈ W c + |σ=σ * for some σ * ≤ 1/r∞ fixed. Any solution with initial datum q ∈ F s (p, σ * ) will tend to solutions with initial datum p ∈ W c + |σ=σ * as r → ∞.
We have now parametrised the centre-stable manifold W cs + | σ=σ * as the union of the stable foliations p F s (p, σ * ) over every point p ∈ W c + | σ=σ * as defined in (5.4). In order to identify the far-field manifold W s + | σ=σ * containing all solutions that decay exponentially as r → ∞, we must isolate the values of the centre coordinates (A 0 , B 0 ) such that the p exponentially decays to zero as r → ∞. Hence, we reduce the problem to the centre manifold.

Reduction to the centre manifold
In order to identify solutions that are exponentially decaying as r → ∞, we need to analyse the dynamics of the centre manifold W c + . To do this, we first derive a normal-form for the extended system restricted to the centre manifold. We need that the hyperbolic part of the system is slaved to the centre modes; a key assumption not proven here. We recall the centre manifold parametrisation introduced in (5.4)-(5.5), and define a centre manifold reduction Then, the reduced vector field on W c + projected onto P 0 is d dr where P 0 is defined in Section 3, and L(σ) is the projection of linear operator L 1 σ restricted to the centre manifold W c + , where σ = 1/r. Then, A 0 , B 0 , σ satisfy a complex ODE system of the form d dr In order to transform (5.8) into a normal-form that we can analyse, we utilise the following result, Lemma 5.1 [28] There exists a change of coordinates where we define the remainder terms R A , R B below. This coordinate change is polynomial in (A, B, σ) and smooth in ε, and T (σ) = O(σ) is linear and upper triangular for each σ, while φ(r) satisfies The remainder terms satisfy This result, with different values of c 0 , c 3 , is proven in [28,Lemma 2], and so we will omit the details. This change of coordinates is called a normal-form transformation, where (5.10) is called the normal-form of the system. As seen in [28] and [38], (5.10) is equivalent to the normal-form of the radial Swift-Hohenberg equation.

Remark 5.2
The radial problem is an O(r −1 ) perturbation from the one-dimensional problem and so, with the added non-autonomous terms, we expect c 0 = c 0 + O(r −1 ) and c 3 = c 3 + O(r −1 ), where c 0 , c 3 are the respective coefficients for the one-dimensional problem. Then, extending the system with respect to σ = r −1 , these O(r −1 ) perturbations of c 0 , c 3 become O(σ) terms and so are lifted to the coefficients of the higher order terms σεA and σ|A| 2 A, respectively. Therefore, the values of c 0 , c 3 are identical to those found in [19, §5] for the one-dimensional ferrohydrostatic problem subject to a linear magnetisation law. That is, for fixed k D ∈ (0, ∞) and M 0 ∈ (0, 1), where k D := kD and M 0 := µ−1 µ+1 , we find and It is straightforward to show that c 0 is positive for all k D and M 0 , one can find a restriction on k D for a given M 0 such that c 3 < 0 ; this is required for the existence of homoclinic solutions in [19] and will be discussed in greater detail in Section 6.3.
It is convenient to extend the transformation in Lemma 5.1 and apply it to our matching coordinates (a, b) 15) such that, in the far-field Using these new coordinates, we can write our core parametrisation (4.8) as where we have used φ(r 0 ) = kr 0 + O(r −2 0 ) + O r0 (|ε| + |d| 1 ), from (5.11). We have also exploited the fact that T (σ) is upper triangular and the coefficient in front of d 2 2 for b(r 0 ) scales with √ r 0 so that it is only affected by the transformation (5.15) at higher order. In the next section, we will be rescaling variables by √ ε, and so it is convenient to define E := √ ε in order to minimise notational complexity.

Geometric Blow-up
So far, we have constructed the core manifold in Section 4 on bounded intervals r ∈ [0, r 0 ], and the far-field stable manifold for r ≥ r ∞ . We choose r 0 large enough such that r ∞ < r 0 , ensuring that W cu − and W s + both exist at the matching point r = r 0 . Following the work of [28,30], we expect decaying solutions on the far-field manifold to have the form A(r) = EA 2 (s), where s = Er, for some function A 2 (s). Then, after this transformation, equations (5.10) have nontrivial solutions that decay exponentially as s → ∞; these are valid for s ≥ δ 0 , i.e. r ≥ δ 0 /E, where δ 0 > 0 is fixed. However, we lose control over bounds on the solution as E → 0. This 'matching gap' is resolved by a second coordinate transformation to maintain control of solutions for r 0 ≤ r ≤ δ 0 /E, as shown in Figure 14. We employ the coordinate transformations used by McCalla & Sandstede [30] for the two-dimensional Swift-Hohenberg equation, as our normal-form (5.10) is equivalent up to leading order.

Rescaling chart
Following the procedure implemented in [30], we introduce the variable z ∼ d dr A /A. This is not to be confused with the vertical coordinate in Section 2. From (5.10), we define and introduce rescaling coordinates In these new coordinates (5.10) becomes where it can be found from (5.12) that Hence, we obtain the system Note that the last equation implies that the variable ε 2 can be arbitrarily fixed; setting ε 2 = 0, we arrive at the system d ds which, for a small fixed ε 2 , has two families of equilibria; Q ± (ε 2 ) = 0, ± √ c 0 , 0, O(|ε 2 | 2 ) with eigenvalues ± √ c 0 , ∓2 √ c 0 , 0, 0 . For both families, the zero eigenvalues have eigenmodes (0, 0, 1, 0) and (0, 0, 0, 1). For the subspace A 2 , z 2 ∈ R, (6.4) can be reduced to the non-autonomous real Ginzburg-Landau equation, By analysing (6.5) linearised about the trivial state for large values of s, one can see that solutions to this equation decay or grow exponentially with rate ± √ c 0 as s → ∞. We see from (6.4) that z 2 = A ′ 2 /A 2 and so, if A 2 ∼ exp(± √ c 0 s), then z 2 = ± √ c 0 . Therefore, to find solutions that are exponentially decaying, we are interested in the centre-stable manifold W cs (Q − ) of the family of equililbria Q − (ε 2 ). Then, we will capture all solutions of (5.10) of size E that decay as r → ∞, for E > 0.

Transition chart
We rewrite (5.10) in terms of coordinates (A 1 , z 1 ) that are obtained by rescaling (A, z) by σ = 1/r, where z is defined in (6.1). In particular, we define Then, we can write (5.10) as where R A , R B are evaluated at (A, B, σ, E) = A 1 , σ 2 1 A 1 z 1 + 1 2 , σ 1 , σ 1 ε 1 . From (5.12), one can find that 1 as shown in [30, §2.4]. Thus, the flow in the transition chart is defined by where (A 1 , z 1 , σ 1 , ε 1 ) ∈ C × C × R × R. For c 3 < 0, (6.7) has precisely two equilibria, P ± := 0, ± 1 2 , 0, 0 . We require a good understanding of the dynamics close to each equilibria, and so we introduce two sets of coordinates that move each respective equilibria to the origin. First, near P + , we use the new variables The linearisation of (6.8) about the origin has the eigenvalues where the linearisation about the origin has the eigenvalues { 1 2 , 1, −1, 1}. Finally, we note that the coordinates in the transition and rescaling charts are related by the transformation and we can therefore transform from one chart to the other in the transverse section ε 1 = σ 2 = 1.

Connecting orbits
We now investigate the dynamics of the system in the transition and rescaling charts. Our goal is to show the existence of heteroclinic orbits q ± (s) that connect the respective equilibria P ± in the transition chart to the equilibrium Q − (ε 2 ) in the rescaling chart in the limit E = 0; see Figure 15 for a schematic diagram. First we note that the subspace σ 1 = ε 2 = 0 is invariant under the flow of the equations in the transition and rescaling charts. Setting σ 1 = 0 in (6.7), we have with (A 1 , z 1 , ε 1 ) ∈ C × C × R. Setting ε 2 = 0 in (6.4), we obtain d ds 12) with (A 2 , z 2 , σ 2 ) ∈ C × C × R. We note that the transformation (6.10) between transition and rescaling charts maps σ 1 = 0 into ε 2 = 0. Also, the real subspaces (A 1 , z 1 , ε 1 ) ∈ R 3 and (A 2 , z 2 , σ 2 ) ∈ R 3 are invariant for (6.11) and (6.12), respectively, and so we will initially restrict our attention to these two subspaces. Figure 15: A sketch of the potential connecting orbits between the core and the equilibrium Q−(ε2) in the transition chart. There exist two connecting orbits: q±(s), connecting P± and Q−, defined in (6.16) and (6.14), respectively. The connections between the core and P+ and P− are u2(r), which corresponds to rings, and u1(r), which corresponds to spots, respectively.

(6.13)
In addition, the linearisation of (6.5) about q(s) does not have a nontrivial solution that is bounded uniformly on R + . If c 3 > 0, then the only bounded solution of (6.5) on R + is A 2 (s) ≡ 0.
This Lemma is proven in [30,Lemma 2.3], and relies on a hypothesis [30, Hypothesis 1] that the general non-autonomous real Ginzburg-Landau equation, where c 0 = 1, and c 3 = −1, has a nontrivial bounded solution q(s) on R + and the linearisation around q(s) does not have a nontrivial bounded solution that is bounded. This hypothesis was later proven with rigorous numerical methods by van den Berg et al. [41,Theorem 1.1]. Next we write the solution q(s) as in the coordinates of the transition and rescaling charts and conclude that the functions , e ρ ∈ R 3 , (6.10) 14) satisfy (6.11) and (6.12), respectively. In [30, Lemma 2.4], McCalla & Sandstede prove that (6.14) is a connecting orbit that forms a transverse intersection of the unstable manifold W u (P + ) of the equilibrium P + = (0, 1/2, 0) of (6.11) and the centre-stable manifold W cs (Q − ) of the equilibrium Q − (0) = (0, − √ c 0 , 0) of (6.12). We note that the real-valued heteroclinic orbit given in (6.14) generates a one-parameter family q + (s) : (A 1 , z 1 , ε 1 ) (ρ) = e iY A + 1 , z + 1 , ε + 1 , (6.10) of heteroclinic orbits of (6.11) and (6.12) between equilibria P + , Q − (0) that are parametrised by Y ∈ R, where orbits related to the choice of Y = 0 and Y = π both lie in the invariant subspace R 3 .
The existence of the connecting orbit (6.15) is reliant on the sign of c 3 being negative, which we will now discuss in more detail. As we alluded to in Remark 5.2, the value of c 3 in the radial problem is equivalent to leading order to that found for the one-dimensional problem in [19], and can be written down explicitly in terms of k D and M 0 , as defined in (5.14). We see from (5.14) that one can find a restriction on k D for a given M 0 such that c 3 < 0, as seen in Figure 5a). For large depth fluids one has a critical magnetic permeability M c ≈ 0.56, or equivalently µ c ≈ 3.54, such that c 3 < 0 for all M 0 > M c ; for experimentally relevant ferrofluids (where 2 < µ < 6) one requires k D > k * D ≈ 1.8. For sufficiently small depths, Figure 5a) appears to suggest that q + (s) exists for all values of µ. However, our free-surface model becomes unphysical for these shallow depths, since one should then also include a magnetic potential below the ferrofluid. Therefore, we will ignore the small-k D shaded region seen in Figure 5a).

Connecting P − and Q −
We start by analyzing (6.11): rewriting the system in (A − , z − , σ − , ε − ) coordinates, where σ − = 0, we find and we obtain an explicit solution given by This solution (6.16) satisfies z 2 (s) → − √ c 0 as s → ∞ and z − (ρ) → 0 as ρ → −∞ (which is equivalent to s → 0), and therefore lies in the intersection of W u (P − ) and W cs (Q − ). We note the real-valued heteroclinic orbit given in (6.16) generates a one-parameter family of heteroclinic orbits of (6.11) and (6.12) between the equilibria P − , Q − (0) that are parametrised by Y ∈ R, where orbits related to the choice of Y = 0 and Y = π both lie in the invariant subspace R 3 .

Spot A
The construction of the up-spot, termed spot A in [29,30], proceeds as follows; we follow the centre-stable manifold W cs (Q − ) near the heteroclinic orbits q − (s) defined in (6.17) for Y = 0, π backward in 'spatial time' ρ past the equilibria P − , as shown in Figure 16. This will define the initial conditions required for solutions to decay to zero as r → ∞, parametrising the centre coordinates of the far-field manifold at r = r 0 . These centre coordinates are then substituted into the foliation parametrisation defined in (5.7) and the final matching between the core and far-field manifolds will take place at r = r 0 .
We begin by tracking the solution in the invariant subspace ε 2 = 0, i.e. at the bifurcation point. We know that the solution q − (s) lies on the intersection W u (P − ) ∩ W cs (Q − [0]). Therefore, we will first parametrise the spot A solution as a small perturbation away from the solution q − (s) at a fixed small value s = δ 0 > 0.

Matching Core and Far-Field Solutions: Spot A
We are now fully equipped to find a nontrivial solution contained in the intersection of the core manifold W cu − and the far-field manifold W s + . We begin by substituting the initial conditions (7.9) into the far-field parametrisation (5.7) following the normal-form transformation (5.15) to obtain a(r 0 ) = A(r 0 ) + O r0 | d| 1 |ε| + |A| + |B| + | d| 1 , The |ε| 1−κ 2 terms are due to the restriction |a| ∈ (ε − κ 2 , ε κ 2 ) for a fixed constant 0 < κ ≪ 1. We also recall the transformed core parametrisation (5.16) where |d| 1 = |d 1 | + |d 2 | + |c 1 | 1 , and |d 2 | 1 = |d 1 | + |c 1 | 1 . Setting these parametrisations equivalent to each other for each respective coordinate ( a, b, a 1 , a −1 , a 2 , a −2 , . . . )(r 0 ) is the same as finding the zeros of the functional and we have defined We introduce the scaling (d 1 , d 2 ) = (ε d 1 , ε 1 2 d 2 ) so that terms in (7.10) and (7.11) scale with the same order in ε. Initially, setting ε = 0, we investigate (7.12) and (7.13) Defining the functional F : Furthermore, the Fréchet derivative DF (0, 0) is invertible, and so, using the implicit function theorem we can solve (7.12) and (7.13) for all values of n ∈ N uniquely for sufficiently small 0 < ε ≪ 1. Matching orders of ε in (7.12) and (7.13), we find that Returning to (7.10) and (7.11), we have 14) Initially setting ε = 0, we obtain the system We formally set ∆ A = 0 and separate (7.15) into real and imaginary parts: this is equivalent to finding zeros of the functional It is apparent that the vector Since c 0 > 0, the Jacobian is invertible and we can, using the implicit function theorem, solve (7.15) uniquely for all sufficiently small ∆, that is, for r 0 large enough and δ 0 small enough, and subsequently (7.14) for all 0 < ε ≪ 1. Reversing the scaling for d, we find that Hence, we have found a spot A solution. We recall that our solution u(r, z) takes the form Here, we have denoted the the vertical coordinate used in Section 2 as z in order to reduce confusion with the geometric blow-up coordinate z. Substituting (7.16) and (7.17) into this form, we can write the spot A solution u A as u A (r, z) = ε for all r ∈ (0, r 0 ). In particular, using the explicit forms of e, f (8.3) defined in Appendix 8.2, the height of the free surface η A (r) has the form for all r ∈ (0, r 0 ). Similarly, for (A − , z − )(ρ) in (7.8), where ρ = ln r − ln(δ 0 ε − 1 2 ), we can invert the coordinate transformations (6.6), (6.1), and (5.9) to find the free surface profile in the transition chart, for r ∈ [r 0 , δ 0 ε − 1 2 ]. Finally, for (A 2 , z 2 )(s) in (7.1), we can invert transformations (6.2), (6.1), and (5.9) to find the free surface profile in the far-field, . This completes the result for spot A.

Spot B
The construction of the secondary spot, termed spot B in [30], proceeds as follows; we follow the centre-stable manifold W cs (Q − ) near the heteroclinic orbits defined in (6.15) for Y = 0, π backward in 'spatial time' ρ past the equilibria P + and P − , as shown in Figure 18. This will define the initial conditions required for solutions to decay to zero as r → ∞, parametrising the centre coordinates of the far-field manifold at r = r 0 . These centre coordinates are then substituted into the foliation parametrisation defined in (5.7) and the final matching between the core and far-field manifolds will take place at r = r 0 .
The construction of spot B in this paper is almost identical to the two-dimensional Swift-Hohenberg equation, which is completed explicitly in [30]. Therefore, we will provide an idea of how to construct spot B, and then state the results which have been proven previously in [30]. The journey in backward ρ taken by the spot B solution through the transition chart can be divided as follows: Section 7.2.1, where the solution moves close to P + ; Section 7.2.2, where the solution travels between P + and P − ; and Section 7.2.3, where the solution moves close to P − . The decomposition of this journey is illustrated in Figure 19.

Dynamics near P +
We begin by tracking the centre-stable manifold W cs (Q − ) for E > 0 backwards in ρ as it passes near to the equilibrium P + . We consider the transition variables (A + , z + , σ + , ε + ) such that P + corresponds to the origin. It is convenient to perform a change of coordinates in order to remove some of the higher order terms from the equations (6.8). Therefore, we apply the following Lemma, Lemma 7.1 [30] There is a smooth coordinate change of the form that transforms (6.8) near the origin into For a proof of this result, see [30,Lemma 3.1]. To find an expression for the manifold W cs (Q − ) near P + , we recall that the solution (6.14) forms a transverse intersection of W cs (Q − ) and W u (P + ). Then, Lemma 7.2 [30] For each sufficiently small δ 0 > 0, there are constants a 0 , ε 0 > 0 such that the following is true. Define Σ 0 := {ε + = δ 0 }; then where η 0 (δ 0 ) = q 0 δ 3/2 0 (1 + O(δ 0 )) is smooth, q 0 > 0 is the constant defined in (6.13), and Y ∈ R is arbitrary.
This result is proven in [30,Lemma 3.4]. Converting back into (A, B) coordinates, by

Matching Core and Far-Field Solutions: Spot B
We are now fully equipped to find a nontrivial solution contained in the intersection of the core manifold W cu − and the far-field manifold W s + . We begin by substituting the initial conditions (7.26) into the far-field parametrisation (5.7) following the normal-form transformation (5.15) a(r 0 ) = A(r 0 ) + O r0 | d| 1 |ε| + |A| + |B| + | d| 1 , We also recall the transformed core parametrisation (5.16) where |d| 1 = |d 1 | + |d 2 | + |c 1 | 1 , and |d 2 | 1 = |d 1 | + |c 1 | 1 . Setting these parametrisations equivalent to each other for each respective coordinate ( a, b, a 1 , a −1 , a 2 , a −2 , . . . )(r 0 ) is the same as finding the zeros of the functional where and we have defined We introduce the scaling (d 1 , d 2 ) = (ε 3 4 d 1 , ε 3 8 d 2 ) so that terms in (7.27) and (7.28) scale with the same order in ε. Initially setting ε = 0, we investigate (7.29) and (7.30) and, defining the functional F : ( Furthermore, the Fréchet derivative DF (0, 0) is invertible, and so we can solve (7.29) and (7.30) for all values of n ∈ N uniquely for sufficiently small 0 < ε ≪ 1. Matching orders of ε in (7.29) and (7.30), we find that Returning to (7.27) and (7.28), we have Initially setting ε = 0, we obtain the system We formally set ∆ B = 0 and separate (7.32) into real and imaginary parts: this is equivalent to finding zeros of the functional The vector Since q 0 > 0, the Jacobian is invertible, and we can therefore solve (7.32) uniquely for all sufficiently small ∆, that is, for r 0 large enough and δ 0 small enough, and subsequently (7.31) for all 0 < ε ≪ 1. Reversing the scaling for d, we find that Hence, we have found a spot B solution. We recall that our solution u(r, z) takes the form Here, we have denoted the the vertical coordinate used in Section 2 as z in order to reduce confusion with the geometric blow-up coordinate z. Substituting for all r ∈ (0, r 0 ), where u B (r, z) decays to zero exponentially as r → ∞. In particular, the height of the free surface η B (r) has the form for all r ∈ (0, r 0 ). Following Lemma 7.4, we can define write the linear flow near P − as where ρ = ln r, for r ∈ r 0 , aε − 3 8 δ2 (1−δ1) (1−δ2) . For the flow between P − and P + near the z 1 -axis we can explicitly solve (7.22) to find . Near the equilibrium P + , we can define solutions by the linear flow, such that we find, to leading order, where ρ = ln r, for r ∈ aε − 3 8 δ1 , δ 0 ε − 1 2 . Finally, for (A 2 , z 2 )(s) in (6.15), we can apply the transformation and so, from (6.13), we see that Inverting the transformations (6.2), (6.1), and (5.9), we find the free surface profile in the far-field in each transition region, which completes the result for spot B.

Rings
The construction of up-ring and down-ring solutions proceeds as follows; we follow the centre-stable manifold W cs (Q − ) near the heteroclinic orbits defined in (6.15) for Y = 0, π backward in 'spatial time' ρ past the equilibria P + , as shown in Figure 20. This will define the initial conditions required for solutions to decay to zero as r → ∞, parametrising the centre coordinates of the far-field manifold at r = r 0 . These centre coordinates are then substituted into the foliation parametrisation defined in (5.7) and the final matching between the core and far-field manifolds will take place at r = r 0 .
Proof. As ∆ 2 (ς) is an even function, we restrict ς = x ∈ R j for j ∈ N; any results for x j ∈ R j are mirrored for x −j = −x j ∈ R −j . Since, for each j ∈ N, ∆ 2 (x) → ±∞ as x → (2j±1)π 2 , we apply a limiting version of the Intermediate Value Theorem: one can find constants a j , b j ∈ R j such that a j < b j , ∆ 2 (a j ) < 0, and ∆ 2 (b j ) > 0. Then, by the IVT, there exists a point x j ∈ [a j , b j ] such that ∆ 2 (x j ) = 0. Therefore, we have proven that ∆ 2 (ς) has at least one root ς = x j ∈ R j , for each j ∈ N. Using the even properties of ∆ 2 (ς), there is an x j ∈ R j for j ∈ Z\{0}, where x −j = −x j . For x ∈ R 0 , one can show that ∆ 2 (x) is convex, and that ∆ 2 (0) is a minimum point. Hence, ∆ 2 (x) ≥ ∆ 2 (0) = Υ 0 > 0. Consequently, ∆ 2 (x) > 0 for all x ∈ R 0 . Lemma 8.3 There exists some fixed y * such that for all y > y * , ∆ 2 (x + iy) = 0.
We define the complex regions K j := R j × i [−y * , y * ], where R j is defined in Lemma 8.2, and y * is some large real number as seen in Lemma 8.3; see Figure 21. Thanks to Lemmas 8.1 and 8.3, we can decompose the kernel of ∆ 2 (ς) as the union j∈Z K j . From now on, we will investigate zeros of ∆ 2 (ς) with ς restricted to K j for some arbitrary j ∈ Z. Lemma 8.4 Define ∆ 3 (ς) : C → C, ς → −ς [M tan(ς) − ς]. Then, ∆ 3 (ς) has four zeros in K 0 and one zero in K j , for all j ∈ Z\{0}.
Proof. We define ς = x + iy and begin with the case when y = 0. Since M tan(x) is a monotonicallyincreasing surjective function for x ∈ R j , ∆ 3 (x) has a repeated zero at x = 0 ∈ K 0 , and a unique non-trivial zero x = x j ∈ K j , for all j ∈ Z\{0}.
Finally, we can explicitly find the kernel of ∆ 2 (ς). We will make use of the symmetric version of Rouché's theorem (See Glicksberg [16]), and so we state it here: Theorem 1 [16] Let K ⊂ C be a bounded region with continuous boundary ∂K. Two holomorphic functions f, g : C → C have the same number of zeros in K, if the strict inequality |f (ς) + g(ς)| < |f (ς)| + |g(ς)| , holds for all values of ς on the boundary ∂K, Proposition 8.5 ∆ 1 (ς) has precisely two purely imaginary zeros ς = ±ikD, each with multiplicity 2, and an infinite set of purely real zeros ς = ± λ j D, j ∈ N, with multiplicity 1.
However, for (M, Υ 0 ) = (M H , Υ H ), we have already four zeros in K 0 , namely ς = ±ikD with double algebraic multiplicity. We define the unique root K j for each j ∈ Z\{0}, namely λ j D := x j as found in Lemma 8.2. By (8.2), the result is proven.

Basis of z -dependent eigenmodes
In Section 3, we introduce a 'spectral' decomposition to reduce the problem to a system of amplitude equations. To do this, we find eigenmodes of the linear operator L ∞ . Explicitly, for the imaginary eigenvalues , if λ n = jπ D , for some j ∈ N. Here, c 1,n and c 2,n are defined such that Ω(e −i , e j ) = δ i,j for any i, j ∈ N.

Core Problem: Quadratic Term
In order to determine the existence of spot A solutions, we need to find the coefficient of the d 2 2 V 3 (r) term in the (4.3) expansion for b(r). To do this, we Taylor expand d 3 and isolate the d 2 2 term. We recall from (4.3) that We set ε = 0 and recall that V * 3 (r) = kπ 2 0, −r [J 0 (kr) + iJ 1 (kr)] , and so we can write The nonlinearity is computed up to quadratic order; we take the quadratic nonlinearity of (2.33)-(2.38), and use the definition of F in (2.42) to write F (u, ε, r) = dG G −1 (u) − 1 L(r)u + F 2 (u, ε, r) + O(|u| 3 ), where the subscript F 2 denotes elements of F with quadratic order. Computing the Fréchet derivative of G (see [22]), we can write the quadratic part of the nonlinearity as