## Abstract

Motivated by the need to understand better the dynamics of non-Fourier fluids, we examine the linear and weakly nonlinear stabilities of a horizontal layer of fluid obeying the Maxwell–Cattaneo relationship of heat flux and temperature using three different forms of the time derivative of the heat flux. Linear stability mode regime diagrams in the parameter plane have been established and used to summarize the linear instabilities. The energy balance of the system is used to identify the mechanism by which the Maxwell–Cattaneo effect (i) introduces overstability, (ii) leads to preferred stationary modes with the critical Rayleigh and wavelengths either both increasing or both decreasing, (iii) gives rise to instabilities in a layer heated from above, and (iv) enhances heat transfer. A formal weakly nonlinear analysis leads to evolution equations for the amplitudes of linear instability modes. It is shown that the amplitude of the stationary mode obeys an equation of the Landau–Stuart type. The two equally excitable overstable modes obey two equations of the same type coupled by an interaction term. The evolution of the different amplitudes leads to supercritical stability, supercritical instability or subcritical instability depending on the model and parameter values. The results are presented in regime diagrams.

## 1. Introduction

The relationship between the temperature of a fluid and the heat flux it produces has been of great importance for centuries because it is a requirement for understanding real-life problems that are influenced by temperature. The empirical relationship between temperature and the heat flux it produces, proposed by Fourier in 1822, leads to a parabolic equation for the evolution of temperature. When an initial value problem is considered for this equation, it produces a solution that immediately extends the temperature to very large distances from the source due to the infinite speed with which the temperature spreads. This shortcoming of the Fourier law motivated studies in search of alternative relationships, which are more acceptable physically. Notable among those early studies is that by Maxwell [1] on kinetic theory of gases. It is now recognized that Fourier law is adequate for some fluid states. However, a vast number of other fluid states, particularly those with small dimensions like nano-fluids and hence associated with relaxation times comparable with their thermal diffusion times, require different relationships. A number of such relationships that replace the parabolic heat equation by a hyperbolic one have been proposed. The hyperbolic equation takes the form of a wave equation and hence provides solutions with finite phase speeds, and waves produced by such an equation are referred to as second sound (or heat) waves [2]. In 1948, Cattaneo [3] introduced a time derivative to Fourier law in solids and in 1950 Oldroyd [4] made an extensive formulation of the constitutive equations of a continuum including a time derivative. In 1969, Fox [5] proposed a time derivative in the form of an objective derivative for solids. In 1972, Carrasi & Morro [6] introduced a time derivative to a Newtonian viscous fluid. The Fourier relation
**q*** is the heat flux, *T** the temperature and *K* the thermal conductivity, was eventually replaced by the relation
*δ*/*δt** here represents the time derivative following the motion so that
*t** is the time and **u*** the velocity (see also Morro [7] and references therein). We shall refer to all fluids for which the heat flux equation (1.2) is relevant, as Maxwell–Cattaneo fluids. These types of fluids form a class of fluids which have been of topical theoretical interest because they are relevant to a wide range of applications, ranging from planetary science (e.g. [8–10]) to liquid helium [11,12]. Four decades ago, it was assumed that the relaxation time is so small that the time derivative in the equation (1.2) is practically insignificant [13], and it can even be of order picoseconds for some fluids [14–16]. For such fluids the Fourier law applies. According to the studies [17,18], there are materials for which the relaxation time is not small, such as sand (21 s) and biological tissue (1–100 s) and these have been the subject of a number of studies (e.g. [19–22]). Recently, Hayat *et al.* [23] studied the impact of Cattaneo–Christov heat flux in a Jeffrey fluid flow with homogeneous–heterogeneous reactions. Neuhauser [24], in her investigation of thermal relaxation in superfluid Helium-3, found that *τ*_{r} has values in the range 30–400 s. Mohammadein [25] in his study of the thermal relaxation time between two-phase bubbly flow found that the relaxation time varied between 10^{−3} and 3 s. Vidal [26] examined the propagation of sound waves near the normal-superfluid transition in 4He liquid and made a comparison between experiments and classical dynamical theories to find that the relaxation time has a range of values in the interval (2×10^{−5}/*π*,2/(10*π*)) s. In addition, nanofluids in small-size industrial products are associated with relaxation times that are comparable to their diffusion time scales and we will see below that the importance of the role of the relaxation time is measured by its ratio with the diffusion time of the system (see [9] for references to a variety of other applications). In general, the relaxation time effect is significant for any material in a situation in which the relaxation time is comparable with the thermal diffusion time. In particular, micro- and nano-fluids involve small dimensions which lead to small diffusion times and hence may be prone to the Maxwell–Cattaneo effect.

A number of suggestions have been made for replacing the time derivative of the Maxwell–Cattaneo equation (1.2) by an objective derivative. Objective derivatives have generally been more popular in elasticity theory and can take a number of forms. There are three well-known objective derivatives due to Jaumann [27], Truesdell & Noll [28] and Green & Naghdi [29]. The Jaumann form of objective derivative is considered more suitable for fluids because it takes into account the local spin (or rotation). However, the introduction of an objective derivative to fluid motions was introduced by Morro, [30] in his study of acceleration waves in thermo-viscous fluids and Straughan & Franchi [31] were the first to use an objective derivative in stability theory. The derivative they used has the form
** ω*** is the vorticity (representing the spin), to study the linear stability of the Benard layer. They found that the stability depends on two material parameters: the Prandtl number

*ν*is the kinematic viscosity and

*κ*the thermal diffusivity, and

*d*is the thickness of the layer and

*τ*

_{κ}(=

*d*

^{2}/

*κ*) is the thermal diffusion time. The parameter

*C*, corresponding to the ratio of half the relaxation time to the thermal diffusion time, is introduced by the Maxwell–Cattaneo (hereinafter referred to as M–C) effect. It was shown in [31] that such a derivative allows instabilities of a layer heated from above in the form of stationary waves and, in addition, it introduces time-dependent instabilities to the layer heated from below provided

*C*is large enough.

The parameter *C* plays the pivotal role in all the instabilities occurring in fluids obeying equation (1.2) whatever the form of the derivative. In the analysis below (see equation (3.11) below), we find it convenient to use the scaled parameter, *C*_{1}, which magnifies *C* by a factor of *π*^{2} (i.e. *C*_{1} is about 10 *C*).

Lebon & Cloot [13] used the derivative in the form
*C*.

Another approach was made by Christov [32] who used a general coordinate transformation method to express the time derivative (*δ*/*δt**) in a form independent of the frame of reference

