Maxwell–Cattaneo double-diffusive convection: limiting cases

Abstract Double-diffusive convection, in which a fluid is acted upon by two fields (such as temperature and salinity) that affect the density, has been widely studied in areas as diverse as the oceans and stellar atmospheres. Assuming classical Fickian diffusion for both heat and salt, the evolution of temperature and salinity are governed by parabolic advection–diffusion equations. In reality, there are small extra terms in these equations that render the equations hyperbolic (the Maxwell–Cattaneo effect). Although these corrections are nominally small, they represent a singular perturbation and hence can lead to significant effects when the underlying differences of salinity and temperature are large. In this paper, we investigate the linear stability of a double-diffusive fluid layer and show that amending Fick's law for the temperature, or the salinity, alone can lead to new modes of oscillation and to very large changes in the preferred wavelength of oscillatory convection at onset. In particular, the salt finger regime of classical double diffusion is here replaced by Maxwell–Cattaneo oscillations when the salt concentration is very high. The more complicated case when both laws are amended is left to a future paper, now in preparation.

involving magnetic fields, double diffusion, etc., use the classical Fick's law to describe the relation between quantities such as, for example, temperature and heat flux. For temperature, this relation (known as Fourier's law in this specific case) takes the form q = −K∇T, where q is the heat flux, T is the temperature and K is the thermal conductivity. The Fourier law predicts an instantaneous response of temperature to the heat flux gradient, leading to a parabolic diffusion equation for the temperature field. The instantaneous response at all points implied by this equation cannot be exactly correct as information must travel at a finite speed. This weakness was recognised by Maxwell (1867) in his study of the theory of gases, who proposed a modified equation incorporating a finite relaxation time. Cattaneo (1948) proposed a similar relation for solids, which was developed further by Oldroyd (1950). Other important contributions were made later, by, for example, Fox (1969) and Carrassi & Morro (1972).
The idea of a finite relaxation time is incorporated into the M-C relation in the temperature equation between the heat flux q and the temperature T, which then takes the form in which τ T is the relaxation time. The operator D/Dt here denotes a generalised Lagrangian time derivative, which should be chosen to give expressions that do not depend on the frame of observation; a specific form of this generalised derivative is discussed in § 2.1. The importance of the thermal relaxation term is typically expressed via the dimensionless M-C coefficient C T , which is defined as the ratio of the thermal relaxation time to twice the thermal diffusion time; i.e. C T = τ T K/(2ρc p d 2 ) = τ T κ/2d 2 , where ρ is the density, c p the specific heat at constant pressure, d a representative length scale and κ the thermal diffusivity. (The factor of two in the denominator of the expression for C T does not seem to be particularly helpful. However, it is the definition of C T used previously in the literature, which we therefore choose to retain for consistency.) Thus the classical Fourier law has C T = 0. The introduction of a finite relaxation time changes the fundamental nature of the parabolic heat equation of Fourier fluids, in which heat diffuses with infinite speed, to a hyperbolic heat equation with a solution in the form of a heat wave that propagates with finite speed (Joseph & Preziosi 1989;Straughan 2011a).
The M-C heat transport effect has been studied in a wide variety of different physical contexts: for example, in solids (Barletta & Zanchini 1997), in fluids (Lebon & Cloot 1984;Straughan & Franchi 1984;Straughan 2009Straughan , 2010Stranges, Khayat & Albaalbaki 2013;Bissell 2015;Stranges, Khayat & deBruyn 2016;Eltayeb 2017), in porous media (Straughan 2013;Haddad 2014), in nanofluids and nanomaterials (Jou, Sellitto & Alvarez 2011;Lebon et al. 2011), in liquid helium (Liepmann & Laguna 1984;Donnelly 2009), in biological tissues (Dai et al. 2008;Tung et al. 2009) and, in the context of magnetoconvection, in stellar interiors and the solar photosphere (Bissell 2016;Eltayeb, Hughes & Proctor 2020). The potential significance of the M-C effect depends on a number of factors, through the definition of the coefficient C T . In gases, the relaxation time τ T can be as small as picoseconds; C T will then only assume O(1) values over very small scales, as can occur for heat transport in nanostructures. Nonetheless, even if C T 1, the M-C effect can still be important if the thermal driving (measured by the Rayleigh number, for example) is extremely high; this is often the case in astrophysical settings. In biological matter, the situation can be quite different, since τ T can be of the order of a second or greater; the M-C effect can then be important on everyday length scales.
In this paper we study the consequences of including the M-C effects on the onset of double-diffusive convection, in which two quantities affect the density of a fluid, but diffuse at different rates. There are now M-C effects to be considered for both temperature and salinity. The most widely studied example of double diffusion is thermohaline convection, in which the competing ingredients are heat and salt, with the diffusion of heat greatly exceeding that of salt. Here we shall stick to the terminology of 'heat' and 'salt', although the equations are relevant in a much wider context. On including the M-C salinity effect, the modified salinity evolution equation can be obtained by analogy with the temperature equation, so can be written as where S is the salt concentration, κ S the saline diffusivity, τ S the relaxation time for salinity and q S the salt flux. We note that the relation between salt flux and salt concentration differs from that between heat and temperature; hence the appearance of K in (1.1) but κ S in (1.2). For later use we define the M-C coefficient for salt as C S = τ S κ S /2d 2 . Double-diffusive convection has been widely studied for many decades in a variety of geophysical, astrophysical and engineering contexts (see, for example, the reviews by Turner (1974), Huppert & Sparks (1984), Turner (1985), Schmitt (1994), Garaud (2018) and the monograph by Radko (2013)), but principally employing the Fickian law for the evolution of both diffusing ingredients. Herrera & Falcón (1995), via a fluid parcel argument, considered the new physics introduced by the thermal M-C effect (but with C S = 0) into double-diffusive convection, motivated by heat transport in neutron stars and compact X-ray sources, where the competing gradients are of temperature and helium concentration. Straughan (2011b) looked at this problem in more depth, but again restricted attention to the case of C S = 0 and for the case of a porous medium rather than a viscous fluid.
Motivated by geophysical and astrophysical considerations, we concentrate here on the case where the M-C coefficients are very small. The modified equations then represent singular perturbations in the time domain, and hence, even when C T , C S 1, new mechanisms for oscillatory instability can arise, provided that the initial gradients of temperature and salinity are very large. The goal of this paper is to investigate the effects of the new terms on the onset of instability in the regime of C T , C S 1. The dependence of the results on C T and C S is very rich, so in this first paper we confine ourselves to the two special cases where one of the M-C coefficients is zero. The general case turns out to merit a second paper, which is in preparation.
The plan of the paper is as follows. In § 2.1 we demonstrate how the M-C effects can be incorporated into the equations describing the evolution of temperature and salinity; in § 2.2 we give the full set of governing equations for M-C double diffusion; § 2.3 presents the linearised stability problem for a basic state with linear gradients of temperature and salinity, with simple boundary conditions. In the two subsequent sections we investigate the special cases of C S = 0 ( § 3) and C T = 0 ( § 4); for each of these we investigate the onset of linear instability for all combinations of the temperature and salinity gradients, concentrating on small values of the M-C coefficients and large values of the gradients. The significance of the results is assessed in the concluding section ( § 5).

