Small-amplitude static periodic patterns at a fluid-ferrofluid interface
Mark D. Groves Jens Horn
math.AP







^1 Fachrichtung Mathematik, Universität des Saarlandes, Postfach 151150, 66041 Saarbrücken, Germany

^2Department of Mathematical Sciences, Loughborough University, Loughborough, LE11 3TU, UK



differential equations, fluid mechanics





M. D. Groves

groves@math.uni-sb.de



We establish the existence of static doubly periodic patterns (in particular rolls, rectangles and hexagons)
on the free surface of a ferrofluid near onset of the Rosensweig instability, assuming a general (nonlinear)
magnetisation law. A novel formulation of the ferrohydrostatic equations in terms of Dirichlet-Neumann operators
for nonlinear elliptic boundary-value problems is presented. We demonstrate the analyticity of these operators in
suitable function spaces and solve the ferrohydrostatic problem using an analytic version of Crandall-Rabinowitz
local bifurcation theory. Criteria are derived for the bifurcations to be sub-, super- or transcritical
with respect to a dimensionless physical parameter.











§ INTRODUCTION

Consider two static immiscible perfect fluids in the regions

    Ω^'    :={ (x,y,z): η(x,z)<y<d}, 
    Ω:     ={ (x,y,z): -d<y<η(x,z)}

separated by the free surface {y=η(x,z)}, where gravity acts in the negative y direction.
The upper fluid is non-magnetisable, while the lower is a ferrofluid with a general nonlinear
magnetisation law

    𝐌=𝐌(𝐇)=m(|𝐇|)𝐇/|𝐇|

expressing the relationship between the magnetisation
𝐌 of the ferrofluid and the strength of the magnetic field 𝐇.
Subjecting the fluids to a vertically directed magnetic field of sufficient strength leads to the emergence of
interfacial patterns (see Figure <ref>).




Small-amplitude static periodic patterns at a fluid-ferrofluid interface
    
M. D. Groves^1,2 and J. Horn^1
    December 30, 2023
========================================================================







In this article we present an existence theory for small amplitude, doubly periodic patterns with 

    η(𝐱 + 𝐥)=η(𝐱)

for every 𝐥∈ℒ, where 𝐱=(x,z) and ℒ is the lattice given by

    ℒ={  m𝐥_1+n𝐥_2: m,n∈ℤ}

with | 𝐥_1|=| 𝐥_2|. 
We are especially interested in three patterns which are observed in experiments (Figure <ref>),
namely rolls, rectangles and hexagons (see Figure <ref>).



  * 

For rolls we seek  functions that are independent of the z-direction   
and choose 𝐥=(2π/ω,0), so that
the periodic base cell  
is given by 
{  x:|x|<π/ω}.

  * For rectangles we choose 𝐥_1=(2π/ω,0), 𝐥_2=(0,2π/ω),
so that the periodic base cell
is given by 
{  (x,z):|  x|, |  z|<π/ω}.

  * For hexagons we choose 𝐥_1=2π/ω(1,-1/√(3)), 
𝐥_2=2π/ω(0,2/√(3)), 
so that we obtain an additional periodic direction 
𝐥_3=𝐥_1+𝐥_2=2π/ω(1,1/√(3))  
and the periodic base cell  
is given by 
{  
(x,z):|  x|<2π/ω, | x-√(3)z|<4π/ω,
| x+√(3)z|<4π/ω}.


Notice that each of these patterns exhibits a rotational symmetry: the shape of the free surface
is invariant under a rotation of the (x,z)-plane through respectively (i) π/2, (ii) π/3 and (iii) π.



In Section <ref><ref> 
we derive the mathematical formulation of this problem from physical principles.
The governing equations (equations (<ref>)–(<ref>)) are formulated
in terms of perturbations ϕ^' and ϕ of magnetic potentials corresponding 
to a uniform vertically directed magnetic field (of strength hμ(h) in the upper fluid and
h in the lower fluid, where
μ is obtained from the magnetisation law by the formula
μ(s)=1+m(s)/s; the potentials
are horizontally doubly periodic, satisfying

    ϕ^'(𝐱 + 𝐥,y)=ϕ^'(𝐱,y)
    ,
        ϕ(𝐱 + 𝐥,y)=ϕ(𝐱,y)

for every 𝐥∈ℒ (with a slight abuse of notation).



This problem was first studied by Cowley and Rosensweig <cit.>. 
Using a linear stability analysis, they found that, 
as the strength h of the magnetic field 
exceeds a critical value h_c, the flat surface destabilises and  a hexagonal pattern of peaks appears. 
This phenomenon is known as the Rosensweig instability. 
A mathematically rigorous treatment of the problem was given by Twombly and Thomas <cit.>, who 
used coordinate transformations to `flatten' the free surface by
transforming  the a priori unknown domains Ω^' and Ω into fixed strips. Applying
Lyapunov-Schmidt reduction reduces these transformed equations for rotationally symmetric patterns (see below)
to a locally equivalent one-dimensional equation 
which is solved using the implicit-function theorem; 
the result is the existence, for values of h near h_c, of rolls and rectangles in addition to the hexagonal pattern.    
Twombly and Thomas's work is however flawed by some miscalculations and mathematical inconsistencies, and is also restricted to
linear magnetisation laws. In this article we present a more systematic approach 
which is motivated by the corresponding study of doubly periodic travelling water waves 
by Craig and Nicholls <cit.>; we also consider general nonlinear magnetisation laws.


We work with dimensionless variables,  in terms of which the problem depends upon
two dimensionless parameters β (whose value β_0 is fixed) and γ (see equation (<ref>)), and `flatten' the equations using Dirichlet-Neumann formalism.
The Dirichlet-Neumann operator G^' for the upper fluid domain 
(given by {η(x,z)<y<1β_0} in dimensionless variables)
is defined as follows.
Fix Φ^'=Φ^'(x,z), solve the linear boundary-value problem

    ϕ^'_xx +ϕ^'_yy+ϕ^'_zz   =0,      η<y<1β_0,
    
    ϕ^'   =Φ^',       y=η, 
    
    ϕ^'_y   =0,       y=1β_0,

and define

    G^'(η,Φ^')    =-(1+η_x^2+η_z^2)^1/2ϕ^'_n|_y=η
     =-(
    ϕ^'_y-η_xϕ^'_x-η_zϕ^'_z
    )|_y=η.

The Dirichlet-Neumann operator G for the lower fluid domain
{-1β_0 < y < η(x,z)} is similarly defined as

    G(η,Φ)    =(1+η_x^2+η_z^2)^1/2μ(|(ϕ+y)|)
    ϕ_n|_y=η
        =μ(|(ϕ+y)|)(
    ϕ_y-η_xϕ_x-η_zϕ_z
    )|_y=η,

where ϕ is the solution of the (in general nonlinear) boundary-value problem

    div(μ(|(ϕ+y)|)(ϕ+y))   =0,       -1β_0<y<η,
    
    ϕ   =Φ,        y=η,
    
    μ(|(ϕ+y)|)(ϕ_y+1)   =μ(1),        y=-1β_0.