The stability of the Benard layer with the heat flux obeying the Christov form of the objective derivative has received considerable attention, Straughan [10,34,35], Tibullo & Zampoli [36], Stranges *et al.* [37], Haddad [38], Bissell [39], Stranges *et al.* [40], Bissell [41]. The paper [40] makes a detailed analysis of the relationship between the M–C number, *C*, and the Knudsen number, *K*_{n} (=the ratio of the mean free path in a rarefied gas to the characteristic length scale of the system, here the thickness of the Benard layer) and obtained the relation

The motivation for this study is multi-fold. First, we discuss further the linear instabilities identified by the relations (1.4) (hereinafter referred to as the SF model) and the relation (1.7) (hereinafter referred to as the LC model) as well as that by the equation (1.8) (referred to as the Christov model). We point out some new results and establish mode regime diagrams, displaying the regions of existence of stationary and overstable modes as well as the regions where they are linearly unstable. We then use these results to address the questions: What is the role played by the M–C effect in enhancing a particular instability? What makes the different objective derivatives give rise to different types of instability? When do the instabilities of a layer heated from above exist and what mechanism drives them? It turns out that the use of the energy balance equation shows that the M–C effect is equivalent to defining an *effective thermal diffusivity*, *κ*_{e}, which depends on the M–C number, the Prandtl number, the wavenumber, the frequency and the spin. It is the versatility of *κ*_{e} that leads to the appearance of these effects.

Secondly, we adopt a formal weakly nonlinear stability analysis to investigate the nonlinear development of the linear instabilities with a view to identifying which ones are stabilized and which ones suffer subcritical instabilities, if any. Now the preferred overstable mode has a critical frequency, *ω*_{c}, that occurs in the form *ω*^{2}_{c} so that there are two equally excitable waves travelling in opposite directions. Sometimes it is argued that on a long time scale, the two waves will be far apart and their interaction can be neglected. Here, we take the interaction effect into account and find that it can play an important part in the evolution of the wave amplitudes. In the case of stationary modes, only one wave is present. In all cases, the amplitudes evolve on a long length scale of the same order of magnitude as that of the slow time scale and as a consequence first order (in both slow time and long length scales) evolution equations of the Landau–Stuart type are obtained. In the case of the overstable mode, a coupled pair of equations is obtained. A detailed analysis of the possible instabilities in the whole parameter space reveals that the nonlinear instabilities that can occur are very much dependent on the model considered and on the relative values of the material parameters.

In §2, we formulate the problem and in §3 we discuss the linear stability, investigate the energy source of the new instabilities introduced by the M–C effect and discuss the influence of the M–C effect on heat transfer. In §4, we examine the nonlinear stability of the three models and in §5 we make some concluding remarks.

## 2. Formulation of the problem

Consider a horizontal layer of viscous non-Fourier fluid of thickness *d*, contained between two horizontal planes. The lower and upper planes are maintained at the uniform temperatures, *T*_{B},*T*_{T}, respectively. We consider a Cartesian system of coordinates *O*(*x*,*y*,*z*) with the origin, *O*, half-way between the two planes, the *z*-axis pointing vertically upwards and *Ox*,*Oy* axes are any two mutually perpendicular horizontal directions.

