## Abstract

On a weakly nonlinear basis, we revisit the pattern formation problem in the Boussinesq convection, for which nonlinear terms of the quadratic order are known to vanish from amplitude equations. It is thus necessary to proceed to the quintic-order approximation in order for the amplitude equations to be generic. By deriving the quintic amplitude equations from the governing PDEs, we examined the bifurcation of steady solutions under rigid–free, rigid–rigid and free–free boundary conditions. Right above the criticality, all the axial solutions are obtained including up- and down-hexagons under the asymmetric boundary conditions and hexagons and regular triangles under the symmetric conditions. Hexagons and regular triangles are unstable whereas rolls are stable as has already been predicted by the cubic-order amplitude equations. Irrespective of the boundary conditions, quintic-order terms stabilize hexagons except near the criticality; rolls and hexagons thus coexist stably in an open region. This suggests that amplitude equations of higher order are possible to predict *re-entrant hexagons*.

## 1. Introduction

At the onset of the Rayleigh–Bénard (hereafter referred to as RB) convection, a flow structure is determined by a combination of eigenmodes whose wavevectors are on a critical circle. Since the eigenvalues have infinite fold degeneracy on the circle, the wavevectors of the critical eigenmodes are usually restricted on a square or a hexagonal lattice to lighten the degeneracy; the governing hydrodynamic equations are reduced to four- or six-dimensional amplitude equations. In the following, we consider a situation in which six points of the hexagonal lattice are exactly on the critical circle. During the 1960s, amplitude equations were derived up to the cubic order on the hexagonal lattice.1 Under the Boussinesq approximation and with the density that is a linear function of the temperature, it is established that rolls are stable right above the criticality. If the Boussinesq approximation is invalid or the density profile ceases to be linear, on the other hand, the hexagons are preferred to the rolls. For details, see Busse (1962, 1967, 1978) for example.

Since the late 1970s, the equivariant bifurcation theory has been applied to pattern formation problems. Under the Boussinesq approximation with a linear density profile and with an up–down reflectional symmetry in the horizontal mid-plane, it is clarified that amplitude equations of the cubic order are not generic; since they cannot give all the branches of axial solutions,2 we cannot distinguish hexagons from regular triangles at the cubic-order approximation (Golubitsky *et al*. 1984, 1988). Nonetheless, stability characteristics of these axial solutions are not altered right above the criticality even if the amplitude equations are truncated at the cubic order. At the lowest order, the signs of the eigenvalues of Jacobian matrices for hexagons and regular triangles are determined from(1.1)where the minus sign is for hexagons and the plus sign for regular triangles. (For , *l*_{3}(0) and *m*_{5}(0), see §§2 and 3.) The *m*_{5}(0) arises at the quintic order and is thus missing at the cubic order. But for the RB convection, *m*_{5}(0) is not responsible for the stability of hexagons and regular triangles since and *l*_{3}(0) are known to have different signs. At least near the criticality, therefore, without evaluating these solutions, they are both concluded to be unstable if they exist.

From the bifurcation point of view, however, at least all the axial solutions need to be distinguished from each other near the criticality. Moreover, the eigenvalue proportional to ∓*m*_{5}(0) tells us which solutions are less unstable between hexagons and regular triangles. If the initial conditions happen to be close to the less unstable solution, it can be realized at least for a certain period of time. To obtain the branches of hexagons and regular triangles and to evaluate *m*_{5}(0), nonlinear terms of at least quintic order are needed in the amplitude equations.

A similar situation arises for hexagons under rigid–free boundary conditions. Since the mid-plane reflection is absent, equilibrium solutions of generic amplitude equations involve rolls and up- and down-hexagons as the axial solutions. Under the Boussinesq approximation, the coefficient of quadratic nonlinear terms vanishes from the amplitude equations, as has been pointed out by Schlüter *et al*. (1965). Up- and down-hexagons can then be distinguished not at the cubic order, but at the quartic order. The signs of the eigenvalues for them are given by(1.2)where the minus sign is for up-hexagons and the plus sign is for down-hexagons. (For , *h*_{3}(0), and *p*_{4}(0), see §§2 and 3.) Again, for the RB convection, the stability of up- and down-hexagons is determined at the cubic order at least near the criticality (see §4). At the quartic order, however, we can find which are less unstable between up- and down-hexagons since arises at the quartic order.

