New Insight into Pseudo-Thermal Convection in Vibrofluidised Granular Systems

Utilising a combination of experimental results obtained via positron emission particle tracking (PEPT) and numerical simulations, we study the influence of a system’s geometric and elastic properties on the convective behaviours of a dilute, vibrofluidised granular assembly. Through the use of a novel, ‘modular’ system geometry, we demonstrate the existence of several previously undocumented convection-inducing mechanisms and compare their relative strengths across a broad, multi-dimensional parameter space, providing criteria through which the dominant mechanism within a given system – and hence its expected dynamics – may be predicted. We demonstrate a range of manners through which the manipulation of a system’s geometry, material properties and imposed motion may be exploited in order to induce, suppress, strengthen, weaken or even invert granular convection. The sum of our results demonstrates that boundary-layer effects due to wall (in)elasticity or directional impulses due to ‘rough’ boundaries exert only a secondary influence on the system’s behaviour. Rather, the direction and strength of convective motion is predominantly determined by the energy flux in the vicinity of the system’s lateral boundaries, demonstrating unequivocally that pseudo-thermal granular convection is decidedly a collective phenomenon.


S1.1: Particle Density Fields
In order to determine a two-dimensional density field for our system, we begin by subdividing the experimental volume into a series of equally-sized 'cells' in the plane of interest. By measuring the time spent by the tracer within each cell and normalising by the total duration of the experiment, we may determine a local 'residence fraction', F R , for each given region of the system.
For an ergodic, steady-state system, the mean particle density within any given cell may be assumed to be directly proportional to F R , allowing us to determine the number density, n, within the i th cell of the system as: where N is the total number of particles in the system as a whole and V c the cell volume. With a known particle diameter, d, the local packing density can then be determined as: η i = n i πd 3 6 (2)

S1.2: Velocity Fields
In order to visualise convection within our system, we must first determine the spatial distribution of velocities within the experimental volume. As with the particle density field determination detailed above, we begin by subdividing the system into a series of cells. For a given cell within the system, the mean local velocity may be determined simply by summing the velocities corresponding to all particle locations falling within the cell and dividing by the total number of data points located within the cell. By performing this process for the x-, yand z-components of velocity individuallty, velocity vector fields such as those shown in the main article may be produced. In order to produce the velocity fields presented in the main manuscript, our experimental system is divided into a series of three-dimensional cells of dimension ∆x × ∆y × ∆z. The precise size of the cells is varied between data sets chosen so as to ensure adequate statistics for the number of, and distribution of, data points acquired which, due to various factors such as run length, system material and tracer activity, may vary between experiments, and cannot be known a priori when performing experiments. For each cell, the velocities in the x-z plane of interest are calculated, and then averaged over all cells in the y-direction, ensuring data that is representative of the whole system. The x and z velocities v x and v z are then normalised by a factor v 2 x + v 2 z such that they may be plotted as unit vectors showing the direction of convection within the system. The plotting as unit vectors ensures that the direction of convection is easily visible -when plotted as unnormalised vectors, the varying lengths of the arrows representing the varying particle velocity can make the direction of convection hard to distinguis in low-velocity regions. Note that the convection velocities (see section 1.4, below) are, however, calculated using unnormalised velocities.

S1.3: Granular Temperature Fields
The granular temperature, also known as the 'fluctuant kinetic energy' of a system may be defined as: where c = |v −v|. Using PEPT data, the granular temperature within the i th cell of our system may be determined as: wherev i is the mean velocity of cell i determined as described in section 1.2, and N i is the total number of particle locations falling within the cell. By repeating this process for all cells within a system, two-dimensional granular temperature profiles such as those shown in the main text may be produced.

S1.4: Convection Velocities
In order to determine a value for the mean convective velocity of our systems, we begin by obtaining (unnormalised) 2D, depth-averaged velocity fields as described in section 1.2. We then determine the vertical centre of convection of the system, which corresponds to the height at which the minimum in |v x | occurs (see Fig. 1). Since, at the vertical centre of convection, we may assume all convective velocity to be in the vertical direction, and we know that the mass flow rate within the system must be conserved, we may compute the mean convective flow rate,v c , of the system, following the methodology of Hsiau and Chen 4 , simply as the mean magnitude of the z-component velocity averaged over all cells lying at this height.
2/5 1.5 S1.5:Additional quantities While the above sections detail the manners in which the quantities relevant to the present work may be determined, a number of other valuable parameters may additionally be extracted from PEPT data. For further information regarding the manners in which these other quantities may be computed, and for more information regarding the PEPT technique in general, please refer to our references 1, 2, 5-8 .