The M-C effect in a moving frame
To derive the governing equations, it is necessary to consider the modifications of the M-C effect to be expected in a fluid moving with velocity u. A number of possible formulations of the nonlinear terms have been proposed (e.g. Lebon & Cloot (1984) commenting on Straughan & Franchi (1984); Christov (2009)). Here we follow the formulation of Christov (2009), who proposed the following frame-invariant equation for the evolution of the heat flux: It is convenient to write this equation as since, on taking the divergence, we obtain where Q T = ∇ · q T and where we have assumed that K is constant. When ∇ · u = 0, as we assume in this paper, the left-hand side of (2.3) can be written as the usual Lagrangian derivative of Q T . Thus the formulation of Christov (2009) is particularly convenient because it is reducible, in that the governing equation can be written in terms of Q T , with q T no longer appearing explicitly.
Although the form of the governing equation for q S has not been addressed previously, it is clear that, at the very least, the diffusion equation is inconsistent with relativity theory and so extra terms are essential. A form analogous to the temperature equation is to be expected, since it should be invariant under space reflection (hence no odd space derivatives) and with the time derivative term as the leading-order expression in an expansion for small τ S . We thus write the evolution equation for Q S as (2.4) 2.2. Governing equations We consider a horizontal layer of an incompressible (Boussinesq) viscous M-C fluid, initially at rest, contained between two planes at z = 0 (bottom) and z = dπ (top). The scaling here with π is helpful in that all factors of π are eliminated from the governing equations. The fluid has kinematic viscosity ν, and thermal and salt diffusivities κ, κ S . The density depends linearly on two components that diffuse at different rates. By analogy with classical thermal convection, we shall denote the faster diffusing component by T (temperature) and the slower by S (salinity). In equilibrium, the fluid is at rest, with temperature and salinity differences across the layer of ΔT and ΔS. The basic state temperature and salinity,T andS, are thus given bȳ where T 0 and S 0 are representative values of temperature and salinity. For the perturbed state, with velocity u = (u, v, w), we express the temperature and salinity by T =T +T, S =S +Ŝ. The density ρ of the fluid obeys a linear relation of the form where the Rayleigh number Ra, the salt Rayleigh number Rs, the Prandtl number σ , and the diffusivity ratio τ are defined by (2.12) With the Rayleigh numbers so defined, positive (negative) Ra is thermally destabilising (stabilising), whereas positive (negative) Rs is solutally stabilising (destabilising).