To our knowledge, amplitude equations with quintic-order terms have not yet been derived from the governing PDEs. The only exception was due to Knobloch (1990). He investigated planforms exhibited by long-wave instability modes; in the RB convection between nearly thermally insulating boundaries or the Bénard–Marangoni convection, for example, long-wave modes give the critical condition. The nonlinear PDEs for a three-dimensional convection are then shown to be reduced to a nonlinear PDE on the horizontal plane; it describes the spatio-temporal evolution of a planform function *f*(*x*, *y*, *t*) and has the formOnce this PDE is derived, it is not difficult to derive the amplitude equations of the quintic order. At the cubic order, in the long-wave problem, a degeneracy *l*_{3}(0)=0 arises. This does not occur in the RB convection between perfectly conducting boundaries with a finite critical wavenumber.

In the present paper, we revisit this very classical problem with a finite critical wavenumber. Our objective is to derive the generic amplitude equations of the quintic order and distinguish all the branches of axial solutions including hexagons and regular triangles under the mid-plane reflection and up- and down-hexagons under rigid–free conditions right above the criticality. In §2, we recall that cubic-order equations are not enough to distinguish all the axial solutions. Section 3 is devoted to explaining how to derive the amplitude equations on a hexagonal lattice. Numerical results are presented in §4. Coefficients of the amplitude equations are evaluated numerically for *P*=0.025, 0.71, 7 and 1000 under rigid–free, rigid–rigid and free–free boundary conditions. Based on the coefficients, we distinguish all the axial solution branches right above the criticality, although they are unstable except for rolls. As a by-product, we show that hexagons are stabilized far above the criticality. Some concluding remarks are given in §5, especially on a relationship between stable hexagons and *re-entrant hexagons*.

## 2. Steady solutions of the amplitude equations

In this section, in order for this paper to be self-contained, we quickly review typical steady solutions of the amplitude equations on a hexagonal lattice. See Buzano & Golubitsky (1983), Golubitsky *et al*. (1984), or Golubitsky *et al*. (1988) for details. In the tables of appendix A, we list the eigenvalues of the Jacobian matrices for those steady solutions.

### (a) Asymmetric case

A hexagonal lattice has the group of symmetry . We consider a situation in which six lattice points are exactly on the critical circle in the wavenumber plane. A disturbance ** ψ**=(

**,**

*v**θ*) added to the conduction state is written aswhere

**=(**

*v**u*,

*v*, ,

*w*) and

*θ*denote the velocity and the temperature components, respectively, and the dots denote higher order terms. The group

*Γ*

_{A}acts on

**=(**

*ζ**ζ*

_{1},

*ζ*

_{2},

*ζ*

_{3})∈

^{3}as and for 0≤

*s*,

*t*<2

*π*. The equation for

*ζ*

_{1}is then generated by the

*Γ*

_{A}-equivariant vector field

*g*

_{1}(

**) as(2.1)where**

*ζ**u*

_{j}=|

*ζ*

_{j}|

^{2}. The

*h*

_{1},

*h*

_{3},

*h*

_{5},

*p*

_{2},

*p*

_{4}and

*p*

_{6}are real-valued functions of four

*Γ*

_{A}-invariants

*σ*

_{1}=

*u*

_{1}+

*u*

_{2}+

*u*

_{3},

*σ*

_{2}=

*u*

_{1}

*u*

_{2}+

*u*

_{2}

*u*

_{3}+

*u*

_{3}

*u*

_{1},

*σ*

_{3}=

*u*

_{1}

*u*

_{2}

*u*

_{3}and and a bifurcation parameter

*δ*(see Buzano & Golubitsky 1983).

As the steady primary solutions, the branches of rolls (R) and hexagons (H) bifurcate from the conduction state (I). The branches of triangles (T) and rectangles (RA) bifurcate as the secondary branches. Solutions exhibiting rolls with the isotropy subgroup are given by ** ζ**=(

*ξ*, 0, 0),

*ξ*∈. Solutions exhibiting rectangles with the isotropy subgroup are given by

**=(**

*ζ**ξ*,

*η*,

*η*),

*ξ*,

*η*∈.