The nonlinearity of  (<ref>)–(<ref>)
is inherited from that of the magnetisation law 𝐌=𝐌(𝐇) (for a linear magnetisation law the value
of μ is constant and (<ref>), (<ref>) are replaced by respectively Laplace's equation and
a linear Neumann boundary condition). In Sections
<ref><ref> and <ref> we show that G^' and G are
analytic functions of respectively (η,Φ^') and (η,Φ) in suitable function spaces and use these operators to
recast the governing
equations in terms of the variables Φ^'=ϕ^'|_y=η and Φ=ϕ|_y=η.
The 
mathematical problem is thus to solve a system of equations of the form

    𝒢(γ,(η,Φ^',Φ))=0,

where 𝒢: ℝ× X_0 → Y_0
is given explicitly by the left-hand sides of equations (<ref>)–(<ref>)
and the function spaces X_0, Y_0 are specified in equation (<ref>).
Observe that this problem exhibits rotational symmetry: it
is invariant under rotations through respectively π, π/2 and π/3 for
rolls, rectangles and hexagons,
and one may therefore replace X_0 and Y_0 by their subspaces of functions that are invariant under these rotations
(denoted by  X_sym and Y_sym). 

In Section <ref> we discuss the existence of small-amplitude solutions
to (<ref>) within the framework of analytic
Crandall-Rabinowitz local bifurcation theory 
(see Buffoni and Toland<cit.>), using γ
as a bifurcation parameter.
According to that theory values γ_0 of
γ at which non-trivial solutions bifurcate from
zero (clearly 𝒢(γ,0)=0 for all values of γ)
necessarily have the property that the kernel of the linear operator
L_0:=d_2𝒢[γ_0, 0]: X_0 → Y_0 is non-trivial.
We show that L_0 is non-trivial if and only if 

    γ_0= r(| k|):=(μ_1(μ_1-1)^2(μ_1|𝐤| |𝐤|/β_0 +S_1 |𝐤| S_1|𝐤|/β_0)^-1-1)|𝐤|^2

for some 𝐤∈ℒ^⋆∖{0},  
where μ_1=μ(1), μ̇_1=μ̇(1) and
S_1=(μ_1/(μ_1+μ̇_1))^1/2. Choosingβ_0<μ_1(μ_1-1)^2/(μ_1+1) and ω so that

    γ_0 = (μ_1(μ_1-1)^2(μ_1ωω/β_0
    +S_1 ωS_1ω/β_0)^-1-1)ω^2

is the unique maximum of the mapping |𝐤| ↦ r(|𝐤|), we find that

    L_0
    =
    ⟨{𝐯sin(𝐤·𝐱),𝐯cos(𝐤·𝐱): 𝐤∈ℒ^⋆ with |𝐤|=ω}⟩,

where v∈ℝ^3 is given by equation (<ref>)
(see the discussion to Figure <ref>); this value
(β_0,γ_0)
of (β,γ) corresponds to the Rosensweig instability. The dimension 
of L_0 is therefore determined by the number of vectors in ℒ^⋆ with length ω;
for rolls, rectangles and hexagons we find that L_0 is respectively 2, 4 and
6 (see Figure <ref> and Sattinger <cit.> for a general discussion of this point).  Because the kernel of L_0 is  multidimensional, one can not use Crandall-Rabinowitz local bifurcation theory directly.
To overcome this problem we replace X_0 and Y_0 by X_sym and Y_sym, thus restricting to
solutions that are invariant under rotations
through respectively π, π/2 and π/3 for rolls, rectangles and hexagons.
These restrictions ensure that L_0=1 
with L_0=⟨v_0⟩, 
where  v_0=𝐯 e_1(x,z) 
and

    e_1(x,z) = {[                                   cosω x                                         ;                            cosω x+cosω z                                         ; cosω x
    +cosω/2(x+√(3)z)
    +cosω/2(x-√(3)z)                                          ].


Verifying the remaining conditions in the analytic Crandall-Rabinowitz local bifurcation theorem yields the following result.
 
The point (γ_0, 0) is a local bifurcation point for (<ref>), 
 that is 
 there exist ε >0, open neighbourhoods W_sym of (γ_0, 0) in ℝ× X_sym and V_sym of 0 in X_sym
and  analytic functions 
w:(-ε,ε)→ V_sym, γ: (-ε,ε)→ℝ with  
γ(0)=γ_0, w(0)=v_0 
such that 𝒢(γ(s),sw(s))=0 for every s∈ (-ε,ε). 
Furthermore   

    W_sym∩ N
    =
    {  (γ(s),sw(s)):0<|s|<ε},

where

    N={  (γ,v)∈ℝ×(V_sym∖{  0}):𝒢(γ,v)=0}.

 

In Section <ref> we examine the bifurcating branches identified in
Theorem <ref>.

 
Branches of small-amplitude doubly periodic solutions to the ferrohydrostatic problem
bifurcate from the trivial solution at γ=γ_0. The bifurcation is



  * transcritical in the case of hexagons,

  * super- or subcritical in the case of rolls and rectangles, depending upon the sign
of a coefficient γ_2 which is determined by μ and ω/β_0.



Explicit formulae for the coefficient γ_2 are given in some special cases
in Section <ref> (such formulae are unwieldy, and it appears in
general more appropriate to calculate them numerically for a specific choice of μ).
We note in particular that for constant μ (corresponding to a linear magnetisation law)
and very deep fluids, rolls bifurcate subcritically for μ<μ_c^1 and
supercritically for μ>μ_c^1, while rectangles bifurcate subcritically for
μ<μ_c^2 and supercritically for μ>μ_c^2, where

    μ_c^1=21/11+8/11√(5),
        μ_c^2=115+160√(2)+8√(184+11√(2))/141+128√(2).

The same values were obtained by Silber and Knobloch <cit.> in a discussion of normal forms for this bifurcation problem and confirmed by Lloyd, Gollwitzer, Rehberg 
and Richter <cit.>  
as part of a wider numerical and experimental investigation.

Finally, we note that supercritical bifurcation of rolls is associated with (supercritical) bifurcation of
spatially localised patterns, whose existence has been established by dynamical-systems
arguments by Groves, Lloyd & Stylianou <cit.>.




§ MATHEMATICAL FORMULATION




 §.§ The physical problem





We consider two static immiscible perfect fluids in the regions

    Ω^':={(x,y,z): η(x,z)<y<d},     Ω:={(x,y,z): -d<y<η(x,z)}

separated by the free interface {y=η(x,z)}.
The upper, non-magnetisable fluid has unit relative permeability and density ρ^',
while the lower is a ferrofluid with density ρ.
The relations between the magnetic fields 𝐇^', 𝐇 and the induction fields 𝐁^', 𝐁  
are given by the identities 

    μ_0𝐇^' =𝐁^',
            μ_0(𝐇+𝐌(𝐇))
    =𝐁

where μ_0 is the vacuum permeability 
and 𝐌 is the 
magnetic intensity of the ferrofluid. (Here, and in the remainder of this paper, equations for `primed' and `non-primed' quantities are supposed to hold in respectively Ω^' and Ω.)  
We suppose that 

    𝐌(𝐇)=m(|𝐇|)𝐇/|𝐇|,

where m is a nonnegative function, so that in particular 𝐌 and 𝐇 are collinear. 
According to Maxwell's equations the magnetic and induction fields are respectively irrotational and solenoidal, 
and introducing 
magnetic potential functions ϕ^', ϕ with

    -ϕ^' =𝐇^',
            
    -ϕ =𝐇,

one finds that these potentials satisfy the equations 

    (ϕ^')=0
    ,
            (μ(| ϕ|) ϕ)=0,

in which  

    μ(s)=1+m(s)/s;

we assume that μ:(0,∞)→ℝ is analytic and satisfies μ(1) + μ̇(1)>0
(so that the linearised version of the equation for ϕ is elliptic). 
Observe that ϕ^' is harmonic while ϕ satisfies a nonlinear elliptic partial differential equation;
this nonlinearity is inherited from that of the magnetisation law 𝐌=𝐌(𝐇)
(for a linear magnetisation law the value of μ is constant and the equation for ϕ reduces
to Laplace's equation).
At the interface we have the magnetic conditions 

    𝐇^'·𝐭_1 =𝐇·𝐭_1,
            𝐇^'·𝐭_2 =𝐇·𝐭_2,
            𝐁^'·𝐧 =𝐁·𝐧,

where 

    𝐭_1   =
    (1,η_x,0)^T√(1+η^2_x)
    ,
            𝐭_2
    =
    (0,η_z,1)^T√(1+η^2_z)
    ,
            𝐧=(-η_x,1,-η_z)^T√(1+η^2_x+η^2_z)

are the tangent and normal vectors to the interface; 
it follows that 

    ϕ^'-ϕ|_y=η =0
    ,
            ϕ^'_n
    -μ(| ϕ|) ϕ_n|_y=η  =0.


The ferrohydrostatic Euler equations are given by 

    -(p'_0+ρ^' gy) =0  
    ,
            μ_0(𝐌(𝐇)·
    (∂_x,∂_z,∂_y)^T) 𝐇
    -(p^⋆+
    ρ gy) =0

(Rosensweig <cit.>),
in which   
g is the acceleration due to gravity, 
p'_0 is the pressure in the upper fluid and
p^⋆ is a composite of 
the magnetostrictive and  fluid-magnetic pressures.
The calculation 

    (𝐌(𝐇)·
    (∂_x,∂_y,∂_z)^T)
    𝐇
    =
    m(|𝐇|)  |𝐇|
    =
    (
    ∫_0^|𝐇|
    m(t)
     d t
    )

shows that these equations are equivalent to

    -(p'_0+ρ^' gy) =C'_0 
    , 
            μ_0∫_0^|𝐇|
    m(t)
     d t
    -(p^⋆
    +ρ gy) =C_0
    ,

where  C'_0,C_0 are constants. 
The ferrohydrostatic boundary condition 
is given by 

    p^⋆
    +μ_0/2
    (𝐌(𝐇)·𝐧)^2= p'_0+2κσ

for y=η(x,z) (Rosensweig <cit.>), 
in which σ>0 is the coefficient of surface tension and 

    2κ    =-η_xx(1+η_z^2)+η_zz(1+η_x^2)
    -2η_xη_zη_xz(1+η_x^2+η_z^2)^3/2

is the mean curvature of the interface. 
Using (<ref>), we find that

    μ_02(𝐌(𝐇)·𝐧)^2
    +μ_0∫_0^|𝐇|
    m(t)
     d t
    +C
    + (ρ^'-ρ)gη
    -2σκ    = 0
    ,

where C=C'_0-C_0, or equivalently  

    C
    +μ_0√(1+η_x^2+η_z^2)
    (ϕ^'_yϕ^'_n
    -μ(| ϕ|)ϕ_yϕ_n) 
    
           
    - (ρ-ρ^') gη
    -2σκ 
    -μ_0(1/2| ϕ^'|^2-M(| ϕ|))
    =0

for y=η(x,z) with

    M(s)
    =
    ∫_0^stμ(t) d t.
 
The requirement that a uniform magnetic field and flat interface solves the physical problem, that is
(η_0,ϕ^'_0,ϕ_0)=(0,μ(h) h y,h y)
is a solution to (<ref>), (<ref>) and (<ref>),  
leads us to chooseC=-μ_0M(h)-μ_0μ(h)(μ(h)-1)h^2/2;
we write (η̃,ϕ̃^',ϕ̃) = (η,ϕ^',ϕ)+(η_0,ϕ^'_0,ϕ_0)
(so that (η̃,ϕ̃^',ϕ̃)=(0,0,0) is the `trivial' solution)
and drop the tildes for notational simplicity.

The next step is to introduce dimensionless variables

    (x̂,ŷ,ẑ)=μ_0h^2/σ(x,y,z),
      η̂:= μ_0h^2/ση,
      ϕ̂^':= μ_0h/σϕ^',
      ϕ̂:= μ_0h/σϕ

and functions

    μ̂(s):= μ(hs),    M̂(s):=M(hs).

We find that

    ϕ^'_xx+ϕ^'_yy+ϕ^'_zz    =0,         η(x,z)<y<1β_0,
    (μ(|(ϕ +y)|)(ϕ +y))     =0,          -1β_0 < y < η(x,z),

with boundary conditions

    ϕ^'_y    = 0,          y=1β_0, 
    
    μ(|(ϕ +y)|)( ϕ_y+1)-μ(1) 
        = 0,          y=-1β_0,
    ϕ^'-ϕ
    +(μ(1)-1)η    =0,        
    y=η(x,z),
    
    
    (ϕ^'+μ(1) y)_n
    -μ(|(ϕ +y)|)(ϕ +y)_n    =0,        
    y=η(x,z),

and

    -M(1)-μ(1)(1/2μ(1)-1)- γη
    -2σκ
         -1/2|(ϕ^'+μ(1) y)|^2+M(|(ϕ+y)|)
         
    +
    √(1+η_x^2+η_z^2)
    (ϕ^'_y+μ(1))
    (ϕ^'+μ(1) y)_n
         
    -√(1+η_x^2+η_z^2)μ(|(ϕ +y)|)
    (ϕ_y+1)
    (ϕ+y)_n
     = 0,         
    y=η(x,z),

where

    α=(ρ-ρ^')gd/μ_0h^2,     β=σ/μ_0h^2d,     γ=αβ

 and the hats have been dropped for notational simplicity.
 
We seek periodic solutions to (<ref>)–(<ref>) satisfying

    η(𝐱 + 𝐥)=η(𝐱)
    ,
        ϕ^'(𝐱 + 𝐥,y)=ϕ^'(𝐱,y)
    ,
        ϕ(𝐱 + 𝐥,y)=ϕ(𝐱,y)

for every 𝐥∈ℒ (with a slight abuse of notation),
where 𝐱=(x,z) and ℒ is the lattice given by

    ℒ={ m𝐥_1+n𝐥_2: m,n∈ℤ}

with |𝐥_1|=|𝐥_2|. 
Choose 𝐤_1,𝐤_2 with 𝐤_i·𝐥_j=2πδ_ij for i,j=1,2 and define the dual lattice ℒ^⋆ to ℒ by

    ℒ^⋆={m𝐤_1+n𝐤_2: m,n∈ℤ},

 so that our periodic functions can be written as 
 
    η(𝐱)
    =
    ∑_𝐤∈ℒ^⋆η__𝐤e^i𝐤·𝐱,
        ϕ^'(𝐱,y)
    =
    ∑_𝐤∈ℒ^⋆ϕ^'__𝐤(y)e^i𝐤·𝐱,
        ϕ(𝐱,y)
    =
    ∑_𝐤∈ℒ^⋆ϕ__𝐤(y)e^i𝐤·𝐱
    ,

where η__-𝐤=η̅__𝐤,ϕ^'__-𝐤=ϕ̅^̅'̅__𝐤,ϕ__-𝐤=ϕ̅__𝐤. 
We are especially interested in three periodic patterns, namely rolls, rectangles and hexagons
(see Figure <ref>).


  * 

For rolls we seek  functions that are independent of the z-direction   
and we choose𝐥=(2π/ω,0),  so that the dual lattice ℒ^⋆ is 
generated by 𝐤=(ω,0) 
and the periodic base cell  
is given by 

    { 
    x:| x|<π/ω}.

Furthermore, the z-independent versions of  equations (<ref>)–(<ref>) 
are invariant under the  reflection x↦ -x (which corresponds to a rotation through π in the (x,z)-plane).

  * For rectangles we choose 𝐥_1=(2π/ω,0) and 𝐥_2=(0,2π/ω),  
so that the dual lattice ℒ^⋆ is generated by 𝐤_1=(ω,0)  and 𝐤_2=(0,ω)  
and the periodic base cell  
is given by 

    { 
    (x,z):| x|, | z|<π/ω}.

Furthermore, equations (<ref>)–(<ref>) 
are invariant under rotations through  π/2 in the (x,z)-plane.

  * For hexagons we choose 𝐥_1=2π/ω(1,-1/√(3)) and 
𝐥_2=2π/ω(0,2/√(3)), 
so that we obtain an additional periodic direction 
𝐥_3=𝐥_1+𝐥_2=2π/ω(1,1/√(3)).
The dual lattice 
ℒ^⋆ is generated by 
𝐤_1=(ω,0)  and 𝐤_2=ω(1/2,√(3)/2)
and the periodic base cell  
is given by

    { 
    (x,z):|x|<2π/ω, |x-√(3)z|<4π/ω|x+√(3)z|<4π/ω}.
 
Furthermore, equations (<ref>)–(<ref>) 
are invariant under rotations through  π/3 in the (x,z)-plane.



The mathematical problem is thus to solve (<ref>)–(<ref>)  
for periodic functions η, ϕ^' and ϕ in the domains Γ, Ω_per^' and Ω_per,     
where 

    Ω_per^':=
    { 
    (x,y,z):(x,z)∈Γ}∩Ω^',
        Ω_per:=
    { 
    (x,y,z):(x,z)∈Γ}∩Ω

and Γ is the parallelogram defined by 𝐥_1 and 𝐥_2 (or by 𝐥 in the case of rolls).






 §.§ Dirichlet-Neumann formalism
 

The Dirichlet-Neumann operator G^' for the upper fluid domain 
{η(x,z)<y<1β_0}
is defined as follows.
Fix Φ^'=Φ^'(x,z), solve the linear boundary-value problem

    ϕ^'_xx +ϕ^'_yy+ϕ^'_zz   =0,      η<y<1β_0,
    
    ϕ^'   =Φ^',       y=η, 
    
    ϕ^'_y   =0,       y=1β_0 
    .

and define

    G^'(η,Φ^')    =-(1+η_x^2+η_z^2)^1/2ϕ^'_n|_y=η
     =-(
    ϕ^'_y-η_xϕ^'_x-η_zϕ^'_z
    )|_y=η.

The Dirichlet-Neumann operator G for the lower fluid domain
{-1β_0 < y < η(x,z)} is similarly defined as

    G(η,Φ)    =(1+η_x^2+η_z^2)^1/2μ(|(ϕ+y)|)
    ϕ_n|_y=η
        =μ(|(ϕ+y)|)(
    ϕ_y-η_xϕ_x-η_zϕ_z
    )|_y=η,

where ϕ is the solution of the (in general nonlinear) boundary-value problem

    div(μ(|(ϕ+y)|)(ϕ+y))   =0,       -1β_0<y<η,
    
    ϕ   =Φ,        y=η,
    
    μ(|(ϕ+y)|)(ϕ_y+1)   =μ(1),        y=-1β_0.

It is also convenient to introduce auxiliary operators H^' and H given by

    H^'(η,Φ^') = ϕ^'_y|_y=η,      H(η,Φ) = ϕ_y|_y=η,

where ϕ^' and ϕ are the solutions to the boundary-value problems (<ref>)–(<ref>)
and (<ref>)–(<ref>).
Using this Dirichlet-Neumann formalism, we can recast the governing
equations (<ref>)–(<ref>) in terms of the variables η,
Φ^'=ϕ^'|_y=η and Φ = ϕ|_y=η as 

    Φ^'- Φ +(μ(1)-1) η    =0,
    
    G^'(η,Φ^') + G(η,Φ)
    +(μ^⋆ - μ(1))    =0

and

    -γη+ 
    ∇·(∇η√(1+ |∇η|^2))
    +1/2(1+|∇η|^2)
    H^'(η,Φ^')^2
    -1/2|∇Φ^'|^2|_y=η
           
    -μ(1)G^'(η,Φ^') - G(η,Φ)-(μ^⋆ - μ(1))
    
           
    +(M^⋆-M(1))-μ^⋆H(η,Φ)-H(η,Φ) G(η,Φ)
    =0
    ,

in which 
∇=(∂_x,∂_z)^T, M(s)=∫_0^s tμ(t) dt,

    μ^⋆   =
    μ(
    (
    |∇Φ|^2
    +2(1-∇η·∇Φ)H(η,Φ)
    +(1+|∇η|^2)H(η,Φ)^2
    +1
    )^1/2),
    
    
    M^⋆   =
    M
    (
    (
    |∇Φ|^2
    +2(1-∇η·∇Φ)H(η,Φ)
    +(1+|∇η|^2)H(η,Φ)^2
    +1
    )^1/2).


We study equations (<ref>)–(<ref>) in the standard Sobolev spaces 

    H_per^r(Γ)
    =
    {ζ
    =
    ∑_𝐤∈ℒ^⋆ζ__𝐤e^i𝐤·𝐱: ζ_-𝐤=ζ_𝐤,
     ‖ζ‖_r
    <
    ∞},     ‖ζ‖_r^2 
    := C(Γ)∑_𝐤∈ℒ^⋆( 1 + |𝐤|^2 )^r 
    |ζ__𝐤|^2,

where ℒ^⋆ is the dual lattice to ℒ, and their subspaces 

    H̅_per^r(Γ) 
       = 
    {ζ∈ H^r_per(Γ) : 
    ζ__0=0
    }
  
consisting of functions with zero mean (the value of the normalisation constant C(Γ) is 2π/ω for rolls, 4(π/ω)^2 for rectangles and
8/√(3)(π/ω)^2 for hexagons).
In Section <ref><ref> below
we establish the following theorem. (A function is `analytic at the origin' if is defined
and analytic in a neighbourhood of the origin; in particular it has a convergent
Maclaurin series.)


Suppose that s>5/2. Formulae (<ref>), (<ref>) and (<ref>) define 
mappings G^', G: H^s_per(Γ) × H^s-1/2_per(Γ) →
 H^s-3/2_per(Γ) and H^', H: H^s_per(Γ) × H^s-1/2_per(Γ) →
 H^s-3/2_per(Γ) 
which are  analytic at the origin.

Define

    X_0 :=
    H^s+1/2_per(Γ)
    ×H̅^s_per(Γ)
    ×
    H^s_per(Γ),
        
     Y_0
    :=
    H^s_per(Γ)
    ×H̅^s-1_per(Γ)
    ×
    H^s-3/2_per(Γ)

for s>5/2. Using Theorem <ref> and the fact that H^r(Γ) is a Banach algebra for r>1, we find that
the left-hand sides of equations (<ref>)–(<ref>)  
define a function 𝒢:ℝ× X_0 → Y_0 which is analytic at the origin. (A straightforward calculation shows that

    ∫_Γ(G^'(η,Φ^') + G(η,Φ)+(μ^⋆ - μ(1)))=0

and explains the choice of functions with zero mean in the second component of Y_0; using functions
with zero mean in the second component of X_0 on the other hand ensures that the kernel of the linear operator
d_2𝒢[γ,0]  does not contain any constant terms for any γ∈ℝ).
The 
mathematical problem is thus to solve
    𝒢(γ,(η,Φ^',Φ))=0,

for (γ,(η,Φ^',Φ)) ∈ℝ× V_0, where V_0 is a neighbourhood of the origin in X_0 and
𝒢(γ, 0)=0 for all γ∈ℝ.  Observe that this problem exhibits rotational symmetry: it
is invariant under rotations through respectively π, π/2 and π/3
for rolls, rectangles and hexagons,
and one may therefore replace X_0 and Y_0 by their subspaces of functions that are invariant under these rotations
(denoted by  X_sym and Y_sym). 




 

 §.§ Analyticity of the Dirichlet-Neumann operators

 
 
We study the boundary-value problems (<ref>)–(<ref>) 
and (<ref>)–(<ref>) by transforming them into equivalent problems
in fixed domains (cf. Nicholls and Reitich <cit.> and 
Twombly and Thomas <cit.>). 
The change of variable

    ỹ
    =y-η1+β_0η,     
    u(x,ỹ,z) =ϕ(x,y,z)

transforms the variable domain Ω_per into the fixed domain   
Σ={(x,y,z): (x,z)∈Γ, y ∈ (-1β_0,0)}
and the boundary-value problem (<ref>)–(<ref>) into

    div(μ^†((u+y)-(F_1(η,u),F_3(η,u),F_2(η,u))^T))
       
    =0,
          -1β_0<y<0,
    
    
    u-Φ   
    =0,
          y=0
    ,
    
    μ^†((u+y).(0,1,0)^T-F_3(η,u))-μ(1)
       
    =0,
          y=-1β_0,

where

    F_1(η, u)  =-η u_x +(1+ β_0 y)η_x u_y,    
     F_2(η, u) =-η u_z +(1+ β_0 y)η_z u_y,


    F_3(η, u) =β_0η u_y/1+β_0η
    +(1+ β_0 y)(η_x u_x+ η_z u_z)
    -(1+β_0 y)^2/1+β_0η
    (η^2_x+η^2_z)u_y

and

    μ^† = μ(| 1/1+β_0η( u - (F_1(η,u),0,F_2(η,u))^T)+(0,1,0)^T|)

(we have again dropped the tildes for notational simplicity).



Suppose that s> 5/2.
There exist open neighbourhoods V and U of the origin in respectively
H_per^s(Γ) × H_per^s-1/2(Γ) and
H_per^s(Σ) such that
the boundary-value problem (<ref>)–(<ref>)
has a unique solution u=u(η,Φ) in U for each (η,Φ) ∈ V. Furthermore u(η,Φ) depends
analytically upon η and Φ.



Write the left-hand sides of equations (<ref>)–(<ref>)
as

    ℋ(u,η,Φ)=0,

and observe that ℋ: H_per^s(Σ) × H_per^s(Γ) ×
H_per^s-1/2(Γ) is analytic at the origin with

    ℋ(0,0,0)=0.

Furthermore, the calculation

    d_1ℋ[0,0,0](u)=[ μ_1 (u_xx+S_1^-2 u_yy+u_zz );                       u|_y=0;     μ_1 S_1^-2 u_y|_y=-1/β_0 ],

where
μ_1=μ(1), μ̇_1=μ̇(1) and
S_1=(μ_1/(μ_1+μ̇_1))^1/2,
and standard existence and regularity theory for elliptic linear boundary-value problems show that
d_1ℋ[0,0,0]: H_per^s(Σ) → H_per^s-2× H_per^s-1/2× H_per^s-3/2(Γ)
is an isomorphism. The stated result now follows from the analytic implicit-function theorem.


The transformation (<ref>) converts the formulae

    G(η,Φ)=μ(|(ϕ+y)|)(ϕ_y-η_xϕ_x-η_zϕ_z)|_y=η,     
    H(η,Φ)=ϕ_y|_y=η

into

    G(η,Φ) = μ^† (u_y-F_3(η,u))|_y=0,      H(η,Φ)=.u_y/1+β_0η|_y=0,

where u=u(η,Φ) is the unique solution to (<ref>)–(<ref>).

Suppose that s>5/2. The formulae (<ref>) define analytic functions G,
H: V → H_per^s-3/2(Γ).



To compute the Taylor-series representations of u and G we begin with the functionν: (H^s-1_per(Σ))^3→ H^s-1_per(Σ) defined by

    ν(T)=μ(|T+(0,1,0)^T|).

Observing that ν is analytic at the origin, we write its Taylor series as

    ν(T)
    =
    ∑_j=0^∞ν^j({T}^(j)),

where ν^j∈ℒ^j_s((H^s-1_per(Σ))^3,H^s-1_per(Σ)) is given by

    ν^j(T_1,…,T_j) = 1/j!d^jν[0](T_1,…,T_j)

and may be computed explicitly from μ (note in particular that ν^0=μ_1).
The functionsu^n ∈ℒ^n_s(H^s_per(Γ)× H^s-1/2_per(Γ),H^s_per(Σ))  
(with u^0=0) in the corresponding series

    u(η,Φ)
    =
    ∑^∞_n=0
    u^n({(η,Φ)}^(n))
  
may be computed recursively by substituting the Ansatz (<ref>), (<ref>)
into equations (<ref>)–(<ref>).
Consistently abbreviating m^n({(η,Φ)}^(n)) to m^n for notational simplicity, one finds after a lengthy but straightforward calculation that  

    (L  u^1)   =0,    (L  u^n-F^n)    =0,          -1β_0<y<0, 
    
    u^1-Φ   =0,     u^n    = 0,          y=0, 
    
    u^1_y   =0,     (L  u^n-F^n)·(0,1,0)^T   =0,          y=-1β_0

for n ≥ 2,  where 
 
    L
       =
    [      1      0      0;      0 S_1^-2      0;      0      0      1 ],
    
    μ_1F^n   =
    μ_1
    (
    (F^n_1,F^n_3,F^n_2)^T-
    μ̇_1/μ_1∑_j=1^n(-β_0η)^ju^n-j_y(0,1,0)^T)
           -
    ∑_j=0^nν^1(T^j)(  u^n-j-(F^n-j_1,F^n-j_3,F^n-j_2)^T)
    
           
    -
    R^n(0,1,0)^T
    -∑_j=0^n
    R^j(  u^n-j-(F^n-j_1,F^n-j_3,F^n-j_2)^T)

and
    F^n_1    =-β_0η u^n-1_x +(1+ β_0 y)η_x u^n-1_y,    
     F^n_2=-β_0η u^n-1_z +(1+ β_0 y)η_z u^n-1_y,
    
     F^n_3    =β_0η∑^n-1_j=0(-β_0η)^ju^n-1-j_y
    +(1+β_0 y)(η_x u^n-1_x+ η_z u^n-1_z)
    
           
    -(1+β_0 y)^2
    (η^2_x+η^2_z)
    ∑^n-2_j=0(-β_0η)^ju^n-2-j_y,
    
    
    T^n   =
    ∑_j=0^n(-β_0 η)^j(
      u^n-j-(F^n-j_1,0,F^n-j_2)^T)
    ,
        
    R^n
    =∑_2≤ j≤ n, 
    
    h_1+…+h_j=nν^j(
    T^h_1,…,T^h_j
    )
    .


The Taylor-series representations of G and H are thus given by  

    G(η,Φ)
    =
    ∑^∞_n=0
    G_n({(η,Φ)}^(n)),
        
    H(η,Φ)
    =
    ∑^∞_n=0
    H_n({(η,Φ)}^(n)),

where 

    G_n
    =
    .
    μ_1
    I^n
    +∑_j=0^n ν^1(T^j)I^n-j+∑_j=0^n R^jI^n-j|_y=0 ,     
    H_n
    =
    ∑_j=0^n
    (-β_0η)^ju^n-j_y|_y=0 ,

and 

    I^n   = 
    ∑_j=0^n(-β_0 η)^j u_y^n-j
    +∑_j=0^n-2(-β_0η)^j 
    (η^2_x+η^2_z)
    u_y^n-j-2
    -(η_x u^n-1_x+ η_z u^n-1_z).
 
For later use we record the formulae

    G_1   =
    μ_1
    u_y^1|_y=0,
    
    
    G_2   =
    μ_1
    (
    ∑_j=0^2(-β_0η)^j 
    u_y^2-j
    -(η_x u^1_x+ η_z u^1_z)
    )
    +
    μ̇_1
    (u_y^1)^2|_y=0,
    
    
    G_3   =
    μ_1
    (
    ∑_j=0^3(-β_0η)^j u_y^3-j
    + 
    (η^2_x+η^2_z)
    u_y^1
    -
    (η_x u^2_x+ η_z u^2_z)
    )
    
           
    +
    μ̇_1
    (
    2∑_j=0^2(-β_0η)^j 
    u_y^2-j
    -(η_x u^1_x+ η_z u^1_z)
    ) u^1_y
    +
    12μ̈_1
    (u_y^1)^3
    +12μ̇_1((u_x^1)^2+(u_z^1)^2)u_y^1
    |_y=0,

where μ̈_1=μ̈(1), and

    H_1=u^1_y|_y=0,
        
    H_2=∑_j=0^2(-β_0η)^ju^2-j_y|_y=0,
        
    H_3=∑_j=0^3(-β_0η)^ju^3-j_y|_y=0

for the first few terms in these series.

The boundary-value problem (<ref>)–(<ref>) is handled in a similar fashion. The change of variable

    ỹ
    =y-η/1-β_0η,     
    u^'(x,ỹ,z) =ϕ^'(x,y,z)

transforms the variable domain Ω_per^' into the fixed domain   
Σ^'={(x,y,z):(x,z)∈Γ, y ∈(0,1β_0)}
and (<ref>)–(<ref>) into

    div ( u^'-(F^'_1(η,u^'),F^'_3(η,u^'),F^'_2(η,u^'))^T)
       
    =0,
              
    0<y<1β_0,
    
    
    u^'-Φ^'   
    =0,
          y=0
    ,
    
    
    u^'_y
       
    =0,
          y=1β_0
    ,

where

    F^'_1(η, u^') =η u^'_x +(1-β_0 y)η_x u^'_y,    
     F^'_2(η, u^')=η u^'_z +(1-β_0 y)η_z u^'_y,


    F^'_3(η, u^') =-β_0η u^'_y/1-β_0η
    +(1-β_0 y)(η_x u^'_x+ η_z u^'_z)
    -(1-β_0y)^2/1-β_0η
    (η^2_x+η^2_z)u^'_y

(we have again dropped the tildes for notational simplicity). The formulae

    G^'(η,Φ^') =-(ϕ^'_y-η_xϕ^'_x-η_zϕ^'_z)|_y=η,     
    H^'(η,Φ^')=ϕ^'_y|_y=η

are converted into

    G^'(η,Φ^') = -u_y^'+F_3^'(η,u^')|_y=0,     
    H^'(η,Φ^') = .u_y^'/1-β_0η|_y=0,

where u^'=u^'(η,Φ^') is the unique solution to (<ref>)–(<ref>).



Suppose that s>5/2.



  (i)
There exist open neighbourhoods V^' and U^' of the origin in respectively
H_per^s(Γ) × H_per^s-1/2(Γ) and
H_per^s(Σ^') such that
the boundary-value problem (<ref>)–(<ref>)
has a unique solutionu^'=u^'(η,Φ^') in U^' for each (η,Φ^') ∈ V^'. Furthermore
u^'(η,Φ^') depends analytically upon η and Φ^'.

  (ii)
The formulae (<ref>)
define analytic functions G^', H^': H^s_per(Γ)× H^s-1/2_per(Γ) →
 H^s-3/2_per(Γ).



The functions u^' n∈ℒ^n_s(H^s_per(Γ)× H^s-1/2_per(Γ),H^s_per(Σ^')) 
(with u^' 0=0) in the Taylor series

    u'(η,Φ^')
    =
    ∑^∞_n=1
    u^' n({(η,Φ^')}^(n))

are computed recursively by substituting this Ansatz into equations (<ref>)–(<ref>);
one finds that

    ( u^' 1)   =0,    ( u^' n-(F^' n_1,F^' n_3,F^' n_2)^T)    =0,          0<y<1β_0, 
    
    u^' 1-Φ^'   =0,     u^' n    = 0,          y=0, 
    
    u^' 1_y   =0,     u^' n_y   =0,          y=1β_0

for n≥ 2, where

    F^' n_1=η u^' n-1_x +(1-β_0 y)η_x u^' n-1_y,    
     F^' n_2 =η u^' n-1_z +(1-β_0 y)η_z u^' n-1_y,


    F^' n_3 =   -β_0η∑^n-1_j=0(β_0η)^ju^' n-1-j_y
    +(1-β_0 y)(η_x u^' n-1_x+ η_z u^' n-1_z) 
           
    -(1-β_0y)^2(η^2_x+η^2_z)∑^n-2_j=0(β_0η)^ju^' n-2-j_y.

The Taylor-series representations of G^' and H^' are given by 

    G^'(η,Φ^')
    =
    ∑_n=1^∞
    G^'_n({(η,Φ^')}^(n)),
        
    H^'(η,Φ^')
    =
    ∑_n=1^∞
    H^'_n({(η,Φ^')}^(n))

with 

    G^'_n   = 
    .-
    ∑_j=0^n-1(
    β_0η)^j u_y^' n-j
    -∑_j=0^n-3(
    β_0η)^j 
    (η^2_x+η^2_z)
    u_y^' n-j-2
    +(η_x u^' n-1_x+ η_z u^' n-1_z)
    |_y=0, 
    
    H^'_n   = 
    ∑_j=0^n-1(
    β_0η)^j u_y^' n-j|_y=0,

and in particular we find that

    G^'_1= -u_y^'1|_y=0,
        
    G^'_2= -.∑_j=0^1(β_0η)^j u_y^'2-j
    +(η_x u^'1_x+ η_z u^'1_z)|_y=0,


    G^'_3= -.∑_j=0^2(β_0η)^j u_y^'3-j
    -(η^2_x+η^2_z)u_y^'1+(η_x u^'2_x+ η_z u^'2_z)|_y=0 .

and

    H^'_1=u_y^'1|_y=0,
        
    H^'_2= .∑_j=0^1(β_0η)^j u_y^'2-j|_y=0,
        
    H^'_3 = .∑_j=0^2(β_0η)^j u_y^'3-j|_y=0.





§ EXISTENCE THEORY


Next we introduce the Crandall-Rabinowitz theorem (cf. Buffoni and Toland <cit.>), 
an application of which yields  a local bifurcation point of the equation

    𝒢(γ,(η,Φ^',Φ))=0,

where the components of 𝒢: ℝ× V_0 → Y_0 are given by
the left-hand sides of (<ref>)–(<ref>). 


Let 
X and Y be Banach spaces, 
V be an open neighbourhood of the origin in X   
and 
ℱ:ℝ× V→ Y be an analytic function 
with ℱ(λ, v) =0 for all λ∈ℝ. 
Suppose also that 




  * L:=d_2ℱ[λ_0, 0]: X → Y is a Fredholm operator of index zero, 

  * L=⟨   v_0⟩ for some v_0∈ X,  

  * the transversality condition
P(d_1d_2ℱ[λ_0,0](1,v_0))≠ 0 
holds, 
where 
P: Y→ Y is a projection 
with L = P.



The point (λ_0,0) is a local bifurcation point, 
that is
there exist ε >0, 
an open neighbourhood W of (λ_0,0) in ℝ× X
and  analytic functions 
w:(-ε,ε)→ V, λ: (-ε,ε)→ℝ with  
λ(0)=λ_0, w(0)=v_0 
such that ℱ(λ(s),sw(s))=0 for every s∈ (-ε,ε). 
Furthermore   

    W∩ N={ (λ(s),sw(s)):0<|s|<ε},

where

    N={(λ,v)∈ℝ×(V ∖{ 0 }):ℱ(λ,v)=0 }.



The first step is to determine the maximal positive  value γ_0 
of the parameter γ 
for which the kernel of the linear operator 
L_0:=d_2𝒢[γ_0,0]:X_0 → Y_0, which is
given by the explicit formula    

    L_0[   η; Φ^';   Φ ]
     =[                                        Φ^'-Φ
    +(μ_1-1)η;                 G_1^'(η,Φ^') + μ G_1(η,Φ)+μ̇_1H_1(η,Φ); η_xx+η_zz-γ_0η-(μ_1G_1^'(η,Φ^')+G_1(η,Φ)+μ̇_1H_1(η,Φ)) ]

with

    G_1^'(η,Φ^')
    =
    ∑_𝐤∈ℒ^⋆
    |𝐤|tanh|𝐤|/β_0Φ^'__𝐤e^i𝐤·𝐱
    ,
        
    G_1(η,Φ)
    =
    μ_1∑_𝐤∈ℒ^⋆
    S_1|𝐤|tanhS_1|𝐤|/β_0Φ__𝐤e^i𝐤·𝐱

and G_1(η,Φ)=μ_1H_1(η,Φ),
is non-trivial.  Writing  v ∈ X_0 
as 

    v(𝐱)=∑_𝐤∈ℒ^⋆𝐯__𝐤e^i𝐤·𝐱

with   
𝐯__𝐤=(η__𝐤,Φ^'__𝐤,Φ__𝐤)^T
and v__-𝐤=v̅__𝐤, we find that    

    L_0v
       =
    ∑_𝐤∈ℒ^⋆
    L_0(|𝐤|)
    𝐯__𝐤e^i𝐤·𝐱,

where 

    L_0(|𝐤|)
    =
    [                         μ_1-1                             1                            -1;                             0                |𝐤|tanh |𝐤|β_0  μ_1S_1^-1  |𝐤| tanhS_1|𝐤|β_0;                    -|𝐤|^2-γ_0            -μ_1|𝐤|tanh |𝐤|β_0 -μ_1S_1^-1  |𝐤| tanhS_1|𝐤|β_0 ]

for 𝐤≠ 0 and  

    L_0(0)
       =
    [ μ_1-1     1    -1;     0     0     0;  -γ_0     0     0 ]

(where we have identified the subspace 
{ 
(η__0,Φ^'__0,Φ__0)^T
:Φ^'__0=0} of ℝ^3 with ℝ^2).

From this observation it follows that 
L_0 is non-trivial if 

    L_0(|𝐤|)
       =μ_1(μ_1-1)^2S_1^-1|𝐤|^2 tanh|𝐤|/β_0tanhS_1|𝐤|/β_0
           -(|𝐤|^2+γ_0)(μ_1|𝐤| tanh|𝐤|/β_0 +S_1 |𝐤| tanhS_1|𝐤|/β_0) =0,

that is

    γ_0=r(| k|):= (μ_1(μ_1-1)^2(μ_1|𝐤| |𝐤|/β_0 +S_1 |𝐤| S_1|𝐤|/β_0)^-1-1)|𝐤|^2

for some 𝐤∈ℒ^⋆∖{0}. The function |𝐤| ↦ r(| k|), which satisfies r(0)=0 and r(|𝐤|) →-∞as |𝐤| →∞, takes only negative values for
β_0>μ_1(μ_1-1)^2/(μ_1+1), while forβ_0<μ_1(μ_1-1)^2/(μ_1+1) it has a unique maximum ω with
r(ω)>0 (see  Figure <ref>); we choose γ_0=r(ω) and
note the relationships

    β_0=μ_1(μ_1-1)^2/2ω̃( h(ω̃)-ω̃ḣ(ω̃)/h(ω̃)^2),
        γ_0=(μ_1(μ_1-1)^2/ω h(ω̃)-1)ω^2,

where ω̃=ω/β_0 and h(ω̃)=μ_1ω̃+S_1 S_1ω̃.



Noting that
L_0(ω) =⟨𝐯⟩, 
where  

    𝐯=([ 1μ_1-1(μ_1S_1^-1tanh S_1ω̃ω̃+1);          -μ_1S_1^-1tanh S_1ω̃ω̃;                               1;                                 ]),

we find that  

    L_0
    =
    ⟨{𝐯sin(𝐤·𝐱),𝐯cos(𝐤·𝐱): 𝐤∈ℒ^⋆ with |𝐤|=ω}⟩

(see Figure <ref>).



  * For rolls the dual lattice ℒ^⋆ is generated by 𝐤=(ω,0), so that |𝐤|,|-𝐤|=ω 
and hence L_0=2.

  * For rectangles the dual lattice ℒ^⋆ is generated by 𝐤_1=(ω,0) and 𝐤_2=(0,ω), 
so that |𝐤_1|,|-𝐤_1|,|𝐤_2|,|-𝐤_2|=ω 
and hence L_0=4.

  * For hexagons the dual lattice ℒ^⋆ 
is generated by 𝐤_1=(ω,0)  and 𝐤_2=ω(1/2,√(3)/2),
so that |𝐤_1|,|-𝐤_1|,|𝐤_2|,|-𝐤_2|, |𝐤_3|,|-𝐤_3|=ω, 
where 𝐤_3=𝐤_2-𝐤_1, 
and hence L_0=6.







Define the projection P_ω:X_0→ X_0, Y_0→ Y_0 
by 

    P_ω
    v(𝐱)
       =
    ∑_𝐤∈ℒ^⋆, 
     |𝐤|=ω𝐯__𝐤e^i𝐤·𝐱,

where v is given by formula (<ref>),
so that

    X_0=X_ω⊕ X_r,
        
    Y_0=Y_ω⊕ Y_r

with

    X_ω=P_ω[X_0],
        
     X_r=(I-P_ω)[X_0],
        
     Y_ω=P_ω[Y_0],
        
     Y_r=(I-P_ω)[Y_0].

Using (<ref>), 
one finds that 
L_0|_ X_ω⊆ Y_ω 
and 
L_0|_ X_r⊆ Y_r  
and we prove 
that L_0:X_0 → Y_0 is a Fredholm operator of index zero in two  steps.

The mapping L_0 
is an isomorphism 
X_r→ Y_r.


The mapping L_0|_ X_r is formally invertible on Y_r with

    L_0^-1(χ,Ψ^',Ψ)^T   =
    ∑_k∈ℒ^⋆ ,
    |𝐤|≠ω
    L_0(|𝐤|)^-1
    (χ__𝐤,Ψ^'__𝐤,Ψ__𝐤)^Te^i𝐤·𝐱,

where 

    L_0(|𝐤|)^-1   =
    μ_1S_1^-1|𝐤|tanhS_1|𝐤|β_0/ L_0(|𝐤|)
    [ (μ_1-1) |𝐤|tanh|𝐤|β_0                     1                     1;         -(γ_0+ |𝐤|^2)              -(μ_1-1)             - (μ_1-1);                     0                     0                     0 ]
    
           +
    
    |𝐤|tanh|𝐤|β_0/ L_0(|𝐤|)
    [          0        μ_1          1;          0          0          0; γ_0 +|𝐤|^2 μ_1(μ_1-1)      μ_1-1 ]
    -γ_0+|𝐤|^2/ L_0(|𝐤|)[   0 0 0  ;   0 1 0  ;   0 1 0   ]

for 𝐤∈ℒ^⋆ with |𝐤|>ω 
and 

    L_0(0)^-1[ χ__0;    0; Ψ__0 ]
    =
    [              -γ_0^-1Ψ__0;                        0; -χ__0
    -(μ_1-1)γ_0^-1Ψ__0 ]


Denoting the right-hand side of equation (<ref>) by (η,Φ^',Φ)^T
and using the estimates

    S_1^-1 |𝐤|tanhS_1 |𝐤|β_0≲ |𝐤|tanh|𝐤|β_0
    ,    |𝐤|tanh|𝐤|β_0/| L_0(|𝐤|)|≲
    |𝐤|^-2

for |𝐤|> ω,
one finds that 

    ‖η‖_ s+1/2^2   ≲(|χ__0|+|Ψ__0|)^2
    +
    ∑_k∈ℒ^⋆, 
    
    |𝐤|≠ω,0
    |𝐤|^2s+1(|𝐤|tanh|𝐤|β_0)^2/( L_0(|𝐤|))^2( 
    |𝐤|tanh|𝐤|β_0|χ__𝐤|
    +|Ψ^'__𝐤|
    +|Ψ__𝐤|
     )^2
       ≲
    |χ__0|^2+|Ψ__0|^2
    +
    ∑_k∈ℒ^⋆, 
    
    |𝐤|≠ω,0
    |𝐤|^2s-3( 
    |𝐤|^2
    |χ__𝐤|^2
    +|Ψ^'__𝐤|^2
    +|Ψ__𝐤|^2)
    
       ≲‖χ‖_ s^2
    +‖Ψ^'‖_ s-1^2
    +‖Ψ‖_ s-3/2^2;

similar calculations yield

    ‖Φ^'‖_ s^2≲‖χ‖_ s^2
    +‖Ψ^'‖_ s-1^2
    +‖Ψ‖_ s-3/2^2
    ,
        ‖Φ‖_ s^2≲‖χ‖_ s^2
    +‖Ψ^'‖_ s-1^2
    +‖Ψ‖_ s-3/2^2
    .

We conclude that L_0^-1: Y_r→ X_r exists and is continuous.


The operator L_0:X_0 → Y_0 is a Fredholm operator of index zero.


A straightforward calculation shows that 𝐯∉L_0(ω) 
(so that (L_0(ω))^2= L_0(ω)) 
and hence 

    Y_ω   =
    ⊕_𝐤∈ℒ^⋆, 
    
    |𝐤|=ω(( L_0(ω)⊕ L_0(ω))sin(𝐤·𝐱)
    ⊕(L_0(ω)⊕ L_0(ω))cos•·𝐱))
    
       
    =
    L_0|_ X_ω⊕ L_0.

Using this decomposition and Lemma <ref>, we find that 

    Y_0 = Y_ω⊕ Y_r = (L_0|_ X_ω⊕ L_0) ⊕ L_0|_U_r
    =
     L_0⊕ L_0.

It follows that 
L_0 is closed and 
 codim  L_0= L_0 , 
 so that L_0:X_0 → Y_0 is a Fredholm operator of index zero. 
   


Because the kernel of L_0 is  multidimensional, we can not use Theorem <ref> directly. 
To overcome this problem, we recall that 𝒢 (and hence L_0)
is invariant under certain rotations (see below)
and seek solutions to (<ref>)  in X_0 that have this rotational symmetry,
denoting the relevant  subspaces of H^r_per(Γ) ,
H̅^r_per(Γ),X_0 and Y_0  
by H^r_sym(Γ),
H̅^r_sym(Γ), X_sym and Y_sym , 
so that 

    X_sym=
    H^s+1/2_sym(Γ)
    ×H̅^s_sym(Γ)
    ×
    H^s_sym(Γ),
        
    Y_sym=
    H^s_sym(Γ)
    ×H̅^s-1_sym(Γ)
    ×
    H^s-3/2_sym(Γ)

for s>5/2. Note that 
X_sym and Y_sym are invariant under P_ω, so that
according to the above analysis L_0:X_sym→ Y_sym is a Fredholm operator of index zero and

    Y_sym= L_0 ⊕ L_0.




  * 

For rolls we consider functions that are independent of the z-coordinate and lie 
in the subspace

    H^r_sym(Γ)
    =
    {ζ∈ H^r_per(Γ):ζ( x )  =  ζ( -x )
     }

of functions which are invariant under rotations through 
π. 

  * For rectangles 
we work with the subspace

    H^r_sym(Γ)
    =
    {ζ∈ H^r_per(Γ):ζ( x,z )  =  ζ( z,-x )
     }

of functions which are invariant under rotations through 
π/2.  

  * For hexagons 
we work with the subspace

    H^r_sym(Γ)
    =
    {ζ∈ H^r_per(Γ):ζ( x,z )=
    ζ( 1/2( x-√(3)z ),
    1/2( √(3)x+z ) )
     }

of functions which are invariant under rotations through 
π/3.  

These restrictions ensure that L_0=1 
with L_0=⟨  v_0⟩, 
where  v_0=𝐯 e_1(x,z) 
with

    e_1(x,z) = {[                                   cosω x                                         ;                            cosω x+cosω z                                         ; cosω x
    +cosω/2(x+√(3)z)
    +cosω/2(x-√(3)z)                                          ].

The projection 
P:Y_sym→ Y_sym onto L_0 along L_0 is given by 

    P(χ, Ψ^', Ψ)^T
    =
    C^⋆((χ__(ω,0), Ψ^'__(ω,0), Ψ__(ω,0))^T·𝐯^⋆)𝐯  e_1(x,z),

where


    𝐯^⋆=([                       γ_0+ω^2μ_1-1; γ_0+ω^2(μ_1-1)^2ω(S_1^-1 S_1ω̃+ω̃);                                  1;                                    ]),


    C^⋆=
    μ_1ω
    S_1^-1tanh S_1ω̃
    -
    μ_1^2
    S_1^-1tanh S_1 ω̃(ω̃
    + S_1 S_1 ω̃)	
    tanhω̃(
    μ_1ω̃
    + S_1 S_1 ω̃)
    +1

(𝐯^⋆∈ℝ^3  solves the equation L_0(ω)^T𝐯^⋆=0
and C^⋆=(𝐯·𝐯^⋆)^-1). 

The transversality condition P(d_1d_2𝒢[γ_0,0](1,v_0))≠ 0 is satisfied.


It follows from the calculation
d_1d_2𝒢[γ_0,0](1,(η,Φ^',Φ)^T)
=(0,0,-η)^T and the formula (<ref>) for P
that

    P(d_1d_2𝒢[γ_0,0](1,v_0))
    =-C^⋆(μ_1-1)^-1S_1^-1tanh S_1ω̃(μ_1ω̃+S_1 S_1ω̃)𝐯  e_1(x,z).


The facts established 
above 
confirm that the hypotheses of 
Theorem <ref> 
are satisfied, 
an application of which yields Theorem <ref>.



§ THE BIFURCATING SOLUTION BRANCHES
 

In this section we examine the bifurcating solution branches 
identified in Theorem <ref> by applying the following
supplement to the Crandall-Rabinowitz theorem.



Suppose that the hypotheses of Theorem <ref> hold. In the
notation of that theorem, 
let Q:X→  X be a projection with Q =  L 
and 
the Taylor series of 
the functions w:(-ε,ε)→ V, λ: (-ε,ε)→ℝ be given by 

    λ(s) = λ_0+s λ_1+s^2λ_2+…,     
    w(s) =  v_0+s w_1+…,

where λ_1,λ_2,…∈ℝ and 
w_1,w_2,…∈ Q. 


  * The coefficient λ_1 satisfies the equation 

    P(1/2!d_2^2ℱ[λ_0,0](v_0,v_0))
    +λ_1P(d_1d_2ℱ[λ_0,0](1,v_0))
    =0

and the bifurcation is transcritical if λ_1 is non-zero (see Figure <ref>). 



  * Suppose that λ_1 is zero. The coefficient λ_2 satisfies the equation 

    P( d_2^2ℱ[λ_0,0](v_0,w_1)
    +1/3!d_2^3ℱ[λ_0,0](v_0,v_0,v_0) )+
    λ_2P( d_1d_2ℱ[λ_0,0](1,v_0) )
    =0,

where w_1∈ Q solves the equation 

    d_2ℱ[λ_0,0](w_1)
    =
    -1/2!d_2^2ℱ[λ_0,0](v_0,v_0).

The bifurcation is supercritical for λ_2>0 
and  subcritical  for λ_2<0 (see Figure <ref>). 






To apply this theorem we write the Taylor series of the analytic functions
w:(-ε,ε)→ V_sym,
γ: (-ε,ε)→ℝ given in
Theorem <ref> as 

    γ(s) = γ_0+s γ_1+s^2γ_2+…,     
    w(s) =  v_0+sw_1+…

with γ_1,γ_2,…∈ℝ and w_1,w_2,…∈ (I-P) 
and introduce the operators  

    L_1=d_1d_2𝒢[γ_0,0](1,·),
        
    Q_0=1/2!d_2^2𝒢[γ_0,0],
        
    C_0=1/3!d_2^3𝒢[γ_0,0].

One finds that L_1(v)=(0,0,-η)^T and

    Q_0   (v,v)
    	
       =
    	[                                                                                           0;                                       G_2^'  
    	+ G_2 +μ̇_1 H_2 
    	+12(μ̈_1 H_1^2+μ̇_1|∇Φ|^2); 12(
    	G^'  2_1-|∇Φ^'|^2
    	+(μ_1-μ̇_1) H_1^2
    	+μ_1|∇Φ|^2
    	-2 G_1 H_1 
    	-μ̈_1 H_1^2-μ̇_1|∇Φ|^2);                                                                -μ_1 G_2^'  
    	-G_2 -μ̇_1 H_2 ]


    C_0   (v,v,v)
    	
        =	[                                                                                                           0;       G_3^' 
    	+ 
    	G_3 +μ̇_1 H_3 
    	+μ̈_1 H_1 H_2 
    	+12(μ̈_1-μ̇_1) |∇Φ|^2H_1 
    	-μ̇_1∇η·∇Φ H_1 
    	+16⃛μ_1 H_1^3; -
    	μ_1 G_3^' 
    	-
    	G_3 -μ̇_1 H_3 
    	-μ̈_1 H_1 H_2 
    	-12(μ̈_1-μ̇_1) |∇Φ|^2H_1 
    	+μ̇_1∇η·∇Φ H_1 
    	-16⃛μ_1 H_1^3;             +
    G_1^' (
    G_2^' -∇η·∇Φ^'
    )
    -H_1 (G_2 +μ_1∇η·∇Φ)
    	+(μ_1-μ̇_1)H_1 H_2 
    	+
    	13
    	(μ̇_1-μ̈_1)
    	H_1^3;                                                    -G_1 H_2+
    	η_x^2η_zz+η_z^2η_xx
    	-2η_xη_zη_xz
    	-32|∇η|^2η ]
    ,

where  v=(η, Φ^', Φ)^T and ⃛μ_1=⃛μ(1).
Theorem <ref> shows that

    γ_1
    =
    -[Q_0(v_0,v_0)]_1·𝐯^⋆/[L_1v_0]_1·𝐯^⋆,

where

    [ζ]_1=ζ__(ω,0)=∫_Γζ e_1

(with componentwise extension),
and

    γ_2   =
    -[2Q_0(v_0,w_1)
    +C_0(v_0,v_0,v_0)]_1·𝐯^⋆/[L_1v_0]_1·𝐯^⋆

for γ_1=0, where 
w_1∈ (I-P) solves the equation

    L_0w_1   =
    -Q_0(v_0,v_0).


A straightforward calculation shows Q_0(v_0, v_0) can be written as a sum
in which each summand is a constant vector multiplied by either e_1^2
or |∇ e_1|^2.
For hexagons we find that γ_1 generally does not vanish,
while 
for rolls

    [e_1^2]_1
    = 
    [1/2
    +1/2cos 2 ω x]_1
    =0
    ,
        [|∇ e_1|^2]_1
    = 
    [1/2
    -1/2cos 2 ω x]_1
    =0,

and for rectangles

    [e_1^2]_1   = 
    [1/2
    +
    cosω (x+z)+cosω (x-z)
    +1/2(cos 2 ω x+cos 2ω z)]_1
    =0,
    
    [|∇ e_1|^2]_1   = 
    [1
    -1/2(cos 2 ω x  +cos 2ω z)]_1
    =0,

so that  in both cases γ_1=0.
Attempting to compute explicit general expressions for γ_2 leads to unwieldy formulae
(it appears more appropriate to calculate them numerically for a specific choice of μ, that is a
specific magnetisation law). Here we confine ourselves to stating the values of the coefficients for two particular special cases.


  * Constant relative permeability μ (corresponding to a linear magnetisation law):
We find that

    γ_2   
    =
    -μμ-1/μ+1ω^2(-μ C_2((8(μ+1)ω^24ω^2+γ_0/t_2-4(μ^2-1)^2ω^3)(1-t_1t_2)^2
       
    +ω^3(μ-1)^4(
    2(
    1+t_1^2)(1-t_1t_2)
    -1/4(1+t_1^2)^2))
    
       
    -μω^2/4γ_0(μ-1)^3/μ+1(1-t_1^2)^2
    -
    (μ+1)^2ω/μ-1t_1(
    3 ω^2/8(γ_0+ω^2)
    -3/2+t_1t_2))

for rolls and

    γ_2   =
    -μμ-1/μ+1ω^2(
    -μ C_√(2)((8(μ+1)ω^22ω^2+γ_0/t_√(2)-2√(2)(μ^2-1)^2ω^3)(1-√(2)t_1t_√(2))^2
       
    -ω^3(μ-1)^4/2√(2)(t_1^4-4t_1^2(1-√(2)t_1t_√(2))))
    
       -μ C_2((8(μ+1)ω^24ω^2+γ_0/t_2-4(μ^2-1)^2ω^3)(1-t_1t_2)^2
       
    +ω^3(μ-1)^4(
    2(
    1+t_1^2)(1-t_1t_2)
    -1/4(1+t_1^2)^2))
    
       
    -ω (μ+1)^2/2(μ-1)t_1(
    5ω^2/4(γ_0+ω^2)-9+2t_1t_2+4√(2)t_1t_√(2))
       
    -ω^2μ(μ-1)^3/2γ_0(μ+1)(1-t_1^2)^2)
 
for rectangles, where t_1=tanhω̃, t__√(2)=tanh√(2)ω̃,
t_2 = tanh 2ω̃, ω̃=ω/β_0 and

    C__√(2)   =
    1/√(2)(μ^2-1)ω(√(2)(ω^2+γ_0)t_√(2)/t_1-γ_0-2ω^2)^-1, 
    
    C_2
       =
    1/2(μ^2-1)ω( 
    2( γ_0+ω^2 )
    t_2/t_1
    -γ_0
    -4ω^2)^-1

The sign of γ_2 clearly depends upon μ and ω̃
(see Figure <ref>).



  * Small values of β_0 (corresponding to deep fluids): Abbreviating 
μ(1),μ̇(1),μ̈(1),⃛μ(1) to respectively μ_1,μ̇_1,μ̈_1,⃛μ_1, one finds that
		
    γ_2   =
    			(
    						t(
    2  μ_1^3⃛μ_1+2  μ_1^2μ̇_1⃛μ_1-2  μ_1^2⃛μ_1-2 μ_1μ̇_1⃛μ_1-112μ_1^2μ̈_1^2+112μ_1μ̈_1^2+42  μ_1^3μ̈_1+49  μ_1^2μ̇_1μ̈_1-8  μ_1^2μ̇_1
               
    			+29  μ_1μ̇_1^2μ̈_1
    			-10  μ_1^2μ̈_1+15 μ_1μ̇_1μ̈_1+3 μ̇_1^2μ̈_1+8  μ̇_1μ_1^3
    +16  μ_1^2μ̇_1^2+8  μ_1μ̇_1^3-16  μ_1μ̇_1^2-8 μ̇_1^3)
    						
           
    +
    16  μ_1^4μ̈_1+32  μ_1^3μ̇_1μ̈_1+16  μ_1^2μ̇_1^2μ̈_1
    +30  μ_1^2μ̇_1μ̈_1+5  μ_1^2μ̈_1^2+10 μ_1μ̇_1^2μ̈_1
    -16  μ_1^2μ̈_1-46 μ_1μ̇_1μ̈_1
           
    -5 μ_1μ̈_1^2
    -10 μ̇_1^2μ̈_1)(μ_1-1)^7μ_1^2/512t^6(μ_1+μ̇_1)
         
    			+
    			(
    			t(
    96  μ_1^6+636  μ_1^5μ̇_1+1129  μ_1^4μ̇_1^2
    +850  μ_1^3μ̇_1^3+217  μ_1^2μ̇_1^4-512  μ_1^5-1720  μ_1^4μ̇_1-1810  μ_1^3μ̇_1^2
                 
    -740  μ_1^2μ̇_1^3-50  μ_1μ̇_1^4+96  μ_1^4+220  μ̇_1μ_1^3+73  μ_1^2μ̇_1^2
    -14 μ_1μ̇_1^3-7μ̇_1^4)	
             					
    +88  μ_1^8+440  μ_1^7μ̇_1
    +880  μ_1^6μ̇_1^2+880  μ_1^5μ̇_1^3
    +440  μ_1^4μ̇_1^4+88  μ_1^3μ̇_1^5-256μ_1^7-928  μ_1^6μ̇_1-1312  μ_1^5μ̇_1^2
             -864  μ_1^4μ̇_1^3
    			-224  μ_1^3μ̇_1^4-80  μ_1^6-416  μ_1^5μ̇_1-614  μ_1^4μ̇_1^2
    -420  μ_1^3μ̇_1^3-102  μ_1^2μ̇_1^4+256 μ_1^5+672  μ_1^4μ̇_1
             
    -8  μ_1^4+460  μ_1^3μ̇_1^2+104 μ_1^2μ̇_1^3-20  μ_1μ̇_1^4+72  μ̇_1μ_1^3+194  μ_1^2μ̇_1^2+84  μ_1μ̇_1^3+10 μ̇_1^4)(μ_1-1)^6μ_1/1024t^6(μ_1+μ̇_1)
          + o(1)

		
as β_0 → 0 for rolls and

    γ_2   =
    			(
    					t	(
    (42 √(2)-60)μ_1⃛μ_1μ̇_1-(42 √(2)-60)μ_1^2⃛μ_1μ̇_1    
    +(411√(2)-289) μ_1μ̈_1^2
    +(4461 √(2)-6190)  μ̈_1μ_1^2μ̇_1
               
    +(3413  √(2)-4766 ) μ̈_1μ_1μ̇_1  ^2-(429√(2)-558)  μ̈_1μ_1μ̇_1
    +(168 √(2)-240)   μ̇_1  ^3
    +(168 √(2)-240) μ_1^2μ̇_1
               
    +(336 √(2)-480 )μ_1μ̇_1  ^2
    +(1870  √(2)-2580  )μ̈_1μ_1^3
    			+(82 √(2)-172  )μ̈_1μ_1^2-(1333 √(2)-1886)  μ̈_1μ̇_1  ^2
               
    +(42 √(2)-60) μ_1^2⃛μ_1-(411 √(2)-289) μ_1^2μ̈_1 ^2-(168 √(2)-240)μ_1μ̇_1^3-(168  √(2)-240)   μ̇_1μ_1^3
               -(336√(2)-480) μ_1^2μ̇_1  ^2-(42 √(2)-60) μ_1^3⃛μ_1)	
    			
           
    			+(976 √(2)-1376) μ_1^4μ̈_1
    			+(1952  √(2)-2752)   μ̇_1μ_1^3μ̈_1
    			+(976 √(2)-1376) μ_1^2μ̇_1  ^2μ̈_1
    			+(445 √(2)-622) μ_1^2μ̈_1 ^2
           
    			+(366 √(2)-468 )μ_1^2μ̇_1μ̈_1-(1414 √(2)-2020) μ_1μ̇_1  ^2μ̈_1-(445√(2)-622) μ_1μ̈_1 ^2-(976 √(2)-1376) μ_1^2μ̈_1
           -(1342 √(2)-1844) μ_1μ̇_1μ̈_1
    			+(1414 √(2)-2020)   μ̇_1  ^2μ̈_1)(μ_1-1)^7μ_1^2/512(10-7√(2))t^6(μ_1+μ̇_1)
         
    			+
    			(
    			t (
    			(3300√(2)-3288 )  μ̇_1μ_1^3
    			+(5745 √(2)-7110 )μ_1^2μ̇_1  ^2
    			+(2410 √(2)-3260) μ_1μ̇_1^3
                 
    			+(6084 √(2)-6168) μ_1^5μ̇_1 
    			+(7377 √(2)-6406 )μ_1^4μ̇_1  ^2-(3190 √(2)-7300) μ_1^3μ̇_1  ^3
                 
    			-(6799 √(2)-10298) μ_1^2μ̇_1  ^4-(58696 √(2)-79024) μ_1^4μ̇_1-(61154 √(2)-81804) μ_1^3μ̇_1  ^2
                 -(11604 √(2)-14456) μ_1^2μ̇_1  ^3
    			+(11550 √(2)-16500) μ_1μ̇_1  ^4
    			+(672 √(2)-448) μ_1^6
                 -(17408 √(2)-23552) μ_1^5-(2351 √(2)-3322)   μ̇_1  ^4
    			+(672 √(2)-448) μ_1^4)
    			
             
    			+(2344√(2)-3056 )μ_1^3μ̇_1  ^5
    			+(11720 √(2)-15280) μ_1^7μ̇_1  
    			+(23440 √(2)-30560) μ_1^6μ̇_1  ^2
             
    			+(23440 √(2)-30560 )μ_1^5μ̇_1  ^3
    			+(11720 √(2)-15280 )μ_1^4μ̇_1^4-(36384 √(2)-49344) μ_1^6μ̇_1
             
    			-(60832√(2)-82880) μ_1^5μ̇_1  ^2-(47328 √(2)-64832 )μ_1^4μ̇_1^3-(14176 √(2)-19520) μ_1^3μ̇_1  ^4
             
    			+(4344 √(2)-6608)   μ̇_1μ_1^3-(478 √(2)-628)μ_1^2μ̇_1  ^2-(4860 √(2)-7016) μ_1μ̇_1  ^3
    			+(2344 √(2)-3056) μ_1^8
             -(8704 √(2)-11776) μ_1^7-(15392 √(2)-21440) μ_1^5μ̇_1
    			-(18422 √(2)-25348) μ_1^4μ̇_1  ^2
             -(8196 √(2)-10904) μ_1^3μ̇_1  ^3
    			+(2410 √(2)-3644) μ_1^2μ̇_1  ^4
    			+(27680 √(2)-37568) μ_1^4μ̇_1
             
    			+(33644 √(2)-46344) μ_1^3μ̇_1  ^2
    			+(18088 √(2)-25328 )μ_1^2μ̇_1  ^3
    			-(3700 √(2)-5176) μ_1μ̇_1  ^4
             
    			-(4016 √(2)-5664) μ_1^6
    			+(8704 √(2)-11776 )μ_1^5
    			+(1850 √(2)-2588)   μ̇_1  ^4
             
    			+(1672 √(2)-2608 )μ_1^4)(μ_1-1)^6μ_1/1024(10-7√(2))t^6(μ_1+μ̇_1)+o(1)


as β_0 → 0 for rectangles, 	where t=√(μ_1(μ_1+μ̇_1))+1.
(Note that

    ω = μ_1(μ_1-1)^2/2(μ_1+S_1)+o(1),     γ_0 = μ_1^2(μ_1-1)^4/4(μ_1+S_1)^2+o(1)

as β_0 → 0.)

We note in particular that for constant μ (corresponding to a linear magnetisation law),
rolls bifurcate subcritically for μ<μ_c^1 and
supercritically for μ>μ_c^1, while rectangles bifurcate subcritically for
μ<μ_c^2 and supercritically for μ>μ_c^2, where

    μ_c^1=21/11+8/11√(5),
        μ_c^2=115+160√(2)+8√(184+11√(2))/141+128√(2).



Figure <ref> shows the sign of γ_2 for the Langevin magnetisation law

    μ(s)=1+M/s((γ s)-1/γ s)

in the limit β_0 → 0, where M and χ_0
are respectively the magnetic saturation and initial susceptibility of the ferrofluid
and γ=3χ_0/M.






This paper has no additional data.

M.D.G. and J.H. collaborated equally on each aspect of this work.

We declare we have no competing interests.

No external funding was involved in this work.



99

CowleyRosensweig67
Cowley MD, Rosensweig RE. 1967  The interfacial stability of a ferromagnetic
  fluid. J. Fluid Mech. 30, 671–688.

TwomblyThomas83
Twombly EE, Thomas JW. 1983  Bifurcating instability of the free surface of a
  ferrofluid. SIAM J. Math. Anal. 14, 736–766.

CraigNicholls00
Craig W, Nicholls DP. 2000  Traveling two and three dimensional capillary
  gravity water waves. SIAM J. Math. Anal. 32, 323–359.

BuffoniToland
Buffoni B, Toland JF. 2003 Analytic Theory of Global Bifurcation.
Princeton, N. J.: Princeton University Press.

Sattinger78
Sattinger DH. 1978  Group representation theory, bifurcation theory and pattern
  formation. J. Func. Anal. 28, 58–101.

SilberKnobloch88
Silber M, Knobloch E. 1988  Pattern selection in ferrofluids. Physica D
  30, 83–98.

LloydGollwitzerRehbergRichter15
Lloyd DJB, Gollwitzer C, Rehberg I, Richter R. 2015  Homoclinc snaking near the
  surface instability of a polarisable fluid. J. Fluid Mech.
  783, 283–305.

GrovesLloydStylianou17
Groves MD, Lloyd DJB, Stylianou A. 2017  Pattern formation on the free surface
  of a ferrofluid: spatial dynamics and homoclinic bifurcation. Physica D
  350, 1–12.

Rosensweig
Rosensweig RE. 1997 Ferrohydrodynamics.
New York: Dover.

NichollsReitich01a
Nicholls DP, Reitich F. 2001  A new approach to analyticity of
  Dirichlet-Neumann operators. Proc. Roy. Soc. Edin. A
  131, 1411–1433.