Linearisation and stability considerations
In this paper, we address the linear stability of the basic state given by (2.5), subject to the standard boundary conditions in which the horizontal boundaries are impermeable and stress-free, and on which the temperature and salinity are fixed. Thus noting that z is now dimensionless. We assume periodicity in the horizontal directions. In general, we may decompose the solenoidal velocity as (2.14) The linearised form of (2.7) shows, however, that T decays for all parameter values, and thus only P is of relevance. Following the usual approach to the classical double-diffusive stability problem, we seek solutions to the linearised versions of (2.7)-(2.11) of the form where the planform function f (x, y) satisfies with ∇ 2 H being the horizontal Laplacian. For the classical problem, with no M-C effects, it is easily shown that the fundamental mode (i.e. m = 1) is the most readily destabilised.
Here we shall also restrict attention to the m = 1 mode, but will discuss this assumption in § 5, in the light of the results.
On substitution from (2.15) into the linearised forms of (2.7)-(2.11), we obtain, after some algebraic manipulation, the following quintic dispersion relation for the growth rate s: a 5 s 5 + a 4 s 4 + a 3 s 3 + a 2 s 2 + a 1 s + a 0 = 0, (2.17) where and where β 2 = k 2 + 1. The third-order system of classical thermohaline convection is recovered by setting C T = C S = 0. The third-order system governing M-C Rayleigh-Bénard convection, studied by Stranges et al. (2013) and Bissell (2015), is recovered by setting C S = Rs = τ = 0. Note that the system governed by (2.17) and (2.18) possesses the following symmetry: To provide a natural link to the classical thermohaline problem, we shall restrict attention to τ < 1; the case of τ > 1 can be recovered through the transformation (2.19).
In this paper, we shall concentrate on determining the conditions for the onset of instability; this may occur either as a direct mode (steady convection), in which case the growth rate s passes through zero, or as an oscillatory mode, in which case, at onset, s = ±iω, with ω ∈ R + . It is traditional in studies of double-diffusive convection to treat Ra as the bifurcation parameter, although, mathematically, there is nothing to favour Ra over Rs. For comparison with the existing literature, we shall maintain this tradition here. We shall refer to the mode that first becomes unstable as Ra is increased as the preferred or favoured mode.
At the onset of steady convection, the coefficient a 0 = 0. Since a 0 has no dependence on either C T or C S , M-C effects therefore have no influence on the onset of steady convection, as is to be expected from the form of the flux equations (2.9) and (2.11). The value of Ra at the onset of steady convection is given by (2.20) The critical value of Ra (s) , which we shall denote by Ra (s) c , is given by the minimum value of Ra (s) over all wavenumbers: Ra (s) is minimised when k 2 = k 2 c = 1/2, thus giving (2.21) We note that for the limiting cases considered here, in which either C T or C S are zero, the coefficient a 5 is zero and the growth rate is governed by a fourth-order equation. Setting s = ±iω, with ω ∈ R + , leads to the coupled equations a 4 ω 4 − a 2 ω 2 + a 0 = 0, a 3 ω 2 − a 1 = 0. (2.22) Since Ra and Rs occur linearly in coefficients a 0 , a 1 , a 2 , a 3 -and are absent in a 4 -then we can, for example, readily combine equations (2.22) either to derive a quadratic equation for ω 2 that does not involve Ra, or eliminate ω 2 to derive an expression quadratic in Ra and Rs. Both approaches turn out to be useful; we shall consider the specific forms of these expressions in the following two sections.
In the classical double-diffusive problem, oscillatory instability can occur only in the first quadrant of the (Rs, Ra) plane (i.e. Rs and Ra both positive). In the absence of M-C effects, the coefficient a 4 = 0; hence, from (2.22), the marginal value of Ra for oscillatory motions (i.e. when s = ±iω, ω ∈ R + ) is given by the expression a 0 a 3 = a 1 a 2 , which becomes (2.23) Using (2.22) and (2.23), the necessary additional constraint of ω 2 > 0 translates to the condition It is straightforward to show that when there is a pair of purely imaginary solutions for s, the third (real) root for s is negative. Thus, since β 6 /k 2 is minimised when k 2 = 1/2, we can see from (2.24) that oscillatory motions are preferred at onset provided that (2.25) from (2.23), the critical Rayleigh number is then given by (2.26) The overall stability boundary for the classical problem is sketched in figure 1. The regime of steady convection in the third quadrant is often referred to as the salt fingering regime; that in the first quadrant in which oscillatory modes are preferred as the diffusive regime. 3. The case of C S = 0 3.1. Stability boundaries As already noted, in this case a 5 = 0, and so the dispersion relation (2.17) reduces to a quartic equation with If we eliminate Ra between equations (2.22) then the frequency ω on the oscillatory boundary is determined by the following quadratic equation for ω 2 : Conversely, eliminating ω 2 gives the following quadratic expression for Ra on the oscillatory boundary (again provided that ω 2 > 0): (3.5c) , then, in comparison with the classical double-diffusive problem, there are only small changes to the critical value of Ra and the corresponding critical k 2 . The question of interest therefore is to ask where small C T makes a fundamental difference. To this end, it is instructive to consider the regimes increasing n thus provides a means of delineating regimes of increasing Rs. As described below, distinct regimes can be identified for 0 < n ≤ 1, 1 < n < 2, n = 2 and n > 2. We have already noted the very different salt fingering and diffusive regimes of the classical problem, shown in figure 1. With the inclusion of M-C effects, there are further significant differences between the first and third quadrants of the (Rs, Ra) plane. The exposition is therefore clearest if we consider these quadrants separately, examining in each the different regimes for the scaling exponent n in (3.6).

First quadrant: Ra
For O(1) values of Rs, the problem is essentially that of classical double-diffusive convection, with oscillatory instability favoured provided that inequality (2.24) is satisfied. Indeed, as we shall see, once this inequality is satisfied then oscillatory motions are preferred for all n. Qualitative changes from the classical problem first arise when n = 1. For classical double diffusion, the critical wavenumber for oscillatory convection (and indeed steady convection also) is given by k 2 = 1/2. To explore how the M-C effect can influence this picture, we consider the regime C T 1, If we write Ra = Ra 0 + C T Ra 1 , the leading-order terms give We note that expression (3.8), which crops up quite frequently, is the wavenumberindependent component of (2.26), the expression for Ra at the onset of oscillatory convection in the classical problem. It can be checked from expression (2.22) that (3.8) leads to ω 2 > 0 only for Rs > 0; thus these considerations apply only to the first quadrant in the (Rs, Ra) plane. At this order, the steady branch is given by Ra = Rs /τ and so, for τ < 1, it lies above the oscillatory branch. We note from (3.8) that at leading order there is no wavenumber dependence of the critical Rayleigh number. This concept, which is similar to the description of the onset of oscillatory magnetoconvection in a strong magnetic field (see, for example, Weiss & Proctor 2014), turns out to be prevalent throughout this entire M-C problem. To determine the critical wavenumber, it is necessary to proceed to the next order. After some algebra, we obtain (3.9) Stationary points (d Ra 1 /dk 2 = 0) are therefore given by the expression which is a quintic polynomial in k 2 . Figure 2(a) shows the evolution of the real parts of the roots of this quintic polynomial as Rs is increased. With Rs = 0, the only real positive root is at k 2 = 1/2 (the classical case). As Rs increases, the existing mode moves to higher k 2 and, at a critical value of Rs , two equal stationary points appear, which then separate with a further increase in Rs ; k 2 remains negative (and thus unphysical) for the other two roots. Figure 2(b) plots the oscillatory stability boundary with no approximations, together with the zeroth-order approximation (3.8), and the first-order correction (3.9); the latter is almost indistinguishable from the boundary of the full system. The three stationary points correspond to the three positive roots for k 2 at Rs = 15 in figure 2(a). Thus the M-C influence is felt for Rs = O(C T −1 ) by the emergence of two new stationary points in the oscillatory stability boundary. As demonstrated in figure 2, the preferred mode is the continuation of the classical mode, moved to higher wavenumbers.  Figure 2. (a) The real parts of the roots of (3.10) for k 2 as a function of Rs for 0 < n ≤ 1; blue lines denote real roots, the red line denotes the real part of a conjugate pair; σ = 1, τ = 0.1. (b) Oscillatory stability boundary (blue solid line), together with the zeroth-order asymptotic expression (3.8) (purple dashed line) and the first-order correction (3.9) (purple squares); C T = 10 −3 , C S = 0, Rs = 1.5 × 10 4 (corresponding to Rs = 15).

