Whirling hexagons and defect chaos in hexagonal non-Boussinesq convection

We study hexagon patterns in non-Boussinesq convection of a thin rotating layer of water. For realistic parameters and boundary conditions we identify various linear instabilities of the pattern. We focus on the dynamics arising from an oscillatory side-band instability that leads to a spatially disordered chaotic state characterized by oscillating (whirling) hexagons. Using triangulation we obtain the distribution functions for the number of pentagonal and heptagonal convection cells. In contrast to the results found for defect chaos in the complex Ginzburg–Landau equation, in inclined-layer convection, and in spiral-defect chaos, the distribution functions can show deviations from a squared Poisson distribution that suggest non-trivial correlations between the defects.

135.3 models, we investigate here the dynamics of hexagonal patterns in the presence of rotation on the basis of the standard hydrodynamical description of non-Boussinesq convection [4], [21]- [23]. Thus, we use the Navier-Stokes equations coupled to the heat diffusion equation (called NONBE henceforth, see below) corresponding to convection in a layer of water near its density maximum. We perform Galerkin stability calculations for the hexagonal pattern and direct numerical simulations of the NONBE and find spatio-temporal chaos that involves an oscillatory instability of the hexagons as well as defects in the hexagonal background pattern. In contrast to the defect chaos of the oscillating hexagons studied in the weakly nonlinear regime [5], here the oscillatory instability does not correspond to spatially homogeneous oscillations but involves a spatial modulation of the oscillation amplitude, and it leads to the nucleation of defects in the hexagonal background pattern.

Basic equations
A horizontal fluid layer of thickness d that is heated from below and cooled from above and that is rotated about a vertical axis with angular velocity ω k is described by the Navier-Stokes equation for the momentum, (1) the continuity equation, and the heat equation, Here u = (u 1 , u 2 , u 3 ) is the fluid velocity, T the temperature, ρ the density, p the pressure, g the acceleration of gravity, ν the viscosity, λ the heat conductivity, c p the specific heat of the fluid, and k is the unit vector in the z-direction. The Kronecker delta is given by δ i j , and i jk is the unit antisymmetric tensor of rank 3. We assume the Einstein summation convention. We have omitted viscous heating, volume viscosity and the centrifugal force. In our analysis we take realistic rigid boundary conditions at the top and bottom plates and keep the top and bottom temperatures fixed, Here T 0 denotes the mean temperature and T > 0 is the temperature difference across the layer. In the lateral directions we choose periodic boundary conditions or introduce a suitable ramp in the imposed temperature gradient (see below).

