format: text

text:
# The electroclinic effect and modulated phases in smectic liquid crystals

## I Introduction

Chiral liquid crystals exhibit a variety of interesting phenomena Kamien . Chiral molecules tend to pack at a slight angle with respect to their neighbors. This preferred packing leads to spontaneous twist in the molecular orientation, producing phases such as the cholesteric and helically twisted smectic–$`C^{}`$ (Sm–$`C^{}`$). In the case of two–dimensional Sm–$`C^{}`$ films, the chiral modulation of the molecular orientation cannot occur in a defect–free fashion. Instead, for large chirality, the system will form a modulated phase consisting of a regular network of defect walls and possibly points. One possible modulated phase is a “striped” phase consisting of a regular array of parallel domain walls (Fig. 1). In between these walls the c director rotates in the sense preferred by the chirality. Another possible modulated phase consists of a hexagonal array of domain walls, with point defects located at the centers and corners of the hexagonal cells (Fig. 2).

Chirality also leads to unusual effects in orientationally disordered phases. As shown by Garoff and Meyer Garoff  using symmetry arguments, if an electric field is applied to a system of chiral molecules in the Sm–A phase, an average molecular tilt will be produced, a response termed the electroclinic effect. These symmetry arguments do not, however, indicate the microscopic basis of the electroclinic effect. To date, two microscopic models have been proposed. In the first, the long axes of the molecules are assumed to be perpendicular to the smectic layers in the Sm–A phase. The electric field causes the molecules to tilt, with an accompanying significant change in the smectic layer spacing. In the second scenario suggested by de Vries deVries , the molecules in the Sm–A phase are assumed to be tilted even in the absence of the electric field, but azimuthally disordered. The electric field then induces azimuthal order, with no significant change in the smectic layer spacing. Experimental studies have shown evidence for both models, at least in terms of the field–induced layer contraction. Recently Selinger et al. Selinger  noted that the optical birefringence of the liquid crystal will also distinguish between the two microscopic pictures. Because the field induces orientational order in the second model while simply tilting the already aligned molecules in the first, the birefringence should change dramatically in the former model, but not the latter. It appears that each of these models has experimental validity, namely there are materials with substantial layer contraction and weak field dependence in the birefringence, and other materials with the opposite properties. However, even a detailed theory based on the second model Selinger  is unable to fully describe the experimental data. Furthermore, many materials show an unusually large electroclinic effect, well into the Sm–A phase. While it is tempting to associate this large effect with strong critical fluctuations Selinger , a completely different molecular mechanism may be relevant.

In this paper we propose a new model for the microscopic origin of the electroclinic effect, one that shares some of the macroscopic characteristics of the de Vries scenario, specifically, no layer contraction and field dependent birefringence. We postulate that what is assumed to be a Sm–A phase in some materials close to the Sm–$`C^{}A`$ transition, could in fact be a modulated version of the smectic–$`C^{}`$ phase, a network of defects, walls and lines, but one where the spacing between the defects is in the subvisible range. This modulated phase would have no long–range azimuthal order or net polarization in the absence of an electric field. Thus, “to the eye”, this phase would appear to be a Sm–A phase. However, the application of an electric field distorts the structure of the defect array and leads to a large electroclinic effect. We now describe our proposed modulated phase and its response to an electric field in detail.

A bulk chiral Sm–$`C^{}`$ can be described by the following elastic free energy deGennes :