n = 2
In § 3.2.1, we saw how, for n = 1, a new pair of stationary points emerges, leading to two minima and one maximum in the Ra (o) versus k 2 curve, all with k 2 = O(1). As Rs is increased, (i.e. n increased), these stationary points separate; at n = 2, they attain distinct asymptotic scalings, with the two minima of Ra (o) having k 2 = O(C T ) and k 2 = O(C T −1 ), and with the maximum having k 2 = O(C T −1/2 ). This is illustrated clearly in figure 3. The value of Ra (o) for the minimum at k 2 = O(C T ) is readily attained analytically, at least to leading order. Noting that β 2 ≈ 1 for small k 2 , the coefficients of (3.4) become, at leading order, The roots have Ra = O(C T −2 ) and Ra = O(C T −3 ), with only the former being admissible (ω 2 > 0). Thus, with k 2 = O(C T ), We see straightaway that steady modes cannot be preferred when n = 2. The wavenumber dependence of Ra (o) has again dropped out at this order, possibly unsurprisingly given the very flat nature of this minimum, as shown in figure 3. If needed, this could be retrieved by retaining the next-order terms. However, numerical calculations suggest that the minimum of Ra (o) at large k 2 is always the smaller, but that the two minima do become arbitrarily close for small values of σ and τ (cf. figures 3a and 3b). In general, it is not possible to obtain an analytic expression for the minimum value of Ra (o) with k 2 = O(C T −1 ). However, we can make progress in the case of C T σ, τ 1, through a distinguished limiting process, in which we first let C T → 0 and then let σ and τ become small. To explore the C T → 0 limit, we introduce the scalings The coefficients of expression (3.4) are almost unchanged, except with Rs replacing Rs, etc., and with C T = 1; the substantive difference is that β 2 = k 2 at this order. We then let σ and τ become small, with the wavenumber scaled ask 2 = σk 2 . The onset of oscillatory instability is then given by (3.14) Hence Ra (o) is minimised when Thus, to this degree of approximation, one cannot distinguish between the minimum at small k 2 (expression (3.12)) and that at large k 2 (expression (3.15)).

n > 2
As for the case of n = 2, there are two minima in Ra as a function of k 2 : the minimum at small wavenumber has k 2 = O(C n−1 T ), that at large wavenumber has k 2 = O(C T −n/2 ).
When k 2 = O(C n−1 T ), the frequency equation (3.2), at leading order, becomes There is one admissible root, given by with the oscillatory stability boundary then given by ( 3.18) For the minimum at large k 2 , with k 2 = O(C T −n/2 ), the frequency equation, at leading order, becomes There is one admissible root, given by (3.20) From the relation ω 2 = a 1 /a 3 we then obtain, to leading order, (3.21) From (3.21), Ra (o) is minimised when From (3.22), we can see that the minimum at large k 2 has Ra = O(C T −(1+n/2) ), whereas from (3.18), that at small k 2 has Ra = O(C T −n ). Thus for n > 2 (and in contrast to the case of n = 2), the two minima have distinct asymptotic scalings: the lower minimum occurs at large k 2 , and is given by (3.22). Furthermore, since the onset of steady convection is given by Ra (s) ≈ Rs/τ = O(C T −n ), oscillatory modes are always preferred for n > 2.