The equations governing the motion of the fluid are the momentum equation, the continuity equation, the conservation of heat flux equation, the heat flux relation and the equation of state
*ν*, and thermal diffusivity, *κ*, have been introduced in §1 above and *ρ*_{0},*T*_{0} are reference values of density and temperature. In writing equation (2.1), we have assumed the Boussinesq approximation in which case the variations in density are neglected except when they occur in the gravity term. The term *δ***q***/*δt** depends on the model taken. Here we will discuss the three models (1.4), (1.7) and (1.8).

We assume a basic state
*β*, is positive for a layer heated from below and negative for a layer heated from above. We assume that the basic state is disturbed slightly so that the variables take the form
*ε* is a small dimensionless amplitude factor. We define characteristic units of distance, time, velocity, temperature, heat flux and pressure by *d*,*d*^{2}/*κ*,*κ*/*d*,|*β*|*d*,|*β*|*κ* and *ρ*_{0}*νκ*/*d*^{2}, respectively. This choice of characteristic quantities reduces to that used for the classical Benard layer and is adopted here in order to facilitate comparison. The dimensionless equations of the perturbations take the form
*α* is defined by
*Ra* is defined by
*y*-component of equation (2.9).

## 3. Linear stability

Here, we establish mode regime diagrams for the existence and preference of stationary and overstable modes for the three models [13,31,34], summarize their main properties and investigate the mechanisms for driving them.

The linear equations are obtained from (2.9)–(2.12) by setting *ε*=0. We assume that the perturbation variables take the form
*x*-axis so that the *y* dependence can be suppressed. If the *y* dependence is required, it can be obtained by modifying the final result without much effort. Here *k* is the wavenumber and *γ* is the growth rate. If *Re*(*γ*)>0, the system is unstable and if *Re*(*γ*)<0, it is stable. If *Re*(*γ*)=0, the system is neutrally stable. This is the state of interest in linear stability theory. If the energy of the least stable neutral state is increased, instability will ensue. We therefore assume that *γ*=i*ω* and our objective is to find the neutral state with the minimum energy (as represented by *Ra*), which is the least stable, and hence the most likely to be unstable once the energy is increased beyond its minimum value. This state is referred to as the preferred (or critical) state, and the associated Rayleigh number, wavenumber and frequency are referred to as the critical or preferred mode.

Substitution of the expressions (3.1) for all the variables in the equations (2.9)–(2.12) with *ε*=0, and after some manipulations, leads to

The application of the boundary conditions to (3.2) leads to the two uncoupled categories of solutions
*z*. Consideration of the dependence of the Rayleigh number, *Ra*, on *n* shows that the fundamental even mode (*n*=1) provides the smallest *Ra*. We therefore take the solution
*A* is an arbitrary constant for the linear theory, but will be determined by the nonlinear theory of §4 below. The expressions for the other variables are

### (a) The linear stability of a layer heated from below

Here, we examine briefly the linear stability properties of the three models, construct the mode regime diagrams for their instabilities and discuss the differences between their stability characteristics. We will then use the energy balance equation for each model to identify the mechanisms driving these instabilities.

#### (i) The SF model

The SF model for a layer heated from below allows both stationary and overstable modes [31]. The stationary mode, obtained from (3.13) by setting *ω*=0, *π*, we have defined the scaled parameters
*m* refers to the minimum value. *The notation* (3.16) *will be used throughout the paper*. We note that the mode (3.15) exists only if *C*_{1}<1 because as *C*_{1} increases to 1,

The overstable mode has the minimum
*C*_{1} exceeds a certain bounding value, *C*_{b}, given by
*C*_{1}<*C*_{b} because the overstable mode does not exist in this range of *C*_{1}, and is given by (3.17) when *C*_{1}≥1 as the stationary mode does not exist in this range of *C*_{1}. In the interval *C*_{b}≤*C*_{1}<1, the preferred mode is determined numerically. For every pair of values in this interval, the minima of the stationary mode (3.15) and of the overstable mode (3.17) are compared and the smallest taken as the preferred mode. We find that as *C*_{1} increases beyond *C*_{b}, the stationary mode remains preferred, although the overstable mode can exist. When *C*_{1} reaches *C*_{p}, the Rayleigh numbers of the stationary and overstable modes become equal and when *C*_{1} increases further, the overstable mode becomes preferred for all *C*_{1}>*C*_{p}. In figure 1, the regions of existence and preference of the stationary and overstable modes are shown in a regime diagram in the (*p*,*C*_{1}) plane.

In complementing the results in [31], numerical computations were performed on the expressions (3.15) and (3.17) for Rayleigh numbers in the region *C*_{p}<*C*<1 where the exchange takes place from the stationary to the overstable mode (figures 2 and 3). We find that

(i)

*C*_{p}decreases as*p*increases,(ii) the critical Rayleigh number increases and the wavenumber decreases as

*C*_{1}increases when the mode is stationary,(iii) at the point where the mode changes from stationary to overstable, the wavenumber experiences a sizable positive jump (figure 3) and the overstable mode becomes preferred even though its wavelength is much smaller, and

(iv) for fixed

*p*, the critical Rayleigh and wavenumbers of the overstable mode decrease as*C*_{1}increases.

The regions II and III in figure 1 are associated with the coexistence of both stationary and overstable modes. When *C*_{1} increases to *C*_{b}, the overstable mode appears, but the energy required for exciting it is more than that required by the stationary mode. As *C*_{1} continues to increase further, the energy required for the stationary mode increases while that of the overstable mode decreases. When *C*_{1} reaches *C*_{p} the two modes require equal amounts of energy and they are equally excitable. Further increase in *C*_{1} reduces the energy required for exciting the overstable mode while that of the stationary mode increases allowing the overstable mode to be preferred for all *C*_{1}>*C*_{p}. (Such a situation is known to occur when the classical Benard layer is subject to constraints like rotation or magnetic fields [44].) The nonlinear analysis of §4 below examines the preference of the two modes for nonlinear instability (figure 7*a*).

These results are consistent with the asymptotic limits for small *C*_{1}
*C*_{1}>*C*_{p}. It is worth noting that *k*_{1c} decreases to 1, *R*_{1c}=*O*(1/*C*_{1}) and *p* is very large we find that *k*_{1c}, *R*_{1c},

When *C*_{1} increases beyond *C*_{p}, the critical Rayleigh number of the preferred overstable mode decreases while the corresponding wavelength increases reducing the influence of diffusion, although it always remains larger than that of the stationary mode which is not preferred for that range of *C*_{1}. On the other hand, when *C*_{1} increases from 0, both the critical Rayleigh number and wavelength of the preferred stationary mode increase with *C*_{1}. Although the wavelength increase is expected to reduce diffusion and promote instability, it is causing the opposite effect and is tending to stabilize the layer. This unconventional behaviour is explained using the energy balance equation in §3c below.

#### (ii) The LC model

The layer heated from below with LC model objective derivative can, in general, support both stationary and overstable modes. The overstable mode is obtained from the dispersion relation (3.13) with *R*_{1} over *k*^{2} yields
*p*,*C*_{1}, and the stationary mode provides the preferred mode of convection with
*C*_{1}=0, the preferred mode is that of the classical Benard layer. Computations of the expressions (3.26) showed that the critical Rayleigh number decreases while the wavenumber increases as *C*_{1} increases. This is consistent with the asymptotic results

#### (iii) The Christov model

The linear stability of this model has been discussed in [34,37,39]. We will therefore restrict our discussion to the part necessary for our nonlinear investigation in §4 below. The resulting dispersion relation (3.14) allows a stationary as well as an overstable mode for a layer heated from below. No instability is possible for a layer heated from above. The stationary mode is the same as the classical Benard layer mode with a minimum Rayleigh number and associated wavenumber given by
*C*_{1}.

The dispersion relation for the overstable mode gives

The preferred mode here is stationary for *C*_{1}<*C*_{b} as the overstable mode does not exist there. The overstable mode coexists with the stationary mode for *C*_{1}≥*C*_{b}, and the critical mode in this region must be determined by comparing the minimum Rayleigh numbers for the two modes by identifying which one is the smaller. This is done by identifying the curve along which the two Rayleigh numbers are equal
*k*_{1c} is given by (3.31), for *C*_{1}. Figure 4 depicts the regions of existence and preference of the two modes. It is clear that *C*_{p}>*C*_{b} so that although *C*_{b} determines the existence of the overstable mode, it is *C*_{p} that is relevant physically.

The critical Rayleigh and wavenumbers remain at the same values (3.29) for all *C*_{1}<*C*_{p}. For *C*_{1}>*C*_{p}, the critical Rayleigh number decreases steadily with *C*_{1} (figure 5) and so does the wavenumber, as pointed out in [39] (figure 5). When *C*_{1}(>*C*_{p}), we find

The change from stationary to overstable mode at *C*_{1}=*C*_{p} is accompanied by a positive jump in the wavenumber, as occurred for the SF model, so that the critical mode changes from a longer wavelength to a shorter one. Further increase in *C*_{1} is accompanied by a steady decrease in the wavenumber. As shorter wavelengths are associated with stronger influence of diffusion which tends to stabilize the system, it is of interest to identify the reason why the system is more unstable for shorter wavelengths when the preferred mode is overstable. This will be clarified when we consider the energy balance of the system below.

### (b) The stability of a layer heated from above

A layer heated from above has been shown to support linear instabilities for the SF and LC models but not for the Christov model. Here, we briefly summarize the properties of these instabilities obtained in [13,31,42] in preparation for the nonlinear investigations of §4 below. We also identify some new characteristics of the instabilities. In particular, we identify the behaviour of the critical Rayleigh number and its associated wavelength on *C*_{1} and *p* with a view to discussing the mechanism for exciting such instabilities.

#### (i) The SF model

A layer heated from above with the SF objective derivative can only support a stationary mode. The minimization of the expression for*R*_{1} with respect to *k* gives
*C*_{1}, we find that
*C*_{1} is small, the critical wavenumber is very large so that the cells are very narrow, resulting in strong diffusion and accordingly the magnitude of the critical Rayleigh number is large. As *C*_{1} increases, both critical Rayleigh and wavenumbers decrease (in magnitude) very rapidly at first and then more slowly. The influence of the M–C effect on this mode is to help the system release the strong influence of diffusion (by enlarging the convection cells) to make the system more unstable. Thus, diffusion plays its usual role on this system.

#### (ii) The LC model

The LC model layer heated from above can support overstable motions only. The critical mode obtained from (3.13) with *C*_{1} is small, all the critical mode parameters are very large
*C*_{1} are associated with the critical wavenumber, *k*_{1c}, of order 1 and the Rayleigh number and frequency both approaching zero
*p*, is unusual. We find
*p*→0, and
*ε*→0, where
*C*_{1}, the critical Rayleigh number increases indefinitely when the Prandtl number approaches zero or close to *α*_{0} when *p* approaches

When we compare the instabilities of the layer heated from above for the SF and LC models we find that they are very different, although the only difference between the two objective derivatives is the sign of *α*, i.e. a change of spin sign. In comparison with the Christov model, it is clear that the instabilities present in the layer heated from above are due to the presence of the spin in the objective derivative. The manner in which this occurs is explained by considering the energy balance of the system below.

### (c) Discussion

The above analysis of the linear stability has shown that the M–C effect introduces three main characteristics to the Benard layer:

the appearance of travelling waves and their preference when the M–C number,

*C*_{1}, exceeds a certain value,*C*_{p}, which depends on the Prandtl number,*p*;instabilities of layers heated from above;

the association of the preferred stationary mode of the layer heated from below with a critical Rayleigh number and wavelength that either both increase or both decrease with the increase in

*C*_{1}.

We will attempt to explain the manner in which the M–C effect gives rise to these characteristics by appealing to the energy balance and the heat flux of the system.

#### (i) Energy considerations

Consider the equations (2.9)–(2.12) with *ε*=0. Multiply (2.9) scalarly by **u**_{1}, integrate over a wavelength volume, *V* , and integrate the pressure term by parts to show that it vanishes and obtain
*θ*_{1} and integrate to get
*E*_{K}, *E*_{T} are the (dimensionless) kinetic and thermal energy densities, respectively,
*E*=*E*_{K}+*E*_{T}
*F*_{B}, viscous diffusion as represented by *F*_{ν}, and thermal diffusion as signified by *F*_{καC}. The terms with *F*_{B} and *F*_{ν} are similar to those obtainable for the classical Benard layer, but the term *F*_{καC} is explicitly influenced by the M–C effect. We have expressed *F*_{ακC} as a sum of two separate terms in order to bring out the differences between the Christov model on the one hand and the SF and LC models, which involve the spin, on the other hand. The term *F*_{κC} is the usual thermal diffusion term here modified by the presence of the M–C number, *C*_{1}, in the denominator. It reduces to its value for the classical Benard layer when *C*_{1}=0. This term is present for all three models. The M–C effect influences *F*_{κC} only for overstable modes, in which case it is reduced by any travelling wave, indicating that *C*_{1} tends to favour overstable modes. The term *F*_{αC} is entirely due to the presence of the spin term in the objective derivatives of the SF and LC models.

Noting that when *C*_{1}=0, *F*_{καC}=*π*^{2}*N*|*A*|^{2}, we see that the M–C influence has effectively replaced the thermal diffusivity, *κ*, by an effective M–C diffusivity *κ*_{e}
*κ* by *κ*_{e}, the energy equation takes the same form as that of the classical Bernard layer. The nonlinear dependence of *F*_{ακC} on *C*_{1},*N*,*p* and *ω*^{2} gives this term a level of versatility that is unusual in diffusive terms. It is this versatility that allows *F*_{ακC} to give rise to instabilities that are not possible for classical fluids. We will now address each of the instabilities identified here.

#### SF model layer heated from below

When *C*_{1}=0, the preferred mode is the stationary mode of the classical Benard layer and *F*_{αC} takes its positive classical value *π*^{2}*N*_{c}|*A*|^{2}. As *C*_{1} increases from 0, noting that *α*,*Ra*>0, the extra term *F*_{καC} becomes potent and takes the value
*F*_{καC} increases thereby enhancing thermal diffusion and consequently the critical Rayleigh number increases. On the other hand, the increase in *C*_{1} tends to promote the appearance of the travelling waves. For any non-zero frequency, *ω*, *F*_{κC} tends to decrease due to the increasing denominator. The spin term *F*_{αC}, given by (3.53), becomes relevant when *C*_{1} exceeds *C*_{b} because it will decrease from its stationary value rapidly at first and then more slowly (see solid line of figure 6 labelled I). The rapid decrease in *F*_{αC} added to the decrease in *F*_{κC} will result in decreasing thermal diffusion contribution to the overstable mode to the energy balance of the system. As *C*_{1} increases beyond *C*_{b}, thermal diffusion of the stationary mode increases, demanding more energy, while that of the overstable mode (as in (3.53)) is decreasing, requiring less energy until *C*_{1}=*C*_{p}, when the energy required by the overstable mode matches that of the stationary mode and further increase in *C*_{1} reduces the energy requirement for the overstable mode, hence remaining preferred for all *C*_{1}>*C*_{p}.

#### SF model layer heated from above

When *C*_{1}=0, the layer is stable. As *C*_{1} increases from 0, the term *F*_{κC} will not be affected while the term *F*_{καC}, noting that *α*>0,*Ra*<0, will take the value (3.53), which is negative for this value of *α*, and hence the M–C effect reverses the role of thermal diffusion of this stationary mode so that it provides energy to help drive the instability. Any non-zero frequency, *ω*, will reduce the term *F*_{καC}. This suggests that the M–C effect tends to reduce thermal diffusion of the stationary mode and tends to increase that of the overstable mode. Thus the preferred mode will remain stationary.

#### LC model layer heated from below

Here, *F*_{αC} vanishes at *C*_{1}=0. As *C*_{1} increases from 0, noting that *α*<0,*Ra*>0, *F*_{αC} is negative and its magnitude will increase with *C*_{1} thus promoting the instability of the stationary mode. The overstable mode is not possible here because the frequency *ω*^{2} will reduce the negative magnitude of *F*_{αC} and the mode will require more energy than the stationary mode.

#### LC layer heated from above

Here we use the same arguments as for the other layers, noting that *αR*_{1c}>0, to find that the M–C effect suppresses the stationary mode in favour of the overstable one. However, while an increase in the Prandtl number leads to a steady decrease in *F*_{καC}, the instability ceases to exist when *p* reaches *F*_{κC},*F*_{αC},*F*_{καC} were made, but no indication of why this occurs has been identified.

#### Christov model diffusion

Here, the term *F*_{αC} vanishes because of the absence of the spin. At *C*_{1}=0, *F*_{κC} takes its classical value when *ω*^{2}=0 and the preferred mode is that of the classical Benard layer. When the M–C effect is potent, the combination *C*_{1} increases from 0, the stationary mode preferred for the classical Benard layer remains preferred until a certain value, *C*_{p}, which depends on the Prandtl number, is reached. Here thermal diffusion has been reduced to the value that matches that of the stationary mode and the two modes are equally excitable. Further increase in *C*_{1} leads to the preference of the overstable mode for all *C*_{1}>*C*_{p}.

#### Influence of the M–C effect on the critical wavenumber

We have observed that the critical wavenumber, *k*_{1c}, of the preferred stationary modes increases (and the associated wavelength decreases) with *C*_{1} when the critical Rayleigh number, *R*_{1c}, is increasing and the reverse can also happen. This does not conform to the general understanding that shorter wavelengths are associated with more diffusion and hence less tendency towards instability. The preferred stationary mode of the SF model layer heated from below is associated with an increasing critical Rayleigh number and an increasing wavelength as *C*_{1} increases. When the expression (3.53) for *F*_{αC} is examined, we find that as *C*_{1} and *R*_{1c}, in the numerator, increase, the M–C effect acts to reduce the effective thermal diffusivity by decreasing the factor *k*_{1c}/*N*_{c} which can occur only for decreasing wavenumber *k*_{1c}, i.e. by increasing the wavelength. The same tendency occurs in the case of the preferred stationary mode of the LC layer heated from below. Here, the expression (3.53) is negative and an increase in the wavelength will increase the value of *k*_{1c}/*N*_{c} and reduce thermal diffusion, thus helping to enhance instability. *Thus the M–C effect acts to adjust the wavenumber in such a way as to reduce the effective thermal diffusion which requires a wavenumber increase in the case of the SF model and a wavenumber decrease in the case of the LC model layers heated from below*.

#### (ii) Heat transfer

The vertical heat flux obtained by the nonlinear analysis below to *O*(*ε*^{2}) is
*C*_{1}=0, and a uniform part which is entirely due to the M–C effect. This has the consequence that whether the layer is heated from below or from above, the M–C effect produces a mean vertical heat flux that is uniform across the layer and hence heat is transferred from hot regions to colder regions thus enhancing instability. We also note that the LC form of the spin (*α*<0) enhances the uniform part of the perturbation mean flux for the layer heated from above more than it does for the layer heated from below. The opposite effect occurs for the SF model.

The Nusselt number is defined by

In the case of the stationary mode, the terms proportional to *C*_{1} in the expressions for the Nusselt number, which are entirely due to the M–C effect, are positive and hence they enhance the heat transfer.

## 4. Nonlinear stability

When the Rayleigh number, *R*_{1}, is increased slightly above its critical value, *R*_{1c}, the neutrally unstable wave will acquire a real growth rate *γ*_{n} which can be expressed in the form
*c* refers to the quantity evaluated at the critical value. The quantities (∂*ω*_{1}/∂*k*_{1})_{c}, (∂*ω*_{1}/∂*R*_{1})_{c} are of order 1 and hence all quantities (*k*_{1}−*k*_{1}_{c}), (*R*_{1}−*R*_{1c}), *γ*_{n} are of the same order. For small values of (*R*_{1}−*R*_{1c}), we let
*u*_{g}=∂*ω*/∂*k* is the group velocity. The length scale here is different from that obtained for a similar problem for a Fourier fluid [46] because of the presence of the term involving the relaxation time. This change in length scale will have a profound effect on the evolution equation for the amplitude *A*.

