Rayleigh-Benard convection in a hard disk system
Pedro L. Garrido























[]garrido@onsager.ugr.es
Instituto Carlos I de Física Teórica y Computacional. Universidad de Granada. E-18071 Granada. Spain 



We do a generic study of the behavior of a hard disk system under the action of a thermal gradient in presence of an uniform gravity field. We observe the conduction-convection transition and measure the main system observables and fields as the thermal current, global pressure, velocity field, temperature field,...
We can highlight two of the main results of this overall work: (1) for large enough thermal gradients and a given  gravity, we show that the hydrodynamic fields (density, temperature and velocity) have a natural scaling form with the gradient. And (2) we show that local equilibrium holds if the mechanical pressure and the thermodynamic one are not equal, that is, the Stoke's assumption does not hold in this case. Moreover we observe that the best fit to the data is obtained when the bulk viscosity  depends on the mechanical pressure.

18-3e
Rayleigh-Benard convection in a hard disk system
    P.L. Garrido
    December 30, 2023
================================================





§ INTRODUCTION


Fluids are one of the subjects in physics that have been studied more intensively during the last centuries <cit.><cit.> <cit.>. Navier-Stokes (NS) equations have been the key ingredient in the understanding of observed fluid's phenomena in nature. Their analytic and numerical solutions are, nowadays, essential to many different scenarios, from the weather prediction up to the aircrafts design. In spite of the overall success obtained by its application to practical problems, there is a lack of a  rigorous connection between the underlying microscopic (typically hamiltonian) dynamics and the macroscopic NS-equations. Therefore there are still several important open issues connected to the NS main hypothesis and key ingredients used to build them. Let us mention a few of them as the local equilibrium hypothesis from which we assume that thermodynamics may be applied (up to some extend) at each macroscopic point of the fluid, the mathematical form of the transport coefficients (heat conduction and viscosities) as a function of the local fields, or the goodness of considering that the local heat conduction follows the Fourier's law or that the stress tensor is linear on the velocity gradients (newtonian fluid). The numerous efforts that have been done so far in the NS derivation show the enormous difficulties on such task. To have some flavour on it, let us mention  the nice review by Esposito et al. about the derivation of NS equations from the Boltzmann equation <cit.>. 
However, there is another way to study the connection between microscopic and macroscopic descriptions: the computer simulation tool. 