135.4
In this paper we focus on weakly non-Boussinesq convection. We therefore keep the temperature dependence of the various fluid properties to leading order and expand them about the mean temperature T 0 in the form [4] ρ(T ) where ρ 0 , ν 0 , λ 0 , and c p0 denote the values of the respective quantities at the mean temperature T 0 . The dots denote higher-order terms to be neglected in the following. To facilitate the analysis of our results in line with the typical experimental procedure we keep T 0 constant when changing the main control parameter T . Then the quantities γ i / T, i = 0, 2, 3, 4, which describe, respectively, the slope of the density, viscosity, heat conductivity, and heat capacity at T 0 , are fixed as well. This implies that the γ i , which give the change of the respective fluid property across the fluid layer, are linear in the main control parameter T and therefore need to be adjusted with increasing T (see below). The usual heat expansion coefficient at T = T 0 is given by α 0 = γ 0 / T , but going beyond the Boussinesq approximation the curvature of ρ(T ) at T 0 , which is proportional to γ 0 γ 1 / T 2 , also comes into play. Note that for constant fluid parameters the volume viscosity term and the centrifugal force can be absorbed into the pressure term. This is not true when the properties depend on the temperature. While the z-dependence of the volume viscosity generates only a higher-order correction, neglecting the centrifugal force is only possible as long as the rotation rates are not too large for a given system size L, i.e. as long as γ 0 Lω 2 /g is small.
Since the fluid velocities are small compared to the sound velocity we make the anelastic approximation [24] (see also [22]) and neglect the time derivative in the continuity equation (2). This simplifies the computation considerably since it reduces the number of evolution equations. Furthermore, v i becomes a solenoidal field, which can be represented in the standard poloidaltoroidal decomposition by two velocity potentials [22]. This automatically enforces the mass conservation.
The conduction solution ( v = 0) of (3) with (5)-(7), (9) and (10) is given by We rewrite the temperature T in terms of the deviation from the conductive profile (11), neglecting its contribution of O(γ 2 3 ), We then obtain as the final dimensionless equations The dimensionless boundary conditions are The nondimensionalized fluid parameters (7)-(10) now read We consider the non-Boussinesq effects to be weak and keep in all material properties only the leading-order temperature dependence beyond the Boussinesq approximation. Therefore the γ 1 -term appears explicitly in (13), while in all other terms it would constitute only a quadratic correction just like the terms omitted in (7)- (10). Correspondingly, we expand the denominators in (13), (15) that contain material properties to leading order in the γ i . In analogy to [4] (equation (6.7)), we further omit non-Boussinesq terms that contain cubic nonlinearities in the amplitudes v i or , as they arise from the expansion of the advection terms v j ∂ j (v i /ρ) and (v j /ρ)∂ j , once the temperature dependence of the density is taken into account. Since we will be considering Rayleigh numbers up to twice the critical value, which implies enhanced non-Boussinesq effects, these approximations may lead to quantitative differences compared to the fully non-Boussinesq system, even though the temperature dependence of the material properties themselves may quite well be described by a linear (or quadratic in the case of the density) approximation.

Computational methods
We solve the NONBE (13)- (15) with top and bottom boundary conditions given by (16) and the material parameters given by (17)- (20) numerically using a number of approaches. As usual, the stability properties of the patterns are determined by Galerkin methods. We assume an infinite extension of the layer in the lateral directions, which is captured by a Fourier expansion on a hexagonal lattice. The Fourier wavevectors q are constructed as linear combinations of the hexagonal basis vectors b 1 = q(1, 0) and b 2 = q(1/2, √ 3/2) as q = mb 1 +nb 2 with the integers m and n in the range |mb 1 +nb 2 | ≤ n q q. The largest wavenumber is then n q q and the number of Fourier modes retained is given by 1+6 n q j=2 j . Typically we used n q = 3. The top and bottom boundary conditions are satisfied by using appropriate combinations of trigonometric and Chandrasekhar functions in z [4,21]. In most of the computations we used n z = 6 modes for each field in equation (13)- (15). The linear analysis yields the critical Rayleigh number R c (γ i ) as well as the critical wavenumber q c (γ i ). In the large-Pr-number regime, to be considered in this paper, the bifurcation to convection is stationary.
To investigate the nonlinear hexagon solutions, we start with the standard weakly nonlinear analysis to determine the coefficients of the three coupled amplitude equations for the modes making up the hexagonal pattern. This allows us to get a first insight into the Küppers-Lortz instability of the rolls, and for the weakly non-Boussinesq case the transition from hexagons to rolls as well as the transition to oscillating hexagons. To obtain the fully nonlinear solutions requires the solution to a set of nonlinear algebraic equations for the expansion coefficients with respect to the Galerkin modes. This is achieved with a Newton solver for which the weakly nonlinear solutions serve as convenient starting solutions. The fully nonlinear solutions are then tested for amplitude and modulational instabilities following the standard methods. Alternatively, we solve the NONBE by direct numerical simulation to describe in general the temporal evolution of the pattern. While this allows for a precise check of the Galerkin results, its main purpose is to capture the complex spatio-temporal dynamics resulting from the instabilities of the hexagonal pattern, which are not accessible in the Galerkin approach. The simulation code involves slight extensions of a previous, well proven version [23,25,26] to include the effect of the non-Boussinesq corrections.