The linearly unstable modes are either stationary or overstable. The stationary modes are standing waves while the overstable modes are travelling waves with phase speeds −*ω*/*k*. The frequency *ω* appears in the dispersion relation for the overstable waves only as *ω*^{2} and the minimization process yields an expression for *ω*/*k*| are equally preferred. As these two waves will be far apart after a long time, previously it was considered possible to neglect any interaction and consider the weakly nonlinear stability of each wave in isolation [47–49]. However, it has been found recently that although the waves may be far apart, it is of interest to consider their interaction. This will provide information on the influence of the modes on each other and how strong the interaction is (e.g. [50]).

The method of solution proceeds formally by expanding all the perturbation variables in the form
*f*=*f*(*x*,*z*,*ξ*,*τ*). We substitute the expansion (4.4) into the full perturbation equations (2.9)–(2.12) and the boundary conditions (2.17). We next equate the coefficients of *ε*^{j}, *j*=1,2,3,… in every equation and boundary condition on the two sides to obtain a set of problems and boundary conditions *j* which are solved successively. The method is straightforward and it may suffice here to say that the problem is closed by considering the problems for *j*=1,2,3. We should, however, note that
*z* is unaffected.

The leading order problem (*j*=1) is the linear problem discussed above. Here the solution has the form
*A*_{1},*A*_{2} are the amplitudes of the two waves to be determined. In the case of the stationary mode, the two equations become identical so that *A*_{1}=*A*_{2}=*A*_{s}. The other variables, which will be required for the evaluation of the coefficients of the evolution equations, can be obtained by adjustment of the expressions (3.8)–(3.12) to the form (4.6).