We know that fluids are typically in nonequilibrium states due to the action of external agents as temperature gradients or driving fields and there appear long range correlations in space and time. Therefore it is expected that their behavior depends strongly on the boundary conditions we impose. That makes a computer simulation of a fluid a very complicated task. It is necessary a priori a very large number of particles to have a reasonable collective description of the fluid and to limit (at some extent) the effects of the details of the boundaries in the bulk system behavior. Also, the long range correlations makes very difficult to make a finite size study of the simulation by trying to extrapolate any macroscopic observable value from a sequence of measurements of it in systems with increasing number of particles.
Finally, the time correlation functions are also power like which implies the need of very long runs to get a large enough set of almost independent configurations to do any observable averaging. In practice we immediately are limited by the finite power of our nowadays computers and we can use a very limited number of particles in our simulations (very far from the Avogrado's number) from which  we can get  some useful data in a computation time that does not exceeds a human life span. Obviously, with such small amount of particles it is difficult to expect any reasonable connection between those simulations and the behavior described by the Navier Stokes equations. In any case, there were some attempts during the eighties with some pioneering computer simulations trying to see macroscopic hydrodynamic behavior on particle systems. Let us mention for instance the work of D.C. Rapaport and E. Clementi <cit.> in which they simulate a fluid flow past a cylindrical obstacle by using a set of particles with a short range pair interaction potential. 

We have recently discovered <cit.> that a bidimensional hard disks system under the action of a temperature gradient, presents a nice property that we call boundary decoupling when studying its stationary state properties. That is, the system rearranges itself in such a way that there is a bulk part of it in which the local hydrodynamic properties hold: the local equilibrium hypothesis and the Fourier's law. The finite size dependences are then restricted to the values of the average heat current crossing the system and the local constant pressure. In other words, for a given stationary state with measured values of the heat current and pressure, the temperature profile in the system bulk is the one corresponding to the solution of the Fourier's law and the local temperature and density are related by the hard-disk equilibrium equation of state. We have checked this property also for a hard disk system under a shear forcing at the boundaries <cit.>. We want to push further the possibility to get some insight on the hydrodynamic behavior of a system by doing computer simulations. The Bénard problem is, in our opinion, the natural next problem to be studied by a computer simulation of a hard disk particle system. 

The Bénard Problem consists in a fluid that is heated from below under the presence of gravity. It is observed that above some critical value of the temperature gradient the convection starts and the fluid move in rolls occupying the full system. This is the Rayleigh-Bénard instability that was observed experimentally by Bénard and afterwards there was interpreted theoretically by Rayleigh.
That is, this nonequilibrium problem includes a transition from conducting to convecting states and a nonzero hydrodynamic velocity field in the convecting stationary state (see for instance some recent books and reviews in <cit.>, <cit.> and <cit.>). 

There are several papers studying the Bénard problem by using simulations of hard disks. Let us mention, for instance, the work from M. Marechal et al. <cit.><cit.> and A. Puhl et al. <cit.>. Mareschal and co-workers observed the transition conducting-convecting for a system of hard disks  even for a very small number of particles (N=5040). However, the very large velocity fluctuations prevented, in that time, to get accurate measurements of local magnitudes. Nevertheless they managed to get some hydrodynamics velocity profiles that they could compare with a numerical resolution of NS equations getting a very reasonable accordance. This system was later studied also by D.C. Rapaport <cit.><cit.> showing the influence of the boundary conditions and the initial configurations on the stationary behavior of the hard disk system. Finally, let us also mention the work from D. Risso and P. Cordero in which the study the onset of the transition <cit.>.

In this paper we show a complete set of computer simulations results whose goal is two fold: first we want to get some insights on the stationary behavior of the system and its dependence on the external parameters (gravity and temperature gradient) by observing the heat current, internal energy, pressure,... and second, we would like to check some of the basic Navier Stokes assumptions like the local equilibrium hypothesis, the Fourier's law, the difference between the mechanical pressure versus the thermodynamic one or the efects of our system compressibility  in the transition <cit.>.



  This paper is divided in several parts. Let us make a summary of the main items and results obtain on each part.
  


  * I. The physical properties of a hard system: First we introduce a catalog of known properties of a hard disks system. We do not know the equilibrium equation of state but we express all the thermodynamic equilibrium quantities as a function of it. We also derive the hydrostatic formula in equilibrium. Finally we write down the  Navier Stokes equations,  the general structure of the transport coefficients and their particular analytical form obtained in the Enskog approximation of the Boltzmann equation for hard disks.


  * II. The model:  We describe the model, the  dynamics of particles and the boundary conditions used in the simulation. We also discuss the election of the parameter values for the simulation. We explain how we do the local measurements of the magnitudes in order to get the hydrodynamic fields. In a computer simulation like this where the ratio between noise and signal is expected to be not small is important to do a careful analysis of errors. We explain the way we have done such analysis. Let us remark that along the paper we use an strict 3σ error interval criterion. Finally we study the relation between the averages of the magnitudes on a local region and its continuous value. That is important to have some clean scheme to define the spatial derivatives of any field.


  * III. Computer simulation results of equilibrium states: We check the program and the way we measure the global and local observables by doing a set of simulations at equilibrium. That is, when the temperatures on both plates are equal but maintaining the gravity field.  We show that the local densities and temperatures follow accurately the theoretical results. Moreover, the local pressure measured by using the virial expression is found to give the correct results. Finally we observe that it is meaningless  to measure local fluctuations of the magnitudes because they are very sensitive to the size of the local cells that in our case is certainly very small.


  * IV. Computer simulation results of nonequilibrium stationary states: global magnitudes. In this part we show the behavior of the global observables as a function of the temperature gradient Δ T and the intensity of the gravity field g. The kinetic energy per particle increases with g in contrast to what happens in equilibrium and for large Δ T it grows like ≃Δ T/3 for all g's. The fluctuations of the kinetic energy increase with Δ T  for all g's and for a fixed Δ T it increases with g. The interesting thing is that the fluctuations divided by the average kinetic energy to the square presents a non-monotonous behavior for a given g value: starting from the equilibrium, it decreases with Δ T, reaches a minimum and then it increases. As we show later in the paper, this minimum is consistent with the transition point between conducting and convecting regimes. Moreover, the critical value of the gradient  at which the transition occurs increases with g. That behavior contrast with the one derived in the Boussinesq approximation of the Navier Stokes equations in which the critical value decreases with g. We explain there that the reason is that our system is compressible and then one should correct the Boussinesq assumptions accordingly.

We also measure the hydrodynamic kinetic energy, that is, the average particle velocity on a cell, squared and summed over all cells. Obviously this magnitude should be zero in the conducting regime and non-zero in the convecting one. We obtain the transition point from the data that is consistent with other methods used in the paper. We also study the hydrodynamic kinetic energy associated with each velocity component. We observe a kind of equipartition independently of the g's for large values of Δ T: e_u,2≃ 3.4 e_u,1 where e_u,i is the hydrodynamic kinetic energy associated to the ith component of the velocity.

The Pressure is computed at the top and bottom boundaries as the moment variation of the particles when colliding with it. We observe that the global barometric formula holds always, that is the difference between the pressures at bottom and top of the system is proprotional to the system weight. This confirms the correctness of the computer simulation. Nevertheless the pressure at each boundary has a non trivial dependece on Δ T and g which we analyze.

Finally we study the energy current that cross each thermal bath boundary. We first check that the current coming into the system is exactly the same to the outgoing one. For large values of Δ T the current goes like Δ T^3/2 for all g's. In general we see that  the gravity hinders the heat current for small values of Δ T while it favours it for large values of it where it correspond to the fully developed convective regime. We define a second critical thermal gradient as the one in which the thermal current for a given g value equals the one when g=0 (pure conducting case). That is, beyond such value, convection becomes “more efficient” that conduction in transporting heat through the system for a given temperature gradient. This second critical gradient is necessary to explain some observed system behavior later in the paper.
 

  * V. Computer simulation results of nonequilibrium stationary states: spatial structures. First we study the distribution of the current line lengths for a set of stationary velocity fields. For conducting states we reasonably fit a gamma distribution to the data. In the convecting regime there is a distribution that depends on the system geometry and we analyze its form. We also analyze  the fraction of velocity vectors whose modulus is larger than their standard deviation, pr. That is, the amount of “real” (or nonfluctuating) hydrodynamic velocity vector. We observe that this magnitude is directly correlated with the average length of the current lines in each case. Moreover, that relation is independent on the value of g. We see that for all g's we obtain pr≃ 0.68 when we are at the second critical gradient. This value coincides with the critical fraction of the covered surface at the percolation threshold of a system of overlapping disks of radius r. That is, the overall picture maybe that when convection begins, the current lines are small. Their length grow with the value of the thermal gradient up to a critical value in which the current lines extends and connect the system boundaries. We conclude that this phenomena is related with a percolation like transition that occurs precisely at the second critical temperature gradient. 

We study the spatial structure of the hydrodynamic velocity components. For large enough temperature gradients it seems that there is  well define limiting velocity profile. Moreover we show that each hydrodynamic velocity component has an universal scaled  field (universal in the sense that, once the field is scaled, it does not depend on Δ T) for large enough values of Δ T. The scaled field, ũ_1,2(x,y) is obtained by

    ũ_1,2=u_1,2(x,y)/σ(u_1,2)

where u_1,2(x,y) is the components 1,2 of the measured hydrodynamic field and σ(u_1,2) are their stardard deviation (the spatial average of the hydrodynamic velocity field is zero). In the process we observe the nice property: σ(u_2)=2σ(u_1) for all fields with any Δ T and g values. We  study  the Inertial Tensor and the fourth moment of scaled fields with different values of Δ T to look for any systematic dependence on the temperature gradient. Moreover, we compute the average mutual euclidean distance between scaled configurations to  confirm that, for large enough values of Δ T such distance tends to zero. The overall analysis is lengthly because the initial data fields have non-negligible fluctuating spatial behavior and intrinsic error bars that should be taken into account to finally obtain a consistent set of results.

We also study the Temperature field. The average y-profile, τ(y) is non-linear and it is convex for g=0 and 5 and it has infection points for g=10 and 15. We see thermal gap in the boundaries, characterize its behavior  and compute the effective thermal gradient that is ≃ 0.82 Δ T for all g-values. Finally we show also the existence of a universal scaled Temperature field that is defined by

    T̃(x,y)=T(x,y)-τ(y)/σ(T)

where σ(t) is the standard deviation with respect τ(y). 
In the same spirit we also analyze  the density and pressure fields. In both cases we observe the existence of universal fields.


  * VI. Connecting with hydrodynamics. In this section we try to check some of the Navier Stokes ansatz and properties at the stationary state. Dynamic properties are beyond the scope of this article. The first thing we check is the continuity equation that express the mass conservation at the hydrodynamic level:
∇(ρ u)=0 from the measured local density and velocity fields. This give us some confidence on the use of the field spatial derivatives in spite that the relative error after derivatives increases up to an average of 10%, large but enough to check some equation consistencies. Second we see if local equilibrium holds at each cell. We observe that it is the case for systems at the conducting regime. However we see a small and noisy systematic deviation from local equilibrium for systems at the convecting regime. We argue that it is so because of the difference between the mechanical pressure (the one we measure) and the thermodyncamic one. Their difference it is known to be proportional to the divergence of the hydrodynamic velocity. The hydrodynamic velocity divergence is small and the data is noisy but nevertheless we manage to discard the Stoke's assumption (∇ u=0) and we find that the best fit occurs when we assume that the bulk viscosity is itself dependent on the mechanical pressure. This is a weak-like result but indicates a way to go deep inside of the properties of the NS equations.
Finally we check stationary NS equations for the non-convective states obtaining a reasonable good coincidence with the computer simulation results.



The overall work has been an effort to extract a coherent and systematic description of the Rayleigh-Benard problem from a microscopic computer simulation of a hard disk system. Of course there are many interesting items and observables we could study. Many of them we honestly tried to analyze (for instance local fluctuations) but we finally applied the strict criterium that anything whose error analysis obscure the average structure should be discarded. Maybe, one may focus on a particular item to design a simulation that gets much better statistics necessary to obtain with good confidence its behavior. We feel that it can be done and it is possible to study hydrodynamic behaviors from microscopic models.



 



§ THE PHYSICAL PROPERTIES OF A HARD DISK SYSTEM


Hard disks is one of the paradigmatic models in many body physics. Having a very simple dynamics it contains many nontrivial  statistical properties that come from collective effects of its individuals. For instance, it has a  phase transition <cit.> or its time correlation functions have the long time tails <cit.>. Moreover its exact Equation of State (EOS) is not yet
derived in the context of the equilibrium statistical physics <cit.>.  Let us state some basic properties that are known and that we are going to use during the computer simulation data analysis.



  * Equilibrium properties:
The canonical partition function for a system with N disks with mass m and radius r enclosed  in a box of sides L_x an L_y at a temperature T is given by
 
    Z(N,S,T)=1/h^2NN!(2π mk_B T )^N Z_c(N,S)
 
 where S is the box surface (S=L_xL_y) and Z_c is the configurational partition function that, in the hard disk case, it doesn't depend on the temperature. Z_c is not known but it can be connected with the equation of state that has been extensively studied numerically and very good analytical expressions have been derived from the extensive computer simulations <cit.>. The equation of state for a hard disk system has the general form:
 
    PS=Nk_BT H(ρ)     ,    ρ=Nπ r^2/S

A well known simple form for H was proposed by Henderson <cit.> that fits  very well the data for de areal density interval ρ∈[0,0.5]:
 
    H(ρ)=1+ρ^2/8/(1-ρ)^2
 
The equation of state is connected with Z_c by the expression:

    P=k_BT∂ Z_c/∂ S  ⇒   NH(ρ)=-ρ∂ Z_c/∂ρ

If we assume that  Z_c≃exp[Nf̅(ρ)] for N large enough, we obtain:

    f̅(ρ)=f̅(ρ^*)-∫_ρ^*^ρdρ̃  H(ρ̃)/ρ̃

where ρ^* is any arbitrary value. We know that in the limit ρ→ 0 the hard disk system behaves like an ideal gas. In such limit f(0)=ln S, substituting  ρ^*=0  above we get

    f̅(ρ)=1/Nln Z_c=ln(π r^2N/ρ)-∫_0^ρdρ̃  H(ρ̃)-1/ρ̃

and the hard disk partition function can be written for N large:

    Z≃ e^Nf(ρ)    ,   f(ρ)=lnT/ρ-∫_0^ρdρ̃  H(ρ̃)-1/ρ̃+C

where C is the constant

    C=ln(2mk_Beπ^2 r^2/h^2)

 From it we can compute all thermodynamic magnitudes, for example, the entropy is given by
 
    s/k_B=f(ρ)+1

the specific heats at constant pressure and volume are respectively 

    c_P/k_B=1+H(ρ)^2 /H(ρ)+ρdH(ρ)/dρ   ,    c_V=k_B

and the coefficients α and κ_T are 

    α=1/TH(ρ)/H(ρ)+ρdH(ρ)/dρ   ,   κ_T=π r^2/k_B Tρ1/H(ρ)+ρdH(ρ)/dρ



  * Hard disks under the action of an external field at equilibrium: the hydrostatic formula

Let us assume now that over the hard disk system, there is an external field acting on the particles. That is, we should add to the Hamiltonian a term of the form:

    H_ext=∑_i=1^NU(r⃗_i)
 
Then, to obtain the equilibrium properties it is convenient to use the grand canonical ensemble. The computation simplifies if we consider a very slow varying external field (we follow here the arguments in Martin-Lof's book <cit.>). That is, U(r⃗)=Ũ(r⃗/L) being Ũ a smooth function in R^2. We can define  a box with a volume that grows with L^2 at a given position x⃗. When L→∞ the particles of the bulk of the box only see an external constant field of value Ũ(x⃗) and then H_ext=NŨ(x⃗). The local thermodynamic properties of the system can be  computed using the grand canonical ensemble just substituting the chemical potential μ→μ-Ũ(x⃗). For instance, the equation of state is given in a parametric form:

    P   =   k_BTa(T,μ-Ũ(x⃗))
    ρ̃   =   k_BTd/dμa(T,μ-Ũ(x⃗))

where ρ̃ is the particle density, a=a(T,μ)=lim_S→∞S^-1log𝒵 and 𝒵 is the gran canonical partition function. We can eliminate the parameter μ by applying the gradient to the pressure and using the definition of ρ̃. Thus, we obtain the hydrostatic formula:

    ∇ P=-ρ̃(x⃗)∇Ũ

Notice that the full argument is general and it is not restricted to a hard disk system.


  * Navier-Stokes Equations for hard disk systems:

Navier Stokes equations describe the dynamics of a macroscopic fluid <cit.> and they are constructed assuming: (1) the local conservation of density, linear momentum and  energy, (2) the local equilibrium property and (3) some local constitutive relations (Fourier's law, Newtonian fluid). The macroscopic evolving fields are: the local mass density, ρ̃(x,t), the local velocity components, u_i(x,t)  ,  i=1,…, d, and the local temperature, T(x,t). Their evolution is determined by:

    ∂_tρ̃+∑_i∂_i(ρ̃u_i)   =   0
    ∂_t u_i+∑_j u_j∂_j u_i   =   1/ρ̃∑_j ∂_jτ_ij+g_i   (i=1,…,d)
    ρ̃c̃_v(∂_t T+∑_i u_i∂_i T )   =   -T∂ P/∂ T|_ρ̃∑_i∂_i u_i+∑_ijτ'_ij∂_i u_j+∑_ij∂_i(κ_ij∂_j T )

 Assuming that we study the so called newtonian fluids then: 

    τ_ij   =   -Pδ_ij+τ'_ij
    τ'_ij   =   η(∂_i u_j+∂_j u_i )+η' δ_ij∑_k ∂_k u_k

Finally, the local equilibrium property also implies that  the equilibrium equation of state applies to the local fields: P=P(ρ̃,T).  Observe that  we also need to know the equilibrium mass specific heat at constant volume c̃_v=c̃_v(ρ̃,T) and  the transport coefficient functions: the thermal conductivity κ=κ(ρ̃,T), the shear viscosity η=η(ρ̃,T) and the second viscosity η'=η'(ρ̃,T) . All these equations plus the value of the external field g_i and the boundary conditions for the fields give us  a complete description of the evolution of a fluid at the macroscopic level. 

For hard disks system we know the EOS and then 

    ∂ P/∂ T|_ρ̃=P/T

and c_v=k_B. Notice that  ρ=Nπ r^2/L_xL_y  and it is related to the mass density by: ρ̃=mρ/π r^2 where m is the mass of one disk and r its radius.

Moreover one can show that the transport coefficients have also a simple dependence on the temperature:  

    η   =   √(k_Bm)/r√(T)E(ρ) 
    η'   =   √(k_Bm)/r√(T)E'(ρ)
    κ   =   √(k_B^3)/r√(m)√(T)K(ρ)


We are going to be interested on the stationary solutions when the external field is g_j=-gδ_j,2. The resulting equations after we substitute these formulas into eqs. (<ref>) are: 

    ∂_1(ρ u_1)+∂_2(ρ u_2)   =   0 
    ρ(u_1∂_1 u_1 +u_2∂_2 u_1 )   =   -∂_1Q+π r∂_1[√(T)(2E∂_1u_1+E'(∂_1u_1+∂_2u_2) )]
       +   π r∂_2[√(T)E(∂_1u_2+∂_2u_1)]
    ρ(u_1∂_1 u_2 +u_2∂_2 u_2 )   =   -∂_2Q+π r∂_2[√(T)(2E∂_2u_2+E'(∂_1u_1+∂_2u_2) )]
       +   π r∂_1[√(T)E(∂_1u_2+∂_2u_1)]-ρ g
    ρ(u_1∂_1 T +u_2∂_2 T )   =   -Q(∂_1u_1+∂_2u_2)+2π r√(T)E((∂_1 u_1)^2+(∂_2 u_2)^2)+π r√(T)E'(∂_1u_1+∂_2u_2)^2
       +   π r√(T)E(∂_1u_2+∂_2u_1)^2+π r[∂_1(√(T)K∂_1T)+∂_2(√(T)K∂_2T)]

where we have defined the reduced pressure, Q=π r^2 P, and  let us remind that E, E' and K are functions of ρ. We have also done, after the computation, the following mapping:

    k_BT→ T   ,   mg→ g   ,  √(m)u_1,2→ u_1,2

With this mapping the equation of state is:

    Q=TF(ρ)=Tρ H(ρ)

Observe that we use in all the equation the areal density ρ instead the mass density ρ̃.

The NS equations for hard disks simplify when looking for non-convective  stationary solutions:

    ρ=ρ^NC(y)  ,    u_x=u_y=0  ,   T=T^NC(y)
 
that are solutions of:

    dQ^NC/dy   =   -gρ^NC
    √(T)K(ρ^NC)dT^NC/dy   =   -J
    
    Q^NC   =   T^NCF(ρ^NC)

with the boundary conditions:

    T(0)=T_0  ,   T(1)=T_1=1  ,  ρ̅=∫_0^1dy ρ^NC(y)

where T_0, g and ρ̅ are the control parameters. Observe that J is a constant that it is be obtained as a function of the external parameters.


  * Transport coefficients (in the Enskog approximation):
Gass  computed the transport coefficients for the hard disk system in the Enskog approximation <cit.>. The shear viscosity η, the bulk viscosity ξ and the thermal conductivity κ are:

    η   =   √(k_Bm)/r√(T)E(ρ)    ,    E(ρ)=0.255/√(π)χ(ρ)[1+2ρχ(ρ)+3.4197ρ^2 χ(ρ)^2]
    ξ   =   √(k_Bm)/r√(T)E_B(ρ)    ,    E_B(ρ)=1.2734 ρ^2χ(ρ)
    κ   =   √(k_B^3)/r√(m)√(T)K(ρ)    ,    K(ρ)=1.029/χ(ρ)[1+3ρχ(ρ)+3.4874ρ^2 χ(ρ)^2]

where

    χ(ρ)=H(ρ)-1/2ρ

and the second viscosity η'=ξ-η.








§ THE MODEL


Our system is composed of a set of N hard disk particles that move in a square box of side L=1. Each disk have unit mass (m=1) and its radius (r) is chosen such that the areal density is ρ̅=Nπ r^2/L^2, that is r=(ρ̅L^2/(Nπ))^1/2. The dynamics of a given disk has two parts: (1) free flight and (2) collisions. 

(1) Free fligth: Between collisions each disk i is under the unique action of a constant external field a⃗=(0,-g). Thus, its equations of motion are given by 
    dv⃗_i/dt=a⃗
 whose trivial solutions are 

    r⃗_i(t)=r⃗_i(0)+v⃗_i(0)t+1/2a⃗ t^2  ,  v⃗_i(t)=v⃗_i(0)+a⃗ t
 
where r⃗_i(0) and v⃗_i(0) are the values of the disk position and velocity after a collision of the disk has taken place.

(2) Collisions: A disk have three types of collisions, with another disk, with the box vertical boundaries and with the box horizontal boundaries. (a) Two  disks collide when their distance equals 2r. When it happens, the linear moment and the kinetic energy is conserved and it is assumed that the velocity components perpendicular to the vector that joints the centers of the disks keep unaltered during the collision. That is, we do not consider the existence of disk internal rotation around its center. (b) When the disk collides with a vertical side of the box at x=0 or x=1 it gets perfectly reflected, that is v_1→ -v_1 and v_2→ v_2. (c) when the disk hits a horizontal boundary at y=0 or y=1 the disks changes its second component of its velocity by getting a random value chosen from a Maxwellian distribution with temperature T_0 or T_1 when hitting at y=0 or y=1 respectively. That is v_1→ v_1 and v_2→ v_2=s_0,1(-2T_0,1log(1-ξ)) being ξ an uniform random variable between [0,1] and s_0=1 and s_1=-1.



Initially the disks are placed regularly on the box with a random initial velocity vector with a modulus that is the square root of the sum of the two thermal bath temperatures. We evolve the system during 5× 10^4 collisions per particle (≃ 5× 10^7 collisions) when we check that the system is in the stationary state. Then we begin to take measures of different observables each 100 collisions per particle. The simulation extends during 10^7 collisions per particle (≃ 10^10 total number of collisions) and we obtain  M=10^5 data for averaging. The observables have typical errors of 3σ/√(M)≃ 0.01 σ that is sufficient to analyze with some detail many of their interesting behavior.


Computer simulations of hard disks systems have many advantages. First, the intrinsic dynamics is very symple and it can be implemented in the computer with very high efficiency without any loss of precision (see for instance <cit.><cit.>). Moreover, its purely kinetic structure makes possible to fix one of the system external parameters (T_0,T_1,g) by just  a time rescaling. In other words: If we rescale time, t=α t', the particle velocities rescale by v=v'/α and then, choosing the temperature of the thermal baths being T'_0,1=α^2 T_0,1 and the gravity field g'=α^2 g, the dynamic evolution of the system with parameters (T_0,T_1,g) is indistinguishable from the one with parameters (T_0',T_1',g'). We can arbitrarily choose α=1/√(T)_1 in order to fix to one the temperature of the the thermal bath at top. Therefore we can do the simulations varying only two parameters  (T_0',T_1'=1,g'). In order to obtain the behavior of any observable for an arbitrary value of T_1 we just should undo the corresponding time rescale on the variables defining the magnitude. 

In order to choose the set of values for the parameters ρ̅ (the average areal density), T_0 and g we should pay attention to the phenomena we want to describe. In this case we would like to study the transient from non-convective regime to the convective one. Then it is reasonable to use the Rayleigh coefficient (Ra) that  is a non-dimensional magnitude that localizes the parameter regions where the fluid is in a convective/non-convective state. It is defined by

    Ra=α gΔ T L_y^3/νκ̅

where g is the intensity of the external field, Δ T=T_0-T_1, ν=η/ρ̃ is the kinematic viscosity, κ̅=κ/ρ̃C̃_P and ρ̃=mρ/π r^2 is the mass density and C_P=c_p/m is the specific heat capacity per unit mass.

For a hard disks case it is known that linearizing NS equations under the Boussinesq approximation there is a critical Rayleigh number Ra_c=27π^4/4≃ 657.51 for the stress free boundary condition case (see for instance S. Chandrasekhar <cit.>). Convection appears when the system has a Rayleigh number above the  critical value.

We have computed Ra by using the Henderson equation of state and the Enskog transport coefficients. We choose the number of particles to be small enough to have a very fast system evolution to get the largest quantity of data to average. One may think that N small would implies very large finite size effects. However,  we already shown in previous studies with hard disks <cit.> that the system have a boundary decoupling property in which the system re adapts its bulk configuration to behave as it was an infinite system with effective boundary conditions composed by the ones we define in the simulation plus a small region with disks around them. For the computer simulation in this paper we chose N=957.
 
 In figure <ref> we show Ra computed using the Enskog transpot coefficients and the Henderson equation of state. In figure <ref> left we see the behavior of Ra for T_0=5 for different values of g=5, 10 and 15. We see that the maximum value of Ra apears for low densities and it is above the critical Rayleigh number. Then, we have chosen that the average system density to be ρ̅=0.2.  Figure <ref> right we plot Ra as a function of T_0 for fix ρ̅=0.2 to have some idea of the critical value of T_0 that separates non-convective to convective regimes. We get the following critical temperatures:
 T_0^c(g=5)=1.6205,  T_0^c(g=10)=1.2233 and  T_0^c(g=15)=1.1376 . That is, when we increase the external field, the convective regime appears for lower temperatures.
 


From now on we will call to the simulated parameters (T_0, T_1=1,g) whose chosen values are: g=0, 5, 10 and 15, and  T_0=1,1.2,1.4,1.6,1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3, 4, …, 19, 20. That is, we have done 112 computer simulations. 

Before begining the analysis of the data we have obtained let us make some comments about notation, how the measurements have been done and the error bars criteria we have used.




  * Notation: In general a vector u have two components u=(u_1,u_2). In order to do local measurements, we divide the unit box into 30× 30 (virtual) square cells of side 1/30. A given cell has the position (n,l) with n,l=1,…,30. We express hydrodynamic magnitudes  at a given macroscopic point (x,y) with x,y∈[0,1] that correspond to the values measured of some microscopic observable at the cell (n,l) such that   x=(n-1)/30+1/60 and y=(l-1)/30+1/60. That is, we use the center of each cell as the macroscopic position of the hydrodynamic fields. Finally, a macroscopic hydrodynamic field is written as u=(u_1(x,y),u_2(x,y))


  * Local Measurements:
The local measurements have been done in the following way. Let a(r,v) be a magnitude that depend on a set of particle positions and velocities r_i and v_i respectively. Say for instance the kinetic energy of a particle: a(r_i,v_i)=v_i^2/2=(v_i,1^2+v_i,2^2)/2 or the local potential energy of a particle: a(r_i,v_i)=g r_i,2 where r_i=(r_i,1,r_i,2). The extensive value of a(r,v) on the cell (n,l)  (n,l=1,…, 30) at time t is given by: 

    A(n,l;t)=∑_i:r_i(t)∈ B(n,l)a(r_i(t),v_i(t))

where B(n,l) is the spatial domain corresponding to the cell (n,l). The total number of particles involved in such sum is given by:

    N(n,l;t)=∑_i:r_i(t)∈ B(n,l) 1

We do M measurement of A(n,l;t) during the system evolution at the stationary state. Then the average value per particle of A is given by:

    a(n,l)=∑_t=1^M A(n,l;t)/∑_t=1^M N(n,l;t)
 
We use this method because it has a better convergence to the limiting value when M→∞ and its fluctuations are smaller than in case we averaged A(n,l;t)/N(n,l;t) over time.


  * Analysis of errors:
Let {A(t)}_t=1^M the set of M measurements of a given observable (local or global one). We assume that the data sequence is time decorrelated and then the law of large numbers apply. That is, the averaged value of A behaves for large M values as:

    A(M)≃ A(∞)+σ(A;M)ξ

where

    A(M)=1/M∑_t=1^M A(t)    ,  σ(A;M)=1/M√(∑_t=1^M(A_t-A(M))^2)

and ξ is a Gaussian random variable with zero mean and variance one. In this paper the error in getting the M→∞ value is given by

    A(∞)=A(M)± 3σ(A;M)

that is, we are assuming that the 99.7% of the data is in the error interval we show in the analysis.

Sometimes, once we have computed the average value of the observable A and its error, we need to study a function of it, say f(A). Let us assume that the error is small and It can be considered as a perturbation of A(∞), then we can do a Taylor expansion:

    f(A(M))=f(A(∞)+σ(A;M)ξ)=f(A(∞))+f'(A(∞))σ(A;M)ξ+1/2f”(A(∞))σ(A;M)^2ξ^2+…

If we average this expression over the random values ξ we obtain that ⟨ f(A(M))⟩=f(A(∞))+1/2f”(A(∞))σ(A;M)^2. There is a small but systematic shift on the observed value f(A(M)) with respect the desired f(A(∞)). Let us define the observable:

    B=f(A(M))-1/2f”(A(∞))σ(A;M)^2

Obviously the average value of B is f(A(∞))  (by construction) up to order σ^4 and the variance of B is:

    σ(B)^2   =   ⟨(B-f(A(∞)))^2⟩=⟨(f'(A(∞))σ(A;M)ξ+1/2f”(A(∞))σ(A;M)^2(ξ^2-1) )^2⟩
       =   σ(A;M)^2[f'(A(∞))^2+1/2f”(A(∞))^2σ(A;M)^2]

The error in computing B is then  B=f(A(∞))± 3σ(B) and we conclude that

    f(A(∞))=f(A(M))-1/2f”(A(∞))σ(A;M)^2± 3σ(B)

Typically the systematic shift is smaller than the error interval and therefore not relevant for the analysis. However, sometimes we obtain small but positive data values for an observable that one knows that it should be zero in average and we should introduce this correction in order to do a correct data analysis. 

Let us consider another typical case: when we compute an observable as an spatial average of local observables with their own error. That is, let A(x,y) a local observable with errors ϵ(x,y). Let

    B=1/N_C∑_(x,y)A(x,y)

If we assume that the local random variables ξ(x,y)  associated to the errors ϵ(x,y) are spatially independent, then we can apply the Lyapunov Central Limit Theorem (see Appendix I) and  the error of B is given by:

    ϵ(B)=1/N_C[ ∑_(x,y)ϵ(x,y)^2]^1/2

Notice that the assumption of spatial independence is not so obvious because the fluid has typically long range correlations that should affect to the data noise behavior. 

However, let us assume that the data error is totally correlated in the sense that the local random variable can be written as ξ(x,y)=ϵ(x,y)ξ  where ξ is now a Gaussian noise with zero mean and variance one. In this case the error of B is given by:

    ϵ(B)=1/N_C∑_(x,y)ϵ(x,y)

Observe that the totally correlated case has an error much larger that the uncorrelated one. It is not clear a priori what scheme to choose to define errors in this case without computing the data correlation function that it is very costly from a computational point of view. In any case the true error on B would have a value in between of the two limiting cases.


  * Differences between Local values and cell averages:

To obtain local values of an observable we have build a set of virtual cells where we have
measured the desired magnitudes. That is an arbitrary choice and we could define other partition and associate a value of the magnitude on each of its elements. Therefore we can assume that at the stationary state, for each microscopic observable, there exists a function defined in ℛ^2 such that its averages over any partition is the value measured on the cell. In other words:
 
    F̃(m,n)=1/Δ^2∫_x̃(m,n)-Δ/2^x̃(m,n)+Δ/2dx ∫_ỹ(m,n)-Δ/2^ỹ(m,n)+Δ/2dy  F(x,y)

  where F̃(m,n) is the average value of the observable that we obtain in the computer simulation and F(x,y) its, assumed existing, underlying continuum density. We consider that 
 x̃(m,n)=(m-1/2)Δ, ỹ(m,n)=(n-1/2)Δ where Δ=1/N_C and m,n=1,…,N_C (N_C=30 in our computer simulations). The question that arises is what is the relation between this local averages and the value of  F and its derivatives at the point (x̃(m,n),ỹ(m,n)). In order to find such relations  let us write the value of F̃ for neighboring cells:
 
    F̃(m+s,n+t)=1/Δ^2∫_x̃(m,n)+sΔ-Δ/2^x̃(m,n)+sΔ+Δ/2dx ∫_ỹ(m,n)+tΔ-Δ/2^ỹ(m,n)+tΔ+Δ/2dy  F(x,y)
 
 We can expand the last expression assuming that F is analytic and Δ is small enough, then:
 
    F̃(m+s,n+t)   =   ∑_n'=0^∞Δ^2n'∑_k'=0^n' F_(2k',2n'-2k')(m,n) a(s,k)a(t,n'-k')
       +   ∑_n'=0^∞Δ^2n'+2∑_k'=0^n' F_(2k'+1,2n'-2k'+1)(m,n) b(s,k)b(t,n'-k')
       +   ∑_n'=0^∞Δ^2n'+1∑_k'=0^n' F_(2k'+1,2n'-2k')(m,n) b(s,k)a(t,n'-k')
       +   ∑_n'=0^∞Δ^2n'+1∑_k'=0^n' F_(2k',2n'-2k'+1)(m,n) a(s,k)b(t,n'-k')

 where F_(n',l')(m,n)=∂_x^n'∂_y^l' F(x,y)|_x̃(m,n),ỹ(m,n)) and
 
    a(s,l)   =   ∑_k=0^l1/(2k)!(2l-2k+1)!s^2k/2^2l-2k
    
      b(s,l)   =   ∑_k=0^l1/(2k+1)!(2l-2k+1)!s^2k+1/2^2l-2k

 Observe that a(-s,l)=a(s,l) and b(-s,l)=-b(s,l). These symmetry properties allows us to break eq. (<ref>) in a set of four equations:
 
    G_1(m,n;s,t)   ≡   1/4[ F̃(m+s,n+t)+ F̃(m+s,n-t)+ F̃(m-s,n+t)+ F̃(m-s,n-t)]
       =   ∑_n'=0^∞Δ^2n'∑_k'=0^n' F_(2k',2n'-2k')(m,n) a(s,k)a(t,n'-k')
    
       G_2(m,n;s,t)   ≡   1/4[ F̃(m+s,n+t)- F̃(m+s,n-t)-F̃(m-s,n+t)+ F̃(m-s,n-t)]
       =   ∑_n'=0^∞Δ^2n'+2∑_k'=0^n' F_(2k'+1,2n'-2k'+1)(m,n) b(s,k)b(t,n'-k')
    
      G_3(m,n;s,t)   ≡   1/4[ F̃(m+s,n+t)+ F̃(m+s,n-t)- F̃(m-s,n+t)- F̃(m-s,n-t)]
       =   ∑_n'=0^∞Δ^2n'+1∑_k'=0^n' F_(2k'+1,2n'-2k')(m,n) b(s,k)a(t,n'-k')
    
      G_4(m,n;s,t)   ≡   1/4[ F̃(m+s,n+t)- F̃(m+s,n-t)+ F̃(m-s,n+t)- F̃(m-s,n-t)]
       =   ∑_n'=0^∞Δ^2n'+1∑_k'=0^n' F_(2k',2n'-2k'+1)(m,n) a(s,k)b(t,n'-k')

 Let us study the G_1 case (the others follow straight forward). First, let us cut the infinite sum on the right hand side by keeping only up to order Δ^2N_0. That is, 
 
    G_1(m,n;s,t)=∑_n'=0^N_0Δ^2n'∑_k'=0^n' F_(2k',2n'-2k')(m,n) a(s,k)a(t,n'-k')+ O(Δ^2N_0+2)

 The unknowns are the set of derivatives F_(n,l) and the data are G_1(m,n;s,t) by varying the values of s and t (being always integers). We have L=(N_0+1)(N_0+2)/2 unknowns and then we should give a set of L-values of pairs (s,t) in order to have a linear set of equations to be solved. We choose the set of natural values (s,t) such that s+t=m with m=0,…,N_0 and s=0,…,m. 
 The tricky issue is to codify the index (s,t) and/or (k',n'-k') in order to create a new index α running from 1,…,L. It is easy to define α as a function of (s,t):
 
    α=m(m+1)/2+s+1    , m=t+s    , m=0,1,…,N_0    , s=0,1,…,m

 where α=1,…,L. We also need the inverse relation, that is, (s,t) as a function of α. Let us assume that s=0 in eq.(<ref>). In this case m=(-1+√(1+8(α-1)))/2. For s>0 the value of m should be the same, therefore
 
    m(α)=Int(-1+√(1+8(α-1))/2)     s(α)=α-1-m(α)(m(α)+1)/2

 Therefore we can write eq.(<ref>):
 
    G_1(m,n;s(α),t(α))=∑_β=1^L  a(s(α),s(β))a(t(α),t(β)) Δ^2n(β) F_(2s(β),2t(β))(m,n) +O(Δ^2N_0+2)

 where t(α)=m(α)-s(α). With this codification we manage to keep all terms in a given perturbation order. By inverting these  set of linear equations we can express F's as a function of G_1 up to order Δ^2N_0.
 Similar strategy can be followed for G_2, G_3 and G_4. The only change is that b(0,t)=0 and we should use a little different set of index (s,t). In particular
 
    G_2   =   G_2(s(α)+1,m(α)-s(α)+1)
    
     G_3   =   G_3(s(α)+1,m(α)-s(α))
    
     G_4   =   G_4(s(α),m(α)-s(α)+1)

 For example, if we choose N_0=1 we get the first derivatives:
 
    F(x̃(m),ỹ(n))   =   1/24[-F̃(m-1,n)-F̃(m,n-1)+28 F̃(m,n)-F̃(m,n+1)-F̃(m+1,n)]
       +    O(Δ^4)
    ∂_x F(x,y)|_x̃(m,n),ỹ(m,n)   =   1/48Δ[5 F̃(m-2,n)+F̃(m-1,n-1)-36 F̃(m-1,n)+F̃(m-1,n+1)
       -   F̃(m+1,n-1)+36 F̃(m+1,n)-F̃(m+1,n+1)-5 F̃(m+2,n)]+O(Δ^4)
    ∂_y F(x,y)|_x̃(m,n),ỹ(m,n)   =   1/48Δ[5 F̃(m,n-2)+F̃(m-1,n-1)-36 F̃(m,n-1)+F̃(m+1,n-1)
       -   F̃(m-1,n+1)+36 F̃(m,n+1)-F̃(m+1,n+1)-5 F̃(m,n+2)]+O(Δ^4)
    ∂_x^2 F(x,y)|_x̃(m,n),ỹ(m,n)   =   1/Δ^2[F̃(m-1,n)-2 F̃(m,n)+F̃(m+1,n)]+O(Δ^2)
    ∂_y^2 F(x,y)|_x̃(m,n),ỹ(m,n)   =   1/Δ^2[F̃(m,n-1)-2 F̃(m,n)+F̃(m,n+1)]+O(Δ^2)
    ∂_x∂_y F(x,y)|_x̃(m,n),ỹ(m,n)   =   1/96Δ^2[44 (F̃(m-1,n-1)-F̃(m-1,n+1)-F̃(m+1,n-1)+F̃(m+1,n+1))
       -   5 (F̃(m-1,n-2)-F̃(m-1,n+2)-F̃(m+1,n-2)+F̃(m+1,n+2))
       -   5(F̃(m-2,n-1)-F̃(m-2,n+1)-F̃(m+2,n-1)+F̃(m+2,n+1))]
       +   O(Δ^4)

 
 We have checked the effect of using F̃(m,n) or F(x̃(m),ỹ(n)) for different observables. In figure <ref> we show as examples the differences between the cell-values and the point values (obtained by the above expressions) for the  temperature and the x-component of the hydrodynamic velocity fields (defined below) for T_0=17 and g=10. We observe how the average value (over the cells) of the differences are -0.0004 (0.0009) and 0.0000 (0.0003) and its standard deviations are 0.008 and 0.003 for T and u_1 respectively. The statistical errors are of order 0.1 and 0.02 in these cases. Moreover, the distribution of the differences are quite random and the normalized distribution (subtracting the average and dividing by the standard deviation) has a third central moment of -0.056 and -0.023 and a kurtosis of 0.667 and  0.589 for T and u_1 cases respectively. That is, we can use F̃(m,n) with the confidence that the 
 corrections do not introduce systematic deviations to the average values because it is just a small noise to be added to the already larger statistical error.  Observe that in <cit.> we measured the global hydrodynamic velocity moments with errors smaller than 10^-3 and it was necessary there to use this corrections to do a correct analysis. 
 
 Nevertheless we have build a well defined scheme to obtain the derivatives of an observable  at any given central cell point. That is going to be very useful when comparing with the Navier-Stokes equations.
 
  
 








§ COMPUTER SIMULATION RESULTS: EQUILIBRIUM


At equilibrium, we know that a system at temperature T,  the pressure and the density are related by the equation of state and the hydrostatic formula, that in our system can be written:

    Q(y)   =   T ρ(y) H(ρ(y))
    dQ(y)/dy   =   -gρ(y)

where Q(y)=P(y)π r^2 and we have considered the mass of the disk to be one: m=1. This expressions  and the boundary conditions, determine the behavior of Q and ρ as a function of the height y. 

The analytic expression of H(ρ) for the hard disk system is unknown. However it is well known that the Henderson expression for H(ρ)=(1+ρ^2/8)/(1-ρ)^2 is a very good approximation of the equation of state (EOS) in the range of ρ∈[0,0.5] with an relative error smaller than the 1%. We  used Henderson's EOS  to solve the differential equation and we got:

    Q   =   Tρ1+ρ^2/8/(1-ρ)^2
    
    y   =   -1/g[logρ -7/8log(1-ρ)+7/8(1-ρ)+9/8(1-ρ)^2+C]

where the constant C is fixed by giving a point: (y_0,ρ_0).
We can check the behavior of our computer program  by measuring the stationary density and pressure profiles for the equilibrium case T_0=T_1=1 and to compare with the exact results.

The local areal density, ρ(n,l),  at cell  (n,l) is computed as:

    ρ(n,l)=N(n,l)π r^2/Δ^2

where N(n,l) is the average number of particles with its center at the cell (n,l) at the stationary state, r is the radius of the particles and Δ is the side length of the cell (Δ=1/30 in all the simulations). This definition is computationally very fast and convenient by its simplicity but it has a drawback when computing the density at the boundary cells: they present a systematic deviation due to the surface exclusion around the walls and therefore, the density is underestimated at these cells. Therefore in our analysis we should exclude the boundary cells and even their neighbor cells to minimize these effects.



We use the virial theorem to compute the pressure field P(x,y). Let us first derive the expression of P for the hard disk case at equilibrium with no external field. Let α be a square spatial region of side L where there are N interacting particles of mass m. The virial theorem states:
 
    ⟨∑_i=1^N r_i· F_i⟩=-2⟨ E_c⟩

 where r_i is the vector position of particle i, F_i is the total force acting on particle i, E_c is the system kinetic energy observable and ⟨·⟩ is the equilibrium ensemble average. A simple proof can be given using the ergodic theorem:
 
    ⟨∑_i=1^N r_i· F_i⟩   =   lim_τ→∞1/τ∫_0^τ dt ∑_i=1^N r_i(t)· F_i(t)
       =   lim_τ→∞m/τ∫_0^τ dt ∑_i=1^N r_i(t)·d^2 r_i(t)/dt^2
       =   -lim_τ→∞m/τ∫_0^τ dt ∑_i=1^N (dr_i/dt)^2=-2⟨ E_c⟩

 The forces acting over each particle are the sum of two terms: internal and external, that is
 
    F_i=∑_j≠ iF_ij+F^ext_i

 Let us include this fact into the left hand side of the virial theorem. First, let us assume that the external force only acts at the system boundaries and its force is perpendicular to the box side with constant magnitude F. Then:
 
    ⟨∑_i=1^N r_i· F^ext_i⟩=(x_R-x_L-y_U+y_D)F̃=-2LF̃=-2SP

  where x_L,D and y_U,D are the coordinates of the box sides and F̃ is the average total force applied on each side. This average force is assumed to be constant in magnitude because it is assumed that the system is in mechanical equilibrium. Finally we have defined the mechanical pressure P=F̃/L.
 
 
 The inter-particle forces give the following contribution:
 
    ⟨∑_i=1^N r_i·∑_j≠ iF_ij⟩=1/2∑_i,j⟨(r_i-r_j)· F_ij⟩=∑_⟨ i,j⟩⟨(r_i-r_j)· F_ij⟩

  where we have use the fact that F_ij=-F_ji ad ⟨ i,j⟩ is the set of different pair of particles. For hard disks, the interaction only occurs when a pair of particles collide. Let us assume that the collision occurs in a time interval 2ϵ very small. Then we can write F_ij
as the variation of the linear moment on such time interval:
  
    F_ij(t_n)≃1/2ϵ(p_i(t_n+ϵ)-p_i(t_n-ϵ))

  where t_n is the time where the collision occurs. Because the collision is elastic, the moment along the vector that connects the center of the disks are exchanged between the particles:
  
    r_ij· p_i(t_n-ϵ)=r_ij· p_j(t_n+ϵ)

 where r_ij=r_i-r_j. Therefore,
 
    r_ij· F_ij=1/2ϵ(r_ijp_i(t_n+ϵ)-r_ijp_i(t_n-ϵ))=1/2ϵr_ij· p_ij

 where p_ij=p_i(t_n+ϵ)-p_j(t_n+ϵ). Finally, we can write:
 
    ⟨∑_i=1^N r_i·∑_j≠ iF_ij⟩=lim_τ→∞1/τ∑_n:t_n∈[0,τ]r_ij· p_ij

 where the sum is over all the pair collisions occurring in the time interval [0,τ].
 
 The virial theorem for hard disks can be written:
 
    -2SP+lim_τ→∞1/τ∑_n:t_n∈[0,τ]r_ij· p_ij=-2NT

 where at equilibrium ⟨ E_c⟩=NT being T the system temperature. Finally:
 
    P=ρ T/π r^2+lim_τ→∞1/2 S τ∑_n:t_n∈[0,τ]r_ij· p_ij

 where ρ=Nπ r^2/S.

When the external field is turned on, and assuming that the conditions to derive the hydrostatic formula hold, we can use the virial expression locally. That is so because we showed that the system is in a kind of local equilibrium with a given temperature and a local mean density and the mechanical pressure obtained with the virial theorem is, in such conditions, equal to the thermodynamic pressure. That is

 
    P(x,y)=ρ(x,y) e_c(x,y)/π r^2+lim_τ→∞1/2 Δ^2 τ∑_n:t_n∈[0,τ]r_ij· p_ij

where Δ=1/30 is the side length of a cell, e_c(x,y)=⟨ E_c(x,y)⟩/⟨ N(x,y)⟩ is the average total kinetic energy in the cell and the sum runs over all particle-particle collisions that occur at the cell (x,y) at the time interval [0,τ] assuming that at time 0 the system is at the stationary state.  



We have measured the cell values for ρ and P for the equilibrium case  (T_0=T_1=1). Then, we have done averages along the x-rows to build the ρ(y) and P(y) functions. We discarded the two first boundary cells in x and y. We show in figure <ref> the obtained results and the comparison with the theory. First we plotted Q(y)/T=π r^2P(y)/T vs ρ(y) for the 26 different values of y and g=0, 5, 10 and 15. We observe how all the data scale in a universal curve. Moreover, the Henderson EOS fits very well such curve. That means that the local equilibrium used in the derivation of the hydrostatic formula is correct and we are capable to capture such property even we are using very small sized cells. That could be due to the strong chaotic disk behavior where, in practice, the spatial and temporal correlations decay very fast. Second, we see how the profiles of ρ(y) and Q(y) for different values of g also follow the solutions obtained with the hydrostatic formula with the Henderson's EOS. 

We also study the relative local fluctuations of the number of particles:

    K(x,y)=⟨ (N(x,y)-⟨ N(x,y)⟩)^2⟩/⟨ N(x,y)⟩

where N(x,y) is the observable number of particles at cell (x,y). We know from equilibrium statistical mechanics that 

    K(x,y)=⟨ N(x,y)⟩/Δ^2T(x,y)κ_T(x,y)

where Δ is the cell side length (Δ=1/30 in our case). κ_T(x,y) is the isothermal compressibility that it may be derive from the EOS and for hard disks case it can be written

    κ_T(x,y)=π r^2/Tρ(x,y)1/H(ρ(x,y))+ρ(x,y) H'(ρ(x,y))

with H'(ρ) being the derivative of H. As in the pressure case above, we have done averages along the x-rows and we have discarded the two first boundary cells. We show in figure <ref> the data obtained for the equilibrium case for different values of g.  We also compare with the theoretical case using the Henderson's EOS. We observe (1) all the data collapse over a common curve which implies that the local equilibrium is again correct and (2) it deviates significantly from the theoretical value. We thing that this is due to the size of the cell in which we measure the fluctuations. In our case Δ≃ 4.09 r which is too small to consistently capture the local pair correlation behavior as it were an infinite system. Therefore we are seeing strong size effects due to the size of our virtual cells. A complete discussion about the behavior of K for hard disks at equilibrium when measuring it on small cells can be found on <cit.>. In order to get reasonable results we would need to simulate systems with at least 10^7 particles which it is far from our actual computer capabilities. We  conclude that there is no reason to study local fluctuations of any magnitudes due to such large finite size effects. In any case our simulation is able to describe with precision the local equilibrium behavior of the hard disk gas and the averages of the main magnitudes as the density, temperature and pressure. 






§ COMPUTER SIMULATION RESULTS: GLOBAL MAGNITUDES


We have measured several global magnitudes like  the averaged kinetic energy (e_c),  the potential energy (e_p), the hydrodynamic kinetic energy (e_u), the pressure (P) and the energy current across the boundaries (J). All of them will give us a general overview of the system behavior. Moreover, we are particularly interested in knowing if the trasition non-convective to convective  appears or not as an abrupt change on their functional behavior and they can give us some clue to determine with precission the transition point. 



  * The averaged kinetic energy: The averaged kinetic energy for a given disk configuration at time t is given by:

    e_c(p,t)=1/N∑_i=1^N p_i(t)^2/2m

We have computed its time average and its second moment:

    e_c=1/M∑_t=1^M e_c(p,t)≡⟨ e_c(p,t)⟩  ,    m(e_c,2)= ⟨ e_c(p,t)^2⟩- ⟨ e_c(p,t)⟩^2

We observe at figure <ref> (left) that the averaged kinetic energy grows monotonically with T_0 in a smooth nonlinear way. The kinetic energy per particle grows with g for any fixed T_0 value. That is, for a fixed external gradient (non zero) any positive variation of the external field increments the system kinetic energy and this hypothetical process  does a net work over the system (speaking in terms of classical equilibrium thermodynamics). Let us remark that variations of g does not affect the value of the averaged kinetic energy at equilibrium  (T_0=T_1=1) that is e_c=1. Therefore, the dependence on g of e_c  is just a pure nonequilibrium property.
The slope of e_c at T_0=1 grows with g:m=0.498 (0.006), 0.852 (0.006), 1.01 (0.01) and 0.997 (0.006) for g=0, 5, 10 and g=15 respectively. 
 For values  T_0>15  e_c  tends to an asymptotic linear behavior with a almost common slope of m'=0.33 for any g value. In fact we get m'=0.330 (0.006), 0.330 (0.006), 0.320 (0.006) and 0.34 (0.02) for the g=0, 5, 10 and 15  cases respectively. That is far from the limiting slope associated to a pure linear temperature profile: e=(T_0+T_1)/2 that has a slope equal to 0.5. It could be though that the thermal resistance at the boundaries makes the system to behave in the bulk with  effective temperatures T_0^* and T_1^*. If that was so we could identify: m' (T_0)=(T_0^*(T_0)+T_1^*(T_0))/2. We'll see later that T_0^*(T_0) grows linearly with T_0: T_0^*≃ 0.17 T_0 and T_1^*(T_0) tends to a constant value. Therefore m' does not match with the linear profile assumption and we can conclude that for large T_0 values there is not a linear profile of temperatures.  

Finally let us to remark that there is no trace of any discontinuity and/or singularity in the data.  Therefore this macroscopic observable does not reflects any singular property on the convective-non-convective transition. 

The second moment relative to e_c: m(e_c,2)/e_c^2  behaves in a more interested manner. At T_0=1 (equilibrium) its value corresponds to 1/N=1/957=0.00104.. as the law of large number predicts. We see at figure <ref> (bottom) how the relative fluctuations, for a fixed g>0 value, decrease as we increment T_0, it reaches a minimum and then it grows. This is in contrast with the equilibrium case (g=0) in which the relative fluctuations grow monotonically with T_0. 
We get an idea on where it is located the minimum by fitting a 9th degree polynomial to the data (plotted by a dotted line in the figure) and equating to zero its derivative. We find that the minimum values are at  T_0^c=2.2, 3.3  and 4.6 for g=5, 10 and 15 respectively. We will see that this singular point where the relative kinetic energy fluctuations  has the smallest value seems to be related to the critical temperature that defines the transition from the non-convective to a convective regime. However these critical temperatures does not coincide with the ones we already computed in the Enskog approximation (dashed vertical lines on figure <ref> (right)). Moreover our critical values increase with g while the ones from the Rayleigh parameter decrease with g. One explanation to this discrepancy is that the Rayleigh critical value is obtained in the Boussinesq approximation that is assuming that the fluid is non-compressible. Some numerical analysis of the linearized Navier Stokes equations found that the critical temperatures for the onset of convection, increase with g <cit.>. In fact, the increment is argumented to be given by gα T L_y^4/c_p. Using the values with the Henderson's EOS and ρ̅=0.2, we find that the critical temperatures are: T_0^c (comp)=3.2, 4.6 and 5.9 for g=5, 10 and 15 respectively. These values are not too far from the ones we have found and the linear dependence on g is similar to ours which is quite remarkable.

Let us finish these arguments by pointing out what would happens it the parameter A(T_0,g)=m(e_c,2)/e_c^2 that localizes the critical value for the temperature were correct. Then one could conclude that:  

    dA/dT_0   <   0    if   T<T_0^c
    dA/dT_0   >   0    if   T>T_0^c

and, at the critical temperature  the following relation should be true:

    d m(e_c,2)/dT_0|_T_0^c=2m(e_c,2)/e_cde_c/dT_0|_T_0^c

It could be interesting to check in some other scenarios this relation to exclude the possibility that we are just seeing a numerical mirage. 




  * The potential energy: It is the energy associated to the constant external acceleration  that is applied to the particles (a⃗=-gj⃗).The potential energy per particle is then defined as:

    e_p=1/N∑_i=1^N mgr_i,2

where r_i,2 is the vertical coordinate associated to the particle i assuming that at r_i,2=0 there is located the thermal bath at temperature T_0. We observe on fig. <ref> that e_p monotonically grows with T_0. That is, the center of mass of the system increases its height when we increase the temperature at the bottom. From the data we can get the asymptotic value fo e_p when T_0→∞ by fitting a polynomial of the form: e_p=e_p^a+a_1/T_0+…+a_5/T_0^5, obtaining: e_p^a=3.21 (0.09), 6.56 (0.09) and 8.9 (0.3) for g=5, 10 and 15 respectively. We cannot exclude a linear dependence of e_p^a on g. One may easely check that the potential energy per particle of a uniform disk system in which all the particles were at the top is e_p^t=g(1-ρ/2ρ^t) and with our data ρ^t≃ 0.3 which is far from the liquid-hexatic phase transition that is at ρ_c≃ 0.70 <cit.> . That is, even in the infinite T_0 limit, the system would still be at a liquid-like phase.

Finally we have found that the slope of the fitted function of the data near the equilibrium point  is found to be m=0.573 (0.006),  1.12 (0.01) and 1.26 (0.02) for g=5, 10 and 15 respectively. Observe that in this case that the slopes are not compatible with a linear dependence on g. Just for completeness, the values at T_0=1 are e_p=1.7365 (0.0003), 2.6263 (0.0005) and 3.343 (0.001).

 
 


  * The hydrodynamic kinetic energy: It is defined as the kinetic energy associated to the velocity of the center of mass at each cell. That is, first lets us define  the average velocity measured at cell (n,l):

    u(n,l;M)=1/N(n,l)M∑_t=1^M ∑_i:r(i,t)∈ B(n,l) v_i(t)    ,    N(n,l)=1/M∑_t=1^M ∑_i:r(i,t)∈ B(n,l)1

where we have follow the local average definition in eq. <ref>.  Then, we define the hydrodynamic kinetic energy:

    e_u(M)=1/2N_C∑_(n,l)ρ(n,l)/ρu(n,l;M)^2

where ρ(n,l) is the mean areal density at cell (n,l), N_C is the total number of cells and the M argument indicates that we have done a large but finite number of measurements.

First observe that the measured values are three orders of magnitude smaller compared to the total kinetic energy or the potential energy ones. This is due to the hydrodynamic separation of scales. This makes the numerical analysis of the convective structures very hard and that's why we need very large time averages to get some clear picture of the system behavior.

In general we expect that u(n,l;M)=0, ∀ (n,l) for the non-convective states, that is for T_0<T_0^c and  u(n,l;M)≠ 0 for the convective states  (T_0>T_0^c).  As we have commented in the error analysis section above, the finite number of measurements makes that u(n,l;M) is typically nonzero due to the small but systematic effect of the noisy data behavior. We should correct this fact using the expressions <ref> and <ref> that in our case are given by:

    e_u≡ e_u(M→∞)≃ e_u(M)- 1/18ρ N_C∑_(n,l)ρ(n,l)ϵ(u(n,l))^2 ±ϵ(e_u)

with

    ϵ(e_u)=1/2ρ N_C[ ∑_(n,l)u(n,l)^2(4ρ(n,l)^2ϵ(u(n,l))^2+ϵ(ρ(n,l))^2u(n,l)^2)  ]^1/2

where
ϵ(u(n,l)) is the measured error of u(n,l)  and ϵ(ρ) is the measured error fo ρ(n,l) (both are three times their standard deviation).

Observe that we subtract to the measured value of e_u(M) a small but monotonous increasing  function of T_0. One checks that the values so obtained have a nice average of zero near T_0=1 and we can try now to fit an appropriate function. We choose  e(u)_fit=0 for T<T_0^c and =(T-T_0^c)^2 (a_0+a_1T_0^2+… a_7 T_0^8) for T_0>T_0^c where the fitting parameters are T_0^c and a's.  We choose (T-T_0^c)^2 as a minimal assumption of asking continuity and first derivative zero both at T_0^c. We fit this function to the points in the intervals T_0∈[2.6,12], T_0∈[4,13]  and T_0∈[5,16] and we obtain: T_0^c=2.1, 3.4 and T_0^c=4.3 for g=5, 10 and 15 respectively. These values change a bit when choosing other polynomials degree and/or the interval of fitting data. However they are consistent with the results obtained from the energy fluctuations and we think that the critical temperature should be around such values. 
 

Finally it is natural to define a kind of anisotropy parameter that measures the relative difference on the x and y components of the hydrodynamic kinetic energy.
We define

    m_u=1/2 e_uN_C∑_(n,l)ρ(n,l)/ρ[u(n,l;M)_2^2-u(n,l;M)_1^2]

First observe that in the non-convective region this observable is singular because the hydrodynamic velocities are zero and this is reflected with a very large set of error bars (that are computed in a similar way as in the e_u case (see <ref> and <ref>)  for the data in that region. The curious result is that this parameter seems to stabilize to a constant value ≃ 0.55  independently (up to error bars) of the T_0 and g parameters. That is, the relative anisotropic behavior seems to be universal with respect the boundary conditions (T_0) and the external forcing (g). This  implies that there is a kind of energy equipartition:

    e_u,2-e_u,1=0.54(e_u_1+e_u,2)     e_u,2≃ 3.4  e_u,1

independently on g and for large enough values of T_0.

 

  * The Pressure: We have computed  the system pressure on each of the thermal bath boundaries as the sum of the moment variation when particles collide with  the bath divided by the absolute time measuring interval times the boundary length:

    P=m/Δ t L_x∑_t_col∈Δ t(v_i,2'(t_col)-v_i,2(t_col))

where v' and v are the particle velocity right after and right before the collision respectively. We first observe that the global barometric formula: 

    Q(0)-Q(1)=g Nπ r^2   ,    Q=Pπ r^2

holds in all cases as figure <ref> shows with high precision. This property is independent on the state of the fluid: non-convective or convective. However, the pressure has a non-trivial dependence with T_0. When g=5 the data is concave over all the T_0 interval. The slope at T_0=1 is 0.1125 (0.0004) and it gradually diminishes until its asymptotic value 0.10 (0.01). However, for g=10, the data has a small slope of 0.1740 (0.0004) at  T_0=1 and it is convex near it, it changes convexity at around T_0=10 and it finally tends to the same asymptotic behavior as in the g=5 case: Q≃ 0.10 (0.01)T_0. Finally, for g=15 the slope at T_0=1 is zero and, again, it tends to the same asymptotic value. Observe that the Q at T_0=1 tend to have the value gρ̅ because the pressure at the top goes to zero. We do not appreciate any singularity in the curve near the transition between non-convective and convective regimes.  
 

  * The current of energy: 
We have computed  the energy current that cross each of the thermal bath boundaries as the sum of the kinetic energy variation when particles collide with  the bath divided by the absolute time measuring interval and the boundary length:

    J̅=1/2Δ t L_x∑_t_col∈Δ t(v_i,2'^2(t_col)-v_i,2^2(t_col))

We see on figure <ref> the averaged behavior of J=J̅ r m^1/2  as a function of T_0 for g=0, 5, 10 and 15. First we observe that J(0)=-J(1) as expected because there isn't any internal energy dissipation mechanism. Second, its behavior is not linear in T_0 implying that one should assume that the thermal conductivity cannot be a constant if we assume that the Fourier's Law applies. Third, for large values of T_0 we see that the data seem to follow the same asymptotic behavior  independently of g:

    J≃   m_1 T_0^3/2

where m_1=0.4726 (0.0003), 0.4870 (0.0009), 0.559 (0.005) and 0.672 (0.002) for g=0, 5, 10 and 15 respectively. 

The current presents a slope at T_0=1 that dependes on g: m=0.81, 0.58, 0.08 and 0.04 for g=0, 5, 10 and 15 respectively.



From the behavior of the heat currents in figure <ref> we observe some interesting physical properties:



  * (1) Any curve J(T_0;g) with a fixed g value have a (non trivial) crossing point with the curve J(T_0;g=0). We call such point T_c,2 and it is solution of the equation J(T_c,2;g)=J(T_c,2;0) (see Table <ref> for computed values). We can conclude this property by seeing the behavior of the slopes of J(T_0;g) near T_0=1: m(g)<m(g=0) for all g's. T_c,2(g) separates neatly two regimes for any fixed g: (1) the temperatures at which the heat current is smaller than the one without gravity, T_0<T_c,2(g). We will call this a bad conducting state because gravity is penalizing heat conduction and (2) the temperatures in which the sytem has a larger current than the corresponding without gravity,  T_0>T_c,2(g). We call it fully convecting state because the action of g through convection is increasing the energy transmission with respect the g=0 case.


  * (2) For any fixed value of T_0, there exists a g_c in which its heat current equals the one of a system with g=0: that is J(T_0;g_c)=J(T_0;0). Again we can observe two different regimes: (1) For g<g_c  J(T_0;g)>J(T_0;0) and (2) g>g_c  where J(T_0;g)<J(T_0;0). Observe this behavior in figure <ref> (left) in a particular case. First we found that the crossing point between the currents with g=10 with the ones with g=0 is at T_c,2=4.96921 with a common heat current of J=5.0315. Then we use this temperature to know the currents for the cases g=5 and 15. Finally we can also conclude that for any fixed value of T_0 there exists a  g^*<g_c  at which the heat current reachs a maximum.


  * (3) We plot at figure <ref> (right) the temperatures at which the hydrodynamic kinetic energy begin to be nonzero (T_c) together with the temperatures at which the heat current equals the one without gravity (T_c,2). We observe an interesting pattern. It seems that for any fixed g>0 there are two critical temperature values that separates three regions: (A) T_0<T_c where there is no macroscopic velocity field and there is only heat conduction (there are bad conducting states). (B) T_0∈[T_c,T_c,2] where there is a small hydrodynamic field but still the heat current is less than the correponding with g=0 and (C) we have hydrodynamic velocity fields and the states are fully convecting ones (in the sense we defined above). If this is the general picture we would like to know if, for instance, there is a coherent convective behavior in region B.


 

 
 



§ SPATIAL STRUCTURES


Once clarified (more or less) the existence of two well defined regimes: non-convective and convective, we can study the spatial structures of the local magnitudes. In general, as Navier-Sokes equations predicts,  we expect that a given local observable A  should be only function of y for the non-convective states, A=A(y), and  A=A(x,y) in the convective states. We are going to analyze among others the behavior of the local  hydrodynamic velocity, temperature, areal density and the virial pressure.



 §.§ THE HYDRODYNAMIC VELOCITY FIELD
 

The velocity field is measured as the average center of mass velocity, eq.(<ref>), at each of the 900 (30× 30) virtual cells in which we have divided our system. The modulus of the hydrodynamic velocity is typically very small compared with the averaged mean particle velocity (for instance, for the case g=10 and T_0=20, the largest cell averaged value is 0.11 compared with the mean particle velocity of around 3) Therefore we needed very large time averages to resolve the structure beyond the underlying chaotic like behavior of the particles. Moreover, as we will see, the field structure can only be analyzed with care far from the critical temperature where we can discriminate the field from fluctuations.  

We see in figure <ref> some snapshots of the configurations for g=5, 10 and 15 cases and five different temperatures.They have been chosen with the following criteria: the first one is below T_c, the second between T_c and T_c,2, and the rest beyond T_c,2. It is clear that the disordered behavior of the configurations below T_c (first row). For temperatures between T_c and T_c,2 some incipient local order is devised but fluctuations seems to dominate along the system. For large temperatures two rolls of convective fluid develops.  For the g=5 the rolls are noisier than the ones of the g=10 or g=15 cases,  in particular near the hot thermal bath. The color of the arrows indicate the velocity magnitude: from the smallest value in the given configuration (dark blue) to the largest value (dark red). We observe how the fluid goes up in the middle of the system with a intermediate velocity and the  fluid goes down following the vertical boundaries with higher velocity. Note that the hot thermal bath gives to the particles a lot of energy but it does not directly contribute to the coherence of the macroscopic movement of the fluid. The external field g is the necessary mechanism that activates such movement: it seems to us that it is essential the acceleration that suffers the particles that goes the way down that, together with the neighbor interactions, tends to align the hydrodynamic velocity vector on the y-direction. This ordered fluid lose part of it acquired momenta when interacting with the disordered hot region but maintains most part of its coherence. Moreover, the created ordered structure is capable to trap and align (in average) the very hot particles from the hot reservoir and the fluid accelerates even going against the external field.   





  §.§.§ Stream lines


In order to analyze the onset of convection from the spatial structure of the velocity vector field we have computed the distribution of its stream lines for each case. We define stream line as the trajectory build from an initial point whose point tangents are the given fixed vector field. That is, let u⃗=(u_1(x,y),u_2(x,y)) a given vector field and (x_0,y_0) an arbitrary initial point. Then, the stream line is solution of the differential equation

    dy/dx=u_2(x,y)/u_1(x,y)

or in parametric form:

    dx/ds=u_1(x(s),y(s))   ,  dy/ds=u_2(x(s),y(s))

 The numerical solutions can be found just by using a Runge-Kutta integrator forward and backwards in time from the initial point. In our case we have a vector field defined over a grid. Then, we can reconstruct the vector field at each integration step and in a given point by linear interpolating it from neighboring grid sites. Because the vector field is not generated by an analytic function and there are noise effects, we should include a stop condition on the integration: u⃗(n+1)·u⃗(n)<0 where u⃗(n) is the vector used in the n-th step integration (from point r⃗(n) to the next one r⃗(n+1)) and u⃗(n+1) is the interpolation vector computed at the point r⃗(n+1) (integration of stream curves from a discrete vector field are used in imaging processing, see for instance ref. <cit.> for some other technicalities). For any given vector field we generate 10000 stream lines with random chosen initial point uniformly distributed inside the vector grid and we compute the statistics of the path-lengths. In Table <ref> (left column) we show a small sample of 100 stream lines corresponding to the g=10 case  and temperatures T_0=1.8 (A), T_0=4 (B), T_0=6 (C), T_0=10 (D) and T_0=20 (E). We choose these points on purpose to show the typical behavior at each of the regions.
 




The averaged stream line length versus T_0 for g=5, g=10 and g=15 is shown in figure <ref>. In the non-convecting region it gets values of about 0.1 and it grows softly until it reaches T_c when it grows fast until it reaches the value 2 at a T_0≃ 10 from which it fluctuates around. Observe that 3 is the perimeter of a rectangle of sides 1× 1/2. The big peak around T_0=10 for g=10 is just due to the discrete structure of our vector field that favors the appearance of large spirals when our algorithm is applied  (think in a closed trajectory with all vectors aligned to it except for one  pointing out). We can explicitly see this effect in the stream trajectories  sample corresponding to the point D in figure <ref>. For large T_0 values (point E) we have mainly closed stream lines and the spirals tend to disappear. We are pretty sure that if we do averaging over other stationary vector field realizations (assuming that there is a fluctuating part on it due to the finite time averaging used) and we simulate much larger systems would erase such peaked behavior. The variance of the stream length distribution is shown in figure <ref> where  we again observe the large variance around T_0=10 for g=10 indicating the strong sensibility on the length to the initial chosen random point.

The distribution of lengths changes its nature from the non-convective to convective regions as we see in Table <ref>. First observe that there are a set of initial seeds whose stream length is less that 1/200. That is, the algorithm stops just at the first steps. Such points correspond typically to strong noisy regions where local hydrodynamic velocity fields are dominated by microscopic fluctuations. We show at figure <ref> the proportion of points with length larger than 1/200. In all cases, the noisy regions tend to diminish when we increase T_0. For T_0<T_c (non-convecting region) there are always about 8% of such points. The proportion diminishes to zero for T_0>>T_c with a sigmoid like of behavior. It is curious that T_c,2 seem to show, in all cases, the typical intermediate point where the proportion of strong noisy points is of about 3% 




For the non-convective region (T_0<T_c) and once discarded the strong noisy data points set, the remainder distribution is reasonable well fitted by a gamma-distribution:

    P(l;α,β)=1/β^αΓ(α)l^α-1e^-l/β

 as we show with the red curve in Table <ref> (A). The α-fitted values are about 2.5 (see figure <ref>). The β fitted values grows with T_0 similarly as ⟨ l⟩ does (in fact it is just a natural length scale). The Gamma distribution fit  fails when we analyze the data corresponding to convecting states (T_0>T_c). The point (B) with T_c<T_0<T_c,2 has a distribution that looks like a shifted-Gamma distribution with a bump that reflects the existence of large stream lines (see for instance the red curve in Table <ref> (B)). For T_0>T_c,2 the Gamma distribution is totally lost and there appears a extended distribution with some nontrivial peaked structure. We see in Table <ref> (C) how there are a collection of “channels”  where the long stream lines focus. They are responsible of the distribution expansion to higher l-values and to the existence of local peaks on it. The distribution at the point (D) confirms our idea about the existence of spirals. We see there how the probability of having stream lines larger than 3 is non-zero. Even there are lines with length around 6  which implies the existence of trajectories wrapping around half of the system several times. We recover a more physical behavior for a fluid  at the point (E) (T_0=20) where the distribution drops fast at l≃ 3. We can construct a simple model to understand, in some way, this limiting stream line distribution. 




Let us assume that there are N closed stream lines turning around the point P:(x_0,y_0). We also assume that along the principal axis with center at P, the lines are uniformly distributed. This implies that their crossing point at the axis with origin at P are mutually separated by (1/2-x_0)/N along the PX positive axis, x_0/N along the PX negative axis, (1-y0)/N along the PY positive axis and y0/N along the PY negative axis. Stream lines should also be perpendicular to the axis at the crossing points. With this set of assumptions we just have to build a family of self avoiding curves connecting the crossing points. The idea is to build them with some structural parameter in such a way that the limiting value of such curves should be the containing rectangle. 

The simplest example is to assume that the stream lines are just rectangles (see figure <ref>). In this case the n-th stream line is built of 8 pieces, two for each of the four quadrants:

    y=y_0+n (1-y_0)/N   x∈[x_0,x_0+n(1/2-x_0)/N]       ;      x=x_0+n(1/2-x_0)/N    y∈[y_0,y_0+n(1-y_0)/N]
    
    y=y_0+n (1-y_0)/N   x∈[x_0-n x_0/N,x_0]       ;      x=x_0-n x_0/N    y∈[y_0,y_0+n(1-y_0)/N]
    
    y=y_0-n y_0/N    x∈[x_0,x_0+n(1/2-x_0)/N]       ;     x=x_0+n(1/2-x_0)/N    y∈[y_0-n y_0/N,y_0]
    
    y=y_0-n y_0/N    x∈[x_0-n x_0/N,x_0]       ;     x=x_0-n x_0/N    y∈[y_0-n y_0/N,y_0]


The n-th rectangle surface is S(n)=n^2/2N^2 with n=1,…, N. The length of the n-th stream line is just the perimeter of the rectangle: l(n)=3n/N. Then, the probability of having a stream line of length l(n) can be built as the surface diference between rectangles n and n-1:

    Q(l(n))=2(S(n)-S(n-1))=2n-1/N^2

and the probability distribution: 

    P(l)=Q(l)/Δ l=2n-1/3N=2/9l-1/3N   l=3/N, 6/N,…, 3

we can do the limit N→∞ and P(l)=2l/9 with l∈[0,3].  It is trivial to compute the average lenth and the variance:

    ⟨ l⟩= 2     ,      m_2(l)=1/2

 We can compare these results with the distribution in the case g=10 and T_0=20 (see right figure E in <ref>).
There the red line is the distribution we have obtained with our model. It is clear that we have captured some of the observed behavior: increasing probability with l, cutoff probability at l=3, average l of around 2 and variance of  around 0.5 (see figure <ref>). Obviously, the rounding structure should be due to the curvy form of the stream lines. 
We tried two other types of stream lines in order to check this point (see figure <ref>)


First we introduced a circular corner of radius R(n) to the above rectangular model of the stream lines (see the figure <ref> bottom-center):

    R(n)=f(n/N) min (L_1,L_2)   ,   f(u)=u(1-u)

where L_1 and L_2 are the length of each quadrant with origin (x_0,y_0): L_1=1/2-x_0, L_2=1-y_0 (first quadrant), L_1=1/2-x_0, L_2=y_0 (second quadrant), L_1=x_0, L_2=y_0 (third quadrant) and L_1=x_0, L_2=1-y_0 (fourth quadrant).
In this case the surface surrounded by the n-th stream line is given by:

    S(n)=n^2/2N^2-(1-π/4)θ(x_0,y_0)f(n/N)^2
 
with

    θ(x_0,y_0)= [min (x_0,y_0)]^2+[min (1/2-x_0,y_0)]^2+ [min (x_0,1-y_0)]^2+[min (1/2-x_0,1-y_0)]^2

Following the same steps we did above for the rectangular stream lines case we obtain in the N→∞ limit

    P(l)=2[u-2θ(x_0,y_0)(1-π/4)f(u)df/du](dl/du)^-1

where l∈[0,3], u∈[0,1],

    l=3u-(2-π/2)w(x_0,y_0)f(u)

and

    w(x_0,y_0)=min (x_0,y_0)+min (1/2-x_0,y_0)+min (x_0,1-y_0)+min (1/2-x_0,1-y_0)

This distribution applied to the case g=10 and T_0=20 give us ⟨ l⟩=1.9392 and m_2(l)=0.5190. Moreover, the distribution form is plotted as a blue line in Table <ref> (E) and there we observe that there is a curvature but still it is far to the observed one.

Finally we tried with a continuous curve at each quadrant. To connect them we ask to the curves be perpendicular to the crossing points at the quadrant axis with lengths (L_1,L_2). Let x_1=nL_1/N  the crossing point of the n-th stream line on the x-axis (with center at (x_0,y_0)) and y_1=nL_2/N the corresponding one with the y-axis. Then we have chosen the function 

    y(x;n/N)=y_1 ỹ(x/x_1;n/N)≡y_1/1-β(1-x/x_1)^β-β y_1/1-β(1-x/x_1)

with β=β_0(1-n/N)^β_1 to compute P(l). Observe that the function is tunned in such a form that y(0)=y_1, y'(0)=0, y(x_1)=0 and y'(x_1)=∞. The distribution obtained in the limit N→∞ is, in this case,

    P(l)=(u^2 dg/du+2ug(u))(dl/du)^-1

where

    l(u)=u h(u)

and

    g(u)   =   ∫_0^1 dzỹ(z;u)
    
    h(u)   =   (1/2-x_0)∫_0^1 dz (1+(1-y_0/1/2-x_0dỹ/dz)^2)^1/2+
    (1/2-x_0)∫_0^1 dz (1+(y_0/1/2-x_0dỹ/dz)^2)^1/2
       +   (1/2-x_0)∫_0^1 dz (1+(y_0/x_0dỹ/dz)^2)^1/2+
    (1/2-x_0)∫_0^1 dz (1+(1-y_0/x_0dỹ/dz)^2)^1/2



This distribution for β_0=0.5 and β_1=0.5 has the form plotted in figure <ref> (bottom-left). In this case we find that ⟨ l⟩=1.8495 and m_2(l)=0.4496. The overall distribution has the form shown by the green curve in Table <ref> (E-green). Observe that the distribution has a similar form compared to the measured one but shifted to the right. This is probably due to the simplified assumption of equidistant crossing points through the axis centered at (x_0,y_0). We cannot go further because the data we have obtained has a limited precision and playing with other distributions would just be a game. We wanted only to show that it is possible to do a systematic study of stream lines statistics  and we can unveil some interesting properties that characterize, from another point of view, the transition between the non-convective regime to the convective one, from a gamma distribution to a peaked distribution build from closed curved stream lines. It is an open issue to understand the transition between those extreme cases but it is beyond the computer effort done in this work. In order to do that study, one should need to focus the simulation in creating many vector configurations at the stationary state, for each one generating stream lines and, finally, do and overall precise study of the resulting distribution.
Nevertheless it could be very interesting to do such hard work. We plot in figure <ref> the ratio of vectors whose modulus are larger than its standard deviation versus the measured stream length for all g and T_0 values. We see how all the points corresponding to different g-values scale in a unique curve. That is, the average length depends on the proportion of “ordered” vector points (the ones whose modulus value is larger to the typical fluctuation interval and then it can be discriminate from the underlying noise). It seems that the stream line distribution maybe related to a kind of percolation phenomena. The overall picture from this point of view could be:


  * Non-convective regime (T_0<T_c): The proportion of non-noisy vector points (pr) are less than ≃ 0.5 and the averaged length seems to to be almost constant with pr.

  * Bad conducting regime (T_c<T_0<T_c,2): We observe that the values for which T_0≃ T_c,2 correspond to  ⟨ l⟩≃ 0.6 and pr≃ 0.68 in all cases. Notice that the 2-d percolation of overlapping random disks of radius r has a critical covered surface at the percolation threshold of ≃ 0.67 that is very near to our observed value. We can naturally connect both phenomena and we can conclude that when the set of “ordered” vectors percolates, the stream lines are long enough to connect both thermal baths and the fully convection regime appears.

  * Fully convecting regime (T_0>T_c,2): the noisy points tend to zero.

We are observing a typical fluctuating phenomena that it could disappear when the number of disks, or lattice points tends to infinity. However there are some open questions that we cannot resolve at this point: There exists the region [T_c,T_c,2] at the hydrodynamic limit? If the answer is positive it would imply the existence of disordered closed stream lines of relative short length. That is, there would be solutions of Navier Stokes equations different from the well known regular convective ones. If  the answer is no, it would imply that T_c=T_c,2 in the hydrodynamic limit and the regular convective solutions would be the unique non-zero ones beyond the pure conducting regime.





  §.§.§ The velocity vector field u=(u_1(x,y),u_2(x,y)) at the fully convective regime


We show in figure <ref> the averaged values obtained for the components of the hydrodynamic velocity field u=(u_1(x,y),u_2(x,y)) for g=10 and T_0 values corresponding to the particular cases studied in figure <ref>. Observe how above T_c there is a well defined spatial ordering for each velocity component that gets clearer as we increment the value of T_0. This favors the idea that the percolating phenomena we described in the above section appears only at the fluctuating level of description and it may disappear in the  hydrodynamic limit.




We analyze with more detail the behavior of the components of the hydrodynamic velocity in some particular case.
Let's first study u_2 as a function of l for the extreme columns n=1 and n=30 (averaged by symmetry)  and the central columns n=15 and n=16  (averaged by symmetry) for the case g=10 and T_0=19. In figures <ref> we see that the fluid goes down from l=30 to l=1 in the columns n=1 and n=30 and goes up from l=1 to l=30 in the columns n=15 and n=16. Let us remind that the center of the cell (n,l) are located at x=(n-1)/30+1/60 and y=(l-1)/30+1/60.

First observe that for such columns the orientation of the hydrodynamic velocity vectors is practically on the y-direction (see right figure on <ref>). That is, these columns  maybe considered as stream lines of the fluid.
The velocity  is almost zero at y=0 and it starts to increase as we move towards smaller y-values. It reaches its maximum value and  decrease the velocity up to zero (orange data points). Just the contrary happens with the central columns: from y=0  it first accelerates and latter it decelerates (green data points). We can fit  to such profiles  the functions:

    u_2(x,y)   =   -y (1 - y) (0.55 + 0.12 y +1.28 y^2)  ,   x=1/60, 59/60    (n=1,30)
    
    u_1(x,y)   =   y (1 - y) (0.22 + 1.00 y)     ,    x=29/60, 31/60   (n=15,16)

 We can define a kind of vertical effective  fluid acceleration as 

    g_eff=1/2∂ u_2^2/∂ y

and from the fitted functions we can get an idea of the acceleration values. For instance: g_eff=0.058 for y=0.25 and g_eff=-0.076 for y=0.75 for columns x=1 and 30 and  g_eff=0.037 for y=0.25 and g_eff=-0.054 for y=0.75 for columns n=15 and 16. All these typical values are far from the external one: g=10. That is, the fluid hydrodynamic velocity field  is having a lot of  friction  from the underlying chaotic movement of particles being at local equilibrium. From Navier-Stokes equations we see that the effective acceleration on the y-direction is the competition of the external forcing, g, the pressure gradient and the viscous terms. We'll try later to study  the role of the last two on the value of g_eff once clarified the behavior of the pressure profile.

We see in figure <ref> the behavior of u_1 as a function of x for the rows l=1 (bottom) and l=30 (top). The data near the hot bath (y=1) is very noisy compared with the data measured on the top. Again the fluid is aligned flowing along the x-direction. The profiles seem to follow a periodic function with modulated amplitude:

    u_1(x,y)= A(x)sin(a_1x+a_2))  ,   y=1

where A(x) is a second order polynomial and a'2 are fitted parameters.

Qualitatively, the general form of the u_1,2 profiles is similiar independently of the values of g or T_0.  However the amplitudes of the profiles increase with T_0 for a fixed g and with g for a fixed T_0 . This is shown in  figure <ref> where we plot the area contained on the u_1 profiles for the top row 30. We observe how the area increase as T_0 increase to a limit value (at least for the g=5 and g=10 cases). Moreover, the values are increased more than twice as we double g which implies a nonlinear behavior as a function of g that we cannot study with the actual set of computer simulations.  The existence of a limiting area implies that it should exists a limiting profile for a fixed g and T_0 large enough. That is, the convective fluid absorbs always a maximum quantity of energy (kinetic). It also indicates that the T_0 values controls a kind of resistance to the macroscopic fluid movement that always exists even for large values of T_0 where its value is minimum. The behavior with g reflects again that the power engine of the convective fluid is the external acceleration.



The idea of the existence of a limiting profile for T_0→∞ inspire us about the existence of a scaling profile at least in the asymptotic regime, that is, when T_0 is large enough. In fact, we are going to try to show below that, up to our computational precision level, there is a hydrodynamic universal profile for T_0 large enough values. We see in figure <ref> that the spatial structure of the velocity field changes with T_0: they have peaks and valleys more and more pronounced when we increase it. The most naive form to see if there is an scaled common field  is to look for a natural measure of the peak-valley differences, rescale each configuration by it and compare them. As scaling parameter we have chosen the configuration standard deviation:

    σ(u_1,2)=[1/N_C∑_(x,y) (u_1,2(x,y)-u̅_1,2)^2]^1/2  ,  u̅_1,2=1/N_C∑_(x,y)u_1,2(x,y)

where N_C is the number of cell points and the sum runs over all (x,y)-cells. We define a scaled configuration as:

    u_1,2^(s)(x,y)=u_1,2(x,y)-u̅_1,2/σ(u_1,2)



In figure <ref> we show the superposition of four different T_0 scaled configurations for g=5, 10 and 15. We see, at a glance, that except for noise effects the configurations seems to be the same. Obviously, these pictures  by themselves don't prove anything but they are a good starting point. The hard work now is to attempt to show that there is not a systematic T_0 dependence on such scaled configurations. From the numerical analysis point of view, we should be very careful. We have already seen that just a variable fluctuation introduces systematic deviations on any function of it. Then we should painfully handle all these effects. The scheme we have followed for this analysis is:


  * Preparation of the scaled fields: We'll take into account the effect of the configuration error bars in the values of σ(u_1,2). From it we'll find the corrections on the scaled fields coming from the fluctuations. 

  * Analysis of the scaled configuration spatial structure and mutual comparison: For each configuration we'll measure the eigenvalues of the inertial moments and the fourth momenta. Finally we also study the mutual distance between configurations. All of this just for looking any systematic T_0 dependence.

  * Obtaining the universal field: We'll average the set of configurations that we assume numerically indistinguishable, Fourier transforms it, we'll discard the noisy modes and we'll do the inverse Fourier transform to get a smooth scaled configuration.
 Finally, we'll study the differences between the universal configuration and each of the original scaled ones. 





  * Preparation of the scaled fields

Let v(x,y;ξ) and σ(x,y) the averaged measured field and its standard deviation respectively. We assume that the v-field depends on an random noise field that prevents it to have a smooth value. 
Then, by assuming that the central limit theorem applies, we can write

    v(x,y;ξ)=v(x,y)+σ(x,y)ξ(x,y)

where v is the field we would like to know and each ξ(x,y) is a Gaussian random variable with zero average and unit variance. Let us mention that one should expect a spatially correlated noise because fluctuations are based on several dynamically conserved quantities like density, momenta and energy, for instance. However, the assumption of independence maybe enough to catch the key ingredients of the errors and the corrections to the measured averaged magnitudes. 

We know that the fields v(x,y) and σ(x,y) have a natural symmetry: v(x,y)=λ v(1-x,y) and σ(x,y)=σ(1-x,y) with λ=± 1. That's because the physics we observe in the system: the magnitudes behave in a mirror like way with respect the axis x=1/2 even in the convective state when two symmetric rolls appear. Therefore, the observables depending on  the u_1 component of the hydrodynamic velocity field should have λ=-1 and any other observable λ=1. Obviously our measured configurations do not strictly respect this symmetry due to the fluctuations. We symmetrize our configurations and their errors to impose such important property: v(x,y;ξ)→ (v(x,y;ξ)+λ v(1-x,y;ξ))/2 and σ(x,y)→ (σ(x,y;ξ)^2+σ(1-x,y;ξ)^2)^1/2/2. With this averaging we also  improve the statistics and therefore the observable error bars.


We can get the averaged field value and its second momenta:

    v̅(ξ)   =   1/N_C∑_(x,y)v(x,y;ξ)=v̅+1/N_C∑_(x,y)σ(x,y)ξ(x,y)
    
    m_2(v;ξ)   =   1/N_C∑_(x,y)(v(x,y;ξ)-v̅(ξ))^2=m_2(v)+2/N_C∑_(x,y)(v(x,y)-v̅)σ(x,y)ξ(x,y)+1/N_C∑_(x,y)σ(x,y)^2ξ(x,y)^2
       +   1/N_C^2∑_(x,y)∑_(x',y')σ(x,y)σ(x',y')ξ(x,y)ξ(x',y')

where 

    v̅=1/N_C∑_(x,y)v(x,y)  ,   m_2(v)=1/N_C∑_(x,y)(v(x,y)-v̅)^2


we follow now the same strategy we explained in the introduction about the error analysis. We define the random variable

    B=m_2(v;ξ)-(1-1/N_C)1/N_C∑_(x,y)σ(x,y)^2

we see that

    ⟨ B⟩   =   m_2(v)
    ⟨ m_2(B)⟩   =   4/N_C^2∑_(x,y)(v(x,y)-v̅)^2σ(x,y)^2

where ⟨·⟩ is the average over the random fields values ξ. Then we can made the strong asumption that the central limit theorem applies to B in the sense that the square root of its second moment give us  the error on the computation of its average. Then

    B=m_2(v)± 3⟨ m_2(B)⟩^1/2⇒ m_2(v)=m_2(v;ξ)-(1-1/N_C)1/N_C∑_(x,y)σ(x,y)^2 ± 3⟨ m_2(B)⟩^1/2

In the computation we have assumed that σ(x,y) are small and we have disregard orders σ^4 in m_2(B) because it is used only on the error bar definition. There, the higher order corrections are negligible because our 3σ criteria is generous enough and it gives us already an idea of the precision level in our computations. 

We are ready to derive the corrections to the scaled configuration that we define as

    v^(s)(x,y;ξ)=v(x,y;ξ)-v̅(ξ)/σ(v;ξ)

where σ(v;ξ)^2=m_2(v;ξ). Then

    v^(s)(x,y)=v(x,y)-v̅/σ(v)=v^(s)(x,y;ξ)[ 1+1/2(1-1/N_C)β(0)-3/2N_Cβ(2)+σ(x,y)^2/N_Cσ(v)^2]-1/N_Cβ(1)± 3σ(x,y)/σ(v)

with

    β(n)=1/N_C∑_(x,y)v^(s)(x,y;ξ)^nσ(x,y)^2/σ(v)^2






We apply the above analysis to the components of the hydrodynamic velocity field. We show in figures <ref> and <ref> how they  average value, u̅_1,2(ξ), and standard deviation, σ(u_1,2;ξ), behave. u̅_1(ξ) is strictly zero because of the symmetrization procedure.  It is interesting to see the systematic queue that appear on σ(u_1,2;ξ) for small T_0 values. It is only due to the local equilibrium fluctuations of the hydrodynamic velocity. If we assume that below the transition the average hydrodynamic velocity at any cell is zero we only see the fluctuations due to the law of large numbers. That is, the modulus of each component is given by u_1,2(x,y)≃√(T(x,y))/√(N_D(x,y)) being T(x,y) the local temperature and N_D(x,y) the total number of data used in the time averaging, both at the cell (x,y). In our case N_D(x,y)=N t ρ(x,y)Δ^2/ρ where N=957 is the number of disks in the simulation, t=100000 is the total amount of measurements, ρ(x,y) is the averaged areal density at the cell, Δ=1/30 is the cell side and ρ is the system areal density. Then

    σ(u_1,2)_le^2=1/N_C∑_(x,y)u_1,2^2(x,y)=ρ/Δ^2NtN_C∑_(x,y)T(x,y)/ρ(x,y)

We plot this result as continuum lines in the center of figures <ref> and <ref> after symmetrizing the data. Observe how the line follows almost perfectly the measured averages. Moreover, we also plot the corrections we have computed above due to the fluctuating character of the data (small dots) and they follow the line and the big dots. That is

    σ(u_1,2)_le^2=(1-1/N_C)1/N_C∑_(x,y)σ(x,y)^2

in the non-convecting region.
We see again in this example the importance to subtract such noise contribution  (even though it is well understood and it  should disappear in the limit for infinite larger averaging or number of disks) because it introduces a systematic deviation in our analysis obscuring, sometimes, the hydrodynamic bulk behavior of the system. 
Finally, the figures in <ref> and <ref>  right show the value of σ(u_1,2) after subtracting the noise term contribution. We couldn't find a simple scaling behavior with g. Nevertheless, we found a striking simple relation:

    σ(u_2)=2σ(u_1)

for every T_0 value of any given g. See for instance in figure <ref> how the points overlap all over the T_0 range for the three g values.
At this moment, we do not know why a parameter that measures the average variation of the field should follow such relation. 

Once we have prepared with care the scaled fields we may start its analysis to see it they are similar or there is a small but systematic dependence on T_0.



  * Analysis of the scaled configuration spatial structure and mutual comparison:

We have measured some structural field parameters. First we have computed the Inertial Tensor with respect the field center of mass. That is,

    I_α,β(u_1,2^(s))=1/N_C∑_(x,y)[r'(x,y;u_1,2^(s))^2δ_α,β-r'_α(x,y;u_1,2^(s))r'_β(x,y;u_1,2^(s))]

where

    r'(x,y;u_1,2^(s))   =   r(x,y;u_1,2^(s))-R(u_1,2^(s))
    
    r(x,y;u_1,2^(s))   =   (x,y,u_1,2^(s)(x,y))
    
    R(u_1,2^(s))   =   1/N_C∑_(x,y)r(x,y;u_1,2^(s))

The Inertial Tensor has the following form:

    I(u_1,2^(s))=
    [ 13/12     0   I_1;     0 13/12   I_2;   I_1   I_2   1/6 ]

where

    I_1=-1/N_C∑_(x,y)(x-1/2)u_1,2^(s)(x,y)  ,   I_2=-1/N_C∑_(x,y)(y-1/2)u_1,2^(s)(x,y)

Observe that I_2=0 for u_1^(s) and I_1=0 for u_2^(s). From it, we can compute its principal axis and the semi principal diameters that are just the eigenvectors and one over the square root of the eigenvalues. We are mainly interested on the set eigenvalues because they contain an averaged information of the topological form of the field. In fact they inverse square root are the main axis of the Poinsot's ellipsoid. The eigenvalues can be analytically computed and are:

    λ_1=13/12  ,  λ_± = 1/2[5/4±√(121/144+4(I_1^2+I_2^2))]

We compute I_1 and I_2 using u_1,2(x,y;ξ) and therefore the eigenvalues are noise-dependent: λ(ξ). Following the same technique above explained we can get the noise-corrected eigenvalues that they are given by

    λ_±=λ_± (ξ)∓4β_2 (I_1^2+I_2^2)/√(121/144+4(I_1^2+I_2^2))

where we show here only the N_C→∞ correction. Nevertheless we use the full corrected form in the results we show in the figures. We observe in the figure <ref> (left) how λ's have a nontrivial dependence on T_0 and the noise term correction doesn't change too much λ(ξ). By other hand for the field u_2^(2) the values of I_2 are so small that the eigenvalues follow the value corresponding to I_2=0, that is, λ_1=13/12=1.08333.... Finally we observe that the scale of variation is of 10^-3 which is of the order of the hydrodynamic velocity resolution. That is, the Inertial Tensor is not giving us a clear answer to our question about the existence of a universal scaled field over a region or, at least, an asymptotic universal one.  




Another global magnitude that capture some of topological structure of a scaled field, v^(s)(x,y), is its n-th momenta:
 
    m_n(v^(s);ξ)=1/N_C∑_(x,y)v^(s)(x,y;ξ)^n

Following the same steps as above we can obtain its noise-corrected form and its error bar:

    m_n(v^(s))=m_n(v^(s);ξ)-Λ_n± 3 ϵ_n

where

    Λ_n   =   -m_n(v^(s);ξ)n/2(1-1/N_C)β(0)+m_n(v^(s);ξ)n(n+2)/2N_Cβ(2)+m_n-1(v^(s);ξ)n^2/N_Cβ(1)
       +   m_n-2(v^(s);ξ)n(n-1)/2N_Cβ(0)+n(n-1)/2(1-2/N_C)β(n-2)-n^2/N_Cβ(n)
    ϵ_n^2   =   n^2/N_C[m_n(v^(s);ξ)^2β(2)-2m_n(v^(s);ξ)β(n)+2m_n(v^(s);ξ)m_n-1(v^(s);ξ)β(1)
       +   β(2n-2)-2m_n-1(v^(s);ξ)β(n-1)+
    m_n-1(v^(s);ξ)^2β(0)
    ]

where β(n) is defined above. We show in figure <ref> the behavior of m_4. We observe that the data is noisy but it is compatible with  a constant behavior o to the existence of an asymptotic one for T_0>13 in all cases. From the u_2^(s) data we see that for the g=15 case the constant regime appears for T_0>14. Moreover, m_4(u_1^(s)) shows that  g=5 and g=10 values have similar behavior while g=15 seems to tend to the same value. However, m_4(u_2^(s)) indicates a clear g-dependence at least between g=5 and the other g-values. We tried to analyze the m_6 cases but the associated error bars obscure any useful analysis.

Finally we have measured the euclidean distance between any two given configurations v^(s)(x,y;T_0,g;ξ) with a fixed g value but at temperatures T_0 and T_0':

    D(v^(s);T_0,T_0',g,ξ)=[1/N_C∑_(x,y)(v^(s)(x,y;T_0,g;ξ)-v^(s)(x,y;T_0',g;ξ))^2]^1/2
 
We plot the set of pair distances for the scaled fields u_1^(s)(x,y;ξ) and u_2^(s)(x,y;ξ) for a given g values in figures <ref> and <ref> (pink surfaces in top figures). We observe that all of surfaces tend to an apparently constant value clearly different from zero. Nevertheless, for any given T_0>13 it seems that the distance is constant (non-zero) for all other T_0'>T_0. Again, the effect is mainly due to the data noise. We compute the noise correction and we find:

    D(v^(s);T_0,T_0',g)^2=D(v^(s);T_0,T_0',g,ξ)^2-Λ± 3Σ

where

    Λ   =   2/N_C∑_(x,y)(v^(s)(x,y;T_0,g;ξ)-v^(s)(x,y;T_0',g;ξ))(α_2(x,y;T_0,g)-α_2(x,y;T_0';g))
       +   1/N_C∑_(x,y)(α_1(x,y;T_0,g)+α_1(x,y;T_0';g))
    Σ^2   =   4/N_C^2∑_(x,y)(v^(s)(x,y;T_0,g;ξ)-v^(s)(x,y;T_0',g;ξ))^2(α_1(x,y;T_0,g)+α_1(x,y;T_0';g))

and

    α_1(x,y;T_0,g)   =   (1-2/N_C)σ(x,y)^2/σ^2+1/N_Cβ(2)v^(s)(x,y)^2+β(0)/N_C-2/N_Cv^(s)(x,y)^2σ(x,y)^2/σ^2+2β(1)/N_Cv^(s)(x,y)
    α_2(x,y;T_0,g)   =   1/2(3/N_Cβ(2)-(1-1/N_C)β(0))v^(s)(x,y)-1/N_Cv^(s)(x,y)σ(x,y)^2/σ^2+β(1)/N_C

where β(n), σ(x,y), σ, ..., are computed for a given configuration with T_0 and g values following the above definitions.
We expect D to be zero or near to zero. Let us define D_±(v^(s);T_0,T_0',g;ξ)=D(v^(s);T_0,T_0',g;ξ)±√(Λ). We expect that 
D_-(v^(s);T_0,T_0',g,ξ)=D_-(v^(s);T_0,T_0',g)± 3Σ_- and D_+(v^(s);T_0,T_0',g,ξ)=D_+(v^(s);T_0,T_0',g)± 3Σ_+. Then, substituting into the D^2 expression we find that the non fluctuating parts follow:  D^2=D_-D_+ and the error parts: Σ_-D_++Σ_+D_-≃Σ. Assuming that D_- is going to be almost zero we can approximate Σ_-≃Σ/D_+. Finally we can write

    D_-(v^(s);T_0,T_0',g)=D(v^(s);T_0,T_0',g;ξ)-√(Λ)± 3Σ/(D(v^(s);T_0,T_0',g;ξ)+√(Λ))

We show in figures <ref> and <ref> the marginal distance D_- for u_1^(s) and u_2^(s) respectively. We see how the bare distance D(v^(s);T_0,T_0',g;ξ) (pink surfaces) have some dependence on T_0 and their limiting values are far from zero. Once we introduce the noise correction √(Λ) the marginal distance D_- tend to zero for T_0>13 which is coherent with the results obtained from other topological parameters. Moreover, we have an estimate of the error bars that we show at the bottom of the figures. Observe the quality of the convergence to zero of D_-for g=10 and 15. g=5 is much noisier but still it is consistent with the zero marginal distance value. 

















  * Obtaining the universal fields:

It is clear now  that the rescaled fields seems to be equal (except for fluctuations) when T_0>13. To obtain the universal scaled field ũ_1,2(x,y) we follow two steps: (i) we average (for a given g-value) all the noise-corrected scaled fields u_1,2^(s)(x,y) with T_0>13 and (ii) we Fourier Transform the averaged field, we disregard modes with fluctuating small amplitude and we inverse Fourier Transform the remaining modes. All of these just to get a smooth field without fluctuating parts. 

In figures <ref> and <ref> we show the modulus of the Fourier Transform for u_1 and u_2 fields respectively. We observe how the dominant modes have ≃ 10  intensity meanwhile the discarded modes have intensities less than 0.4. In the table <ref> we show the cut-off values used in each case. In general we choose such values as the minimum one in which we have a compact and continuum set of modes.  In Figure <ref> we show the final smoothed universal fields ũ_1(x,y) and ũ_2(x,y) for g=5, 10 and 15. One observes at a glance that the universal fields for different g-values are very similar. When subtracting we see systematic small differences that have spatial structure. However, such differences are small and of the order of the error bars. At this point and with our set of simulations we cannot conclude about the existence of a g-independent universal profile.

The final check is to see if the difference between any T_0>13 scaled configuration and the universal field gives place to a random fluctuating field with amplitudes smaller than the error bars of the scaled field. The error for a point (x,y) of any given rescaled configuration have the form of 3σ(x,y)/σ assuming that the central limit theorem holds. Figure <ref> shows, as an example, the difference between u_1,2^(s)(x,y;T_0=17) and the corresponding ũ_1,2(x,y). We observe how practically all the fluctuating data is between the error bar interval. Moreover, there is a spatial dependence that is similar for the differences and the error bars. From the point of view of the data errors we cannot do better than this. However, we can try to see if there is any systematic spatial regularity in the differences. 
We define the difference field:

    w_d(x,y;T_0,ξ)=σ(T_0)u_1,2^(s)(x,y;T_0;ξ)-ũ_1,2(x,y)/σ(x,y;T_0)

In this form, if u_1,2^(s) follows exactly the noise ansatz: u_1,2^(s)(x,y)=ũ_1,2(x,y)+σ(x,y)ξ(x,y)/σ, then w_d(x,y;T_0;ξ)=ξ(x,y), that is, we should see a pure white noise. In particular we defined the average of its absolute value:

    Dif=1/N_C∑_(x,y)| w_d(x,y;T_0,ξ) |

In case that w_d was a pure white noise the value of Diff would be 2/√(2π)≃ 0.798...
We see in figure <ref> this magnitude as a function of T_0 (≥ 13). We observe how this magnitude goes around the white noise value. The small deviations observed could be due to the cutoff used in which we could discard some weak structural behavior in modes that are superimposed to the noise and/or some small system correlations. Nevertheless, the white noise average behavior seem to us quite remarkable because it implies that our analysis is capturing the principal parts of the system behavior. Finally, we show in figure <ref> the Fourier Transform of w_d for T_0=17 and g=5, 10 and 15 just to see if there is any clear spatial structure. We do not see any clear regularity that would mean the existence of a systematic deviation between the universal field and the scaled one. 
 All of this support the idea that our universal field is a good average for each T_0>13 velocity field.


We may conclude, after this extent statistical analysis, that the hydrodynamic velocity field components have an universal scaling property, that is, we can write them as:

    u_1,2(x,y)= σ(u_1,2(T_0,g))ũ_1,2(x,y;g)

That is, the T_0 dependence can be explicitely extracted from the fields. As we commented above, we cannot exclude the possibility of having ũ_1,2 independent on g. Another questions that remains also unsolved is: if we could improve data statistics  for the T_0<14 configurations in such a way that their errors were smaller than the typical profile variation, would we still find the same scaling behavior? Is then, the scaling behavior typical of all convective states? or  is it just an asymptotic like property? Those are open questions that are beyond the computational effort presented in this paper.




  * Temperature field: 

We define local temperature as the kinetic energy measured with respect the local center of mass reference frame. 
In other words, the local temperature at cell (n,l) is:

    T(n,l)=K(n,l) -1/2u(n,l)^2

where K(n,l) is the averaged one particle kinetic energy at box (n,l):

    K(n,l)=1/2N(n,l)M∑_t=1^M∑_i:r_i(t)∈ B(n,l)v_i(t)^2   ,   N(n,l)=1/M∑_t=1^M∑_i:r_i(t)∈ B(n,l)1

and u(n,l) is the modulus of the hydrodynamic velocity studied in the previous section. We see that the temperature definition is composed of two terms that they are independently measured. We consider that the error is the sum of the errors:

    ϵ(T)=ϵ(K)+| u_1|ϵ(u_1)+| u_2|ϵ(u_2)








We show in figure <ref> two typical cases: one without convection (T_0=2 and g=10) and the other with convection (T_0=18 and g=10). Observe that both figures present some general features: (1) The profile is nonlinear on the vertical direction y, (2) there is no spatial structure on the x direction when we are at a non-convective state and there is a small but systematic x non-constant profile when the system is in a convective state and (3) there is a  thermal gap on the boundaries, that is, the extrapolated profile for y=0 and y=1 doesn't match the thermal bath temperatures.

In order to characterize the temperature profile we first study the x-averaged profiles:

    τ(y)=1/30∑_n=1^30 T(n,l)   ,   y=(l-1)/30+1/60



We see in figure <ref> the averaged temperature profiles τ. We manage to fit all the data points with a very simple fourth order polynomial that permits us to refine the analysis. First, we observe for g=10 and T_0>12 that the profile has a small change of curvature at the middle of the profile. The same happens for g=15 and T_0>15. In order to characterize these changes, we subtract to the data the tangent to the profile at y=11/60 (which is the coordinate of the 6th cell) which is computed using the fitted curve. If  the difference is positive means that the slope grows (is negative) and the curve is concave. In contrast, if the difference is negative the curve is convex.
We plot such differences on figures <ref>. We can see that for g=0 and 5 and any value of T_0 the profile is convex. Moreover, the convexity grows with T_0. In contrast, for  g=10  and T_0>12 the τ(y) profile has a changing convexity: it is almost flat near the hot thermal bath and it is concave for a while and again convex near the cold thermal bath. Also, this changing behavior is reinforced for large values of T_0. The same happens for g=15 with a concave behavior stronger than g=10 but appearing for T_0>15. From the phenomenological fit we can get the inflection point by doing the second derivative and looking the zeros of the resulting polynomial. In figure <ref> it is shown the inflection points x^* as a function of T_0 for g=10 and 15. They grow with T_0 but all of them are around x≃ 0.5, far from the boundaries. All of this indicate a real strong nonlinear behavior of the y-temperature profile. The increment of g emphasizes the nonlinear character of τ(y). 




Finally we have studied the behavior of the thermal gaps at the boundaries. We see that the profiles do not extend to the thermal bath values, that is, τ(y=0)≠ T_0 and τ(y=1)≠ 1. This is a well known phenomena called thermal resistance that appear when we there are currents crossing a singular boundary that, in our case is due to the use of thermal baths that act to the particles when they hit the walls. Typically it should disappear at the thermodynamic limit. We define the thermal gap as

    G=|τ(y=0,1)-T_0,1|

In figure <ref> we show how the gaps G behave with T_0 for g=0, 5, 10 and 15.  The thermal gap for the hot reservoir grows monotonically with T_0. For large values of T_0 it seems to growth linearly. In this region, we have done fits to the data of the form G=a_0T_0+a_1+a_2/T_0+a_3/T_0^2+… and we got that  G≃ 0.17(0.02)T_0, 0.17(0.01)T_0, 0.17(0.01)T_0 and 0.18(0.02)T_0 for g=0, 5, 10 and 15 respectively. We cannot discard a g-independent behavior of G≃ 0.17T_0. By other hand, the thermal gap behavior for the cold reservoir has a richer structure (maybe due to the fits used and/or the extrapolated data goodness). For g=0 and 5 it seems to monotonically grow up to a limiting value which is computed by using a fit of the form G=a_1+a_2/T_0+a_3/T_0^2+… to the T_0>9 points: a_1=1.10(0.07) and a_1=1.13(0.06) for g=0 and 5 respectively. For g=10 and 15 the gap grows fast up to the critical temperature T_c, then stays more or less constant and it smoothly grows again up to a limiting value of a_1=1.6(0.1) and 3.0(0.1) for g=10 and 15 respectively. We cannot clearly relate the constant region with the T_c,2 critical temperature defined above. 

However, independently of the different cold and hot gap behaviors, the effective thermal gradient, Δτ=τ(y=0)-τ(y=1), present a very smooth behavior as a function of the imposed external gradient, Δ T=T_0-1. We see in figure <ref> its behavior with T_0 and for the three g-values. Except for a small region near the equilibrium, the effective gradient is proportional to Δ T: Δτ≃ bΔ T, with b=0.86(0.05), 0.81(0.09), 0.88(0.06) and 0.79(0.06) for g=0, 5, 10 and 15 respectively. We cannot discard the existence of a common value for all g's. It looks like the system re adapts itself to such effective gradient once we fix the values of the external thermal baths. This could be a characteristic of hard disk systems where the properties of the system are invariant under a global rescale of the bath temperatures. That's why we can fix T_1=1 without any loss of generality.

Once we have study the x-averaged behavior of the temperature field, we think it is natural to subtract this x-average, τ(y), to the field in order to study the spatial structure of the deviations that appear in the convective state. Again, we think that this field may also have the scaling behavior that we discovered in the hydrodynamic velocity field. Therefore we define the scaled excess of temperature, T^(s)(x,y), as:

    T^(s)(x,y)=T(x,y)-τ(y)/σ(T)

where

    σ(T)^2=1/N_C∑_(x,y)[T(x,y)-τ(y)]^2   ,  τ(y)=1/L_C∑_x T(x,y)

where L_C=30 in our case.



We have plotted in figure <ref> the superposition of the scaled excess of temperature fields for T_0=17, 18, 19 and 20 and different g values. We observe that, again, it seems that there are common surfaces for each g value. It is also clear that there are the fields depend on the g-value. between g values. We can follow then the same items we did for the hydrodynamic velocity field. Let us go through them.



  * Preparation of scaled excess field:
Let T(x,y;ξ) and σ(x,y) the averaged measured field abd its standard deviation respectively where we already have symmetrized the data with respect the x=1/2 axis. As before we assume that we get a value that is a particular realization of an inherent Gaussian white noise and therefore:

    T(x,y;ξ)=T(x,y)+σ(x,y)ξ(x,y)

where ξ(x,y) is an uncorrelated Gaussian random variable with zero mean and unit variance. The y-profile is then given by:

    T(y;ξ)=T(y)+1/L_C∑_x σ(x,y)ξ(x,y)  ,   T(y)=1/L_C∑_x T(x,y)

and the excess of the v-field is defined by:

    v(x,y;ξ)=T(x,y;ξ)-T(y;ξ)=v(x,y)+σ(x,y)ξ(x,y)-1/L_C∑_x σ(x,y)ξ(x,y)

where v(x,y)=T(x,y)-T(y). Observe that there is a difference with respect the hydrodynamic velocity field case: we subtracted there the overall field average and here just the y-averaged profile. This makes the computation rather different: for instance ∑_x v(x,y;ξ)=0 and ∑_x v(x,y)=0. The second momenta is now given by:

    m_2(v;ξ)=1/N_C∑_(x,y)v(x,y;ξ)^2=m_2(v)+A_1(v;ξ)+A_2(v;ξ)  ,   m_2(v)=1/N_C∑_(x,y)v(x,y)^2

where

    A_1   =   2/N_C∑_(x,y)v(x,y)σ(x,y)ξ(x,y)
    
    A_2   =   1/N_C∑_(x,y)σ(x,y)^2ξ(x,y)^2-1/L_C∑_y[1/L_C∑_x'σ(x',y)ξ(x',y)]^2

following the techniques already explained we find that

    m_2(v)=m_2(v;ξ)-⟨ A_2(ξ)⟩± 3 √(⟨ A_1(ξ)^2⟩)

where ⟨·⟩ the the average over the uncorrelated white noise field. We show in figure <ref> the behaviors of m_2(v;ξ) and m_2(v).
First, we observe how the m_2(v;ξ) have a queue for small T_0 values that it is due to the noisy data behavior. The computed ⟨ A_2(ξ)⟩ term is sensibly larger that the measured value. This effect could be due to our assumption that the σ(x,y) we have obtained is not the intensity of a pure white noise, that is, there are correlations effects not taked into account in the error theory we have developed. Another related possibility is that we have build  σ(x,y) as the straight superposition of the K, u_1 and u_2 errors assuming that they were mutually independent. We know that there are some correlations between the fluctuations in the total kinetic energy of a cell and the fluctuations of the average values of the cell velocity. That is, we are overestimating the error in T by definition.  With the set of data we have obtained we cannot make a precise and independent computation of the T's fluctuation to study its behavior. In any case we observe that multiplying the σ(x,y) by a factor f=0.87 we manage to fit the observed queue with the computed A_2 correction in all cases. We are going to use this effective factor in the error definition from now on. For instance, in figure <ref> (right) we show m_2(v) as a function of T_0 and  we already included this error correction. In this case, observe that the values of m_2(v) are very small for T_0<13 in contrast with the behavior of the hydrodynamic velocity field. Moreover, we see there that m_2(v) grows and it doesn't seem to reach any finite limit. 

Finally, the measured scaled excess field is given by

    T^(s)(x,y;ξ)=T(x,y;ξ)-T(y;ξ)/√(m_2(T(x,y;ξ)-T(y;ξ)))=v(x,y;ξ)/√(m_2(v;ξ))
 

It can be expanded with respect of σ(x,y) which is assumed to be small:

    T^(s)(x,y;ξ)=T^(s)(x,y)+α_1(x,y;ξ)+α_2(x,y;ξ)+O(σ^3)

where

    α_1(x,y;ξ)   =   σ̃(x,y)ξ(x,y)-1/L_C∑_x'σ̃(x',y)ξ(x',y)-1/2A̅_1 T^(s)(x,y)
    α_2(x,y;ξ)   =   -1/2A̅_1σ̃(x,y)ξ(x,y)+1/2A̅_1 1/L_C∑_x'σ̃(x',y)ξ(x',y)-1/2A̅_2 T^(s)(x,y)+3/8A̅_1^2 T^(s)(x,y)

where σ̃(x,y)=σ(x,y)/σ(v), A̅_i=A_i/σ(v)^2 (i=1,2) and σ(v)=m_2(v)^1/2.



  * Analysis of the scaled configuration:
The Inertial Tensor of the scaled field is diagonal by construction and it doesn't depend on any T^(s) property. By other hand we have studied the 4'th momenta of it: 

    m_4(T^(s);ξ)=1/N_C∑_(x,y)T^(s)(x,y;ξ)^4

Its corrected noise form as a function of the observed one is given by:

    m_4(T^(s))=m_4(T^(s);ξ)-Λ± 3 ϵ

with

    Λ   =   4/N_C∑_(x,y)T^(s)(x,y)^3⟨α_2(x,y;ξ)⟩+6/N_C∑_(x,y)T^(s)(x,y)^2⟨α_1^2⟩
    ϵ   =   16/N_C^2∑_(x,y)∑_(x',y')T^(s)(x,y)^3T^(s)(x',y')^3⟨α_1(x,y;ξ)α_1(x',y';ξ)⟩




In figure <ref> we observe again how the m_4'th moment tends for each g to a constant value independent on T_0. The g dependence of such constants is small.

We have also studied the averaged euclidean distance between any two configurations T^(s)(x,y;T_0,g;ξ) with a fixed g value and different T_0 temperatures. As in the case of the hydrodynamic velocity field, we can write:

    D(T^(s);T_0,T_0',g)^2=D_-(T^(s);T_0,T_0',g)D_+(T^(s);T_0,T_0',g)

where

    D_-(T^(s);T_0,T_0',g)   =   D(T^(s);T_0,T_0',g;ξ)- Λ± 3ϵ(D_-)
    
    D_+(T^(s);T_0,T_0',g)   =   D(T^(s);T_0,T_0',g;ξ)+ Λ± 3ϵ(D_+)


    Λ^2   =   1/N_C∑_(x,y)[⟨α_1(x,y;ξ;T_0)^2⟩+ ⟨α_1(x,y;ξ;T_0')^2⟩
       +   2(T^(s)(x,y;T_0)-T^(s)(x,y;T_0'))(⟨α_2(x,y;ξ;T_0)⟩-⟨α_2(x,y;ξ;T_0')⟩)]

and

    ϵ(D_-)=√(⟨ a(ξ)^2⟩)/| D(T^(s);T_0,T_0',g;ξ)- Λ|+D(T^(s);T_0,T_0',g;ξ)+ Λ

with

    a(ξ)=2/N_C∑_(x,y)(T^(s)(x,y;T_0)- T^(s)(x,y;T_0'))(α_1(x,y;ξ;T_0)-α_1(x,y;ξ;T_0'))


We plot in figure <ref> the mutual marginal distance, D_-, between configurations. We observe how the raw distance tend to a non-zero constant that, once we include the noise corrections, is zero for T_0>11 in all cases. Let us remark that we have included the correcting factor in the T's errors. Without including it, the corrected distances would be constant but negative. That is consistent with our observation that we overestimated the error in T.




  * The universal fields: At this point we follow the same path we did with the hydrodynamic velocity field. First, we averaged the scaled excess temperature fields with T_0=14, … 20 just to have better statistics. Then we do a discrete Fourier Transform of the resulting field (see figure <ref>) discarding the noisy wave numbers (we show in table <ref> the limiting intensity below which the modes are neglected) and, finally, we do the inverse Fourier Transform obtaining the universal scaled excess of temperature field (see figure <ref>). 


 Finally we have also compared the universal scaled field with the one with T_0=17. We again see in figure <ref> how such differences are constrained to the 3σ(x,y) error bars of the measured field. Morever, we again see the average of the absolute value of w_d-like field (as we defined it in the hydrodynamic velocity section) follows remarkably well the white noise value. (see figure <ref>). The Fourier Transform of the w_d field is shown in figure <ref>. We observe in this case a very weak structure around the (0,0) mode. It could be just our effective error definition that it does not capture correctly the real spatial correlations of the noise. 















We can conclude that there exists an universal field T̃(x,y) in such a way that

    T(x,y;T_0,g)= T(y;T_0,g)+σ(T(x,y)-T(y);T_0,g))T̃(x,y;g)

for the configurations with T_0∈[14,20]. 



  * Density profiles: 




The local density is defined by eq. (<ref>). In figure <ref> we present two typical measured density  fields for g=10: T_0=2 (left figures) and T_0=18 (right figures).  For T_0=2 we are at a non-convective state and we observe that a barometric type of density profile is present: due to the effect of the driven external field the particles tend to accumulate at the bottom of the system. We also see that (as in the temperature case) ρ(x,y) does not present structure on the x-direction. However, for T_0=18 the density profile is inverted and now the low density region is at the bottom of the cage and the high density is at the top of it. Moreover, there is a clear spatial structure on the x-direction. We can see in both cases the boundary effects we already commented and that it results in a decrease of their value near the walls.

It seems reasonable to try to characterize the density marginal ω(y) defined by the average of its x values:

    ω(y)=1/N_A∑_x∈ Aρ(x,y)

where A is the set of cells we use for the averaging and N_A is the total number of cells. We discard the two first columns and rows near the boundaries. In our analysis we choose A={(2n-1)/60 | n=3,…,28}  and N_A=26.


We observe some interesting properties (see figure <ref>): (1) For T_0=1 the density profile presents a typical barometric equilibrium density distribution we already studied.  When the external field is on the bottom density is larger than the top density. (2) The barometric-like profile tends to disappear as we increase the value of T_0. (3) There is an inversion of the density for large values of T_0: the large density region is at the top of the cage. (4) For g=0 the density is always a monotonous increasing function for T_0>1.
(5) We do not see relation between being at the non-convective or convective state and the profile type. In fact, for g=5 and T_0>4 the profiles are monotonous increasing function of y and for T_0<4 are decreasing functions of y, far from the appearance of the convective state which occurs at T_c≃ 2.1 in this case. Similar thing occurs for g=10 and g=15. However, there is a clear relation with T_c,2: w(y) is a monotonous decreasing function for T_0<T_c,2 and it looses such monotonicity for T_0>T_c,2 appearing first some increasing parts (positive slopes) for y≃ 1 when T_0 has values near T_c,2 and w being a monotonous increasing function for large enough values of T_0. Again it seems that the strong convecting state has the properties that a fluid in a convective state should have: low density fluid is at the bottom of the container and it tends to go up and the cold high density fluid is at the top and it tends to go down. The unclear issue is to fit this image for the region T_c<T_0<T_c,2 where we see convection and there is no density inversion.
(6) It is curious that for T_0<T_c,2 the best fitting function is the equilibrium like one with varying constants: y=a(0)+a(1)log w+a(2)log(1-w)+a(3)/(1-w)+a(4)/(1-w)^2. For T_0>T_c,2 the best fit is a fourth order polynomial.























Again we can study now the scaled excess density field: :

    ρ^(s)(x,y)=ρ(x,y)-w(y)/σ(ρ)

where

    σ(ρ)^2=1/N_C∑_(x,y)[ρ(x,y)-w(y)]^2   ,   w(y)=1/L_C∑_x ρ(x,y)

and N_C=L_C^2 and L_C=26  in our case and the sums exclude the two nearest boundary rows and columns.
At this point the analysis follows the same steps we did with the temperature field. We show in figures <ref> the superposition of the scaled density fields ρ^(s) for T_0=17, 18, 19 and 20 we see again that all of them seem to follow the same surface without any systematic deviation. In figures <ref> we present the behavior of σ(ρ)^2. We observe  that here there is not a systematic tail on σ(ρ)^2 for small values of T_0. Moreover it seems that the values are non-zero and small.  The correction we compute assuming a random white noise, ⟨ A_2⟩, (small dots with red color in the figure) is exactly the part to be corrected.
Finally, in figure <ref> we show how the fourth momenta of ρ^(s) tends to a constant for large values of T_0. We see in figure <ref> how the mutual marginal distance tend to zero once we subtract the noise effect. We average the  scaled configurations with T_0∈ [14,20] values for g=5 and g=10 and with T_0∈ [16,20]  for g=15 (because we detected in this case some small but systematic deviations when we included the ones with T_0=14 and 15). 
In order to get rid of the noise effects, we do the Fourier Transform  of the averaged configurations (figure <ref>) and we discard the modes below the cutoff showed in Table <ref>. After doing the inverse Fourier Transform to the remaining modes we get the universal scaled configuration for the excess of the density field ρ̃(x,y) (figure <ref>). Finally, the differences between these universal configurations and the scaled ones behave like  a ramdom white noise that is bounded by our measured error (figures <ref>, <ref> and <ref> and see comments in the temperature field section).

We conclude this analysis by remarking that we have observed that there exists  an universal field ρ̃(x,y) such that:

    ρ(x,y:T_0,g)=w(y;T_0,g)+σ(ρ(T_0,g)-w(T_0,g))ρ̃(x,y;g)

for T_0's large enough. 
The universal profiles seem to be equal. In fact we see that the difference between different g-values is, at most of 5% and there is a small but systematic spatial deviations on the differences. However we cannot discard from the numerics the possible existence of a common universal profile independently of the value of g. 

NOTE: We have discarded two rows and columns near the boundary to study the scaling of the excess density field. That changes the σ(ρ(T_0,g)-w(T_0,g)) with respect to the one we would get if we used all the points in the system and it maybe could change the universal character of the scaled density field.  Let us find the changes that we get when we change the domain to compute the scaled field.
Let us assume that in a given square domain D× D we have a field that can be written as

    f(x,y;T_0)=f(y;T_0)+σ(T_0)f̃(x,y)

where

    f(y;T_0)=1/| D|∫_Ddx f(x,y;T_0)  ,  σ(T_0)^2=1/| D|^2∫_Ddx∫_Ddy [f(x,y;T_0)-f(y;T_0)]^2

Let us now consider that we compute the scaled field in a different domain D'× D':

    f^(s)(x,y;T_0,D')=f(x,y;T_0)-f(y;T_0,D')/σ(T_0,D')
 
where now

    f(y;T_0,D')=1/| D'|∫_D'dx f(x,y;T_0)  ,  σ(T_0,D')^2=1/| D'|^2∫_D'dx∫_D'dy [f(x,y;T_0)-f(y;T_0,D')]^2

Is it f^(s)(x,y;T_0,D') universal?, that is, f^(s)(x,y;T_0,D')=f̃(x,y;D')? It is a matter of simple algebra to see that:

    f(y;T_0,D')   =   f(y;T_0)+σ(T_0)/| D'|∫_D'dxf̃(x,y)
    σ(T_0,D')   =   σ(T_0)Σ(D')

where

    Σ(D')^2=1/| D'|^2∫_D'dx∫_D'dy[f̃(x,y)-1/| D'|∫_D'dx'f̃(x',y)  ]^2

We can substitute these relations into the f^(s)(x,y;T_0,D') definition and we get

    f^(s)(x,y;T_0,D')=1/Σ(D')[f̃(x,y)-1/| D'|∫_D'dx'f̃(x',y) ]≡f̃(x,y;D')

that is, the scaled field in the new domain is still universal (it doesn't depend on T_0). However we pay the price of having a different scaled universal function because, in any case, it depends on the domain in which we define the scale.






  * Pressure profiles: 
We use the virial theorem to compute the pressure field P(x,y) even in the nonequlibrium stationary state:
 
    P(x,y)=ρ(x,y) e_c(x,y)/π r^2+lim_τ→∞1/2 Δ^2 τ∑_n:t_n∈[0,τ]r_ij· p_ij

where Δ=1/30 is the side length of a cell, e_c(x,y)=⟨ E_c(x,y)⟩/⟨ N(x,y)⟩ is the average total kinetic energy in the cell and the sum runs over all particle-particle collisions that occur at the cell (x,y) at the time interval [0,τ] assuming that at time zero the system is at the stationary state. From the derivation of this expression we can call P(x,y) to be the local mechanical pressure that, as we will see, it is non trivially  related with the thermodynamic pressure when there is a nonzero hydrodynamic velocity field. From a computational point of view we already checked this expression for the equilibrium case T_0=T_1=1 and we showed the goodness of our simulation when we compared the data obtained with the exact result (assuming the Henderson's EOS).  

At figure <ref> we present two typical measured pressure fields for g=10: T_0=2 (left figures) and T_0=18 (right figures).  For T_0=2 we are at a non-convective state and we observe that the pressure decreases with y and it has some structure on x when we are at the convective state (right figure in <ref>).  


Again we first study the marginal q(y) defined by the average of its x values:

    q(y)=1/N_A∑_x∈ AQ(x,y)

where A is the set of cells we use for the averaging and N_A is the total number of cells. We discard only the two rows near the thermal walls. However, we do not discard the columns near the vertical boundaries in order to compare with the measured pressure on the top and bottom boundaries (see section Global Magnitudes: the pressure). In our analysis we choose A={(2n-1)/60 | n=1,…,30}  and N_A=30. 

 In figure <ref>  we see: (1) the pressure is a decreasing monotone function on y,  (2) there is a smooth change of curvature of the pressure profiles with the temperature gradient from concave to convex at around T_0=T_c,2. However we cannot precise the transition temperature because the second derivatives of the polynomial fitted functions present some wavy behavior that obscure the convexity property, (3) we see in eq. <ref> that the effect of g into the equilibrium pressure profiles is through a trivial rescaling of the y coordinate by g: y̅=gy. However, for any nonequlibrium state this doesn't happen and the pressure profiles have a nontrivial dependence on g and (4) the naive extrapolation of the fitted 6th order polynomial to the extreme points y=0,1 (Q_v) coincide with the measured values of the pressure at the boundaries (Q_w). In fact, for the g=0 case the relative error is, at most, of about 0.15%. For the g=5, 10 and 15 cases the relative errors for the hot boundary (y=0) is about 1%. For the cold boundary (y=1) the relative errors fluctuate much more as we see in figure <ref>. That is due to the fitting function we are using that doesn't represent the fast decay that is present for small values of T_0. We already had some similar troubles in the study of the density field for such values. There we used the equilibrium form for the fitting function in order to get very good fits.   















We define the scaled excess pressure field:

    Q^(s)(x,y)=Q(x,y)-q(y)/σ(Q)

where

    σ(Q)^2=1/N_C∑_(x,y)[Q(x,y)-q(y)]^2   ,   q(y)=1/L_C∑_x Q(x,y)

and N_C=L_C^2 and L_C=26  in our case and the sums exclude the two nearest boundary rows and columns.

We show in figures <ref> the superposition of the scaled pressure fields Q^(s) for T_0=17, 18, 19 and 20. We see how the g=5 is so noisy that we cannot discriminate, at a glance, a regular surface. For g=10 and 15 we see again that all of them seem to follow the same surface without any systematic deviation. In figures <ref> we present the behavior of σ(Q)^2 and  the fourth momenta of Q^(s).
Due to computational problems we didn't get the errors in the local virial pressure and we are unable to get the corrections to the field due to the noise of the data. In any case we see that their behavior is very similar to the ones of the temperature or the density fields. In particular the fourth momenta tend to a constant value independent of T_0. The rest of the analysis follows the same path of the other already analyzed magnitudes.

 We see in figure <ref> how the mutual marginal distance tend to a constant. The constant is larger in the g=5 case due to the intrinsic fluctuations. We average the  scaled configurations with T_0∈ [14,20] values for all g values. 
In order to get rid of the noise effects, we do the Fourier Transform  of the averaged configurations (figure <ref>) and we discard the modes below the cutoff showed in Table <ref>. After doing the inverse Fourier Transform to the remaining modes we get the universal scaled configuration for the excess of the pressure field Q̃(x,y) (figure <ref>). Finally, the differences between these universal configurations are shown in figure <ref>.  .


We again observe that, for the configurations with T_0∈[14,20] exists an universal field Q̃(x,y) such that 

    Q(x,y;T_0,g)= q(y;T_0,g)+σ(Q(T_0,g)-q(T_0,g))Q̃(x,y;g)






§ CONNECTING WITH HYDRODYNAMICS


Once we have studied global and local magnitudes from the microscopic model, we would like to compare the observed system behavior with the one we would obtain from stationary Navier-Stokes equations applied to a hard disk system. 
 We already explained that NS equations are based on local conservation laws and in two relevant assumptions:


  * Local equilibrium hypothesis: It is assumed that thermodynamics holds locally. That is, we can define local thermodynamic quantities like density, temperature, pressure, ... and they are assumed to be related via the equilibrium thermodynamic relations. For instance by the equation of state.

  * The constitutive relations: The local energy current is assumed just to be proportional to the temperature gradient (Fourier law) and the pressure tensor has only linear dependence on the hydrodynamic velocity gradients (newtonian fluid). Moreover, there are transport coefficient functions like the thermal conductivity, the shear viscosity and the second viscosity that they are assumed to depend only on the local values of density and temperature.

Some questions may arise at this point: Up to what extend the computer simulation is going to give us results to be compared with hydrodynamics?  Is local equilibrium hypothesis valid always or only for (relatively) small external temperature gradients? and the Fourier law or the newtonian fluid hypothesis? Also there are other interesting questions to discuss, for instance What is the relation between the mechanical pressure and the thermodynamic pressure in the convective regime?, or
Is it correct to use the Boussinesq approximation to study the non-convective-convective transition? We'll intent to answer some of these questions. 



  * 0. The Problem of the numerical derivatives

We have measured in our computer simulation a set of observables  in  each of the 30× 30 virtual cells in which  we have divided our system. The obtained averaged values are, in most cases, smooth enough to look like a spatial regular field. In order to compare our results with the NS-equations (without solving them) we need to do the derivatives to such computed discrete fields. There are many strategies to do derivatives of a discrete set of variables, for instance discrete local schemes or data polynomial interpolations before doing analytical derivatives. However, if we want to get the goodness or consistency of our results compared with the NS equations, we need to control the propagation of the errors on the derivatives. They  are coming from two sources: the intrinsic error bars in each measured datum and the chosen derivative scheme. After playing with several possibilities we decided use the following strategy. Let (a(n,l),σ(n,l)) the computed value of the observable and its statistical standard deviation at cell (n,l) then: 


  * (a) We create a new data set by perturbing a(n,l) by a gaussian random number with zero mean and σ(n,l) standard deviation: ã(n,l)=a(n,l)+σ(n,l)ξ(n,l) with ξ(n,l)∈ N(0,1).

  * (b)We fit the data to a seventh order polynomial on each variable. That is,  the fitted function to the data, A(x,y), is of the form:

    A(x,y)=∑_n=1^N_x∑_l=1^N_yc(n,l)(x-1/2)^n(y-1/2)^l

where typically N_x=N_y=7.
We can include boundary terms if necessary. For instance, if we know that the observable should be zero at x=0,1 we multiply the polynomial by x(1-x).  

  * (c) We do an analytic derivative to the fitted function A(x,y). 

  * (d) We compute the values of the derivative and its square at the (x,y) points corresponding to the center of the virtual cells. 

  * (e) We iterate steps (a)-(d) 1000 times and average at each cell point the derivative and its square. 

  * (f) the local derivative is the average value and its dispersion is the standard deviation of the thousand data values obtained. 

We show in figures <ref> and <ref> the results of this strategy by plotting the partial derivatives of the density, ρ, pressure, Q, temperature T and the hydrodynamic velocity components u_1 and u_2 for g=15 and T_0=20. Observe how the derivatives present a larger error interval on the boundaries where there are no constraints which it is due to the variability of the fitted function for values on the boundary (remind that we take the center of the box, x=1/2 and y=1/2, as the seed of the polynomial expansion). It is clear that the error is reasonably small  and it gives us some hope into obtain some reasonable info about the NS-equations behavior.




  * 1. The continuity equation

The continuity equation express the hydrodynamic mass conservation:

    ∂_1(ρ u_1)+∂_2(ρ u_2)=0

We may compute each of the derivatives of our discrete fields ρ u_1 and ρ u_2 obtained in the computer simulation and check that its sum gives zero. Obviously, from our polynomial fitting method (5 or 7th order polynomial) we cannot expect a perfect matching (that would only occur if we used tha exact analytic solution). However we should see that the answer is contained inside the computed error bars. In order to convince the reader that the continuity equation holds let us use two complementary points of view. First we show in figure <ref> the numerically obtained ∂_1(ρ u_1)  and ∂_2(ρ u_2) fields, and its corresponding sum for the case with T_0=20 and g=5, 10 and 15. We observe that in all cases the sum of each component is compatible with the zero value. The relative error is smaller as we increase the g value. For g=15, the variability of the derivatives is of order 0.4 and the error of the sum is bounded between ± 0.03, that is a 0.06 interval of variation. That is, we are seeing a true cancellation when summing up the derivatives of two different functions more beyond the error noise.

The second strategy is to plot ∂_1(ρ u_1) vs ∂_2(ρ u_2) for each virtual cell point where the derivative is computed. If continuity equation holds the data should follow the line y=-x for any value of T_0 and g. We show in figure <ref> the data for T_0=14, 17 and 20 for g=5, 10 and 15. We observe how, beyond error bars, the data follows the expected line. We can conclude that the numerical analysis permits us to show that the basic hydrodynamic continuity equation is correct This is the minimum requirement we had to pass to continue, with some confidence, the analysis of the NS hydrodynamic equations from our computer simulation.



 


  * 2. Checking the Local Equilibrium hypothesis

In order to close the conservative equations in Navier-Stokes we need to relate the local temperatures, pressures and densities assuming that equilibrium  applies at each macroscopic spatial point. This is a well assumed property that is called local equilibrium hypothesis. Let us emphasize that it doesn't mean that the fluid is at equilibrium locally. That is, the local probability measures is not a Gibbs measure because there are currents of energy, momentum and matter that breaks the symmetry properties of the Hamiltonian. Local equilibrium just means that the equation of state of equilibrium applies at each macroscopic point. This could be so, probably because the nonequilbrium measure distribution shape near the most probable configurations is similar to the corresponding Gibbs distribution but the tails deviate from it in a highly nontrivial way.

We know that the equilibrium equation of state (EOS) for hard disks is unknown. However, there are a set of analytic expressions that approximate the measured results with relative errors smaller than 1% (for instance the Henderson EOS). These expressions may help us to check, up to some extend, if local equilibrium holds for the hard disk system. In any case, we have an important advantage in the hard disk case. We know from the equilibrium Gibbs probability measure that the EOS has the structure:
 
    Q=T F(ρ)=Tρ H(ρ)

because the interaction potential has only two values: zero (when disks move freely) or infinity (when disks are in contact).
If local equilibrium applies,  the values of density, temperature and  pressure measured at each cell should follow the same relation given by eq.(<ref>) and it should be independent on the external gradient T_0 or the value of the external driving field g. This is a very exigent assumption. In figure <ref> we plot the 75712 points of (ρ(x,y),T(x,y),Q(x,y)) that we have obtained in the set of all simulations shown in this paper. We have discarded the points near the boundaries (the nearest two rows and columns) to eliminate the exclusion surface effects due the existence of the rigid boundaries. We observe in the figures the variety of profiles in the (ρ,T,Q) space and how all follow (apparently) the same surface in such space. 



We can be more precise. If we plot (ρ(x,y),Q(x,y)/ρ(x,y)T(x,y)), the EOS for hard disks, eq.(<ref>), predicts that all points should follow exactly the same function: H(ρ). 
In figure <ref> we show the   75712 data points in such graph. The dotted line is the function proposed by Henderson that is assumed to be correct (up to 1% of relative error):

    H_H(ρ)=1+ρ^2/8/(1-ρ)^2







We observe how the complete set of data follows the same curve that coincides apparently with the Henderson's proposal. However, we appreciated two small but systematic deviations: 


  * (1) The points corresponding to the non-convective regime deviate, in average, from the Henderson's EOS: ⟨Δ H⟩=∑_i=1^N_D(H_i-H_H(ρ_i))/N_D=-0.00098 where i runs over all data cells pertaining to configurations at non-convective states, the total number of such data is N_D=18928, H_i=π r^2 P_i/ρ_iT_i, and P_i, ρ_i and T_i are the measured pressure, density and temperature at each cell i. The data mean standard deviation of such differences is σ(Δ H)=0.00402297... and the statistical error of the average is 3σ(Δ H)/√(N_D)=0.0000877 which is almost 10 times smaller that the observed average. That is, we are seeing some deviation from the Henderson's EOS that is small but systematic. Moreover, if we normalize the differences: ξ=(Δ H-⟨Δ H⟩)/σ(Δ H) we obtain a non-Gaussian behavior (see figure <ref>). This implies that H_H(ρ) is not the real EOS (as we knew already). It is surprising to us how our data detects such fact in a statistic way. Nevertheless the data set is in the 1% relative error taking as reference the Henderson's EOS as it was expected (see figure <ref> left).

  * (b) Some points in the range ρ∈[0.2,0.3] deviated (more than an usual fluctuation) from the rest. After some analysis we found that part of the data corresponding to the convective regime (see figure <ref> right) are the one that deviates systematically from the scaling behavior. In order to study such observed deviations we show in figure <ref>  the relative differences between the data points with the Henderson's proposal. 
We observe how, for all the data pertaining to states in non-convecting states, the points are distributed around zero and they do not cross the limiting errors lines of 1%. However, the data from convective states are distributed in a non symmetric manner around zero and there is a set of points deviating far from the limiting error lines. We think that this is not casual and there should be some reason beyond fluctuating data and/or errors. That is, it seems that the existence of a nonzero hydrodynamic velocity field affects to the existence of a unique H(ρ) function. In fact, the observed deviation grows systematically with g which is correlated with the fact that the modulus of the observed velocity field also grows with g. 


 

In order to understand the apparent violation of the local equilibrium from the data corresponding  to systems in the convective regime, we should remind some definitions.
In Navier-Stokes equations, the pressure appears as a part of the stress tensor:

    τ_ij=-Pδ_ij+η(∂_iu_j+∂_ju_i)+η'δ_ij∑_k∂_ku_k

where P is the Pressure that appears in the EOS (also called Thermodynamic Pressure). One can also define  the mean normal stress:

    P̅=-1/dTr(τ)

where d is the system dimension and it is written is this way to guarantee its invariance with respect any change of coordinates. Let us recall that the stress tensor is defined from the vector small  forcing (δ f) acting over a small surface (nδ S):

    δ f_i=∑_jτ_ijn_jδ S

For a fluid at rest (u=0) we find that δ f=-P nδ S. That is, the mean normal stress coincides with the Thermodynamic Pressure which is the force applied by unit surface and perpendicular to it. When the fluid is moving, P̅ and P differ and they are related by:

    P̅=P-(η'+2/dη)∇· u

where η'+2/dη=ξ is known as the bulk viscosity.
P̅ is often called Mechanical Pressure because it measures the total averaging forcing at a given fluid point. Let us mention that typically it is assumed that the bulk viscosity is a function of the local variables ρ(x) and T(x). However there is also the possibility that ξ being a function of P̅, T and ρ implying that P̅ could be a nonlinear function of ∇· u (see for instance <cit.>).

Observe that when the fluid is at rest or the fluid is incompressible (∇· u=0) both magnitudes coincide P̅=P. Stokes in 1845 assumed that, in general, the fluids have ξ=0 (Stokes assumption). This assumption is still under discussion but it seems that the actual theoretical and experimental evidence is against it (see for instance <cit.>).

At this point we should think about what we have really measured by using the virial pressure expression (<ref>). Looking with care its derivation we can conclude that the virial pressure fields P(x,y) analyzed above are in fact the Mechanical Pressure fields P̅(x,y) because we are measuring the TOTAL external forcing applied to a given region. It this is so, the observed deviations to the expected EOS from the data pertaining to convecting states should be due to the use in the EOS of P̅ instead to the correct P fields. Therefore, we would like to get the thermodynamic pressure from our data by using eq. (<ref>) and afterwards to check if the systematic observed deviations disappear. Obviously the first question to be solved is to see how P̅ depends on the local ρ, T and ∇· u. In particular we would like to check if it is just a linear function on ∇· u.

Let us first write eq. (<ref>) for hard disk systems. We assume the general form:

    P̅(x,y)=ρ(x,y) T(x,y)H(ρ(x,y))/π r^2-ξ(ρ(x,y),T(x,y),P̅(x,y))∇· u

where H(ρ)=π r^2 P/ρ T (see eq.(<ref>)). Dimensional arguments imply that the bulk viscosity should be of the form

    ξ(ρ,T,P̅)=√(T)/rE̅_B(ρ,π r^2P̅/T)

where the arguments of the function E̅_B are dimensionless (see for instance the expression in Enskog's approximation). Substituting into eq.(<ref>) we find:

    H̅(ρ)=H(ρ)- Ẽ_B(ρ,H̅) ∇· u/√(T)

where H̅=π r^2P̅/ρ T and Ẽ_B=π rE̅_B/ρ. That is, H̅ is solution of an implicit equation and it has the general form of: 

    H̅=G(ρ, ∇· u/√(T))

 We know that for u=0, H̅=H and then G(ρ,0)= H(ρ). Then it is convenient to write:

    H̅(ρ,∇· u/√(T))=H(ρ)+G_1(ρ,∇· u/√(T))

where we have assumed that G(x,y) is an analytic function. 

Observe that if G_1(x,y)=g(x)y we are in the case where the bulk viscosity only depends on the local temperature and local density (ξ=ξ(ρ,T)=√(T)E̅_B(ρ)/r) and g(ρ)=-π r E̅_B(ρ)/ρ. In a general case, G_1(x,y)= g_1(x)y+O(y^2) when y→ 0. We expect also that for very low density the system tends to the ideal gas behavior and then H̅≃ 1+O(ρ). But H≃ 1+O(ρ) and then G_1 should have the form: G_1(x,y)=xG̅_1(x,y) with G̅_1 being an analytic function on both arguments.

Then, by studying eq.(<ref>) using our data for H̅, ρ, ∇· u and T, we can try to fit H and G_1 functions. However we deal with several numerical problems: (a) We see that H̅ does not differ too much from H and the data obtained for H̅ fluctuates around H, therefore it is going to be very difficult to obtain a clear picture of the G_1 function and (b) the numerical derivatives contain some rough structure and errors that causes a noisy result that obscures the system behavior. Therefore, we can probably discard or confirm tendencies from a statistical analysis of a set of structured cloud of points.



We already know that, in general, numerical derivatives contain an extra noise due to the necessity of doing some fit and/or interpolation to the discrete data field. When computing the divergence of the velocity field we have two ways of getting it: (1) obtaining the derivatives  ∂_1u_1 and ∂_2u_2 (using the method above explained) and summing up both fields or (2) by assuming that the continuity equation is correct we obtain from it the relation ∇· u=-u·∇ρ/ρ, then we just need to compute the derivatives of the ρ field. Each method depend on the numerical computation of different field derivatives. We show ∇· u in figure <ref> computed by the two methods for T_0=17 and g=15. Observe that the direct method (1)  give us a wavy field that seems to follow, in average, the one we obtain by the method (2) which is much smoother and its error bars are smaller.  That is due to the less fluctuating behavior of the density field compared with the hydrodynamic velocity field that also affects to the error propagation when doing the derivatives. It is clear that the second method is the one we should use in the EOS analysis.



Once we solved the problem of how we measure the ∇· u, we can continue the analysis of the mechanical pressure. First we plot in figure left in (<ref>) the measured local pressure H̅=Q̅/ρ T=π r^2P̅/ρ T as a function of local (ρ,∇· u/√(T)) for all the cells (discarding  two rows and columns nearest to the four boundaries) from stationary configurations we have obtained in this numerical experiment (T_0∈[1,20] and g=0,5,10 and 15). Let us point out that we did the full analysis separately for the different g values and we did not see any clear dependence on g, thats why we join all data sets. There are 75712 data points in total. Observe in the figure how the data from non-convective configurations  accumulate into a curve with ∇· u=0. The rest of the data extends over the space following (apparently) a smooth surface.

We have to fit a regular surface of the form given by eq.(<ref>) to all these points without any theoretical help to build a test function. All we know is: (a) H(ρ) is the real EOS for hard disks and that the Henderson's EOS is a good approximation to the real one but there are other well known EOS approximations having good behaviors on different density ranges (see for instance the Chapter 3 in <cit.>) and (b) the function G_1 is assumed to be analytic and it is of the form G_1(x,y)= xyG_2(x,y) due to physics constraints (see above).

Therefore, after many trials, we have chosen to fit functions of the form:

    H̅(x,y)=1/(1-x)^2[1+∑_n=1^M_1c_2(n)x^n+1]+∑_n=1^M_2∑_l=1^M_3 c_1(n+M_2(l-1))x^ny^l
 
where M=M_1+M_2*M_3 is the total number of parameters we use to fit the surface. Let's point out that many  coefficients that measure the goodness of a fit improve with the total number of parameters. Therefore if we want to compare different approaches we need to fix M to do not introduce spurious improvement in the fit.
Let us show explicitly the analysis over three particular structures with M=8 and M=10 as typical examples that confirms our final answer. We have tested many other M's that give us the same final results. These trial functions are for M=8: A(M_1=3,M_2=1,M_3=5) that permits a nonlinear behavior in ∇· u, B(8,0,0) that assumes that the Stoke's assumption holds and tries to fit only the EOS, and C(3,5,1) that assumes that the bulk viscosity depends only on the local values of ρ and T. For M=10 the corresponding cases are A(4,2,3), B(10,0,0) and C(4,6,1).

In table <ref> we show the obtained parameter values with their standard deviations.


We should develope a criteria for a given M value to elect which one, A, B or C, is the best fit. We focus on the distribution of the differences between the data and each of the fitted surfaces:
Δ H_i=H̅_i-H̅(ρ_i,∇· u_i/√(T_i) for i=1,…, N_D=75712. 
 In particular we have computed the statistics of the differences that includes the distribution of values and the first four central moments:
 
    ⟨Δ H⟩   =   1/N_D∑_i=1^N_DΔ H_i 
    σ(Δ H)^2   =   1/N_D∑_i=1^N_D(Δ H_i-⟨Δ H⟩)^2
    ϵ(Δ H)   =   3σ(Δ H)/√(N_D)
    
     s(Δ H)   =   1/N_Dσ(Δ H)^3∑_i=1^N_D(Δ H_i-⟨Δ H⟩)^3
    κ(Δ H)   =   1/N_Dσ(Δ H)^4∑_i=1^N_D(Δ H_i-⟨Δ H⟩)^4-3

 ϵ(Δ H) is a measure of the statistical error of ⟨Δ H⟩.
 Moreover, we have also studied the corresponding spatial point distribution by computing first the center of mass:
  
    R=1/N_D∑_i=1^N_Dr_i     r_i ≡ (Δ H_i,ρ_i,∇· u_i/√(T_i))

and then obtaining the inertial tensor with respect the center of mass: 

    I_α,β=1/N_D∑_i=1^N_D[r'_i^2δ_α,β-r'_i,αr'_i,β]

where r'_i=r_i-R. From the inertial tensor we can get the 
principal axes of inertia that are defined by its eigenvectors which give us some averaged idea of the spatial point distribution. 


In order to interpret all these parameters let us think what would happened to the data if the fitted function H̅ was the exact one. If we neglect the mutual correlations from the data (they come from fields whose values at each cell are not independent from the other cells), then we assume that the points Δ H_i should be distributed, in this case, evenly around the plane (x,y,0). We do not expect a total homogeneous distribution in space because the sampling is, by construction, not homogeneous (it depends on temperature gradients, external fields and the resulting fields have non-linear and non-uniform values). Moreover, we also assume that the probability distribution of Δ H should be a Gaussian and/or a near Gaussian distribution with large kurtosis that would reflect the correlation of the  data.  In fact, it is simple to show that the moments of the distribution ξ=(Δ H-⟨Δ H⟩)/σ(Δ H) are proportional to the moments of the differences between the fitted function and the correct one (assuming that the data fluctuates with a Gaussian distribution around the correct values of the function). That is, we would get s=0 and κ=0 if the fitted function was exactly the correct one (except for errors due to the finite size of the sampling).

That is the picture we have in mind. Therefore, a given fit would be better than other if: (a) ⟨Δ H⟩=R_3 (the third coordinate of the center of mass) is nearest to zero, (b) s(Δ H) and κ(Δ H) have the smaller possible value and (c) one of the principal axes of inertia is nearest to the vector (0,0,1) because it means that the points are more evenly distributed in space around the plane (x,y,0).
A final element  is that we can compute the error of the center of mass third coordinate R_3=⟨Δ H⟩: ϵ(Δ H)=3σ(Δ H)/√(N_D). The fit is good if the error bars around the value of R_3 contains the zero value: 0∈[R_3-ϵ(Δ H),R_3+ϵ(Δ H)]. If that doesn't happen the fit is considered wrong.

In figure <ref> we show, for the case M=8, the point spatial distributions with the center of mass and the principal axes of inertia and the probability distribution of ξ=(Δ H-⟨Δ H⟩)/σ(Δ H) that it is compared with a Gaussian N(0,1).
We show in Table <ref> the values obtained for the above parameters for the three reference fits A,B and C and M=8 and 10.  We highlighted the fit with the best value.






 There are some comments to make the best fit's final choice:


  (a) The fit B (Stokes assumption) is the worse one always. None of the parameters in table <ref> does better than the fits A or C for any M value. We can conclude that the Stokes assumption is not followed by our system. 

  (b) The fit C (Mechanical pressure proportional to ∇· u) is,in general, worse than the fit A. For instance, for the M=8 case we observe that the R_3 value is so large that the zero value is barely included in the error interval of ϵ(Δ H) and for the M=10 case the zero is out of the error interval meaning that the fit is useless.  Moreover the skewness of the data with fit C is always larger  the one's with A. C does better in the kurtosis value for M=8 but no with M=10 (but it just reflects the existence of data correlation that we cannot evaluate its real magnitude). C also does better that A in the  the third component of the principal inertial axis for M=8  but this is maybe due to a better description by the fitted function of the density dependence of G_1. This is confirmed when we use M=10 where we see that  v_3 for A is larger than the corresponding one with the C fit.  

  (c) We have observed that the overall behavior of fits type A, B or C is conserved when we do other fits with different total number of parameters. 

We can conclude that, from a statistical point of view, the fit type A (for most of the M values checked)  is  the most consistent with the data. That is, from our analysis we can clearly discard the Stokes assumption in our system and that the mechanical pressure is a nonlinear function on ∇· u.  Therefore, we are quite convinced  that 
the bulk viscosity is a function of the local density, ρ, local temperature T and the local mechanical pressure P̅.



In figure <ref> we show the row data (ρ_i,H̅_i) (first row) and the corrected mechanical pressure using the G_1 function obtained with the A,B and C-fits, (ρ_i,H̅_i-G_1(ρ_i,∇· u_i/√(T_i)) in the second. third and fourth rows respectively. We appreciate with some detail at the center figures at each row how the data is corrected and they scale in a more homogeneous way for the fit A compared with the original set of data. The correction is worse for fits B and C.



We show in figure <ref> the A fit overlapped to the data. The fitted function can be written:

    H̅=H_A(ρ)-π r^2/ρ Tξ(ρ,∇· u/√(T))∇· u

with

    H_A(ρ)=1-0.0726 ρ^2 +1.6265ρ^3 -4.9088 ρ^4+4.9040ρ^5/(1 - ρ)^2

that is the A-fit EOS and

    ξ(ρ,∇· u/√(T))=√(T)/rρ^2[91.08-386.40ρ+78.45ỹ-225.82ρỹ+18.38ỹ^2-70.02ρỹ^2]

where  ỹ=∇· u /(π r√(T)) anr π r= 0.02562.. in our simulations. We can compare with the result at the Enskog approximation. There the leading term is 49.72 that is almost half of our result.  We do not pretend that these fitted functions have the correct mathematical structure. Nevertheless, there are reasonable smooth fits and they are compatible with the already known approximations (for instance with the hard disk EOS) and, more important, there are necessary to rescale the data to obtain a better distribution of them around the fitted functions.  After all these analysis we may conclude that: (1) we detect the existence of the  mechanical pressure concept and its different behavior at the convective regime with respect the thermodynamic pressure, (2) Assuming that the average stress tensor is the mechanical pressure that we have measured, we conclude that the Stokes relation cannot be correct (in this case), (3) the bulk viscosity seems to be a function of the local density, temperature and the mechanical pressure and then the mechanical pressure cannot be a linear function of the divergence of the hydrodynamic velocity vector field.


  * 3. Checking the stationary Navier-Stokes equations for non-convective states

We already know that the stationary Navier Stokes equations for non-convective states are:

    dQ/dy   =   -gρ
    √(T)K(ρ)dT/dy   =   -J

where it is assumed the functions only depend on the y-coordinate. We should add the local equilibrium hypothesis to these equations and the temperature boundary conditions. Our aim is to check that the data is compatible with these two equations. Let us study first the equation for the pressure.

We can integrate the pressure equation:

    Q(y)-Q(y_0)=-g∫_y_0^y dy̅ρ(y̅)

If we want to get a precise check, it is convenient to take into account the differences between the measured variables in a cell and the above continuum expression.  First, we can define the pressure value over a one-dimensional cell, Q̅(n):

    Q̅(n)=1/N_C∑_m=3^N_C-2 Q(m,n)

where Q(m,n) is the pressure measured in the cell (m,n) and we have discarded the columns near the vertical boundaries.
Assuming the existence of the continuum functions Q(y) and ρ(y) we can write

    Q̅(n)=1/Δ∫_ỹ(n)-Δ/2^ỹ(n)+Δ/2 dy  Q(y)

where ỹ(n)=(n-1/2)Δ. 
Substituting equation (<ref>) into eq. (<ref>) we get:

    Q̅(n)=Q(y_0)-g/Δ∫_ỹ(n)-Δ/2^ỹ(n)+Δ/2 dy ∫_y_0^ydy̅ρ(y̅)


The integration variables of the double integral can be exchanged with the corresponding change of the integral limits. In this way we can do one of the integrals.The result, when y_0<ỹ(n)-Δ/2, is:

    Q̅(n)=Q(y_0)-g∫_y_0^ỹ(n)-Δ/2dyρ(y)-g/Δ∫_ỹ(m,n)-Δ/2^ỹ(m,n)+Δ/2 dyρ(y)(ỹ(n)+Δ/2-y)
 
We now choose y_0=ỹ(n_0)-Δ/2=(n_0-1)Δ (the bottom of the cell n_0) and we can write

    Q̅(n)=Q(y_0)-gΔ∑_l=n_0^n-1ρ̅(l)-g/Δ∫_(n-1)Δ^nΔdyρ(y)(nΔ-y)

where ρ̅(l) is the average density at a cell at height l. The ρ(y) inside the integral can be Taylor expanded around the ỹ(n) coordinate and we integrate order by order in the expansion. Afterwards, we obtain:

    Q̅(n)=Q(y_0)-gΔ∑_l=n_0^n-1ρ̅(l)-g/Δ∑_k=0^∞(-1)^kΔ^k+1/k! 2^k+2(1-(-1)^k/k+2+1+(-1)^k/k+1)d^kρ/dy^k|_y=ỹ(n)

To find Q(y_0) as a function of Q̃(n_0) it is enough to substitute n=n_0 in the last expression. After its substitution we get:

    Q̅(n)-Q̅(n_0)   =   -gΔ∑_l=n_0^n-1ρ̅(l)
       -   g/Δ∑_k=0^∞(-1)^kΔ^k+1/k! 2^k+2(1-(-1)^k/k+2+1+(-1)^k/k+1)[d^kρ/dy^k|_y=ỹ(n)-d^kρ/dy^k|_y=ỹ(l_0)]

Finally, the derivatives at a central point can be expressed as variations with cell functions ρ̅ as it was explained in the introduction. Then we get the cell version of the barometric equation:

    Q̅(n)-Q̅(n_0)   =   -gΔ∑_l=n_0^n-1ρ̅(l)-gΔ/2(ρ̅(n)-ρ̅(n_0))
       +   gΔ/24[ρ̅(n+1)-ρ̅(n-1)-ρ̅(n_0+1)+ρ̅(n_0-1)]+O(Δ^3)≡ -D(n,n_0)

We use this last expression to check the barometric equation. We measure by one hand
Δ Q(n,n_0)≡Q̅(n_0)-Q̅(n) and by the other hand D(n,n_0) defined in eq. (<ref>). We choose the values n_0=4,…,N_C-4 and n=n_0+1,…,N_C-3 to not to use the data from the boundary cells.


We study the barometric equation in three temperature regions for all the g's we have simulated: (I) T<T_c(g) (non-convecting states), (II) T_c(g)<T<T_c,2(g) (weak convecting states) and (III) T>T_c(g) (fully convecting states). Obviously we expect only good agreement for data in region (I).In figure <ref> we show  the results obtained.
We see in the first row how the 8004 data of region (I) behave. Δ Q(n,l) follows a clean straight line for all distance differences and g's studied. We do a linear fit to the data 
Δ Q(n,l)_F=a_0+a_1 D(n,l) whose results are in table <ref>. In order to see small deviations we study the remainders: Δ Q(n,l)-Δ Q(n,l)_F and their statistical distribution and central moments. We observe that the data follows a deformed Gaussian distribution with a high kurtosis value. That may indicate that the data is somehow correlated. That's could be because the data differences are build over all pairs of a given set of pressures and densities. In region (II) we cannot see any obvious deviation to the barometric equation.
Finally, we see in region (III) how the data presents an obvious structure indicating the they do not follow the barometric equation. The data structure grows with the g value in accordance with the strong convection behavior.
It is somehow curious how in regions (II) and (III) the linear fit still is very good. That is, in average the data seems to follow the barometric equation with some small dispersion even though there are convecting structures on the system. That is consitent with the fact that the overall difference between the pressures on top and bottom of the system always follow the barometric relation (see above).

Finally, we can check if the Fourier's Law  is consistent with our data:

    √(T)K(ρ)dT/dy=-J

in order to do that we should show that exists a function K(ρ) that is independent on g, T_0 and ρ̅) and  has the form:

    K(ρ̅)=-J√(T(y))[dT/dy]^-1|_y=ρ^-1(ρ̅)

In order to get K we should compute the local spatial derivatives of the T(y) profiles at a given cell, its local density and the heat current passing through the system. The main problem is to get a precise value of the derivative of a function build with only 30 points. The measured profiles are quite smooth and with relative small error bars then the appropriate strategy is to fit a function to the data and doing its analytic derivative. There are several problems however: (1)  Any function we use to interpolate the data to the data minimize the distance to the data points but, typically, they do not cross through the points. That is, even though we had a set of discrete data points obtained from a given function, the fitted one will cross the real function by above and below it. Therefore, the derivatives at a given point have a systematic and rather large dispersion due to such undulations. This effect increases if, like in our case, the data have error bars, and (2) for small gradients and/or large g values, the temperature profiles are very flat and the heat currents very small. That is, the ratio current/derivative is in this case very sensitive to the errors in the data and in the fit. Typically, we will discard these profiles.   

We observe that the temperature profiles have a well defined convexity: d^2T/dy^2<0. Therefore we have used a fitting function with such property:

    p(y)=-∑_n=1^s∑_l=1^s a(n)a(l)/(n+l+2)(n+l+1)y^n+l+2+a(s+1)y+a(s+2)

with a(n), n=1,… , s+2 being the free parameters to be fitted.

In figure <ref> we show the fits for the values of T_0 being at the non-convective regime for g=0, 5, 10 and 15. The fits are reasonable good with regression coefficients near one. Observe that some data for low gradients present the commented smaller values of R due to the flat profile problem we already commented.
In figure <ref> we show the values of K obtained by doing the derivative of T(y) of the fitted function at the point y for states at the non-convecting regime (left figure) and convecting one (right figure) for most of the gradients and all the g's we have simulated. We do not show the cases having large dispersions due to having near flat profiles and/or fits with low values of the regression coefficient.  We observe that the data scale with still some dispersion. Moreover we also see the already observed deviation from Enskog result. 

Finally we checked the possibility that the y profile for T_0 values at the convective regime, might follow the same Fourier's Law and/or some effective one. In figure <ref> (right)  we show that there is no trace of any scaling. 







  * 4. Some comments about the boundary conditions

In the computer simulation the boundary conditions choosed (reflecting vertical walls at x=0,1 and thermal walls at y=0,1) condition the particle behavior and the fields measured at the stationary state. In particular we observe:

    u_1(0,y)=0    u_1(1,y)=0   u_2(x,0)=0    u_2(x,1)=0

We also showed above that the difference of pressures measured at the top and bottom of the system follow the barometric formula with high precision (see figure <ref>):

    P(1)-P(0)=-gρ

where ρ is the system average density. We can use this information to know the kind of boundary conditions we should use on the Navier-Stokes equations to obtain this property. We know that the momentum equations in NS can be writen: 

    ∂_1P   =   ∂_1(τ_11'-ρ u_1^2)+∂_2(τ_12'-ρ u_1u_2)
    ∂_2P   =   ∂_1(τ_12'-ρ u_1u_2)+∂_2(τ_22'-ρ u_2^2)-ρ g

We can integrate both equations on x and y over all the domain. The left hand sides are:

    ∫_0^1dx∫_0^1dy∂_1 P(x,y)   =   ∫_0^1dy[P(1,y)-P(0,y)]=0
    ∫_0^1dx∫_0^1dy∂_2 P(x,y)   =   ∫_0^1dx[P(x,1)-P(x,0)]≡ P(1)-P(0)
 
The first average is zero by the boundary conditions symmetry. The complete set of pressure equations are, after integration,

    0   =   ∫_0^1dy[τ_11'(1,y)-τ_11'(0,y)]+∫_0^1 dx[τ_12'(x,1)-τ_12'(x,0) ]
    
    P(1)-P(0)   =   -gρ+ ∫_0^1dy[τ_12'(1,y)-τ_12'(0,y)]+∫_0^1 dx[τ_22'(x,1)-τ_22'(x,0) ]

Therefore, the conditions to get the observed barometric formula are:

    ∫_0^1dy[τ_11'(1,y)-τ_11'(0,y)]+∫_0^1 dx[τ_12'(x,1)-τ_12'(x,0) ]   =   0
    ∫_0^1dy[τ_12'(1,y)-τ_12'(0,y)]+∫_0^1 dx[τ_22'(x,1)-τ_22'(x,0) ]   =   0

In hydrodynamics  one can use many different boundary conditions for the hydrodynamic velocities. For instance, the no-slip boundary conditions consider u(x,y)=(0,0) ∀ (x,y)∈∂Λ that is incompatible with the results we have obtained in our simulation. On the other hand, the stress-free boundary conditions are 

    u(x,y)· n=0     ,  τ_12'(x,y)=0      ∀ (x,y)∈∂Λ

where n is a unitary perpendicular vector to the boundary. Observe that this boundary conditions are compatible with our observations if

    ∫_0^1dyτ_11'(1,y)   =   ∫_0^1dyτ_11'(0,y)
    ∫_0^1 dxτ_22'(x,1)   =   ∫_0^1 dxτ_22'(x,0)

In conclusion, our computer simulation is compatible with the well known in hydrodynamics slip boundary conditions if we add the properties given by equations (<ref>).






§ APPENDIX I: LYAPUNOV CENTRAL LIMIT THEOREM


Let X_1, X_2, …, X_N independent random variables such that 

    ⟨ X_i⟩=X̅_i   ,   ⟨ (X_i-X̅_i)^2⟩=σ_i < ∞

If for any δ >0 the following condition holds (Lyapunov condition):

    lim_N→∞1/s_N^2+δ∑_i=1^N⟨| X_i-X̅_i|^2+δ⟩=0

where

    s_N^2=∑_i=1^Nσ_i^2

Then 

    lim_N→∞1/s_N∑_i=1^N(X_i-X̅_i) d→  N(0,1)

where d means a convergence on distribution and N(0,1) is a gaussian distribution with 0 average and σ=1. We can formally rewrite the later property:

    1/N∑_i=1^NX_i N→∞≃1/N∑_i=1^NX̅_i+ξ/N√(∑_i=1^Nσ_i^2)

where ξ is a gaussian random variable with N(0,1) distribution.

As an example  let us assume that the X_i random variables are gaussian distributed with N(X̅_i,σ_i) where
σ_i∈ [σ_0,σ_1]. In this case the Lyapunov condition holds. Let us show it.

First we compute  the average:

    ⟨| X_i-X̅_i|^2+δ⟩=1/√(2π)σ_i∫_-∞^∞dx | x-X̅_i|^2+δexp[-(x-X̅_i)^2/2σ_i^2]=C(δ)σ_i^2+δ

where C(δ)=Γ((3+δ)/2)2^(3+δ)/2/√(2π) is a positive constant.
 Then, the Lyapunov condition can be written:
 
    ∑_i=1^N σ_i^2+δ/(∑_i=1^N σ_i^2)^1+δ/2N→∞→ 0

 For δ=2 we know that
  
    ∑_i=1^N σ_i^4/(∑_i=1^N σ_i^2)^2 <1/1+(N-1)σ_0^4/σ_1^4

  and it goes to zero when N→∞ whenever σ_0>0 and σ_1<∞.
  


§ ACKNOWLEDGEMENTS


We want to acknowledge the support from  National Science Foundation [grant DMR1104500] and AFOSR [grant F49620-01-0154] and FIS2017-84256-P Spanish Ministry MINECO. 



99
Alder62 Alder, B.J. and Wainwright, T.E., Phase Transitions in Elastic Disks, Physical Review, 127,  359-361 (1962).
Batchelor Batchelor, G.K., An introduction to fluid dynamics, Cambridge University Press (2000).
Boden Bodenschatz, E., Pesch, W. and Ahlers, G., Recent developments in Rayleigh-Bénard convection, Annu. Rev. Fluid Mech. 32, 709-778 (2000).
Bormann Bormann, Andreas S., The onset of convection in the Rayleigh-Bénard problem for compressible fluids, Continuum Mechanics and Thermodynamics, 13, 9–23 (2001). Numerical Linear Stability Analysis for Compressible Fluids in Analysis and Numerics for Conservation Laws, Warnecke, Gerald (ed.), pp. 93–105,  Springer Berlin Heidelberg (2005).
Buresti Buresti, G., A note on Stokes' hypothesis, Acta Mechanica, 226, 3555-3559 (2015).
Chandra Chandrasekhar, S., Hydrodynamic and Hydromagnetic stability, Dover, New York (1981).
Cordero95 Cordero, P., Marín, M. and Risso, D., Efficient Simulations of Microscopic Fluids: Algorithm and Experiments., Chaos, Solitons & Fractals,  6,  95-104, (1995). 
delPozo del Pozo, J., Garrido, P.L. and Hurtado P.I., Probing local equilibrium in nonequilibrium fluids, Physical Review E, 92, 022117 (2015); Scaling laws and bulk-boundary decoupling in heat flow, Physical Review E, 91, 032116 (2015). 
delPozo2 del Pozo, J., Garrido, P.L. and Hurtado P.I. in preparation.
Engel Engel, Michael and Anderson, Joshua A. and Glotzer, Sharon C. and Isobe, Masaharu and Bernard, Etienne P. and Krauth, Werner, Hard-disk equation of state: First-order liquid-hexatic transition in two dimensions with three simulation methods, Phys. Rev E, 87, 042134, (2013).
Esposito Esposito, R., Lebowitz, J.L. and Marra, R., On the derivation of hydrodynamics from the Boltzmann equation, Physics of Fluids, 11, 2354-2366 (1999).
Gad Gad-el-Hak, M., Stokes's hypothesis for a Newtonian, Isotropic Fluid,  Journal of Fluids Engineering 117, 3-5 (1995).
Gallavotti Gallavotti, G., Foundations of Fluid Dynamics, Springer, Berlin (2002).
Gass71 Gass, D.M., Enskog Theory for a Rigid Disk Fluid, Journal of Chemical Physics 54, 1898-1902 (1971).
Hen75 Henderson, D., A Simple Equation of State for Hard Discs, Mol. Phys. 30, 971 (1975).
Lappa Lappa, M. Thermal Convection: Patterns, Evolution and Stability, Wiley 2010. 
Manela Manela, A. and Frankel, I., On the Rayleigh-Bénard problem: dominant compressibility effects, J. Fluid Mech. 565, 461-475 (2006).
Mar87 Mareschal, M. and Kestemont, E., Order and Fluctuations in Nonequilibrium Molecular Dynamics Simulations of Two-Dimensional Fluids, Journal of Statistical Physics 48, 1187-1201 (1987).
Mar88 Mareschal, M., Malek Mansour, M., Puhl, A. and Kestermont, E., Molecular Dynamics versus Hydrodynamics in a two dymensional Rayleigh-Bénard System, Phys. Rev. Lett. 61 2550-2553 (1988).
MartinLof Martin-Lof, A., Statistical Mechanics and the Foundations of Thermodynamics, Lecture Notes in Physics 101, Springer-Verlag, Berlin (1979).
Mulero Mulero, A. (Ed.), Theory and Simulation of Hard-Sphere Fluids and Related Systems,
Lect. Notes Phys. 753 (Springer, Berlin Heidelberg 2008).
Mutabazi Mutabazi, I.,  Wesfreid, J.E. and Guyon, E. (eds.) Dynamics of Spatio-Temporal Cellular Structures: Henri Benard Centenary Review, in Springer Tracts in Modern Physics 207 (Springer, 2006). 
Puhl Puhl, A., Malek Mansour, M. and Mareschal, M., Quatitative comparison of molecular dynamics with hydrodynamics in Rayleigh.Bénard convection, Physical Review A, 40, 1999-2012 (1989).
Raja Rajagopal, K.R., A new development and interpretation of the Navier Stokes fluid which
reveals why the ‘‘Stokes assumption’’ is inapt, International Journal of Non-Linear Mechanics 50, 141–151  (2013).
Raja2 Rajagopal, K.R., Remarks on the notion of “pressure”, International Journal of Non-Linear Mechanics 71, 165–172  (2015). 
Rapaport_1986 Rapaport, D.C., Eddy formation in obstructed fluid flow: a molecular dynamics study, Physical Review Letters, 57, 695-698 (1986).
Rapaport_1988 Rapaport, D.C., Molecular Dynamics study of Rayleigh-Bénard Convection, Physical Review Letters, 60, 2480-2483 (1988).
Rapaport_1992 Rapaport, D.C., Unpredictable convection in a small box: Molecular dynamics experiments, Physical Review A, 46, 1971-1984 (1992).
Rapa09 Rappaport, D.C., The Event-Driven Approach to N-Body Simulation
, Progress of Theoretical Physics Supplement, 178, 5-14 (2009).
Risso Risso, D. and Cordero, P., Empirical determination of the onset of convection for a hard disk system, in Instabilities
and Nonequilibrium Structures IV, Tirapegui, E. and Zeller, W. (eds.). Kluwer Academic Publishers (1993).
Roman Román, F.L., White, J.A.  and Velasco, S., Fluctuations in an equilibrium hard-disk fluid: Explicit size effects. Journal of  Chemical Physics 107 (12), 4635-4641 (1997).
Stalling Stalling, D. and  Hege, H.C., Fast and Resolution Independent Line Integral Convolution, in Proceeding
SIGGRAPH '95 Proceedings of the 22nd annual conference on Computer graphics and interactive techniques,  Pages 249-256, (1995) ISBN:0-89791-701-4. Cabral, B. and Leedom, L.  Imaging vector fields using line integral convolution. Computer Graphics, 27, 263–272 (1993).
Wain71 Wainwright, T.E., Alder, B.J., and Gass, D.M., Decay of Time Correlations in Two Dimensions,  Physical Review A, 4, 233-237 (1971).