Overall stability boundary
It is important to put together the above ideas in order to determine the critical Rayleigh number and the associated critical wavenumber for a wide range of Rs (i.e. for a range of n). Figure 4(a) shows Ra c for the range 10 −5 ≤ Rs ≤ 10 15 (corresponding to −5/3 ≤ n ≤ 5) for fixed values of C T = 10 −3 , σ = 1, τ = 0.1, together with the asymptotic results for n > 2. Figure 4(b) shows k 2 c at the onset of instability. For small enough Rs, the onset is always steady, with k 2 c = 1/2. As Rs is increased, the preferred mode becomes oscillatory, initially with little change in the wavenumber of the critical mode. However, for n 1, the critical wavenumber increases with increasing Rs: for the n = 2 regime, Ra c and k 2 c both depend linearly on Rs (e.g. expression (3.15) for small σ and τ ), whereas for n > 2, Ra c and k 2 c both vary as Rs 1/2 (expression (3.22)). Also shown in figure 4(a) is the onset of oscillatory instability for the classical problem. In terms of the critical Rayleigh number (but, as noted above, not the preferred wavenumber), it can be seen that the M-C effect only comes into play strongly for n 2, where it leads to a clear enhanced destabilisation of the oscillatory modes.
3.3. Third quadrant: Ra < 0, Rs < 0 3.3.1. n ≤ 2 As discussed in § 2.3, for classical double-diffusive convection the onset of instability in the third quadrant is always via a steady bifurcation (salt fingers). Therefore, in comparison with the first quadrant (for which the classical problem has an oscillatory instability), larger values of |Rs| are needed before there are substantive variations from classical double-diffusive convection: M-C effects are not felt until n = 2 (i.e. Rs = O(C T −2 )).
For n = 2, the minimum in the Ra (o) curve as a function of k 2 is located in the regime of k 2 = O(C T −1/2 ) (which in the first quadrant is the location of a local maximum). The behaviour in this regime is captured by the orderings k 2 ∼ β 2 = O(C T −1/2 ), and At leading order, (3.4), governing the onset of oscillatory convection, becomes To this order there is no wavenumber dependence of Ra and thus to determine the dependence on k 2 , we need to retain terms of the next order; we thus write Ra = Ra 0 + C T 1/2 Ra 1 . The leading-order terms determine Ra 0 via (3.23); Ra 1 is determined at the next order through the expression If F, A and B are of the same sign, then Ra 1 is minimised wheñ (3.26) Figure 5 shows the stability boundaries for oscillatory and steady modes as functions of k 2 , together with the associated values of ω 2 ; also shown are the asymptotic results (3.23) and (3.24). It can be seen that oscillatory instability occurs within a closed loop in the (Ra, k 2 ) plane. As |Rs| decreases (i.e. Rs becomes less negative), the loop shrinks, eventually disappearing altogether; equivalently, and maybe what is a more natural way of thinking, in terms of an experiment, the loop appears at a critical value of Rs as |Rs| is increased. An important feature is that this spontaneous appearance may occur either above or below the steady branch. To determine the least negative value of Rs for which oscillations are possible thus requires determination of the critical value of Rs for which the loop first appears, calculation of Ra at that value, and a comparison with the value of Ra on the steady branch at the same value of Rs.
The loop of oscillatory instability appears when the two roots of (3.23) coincide. After a little algebra, it can be shown that the critical value of Rs ( Rs , say) is then given by For comparison, we must also calculate the values of Rs at which the stability boundaries for the steady and oscillatory modes intersect ( Rs so , say). At leading order, the steady boundary (2.21) is approximated by Ra = Rs /τ . Substituting into (3.23), the leading-order expression for the oscillatory stability boundary, gives . (3.28) The comparison to be made is between Ra evaluated at Rs , when the loop first appears, and Rs /τ , which is the value of Ra on the steady branch. If Ra( Rs ) < Rs /τ (both negative recall), then the loop of oscillations appears beneath the steady branch (i.e. at higher | Ra|); the critical value of Rs for the onset of oscillations is then Rs . If, on the other hand, Ra( Rs ) > Rs /τ , then the loop of oscillations appears above the steady branch (i.e. at lower | Ra|); the critical value of Rs for the onset of oscillations is then Rs so . It can be shown that Rs = Rs so when Figure 6(a) shows the critical value of Rs as a function of τ for the case of σ = 0.3. For this value of σ , the critical value of τ distinguishing the two regimes, given by expression (3.29), is τ = 0.649. We note also, from (3.28), that Rs so → −∞ as (σ − σ τ − τ 2 ) → 0; thus oscillations can be preferred in the third quadrant only if the following inequality is satisfied: (3.30) for σ = 0.3, this gives the lower bound of τ = 0.418. Figure 6(b) shows the corresponding critical value ofk 2 , as determined by expression (3.26), at the onset of oscillatory instability.

n > 2
In the third quadrant, the oscillatory mode remains minimised when k 2 = O(C T −1/2 ) for n > 2. This property thus allows us readily to determine the oscillatory stability boundary from (3.23) (valid for n = 2) simply by considering the case of | Rs |, | Ra| 1. In this limit, (3.23) factorises to give the two leading-order solutions Recall that at leading order, the steady branch is given by Ra = Rs /τ . Thus, since τ < 1, the root Ra = Rs always lies above the steady branch and hence is of no physical relevance. However, the root Ra = (σ + τ ) Rs /σ lies below the steady branch provided that inequality (3.30) holds; i.e. the same condition applies for n > 2 as for n = 2.
Since there is cancellation at leading order in expression (3.25b) for the coefficient A, in order to determine the preferred value of the wavenumber at onset it is necessary to include the first-order correction to Ra, which gives (3.32) 3. For small enough |Rs|, the preferred mode at onset is always steady, with k 2 = 1/2, as shown in figure 7(b). For n > 2, since inequality (3.30) is satisfied with this choice of parameters, oscillatory modes are preferred; figure 7(b) confirms that the critical wavenumber for the oscillatory modes has k 2 = O(C T −1/2 ) for all n > 2. If there is a transition between steady and oscillatory modes, as here, it occurs when n ≈ 2; for the parameter values in figure 7, the transition occurs at Rs = −7.9 × 10 6 . For both the steady and oscillatory branches, Ra c varies linearly with Rs; the offset between the two is made clear in the inset, which plots Ra c /|Rs| against Rs.

Stability boundaries
As for the case of C S = 0, the coefficient a 5 defined by (2.18a) vanishes; the dispersion relation (2.17) reduces to a quartic equation with coefficients Again, for an oscillatory root we have the conditions Eliminating Ra gives the following quadratic equation for ω 2 : Conversely, we could eliminate ω 2 to obtain the following quadratic expression for Ra on the oscillatory boundary (provided that ω 2 > 0): We first recall that once Rs is such that inequality (2.24) is satisfied then oscillatory modes are preferred. In an analogous fashion to the case of C S = 0, we may investigate the changes from the classical problem when n = 1. We consider the regime C S 1,  Figure 8. (a) The real parts of the roots of (4.9) for k 2 as a function of Rs for the n = 1 regime; the blue line denotes a real root, the red lines denote the real part of a conjugate pair; σ = 1, τ = 0.5. (b) Oscillatory stability boundary (blue solid line) for C S = 10 −3 , C T = 0, Rs = 2 × 10 4 (corresponding to Rs = C S Rs = 20), together with the zeroth-order asymptotic expression (4.7) (purple dashed line) and the first-order correction (4.8) (purple squares). (1), and write Ra = Ra 0 + C S Ra 1 . On substitution into (4.5), the leading-order terms give

