The effect of temperature-dependent viscosity and thermal conductivity on the onset of compressible convection

The linear equations of thermal convection in a compressible fluid with non-constant transport coefficients are derived. The criterion for the onset of convection is established, based on linear stability analysis, for a range of different temperature-dependent profiles of thermal conductivity and viscosity. Temperature-dependent transport coefficients are shown to lead to a more complex behaviour than their constant counterparts, and modifies the stability condition of the fluid. When the Rayleigh number is defined in terms of the mid-layer physical properties and the temperature gradient at the top is held constant, increasing the temperature-dependence of thermal conductivity is found to raise the critical Rayleigh number dramatically, as the convective disturbance is then concentrated mainly at the top of the layer. In contrast, for viscosity a more subtle effect on stability is identified.


Introduction
High-resolution observations of the solar surface are continually revealing a variety of multi-scale magnetic features (e.g. Harrison 2008, Wiegelmann et al. 2014. In order to explain these observed structures and activities, research has been conducted over many years, and it is known that the dynamics in the convection zone significantly contribute to what is observed on, and above, the surface of the Sun (e.g. Galloway and Weiss 1981, Cattaneo et al. 2003, Weiss and Proctor 2014. In order to understand convection in stellar interiors, such as in the Sun, extensive studies have been carried out to explore convective instabilities (see, e.g. Goody 1956, Chandrasekhar 1961, Jones and Moore 1978, Zappoli et al. 2014. However, the formalism of most models introduced simplifying mathematical assumptions to reduce the complexity of the problem considerably and employ analytical methods, in order to provide a preliminary exploration of the physical processes. Many have utilised the Boussinesq approximation, in which density variations are neglected in the governing equations, except when the gravitational force is considered (see particularly Boussinesq 1903, Spiegel andVeronis 1960).
In some astrophysical and geophysical systems, particularly in stars like the Sun, the inevitable large length-scales result in substantial density variations through the deep convective layer. Therefore, the Boussinesq approach is an oversimplication to explore the many scale heights of solar convection. Spiegel (1965) constructed the linearised system for a fully compressible medium, to determine the onset of steady convection, and provide an insight into the non-linear evolution of compressible convection. However despite incorporating compressibility, simplifying assumptions regarding the transport coefficients still remained. In this paper, as in others (see, e.g. Gough et al. 1976, Calkins et al. 2014, Liu and Sun 2019, the thermal conductivity and viscosity are taken spatially uniform, and independent of the thermodynamic variables. In some situations, it may be acceptable to consider such approximations, e.g. if the vertical extent of temperature is sufficiently small. But generally, the transport coefficients are not constant and depend on both the magnetic field and temperature (Priest 1982, Spitzer 2006. For a fully ionised plasma, both viscosity and thermal conductivity are dominantly proportional to T 5/2 . Thus, dependency will be important for calculations where there is a large temperature difference between the top and bottom of the domain as in the solar convection zone, and there is the potential for complex local dynamics. Few attempts have been undertaken in order to understand the effect of non-constant transport coefficients in linear studies of convection, in a compressible atmosphere. Vickers (1971) considered a position-dependent viscosity and a conductivity that was a function of temperature, which was later extended by Graham and Moore (1978). Appropriate nonconstant transport coefficients have been shown in some contexts to significantly alter the stability threshold of a convecting fluid (Glatzmaier andGilman 1981, Drew et al. 1995). Here we will explore the effect of non-constant transport coefficients for a fully compressible fluid by deriving the hydrodynamic equations for marginal stability to incorporate the temperature-dependent thermal conductivity and viscosity in an arbitrary form that will allow employing the temperature relations outlined in Spitzer (2006). Linear stability analysis will be conducted to investigate the stability threshold of a fully compressible stratified system. This paper will proceed as follows: Sections 2 and 3 outline the model, numerical approach, and parameter selection. Section 4 discusses the results and this is followed by the conclusions in section 5.

Derivation of linear stability equations
A three-dimensional domain in a Cartesian geometry is considered, with the x and y coordinates representing the horizontal directions, and the z-axis pointing vertically downwards, parallel to the constant gravitational force, g. In dimensional form, the set of equations for the evolution of a compressible hydrodynamic system are We also have where R * is the gas constant, and the viscous stress tensor In the above equations, ρ is the fluid density, p is the pressure, and T is the temperature. The transport coefficients, K and μ are assumed functions of temperature, and take the following form where T * , K 0 , and μ 0 are the reference temperature, thermal conductivity, and viscosity respectively, taken to be the value at the top of the domain. The indices q and r for a fully ionised hydrogen plasma are q = 2.5 and r = 2.5 (Priest 1982, Spitzer 2006. Here, the indices are assumed constants that will be varied. Let ρ = ρ 0 + ρ , T = T 0 + T , u = u 0 + u , and p = p 0 + p , where ρ 0 , T 0 , u 0 , and p 0 are the static solutions of the governing equations and ρ , T , u , and p are to have infinitesimal amplitudes. Assuming a basic state where u 0 = 0, T = T 0 (z), and ρ = ρ 0 (z), we obtain dp 0 dz The evolution equations for the perturbed quantities are derived by inserting the decomposed variables and linearising. Hence, we obtain where τ is the stress tensor for the perturbed component of the velocity.
To obtain non-dimensional form of the equations, the unit of length is scaled by the depth of the layer, d, density and temperature are scaled by their initial values at the upper surface, ρ * and T * respectively. Velocity is scaled by the sound travel time across the layer in terms of the isothermal sound speed, √ R * T * , and is related to the unit of time d/ √ R * T * . Thus, Each perturbed quantity can be expressed in the form f (z) exp(ikx + ily + st), where k, l ∈ R are the horizontal wavenumbers, s = s r + is i ∈ C with s r being the growth rate of the instability, and f describes the variation of the disturbances across the layer. To simplify the problem further, Squire's transformation can be applied to reduce the three-dimensional problem to an equivalent two-dimensional problem (see, e.g. Squire 1933, Drazin andReid 2004). Thus, (6a-c) are reduced to an equivalent two-dimensional problem in which D ≡ d/dz. The system (7a-d) can be written in the form sf = Af, where the vector f = [ρ , u , w , T ] T is the solution vector containing the eigenfunctions, and A is a 4 × 4 matrix that consists of linear differential operators in z. Equations (7a-d) are solved numerically by dividing the layer depth 0 ≤ z ≤ d into n uniformly distributed points, and the differential operators are approximated using a central fourth-order finite difference approximation. The method we use is adapted from the general method discussed in Favier et al. (2012) and Witzke et al. (2015), where the Schur factorisation is utilised to determine the eigenvalues and eigenvectors.

Boundary conditions, initial conditions, and parameter choices
We will employ impermeable, stress-free velocity and constant temperature boundary conditions on the top and bottom boundaries where θ = T/T * is the temperature stratification. The equilibrium temperature and density distributions are given by respectively. The background density and temperature profiles vary depending on the values of q and is unaffected by r since viscosity, in this case, does not modify the basic state. Figure 1 illustrates the differences in the initial temperature and density profiles for two q values. For this preliminary investigation, our principal objective is to understand the general influence of non-constant transport coefficients on stability. We choose to present our survey of critical values of instability by setting the thermal diffusivity, C k , as the control parameter and fixing the thermal stratification, θ, and the Prandtl number, σ , to unity in all cases while the variable parameters are the polytropic index, m, thermal conductivity index, q, and viscosity index, r. A summary of the input parameters is shown in table 1.
The resolution for this problem has been carefully selected and tested for the results quoted in this work. Due to the large density contrast, as q increases, we are limited by our choices of q since it requires more n-points than computational resources allow. For example, for q = 2.5 and m = 1.0, the system remains under-resolved for n = 3000.

Results
Here we focus attention on the stability threshold while we vary the characteristics of the transport coefficients. Equations (7a-d) are solved numerically for a range of wavenumbers, with stability threshold solutions discussed in terms of the critical Rayleigh number, in addition to the thermal diffusivity, for varying polytropic index, thermal conductivity index, and viscosity index. The definition of Vickers (1971) for the Rayleigh number is where the subscript 1 refers to the value of the variable evaluated at a level within the layer. Note that the form (10) of the Rayleigh number differs slightly from Spiegel (1965), where constant transport coefficients are considered. However, to enable direct comparisons between constant and non-constant transport coefficients, through varying q and r, the definition (10) of the Rayleigh number is modified to as in Spiegel's definition of the Rayleigh number, such that the superadiabatic temperature gradient is included. The Rayleigh number is depth-dependent and so, we follow the common convention of expressing Ra in terms of the mid-layer value for all cases (see, e.g. Brandenburg et al. 1990, Hurlburt et al. 1994. To determine the critical Rayleigh number, Ra c , the eigenvalue problem is solved, such that the most unstable mode is found for each parameter regime. By setting the indices of  the transport coefficients to zero, we reduce the problem to the constant thermal conductivity and viscosity case. Therefore, part of our results can be compared to Gough et al. (1976). For each fixed m, the critical Rayleigh number Ra c and associated critical wavenumber k c are found and presented in table 2 for varying q and fixed r = 0, and in table 3 for varying r and fixed q = 0, such that for Ra > Ra c convection will ensue. As values of Ra become larger than Ra c , the band of wavenumbers k, in which perturbation grows, expand. Note, the system initiates close to the Boussinesq limit for stability of Ra c = 27π 4 /4, since the density variation is minimal for small m, q, and r.
Here we are focused on discovering how Ra c behaves with varying q and r. The quantitative indication of our results reveals the stabilisation of the compressible medium as both q reaches the Spitzer q value, 2.5, in table 2, and r reaches the Spitzer r value, 2.5, in table 3. However, stabilisation is more effective as q is increased. This remark mainly lies in the effect of the large density stratifications created by increasing q. In other words, the substantial increase in Ra is primarily due to changes in the depth-varying background temperature profile as the exponent q is varied. Increasing q greatly reduces the superadiabatic gradient near the bottom boundary which consequently enhances density stratification. In agreement with the results of Gough et al. (1976) and Calkins et al. (2014), increasing stratification is found to lead to a growth in the critical Rayleigh number.
Although we chose to fix the thermal stratification to unity as we first and foremost want to focus on q and r variation, we have conducted some additional cases to examine how varying the thermal stratification influences the onset of convection. Given the particular choice of m = 0.5, q = 0.1, and r = 0, the critical Rayleigh numbers (wavenumbers) are 708.47 (2.2), 896.33 (2.2), 1267.54 (2.39), and 1826.57 (2.39) for θ = 0.5, 2, 5, and 10, respectively and so we find that increasing the thermal stratification is found to enhance the stability of the system.
The plots in figure 2 illustrates how variations in q (for fixed r) and r (for fixed q) influence the marginal stability for convection in terms of the thermal diffusivity, C k , for a range of polytropic indices, m. By employing a logarithmic scale, we notice a decrease in C k with increasing q. The onset of instability becomes significantly smaller as q grows, which leads to an increase in the critical Rayleigh values (as shown in table 2). It was not possible to determine the critical Ra values for large q at large m due to computational limitations.
When varying viscosity through the parameter r, an almost linear decrease in the value of C k for onset of instability was found in figure 2 for different values of the polytropic index. Generally, as the polytropic index reaches the adiabatic limit, the system becomes more stable, and so the thermal diffusivity required for onset of the instability is expected to decrease. Interestingly, this is not true until m reaches 0.7. This behaviour in the polytropic index was also observed for varying thermal conductivity when q 0.5. The thermal diffusivity transition, with respect to the polytropic index, cannot be extracted directly from tables 2 and 3, given our definition of the Rayleigh number. This is because both the reduction of C k and the increase of the polytropic index promotes system stability.
Profiles of the marginal eigenfunctions for the vertical velocity, temperature, and density are shown in figure 3, to provide a better insight of the local dynamical properties. The eigenfunctions displayed in this figure are for two values of the polytropic index, m = 0.1 and m = 1, for different values of q. In the constant transport coefficient regime, where q = 0, the shape of the temperature eigenfunctions can be seen to be parabolic in nature with a symmetry around z = 0.5 for both m = 0.1 and m = 1. However, by increasing q, changes in the general form of the eigenfunctions are revealed. The symmetry of the eigenfunctions breaks and becomes skewed towards the upper boundary. This is because the temperature gradient, as q increases, becomes steeper near the top of the layer which therefore is prone convective instability and produces this deviation in the eigenfunctions. Convection is driven in superadiabatic regions though it can penetrate to some extent into the subadiabatic regions below. Using our expression for the equilibrium temperature (equation (9a)) we can write the superadiabatic temperature gradient as From this expression we see that at m = 1 and q = 2 (recall θ = 1) the superadiabatic region is confined to a thin layer near z = 0, but that for lower m and lower q it reaches deeper into the layer. Generally, growth in temperature produces a decline in density, which consequently causes fluid motion due to pressure and other forces when differences in density occur under the influence of gravity (Böhm-Vitense 1992, Subramanian andBalasubramaniam 2001). This pattern can be seen in both eigenfunctions for density and velocity, where the peak height increases with increasing q. By comparing the eigenfunctions for m = 0.1 and m = 1, we examine that the skewness to the top boundary has sharper peaks for large m. Further, the eigenfunctions are zero at the bottom, and the onset of convection takes place at a sublayer near the upper domain. This is again due to the nature of the background density, where the system is less dense at the upper boundary (as imposed by the boundary conditions). As the fluid becomes more dense when moving downwards, the effect of buoyancy braking overtakes as the dominant mechanism driving the flow, where buoyancy forces act as a brake (see, e.g. Hurlburt et al. 1986), and so leading to enhanced stability. Massaguer and Zahn (1980) and Hurlburt (1983) have reported that buoyancy braking is responsible for the enhancing effect of stabilisation in layers with large density. The perturbed temperature, density, and vertical velocity eigenfunctions for constant and temperature-dependent viscosities, for various values of r, are depicted in figure 4. As r is varied, the symmetry is largely maintained with increasing r and is only mildly skewed. This is predominantly because the basic state remains almost unchanged for varying r, as opposed to varying q. The eigenfunctions at onset indicate that the temperaturedependancy of viscosity through varying r enhances the temperature perturbation about the midpoint of the layer, which accordingly attenuates density perturbations. However, as the viscosity increases with temperature, the effect on convective instability becomes evident. The vertical velocity eigenfunction reveals that increasing r is to inhibit the fluid flow, which therefore contributes to delaying the onset of convection. For larger m, a qualitatively comparable behaviour is observed, with the addition of the enhanced impact of stabilisation due to density stratification.

Conclusions
The form of transport coefficients, dynamic viscosity and thermal conductivity, are often simplified to be constant in stellar modelling. To explore the impact of non-constant transport coefficients in the modelling of convective instabilities, a general form of the equations governing thermal convection in a compressible polytropic atmosphere, using the Spitzer relations for temperature-varying thermal conductivity and viscosity, were derived for the first time and the stability of the system was examined using linear stability analysis for each non-constant transport coefficient separately. The linear equations were solved numerically to determine the nature of the unstable modes, together with the structure of the eigenfunctions.
Our investigation provides insight into the complexity introduced by non-constant transport coefficients. Here, we focused on the conditions for instability. In all cases, the calculated values of the critical Rayleigh number for marginal stability was found to be consistently higher as the indices q and r were increased. For large q, the onset of convection was found to require huge computational efforts, which restricted the parameter space. Nevertheless, given that the structure of the static atmosphere was to a large extent determined by q, increasing q was found to lead to large density stratifications, which led to stabilising the onset of convection significantly.
For fluids with high viscosity, resulting from increasing the parameter r, the stabilising effect was found to be weaker, in comparison to varying q. This is due to the small changes in the distribution of the basic temperature and density profiles, as r was varied. The polytropic index in all of the cases investigated in this study was shown to stabilise the system as it reached the adiabatic limit. An interesting result is the behaviour of the thermal diffusivity for m < 0.7, where the thermal diffusivity was observed to increase as the polytropic index increases.
The general behaviour of the eigenfunctions for the temperature-dependent transport coefficients revealed a far more interesting impact on the convecting fluid. This reinforces the notion that the choice of transport coefficients can have a great influence on the overall dynamics of convection, and thus appropriate form of coefficients should be carefully selected when modelling stellar interiors.