The problem for *j*=2, yields an equation similar to that for the leading order problem but with *k*,*ω* replaced by 2*k*,2*ω*, respectively. This is readily solved. The third order problem leads to an equation of the form
*L* is given by (3.3) for the SF and LC models, while (3.5) applies to the Christov model. The quantity *P*_{0} is a function of the variables of the two previous problems *j*=1,2. This equation is the non-homogeneous form of the equation for *j*=1. The existence of a solution then requires that the r.h.s. be orthogonal to the the solution, *j*=1. Thus,
*a*_{0},*a*_{1},*a*_{2},*a*_{3},*a*_{4} are functions of *k*_{1c},*ω*_{1c},*R*_{1c},*p*,*C*_{1} and the overbar denotes the complex conjugate. They are given in appendix A. Note that the coefficients of the second equation are the complex conjugates of the first equation, except the coefficient of ∂*A*_{i}/∂*ξ* (*i*=1,2) because of the appearance of the expression (i*k*), which is the same for all the waves as stated in the expressions for *E*_{1},*E*_{2} of equation (4.7).

The *ξ*-dependence can be removed from the equation by a linear transformation in *ξ* and *τ* [49]. We now analyse the evolution of the amplitudes as the slow time increases. We will not include the details of the analysis, and the reader is referred to excellent detailed descriptions of the basic aspects of the analysis in [51], chapter 1, [52], chapter 5, [53], chapter 7.