Rs
The wavenumber dependence appears at the next order, where Ra 1 is given by (4.8) Stationary points (d Ra 1 /dk 2 = 0) are therefore given by the expression which is a quintic polynomial in k 2 . Figure 8(a) shows the evolution of the real parts of the roots of this polynomial as Rs is increased. By contrast with the case of C S = 0 (cf. figure 2a), no new positive stationary points arise. As Rs is increased, and contrary to the behaviour in the first quadrant when C S = 0 (see § 3.2.1), the minimum at k 2 = 1/2 of the classical problem first moves to slightly lower wavenumbers, for all σ and τ < 1, before then moving gradually to higher wavenumbers. The smallest value of k 2 c , taken over Rs , is given by the positive root of the quadratic equation (4.10) Figure 8(b) plots the oscillatory stability boundary, together with the zeroth-order approximation (4.7), and the first-order correction (4.8); the latter is almost indistinguishable from the boundary of the full system.

n = 2
In the first quadrant, the critical Ra for oscillatory modes is minimised when k 2 = O(C S −1/2 ). At leading order, the coefficients of (4.5) become In terms of the scaled variables defined by Ra = C S −2 Ra, Rs = C S −2 Rs , this gives the following quadratic equation: It may be verified that in this regime, oscillations are favoured for all values of the parameters. To see this, we note that the marginal curve for oscillations (4.12) is a hyperbola passing through the origin, which can be written in the form (4.13) This hyperbola intersects the line Rs = 0 at Ra = 0, Ra = 1/4σ . Thus it is the lower branch of the hyperbola that passes through the origin, and the gradient along this branch decreases monotonically along the branch from (σ + τ )/(1 + σ ) to σ/(1 + σ ); hence the gradient is always less than 1/τ , and oscillations are therefore always preferred in the first quadrant when Rs = O(C S −n ), n ≥ 2. As in § 3.3, we must take the analysis to the next order to determine the wavenumber dependence of the critical Rayleigh number; we thus write Ra = Ra 0 + C S 1/2 Ra 1 . The leading-order terms determine Ra 0 via (4.12); Ra 1 is determined at the next order through the expression F Ra 1 = Ak 2 + B k 2 , (4.14) where If F, A and B are of the same sign, then Ra 1 is minimised wheñ (4.16) Figure 9 shows Ra (o) as a function of k 2 , together with the asymptotic results (4.12) and (4.14), for C S = 10 −4 and C S = 10 −8 . ; the purple squares include the first-order correction given by (4.14). At the smaller value of C S , the numerical results for the full system and the first-order asymptotic results are essentially indistinguishable.

n > 2
As in the third quadrant with C S = 0 ( § 3.3.2), for n > 2 the oscillatory mode remains minimised when k 2 = O(C S −1/2 ). The ensuing analysis to determine the critical wavenumber is thus entirely analogous, with the oscillatory stability boundary determined from (4.12) (valid for n = 2) by considering the case of | Rs |, | Ra| 1. In this limit, (4.12) factorises to give the following two leading-order solutions: where the latter is the critical value for the onset of oscillatory instability (which is always the preferred mode). To determine the critical value of the wavenumber at the onset of oscillatory instability, it is again necessary to include the first-order correction to Ra, which gives (4.18) From (4.16), having made use of (4.15) and (4.18), we can determine the critical value of k 2 at the onset of instability ask The onset of instability in the first quadrant for C T = 0 is related to that in the third quadrant for C S = 0. The leading-order expressions (3.32) and (4.18) are connected through the transformation (2.19). The critical wavenumbers are identical in the two cases, although this is not readily shown via the transformation (2.19) -the symmetry is lost since both analyses involve perturbations in Ra but not in Rs. Figure 10 shows Ra c and k 2 c for 10 −5 ≤ Rs ≤ 10 15 with C T = 10 −3 , σ = 1, τ = 0.5. The steady mode, with k 2 c = 1/2, gives way to the oscillatory mode when Rs = O(1). When it first becomes dominant, the oscillatory mode also has k 2 c = 1/2. As Rs is increased, 927 A13-21 in the n = 1 regime, k 2 c decreases slightly, as explained in § 4.2.1, before increasing sharply in the n = 2 regime. For n ≥ 2, k 2 c = O(C S −1/2 ), with the precise asymptotic value (4.19) for n > 2. Also shown in figure 10(a) is the boundary for the onset of oscillatory motion for the classical problem. In terms of Ra c , it can be seen that the M-C effect comes into play for n 2, although, since Ra c depends linearly on Rs for both these boundaries, the difference is small on a logarithmic scale; it is made clearer in the inset, which brings out the difference in the multiplicative factors in (2.23) and (4.17). It is though worth reiterating that the difference in the spatial structure of the preferred modes between the case of C S = 0 and that of C S 1 is significant; as illustrated in figure 10(b), the former has k 2 c = 1/2, the latter k 2 c = O(C S −1/2 ).