S1.6: Reliability and Repeatability of Experiments
In order to ensure the reliability, generality and repeatability of our experimental results, a two-step process is undertaken. Firstly, all experimental data sets are repeated. Repeat experiments are conducted on different days (to ensure an independence from precise humidity values and other imperfectly controllable external variables which may vary from day to day * ) and in inverse order (to ensure that wear on the system and/or particles does not influence our results).
Secondly, in order to ensure that our experimental runs are long enough to produce suitable statistics, each data set of duration t f ull is divided into two non-overlapping partial data sets each of length t partial = 1 2 t f ull . Three one-dimensional packing density profiles (the simplest spatially-variant whole-field parameter extractable from PEPT data) are then calculated from the three data sets -one each for the full data set and each partial set -using the same cell size (see section 1.1). The full-data and partial-data profiles are then compared against one another point-by-point, allowing a mean percentile error to be calculated. Note that only points falling within the bulk of the medium (here defined as twice the vertical centre of mass of the system -a measure shown viable in previous studies 7 ) are considered. If the mean percentile error calculated in this manner lies within 5% for both partial data sets, the data is assumed consistent, i.e. the experiment has been run for a suitably long time.
Once adequate statistics within each repeated experiment have been established, a similar cross-comparison of said repeated experiments is then conducted to ensure reproducibility.

S2: Discrete Element Method Simulations
In this work, we employ discrete element method (also known as 'discrete particle method', DPM) simulations to computationally model our system. While the system modelled is described in the main text, we discuss here in greater detail the specific algorithms used to model the behaviour of the simulated particles and housing geometry, including the relevant normal and frictional force laws employed.

S2.1: Simulation Details
The simulation data utilised in this paper is acquired using the MercuryDPM discrete element method software package 9-12 . In the current study, interactions between particles are modelled through the implementation of a frictional spring-dashpot model * Note that, due to the large particle size used in our experiments, and the temperature-controlled nature of the laboratory used to perform experiments, such matters are highly unlikely to affect our results in any significant manner.

3/5
with linear elastic and dissipative contributions, the relevant normal and tangential forces being calculated, respectively, as 13,14 : and Here, the symbols v n i j and v t i j represent, respectively, the normal and tangential components of (relative) velocity between a given pair of interacting particles; δ t i j is the elastic tangential displacement, whose definition and significance are described in detail in reference 15 ); k n and ζ n represent, respectively, the relevant spring constant and damping constants, whose values are determined as: In the above, ε is a (user-defined) restitution coefficient, and t c is the contact time 16 , whose value is also determined by the user. In the present work, t c is chosen so as to give a particle stiffness of 1.5 × 10 5 N/m, a value found to provide a realistic degree of overlap (δ n i j ) between colliding particles. In the above equations, the variable m i j corresponds to the reduced mass of a pair of colliding particles i and j, and can be calculated as: The tangential spring and damping constants are calculated from their normal counterparts as k t = 2 7 k n and ζ t = ζ n . We implement the Coulomb friction law by applying a static yield criterion which acts to truncate the magnitude of δ t i j at a value equal to µ f n i j , such that the inequality f t i j ≤ µ f n i j is fulfilled Finally, we consider the force of gravity, which acts uniformly on all particles within the system in a direction along the negative z-axis as defined in the main text.
The force relations obtained above are integrated in time using a Velocity-Verlet time-stepping algorithm 17 with a step size δt = tc 50 to model the evolution of the velocities and positions of the particles within the system.

S2.1: Validation of Simulation Model
For all experimental cases studied, simulations were performed implementing identical system parameters, in order to assure a qualitative match in the behaviours of interest between simulation and experiment -i.e. that our numerical systems are representative of our 'real', physical experiments. Specifically, it was ensured that the number and orientation and positioning of the convection rolls, where present (and the absence of convection where relevant) showed a direct correspondence for an identical parameter set. In Fig. 3, we show equivalent numerical and experimental data sets for the 4 main experimental cases discussed in the main text: flat, dissipative walls (normal convection), flat walls with ε wall = ε particle (zero convection), oscillating sawteeth with upward orientation (inverse convection) and oscillating sawteeth with downward orientation (inverse convection). As can be seen from the images provided, our simulations indeed behave as expected based on our experimental results.  Figure 3. Comparison of experimental (above, black arrows) and numerical (below, red arrows) results for equivalent systems, namely -from left to right -bounded by dissipative PMMA sidewalls, relatively elastic steel sidewalls, upward-oriented oscillating sawteeth and downward-oriented oscillating sawteeth.