Next we let

The coupled amplitude equations (4.13) and (4.14) govern the evolution of the amplitudes of the two waves travelling in opposite directions. This is indicated in the phase equations (4.15), (4.16) by the different signs on the r.h.s. In the case of the stationary mode, *ω*=0, and the coefficients are real. The two equations become identical and give the Landau–Stuart equation

The coupled amplitude equations (4.13)–(4.14) cannot be readily solved analytically. However, the evolution of the amplitudes can be described by consideration of their equilibrium points and the computation of the coefficients *d*,ℓ_{1},ℓ_{2}. It is convenient to introduce the two quantities
*r*_{1}=*r*_{2}=0. This is referred to as the first equilibrium point when another equilibrium point occurs at
*d*_{r} (i.e. it behaves like *i*=1,2) so that the solutions are growing if *d*_{r}>0 (or *d*_{s}>0) or decaying if *d*_{r}<0 (or *d*_{s}<0). At the second equilibrium point, when it exists, the behaviour of the solution depends on the two quantities

Now the computations have shown that except for the SF layer heated from above, *d*_{r}>0 so that the second equilibrium point (4.20) exists only if *S*<0. This makes *s*_{2}<0. If *D*>0, then both *s*_{1},*s*_{2} are negative, in which case the point is stable. In this situation, the amplitude increases from its value at the bifurcation point to *r*_{0} at the second equilibrium point where it remains for all further times. The system is here said to be supercritically stable. This is referred to as region III (defined below) in figure 7. If, however, *D*<0, then *s*_{1}>0,*s*_{2}<0, in which case the point is a saddle point. In this case, the amplitude approaches *r*_{0} along a hyperbolic path in the phase space but maintains its increase past it (i.e. the amplitude as a function of time experiences an inflection point). Here the system is supercritically unstable despite the occurrence of a second equilibrium point. This case is referred to as region II (defined below) in figure 7.