Let us concentrate on solutions satisfying *ζ*_{1}=*ζ*_{2}=*ζ*_{3}; the solutions exhibit hexagons or triangles. Setting (*j*=1, 2, 3) and *ϑ*_{1}+*ϑ*_{2}+*ϑ*_{3}=*Θ*(*t*), we write the branching equations as(2.2)Solutions having the isotropy subgroup D_{6} exhibit hexagons. They satisfy cos *Θ*=±1 and are given by ** ζ**=(

*ξ*,

*ξ*,

*ξ*),

*ξ*∈. The plus sign indicates up-hexagons (or

*l*-hexagons) while the minus sign indicates down-hexagons (or

*g*-hexagons). We denote these solutions by H

_{+}and H

_{−}, respectively. The

*r*for H

_{±}is governed by(2.3)Schlüter

*et al*. (1965) point out that

*p*

_{2}(0) vanishes if linear operators in the nonlinear governing PDEs are all self-adjoint. In this case, at the lowest order, a term proportional to

*q*is needed in

*h*to extract the difference between cos

*Θ*=1 and −1. We thus need to proceed to the quartic order, at least, in order for the Taylor expansion of

*h*to involve the cubic-order nonlinear term proportional to cos

*Θ*.

Solutions having the isotropy subgroup D_{3} exhibit triangles. They are given by cos *Θ*≠±1, i.e. ** ζ**=(

*ζ*,

*ζ*,

*ζ*),

*ζ*∈. From the equations(2.4)it is obvious that, at least, quintic-order terms are needed in the amplitude equations to specify

*r*and

*Θ*(≠

*nπ*) uniquely.

In summary, to distinguish H_{+} from H_{−}, we need to derive the amplitude equations of the quartic order, whereas to resolve the solutions T bridging between H_{+} and H_{−}, we need to derive the amplitude equations of the quintic order, at least.

### (b) Symmetric case

Suppose that a horizontal fluid layer is confined between *z*=0 and *d*. In the presence of an up–down reflectional symmetry about the horizontal mid-plane *z*=*d*/2, we require *z*→*d*−*z* : *w*→−*w*, *θ*→−*θ*. This symmetry causes an additional group Z_{2} that acts on ** ζ**=(

*ζ*

_{1},

*ζ*

_{2},

*ζ*

_{3}),

*ζ*∈

^{3}as

*σ*

_{h}:

**→(−**

*ζ**ζ*

_{1}, −

*ζ*

_{2},−

*ζ*

_{3}). The equation for

*ζ*

_{1}generated by the -equivariant vector field is given by(2.5)where

*l*

_{1},

*l*

_{3},

*l*

_{5},

*m*

_{5},

*m*

_{7}and

*m*

_{9}are real-valued functions of

*Γ*

_{S}-invariants

*σ*

_{1},

*σ*

_{2},

*σ*

_{3}and

*q*

^{2}and a bifurcation parameter

*δ*. See Golubitsky

*et al*. (1984) or (1988) for further details.

In addition to the branches of rolls (R) and hexagons (H), the branches of solutions exhibiting patchwork quilts (PQ) and regular triangles (RT) bifurcate as the primary solution branches. Four secondary solutions possibly to exist, i.e. triangles (T), rectangles (RA), imaginary rectangles and bimodals. Solutions PQ having the isotropy subgroup are given by ** ζ**=(0,

*ξ*,

*ξ*),

*ξ*∈.

Let us now focus on solutions satisfying *ζ*_{1}=*ζ*_{2}=*ζ*_{3}. Again, setting and *ϑ*_{1}+*ϑ*_{2}+*ϑ*_{3}=*Θ*(*t*), we write the branching equations as(2.6)Solutions exhibiting hexagons with the isotropy subgroup D_{6} are given by cos *Θ*=±1, i.e. ** ζ**=(

*ξ*,

*ξ*,

*ξ*),

*ξ*∈. Since the arguments of

*l*and

*m*involve

*q*always in the form

*q*

^{2}, the

*r*for solutions H satisfies(2.7)Therefore, the amplitude equation (2.5) generated by the

*Γ*

_{S}-equivariant vector field cannot distinguish between cos

*Θ*=−1 and 1.

The isotropy subgroup of solutions RT is D_{3}⊕Z_{2}. RT are given by cos *Θ*=0, i.e. ** ζ**=(i

*ξ*, i

*ξ*, i

*ξ*),

*ξ*∈. The

*r*for RT is governed by(2.8)Solutions H are distinguished from those of RT by the existence of the second term in (2.7). This is at the quintic order in (2.5).