Results
Instead of extensive parameter studies, we concentrate in this paper on a specific, interesting scenario that should be experimentally realizable. We focus our investigation on a fluid with a fairly large Prandtl number (water). For smaller Prandtl numbers our previous results using a Swift-Hohenberg equation and also some initial test runs using the Navier-Stokes equations indicate that rolls [4], which may be nucleated in defects or at the wall, take over the hexagon pattern even for relatively small values of the Rayleigh number. Furthermore, we have confined ourselves in this work to one representative rotation rate ( = 65, see below).
To obtain strong non-Boussinesq effects the layer would have to be taken quite thin. However, then the dimensionless rotation rates obtained for a given physical rotation rate are relatively low. Therefore, we decided to consider water near its density maximum as the convecting fluid. Strong non-Boussinesq effects at larger layer thicknesses are also obtained with Stability limits of stationary hexagon pattern for the parameters given in table 1 (n q = 3, n z = 6, except for dark green diamond, which is for n q = 5, n z = 8): long-wave steady (black circle), long-wave oscillatory (red squares), oscillatory short-wave instability (green diamonds), steady short-wave instability (blue triangles).  glycerin. However, since the rotation rate is made dimensionless using the viscous diffusion time d 2 /ν, the extremely large viscosity ν of glycerin yields only small dimensionless rotation rates. The second row in table 1 gives the non-Boussinesq parameters used in all computations in this paper. For = 65 they correspond to a critical temperature difference of 10 • C across the layer at a mean temperature of 12 • C. Note that for a given thickness of the fluid layer the critical Rayleigh number increases significantly with the rotation rate. Therefore, the coefficients γ (c) i , which give the non-Boussinesq effects at onset, also increase strongly with the rotation rate. Consequently, while in the absence of rotation the system is essentially Boussinesq for the chosen layer thickness of d = 0.46 cm (first line in table 1), it is strongly non-Boussinesq for = 65. Table 2 gives the critical Rayleigh number and the critical wavenumber for the investigated layer of water with thickness d = 0.46 cm with and without rotation.
We first determine numerically the linear stability of the stationary hexagon patterns as a function of the wavenumber and the control parameter = (R − R c (γ i ))/R c (γ i ) for = 65 and the corresponding non-Boussinesq coefficients γ (c) i as given in table 1. As indicated before, the γ i are linear in the temperature difference T and are therefore given by γ i = γ (c) i (1 + ). All results were obtained with n q = 3 and n z = 6. In the parameter regime investigated these Table 2. Dependence of the critical Rayleigh number and of the critical wavenumber on the rotation rate for a layer thickness of d = 0.46 cm. The critical Rayleigh numbers in the Boussinesq approximation (γ i = 0) and for γ i = γ i ( T = T c ) are denoted by R c0 and R cγ , respectively. Similarly for the critical wavenumber q c0,γ . The physical rotation rate for this cell thickness is given by ω and the vertical thermal diffusion time by τ th . cut-off parameters are sufficient. As shown in figure 1 the hexagon pattern can undergo a variety of linear instabilities in the nonlinear regime. Close to threshold the dominant instability is a long-wave instability, i.e. its growth rate vanishes for vanishing modulation wavenumber β (Floquet parameter) but grows quadratically with β for small β. Very close to onset it can be steady (black circles). For slightly larger , however, it becomes oscillatory (red squares). As the control parameter is increased a steady instability at finite Floquet parameter limits the stability region of the hexagons for large wavenumbers (blue triangles). Increasing the control parameter, on the low-wavenumber side the stability region becomes limited by an oscillatory instability that sets in with a finite Floquet parameter (green diamonds). For non-Boussinesq coefficients that differ slightly from those of table 1 (by O(20%)) an oscillatory instability with vanishing Floquet parameter determines the stability limit for low q and ≈ 1. For the parameters in table 1 it is, however, pre-empted by the Hopf bifurcation at finite modulation wavenumber. The spatially homogeneous and the short-wave oscillatory instabilities are similar to the oscillatory instability of hexagons that is found within the framework of three coupled Ginzburg-Landau equations [5,6,10,11]. There it replaces the steady instability of hexagons that in the absence of rotation leads to an unstable mixed-mode solution. With rotation it turns into a supercritical Hopf bifurcation at which the oscillating-hexagon solution branches off the stationary hexagons.
Numerical simulations support the expectation that the modes corresponding to the Hopf bifurcation are similar to the oscillating hexagons. In a small computational domain of length L = 2 × 2π/q, which contains only four convection cells, the oscillatory side-band instability is suppressed and the homogeneous Hopf bifurcation arises when the parameters are chosen sufficiently far beyond the stability limit corresponding to the short-wave oscillatory instability (green diamonds in figure 1). This instability leads to an elliptic deformation of the convection cells which oscillates in time, giving the cells an appearance as if they were rotating or in fact whirling. In sufficiently large domains the side-band instability leads to a spatial modulation of these oscillations as shown in figure 2(a) (L = 8 × 2π/q with 64 × 64 Fourier modes).
As the instability develops and the oscillation amplitude grows the regular spatial structure is destroyed and defects arise in the pattern as pictured in figure 2(c). This temporal evolution is illustrated in movie2. While the apparently chaotic dynamics appear to be reaching a statistically steady state, they do not, however, persist, and around t = 130 a relative rapid transition to an ordered, stationary hexagon pattern occurs. It has a slightly larger wavenumber q, which is in the stable regime (cf figure 1). During this chaotic transient the overall disordered pattern  intermittently exhibits relatively large domains of ordered hexagonal patterns. The duration of the transient increases with as illustrated in figure 3, presumably reflecting the narrowing of the band of stable wavenumbers with increasing . To monitor the activity of the pattern we consider as a representative measure the temperature field in the mid-plane of the layer and plot the 'whirling' activity W given by the L 1 -norm of the time derivative of its Fourier modes 4 normalized by the L 1 -norm of the Fourier modes themselves. As is increased the fluctuations in this quantity decrease, which is expected to reduce the probability of excursions to one of the linearly stable stationary states, and eventually for = 1.1 our simulations did not end up in a stationary pattern for times as large as t = 400.
The duration of the transients is presumably related to the size of the system. We therefore also performed simulations in larger systems of size 16 × 2π/q (128 × 128 Fourier modes). For this system size the dynamic state persists for quite a long time (t max ≈ 300τ th ≈ 5 h) even 135.10 for control parameters as low as = 0.87 (compare with the transients for the smaller system shown in figure 3). This increased duration of the transient is due to the longer time it takes for the defects in the pattern to annihilate. For = 1.0, a snapshot of the pattern and of the corresponding temporal derivative is shown in figures 4(a) and (b). The latter makes the spatially intermittent behaviour of the dynamics more apparent. The dynamics are illustrated in the two corresponding movies (movie4A, movie4B).
For the characterization of weakly disordered hexagon patterns it is often advantageous to locate defects in the pattern by demodulating the pattern using the (three) main wavevectors comprising the pattern (e.g. [13,8]). This allows the identification of the penta-hepta defects that are characteristic of weakly disordered hexagon patterns and which consist of correlated convection cells with five and seven neighbours, respectively. Obviously, this approach requires that the Fourier spectrum of the pattern exhibits well defined peaks at the three hexagon modes.
In the strongly disordered patterns obtained in our simulations the spectrum does not exhibit clear peaks. Demodulation is therefore not possible and other techniques for the analysis of the patterns have to be brought to bear. We use a triangulation of the pattern to locate the defects. To that end we map the convection pattern to a point lattice by identifying the local minima in the temperature with nodes in a lattice. To avoid spurious minima we require that the minima fall below a certain threshold. The coordination number of each node is then obtained using a Delaunay triangulation, which is the unique triangulation obtained by the requirement that the circle through the corners of each triangle does not contain any other node of the pattern. We use the program triangle [27] to perform the triangulation. three modes that make up the hexagonal pattern (see e.g. figure 25 in [23]). Isolated pentagons or heptagons, however, do not correspond to single dislocations.
In various previous analyses of defect-dominated spatio-temporal chaos it was found that the distribution function of the number of defects in the patterns gives some insight into the dynamics of the system. More specifically, in the defect chaos obtained in the complex Ginzburg-Landau equation [14], in electroconvection of nematic liquid crystals [17], in undulation chaos in inclined-layer convection [18,19] as well as in spiral-defect chaos in convection at low Prandtl numbers [20] the statistics were found to be close to a Poisson-type distribution. Such a distribution arises if defects of opposite topological charge are created randomly in pairs with a fixed probability and are annihilated with a rate that is proportional to the square of the density of the defects. The Poisson statistics therefore suggest that the defects move like uncorrelated random walkers. In contrast, for the penta-hepta defect chaos that was found in hexagon patterns in the presence of rotation within the framework of a Swift-Hohenberg-type equation the distribution function was found to be considerably broader [8]. The origin of this broadening was identified to be the induced nucleation of dislocations. In this process an already existing penta-hepta defect induces the nearby nucleation of a dislocation pair, which then combines with the penta-hepta defect to form two penta-hepta defects. As a result, the creation rate for dislocations and the ensuing penta-hepta defects depends on the defect density [8,28].
Using the triangulation procedure, we have measured the distribution functions for the number of pentagons and heptagons. For the state corresponding to figure 4 they are presented in figures 5(a) and (b), where the blue circles and red triangles denote the distribution functions for pentagons and heptagons, respectively. As mentioned previously, in strongly disordered patterns pentagons and hexagons are not necessarily bound in penta-hepta defects. Correspondingly, the number of pentagons does not always agree with the number of heptagons. The distribution function for their difference is given by the solid line in figure 5(b) along with the distribution function for rectangles (blue circles) and octagons (red triangles).
As a first step it appears reasonable to assume that the dynamics of the pentagons and heptagons are captured by a kinetic model that is analogous to the description of dislocations [14]. For periodic boundary conditions [18] this leads then to the standard Poisson distribution for each species, where α ≡ N 2 and I (2 √ α) is the modified Bessel function. In some simulations we observe that a given whirling cell creates a new convection cell nearby, which implies the creation of defects. The newly created cell then often merges with the whirling cell. Conversely, defects in the pattern induce whirling activity. This mutual nonlinear enhancement of defects and whirling activity suggests that the existence of defects may lead to the creation of further defects, somewhat similar to the induced nucleation in penta-hepta defect chaos [8]. As shown by the dashed curves in figure 5(a), for the state of figure 4, the Poisson distribution (21) does not fit the pentagon and the heptagon distributions very well. The distributions therefore suggest that the creation rates for the pentagons and heptagons are not independent of their respective densities. In some preliminary parameter studies involving non-Boussinesq coefficients that may not correspond to a realistic convection cell with water as fluid, we have found distributions that are fitted quite well by the Poisson distribution and others that are even broader than the distributions obtained in figure 5. This suggests that a mechanism akin to induced nucleation may be operative in this system, with its significance depending on the specific parameter values [29]. A study of the dependence of the distribution on the physical parameters of the fluid is beyond the scope of the present paper. Figure 4(b) shows that in the disordered state not all cells are active and at times sizable patches of relatively ordered quiescent domains arise. To characterize this intermittent behaviour we consider the distribution function for the time derivative of the temperature as obtained from all snapshots of the type shown in figure 4(b). The logarithmic scale of figure 6 reveals an exponential decay of the distribution function for = 1.0. This exponential dependence is consistent with the visually apparent spatio-temporal intermittency of the activity of the pattern. The pronounced peak exhibited for = 0.87 represents the large spatial domains in which the hexagons are ordered and quiescent.
In all simulations so far we have employed periodic boundary conditions. In experiments boundary effects may strongly perturb the patterns. We therefore have also performed some simulations in which we mimic a circular container by applying a strong radial subcritical ramp in the Rayleigh number that suppresses any convection outside a certain radius. In previous simulations [23] this was seen to provide a reasonable idea of the impact of boundaries. Experimentally, it has been found that boundaries tend to induce defects and may also promote the formation of rolls, which could invade the system and replace the hexagonal pattern [30]. For the parameter regimes investigated here we find that no patches of rolls invading from the side replace the hexagons. Instead the precession of the cells near the boundary and the associated formation of defects triggers persistent dynamics in the whole convection cell even for parameter values for which in the absence of the boundaries the pattern eventually becomes ordered and stationary. Figure 7(a) shows a snapshot for = 0.87 with the associated temporal evolution displayed in movie7A. A snapshot for a larger system size and = 1.0 is shown in figure 7(b) (movie7B). We therefore expect that the defect-dominated spatio-temporal chaos will also be accessible in experiments.
We have not investigated in detail the nonlinear evolution ensuing from the steady instability at finite wavenumber that limits the wavenumber band on the high-q side. In preliminary simulations the instability led to the formation of defects in the pattern, which eventually, however, disappeared, yielding a stationary, regular hexagonal pattern.