If *S* is positive, then there is no second equilibrium point and the amplitudes will grow with time giving rise to supercritical instability. Here the sign of *D* is not relevant. This case is referred to as region I (defined below) in figure 7.

The results for the stationary mode can be deduced from those of the overstable mode by replacing *d*_{r},*S* by *d*_{s},ℓ_{s}, respectively, while *D* is not relevant.

Considering the signs of the quantities *d*_{r},*D*,*S* together with the observation that the two equations remain the same if we exchange the two amplitudes, we find that there are, in general, four possible situations that can arise. We will refer to them as regions I–IV in the parameter space (*p*,*C*_{1}) represented in figure 7 and in the discussion below.

(i) Region I:

*d*_{r},*S*>0. Here both linear and nonlinear contributions to the amplitudes are positive and there is only one equilibrium point, which is unstable. The amplitudes grow with time.(ii) Region II:

*d*_{r}>0,*S*<0,*D*<0. Here there are two equilibrium points. The bifurcation point is unstable while the second equilibrium point is a saddle point. The amplitudes approach an inflection point at*r*=*r*_{0}but continue to increase monotonically with time. The waves are supercritically unstable.(iii) Region III:

*d*_{r}>0,*S*<0,*D*>0 . The second equilibrium point at*r*_{1}=*r*_{2}=*r*_{0}is stable. The amplitudes increase from their value at the bifurcation point to the value*r*_{0}where they remain as a second steady solution. The waves are referred to as supercritically stable.(iv) Region IV:

*d*_{r}<0,*S*>0. The amplitude of the stationary mode will grow nonlinearly but decay linearly. As time increases, the two effects will balance and the waves will acquire the amplitude*r*_{0}and remain steady. Here, the wave suffers subcritical instability because it will be unstable for*R*_{1}<*R*_{1c}(because*d*_{r}<0).

In general, there is another possibility in which both *d*_{r},*S* take negative values indicating that both linear and nonlinear terms on the r.h.s. of equations (4.13) and (4.14) are negative. This situation corresponds to a nonlinear stabilization of a preferred mode. However, this case does not occur for any of the five layers considered here.

We should mention here that the equations (4.13) and (4.14) can possess two other equilibrium points at *r*_{2}=0,*r*_{1}=*r*_{0} and *r*_{1}=0,*r*_{2}=*r*_{0}, with *r*_{0}=−*d*_{r}/ℓ_{1r}. These two cases refer to the situations when only one overstable mode is present. They have identical behaviour because the coefficients of the evolution equation are identical. Computations showed that stability of a single overstable mode is similar to that of the two interacting waves in the case of the Christov model but the LC layer heated from above is supercritically unstable for all *C*_{1} and *C*_{1}.

The coefficients *d*_{r},*D*,*S* have been computed in the parameter space (*p*,*C*_{1}) for the critical modes of the different layers. The results are summarized in nonlinear stability regime diagrams in figure 7, except the SF model layer heated from above. We make the following observations:

The preferred stationary mode of the SF layer heated from above is independent of the Prandtl number and