Solutions exhibiting triangles have the isotropy subgroup D_{3}. They are given by cos *Θ*≠±1,0, i.e. ** ζ**=(

*ζ*,

*ζ*,

*ζ*),

*ζ*∈. The

*r*and

*Θ*for T are governed by(2.9)Therefore, to specify

*r*and

*Θ*(≠

*nπ*/2), we need to proceed to the seventh-order approximation in the expansion of (2.5).

As mentioned previously, to distinguish four primary solution branches, we need to involve the terms of the quintic order in the amplitude equations.

## 3. Centre manifold reduction

Let us denote the horizontal coordinates by *x*^{*} and *y*^{*}, and the vertical one opposite to the gravity by *z*^{*}, where an asterisk denotes a dimensional quantity. We consider a horizontal fluid layer with an infinite extent that is confined between two boundaries located at *z*^{*}=0 and *d*. Temperatures on the top and bottom boundaries are kept to be constant at *T*^{*}=*T*_{t} and *T*_{b}(>*T*_{t}), respectively. Under the Boussinesq approximation, the velocity *v*^{*}, the pressure *π*^{*} and the temperature *T*^{*} are governed by(3.1)Here, as usual, *ρ* denotes the density; *g* is the acceleration due to the gravity; *μ* is the viscous coefficient; and *κ* is the thermal diffusivity. The *ρ*_{0} is the density at a reference temperature *T*_{0}. In §4*c*, we perturb the self-adjointness of the linear operators. For this purpose, we assume the density of the fluid to be a weak quadratic function of the temperature,(3.2)After an appropriate non-dimensionalization, the Rayleigh number *R*, the Prandtl number *P* and a parameter *ϵ* arise,(3.3)Here, *ϵ* measures a deviation of *ρ* from the linear profile. In this paper, except for §4*c*, we consider a situation where *ϵ*=0 holds. We regard *δ*=(*R*−*R*_{c})/*R*_{c} as the bifurcation parameter, where *R*_{c} denotes the critical value of *R*. Since it is not easy to change physical properties continuously during laboratory experiments, we prescribe the Prandtl number *P*. We also fix the wavenumber at the critical value *k*_{c}; even though a thermal imprinting technique is introduced, it is difficult to control the excited wavenumber continuously during experiments.

First, we split the flow field into the heat conduction state and the disturbance. By (** v**,

*θ*,

*π*), we denote disturbance components of the velocity, temperature and the pressure. We then eliminate the pressure

*π*from a non-dimensional version of the nonlinear PDEs (3.1) and set

**=(**

*ψ***,**

*v**θ*), for brevity. Governing equations and homogeneous boundary conditions for

**are compactly written as(3.4)where , and denote linear operators and denotes quadratic nonlinear terms. We impose three kinds of boundary conditions: two rigid boundaries at**

*ψ**z*=0, 1; a rigid boundary at the bottom

*z*=0 and a free boundary at the top

*z*=1; and two free boundaries at

*z*=0, 1. For brevity, we denote them by

*=0. We note that and are self-adjoint. For vanishing*

**ψ***ϵ*, is self-adjoint. Otherwise, is non-self-adjoint.

We triply expand ** ψ** in terms of the double Fourier series and linear eigenfunctions,(3.5)There is no need to assume

*ϵ*=0 in the course of our derivation of amplitude equations. Using the biorthogonality relationship, we obtain infinite-dimensional ODEs for ,(3.6)The and span the centre manifold (see Carr 1981)(3.7)At the quintic-order approximation, a straightforward but lengthy manipulation yields the equation for

*ζ*

_{1}in the form(3.8)where notations have been changed for simplicity as , , , , and .3 For some details of the centre manifold reduction, see Fujimura (1997).

For *ϵ*=0, *λ*=0 holds, as had been pointed out by Schlüter *et al*. (1965). In the absence of the up–down reflectional symmetry, amplitude equations are generated by (2.1). The Taylor expansion of the functions *h*_{j} and *p*_{j} about the origin and a truncation of the resultant equation yield, up to the quintic order,(3.9)Coefficients in (3.9) are correlated with those in (3.8) byFor *ϵ*=0, *p*_{2}(0) has to vanish. An *O*(*ϵ*) perturbation yields a non-vanishing *p*_{2}(0), which can be written as *ϵp*_{2,ϵ}(0).

In the presence of the reflectional mid-plane symmetry Z_{2}, *λ*, *ν*_{1}, *ν*_{2} and *ν*_{3} vanish. Equations (2.1) and (2.5) imply that if the symmetry Z_{2} is slightly broken by an *O*(*ϵ*) perturbation, the equation for *ζ*_{1} that is generated by the *Γ*_{A}-equivariant vector field is written as(3.10)The Taylor expansion of the functions *l*_{j}, *m*_{j}, and about the origin and a truncation of the resultant equation yield, up to the quintic order,(3.11)Here, we retained only the most essential contributions out of the *O*(*ϵ*) perturbation terms. Coefficients in (3.11) are correlated with those in (3.8) by

## 4. Numerical results based on generic equations

This section is divided into three subsections. In §4*a*, we show numerical values of all the coefficients involved in the amplitude equations (3.9) and (3.11). Section 4*b* is devoted to describing our numerical results for water (*P*=7) with a linear density profile (*ϵ*=0). Section 4*c* describes an effect of the *O*(*ϵ*) perturbation on the bifurcation for water. To unfold the degeneracy in (2.5), we perturb the mid-plane symmetry by assuming the weak quadratic density profile (3.2) that is known to be relevant for water. Recall that the density of water has the maximum at *T*=4°C. It is thus well approximated by the quadratic profile near 4°C. If we assume *T*_{b}−*T*_{t}=1°C, we have *ϵ*≃0.04 at *T*_{0}=10°C. Therefore, *ϵ*≃0.04 is considered to be achievable in laboratory experiments. We note that *ϵ*≥0 for liquid and *ϵ*≤0 for gas. The neutral stability curves for different *ϵ* show that, with a slight increase of *ϵ*, the critical Rayleigh number decreases slightly without a significant change in the critical wavenumber.

### (a) Coefficients in the amplitude equations (3.9) and (3.11)

In the expansions (3.5), 37 Fourier modes with −3≤*m*≤3 and −3≤*n*≤3 are included. All the Fourier coefficients are expanded in 20 eigenfunctions.4 Therefore, our centre manifold reduction is from a 740-dimensional system to a six-dimensional one. To evaluate the coefficients in (3.8), we first need to determine all the coefficients in (3.6). For this purpose, we solve the linear eigenvalue problems and the corresponding adjoint problems numerically by means of expansions in Chebyshev polynomials. The eigenfunctions are normalized, such that, at *z*=1/2, for and for . Numerical integrations necessary in computing and in (3.6) are due to the Gauss–Legendre quadrature. The values of coefficients in (3.9) and (3.11) are listed in tables 1–3, respectively, under rigid–free, rigid–rigid and free–free boundary conditions.

### (b) Bifurcation diagram for *ϵ*=0

In this section, we assume that the density is a linear function of the temperature, i.e. .

Figure 1*a*,*b* shows steady solution branches under rigid–free boundary conditions. Figure 1*a* is due to the cubic truncation of (3.9). This bifurcation diagram is qualitatively the same as that given by the cubic truncation of (3.11) under the mid-plane reflectional symmetry. The solutions RA degenerate into solutions PQ, although the latter never exist as solutions of the generic amplitude equations generated by (2.1) without the mid-plane symmetry. Solutions of rolls (R) are found to be orbitally stable. The stability assignment attached to a label of a branch is such that the plus and minus signs denote the real part of the eigenvalue to be positive and negative, respectively. A multiplicity of an eigenvalue and the zero eigenvalue are not listed. The assignments are in order consistent with the entries in tables 4 and 5 of appendix A. At the cubic order, the solutions of H_{+} and H_{−} are not distinguished from the solutions T. Besides two zero eigenvalues forced by the symmetry, the Jacobian matrices of H_{±} and T have one negative eigenvalue, positive eigenvalues with multiplicity two and one zero eigenvalue corresponding to the absent (see equation (1.1)).

At the quartic order, the solution branch of H_{+} is well separated from the branch of H_{−}. The branch of PQ is now replaced with the branch of RA.5 Figure 1*b* is due to the quintic equation (3.9). The solution branch of T bifurcates from the branch of H_{+}; its bifurcation point *δ*=0.0367 is indicated by the filled circle. The solutions T are unstable. Beyond the bifurcation points of the branches of RA, the solutions of both H_{+} and H_{−} are stabilized.

Let us now examine steady solution branches under the mid-plane symmetry. Figure 2*a* shows the branches under rigid–rigid boundary conditions. The branch of solutions RT bifurcates as the primary solution, but is unstable. Solutions of rolls are always stable. Moreover, for *δ*>0.28, solutions H are stabilized. A qualitatively similar bifurcation diagram is obtained under free–free boundary conditions. Note that the stability boundary of H, *δ*≃0.28, is far above the criticality. If we proceed to the seventh order or higher, we might have a different bifurcation diagram from ours for *δ*≃0.28. In this respect, it is worth noting the work of Nishida *et al*. (2002). They analysed the bifurcation of steady solutions fully numerically under free–free boundary conditions. On a rectangular lattice with the aspect ratio on the *xy*-plane, they obtained a bifurcation diagram that is qualitatively consistent with figure 2*a* except for the stability of solutions PQ and RA. The inconsistencies upon the stability are considered to be due to the functional space they took. Therefore, we expect that, over the range of *δ* we examined, the effect of higher order terms does not change our bifurcation diagrams qualitatively.

### (c) Effect of the *O*(*ϵ*) perturbation

Figure 2*b*–*d* shows the bifurcation diagrams based on (3.11) for *ϵ*=0.002, 0.01 and 0.04, respectively. The solutions RT for *ϵ*=0 are now replaced with the solutions T. Figure 2*e* shows that the branch of T bifurcates from the branch of H_{−}. Figure 2*f* indicates the bifurcation structures that are consistent with those based on the cubic-order equations for *ϵ*=*O*(1) (see Busse (1978) for an example). In figure 2*b*–*f*, the solutions T are found to be unstable. For *ϵ*≠0, both H_{+} and H_{−} are stabilized far above the criticality.

In figure 3, we depict the boundaries of existence regions of steady solutions. Although the boundaries of the regions branching from the origin are not straight, the bifurcation characteristics for all the different but relatively small *ϵ* are qualitatively the same. In figure 3*c*, the solutions of H_{+} are stable between two curves indicated by H_{+}(LP) and H_{+}→RA. The solutions of rolls are stable on the r.h.s of the curve indicated by R→RA. These are consistent with the stability diagram obtained by Busse (1978) for *ϵ*≠0. Figure 3*a* shows that the stability boundaries of up- and down-hexagons are at approximately *δ*≃0.28 and almost independent of *ϵ*.

## 5. Discussion

In this paper, we derived the amplitude equations of the quintic order by restricting the wavevectors of the critical eigenmodes on the hexagonal lattice. The bifurcation of steady solutions was examined under three kinds of boundary conditions. All the axial solutions were distinguished in the neighbourhood of the criticality. Right above the criticality, down-hexagons are less unstable than up-hexagons in the absence of the mid-plane symmetry, whereas hexagons are less unstable than regular triangles under the symmetry. In the latter case, an introduction of the perturbation on the linear density profile breaks the mid-plane symmetry as well as the self-adjointness. The perturbation thus forces unstable regular triangles to be replaced with unstable triangles under the symmetric boundary conditions.

For air (*P*=0.71) and a fluid with *P*=1000, bifurcation structures are qualitatively the same as those for water. For mercury (*P*=0.025) under rigid–rigid boundary conditions, on the other hand, we can conclude that solutions of rolls are stable whereas solutions H and RT are unstable only very close to *δ*=0. The bifurcation diagram for mercury is not reliable except very close to *δ*=0 since the radius of convergence of the amplitude equations is unexpectedly narrow for small Prandtl number fluids. For a two-dimensional RB problem, the radius of convergence has been examined. It was Kuo (1961) who proceeded to the eighth-order approximation in the weakly nonlinear expansion under free–free boundary conditions and pointed out that the parameter expansion based on *δ* yields alternating series and the radius of convergence is narrow. It is reported that, instead of the *δ*=(*R*−*R*_{c})/*R*_{c}, an introduction of a new parameter (*R*−*R*_{c})/*R* (=*δ*_{Kuo}, say) together with suppresses the alternating feature of the series and gives a much wider range of validity. According to his eqn (5.16), this is the case for *P* larger than *O*(1). If the value of *P* decreases far below *O*(1), on the other hand, the alternating feature recovers and the radius of convergence is getting narrower. The introduction of the new parameter *δ*_{Kuo} loses its advantage for small Prandtl number fluids. By proceeding to at most the 25th order, under rigid–rigid boundary conditions, Generalis & Fujimura (2008) numerically examined the range of validity by comparing weakly nonlinear results with fully numerical solutions for rolls based on the Newton–Raphson iteration.6 They concluded that, for *P*<0.25, the r.h.s. of the Stuart–Landau equation is almost alternating and the range of convergence is quite narrow, whereas, for *P*>0.25, the series no longer exhibits the alternating feature. This is consistent with Kuo's result based on the parameter *δ*_{Kuo}=(*R*−*R*_{c})/*R*.

In all the cases for *P*≥0.71, solutions of rolls and hexagons are simultaneously stable in an open region with relatively large *δ*. In figure 4, we depict the stability boundary of the hexagons together with the neutral stability curve, a part of the Busse balloon and the stability boundary of hexagons obtained by Clever & Busse (1996).

Above the neutral curve, the conduction state is unstable with respect to a disturbance that is composed of steady plane waves having the same wavenumbers but with different wavevectors. The planform of such a disturbance is undetermined on the linear basis. Inside the Busse balloon, solutions of rolls are stable. If the parameter crosses the Busse balloon from inside, solutions of rolls are destabilized due to the secondary instability mechanisms.

Clever & Busse (1996) numerically examined the secondary instability of hexagons under the mid-plane reflectional symmetry for Boussinesq fluids. Inside the stability boundary shown by the short/long-dashed line (figure 4), hexagons are stable with respect to three-dimensional disturbances whose wavevectors are on the hexagonal lattice. This stable region of hexagons is inside the Busse balloon. Therefore, inside the stability boundary of hexagons, either hexagons or rolls are realized depending on the initial conditions. Their result supports the coexistence of stable rolls and stable hexagons (*re-entrant hexagons*), which was experimentally observed for nearly Boussinesq fluids (Assenheimer & Steinberg 1996).

Since the early 1990s, re-entrant hexagons have repeatedly been found as stable Turing patterns in reaction–diffusion systems and as stable convection patterns in Boussinesq and non-Boussinesq fluids. As stable Turing patterns, Ouyang & Swinney (1991*a*,*b*) first observed re-entrant hexagons in a chlorite–iodide–malonic acid (CIMA) reaction. Borckmans *et al*. (1992) and Verdasca *et al*. (1992) obtained the re-entrant hexagons for the two-dimensional Brusselator model, Dufiet & Boissonade (1992*a*) for the Schnackenberg model, Dufiet & Boissonade (1992*b*) for an activator–substrate depletion model, Dewel *et al*. (1996) for a FitzHugh–Nagumo-type model and Mosekilde *et al*. (1998) for the Lengyel–Epstein model related to the CIMA reaction, for example. As stable convection patterns, on the other hand, Dewel *et al*. (1995) and Hilali *et al*. (1995) obtained re-entrant hexagons in the generalized Swift–Hohenberg equation. Recently, Madruga *et al*. (2006) and Madruga & Riecke (2007) carried out fully numerical analyses based on the Galerkin approximation. They gave an explanation for the basic mechanism of experimental results for strongly non-Boussinesq fluids (Roy & Steinberg 2002).

Among the works mentioned above, re-entrant hexagons are predicted by means of the amplitude equations of the cubic order (Verdasca *et al*. 1992) and by means of the quintic equations (Hilali *et al*. 1995). In the latter case, nonlinear saturation never occurs at the cubic order since the primary bifurcation is subcritical. Nonlinear terms of the quintic order are thus needed to realize the saturation. Note that, in both of the cases, the equations involve non-vanishing quadratic nonlinear terms. By contrast, in the present paper, our amplitude equations for *ϵ*=0 lack the quadratic terms. We proceeded to the quintic order since amplitude equations of the cubic order are not generic. Under the symmetry Z_{2} : * ζ*→−

*, Dewel*

**ζ***et al*. (1995) considered a situation where the

*k*=0 mode as well as the six modes on the critical circle are altogether critical. The resultant seven-dimensional amplitude equations of the cubic order, which lack quadratic terms, are shown to possess stable solutions of hexagons for the Swift–Hohenberg equation. In our RB problem with

*k*

_{c}=

*O*(1), the eigenmode with

*k*=0 is strongly damping so that the

*k*=0 mode cannot be regarded as the critical one. Madruga

*et al*. (2006) extended amplitude equations of the cubic order with non-vanishing quadratic terms so as to take account of the spatial modulation.

Returning back to figure 4, we note that above our stability boundary, hexagons are stable with respect to three-dimensional disturbances whose wavevectors are exactly on the critical circle connecting six hexagonal lattice points. At *k*=*k*_{c}, Clever and Busse reported no stable hexagons whereas our result shows stable hexagons for *δ*>0.28. The discrepancy between our result and theirs is remarkable and is considered to be partly due to the truncated nonlinear terms in the amplitude equations higher than the quintic order. It might also be worth noting that the wavevectors of the secondary disturbance in Clever and Busse are on the hexagonal lattice, but are not necessarily on the critical circle. Therefore, the stable region of hexagons they obtained needs to be the same as or narrower than our stable region. A similar situation arises in comparison with Nishida *et al*. (2002). They reported that, for *P*=10, the hexagons are stabilized for *δ*>0.845 at *k*=*k*_{c} under free–free boundary conditions. The lower limit of *δ* is much larger than our value 0.191. To resolve the discrepancy, one possible way is to rederive the amplitude equations by assuming that the centre manifold in (3.7) is spanned by , , , , , and *δ*_{Kuo}, and supplying the equations for the centre modes with the additional equation . Rederiving the equations is straightforward, but is so complicated that it will form future work. As is seen in figure 4, we find that the Rayleigh number of our stability boundary decreases with the decrease of *k*. The boundary ends up on the neutral curve at approximately *k*≃2.3, where the weakly nonlinear theory should be the most reliable. We have to note, however, that the 1 : resonance takes place at the cubic order for *k*=2.341 and the 1 : 2 resonance takes place at the quintic order for *k*=2.165. The interaction points are marked in figure 4 by open and filled circles.7 Before emphasizing the validity of our quintic-order amplitude equations near *k*≃2.3, therefore, we have to know whether or not these resonant interactions alter the stability of hexagons. The effect of the resonance on the stability of hexagons will thus be the subject of future work.

## Acknowledgments

The authors thank Edgar Knobloch for pointing out the reference of Clever & Busse. They also thank an anonymous referee for constructive comments and pointing out the reference of Kuo.

## Footnotes

↵The exception was due to Busse (1962); in his extensive work on the role of the temperature dependence of physical properties, he examined the effect of quartic-order terms on the stability of rolls and hexagons under asymmetric boundary conditions.

↵We call a solution axial if its isotropy subgroup is axial with one-dimensional fixed-point subspace. See Golubitsky & Stewart (2002) for example.

↵Since our problem is parametrized by

*δ*, the centre manifold should also depend on*δ*, i.e. . The equation needs to be added to the equations for the centre modes (3.6). An increase in the dimension of the centre manifold makes the derivation much more complicated. We thus avoid such a situation by simply retaining*σζ*_{1}term in (3.8) without adding . We evaluate the*δ*dependence of*σ*from the linear stability characteristics and set*σ*≃*σ*_{δ}*δ*for small*δ*.↵Changing the number of eigenfunctions, we checked the convergence of the expansion in eigenfunctions. It is concluded that, with 20 eigenfunctions, all the coefficients converged within five decimal digits (see also Fujimura 1997).

↵Busse (1962) showed that, for

*P*≫1, the quartic-order terms exert a stabilizing effect on hexagons for*δ*>2.4 whereas they destabilize rolls for*δ*>6.2. Qualitatively consistent results are obtained based on our quartic truncation, although the parameter range is outside the range of validity of the amplitude equations, i.e. for*P*=1000, hexagons are stabilized for*δ*>1.87 whereas rolls are destabilized for*δ*>5.95.↵The weakly nonlinear analysis is by means of the amplitude expansion scheme instead of the parameter expansion used by Kuo.

↵On the hexagonal lattice, the 1 : resonance at the quadratic order was investigated by Daumont

*et al*. (1997) for the Swift–Hohenberg equation, whereas the 1 : 2 resonance at the quadratic order was examined by Fujimura (2008) for a two-layered RB problem.- Received November 23, 2007.
- Accepted May 6, 2008.

- © 2008 The Royal Society