4.3.
Third quadrant: Ra < 0, Rs < 0 4.3.1. n = 2 As for the third quadrant in the case of C S = 0, discussed in § 3.3, steady convection in the form of salt fingers is always preferred for n < 2. For n = 2, oscillatory solutions arise for a range of wavenumbers lying between two Takens-Bogdanov points, defined as where the steady and oscillatory branches coincide, with the oscillatory branch having two zero frequencies. At such points, the coefficients a 0 and a 1 of the dispersion relation (given by expressions (4.1d,e)) must both vanish. For large |Rs| (n ≥ 2), simultaneous solution of the equations a 0 = 0 and a 1 = 0 requires large k 2 . Thus, making only the approximation that β 2 ≈ k 2 , eliminating Ra gives the following condition for Takens-Bogdanov points:  Figure 11. Critical Ra for the steady mode (red line) and oscillatory mode (blue line) for σ = 1, τ = 0.2, C S = 10 −3 , C T = 0 with (a) Rs = −2.5 × 10 6 and (b) Rs = −6.5 × 10 6 . In (a), steady convection is favoured, with the preferred mode having k 2 = 0.5; in (b), oscillatory convection is favoured, with the preferred mode having k 2 = 5954. roots for k 2 in (4.20); this occurs when When Rs = Rs TB , Ra is given by Thus, oscillatory solutions first appear in the third quadrant if τ > σ/(1 + 2σ ); otherwise they first appear in the second quadrant. When oscillatory solutions first appear, at Rs = Rs TB , they are not favoured, even when (Rs TB , Ra TB ) lies in the third quadrant. As Rs is decreased below Rs TB , the extent of k 2 for which oscillatory solutions exist increases, but initially without oscillations being preferred, as shown in figure 11(a). However, as Rs is decreased further -but still within the n = 2 regime -the minimum of the oscillatory loop lies below the minimum of the steady branch and oscillations are then favoured; this is the situation shown in figure 11(b). Analytically, it is not possible to pin down the transition values of Rs and k 2 precisely since, save for the approximation β 2 ≈ k 2 , no additional simplification is possible of the terms in (4.5), which are all of the same order.