$`F`$ $`=`$ $`{\displaystyle }d^3x[{\displaystyle \frac{1}{2}}K_1(\mathbf{}\mathbf{}𝒄)^2+{\displaystyle \frac{1}{2}}K_2(𝒄\mathbf{}\mathbf{}\mathbf{\times }𝒄+q)^2+{\displaystyle \frac{1}{2}}K_{13}(\widehat{𝒛}\mathbf{}\mathbf{}\mathbf{\times }𝒄)(𝒄\mathbf{}\mathbf{}\mathbf{\times }𝒄)`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}K_3(\widehat{𝒛}\mathbf{}\mathbf{}\mathbf{\times }𝒄)^2w\widehat{𝒛}\mathbf{}\mathbf{}\mathbf{\times }𝒄]`$
where $`𝒄`$ is the projection of the nematic director on the plane of the layers. We assume that the molecular tilt is uniform and thus $`𝒄`$ is assumed to be of fixed length. Chirality appears in two terms in this free energy. First, chirality induces a spontaneous twist of the $`𝒄`$ director as indicated in the second term on the right hand side of Eq. (1), which is minimized in a helical state with pitch $`2\pi /q`$. Second, chirality induces spontaneous bend within each smectic layer as indicated by the last term in the free energy, where we assume $`w>0`$. The free energy as written is zero for the perfect helical state with pitch $`2\pi /q`$. Because the smectic–$`C^{}`$ phase is ferroelectric, with the spontaneous polarization perpendicular to $`𝒄`$, we must add to this free energy an electrostatic term arising from space charge density proportional to $`\widehat{𝒛}\mathbf{}\mathbf{}\mathbf{\times }𝒄`$. We ignore this term for now, since we will find that in our proposed modulated phase, it adds only a constant to the free energy density, independent of the detailed structure of our model.

Considerable work has been done on analyzing Eq. (1) in two–dimensional chiral smectic liquid crystals, where the terms proportional to $`K_2`$ and $`K_{13}`$ in the free energy above are absent, and thus, there will be no competition between the two chiral terms. While only one chiral term is present, it is impossible to impose a spontaneous bend throughout the system. Defects must be introduced which separate domains of spontaneous bend, thus leading to a modulated phase. Two types of modulated phases have been proposed, a “striped” phase (Fig. 1) with a one–dimensional modulation Langer ; Petschek1 ; Petschek2 , and a hexagonal phase (Fig. 2) with a two–dimensional modulation. In the former case, parallel domain walls separate regions of spontaneous bend, while the latter consists of a two–dimensional array of domain walls and point disclinations. A $`+1`$ point disclination lies at the center of each hexagonal cell, while the boundaries of the cell are domain walls where the $`𝒄`$ director rotates by $`\pi `$ (a description of internal structures of these walls in terms of a nonfixed length $`𝒄`$ director field is given in Petschek1 ; Petschek2 ). A $`\frac{1}{2}`$ point disclination is located at each corner of the cell. The hexagonal phase which involves nearly pure bend is favored over the striped phase when $`K_1K_3`$. The size of the unit cell is inversely proportional to the strength of the chirality parameter $`w`$. The free energy of the hexagonal lattice includes (see Eq. (5) below) an attractive logarithmic interaction between the $`+1`$ disclination at the center of the cell and the $`\frac{1}{2}`$ disclinations at the corners. Thus, strongly chiral materials, where the unit cell is small, will favor this interaction, and hence the hexagonal phase Kamien .

Returning to the full three–dimensional free energy, Eq. (1), we consider the possibility of a modulated phase. A bulk modulated phase must somehow compromise between the two competing chirality–induced tendencies: spontaneous in–layer bend and spontaneous twist about the layer normal. To describe the electroclinic effect in our proposed scenario the bulk modulated phase must have zero net polarization. Thus, a bulk modulated phase produced by simply extending the two–dimensional striped phase in the third direction is not acceptable, even if the pitch of the spontaneous twist were assumed to be infinitely large. A three–dimensional model for a thick chiral smectic film was proposed in Ref. Gorecka , where the in–plane modulation and domain walls are confined to surface regions, and the director twists in the interior of the film. Here we consider a space–filling modulated structure, with no need for surface structures.

Our starting point in constructing a suitable bulk modulated phase is the director configuration of the $`+1`$ point disclination at the center of the two–dimensional hexagonal cell (Fig. 2). Expressing the $`𝒄`$ director as $`𝒄=(\mathrm{cos}\theta (x,y),\mathrm{sin}\theta (x,y))`$, the angle $`\theta `$ is given by:

$$\theta (x,y)=\varphi (x,y)+\varphi _0,$$
(2)
where the polar angle $`\varphi (x,y)`$ is given by $`\varphi (x,y)=\mathrm{arctan}(y/x)`$, and $`\varphi _0=\pi /2`$ (for the sign of the chirality appropriate to the structure shown in Fig. 2; for the opposite sign of chirality, $`\varphi _0=\pi /2`$, and the director circulation will be clockwise). This configuration represents pure bend and describes the director configuration away from the domain walls which bound the unit hexagonal cell. We now form a three–dimensional simple hexagonal lattice, with a unit cell of height $`h`$, and cross–sectional area as in the two–dimensional case. We specify the director configuration as follows. In the midplane of the unit cell (taken to be the $`z=0`$ plane) we assume the perfect bend director configuration of Eq. (2). We then introduce twist in the $`z`$ direction, by writing the angle $`\theta (x,y,z)`$ of the $`𝒄`$ director as follows:

$$\theta (x,y,z)=\varphi (x,y)+\varphi _0(z),$$
(3)
with

$$\varphi _0(z)=\pi /2+\pi z/\mathrm{},h/2<z<h/2.$$
(4)
The pitch $`\mathrm{}`$ of this helical structure will be determined by minimizing the free energy density of the unit cell. With this model for the modulated phase, the $`K_{13}`$ and $`w`$ terms in the free energy, Eq. (1), can be combined into a term like the $`w`$ term with $`w_{eff}=w+\pi /\mathrm{}`$. We find that $`\mathrm{}`$ is so large that for simplicity we set $`w_{eff}=w`$.

Thus the unit cell has a $`+1`$ disclination line running along the central $`z`$ axis of the cell and parallel $`1/2`$ disclination lines along the edges of the cell. Within the cell the director twists along the $`z`$ axis, in the sense preferred by the chirality (see Fig. 3). The top and bottom faces of the cell (perpendicular to the $`z`$ axis) are domain walls where the director undergoes a rotation by $`2\pi h/\mathrm{}`$. The six other faces of the cell are walls where the $`𝒄`$ director rotates by $`\pi `$.

The free energy density of the hexagonal unit cell (of side $`R`$) relative to the helical state is given by:

$$\text{f}=\frac{2\pi K\mathrm{ln}(R/R_o)+K^{}}{3\sqrt{3}R^2/2}+K_2\left(q\frac{2\pi }{\mathrm{}}\right)^2\frac{2\mathrm{ln}3w\mathrm{}}{\pi hR}\mathrm{sin}(\pi h/2\mathrm{})+\frac{2A_T}{h}+\frac{2A_S}{\sqrt{3}R},$$
(5)
where we have assumed that $`K_1=K_3=K`$ to simplify the elastic energy expression associated with the disclination lines. The core size of the disclination lines is denoted by $`R_o`$, and we represent the domain wall energy per unit area of the top and bottom faces of the cell by $`A_T`$, and the wall energy per unit area of each side face by $`A_S`$. The term proportional to $`\mathrm{ln}(R/a)`$ includes the elastic energy per unit cell associated with the $`+1`$ disclination line at the center of the cell and the surrounding $`1/2`$ disclination Petschek2 . The term $`K^{}`$ represents the remainder of the splay–bend elastic energy per unit cell. The term proportional to $`w`$ is obtained by integrating the corresponding term in Eq. (1) over the hexagonal unit cell using the director profile Eq. (3), thus ignoring the distortion of this profile near the edges of the cell.

Before minimizing this free energy density we consider the response of the modulated state to an electric field, i.e., the electroclinic effect. We assume that the direction of the molecular dipole moment $`𝒑`$ is given by $`𝒄\mathbf{\times }\widehat{𝒛}`$. In the two–dimensional hexagonal modulated phase (see Fig. 2), the molecular dipole moments near the center of the cell will point radially outward, given the bend configuration of the $`𝒄`$ director. Thus, an electric field applied in the $`+x`$ direction will shift each $`+1`$ disclination in the opposite direction, inducing a net polarization in the $`+x`$ direction, as shown in Fig. 4.

For the three–dimensional cell we assume that the $`+1`$ disclination lines at the cell centers also shift rigidly in the opposite direction to the field, ignoring any possible bending of the line due to the twist of the director in each cell (this should be a higher–order effect in the electric susceptibility). Assuming that the disclination lines shift by a distance $`d`$ in the $`x`$ direction, we find that the polarization (the net dipole moment of a cell divided by its volume $`V`$) is given to leading order in $`d`$ by:

$$P=\frac{2dp\mathrm{}R}{Va^3}\mathrm{sin}(\pi h/2\mathrm{}),$$
(6)
where $`p`$ is the molecular dipole moment and $`a`$ is the average intermolecular spacing. In deriving $`P`$ we have assumed that the unit cell is a cylinder of radius $`R`$ centered on the original position of the line; ignoring the hexagonal shape of the cell merely changes the numerical prefactor appearing in the above expression.

The distance $`d`$ is determined by minimizing the change in the free energy density due to the presence of the field. Displacing the disclination line by a distance $`d`$ costs elastic energy proportional to $`Kh(d/R)^2`$. Thus, $`d`$ is given by:

$$d=\frac{p\mathrm{}ER^3}{Kha^3}\mathrm{sin}(\pi h/2\mathrm{}),$$
(7)
the polarization by:

$$P\frac{p^2\mathrm{}^2R^2E}{Ka^6h^2}\mathrm{sin}^2(\pi h/2\mathrm{}),$$
(8)
and the susceptibility by:

$$\chi _P\frac{dP}{dE}=\frac{p^2\mathrm{}^2R^2}{Ka^6h^2}\mathrm{sin}^2(\pi h/2\mathrm{}).$$
(9)
We now consider the minimization of $`f`$, which is a function of $`R,\mathrm{},`$ and $`h`$. Minimizing first with respect to $`R`$, we find the following condition:

$$R=\frac{4\pi K}{3f_W}\mathrm{ln}(R/R_o^{}),$$
(10)
where

$$f_W=\frac{w\mathrm{}\sqrt{3}\mathrm{ln}3}{\pi h}\mathrm{sin}(\pi h/2\mathrm{})A_S,$$
(11)
and

$$R_o^{}=R_o\mathrm{exp}\left(\frac{1}{2}\frac{K^{}}{2\pi K}\right).$$
(12)
Eq.(10) has a solution for $`R>R_o^{}`$ providing $`f_W>0`$. This condition is physically equivalent Petschek1  to requiring that the domain wall energy of the side walls of the cell be negative in a non–fixed length director model, where the wall energy includes surface contributions, i.e., contributions from the total gradient term $`\mathbf{}\mathbf{\times }𝒄`$. Numerical minimization of the free energy density, Eq. (5), with respect to $`R,\mathrm{},`$ and $`h`$, (choosing $`K=K_2=1,K^{}=A_T=A_S=0.1`$) indicates that for values of $`w`$ that yield $`f<0`$, (i.e., when the modulated state is favored over the helical state), $`h10^5,\mathrm{}/h100`$, and:

$$\frac{2\mathrm{}}{\pi h}\mathrm{sin}(\pi h/2\mathrm{})1.$$
(13)
The condition $`\mathrm{}h`$ is sensible, because in that case the director pattern in each cell will be nearly pure bend. In this limit the condition $`f_W>0`$ implies that $`A_S<w\sqrt{3}\mathrm{ln}3/2w`$.

With $`\mathrm{}`$ very large, the polarization space charge density, proportional to $`\widehat{𝒛}\mathbf{}\mathbf{}\mathbf{\times }𝒄`$ varies in each cell as $`1/r`$, with $`r`$ the distance from the central disclination, and is almost independent of $`z`$. This produces an electric field $`𝑬`$ which is almost purely radial, and constant in magnitude, considering the cell to be approximated by a cylinder. At the outer boundary of the cell there is a surface charge which precisely cancels the charge in the interior of the cell, so that there is no electric field outside each cell due to the charge distribution within the cell and on its boundaries. Thus the electrostatic energy density can be evaluated in a single cell, and is proportional to $`E^2`$, independent of the radius or height of the cell. Although there are also ionic charges present in liquid crystals, typical estimates of their screening length indicate that they would not alter the electric fields at the length scales we are considering.

Assuming that $`\mathrm{}`$ is large, as indicated by the numerical minimization of Eq. (5), (specifically, $`\mathrm{}q^1`$), and using Eq. (13), the free energy density Eq. (5) can be approximated as:

$$fK_2q^2+\frac{P_0^2}{8\pi ϵ}\frac{\sqrt{3}f_W^2(2\mathrm{ln}(R/R_o^{})1)}{8\pi K\mathrm{ln}^2(R/R_o^{})}.$$
(14)
in which the electrostatic term, with spontaneous polarization density $`P_0`$ and dielectric constant $`ϵ`$, has been included.

We determine $`R`$ by numerically solving Eq. (10), using Eq. (13), and choosing $`K=1,K^{}=A_S=0.1`$ and $`K_2q^2+P_0^2/(8\pi ϵ)=0.0025`$ ; the results are shown in Fig. 5.

Using these results for $`R`$ we obtain the approximate free energy density, Eq. (14), as a function of $`w`$ (as shown in Fig. 6, choosing $`K_2=1,A_T=0.1`$). We see that if $`w0.39`$, the hexagonal phase is favored over the helical state. The numerical minimization of the full expression for the free energy, Eq. (5), indicates that in this phase $`h10^5`$, in units of an intermolecular spacing. Thus, the hexagonal phase in a real system would likely consist of a single lattice plane of hexagonal cells, i.e., the system is effectively two–dimensional.

From Eqs. (9) and (13) we see that:

$$\chi _PR^2.$$
(15)
We plot $`\chi _P`$ as a function of $`w`$ in Fig. 7.

## II Conclusions

Returning now to the observations that motivated this work, we can see that if $`R`$ is a sub-visible length, then the modulated phase would resemble a Sm–A phase, optically. In our model, the apparent Sm–A to Sm–$`C^{}`$ transition, as temperature is lowered, would actually be mimicked by the gradual change of physical parameters, so that below the transition temperature, the helical Sm–$`C^{}`$ phase would have lower free energy than the modulated phase.

In our model, this would be a first order transition. Since the driving term in the free energy density, Eq. (5), favoring the modulated phase is the term proportional to $`w`$, one could say that if $`w`$ decreases as temperature is lowered while all other parameters remain constant, then when $`w`$ passes through a particular value (about 0.39 in our calculated example) a first order transition to the helical Sm–$`C^{}`$ phase would occur. Approaching this transition, $`R`$ would be increasing, as would the susceptibility, $`\chi _P`$, as seen in the figures above.

The value of $`R`$ could well remain sub-visible throughout the range of the modulated phase. The electroclinic effect could be large throughout the temperature range of the modulated phase, with a rapid increase as the transition to the helical Sm–$`C^{}`$ phase is approached. The mean smectic layer spacing would be almost constant through the phase transition. In the modulated phase it would be governed by an average between the layer spacing in the bulk and that in the side walls of the hexagonal cells, in which it is probably larger, due to a smaller molecular tilt angle. Therefore, at the phase transition, when the walls disappear, a very slight decrease in mean layer spacing would result. Qualitatively, these features agree with observations.

One could try to orchestrate a more physically appealing model of how the parameters of the free energy might change with temperature. For instance, $`A_S`$ steadily increasing would have essentially the same result as $`w`$ steadily decreasing. More generally, any combination of $`K_2q^2+P_0^2/(8\pi ϵ)`$ increasing and/or $`f_W`$ decreasing with falling temperature could produce the qualitative features of the apparent Sm–A phase and the phase transition.

The most valuable test for the validity of this model in real materials would be the direct observation of the periodic structure of the modulated phase. The first appropriate sub-optical probe that comes to mind is small angle X-ray diffraction, especially resonant X-ray diffraction that is most sensitive to variations of molecular orientation in a structure resonant . Hiroshi Yokoyama HY  has pointed out to us that atomic force microscopy might also be a suitable probe, if the modulated phase structure is preserved in a thin layer deposited on a substrate.

Remarkably, in a recent paper bookshelf  on one of the materials whose properties stimulated this work, optical observations consistent with a modulated phase, just at the transition from the Sm–A to the Sm–$`C^{}`$ phase, are reported. In a thin layer of material in the “bookshelf” geometry, i.e., the smectic layers viewed edge-on, thin parallel stripes are observed running perpendicular to the smectic layers. These stripes are possibly consistent with the long thin hexagonal cells that we propose here.

## Acknowledgments

We thank R. Pindak for helpful discussions. This work was supported by the National Science Foundation under grants DMR–9873849 and DMR–9974388.


text:
# Convective Heat Transport in Compressible Fluids

## I Introduction

Recently much attention has been paid to organized fluid motion in turbulent convection in the Rayleigh-B$`\stackrel{´}{\mathrm{e}}`$nard geometry . Although the conventional hydrodynamic equations are constructed for (nearly) incompressible fluids , we may mention a number of convection experiments in compressible one-component fluids in the supercritical region , together with those in non-critical fluids such as water or Hg . In these studies the Nusselt number $`Nu`$ representing the efficiency of convective heat transport has been measured at large values of the Rayleigh number $`Ra`$ defined by

$$Ra=\alpha _p\rho gL^3\mathrm{\Delta }T/\eta D.$$
(1.1)
Here $`g`$ is the gravity constant, $`\mathrm{\Delta }T=T_{\mathrm{bot}}T_{\mathrm{top}}`$ is the difference between the bottom and top temperatures, and $`L`$ is the cell height. As the critical point is approached in one-component fluids, the thermal expansion coefficient $`\alpha _p=(\rho /T)_p/\rho `$ grows strongly as $`\xi ^{\gamma /\nu }`$ (in the same manner as the isothermal expansion coefficient $`K_T=(\rho /p)_T/\rho `$ and the isobaric specific heat $`C_p`$), the thermal diffusivity $`D`$ decreases as $`\xi ^1`$, and the shear viscosity $`\eta `$ is nearly a constant. Here $`\xi `$ is the thermal correlation length growing as $`(T/T_c1)^\nu `$ on the critical isochore with $`\gamma 1.24`$ and $`\nu 0.625`$. Hence, in the critical region, $`Ra`$ can be extremely large; for example, $`Ra10^{13}`$ even for not very long $`L`$($`10`$ cm). The Prandtl number $`Pr=\eta /\rho D`$ was in the range of 1-100.

High compressibility of supercritical fluids gives rise to some unique features not encountered in incompressible fluids. (i) First, the transient behavior after application of a heat flux from the bottom is strongly influenced by the so-called piston effect , as revealed by recent high-precision experiments on 3He  and reproduced by subsequent simulation . (ii) Second, as $`\alpha _p`$ grows, the usual mechanism of convection onset $`Ra>Ra_\mathrm{c}(1708)`$ is replaced by that of the Schwarzschild criterion . That is, for large compressible fluid columns (even far from the critical point), convection sets in when thermal plumes continue to rise upward adiabatically. This occurs when the applied temperature gradient $`|dT/dz|`$ is larger than the adiabatic gradient ,

$$a_\mathrm{g}=(T/p)_s\rho g,$$
(1.2)
which is equal to $`0.034`$ mK$`/`$cm for 3He and $`0.27`$ mK$`/`$cm for CO2. This is the condition that the entropy $`s=s(T,p)`$ per unit mass decreases with height as $`ds/dz=(C_p/T)[dT/dz+a_\mathrm{g}]<0`$, under which the entropy of fluid elements adiabatically convected upward is larger than that of the ambient fluid. More precisely, Gitterman and Steinberg  found that the convection onset for compressible fluids is given by $`Ra^{\mathrm{corr}}>Ra_\mathrm{c}`$, where $`Ra^{\mathrm{corr}}`$ is a corrected Rayleigh number defined by

$$\begin{array}{ccc}\hfill Ra^{\mathrm{corr}}& =& (\alpha _p\rho gL^3/\eta D)(\mathrm{\Delta }Ta_\mathrm{g}L)\hfill \\ & =& Ra(1a_\mathrm{g}L/\mathrm{\Delta }T).\hfill \end{array}$$
(1.3)
This is a natural consequence because the effective temperature gradient seen by the raising plumes is given by $`\mathrm{\Delta }T/La_\mathrm{g}`$. At the convection onset we thus have

$$(\mathrm{\Delta }T)_{\mathrm{on}}=a_\mathrm{g}L+Ra_\mathrm{c}D\eta _0/(g\rho \alpha _pL^3).$$
(1.4)
where the second term behaves as $`(T/T_c1)^{\gamma +\nu }L^3`$ and can be much smaller than the first term even for small $`L(1`$mm ) as $`TT_c`$. The relation (1.4) has been confirmed in SF6 , and in 3He . (iii) Third, in steady convective states, experimental curves of $`Ra(Nu1)`$ vs $`Ra^{\mathrm{corr}}`$ were collapsed onto a single universal curve for various densities above $`T_\mathrm{c}`$  and for various average reduced temperatures on the critical isochore . These empirical results are highly nontrivial, because $`Nu`$ can in principle depend on $`Ra`$, $`Pr`$, and $`a_\mathrm{g}L/\mathrm{\Delta }T`$, while $`Nu`$ is a function of $`Ra`$ and $`Pr`$ for incompressible fluids.

For various fluids under relatively large $`\mathrm{\Delta }Ta_\mathrm{g}L`$ (where $`Ra^{\mathrm{corr}}Ra`$), data of $`Nu`$ have been fitted to a simple scaling law,

$$NuRa^a.$$
(1.5)
The exponent $`a`$ has been in a range from 0.28 to 0.31 and, in particular, a theoretical value $`2/7`$  was generally consistent with data for $`Ra10^{12}`$ . Moreover, measurements of the patterns of isothermal surfaces  and the velocity  have been informative on plume motion and a large-scale circulating shear flow in small-aspect-ratio cells . Several authors have also performed numerical analysis of convection at large $`Ra`$ in two dimensions (2D)  and in three dimensions (3D) . Even in 2D salient features in the experiments have been reproduced. In these simulations, if the temperature is averaged over a long time, the temperature gradient is localized in thin boundary layers with thickness $`\mathrm{}_T`$ related to $`Nu`$ by

$$Nu=L/2\mathrm{}_T.$$
(1.6)
Both in 2D and in 3D (if visualized from side), the plumes tend to be connected from bottom to top for large $`Pr`$ because of slow thermal diffusion, while they become diffuse far from the boundaries for small $`Pr`$. In the 3D simulations with periodic or free-slip sidewalls , local boundary shear flows were observed between incoming plumes and outgoing networks of buoyant sheets in horizontal planes close to the boundaries.

In this paper we will derive and examine hydrodynamic equations for compressible fluids under gravity in the supercritical region, in which the oscillatory motion of sound has been averaged out . Since the time scale of convective motions is much longer than that of the acoustic wave $`t_{\mathrm{ac}}=L/c`$ (typically of order $`10^4`$s for $`L1`$ cm), such a description is convenient theoretically and is even indispensable for numerical analysis. Our dynamic equations are a natural generalization of the usual hydrodynamic equations . Our new predictions are unique particularly when the piston effect comes into play, as has been demonstrated in the previous simulation  relatively close to the convection onset. This paper will present 2D simulations of our hydrodynamic equations for much larger $`Ra^{\mathrm{corr}}`$ both in transient and (dynamical) steady states. Even in steady states, we will find some characteristic features of turbulent states, which have not been reported in the previous simulations , such as the logarithmic velocity profile of the velocity near the boundary  and random reversal of the large-scale circulating flow in small-aspect-ratio cells . We will also point out that individual arrivals of plumes at the boundaries cause global temperature fluctuations in the cell via the piston effect. The resultant noise level of the temperature fluctuations grows as the critical point is approached.

## II Theoretical Background

### A Hydrodynamic equations

We consider a supercritical fluid on the critical isochore in a cell with the bottom plate at $`z=0`$ and the top plate at $`z=L`$. The $`z`$ axis is taken in the upward direction and the total fluid volume is fixed at $`V`$. The temperature disturbance $`\delta T(𝒓,t)=T(𝒓,t)T_{\mathrm{top}}`$ measured from the temperature $`T_{\mathrm{top}}`$ at the top boundary is much smaller in magnitude than $`T_{\mathrm{top}}T_\mathrm{c}`$. Hereafter $`ϵ`$ will be used to denote the reduced temperature at the top boundary, which satisfies

$$ϵ=T_{\mathrm{top}}/T_\mathrm{c}1\mathrm{\Delta }T/T_c.$$
(2.1)
We assume that the gravity-induced density stratification is not too severe such that the thermodynamic derivatives are nearly homogeneous in the cell. This is satisfied when $`|\rho /\rho _c1|(\rho /p)_TgLϵ^\beta `$ with $`\beta 0.33`$ . This condition is rewritten as

$$ϵ^{\beta +\gamma }a_\mathrm{g}L/T_\mathrm{c}.$$
(2.2)
In the theoretical literature on convection , the top and bottom temperatures $`T_{\mathrm{bot}}`$ and $`T_{\mathrm{top}}`$ are constant parameters. However, in most of recent convection experiments, particularly in cryogenic ones, the heat flux at the bottom $`Q=\lambda (dT/dz)_{z=0}`$ and $`T_{\mathrm{top}}`$ have been fixed. The $`\lambda `$ is the thermal conductivity. Furthermore, if the top and bottom walls are made of a metal with high thermal conductivity, the boundary temperatures become homogeneous in the lateral directions (unless local temperature changes are too fast). Then $`T_{\mathrm{bot}}(t)`$ and hence $`\mathrm{\Delta }T(t)`$ are functions of time only. Metcalfe and Behringer  performed linear stability analysis of convection onset under this cryogenic boundary condition. In the nonlinear regime it is of great interest how the boundary layer thickness and the plume generation depend on the boundary condition.

In equilibrium the pressure gradient is given by $`\rho g\rho _\mathrm{c}g`$. In nonequilibrium we set

$$p(𝒓,t)=p_0\rho _\mathrm{c}gz+p_1(t)+p_{\mathrm{inh}}(𝒓,t),$$
(2.3)
where $`p_0`$ is a constant, $`p_1(t)`$ and $`p_{\mathrm{inh}}`$ are the homogeneous and inhomogeneous parts induced by $`\delta T`$, respectively. That is, we assume $`p_{\mathrm{inh}}=0`$, where $`\mathrm{}𝑑𝒓(\mathrm{})/V`$ represents the space average. Then $`p_1`$ is related to the space average of $`\delta T`$ by

$$p_1(t)=(p/T)_\rho \delta T(t),$$
(2.4)
which follows from the thermodynamic relation $`dp=(p/T)_\rho dT+(p/\rho )_Td\rho `$ and the condition that the space average of the density deviation vanishes ($`\delta \rho =0)`$. It is important that the combination $`p(𝒓,t)+\rho _\mathrm{c}gz`$ is nearly homogeneous or $`|p_1(t)||p_{\mathrm{inh}}(𝒓,t)|`$ for fluid motions much slower than the acoustic time $`t_{\mathrm{ac}}=L/c`$ ($`c10^4`$ cm$`/`$s being the sound velocity)

Now we derive the equation for $`\delta T`$ from the heat conduction equation

$$\rho T\left(\frac{}{t}+𝒗\right)s=\lambda ^2\delta T,$$
(2.5)
where $`s(𝒓,t)`$ is the entropy per unit mass. Here $`\delta s`$ consists of the equilibrium part $`s_{\mathrm{eq}}(z)`$ with

$$\frac{d}{dz}s_{\mathrm{eq}}(z)=\left(\frac{s}{p}\right)_T\rho g=T^1C_pa_\mathrm{g}$$
(2.6)
and the nonequilibrium deviation,

$$\delta s(𝒓,t)=T^1C_p\left[\delta T(𝒓,t)\left(\frac{T}{p}\right)_sp_1(t)\right].$$
(2.7)
With the aid of the thermodynamic identity $`(T/p)_s=(T/p)_\rho (11/\gamma _s)`$, we rewrite (2.5) to obtain the desired equation for $`\delta T`$,

$$\left(\frac{}{t}+𝒗D^2\right)\delta T=a_\mathrm{g}v_z+\alpha _s\frac{d}{dt}\delta T,$$
(2.8)
where $`D=\lambda /C_p`$ is the thermal diffusivity and

$$\alpha _s=1\gamma _s^1.$$
(2.9)
The specific heat ratio $`\gamma _s`$ behaves as

$$\gamma _s=C_p/C_Vϵ^{\gamma +\alpha }1,$$
(2.10)
where $`C_pϵ^\gamma `$ and $`C_Vϵ^\alpha `$ are the specific heats (per unit volume) at constant $`p`$ and $`V`$, respectively, with $`\alpha 0.1`$. The first term on the right hand side of (2.8) arises from $`ds_{\mathrm{eq}}/dz`$. Inside plumes the temperature is adiabatically cooled if they go upward ($`v_z>0`$), or adiabatically warmed if they go downward $`(v_z<0)`$. In this way this term suppresses upward motion of warmer plumes from the bottom and downward motion of cooler plumes from the top, resulting in the Schwarzschild criterion of convection onset (the adiabatic temperature gradient effect). On the other hand, the second term arises from $`p_1(t)`$, leading to the piston effect . It is worth noting that the space integral of (2.8) in the cell becomes

$$VC_V\frac{d}{dt}\delta T=\lambda 𝑑a𝒏\delta T,$$
(2.11)
where use has been made of $`𝒗=\mathrm{𝟎}`$. The right hand side represents the rate of heat supply from the boundary surface, where $`da`$ is surface element and $`𝒏`$ is the outward surface normal. Its time-integration is the total heat supply expressed as $`V\rho \delta s`$, resulting in

$$C_V\delta T(t)=\rho \delta s(t),$$
(2.12)
which also follows (2.4) and the space average of (2.7). The appearance of $`C_V`$ in on the left hand side of (2.12) is a natural consequence under the fixed volume condition. Notice that (2.7) can also be written as

$$\delta T(𝒓,t)=\frac{\rho T}{C_p}\delta s(𝒓,t)+\rho T\left[\frac{1}{C_V}\frac{1}{C_p}\right]\delta s(t).$$
(2.13)
This relation holds even in gravity if $`\delta s`$ is the deviation of $`ss_{\mathrm{eq}}(z)`$ as in (2.7). In addition, the density deviation $`\delta \rho =\rho \rho `$ is written in our approximation as

$$\delta \rho =\rho K_Tg(zL/2)\rho \alpha _p(\delta T\delta T),$$
(2.14)
where $`K_T=(\rho /p)_T/\rho `$ and we have set $`\delta \rho =0`$.

Since $`C_pC_V`$ near the critical point, the homogeneous part of $`\delta T`$ (second term) in (2.13) can easily dominate over the inhomogeneous part (first term) even when $`\delta s`$ is localized near a heated wall. Indeed, if a thermal disturbance is produced within a thermal boundary layer with thickness $`\mathrm{}`$ near the boundary, the ratio of the homogeneous part ($`\delta s`$) to the localized inhomogeneous part ($`\delta s`$) in (2.13) is of order $`(\gamma _s1)\mathrm{}/L`$ where $`L`$ is the characteristic system length. Temperature homogenization is achieved when $`(\gamma _s1)\mathrm{}L`$. By setting $`\mathrm{}=(Dt_1)^{1/2}`$ we obtain the time constant of this thermal equilibration (the piston time) in the form,

$$t_1=L^2/D(\gamma _s1)^2.$$
(2.15)
Next we consider the momentum equation for the velocity field $`𝒗(𝒓,t)`$. On long time scales, sound waves decay to zero and the incompressibility condition

$$𝒗=0$$
(2.16)
becomes nearly satisfied ($`t_{\mathrm{homo}}`$) . Then the dissipation of $`𝒗`$ is produced by the shear viscosity $`\eta `$ and the usual Navier-Stokes equation in the Boussinesq approximation may be set up in the form ,

$$\left(\frac{}{t}+𝒗\right)𝒗=\frac{p_{\mathrm{inh}}}{\rho }+\alpha _pg\delta T𝒆_z+\frac{\eta }{\rho }^2𝒗,$$
(2.17)
where the inhomogeneous part $`p_{\mathrm{inh}}`$ ensures (2.16), $`𝒆_z`$ is the unit vector along the $`z`$ axis, and $`\rho `$ is the average density. The two equations (2.8) and (2.17) are our fundamental dynamic equations closed under (2.16). In the conventional theory , (2.17) has been used, but the right hand side of (2.8) vanishes.

As another characteristic feature near the critical point, the Prandtl number behaves as

$$Pr=\eta /\rho Dϵ^\nu .$$
(2.18)
For example, $`Pr=350`$ at $`T/T_\mathrm{c}1=10^3`$ in 3He. This means that the time scale of the thermal diffusion is much slower than that of the velocity in the critical region. Based on this fact, the simulation in Ref.  was performed using the Stokes approximation in which the left hand side of (2.17) is set equal to zero. Good agreement with the experiments  was then obtained for $`Ra^{\mathrm{corr}}/Ra_\mathrm{c}15`$ at $`ϵ=0.05`$.

For $`Pr1`$, let us estimate the upper bound of $`Ra^{\mathrm{corr}}`$ below which the Reynolds number $`Re`$ is smaller than 1 or the Stokes approximation is allowable. The characteristic temperature variation $`(\delta T)_{}`$ changing perpendicularly to the $`z`$ axis and the characteristic velocity field $`v_{\mathrm{pl}}`$ are related by

$$v_{\mathrm{pl}}(\alpha _p\rho _\mathrm{c}g/\eta k^2)(\delta T)_{},$$
(2.19)
where $`k2\pi /L`$ for roll patterns. If $`Ra^{\mathrm{corr}}/Ra_\mathrm{c}`$ is considerably (but not much) larger than 1, $`(\delta T)_{}/\mathrm{\Delta }T`$ is of order 1 (but somewhat smaller than 1). Then we obtain

$$v_{\mathrm{pl}}(Ra^{\mathrm{corr}}/Ra_\mathrm{c})D/L.$$
(2.20)
Thus the small Reynolds number regime is written as

$$Ra^{\mathrm{corr}}/Ra_\mathrm{c}Pr,$$
(2.21)
where use has been made of $`Rev_{\mathrm{pl}}L\rho /\eta `$. For $`Pr1`$ there is a sizable range of $`Ra^{\mathrm{corr}}`$ in which the Stokes approximation is justified. In passing, for $`0<Ra^{\mathrm{corr}}/Ra_\mathrm{c}11`$, the theory of the amplitude equation  predicts

$$v_{\mathrm{pl}}L/D(\delta T)_{}/\mathrm{\Delta }T(Ra^{\mathrm{corr}}/Ra_\mathrm{c}1)^{1/2},$$
(2.22)
from which we have $`Nu1Ra^{\mathrm{corr}}/Ra_\mathrm{c}1`$ because the convective heat current is of order $`C_p(\delta T)_{}v_{\mathrm{pl}}`$. In the next section we will estimate $`v_{\mathrm{pl}}`$ for much larger $`Ra`$.

Analogously to (2.19), the inhomogeneous pressure deviation $`p_{\mathrm{inh}}`$ is estimated as $`p_{\mathrm{inh}}(\alpha _p\rho _\mathrm{c}g/k)(\delta T)_{}`$. If we assume $`p_1(t)(p/T)_s\mathrm{\Delta }T`$ from (2.4) and $`\mathrm{\Delta }T(\delta T)_{}`$ as in (2.20), we find that $`p_{\mathrm{inh}}/p_1(t)`$ is of order $`ϵ^\gamma a_\mathrm{g}/T_ck`$ and is much smaller than 1 from (2.2). This estimation justifies the assumption of the homogeneity of $`\delta p(𝒓,t)+\rho _cgz`$ made below (2.4).

### B Free energy and heat production rate

In the presence of small deviations of the temperature and the density, $`\delta T`$ and $`\delta \rho `$, around an reference equilibrium state, we have an increase of the free energy functional $`\delta F`$. Up to the bilinear order of the deviations, it is of the form ,

$$\delta F=𝑑𝒓\left[\frac{C_V}{2T}(\delta T)^2+\frac{1}{2\rho ^2K_T}(\delta \rho )^2+gz\delta \rho \right],$$
(2.23)
where the third term is the potential energy in gravity. All the deviations are assumed to change slowly in space compared with the thermal correlation length $`\xi `$. If we express $`\delta \rho `$ in terms of $`\delta T`$ as in (2.14), we obtain

$$\delta F=\frac{1}{2T}𝑑𝒓\left[C_p(\delta T\delta T)^2+C_V\delta T^2\right],$$
(2.24)
where the constant term is omitted. We notice that $`\delta F`$ decreases dramatically for $`\gamma _s1`$ in the process of adiabatic temperature homogenization. Furthermore, in the presence of velocity field, the total free energy change is the sum of $`\delta F`$ and the kinetic energy of the velocity field,

$$F_\mathrm{K}=\frac{1}{2}𝑑𝒓\rho 𝒗^2.$$
(2.25)
Its time derivative is calculated from our dynamic equations (2.7) and (2.17) in the form,

$$\begin{array}{ccc}\hfill \frac{d}{dt}(\delta F+F_\mathrm{K})& =& 𝑑𝒓(ϵ_{\mathrm{th}}+ϵ_{\mathrm{vis}})\hfill \\ & +& \lambda T^1𝑑a[\delta T(𝒏\delta T)],\hfill \end{array}$$
(2.26)
where $`ϵ_{\mathrm{th}}`$ and $`ϵ_{\mathrm{vis}}`$ are the thermal and viscous heat production rates (per unit volume) , respectively, defined by

$$ϵ_{\mathrm{th}}=\lambda T^1|\delta T|^2,$$
(2.27)
$$ϵ_{\mathrm{vis}}=\eta \underset{ij}{}(v_i/x_j)^2.$$
(2.28)
In the second term of (2.26) the surface integral is over the boundary of the cell, $`𝒏`$ being the outward unit vector. In terms of the heat flux from the bottom $`Q`$, it is expressed as $`VQ\mathrm{\Delta }T/TL`$ if the top temperature is fixed.

### C Basic relations in steady states

We consider steady convective states in the Rayleigh-B$`\stackrel{´}{\mathrm{e}}`$nard geometry, in which the flow pattern is either time-independent not far above the convection onset or chaotic at larger $`Ra`$. We treat $`\mathrm{\Delta }T`$ as a constant parameter. Under the condition of fixed heat flux at the bottom, however, $`\mathrm{\Delta }T(t)`$ exhibits rapidly-varying fluctuations in chaotic states. In this case $`\mathrm{\Delta }T`$ in the following relations represents the time-average of $`\mathrm{\Delta }T(t)`$. The steady state averages (over space and time) will be denoted by $`\mathrm{}_\mathrm{s}`$ to distinguish them from the space averages $`\mathrm{}`$ used so far.

We make (2.8) and (2.17) dimensionless by measuring space and time in units of $`L`$ and $`L^2/D`$ and setting $`\stackrel{~}{𝒓}=L^1𝒓`$ and $`\stackrel{~}{t}=DL^2t`$. The temperature deviation is written as

$$\delta T(𝒓,t)/\mathrm{\Delta }T=1\stackrel{~}{z}+Ra^1(\stackrel{~}{𝒓},\stackrel{~}{t}),$$
(2.29)
where $`\stackrel{~}{z}=z/L`$. The dimensionless function $``$ becomes nonvanishing in convective states and obeys

$$\left(\frac{}{\stackrel{~}{t}}+𝑽\stackrel{~}{}\stackrel{~}{}^2\right)=Ra^{\mathrm{corr}}V_z+\alpha _s\frac{d}{d\stackrel{~}{t}},$$
(2.30)
where $`\stackrel{~}{}=L`$ is the space derivative in units of $`L`$. Then the (average) heat flux at the bottom is written as $`Q=(\lambda \mathrm{\Delta }T/L)[1+Ra^1f_\lambda ]`$, where

$$f_\lambda =(/\stackrel{~}{z})_{\stackrel{~}{z}=0}_\mathrm{s}.$$
(2.31)
The $`f_\lambda `$ is a function of $`Ra^{\mathrm{corr}}`$ and $`Pr`$. The Nusselt number $`Nu=QL/\lambda \mathrm{\Delta }T`$ is expressed as

$$Nu=1+Ra^1f_\lambda .$$
(2.32)
As the boundary condition of $``$ we require $`=0`$ at $`\stackrel{~}{z}=0`$ and $`1`$ if $`T_{\mathrm{top}}`$ and $`T_{\mathrm{bot}}`$ are fixed. However, If $`T_{\mathrm{top}}`$ and $`Q`$ at the bottom are fixed, we have $`=0`$ at $`\stackrel{~}{z}=0`$ and $`/\stackrel{~}{z}=Ra(Nu1)`$ at $`\stackrel{~}{z}=0`$. The dimensionless velocity $`𝑽(\stackrel{~}{𝒓},\stackrel{~}{t})=(L/D)𝒗`$ obeys

$$\frac{1}{Pr}\left(\frac{}{\stackrel{~}{t}}+𝑽\stackrel{~}{}\right)𝑽=\stackrel{~}{}P_{\mathrm{inh}}+𝒆_z+\stackrel{~}{}^2𝑽,$$
(2.33)
where $`P_{\mathrm{inh}}`$ ensures $`\stackrel{~}{}𝑽=0`$.

Here we assume that the piston term, the second term on the right hand side of (2.30), can be neglected in steady states. For $`ϵ=0.05`$, the piston term in steady states is less than a few percents of the convection term $`𝒗\stackrel{~}{}`$ in (2.30) except at the boundaries. It thus produces no significant effects on steady state heat transport (on $`Nu`$), while it can be crucial in the initial transient stage . Then, if the piston term in (2.30) is neglected, (2.30) and (2.33) become of the same form as those of usual incompressible fluids except that $`Ra^{\mathrm{corr}}`$ appears in place of $`Ra`$. At much smaller $`ϵ`$, however, this assumption is questionable, because the noise part of $``$ grows as $`ϵ0`$, as will be discussed later in the next section. We may conclude the following (at least for $`ϵ=0.05`$). (i) It follows the Gitterman-Steinberg criterion $`Ra^{\mathrm{corr}}>Ra_\mathrm{c}`$ in convective states in the compressible case . (ii) It is more nontrivial that the combination

$$Ra(Nu1)=f_\lambda (Ra^{\mathrm{corr}},Pr)$$
(2.34)
should be a universal function of $`Ra^{\mathrm{corr}}`$ and $`Pr`$ from (2.32) in agreement with the experiments . Notice that $`Ra(Nu1)=f_\lambda (Ra,Pr)`$ holds for incompressible fluids in terms of the same $`f_\lambda `$. These experiments and more decisively that by Ahlers and Xu  indicate that $`f_\lambda `$ should be nearly independent of $`Pr`$ once $`Pr`$ considerably exceeds 1. In the 3D simulation by Verzicco and Camussi , $`Nu`$ became independent of $`Pr`$ for $`Pr0.5`$. Theoretical support of this behavior using scaling arguments was presented in Ref..

In steady states we may also derive some exact relations for variances among $`\delta T`$ and $`𝒗`$. Using the dynamic equations (2.8) and (2.17) we calculate the averages of $`(\delta T)^2/t`$, $`𝒗^2/t`$, and $`(z\delta T)/t`$ to obtain

$$|\delta T|^2_\mathrm{s}=a_{\mathrm{th}}^2+a_{\mathrm{th}}(a_{\mathrm{th}}a_\mathrm{g})(Nu1),$$
(2.35)
$$\underset{ij}{}(v_i/x_j)^2_\mathrm{s}=Ra(D/L^2)^2(Nu1).$$
(2.36)
We also obtain a cross correlation,

$$v_z\delta T_\mathrm{s}=a_{\mathrm{th}}D(Nu1),$$
(2.37)
which is nothing but the average convective heat flux (if $`\lambda C_p`$ is multiplied). Here $`a_{\mathrm{th}}\mathrm{\Delta }T/L=d\delta T/dz_\mathrm{s}`$ is the average temperature gradient and $`a_\mathrm{g}`$ is the adiabatic temperature gradient defined by (1.2). If we use the usual hydrodynamic equations for incompressible fluids, the right hand side of (2.35) becomes $`a_{\mathrm{th}}^2Nu`$, while (2.36) and (2.37) remain the same . In addition, (2.35) indicates $`a_{\mathrm{th}}>a_\mathrm{g}`$ in convective states in which $`Nu>1`$. This is consistent with the convection criterion $`Ra^{\mathrm{corr}}>Ra_\mathrm{c}`$. We obtain the averages of the two dissipation rates in (2.27) and (2.28) by multiplying $`\lambda /T`$ and $`\eta `$ to (2.35) and (2.36), respectively. Using the thermodynamic identity $`T\alpha _p=C_p(T/p)_s`$, we obtain

$$ϵ_{\mathrm{th}}_\mathrm{s}+ϵ_{\mathrm{vis}}_\mathrm{s}=T^1\lambda a_{\mathrm{th}}^2Nu,$$
(2.38)
$$(ϵ_{\mathrm{th}}_\mathrm{s}T^1\lambda a_{\mathrm{th}}^2)/ϵ_{\mathrm{vis}}_\mathrm{s}=a_{\mathrm{th}}/a_\mathrm{g}1,$$
(2.39)
The first relation (2.38) also follows from the average of (2.26). The second relation (2.39) holds only in convective states ($`Nu>1`$), while the right hand side is replaced by $`C_pa_{\mathrm{th}}/T\alpha _p\rho g=a_{\mathrm{th}}/a_\mathrm{g}`$ for the usual hydrodynamic equations of incompressible fluids.

## III Simulation results

We perform numerical analysis of (2.8) and (2.17) in 2D using parameters of 3He in a cell with $`L=1.06`$ mm. The reduced temperature is $`ϵ=0.05`$ (except in Fig.13), where $`\gamma _s=22.8,T\alpha _p=26.9`$, $`\lambda =1.88\times 10^4`$ erg/(cm2s K), $`D=5.42\times 10^5`$ cm$`{}_{}{}^{2}/`$s, and $`Pr=7.4`$ . The condition (2.2) is well satisfied. The piston time $`t_1`$ in (2.15) is given by $`0.42`$ s. We apply a constant heat flux $`Q`$ at the bottom $`z=0`$ for $`t>0`$ with a fixed top temperature $`T_{\mathrm{top}}`$ at $`z=L`$. In steady states we have $`Ra^{\mathrm{corr}}/Ra_\mathrm{c}=0.90[\mathrm{\Delta }T/a_\mathrm{g}L1],`$ where $`a_\mathrm{g}L=3.57`$ $`\mu `$K. Thus $`(\mathrm{\Delta }T)_{\mathrm{on}}=7.6`$ $`\mu `$K and $`Q_{\mathrm{on}}=13.5`$ nW$`/`$s at the convection onset. We assume homogeneity of the boundary temperatures, $`T_{\mathrm{top}}`$ and $`T_{\mathrm{bot}}`$, in the lateral $`x`$ direction.

In the experiments the aspect ratio was 57, so in the simulation  the periodic boundary condition was imposed in the $`x`$ direction with period $`4L`$. This period was chosen because the roll period is close to $`2L`$ slightly above the onset for infinite lateral dimension . Then, in steady states in the region $`1<Q/Q_{\mathrm{on}}5`$, the linear relation

$$Q/Q_{\mathrm{on}}1A_0[\mathrm{\Delta }T/(\mathrm{\Delta }T)_{\mathrm{on}}1]$$
(3.1)
was numerically obtained with $`A_02.2`$ in good agreement with the experiments. From $`Nu=[Q/\mathrm{\Delta }T]/[Q_{\mathrm{on}}/(\mathrm{\Delta }T)_{\mathrm{on}}]`$, the behavior of $`Nu`$ is known from (3.2) in the range $`1<Q/Q_{\mathrm{on}}5`$. In particular, slightly above the onset, we have

$$Nu1A_1(Ra^{\mathrm{corr}}/Ra_\mathrm{c}1)+\mathrm{}.$$
(3.2)
where $`A_10.64`$ in fair agreement with the theoretical value ($`A_10.70`$ for $`Pr=7.4`$) . This behavior is also consistent with (2.22).

In this work we are interested in fluid motion for relatively large $`Ra`$ up to $`3\times 10^6`$. In the following we show two sets of the numerical results. In the first set, periodic sidewalls are assumed at $`x=0`$ and $`x=L_{}`$ with period $`L_{}=4L`$ as in Ref.. In Table 1 the steady state values of $`\mathrm{\Delta }T`$, $`Ra^{\mathrm{corr}}`$, $`Ra`$, $`Nu`$, and $`\overline{R}e`$ are written, where $`\overline{R}e`$ is a Reynolds number to be defined in (3.11). They are obtained for $`Q=0.0458`$ $`\mu `$W$`/`$cm$`{}_{}{}^{2}(3.4Q_{\mathrm{on}}`$), $`0.965`$$`\mu `$W$`/`$cm$`{}_{}{}^{2}(71Q_{\mathrm{on}}`$), and $`122.2`$$`\mu `$W$`/`$cm$`{}_{}{}^{2}(9\times 10^3Q_{\mathrm{on}}`$). For the smallest $`Q`$ the system tends to a time-independent convective state, as already studied in Ref., while for the other values of $`Q`$ the system tends to a chaotic state without macroscopic boundary shear flow. In the second set, we perform simulations for $`A=1,2,`$ and 3 with insulating and rigid sidewalls at $`x=0`$ and $`AL`$, at which $`𝒗=0`$ and through which there is no heat flux ($`\delta T/x=0`$), as will be presented in Figs.4, 12, and 13.

In addition, if the temperature difference will be simply written as $`\mathrm{\Delta }T`$, it should be taken as the time average of $`\mathrm{\Delta }T(t)`$ in a steady state. We also assume that $`Pr`$ is considerably larger than 1 in the following arguments.

### A Transient behavior

In Fig.1 we show numerically calculated $`\mathrm{\Delta }T(t)=T_{\mathrm{bot}}(t)T_{\mathrm{top}}`$ for $`Q=0.965`$$`\mu `$W$`/`$cm2 in (a) and for $`Q=122.2`$$`\mu `$W$`/`$cm2 in (b). They nearly coincide with the upper broken curve without convection ($`𝒗=\mathrm{𝟎}`$) in the initial stage before the maximum is attained. The latter curve is calculated from (2.8) as

$$[\mathrm{\Delta }T(t)]_0=\frac{Q}{\lambda }\sqrt{\frac{Dt}{\pi }}\left[4_0^{\mathrm{}}\frac{ds}{\sqrt{\pi s}}\frac{1e^s}{s+t/t_1}\right],$$
(3.3)
where $`t_1`$ is defined by (2.15) and the integral in the brackets behaves as $`(\pi t_1/t)^{1/2}`$ for $`tt_1`$ . If the piston term is absent and $`𝒗=\mathrm{𝟎}`$, (2.8) becomes the simple diffusion equation, yielding $`[\mathrm{\Delta }T(t)]_0=(2Q/\lambda )(Dt/\pi )^{1/2}`$, which is about half of $`[\mathrm{\Delta }T(t)]_0`$ in (3.3) for $`tt_1`$ (see Fig.3 in Ref.). We also show the numerically calculated $`\mathrm{\Delta }T(t)`$ at fixed pressure where the piston term is absent ($`\alpha _s=0`$ in (2.8)) but $`𝒗\mathrm{𝟎}`$. In (a) the experimental curve is shown to have a lower peak and overdamp more slowly than in our simulation. In (b) the selected value of $`Q`$ is in the region where no overshoot was observed in the experiment. See also Fig.11, where the numerical curves of $`\mathrm{\Delta }T(t)`$ will be given for other choices of the parameters.

In Fig.2 we show time evolution of the temperature profile at $`Q=122.2`$ $`\mu `$W$`/`$cm for periodic sidewalls. In (A) and (B) small-scale mushroom-like plumes are ejected from the bottom. In (C) and (D) they reach the top and are flattened there. In this initial stage the typical raising speed $`v_{\mathrm{pl}}`$ is estimated as $`L/t_{\mathrm{tr}}`$ where $`t_{\mathrm{tr}}`$ is the traversing time. From (A)-(C) we find that it is nearly equal to the free-fall velocity $`v_\mathrm{g}`$ defined by

$$v_\mathrm{g}=(Lg\alpha _p\mathrm{\Delta }T)^{1/2}=(RaPr)^{1/2}D/L,$$
(3.4)
which is $`2.37`$ cm$`/`$s. In this case the plumes leave the bottom at zero velocity and go upward with their velocity roughly of the form,

$$v_{\mathrm{pl}}(t)=v_{\mathrm{}}\left[1\mathrm{exp}[(tt_0)/t_{\mathrm{vis}}]\right],$$
(3.5)
where $`t_0`$ is the departure time, $`t_{\mathrm{vis}}\rho R^2/\eta `$ is the viscous relaxation time with $`R`$ being the plume size, and

$$v_{\mathrm{}}R^2g\rho \alpha _p\mathrm{\Delta }T/\eta $$
(3.6)
is the terminal velocity achieved by balance between the buoyancy and the viscous drag. For $`t_{\mathrm{tr}}t_{\mathrm{vis}}`$ the viscous drag is negligible and we have $`v_{\mathrm{pl}}(t)v_\mathrm{g}^2(tt_0)/L`$ and $`t_{\mathrm{tr}}L/v_\mathrm{g}`$. Thus, if the initial velocity is much less than $`v_\mathrm{g}`$, the free-fall condition becomes

$$R/L(Pr/Ra)^{1/4},$$
(3.7)
under which $`v_{\mathrm{}}=(R/L)^2(Ra/Pr)^{1/2}v_\mathrm{g}v_\mathrm{g}`$. In Fig.2, $`R/L1/3`$ and $`(Pr/Ra)^{1/4}0.04`$, so the above condition is satisfied.

With the arrival of the plumes the heat current increases at the top, because $`T_{\mathrm{top}}`$ is fixed, and a negative deviation of $`\delta s`$ is produced in a layer near the top. As can be known from (2.13), the piston effect is then operative, resulting in a homogeneous lowering of the temperature in the whole cell. In the time region around (E) the fluid is vigorously mixed with high Reynolds numbers. More precisely, the height-dependent Reynolds number $`\widehat{R}e(z,t)`$ to be defined in (3.12) below is about $`20`$ except in the vicinity of the boundaries. A downward flow of cooler fluid regions is then produced from the top. In the steady state (F), the temperature deviation becomes considerably smaller than in the transient states, and the localized boundary shear flows are produced between outgoing and incoming plumes with thickness $`\mathrm{}_v`$ much smaller than $`L`$.

The overshoot is more clearly illustrated in Fig.3, which displays the average of $`\delta T(x,z,t)`$ taken in the $`x`$ direction,

$$\overline{\delta T}(z,t)_0^L_{}\frac{dx}{L_{}}\delta T(x,z,t),$$
(3.8)
for the points (A), (C), (E), and (F) in Fig.1b. As a characteristic feature, the temperature in the interior consists of global changes due to the piston effect and bumps due to localized plumes. In (E) the cooler layer becomes thicker temporarily near the top due to the excess heat flow.

In our simulation the raising plumes leave the bottom and reach the top nearly simultaneously, resulting in a homogeneous temperature change. (i) Not far above the onset this mechanism is the main cause of the overshoot in compressible fluids. Note that a small peak appears in $`\mathrm{\Delta }T(t)`$ even in the fixed pressure case ($`\alpha _s=1`$) as shown in Fig.2 of Ref. and as was observed by Behringer and Ahlers . Furthermore, in Ref., the time scale of the overshoot (from the maximum to the minimum of $`\mathrm{\Delta }T(t)`$) due to the piston effect was predicted to be of order $`t_D/(Ra^{\mathrm{corr}}/Ra_\mathrm{c}1)`$, where $`t_D=L^2/4D`$($`50`$ s) is the diffusion time. This fairly agrees with later analysis of the experimental data . (ii) For much larger $`Q`$ such as those in Figs.1a and 1b, however, the downward flow from the top is also rapid enough to produce large overshoot, as demonstrated by the curves at fixed (height-dependent) pressure. Whether fixed is the volume or the pressure, the time scale of the overshoot is of the order of the traversing time $`L/v_\mathrm{g}`$ of the plumes due to gravity.

As regards the overshoot behavior of $`\mathrm{\Delta }T(t)`$, agreement between our simulation and the experiment  becomes worse with increasing $`Q`$. We point out the possibility that in the experiment a synchronous arrival of plumes at the top might have not been realized for very large $`Q`$ or for very short $`L/v_\mathrm{g}`$ because of large lateral dimensions of the cell used. That is, if some plumes arrive at the top and others leave the bottom at the same time, negative interference between currents up and down will suppress overshoot.

### B Steady state behavior

Now we discuss the Nusselt number $`Nu`$ in steady states. Fig.4 shows the combination $`Ra^{\mathrm{corr}}(Nu1)/(Ra^{\mathrm{corr}}Ra_\mathrm{c})`$ vs $`Ra^{\mathrm{corr}}/Ra_\mathrm{c}1`$ for periodic sidewalls and for $`A=1,2,`$ and 3. This combination depends on $`Ra^{\mathrm{corr}}`$ and $`A`$ from (2.34) in steady states. The data (solid line)  excellently agree with the numerical results for periodic sidewalls. We find that the scaling relation (1.5) nicely holds for $`Ra^{\mathrm{corr}}/Ra_\mathrm{c}10`$ for periodic sidewalls, while it holds only for $`Ra^{\mathrm{corr}}/Ra_\mathrm{c}10^3`$ at $`A=1`$. The exponent $`a`$ in (1.5) is close to $`2/7`$, but $`a=1/4`$ is also consistent with our numerical data. If $`A1`$ and $`Ra^{\mathrm{corr}}`$ is not very large such that the plume size is of order $`L`$, large-scale fluid motions are suppressed by the rigid sidewalls. This marked tendency of the $`A`$ dependent crossover of $`Nu`$ was already reported in measurements for $`A=0.5,1`$, and $`6.7`$ . It was also confirmed in the 3D simulation by Kerr  for periodic sidewalls.

In Fig.5 we show the steady-state temperature deviation $`\overline{\delta T}(z)`$ averaged in the $`x`$ direction as in (3.8) and in time for the three values of $`Q`$ in Table 1 for periodic sidewalls with period $`L_{}=4L`$. The averages taken along the $`x`$ direction become only weakly fluctuating in time in steady chaotic states (the relative fluctuations being of order $`10\%`$ for the largest $`Q`$). As has been observed ubiquitously in the previous simulations, the temperature gradient becomes localized within thermal boundary layers with thickness $`\mathrm{}_T`$. Because $`\mathrm{\Delta }T2\mathrm{}_TQ/\lambda `$ for $`\mathrm{}_TL`$, it is related to $`Nu`$ by (1.6). The arrows in Fig.5 represent the maximum points, $`z=\mathrm{}_v`$ and $`L\mathrm{}_v`$, of the variance of the horizontal velocity defined by

$$v_x^{}(z)=\left[_0^L_{}\frac{dx}{L_{}}v_x(x,z,t)^2\right]^{1/2}.$$
(3.9)
In Fig.6 we plot the normalized velocity variances, $`v_x^{}(z)/v_\mathrm{g}`$ in (a) and $`v_z^{}(z)/v_\mathrm{g}`$ in (b), where $`v_\mathrm{g}`$ is defined by (3.4) and

$$v_z^{}(z)=\left[_0^L_{}\frac{dx}{L_{}}v_z(x,z,t)^2\right]^{1/2}.$$
(3.10)
The time average of $`v_x^2`$ and $`v_z^2`$ in the brackets is also taken in these figures. On one hand, $`v_x^{}`$ take maxima at $`z=\mathrm{}_v`$ and $`L\mathrm{}_v`$, where $`\mathrm{}_v`$ is hardly distinguishable from $`\mathrm{}_T`$. On the other hand, $`v_z^{}`$ is largest at the middle of the cell. We also find that the sum (the kinetic-energy variance) $`(v_x^{})^2+(v_z^{})^2`$ is nearly constant in the interior, which was a finding reported in Ref.. At large $`Ra`$ the maxima of $`v_x^{}`$ and $`v_z^{}`$ are of the same order and will be identified as the typical plume velocity $`v_{\mathrm{pl}}`$. In our simulation we have $`v_{\mathrm{pl}}0.1v_\mathrm{g}(Ra^{1/2})`$, which is consistent with velocity measurements .

Kerr and Herring  made similar plots of the height-dependent velocity variances in their 3D simulations for free-slip sidewalls. They found that the characteristic length $`\mathrm{}_v`$ defined by the peak positions of $`v_x^{}(z)`$ becomes longer than $`\mathrm{}_T=L/2Nu`$ with increasing $`Ra`$; for example, for $`Pr=7`$ they obtained $`\mathrm{}_v/\mathrm{}_T1`$ at $`Ra=10^4`$ and $`\mathrm{}_v/\mathrm{}_T3`$ at $`Ra=10^7`$. Verzicco and Camussi obtained a similar slow growing of $`\mathrm{}_v/\mathrm{}_T`$ at large $`Ra`$ for $`Pr>1`$ in their 3D simulation with $`A=1`$ . Also similarly, our 2D simulation with $`Pr=7.4`$ gives $`\mathrm{}_v/\mathrm{}_T=2.54`$ and $`1.1`$ for $`Q=122.2`$ and $`0.965`$ $`\mu `$W$`/`$cm, respectively, but we cannot draw a definite conclusion because of our limited range of $`Ra`$.

In Fig.7 we plot an overall Reynolds number $`\overline{R}e`$ vs $`R^{\mathrm{corr}}/R_\mathrm{c}1`$ in the simulation for periodic sidewalls. It is defined by

$$\overline{R}e=\frac{\rho }{\eta }\left[|𝒗𝒗|^2/|^2𝒗|^2\right]^{1/2},$$
(3.11)
where the averages are taken in the whole space region. The $`\overline{R}e`$ is smaller than 1 for $`R^{\mathrm{corr}}/R_\mathrm{c}5`$ . For larger values of $`R^{\mathrm{corr}}`$, it exceeds 1 and the effective exponent $`(\mathrm{ln}\overline{R}e)/(\mathrm{ln}Ra^{\mathrm{corr}})`$ is from $`1/4`$ to $`1/3`$. However, as suggested by Fig.6, the strength of the velocity fluctuations strongly depends on the distance from the boundary, so it is more informative to introduce a height-dependent Reynolds number,

$$\widehat{R}e(z)=\frac{\rho }{\eta }\left[_0^L_{}𝑑x|𝒗𝒗|^2/_0^L_{}𝑑x|^2𝒗|^2\right]^{1/2},$$
(3.12)
where the time averages of the integrands are taken. As shown in Fig.8, $`\widehat{R}e(z)`$ takes maxima at $`z\mathrm{}_v`$ and $`L\mathrm{}_v`$ of order

$$\widehat{R}e(\mathrm{}_v)\mathrm{}_vv_{\mathrm{pl}}\rho /\eta ,$$
(3.13)
where $`v_{\mathrm{pl}}v_x^{}(\mathrm{}_v)`$. This relation indicates $`\widehat{R}e(\mathrm{}_v)Ra^{1/2a}`$ with $`a2/7`$ from $`v_{\mathrm{pl}}0.1v_\mathrm{g}`$ and $`\mathrm{}_v\mathrm{}_T`$. The $`\widehat{R}e(z)`$ becomes considerably smaller in the interior than at $`z\mathrm{}_v`$, whose origin is the sparseness of the plumes in the interior (see (3.20) below). We confirm that $`\overline{Re}`$ is of the order of the space average $`_0^L𝑑z\widehat{R}e(z)/L`$. In the literature , however, the (large-scale) Reynolds number has been identified as $`Re=v_{\mathrm{pl}}L\rho /\eta `$, which is much larger than $`\widehat{R}e(\mathrm{}_v)`$ in (3.13) by $`L/\mathrm{}_v`$. (For roll patterns, as was discussed below (2.21), we uniquely have $`Re=v_{\mathrm{pl}}L\rho /\eta `$.)

At very large $`Ra`$ the boundary layers should gradually crossover from a laminar state to a turbulent state except within thin viscous sublayers with thickness $`z_0`$ much shorter than $`\mathrm{}_v`$. In the inertial region $`z_0z\mathrm{}_v`$ of the boundary layer, it is natural to expect the logarithmic velocity profile ,

$$v_x^{}(z)=b_0^1(\sigma _0/\rho )^{1/2}[\mathrm{ln}(z/z_0)+c_0],$$
(3.14)
where $`\sigma _0`$ is the amplitude of the shear stress at the boundary with $`b_0`$ and $`c_0`$ being dimensionless numbers of order 1. We may set $`\sigma _0=\eta lim_{z0}D_{xz}(z)`$, where $`D_{xz}(z)`$ is the variance of the velocity gradient,

$$D_{xz}(z)=\left[_0^L_{}\frac{dx}{L_{}}\left(\frac{}{z}v_x(x,z,t)\right)^2\right]^{1/2}.$$
(3.15)
Then $`v_x^{}(z)(\sigma _0/\eta )z`$ as $`z0`$. It is appropriate to define $`z_0`$ by

$$z_0=\eta /(\rho \sigma _0)^{1/2},$$
(3.16)
which ensures $`\widehat{R}e(z_0)1`$. The size of $`\sigma _0`$ should be equal to the typical size of $`\rho v_xv_z`$ at $`z=\mathrm{}_v`$ even if we consider localized shear flows, so we also have

$$\sigma _0\rho v_{\mathrm{pl}}^2.$$
(3.17)
The ratio of the two lengths $`z_0`$ and $`\mathrm{}_v`$ is given by

$$\mathrm{}_v/z_0v_{\mathrm{pl}}\mathrm{}_v\rho /\eta \widehat{R}e(\mathrm{}_v),$$
(3.18)
which grows with increasing $`Ra`$. In Fig.9a, $`v_x^{}(z)`$ is fitted to the above logarithmic form in the inertial region for $`Q=122.2`$ $`\mu `$W$`/`$cm, where $`(\sigma _0/\rho )^{1/2}=0.067v_\mathrm{g}=0.16`$ cm$`/`$s, $`b_0=1.2,c_0=0.97,`$ and $`z_0=0.025L`$. In Fig.9b, we plot $`v_x^{}(z)`$ and $`zD_{xz}(z)`$ on a logarithmic scale. We may conclude that these quantities do not behave as $`z`$ in the inertial region of the boundary layers, although our $`Ra`$ is not large enough to unambiguously demonstrate the logarithmic velocity profile. Here we point out that our results are not consistent with Shiraiman and Siggia’s primary assumptions of $`\mathrm{}_T<\mathrm{}_v`$ and the linear profile of the mean shear flow, $`v_xz`$, in the region $`z<\mathrm{}_T`$ .

In contrast to the averages taken along the $`x`$ direction, those taken along the $`z`$ direction are rapidly varying functions of time at large $`Ra`$ due to the random plume motions. We consider the vertical velocity variance defined by

$$v_z^{}(x,t)=\left[_0^L\frac{dz}{L}v_z(x,z,t)^2\right]^{1/2}.$$
(3.19)
In Fig.10 we display snapshots of $`v_z^{}(x,t)`$, where the time average is not taken and peaks arising from the plumes become more apparent with increasing $`Q`$. For our $`Ra`$ realized, the space regions occupied by the plumes become more sparse with increasing $`Ra`$ in the interior. As the plumes move through the cell, they remain distinguishable from the ambient fluid because the thermal diffusion length $`(DL/v_{\mathrm{pl}})^{1/2}`$ does not much exceeds $`\mathrm{}_v`$. So we may define the volume fraction of the plumes $`\varphi _{\mathrm{pl}}`$. The convective heat current is of order $`\varphi _{\mathrm{pl}}v_{\mathrm{pl}}C_p\mathrm{\Delta }T\lambda Nu\mathrm{\Delta }T/L`$, leading to

$$\varphi _{\mathrm{pl}}D/\mathrm{}_Tv_{\mathrm{pl}},$$
(3.20)
which is of order $`Pr^1\widehat{R}e(\mathrm{}_v)^11`$ from (3.13). For much larger $`Ra`$, the plumes will generate smaller scale eddies, ultimately leading to fully developed turbulence in the interior, as will be discussed in Section 4.

### C Overall Temperature Fluctuations

When a plume with a volume $`V_{\mathrm{fl}}`$ reaches the boundary, it transfers a heat of order $`C_p\mathrm{\Delta }TV_{\mathrm{fl}}`$ to the boundary wall. As indicated by (2.13), the piston effect then gives rise to a homogeneous change in $`\delta T(t)`$ of order

$$(\delta T)_{\mathrm{fl}}\gamma _s(V_{\mathrm{fl}}/V)\mathrm{\Delta }T.$$
(3.21)
Of course, the real plumes are extended objects and are continuously arriving at the boundary in high $`Ra`$ convection. Thus $`V_{\mathrm{pl}}/V`$ in the above formula should be regarded as the fluctuation amplitude of the plume volume fraction $`\varphi _{\mathrm{pl}}`$ in the interior, although we do not know its dependence on $`Ra`$ etc at present. If $`T_{\mathrm{top}}`$ and $`Q`$ at the bottom are fixed as in our simulation, $`\mathrm{\Delta }T(t)`$ should also consist of fluctuations of the same origin. Because of the strong critical divergence of $`\gamma _s`$, we expect that the relative amplitude $`(\delta T)_{\mathrm{fl}}/\mathrm{\Delta }T`$ would increase as $`ϵ`$ is decreased with a fixed size of $`\mathrm{\Delta }T`$.

Fig.11 displays time sequences of $`\delta T(t)`$ and $`\mathrm{\Delta }T(t)`$ at fixed volume and pressure for periodic sidewalls with $`L_{}=4L`$, which demonstrates strong correlations between these two deviations at fixed volume. In case (a) (upper figure) we set $`ϵ=0.05`$ ($`\gamma _s=22.8)`$, $`\mathrm{\Delta }T=0.17`$mK, $`Pr=7.4`$, $`Ra=7.38\times 10^4`$, and $`Nu=4.06`$, while in case (b) (lower figure) we set $`ϵ=0.01`$ ($`\gamma _s=119)`$, $`\mathrm{\Delta }T=0.19`$mK, $`Pr=37.7`$, $`Ra=4.14\times 10^5`$, and $`Nu=6.04`$. The steady state values of $`\mathrm{\Delta }T`$ in the two cases are chosen to be only slightly different. At fixed volume, the fluctuations of $`\delta T(t)`$ and $`\mathrm{\Delta }T(t)`$ are strongly correlated, and are larger and slower for (b) than for (a) in steady states $`(t100)`$. This is because of the critical enhancement of the piston effect and the critical decrease of $`D`$ with decreasing $`ϵ`$. At fixed pressure, where the piston effect is absent, $`\mathrm{\Delta }T(t)`$ exhibits noises much smaller than those at fixed volume and $`\delta T(t)`$ smoothly changes in time. It is worth noting that this noise increase at fixed volume accompanied with an increase of $`Ra`$ is contrary to the usually measured noise behavior of the temperature. For non-critical fluids, if the temperature is measured at the center of a cell, its fluctuation amplitude divided by $`\mathrm{\Delta }T`$ is known to decrease with increasing $`Ra`$ as $`Ra^{\beta _n}`$. The exponent $`\beta _n`$ was about $`0.15`$ in a cell with $`A=1`$ .

### D Random Reversal of Macroscopic Flow

For a convection cell with $`A1`$, it is well-known that large-scale shear flow develops near the boundary of the cell for large enough $`Ra`$ . Moreover, it has also been observed that the global circulation changes its orientation over long time scales . For the case of $`A=1`$, $`ϵ=0.05`$, $`Q=40.7`$ $`\mu `$W$`/`$cm2, $`Ra=1.68\times 10^6(Ra^{\mathrm{corr}})`$, and $`Nu=5.97`$, we plot a numerical time sequence of a circulation $`\mathrm{\Gamma }(t)`$ in Fig.12. Here,

$$\begin{array}{ccc}\hfill \mathrm{\Gamma }(t)& =& _d^{Ld}𝑑x[v_x(x,Ld,t)v_x(x,d,t)]/L\hfill \\ & +& _d^{Ld}𝑑z[v_z(Ld,z,t)v_z(d,z,t)]/L,\hfill \end{array}$$
(3.22)
where the integration is along a square contour with distance $`d=0.05L`$ from the cell boundary. This quantity is positive for clockwise circulation and negative for counterclockwise circulation. In Fig.12, $`\mathrm{\Delta }T(t)`$ is also plotted, which exhibits particularly large fluctuations on the occasion of orientation changes. This is a natural result because large-scale reorganization of the flow pattern is needed for an orientation change. Fig.13 illustrates the velocity patterns at $`t=228,269`$, and $`311`$s in Fig.12. They closely resemble a picture of the measured velocity pattern in Ref..

## IV Scaling theory

Rayleigh numbers realized in the existing simulations are still moderate in the sense that the plumes do not have enough kinetic energies such that they do not generate fully developed turbulence in the interior. In this pre-asymptotic regime of steady states, we may understand the numerical and experimental data using a very simple zeroth-order theory. First, we set $`\mathrm{}=\mathrm{}_T=\mathrm{}_v`$ neglecting the possible small difference between $`\mathrm{}_T`$ and $`\mathrm{}_v`$ mentioned below (3.10). The plume sizes in the horizontal direction are also of order $`\mathrm{}`$. Second, in our simulation the plumes are ejected into the interior with a velocity $`v_{\mathrm{pl}}`$, for which the viscous drag and the buoyancy are balanced or

$$\eta \mathrm{}^2v_{\mathrm{pl}}g\alpha _p\mathrm{\Delta }T.$$
(4.1)
Thus $`v_{\mathrm{pl}}`$ is of the order of the terminal velocity $`v_{\mathrm{}}RaD\mathrm{}^2/L^3`$ in (3.6) with $`R\mathrm{}`$. In the interior we find that (i) gravity-induced acceleration of the plumes is suppressed by the viscous drag, (ii) $`(v_x^{})^2+(v_z^{})^2`$ is nearly independent of $`z`$ as stated below (3.10), and (iii) the last two terms on the left hand side of (2.17) are numerically of the same order. For example, the ratio of the average of $`(\alpha _pg\delta T)^2`$ in the $`x`$ direction to that of $`|(\eta /\rho )^2𝒗|^2`$ is about 4 at $`z\mathrm{}`$ and is fluctuating around 1 in the interior for the largest $`Q`$ in Table 1. These support $`v_{\mathrm{pl}}v_{\mathrm{}}`$ in the interior. Third, to the sum rule (2.36) for the velocity gradients, the contribution from the boundary layers is of order $`v_{\mathrm{pl}}^2/\mathrm{}L`$, while that from the interior is of order $`\varphi _{\mathrm{pl}}v_{\mathrm{pl}}^2/\mathrm{}^2Dv_{\mathrm{pl}}/\mathrm{}^3`$ from (3.20). If use is made of (4.1) and the sum rule (2.36), these boundary-layer and bulk contributions become both of order $`RaNu(D/L^2)^2`$, which has also been confirmed numerically. Thus,

$$v_{\mathrm{pl}}Ra^{1/2}D/L,$$
(4.2)
$$NuL/\mathrm{}1/\varphi _{\mathrm{pl}}Ra^{1/4}.$$
(4.3)
These quantities are independent of $`Pr`$. In particular, the independence of $`Nu`$ on $`Pr`$ is consistent with the experiments . Note that $`v_{\mathrm{pl}}(v_{\mathrm{}})`$ is smaller than $`v_\mathrm{g}`$ in (3.4) by $`Pr^{1/2}`$.

Our height-dependent Reynolds number at $`z=\mathrm{}`$ in (3.13) becomes

$$\widehat{R}e(\mathrm{})Ra^{1/4}/Pr.$$
(4.4)
The usual large-scale Reynolds number is given by $`Rev_{\mathrm{pl}}L\rho /\eta Ra^{1/2}/Pr`$. As $`\widehat{R}e(\mathrm{})`$ exceeds a crossover value $`Re^{}`$, plumes will induce turbulence in the interior. Our simple scaling theory is valid for $`Ra(Re^{}Pr)^4`$. In our simulation we have $`\widehat{R}e(\mathrm{})0.38Ra^{1/4}`$, so that if we set $`Re^{}10^3`$ (regarding plumes as jets ), the upper bound is crudely estimated as $`5\times 10^{13}`$. The transition from the scaling (4.3) to the asymptotic scaling occurs over a very wide range of $`Ra`$. Similarly, Grossmann and Lohse  considered a transition of a laminar boundary-layer flow to a turbulent boundary layer when the local Reynolds number on the scale of $`\mathrm{}_v`$ at $`z\mathrm{}_v`$ exceeds a value of order $`420`$. Then $`Nu`$ was claimed to be better expressed by

$$NuRa^{1/4}(1+C_1Ra^b)$$
(4.5)
than the single power-law form, where $`C_1`$ and $`b`$ are small coefficient and exponent, respectively, dependent on $`Ra`$ and $`Pr`$ under investigation (both being of order $`0.1`$). This proposed form of $`Nu`$ was later claimed to be in good agreement with data .

Here it would be informative to add more supplementary explanations of the previous scaling theories. (i) Shiraiman and Siggia  assumed fully developed turbulence in the interior. Then the maximum of the turbulent velocity gradient is of order $`S_\mathrm{d}=\eta k_d^2/\rho =(v_{\mathrm{pl}}^3\rho /L\eta )^{1/2}`$ at the smallest eddy size $`k_d^1`$ ($`(\eta /\rho v_{\mathrm{pl}})^{3/4}L^{1/4})`$ if the Kolmogorov cascade is assumed with the energy dissipation rate $`v_{\mathrm{pl}}^3/L`$  (the sparseness of ejected plumes being neglected). If the left hand side of the sum rule (2.36) is estimated as $`S_d^2`$, it follows the relation,

$$v_{\mathrm{pl}}(PrNuRa)^{1/3}D/L.$$
(4.6)
Furthermore, they assumed the linear horizontal velocity profile $`v_x(\sigma _0/\eta )z`$ in the region $`z\mathrm{}_TL/Nu`$, where $`\sigma _0`$ is given by (3.17). From the thermal diffusion equation $`v_x\delta T/x=D^2\delta T`$ (the time-dependent fluctuations being neglected), they obtained the scaling,

$$\mathrm{}_T^3\rho v_{\mathrm{pl}}^2/\eta DL,$$
(4.7)
by setting $`/xL^1`$ and $`^2(/z)^2\mathrm{}_T^2`$ for a cell with $`A1`$. From (4.6) and (4.7) they found

$$NuPr^{1/7}Ra^{2/7}.$$
(4.8)
However, as discussed below (3.14), our simulation suggests that the velocity deviates significantly from the linear profile in the boundary layers. (ii) Castaing et al.  assumed the balance (4.1) at the length $`\mathrm{}_T`$,

$$v_{\mathrm{pl}}\mathrm{}_T^2g\alpha _p\mathrm{\Delta }T/\eta RaNu^2D/L.$$
(4.9)
They furthermore assumed that the typical temperature scale in the interior is $`(\delta T)_cv_{\mathrm{pl}}^2/\alpha _pgL`$ and that the average heat current ($`Nu\lambda \mathrm{\Delta }T/L`$) is of order $`C_p(\delta T)_cv_{\mathrm{pl}}`$. From these relations $`(\delta T)_c`$ may be eliminated to give (4.6). If we combine (4.6) and (4.9), we are again led to (4.8). Therefore, to justify their arguments, the presence of fully developed turbulence in the interior seems to be required. (iii) Grossmann and Lohse  estimated the bulk and boundary-layer contributions to the sum rules for the temperature gradient and the velocity gradients, the incompressible version of (2.35) and (2.36). Their primary assumption is that the boundary layer thickness for the velocity is given by $`\mathrm{}_vL/Re^{1/2}`$ in terms of the large-scale Reynolds number $`Re`$ . Note that this assumption is not consistent with our zeroth-order scaling theory with respect to the $`Pr`$ dependence. In particular, in the case where $`Pr>1`$ and the boundary-layer contributions are dominant both for the temperature and the velocity, they obtained $`NuPr^{1/12}Ra^{1/4}`$. In this case we also find $`\mathrm{}_v/\mathrm{}_TPr^{1/3}`$ from their theory. They predicted that this pre-turbulent scaling crossovers to the asymptotic turbulent scaling very slowly as in (4.5).

## V Concluding Remarks

We have presented a hydrodynamic model of compressible fluids properly taking into account the piston effect and the adiabatic temperature gradient effect. Though performed in two dimensions, our simulation has revealed some new effects unique in near-critical fluids, such as the overshoot behavior and the amplification of the overall temperature fluctuations as $`TT_c`$. It generally explains the experimental findings , but a discrepancy remains in the overshoot behavior at high heat flux $`Q`$ as discussed in Section 3. It is desirable to extend simulation to smaller $`ϵ`$ and higher $`Ra`$. Also more experiments on the overshoot and the temperature noises etc. are needed to resolve the discrepancy and to confirm the new predictions. As by-products, we have numerically examined steady state properties not treated in the previous simulations, such as the logarithmic velocity profile and the random reversal of macroscopic shear flow. They are universal aspects present both in compressible and incompressible fluids.

We have assumed that the fluid is in the supercritical region not very close to the critical point such that the conditions (2.1) and (2.2) are satisfied. However, if $`\mathrm{\Delta }T`$ exceeds $`TT_\mathrm{c}`$ or if $`T_{\mathrm{top}}`$ is below $`T_\mathrm{c}`$, we encounter a variety of new effects such as boiling and wetting under heat flow and gravity . We believe that such problems should provide us a new challenging field in which nonlinear dynamics and phase transition dynamics are coupled. These problems are beyond the scope of this paper.

We thank H. Meyer for valuable suggestions and comments. Thanks are also due to P. Tong for informative correspondence. This work is supported by Japan Space Forum grant H13-264.

Fig.1. $`\mathrm{\Delta }T(t)`$ vs time (solid line) calculated from (2.8) and (2.17) for (a) $`Q=0.965`$ $`\mu `$W$`/`$cm2 and (b) $`Q=122.2`$ $`\mu `$W$`/`$cm2. The temperature profiles for the points ($`\mathrm{}`$) on the curve in (b) are given in Fig.2. The experimental data ($`+`$)  are shown in (b). The upper broken curves in (a) and (b) represent the theoretical result (3.3) obtained from integration of (2.8) with $`𝒗=\mathrm{𝟎}`$. The dotted curves represent the numerical ones in the fixed pressure condition without the piston effect.

Fig.2. Temperature profiles at (A), (B), (C), (D), (E), and (F) on the curve of $`Q=122.2`$ $`\mu `$W$`/`$cm2 in Fig.1b ($`\mathrm{}`$). The temperature (and velocity) deviations are more enhanced in the transient states (A)$``$(E) than in a steady state (F). The $`\delta T`$ at the bottom boundary $`z=0`$ is equal to $`\mathrm{\Delta }T(t)`$ in Fig.1b. The plumes tend to be connected between bottom and top because $`Pr=7.4`$.

Fig.3. Time evolution of $`\overline{\delta T}(z,t)`$ defined by (3.5) at the points (A), (C), (E), and (F) in Fig.1b for $`Q=122.2`$ $`\mu `$W$`/`$cm2.

Fig.4. Numerical results of $`Ra(Nu1)/(Ra^{\mathrm{corr}}Ra_\mathrm{c})`$ vs $`Ra^{\mathrm{corr}}/Ra_\mathrm{c}1`$ in steady states, obtained under the periodic boundary condition ($`+`$) and for $`A=3(\mathrm{}),2(),`$ and $`1(\times )`$. The first curve ($`+`$) is close to the experimental results for $`A=57`$  (solid line) and is well fitted to the scaling form (1.5) with $`a2/7`$ for $`Ra^{\mathrm{corr}}/Ra_\mathrm{c}10`$ . With decreasing the aspect ratio $`A`$, crossover to the scaling occurs at much larger $`Ra^{\mathrm{corr}}`$.

Fig.5. Height-dependent average temperature profiles $`\overline{\delta T}(z)`$ divided by $`\mathrm{\Delta }T`$ in steady states for the three $`Q`$ values in Table 1. The arrows represent the maxima of $`v_x^{}(z)`$ in Fig.6a.

Fig.6. Normalized height-dependent variances, $`v_x^{}(z)/v_\mathrm{g}`$ for the horizontal velocity in (a) and $`v_z^{}(z)/v_\mathrm{g}`$ for the vertical velocity in (b) in steady states for the three $`Q`$ values in Table 1.

Fig.7. Oveall Reynolds number $`\overline{R}e`$ defined by (3.11) as a function of $`Ra^{\mathrm{corr}}/Ra_\mathrm{c}1`$ in steady states for $`Q=122.2`$ $`\mu `$W$`/`$cm2.

Fig.8. Height-dependent Reynolds number $`\widehat{R}e(z)`$ defined by (3.9) in steady states for the three $`Q`$ values in Table 1.

Fig.9. (a) Height-dependent velocity variance $`v_x^{}(z)`$ defined by (3.9) on a semi-logarithmic scale in steady states for $`Q=122.2`$ $`\mu `$W$`/`$cm2. (b) $`v_x^{}(z)`$ (upperline) and velocity gradient variance $`zD_{xz}(z)`$ defined by (3.15) (lower line) on a logarithmic scale.

Fig.10. Snapshots of the normalized velocity variance $`v_z^{}(x,t)/v_\mathrm{g}`$ averaged in the $`z`$ direction defined by (3.17) for the three values of $`Q`$ in Table 1. The system is periodic with period $`4L`$ in the $`x`$ direction. The peak heights increase with increasing $`Q`$. For the largest $`Q`$ this quantity changes in time as the plumes move in the cell, while for the other $`Q`$ it is weakly dependent on or independent of time.

Fig.11. $`\delta T(t)`$ and $`\mathrm{\Delta }T(t)`$ at fixed volume (solid line) and at fixed pressure (broken line) for $`ϵ=0.05`$ (upper figure) and $`0.01`$ (lower figure). The noises of these quantities at fixed volume increase as the reduced temperature $`ϵ`$ is decreased.

Fig.12. Time evolution of the circulation $`\mathrm{\Gamma }(t)`$ defined by (3.22) (upper figure) and $`\mathrm{\Delta }T(t)`$ (lower figure) for $`Q=122.2`$ $`\mu `$W$`/`$cm2 in a cell with $`A=1`$. The orientation of the macroscopic flow changes on a time scale of 50 s. The sign of $`\mathrm{\Gamma }(t)`$ represents the orientation of the macroscopic circulation, while the fluctuations of $`\mathrm{\Delta }T(t)`$ become large when the orientation changes.

Fig.13. Velocity patterns at $`t=228,269`$, and $`311`$ s for the run in Fig.12. At $`t=228`$ s the orientation is counterclockwise, while at $`t=311`$ s it is clockwise. At $`t=269`$ s two large eddies with different orientations can be seen.