Conclusion
In this paper we have considered hexagonal non-Boussinesq convection focusing on rotating systems, in which the chiral symmetry is broken. The main motivation for this work stems from the interesting chaotic states that had been identified in previous investigations of the stability of hexagonal patterns and of the dynamics resulting from their instabilities. The chaotic states were associated with the secondary instability to oscillating hexagons [5,6] and with the nucleation of dislocations [8,13,28], respectively. These results were obtained within order-parameter and Ginzburg-Landau models. To make closer contact with experimental systems we have investigated in this paper the effect of rotation on realistic hexagonal convection in a thin layer of water in a temperature regime near the anomalous density maximum.
Using a numerical stability analysis and direct numerical simulations of the Navier-Stokes equations we have identified an oscillatory side-band instability, which leads to spatially 135.14 modulated whirling hexagons that in turn become disordered in space and exhibit irregular persistent dynamics, which coexist with linearly stable regular hexagon patterns.
Coexistence of spatio-temporal chaos with ordered patterns has been observed previously in roll convection at low Prandtl numbers [1,31]. There the spiral defect chaos is sustained by the mean flow, which compresses the rolls so as to push them locally across their stability limit with respect to the skew-varicose instability, which then leads to the formation of defects [32]. This mechanism for sustaining dynamics is reminiscent of that observed earlier in a much smaller system [33,34].
The method employed in [32] to analyse the disordered pattern, which extracts the local wavenumber and orientation of the rolls, is not applicable to the strongly disordered hexagonal patterns found in our simulations. We therefore use Delaunay triangulation to measure the distribution functions for the number of pentagons and heptagons. Both show deviations from a squared Poisson distribution, which suggests that their dynamics (including their creation and annihilation) exhibit correlations that are absent in the dynamics of dislocations in the complex Ginzburg-Landau equation [14] and in inclined-layer convection [18,19] as well as in the dynamics of spiral defects in spiral-defect chaos [20]. It is, however, not quite clear whether a kinetic description of the pentagons and heptagons similar to that of dislocations or of pentahepta defects [8,14,18] is in fact appropriate. To characterize the intermittent appearance of ordered domains in the chaotic states we have determined the distribution function for the temporal derivative of the temperature field. It exhibits an exponential decay reflecting the spatio-temporal intermittency. For smaller Rayleigh number a pronounced peak appears, which reflects the large ordered, quiescent domains.
In this paper we have focused on one set of fluid parameters (cf table 1). For these values the oscillatory instability with finite modulation wavenumber pre-empts a spatially homogeneous oscillatory instability. We found, however, that the situation can be reversed if the non-Boussinesq parameters are changed only slightly. In these cases, which may not correspond to a realistic convection cell with water as fluid, we have found defect-chaotic states with distribution functions that are significantly wider than the squared Poisson distribution suggesting again a nonlinear interaction between the whirling of the hexagons and the formation of defects, as well as distribution functions that are well described by the squared Poisson distribution. In addition, we have identified a state in which the whirling is spatially and temporally intermittent in the form of bursts while the underlying hexagonal planform has essentially no defects [29].