n > 2
For n > 2, expression (4.20) for the Takens-Bogdanov points, being based only on the condition k 2 1, is still valid, with roots (4.23) Thus, in comparison with the case of n = 2, the arc of oscillatory instability is extended; the lower value of k 2 remains O(C S −1 ), with the upper value now O(C S 1−n ). The minimum of the oscillatory stability boundary is captured in the same scaling for k 2 as the upper Takens-Bogdanov point. Thus, with k 2 ∼ C S 1−n , the frequency equation ( at leading order, which factorises as The relevant root for the stability boundary is given by ω 2 = τ k 2 /2C S , which is always admissible. The oscillatory stability boundary is then determined from the leading-order terms in expression (2.22) to be Since the onset of steady convection occurs when Ra = Rs/τ = O(C S −n ), then, for n > 2, oscillatory convection, with onset given by (4.27), is always preferred. Figure 12 shows Ra c and k 2 c for −10 15 ≤ Rs ≤ −4, for fixed values of C S = 10 −3 , C T = 0, τ = 0.5, σ = 0.3. For small enough |Rs|, the preferred mode at onset is always steady, whereas for n 2, oscillatory modes are preferred. There is an abrupt change in scale of the preferred mode, from k 2 c = 1/2 to k 2 c = O(C S −1 ) in the n = 2 regime. For n > 2, oscillatory modes are vastly preferred, with Ra c scaling as Rs 2 ; k 2 c increases linearly with Rs.

Conclusion and discussion
We have investigated the onset of linear instability for double-diffusive convection, incorporating, separately, the M-C effects for heat and salinity. In most physical situations, the M-C effects are very small, and we have therefore concentrated our attention on the case of C 1, where here in the discussion, we use C to denote either C T or C S . The linear theory for classical double-diffusive convection is fairly straightforward, chiefly because the most favourable horizontal wavenumber for instability is independent of the governing parameters of the problem -in contrast to the somewhat related problems of rotating convection or magnetoconvection. However, for two reasons, the introduction of either one of the M-C effects -although 'just' an extra time derivative -leads to a vastly richer problem. The first is that the additional time derivative is multiplied by the M-C coefficient and hence represents a singular perturbation in the time domain. The second is that the dependence on wavenumber becomes decidedly non-trivial.
Since C is small, the influence of the M-C effect is felt only for concomitantly high thermal and salinity gradients. This competition is most elegantly treated by determining the critical thermal Rayleigh number Ra c for the regimes defined by Rs = O(C −n ) (n > 0). The only symmetry of the system is given by (2.19); thus the two problems of C S = 0 and C T = 0 require separate analyses. The chief influence of both M-C effects is to favour oscillatory convection, with the preferred mode having a large wavenumber. Oscillatory instabilities in double-diffusive convection may be regarded as capitalising on a favourable phase lag in the two components that contribute to the density, allowing a fluid parcel displaced upwards to return to its original position with a higher density than it had initially; it will then overshoot, with repetition of the process leading to growing oscillations. The M-C effects introduce a new fast time scale -this can then be exploited by modes of small horizontal scale to lead to new oscillatory modes of instability.
For classical double-diffusive convection, the onset of instability with Rs and Ra both positive (the first quadrant) is steady for small Rs, but oscillatory once inequality (2.24) is satisfied (the diffusive regime). By contrast, the onset of instability with Rs and Ra both negative (the third quadrant) is steady for all values of Rs (the salt fingering regime).
In the first quadrant, for C S = 0 or C T = 0, there is some influence of the M-C effect when n ≈ 1, but this represents only a small change in the critical value of k 2 from its classical value of k 2 c = 1/2. Much more drastic changes arise when n ≈ 2. For C S = 0, with n ≥ 2, the preferred mode becomes small scale, with k 2 c = O(C T −n/2 ); the critical Rayleigh number is given by Ra c = O(C T −(1+n/2) ) = O(C T n/2−1 Rs). For n = 2, Ra c = O(Rs) and hence is of the same order of magnitude as the critical Ra of the classical problem -though the preferred mode is characterised by k 2 c = O(C T −1 ) rather than k 2 = O(1). For n > 2, Ra c is asymptotically smaller than Rs (and hence than the critical Ra of the classical problem). These results are captured in figure 4. For C T = 0, as shown in figure 10, the preferred mode has k 2 c = O(C S −1/2 ) for n ≥ 2, with Ra c = O(Rs); here Ra c is of the same order of magnitude, but less than the critical value for the classical problem.
In the third quadrant, in which there is no oscillatory instability in the classical problem, M-C effects are not felt until n = 2. For C S = 0, oscillatory instabilities are favoured only for sufficiently small σ , given by inequality  figure 12, for n > 2, the oscillatory branch then has a distinct asymptotic scaling to that of the steady branch.
The preference for small horizontal scales (k 1) of the preferred modes at large Rs has two interesting consequences. The first concerns the vertical structure of the perturbed quantities, as described in expression (2.15). In the leading-order analysis of § § 3 and 4, we make use of the approximation β 2 = k 2 + 1 ≈ k 2 . It is, though, of course the case that for modes with vertical wavenumber m > 1, the approximation β 2 = k 2 + m 2 ≈ k 2 still holds: leading-order expressions for k 2 1 are therefore degenerate in the vertical wavenumber. Unlike for classical double-diffusive convection, determining the preferred value of m is thus not immediately apparent, and requires examination of the higher-order terms. Here, to avoid further complications, we have considered only the m = 1 mode. However, numerically, we have also considered the case of m > 1 and found that in all of the examples considered, the m = 1 mode is preferred. The second consequence of the preference for k 1 modes concerns the influence of the boundary conditions on the planes z = 0, π. In this paper, for algebraic simplicity, we have assumed that the boundaries are stress-free. For no-slip boundaries, it can be shown that the dominant terms in the interior give the same solutions as in the stress-free case, with the new boundary conditions accommodated by viscous boundary layers. This is reminiscent of the behaviour for convection in an imposed strong magnetic field, discussed in more detail in Eltayeb et al. (2020).
We have concentrated in this paper on determining the conditions for the onset of instability and the nature of the preferred mode at onset. An alternative approach to the linear problem is to consider the mode of maximum growth rate in the unstable regime. Indeed, sufficiently far into the unstable regime, close to the line Ra = Rs, we expect steady convection to dominate and for the modes to resemble those of classical double-diffusive convection, perturbed slightly by the M-C effect. What our analysis has revealed, however, is that there are regions in (Rs, Ra) space that are stable in the classical problem, but that can be destabilised by a new instability mechanism arising from the M-C effect. This feature is particularly marked in the first quadrant of the (Rs, Ra) plane for C S = 0 and the third quadrant for C T = 0, for which the power-law dependence of Ra c on Rs is different for the classical and M-C problems.
It is important to see where M-C effects may play a role in double-diffusive convection in physical systems. As we have shown, even very small values of C can be significant, provided that the Rayleigh numbers are sufficiently high. In an astrophysical context, double diffusion in the diffusive regime (often referred to as semiconvection in the astrophysical literature) can be of importance for massive stars and stars on the horizontal branch. Here, one indeed expects very high values of Ra and Rs. Although estimating the Rayleigh numbers is not straightforward, one might assume that Ra and Rs are comparable to Ra in the Sun, which is O(10 20 ) at the base of the convection zone, with the length d taken to be a scale height H (Ossendrijver 2003). We may then ask what value of τ T (or τ S ) will give C = O(Rs −1 ), i.e. the n = 1 regime? With the solar values of H = O(10 8 m) and κ = O(10 3 m 2 s −1 ), we obtain τ T = O(10 −7 s). Direct calculations of the relaxation times for gases in stellar interiors are elusive, but for comparison, we note that kinetic calculations of τ T for a rarefied gas yield τ T = O(10 −9 s) (Carrassi & Morro 1972;Stranges et al. 2016). Thus, in 'normal' stellar interiors, the n = 1 regime may be accessible, but probably not the regimes for larger n, which would require unfeasibly large values of the relaxation time. However, as pointed out by Herrera & Falcón (1995), in dense degenerate stellar material, where thermal conductivity is dominated by electrons, one does indeed expect a much longer relaxation time, owing to the larger mean free path of electrons.
At the other extreme of very small length scales -as found in microfluidic devices, for example -values of C are readily obtained that are small, but not outrageously so. For example, for hydrogen at room temperature and atmospheric pressure, C T = O(10 −2 ) when d = 10 −6 m. Under such conditions, the regimes for n 2 are attainable at fairly modest values of Ra and Rs. Microfluidic platforms offer precise capabilities in controlling mass transport for biological and medical studies (Kuo & Chiu 2011). Some form of double-diffusive convection, influenced by M-C effects, may then play a role in manipulating chemical concentrations at the microscale. It does not appear that there has been any attempt either to calculate τ S theoretically or measure it experimentally; we hope that our study may trigger work in this direction.
The natural extension of the work described in this paper, which we are currently pursuing, is to consider the linear stability problem in which C T and C S are small but both non-zero; this clearly adds another level of complexity, particularly since the relative sizes of C T and C S will also be an important factor. In terms of investigating the nonlinear evolution, there are several features of interest that are contingent on the presence of the M-C effects. As shown in both § § 3 and 4, and unlike in classical double diffusion, steady and oscillatory modes of vastly different horizontal scales can be equally readily excited, thus prompting the question of the emergent dominant nonlinear scale. Nonlinear simulations of double-diffusive convection often employ periodic boundary conditions for the perturbations in the vertical as well as horizontal directions. Physically, this reduces the impact of impermeable boundaries and is thus more relevant in geophysical and astrophysical settings. Mathematically, it allows the presence of 'elevator' modes (i.e. modes with no vertical structure) and also couches the problem in terms of the density ratio R 0 = |Ra|/|Rs|, rather than Ra and Rs individually (see, for example, Hughes & Brummell 2021). This approach has proved particularly fruitful in the study of the formation and maintenance of layers or 'staircases', which can appear in both the diffusive and fingering regimes (e.g. Radko 2003;Stellmach et al. 2011;Mirouh et al. 2012). One of the most important consequence of the formation of double-diffusive staircases is the marked increase in turbulent transport, in comparison with the unlayered state. It would therefore certainly be of interest to pursue the nonlinear development of the instabilities studied here, but with such modified boundary conditions, in order to investigate how the layering process and associated turbulent transport may be influenced by the new time scale and very different dynamics introduced by the M-C effects.