*d*_{r}is positive for all*C*_{1}. The range of*C*_{1}is divided into two parts. In the interval 0<*C*_{1}≤13.01, ℓ_{s}is negative, and a stable second equilibrium point exists. The solution here is supercritically stable. In the range*C*_{1}>13.01, ℓ_{s}is positive and the layer is supercritically unstable.Supercritical instability (regions I and II) is present for moderate values of

*C*_{1}in the case of a SF layer heated from below and for large values of*C*_{1}in the cases of a LC layer heated from below and the SF layer heated from above. For the LC layer heated from above, it is present for a bounded region as shown in figure 7*d*. Supercritical stability does not occur for the Christov layer.Supercritical stability (region III) is present in the SF layer heated from below, the LC layer heated from below and the LC layer heated from above. However, the regions of the parameter space where they occur differ from layer to layer (figure 7).

Subcritical instability occurs for the LC layer heated from below in a strip of the (

*p*,*C*_{1}) plane adjacent to the axes (region IV in figure 7c).

## 5. Conclusion

The linear instabilities identified in [13,31,34] referred to as SF, LC and Christov models, respectively, for fluids obeying the M–C heat flux relationship, have been investigated with a view to identifying the mechanisms responsible for the appearance of the instabilities absent in Fourier fluids and those responsible for the choice of type of instability (stationary or travelling) that occurs. The instabilities of all the models depend on two fluid parameters: the Prandtl number, *p*, and the M–C number, *C*_{1}. The linear stability regime diagram established for each layer displays the regions where the different types of mode can exist and where they are preferred. Two situations, namely an LC model layer heated from below and a SF model layer heated from above, can support stationary modes only. For the other cases, the plane (*p*,*C*_{1}) is divided into regions where stationary and overstable modes can exist, with both stationary and overstable modes coexisting in some regions of the parameter plane. The stationary mode always exists for smaller values of *C*_{1} while the overstable mode occupies the region for larger *C*_{1}. Both stationary and overstable modes are possible for moderate values of *C*_{1} with the boundary between them varying with the Prandtl number. However, the SF overstable mode cannot exist for *C*_{1}=0 apply.

The mechanism driving the overstable modes in layers heated from below and destabilizing layers heated from above introduced by the M–C effect has been investigated using energy considerations. It is found that the M–C effect tends to enhance travelling waves and interacts with fluid spin to influence the effect of thermal diffusion on the wavenumber and frequency of the possible waves. This action manifests itself in the expression for the thermal diffusion contribution to the energy balance of the system. This expression is found to reduce to that obtained in the absence of the M–C effect if we define an *effective thermal diffusivity* *κ*_{e} as given by (3.53). Now *κ*_{e}, in general, depends on the M–C number, the wavenumber, the frequency, the spin and the Prandtl number. The dependence of the expression on these quantities is such that it is reduced in value for certain wavenumbers and frequencies, and sometimes reversed (from positive to negative) so that it supports instability. As discussed in detail in §3c above, it has been shown in particular that (i) the nonlinear dependence on the combination *αRa*) determines the existence of instabilities in layers heated from above, and (iii) the dependence of *κ*_{e} on the wavenumber explains why preferred stationary modes with increasing Rayleigh number are associated with increasing wavelengths and those with decreasing critical Rayleigh number are associated with decreasing wavelengths.

The M–C flux is found to give rise to contributions to the mean vertical flux and Nusselt number that are uniform across the layer thereby promoting the transfer of heat across the layer even when the layer is heated from above.

A formal nonlinear analysis of the least stable neutral modes of the linear theory taking full account of the interaction of the two equally preferred overstable modes travelling in opposite directions yielded a pair of coupled evolution equations for the amplitudes of the two overstable modes. The two equations coincide when the mode is stationary. The amplitudes evolve on a length scale that is long compared to that occurring for a Fourier fluid. Consequently, the two equations take the form of a Landau–Stuart equation with a term representing the interaction between the two waves. This resemblance of the equations to the well-known Landau–Stuart equation has made it possible to infer the general stability properties of the waves using the signs of their coefficients. The results are summarized in §4 and presented in figure 7 above. In addition, it can be mentioned here that

The Christov layer, which does not possess the spin factor (i.e.

*α*=0), is supercritically unstable for the whole parameter space (*p*,*C*_{1}). When*C*_{1}is large, a second equilibrium point exists, but as a saddle point so that the amplitude suffers a point of inflection, and remains supercritically unstable (region II).The LC layer heated from below is supercritically stable for large

*C*_{1}and subcritically unstable for small values of*C*_{1}whatever the value of*p*, while the LC layer heated from above can suffer supercritical stability and supercritical instability depending on the relative values of*C*_{1}and*p*.None of the five layers (SF layers heated from below or from above, LC layers heated from below or from above and the Christov layer heated from below), considered above can give rise to all four types of possible behaviour (supercritical stability, supercritical instability, subcritical instability and stability) in the parameter space.

None of the preferred modes of the linear theory is stabilized by the nonlinear interactions.

## Competing interests

I declare I have no competing interests.

## Funding

This work did not involve any formal funding support.

## Acknowledgements

I wish to thank the anonymous referees for some helpful suggestions and Sultan Qaboos University for the use of general academic facilities.

## Appendix A

Here, we give the expressions for the coefficients of the evolution equations (4.13) and (4.14). The coefficient *a*_{1} is not given because it is not required for the calculations.

For the SF (

- Received October 20, 2016.
- Accepted April 13, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.