Ordering in Granular Rod Monolayers Driven Far from Thermodynamic Equilibrium

Ordering in Granular Rod Monolayers Driven Far from Thermodynamic Equilibrium

Thomas Müller Experimentalphysik V, Universität Bayreuth, D-95440 Bayreuth, Germany    Daniel de las Heras Theoretische Physik II, Universität Bayreuth, D-95440 Bayreuth, Germany    Ingo Rehberg    Kai Huang kai.huang@uni-bayreuth.de Experimentalphysik V, Universität Bayreuth, D-95440 Bayreuth, Germany
September 8, 2019

The orientational order in vertically agitated granular rod monolayers is investigated experimentally and compared quantitatively with equilibrium Monte Carlo simulations and density functional theory. At sufficiently high number density, short rods form a tetratic state and long rods form a uniaxial nematic state. The length-to-width ratio at which the order changes from tetratic to uniaxial is around in both experiments and simulations. This agreement illustrates the universal aspects of the ordering of rod-shaped particles across equilibrium and nonequilibrium systems. Moreover, the assembly of granular rods into ordered states is found to be independent of the agitation frequency and strength, suggesting that the detailed nature of energy injection into such a nonequilibrium system does not play a crucial role.

45.70.-n, 05.70.Ln, 64.70.M-

I Introduction

The ordering of anisotropic particles is a universal phenomenon appearing widely in nature, ranging from thermally driven molecules or colloids Stephen and Straley (1974); Gennes and Prost (1995); Onsager (1949); Vroege and Lekkerkerker (1992) to active particles such as bacteria colonies Zhang et al. (2010), actin filaments Schaller et al. (2010); Sanchez et al. (2012), animal groups Buhl (2006); Couzin et al. (2005); Liu et al. (2013), and living liquid crystals Zhou et al. (2014). In equilibrium lyotropic systems, such as hard rods interacting only through excluded volume interactions, the transition of sufficiently anisotropic particles into various ordered states is entropy driven. The loss in rotational degrees of freedom in the ordered state is compensated by the gain in the translational ones Onsager (1949); Frenkel (1999); Vroege and Lekkerkerker (1992). Taking a two-dimensional system of hard rectangles as an example, a tetratic state with four-fold rotational symmetry was discovered in Monte Carlo (MC) simulations Wojciechowski and Frenkel (2004); Donev et al. (2006), and studied theoretically with density functional theory (DFT) Martínez-Ratón et al. (2005); Martínez-Ratón and Velasco (2009); Geng and Selinger (2009). The number density and the length-to-width ratio (aspect ratio) of the particles were found to be the key parameters determining the ordered states of hard rectangles with only excluded volume interactions Martínez-Ratón et al. (2005). Given the ubiquity of ordering transitions in nature, it is important to ask how well the existing knowledge about such transitions in equilibrium (thermal) systems can be extended to nonequilibrium (athermal) systems.

Due to the dissipative interactions between particles, agitated granular matter has been frequently used as a nonequilibrium model system for phase transitions Jaeger et al. (1996); Ristow et al. (1997); Götzendorfer et al. (2006); Eshuis et al. (2007); Fingerle et al. (2008); Huang et al. (2010); May et al. (2013). Rich and often counterintuitive dynamical behavior Börzsönyi and Stannarius (2013) has been discovered for granular rods, including vortex patterns Blair et al. (2003), collective swirling motions Aranson et al. (2007), giant number fluctuations Narayan et al. (2007); Aranson et al. (2008), violation of the equipartition theorem Harth et al. (2013), and an enhanced ordering transition in an effective ‘thermal’ bath of spherical particles Kumar et al. (2014). Reminiscent to equilibrium systems, ordering transitions of vertically agitated granular rods were investigated in three-dimensional (3D) and quasi-two-dimensional systems. In 3D, the aspect ratio of the rods was found to influence the ordered states of cylindrical rods Yadav et al. (2013). In quasi-two-dimensional systems, a bulk isotropic-uniaxial nematic (I-U) transition was observed for cylindrical rods with large aspect ratios Galanis et al. (2006) and an effective elastic constant was characterized quantitatively Galanis et al. (2010a). Particularly in strict monolayer systems, the shape of the rods was found to play an important role in determining the ordered states: Tetratic, nematic or smectic order was found for cylindrical rods, tapered rods or rice particles, respectively Narayan et al. (2006). Moreover, tetratic order was also found for tubular shaped particles and the influence of the container shape was discussed in Sánchez and Huerta (2014).

Despite all of these progresses, it is still unclear to which extent one can draw quantitative connections between systems in and out of thermodynamic equilibrium. More specifically, a quantitative comparison between the state diagram of dissipative granular rods and that of the corresponding equilibrium system is still lacking. This quantitative comparison is the purpose of the present work. Here we investigate experimentally the orientational order in monolayers of cylindrical granular rods driven far from thermodynamic equilibrium, and compare the results to MC simulations as well as DFT of the analogous equilibrium system. Focusing on the bulk region of the system, we detect both tetratic and uniaxial nematic states by varying the aspect ratio of the rods. We demonstrate that the aspect ratio and the number density of rods are the key parameters determining the state diagram in both systems. We find a common aspect ratio that separates tetratic and uniaxial nematic states in both experiments and MC simulations. Such an agreement illustrates the universal aspects of the ordering of rod-shaped particles.

Ii Methods

ii.1 Experiments

Figure 1: (Color online) Sketch of the experimental set-up. The closed cylindrical container of height and radius is driven sinusoidally against gravity with an electromagnetic shaker. The rods have a length and a diameter . The embedded image shows a close view of the detected particles.

A sketch of the experimental set-up is shown in Fig. 1. Monodisperse polyvinyl chloride (PVC) rods of diameter and length , cut from welding wires of  mm (aspect ratio ) or  mm (), are confined in a cylindrical container of height and radius  cm. The ratio is chosen for both diameters to ensure a monolayer of particles; that is, no rods can cross or jump over each other. The inner surface of the container is covered with antistatic spray (Kontakt Chemie, Antistatik 100) to minimize electrostatic forces. An electromagnetic shaker (Tira TV50350) is employed to drive the sample sinusoidally against gravity with frequency  Hz and peak acceleration , where is the peak vibrational amplitude and is the gravitational acceleration. The acceleration is monitored with an accelerometer (Dytran 3035B2). We capture high contrast images of the rods using backlight LED illumination and a camera (IDT MotionScope M3) mounted above the container. The camera is synchronized with the shaker so as to capture images at a fixed phase of each vibration cycle. The images are subjected to an analysis algorithm that determines the center of mass and the orientation of the -rod with . is the angle of the main rod axis with respect to a fixed laboratory axis, and is the total number of rods in the container. The detection rate is  % for  mm and  % for  mm.

To systematically study the collective behavior of the rods, we vary the global area fraction between and , and the aspect ratio between and . For each and , we vary the peak acceleration with a step of from to and back. The waiting time between each step is fixed at  minutes. We repeat the whole cycle at least times.

ii.2 Monte Carlo Simulations

Correspondingly, we model the particles as two-dimensional hard rectangles of length and width that interact through excluded volume interactions. of such particles are placed in a box with dimensions and along the - and -axes, respectively. We use periodic boundary conditions along both axes and study the equilibrium bulk configurations by means of standard MC simulations Allen and Tildesley (1987) in the canonical ensemble. That is, we fix the number of particles and the system area (the temperature is irrelevant in hard models). The number of particles is similar to that in the experiments, . We use simulation boxes with rectangular and square shapes. No difference has been found between the two geometries.

The simulation method is as follows. In order to equilibrate the system we start at very high area fractions, . We place the particles, with their main axes pointing in the same direction, in a rectangular lattice. Next we run Monte Carlo steps (MCSs). Each MCS is an attempt to move and rotate all the particles in the system. The maximum displacement and maximum rotation that each particle is allowed to perform in a MCS is determined such that the acceptance probability is . Then we remove a few randomly chosen particles, recalculate and , and start a new simulation. The number of removed particles is such that the change in area fraction is . In order to rule out metastable configurations related to the preparation of the initial state, we discard simulations with . When the area fraction is below that limit we start the proper simulation. For each simulation we first run MCSs to equilibrate the system and then accumulate data over MCSs. For selected we have also simulated the system by increasing the number of particles, i.e., by adding particles instead of removing them. We have found no differences between both methods.

ii.3 Density functional theory

We use an Onsager-like DFT with Parsons-Lee rescaling. A similar DFT was previously used to analyze the state diagram of two-dimensional rods confined in a circular cavity de las Heras et al. (2009). We are interested in the behavior of fluid states in which the density is spatially homogeneous. Hence we can write, without loss of generality, the one body density distribution as


where is the number density and is the orientational distribution function. Here is the angle with respect to the director. is normalized such that


We split the free energy into two parts


where is the ideal gas part and is the excess part accounting for the excluded volume interactions. The ideal free energy per unit of area is given exactly by


where with the Boltzmann’s constant and the absolute temperature. is the (irrelevant) thermal volume that we set to one. The excess part is approximated by


is the excluded area between two rectangles with relative orientation :


and is the excess free energy per particle of a reference system of hard disks at the same area fraction as our system of hard rectangles. The diameter of the disks is selected such that both disks and rectangles have the same area. Following Baus and Colot Baus and Colot (1987) we approximate by:


with . Eq. (5) recovers the Onsager approximation in the low density limit.

Finally, the grand potential is given by


with the chemical potential. We minimize with respect to and in order to find the equilibrium states. We use a standard conjugated gradient method to minimize the functional. We use a truncated Fourier expansion to describe . We truncate the expansion such that the absolute value of the last coefficient in the expansion is smaller than .

Iii Results and Discussion

This section is organized as follows: We first introduce the ordered states observed in experiments and MC simulations in section III.1. In section III.2, we analyze the influence of the container walls and the driving conditions in the experiments. Finally in section III.3, we quantify the ordering transition threshold for various aspect ratios and compare the state diagrams obtained experimentally, via MC simulations and with DFT.

iii.1 Ordered states

Figure 2: (Color online) Raw experimental images (topview) showing typical configurations of rods with aspect ratio  (a) and  (b) at high global area fractions. The yellow (light gray) dashed circle indicates the region of interest.

Figure 2 shows typical snapshots of the ordered states obtained experimentally. Short rods (a) tend to develop tetratic order with two alignment directions perpendicular to each other. Long rods (b) form uniaxial nematic order with only one preferred alignment direction. In both cases, the container promotes either homeotropic (perpendicular) or planar (parallel) anchoring of the rods close to the boundary. To minimize the boundary effects, we consider only those particles located in the central region of the container, as marked in Fig. 2. A quantitative justification of this region of interest (ROI) will be given in section III.2. Sometimes during the experiments, especially at low global area fractions, we observe regions with very low number density of rods (almost empty regions). As we are interested in the bulk behavior, we discard those configurations in which the “empty regions” and the ROI overlap.

Figure 3: (Color online) Typical snapshots of tetratic (left column, ) and uniaxial nematic (right column, ) states. (a) and (d) are reconstructed from the positions and orientations of the particles detected in the center region of the container. (b) and (e) are from MC simulations with periodic boundary conditions. The particles are color coded according to their orientations. (c) and (f) show the orientational distribution functions of the particles in experiments (gray bars) and simulations (solid line). is the angle with respect to the director.

Figure 3 shows a direct comparison of the ordered states obtained in both experiments and MC simulations. The color coded particle configurations are reconstructed from granular rods in the ROI (upper panels) and from MC simulations (middle panels) with periodic boundary conditions. In the tetratic state with fourfold rotational symmetry (left column), the orientational distribution function , where is the angle with respect to the director , has two peaks at and (c). In contrast, in the uniaxial nematic state (right column), the elongated particles are oriented on average along the director, yielding only one peak at (f). The director is calculated as the normalized eigenvector of the largest eigenvalue of the tensorial order parameter . Here is the th Cartesian coordinate of the unit vector , is the Kronecker delta, and denotes an average over the rods Cuesta and Frenkel (1990); de las Heras and Velasco (2014). To quantify the orientational order we measure


where and are the uniaxial and tetratic order parameters, respectively. In an isotropic state (no orientational order) and vanish. In a uniaxial nematic state and . Finally in a tetratic state and . The states in Fig. 3 are selected such that and are comparable in both experiments and MC simulations.

iii.2 Experiments: The influence of boundary and driving

Figure 4: (Color online) Wall-rod angular correlation function as a function of the rescaled distance to the wall for various . A sketch with various definitions is shown in the inset. The data are obtained through an average over all , global area fractions and experimental runs. The typical error for () is comparable to the size of the symbols for . Note the logarithmic scale on the y-axis.

Experiments Galanis et al. (2006); Narayan et al. (2006); Galanis et al. (2010b); Kudrolli et al. (2008) and MC simulations de las Heras and Velasco (2014) have shown that the container induces a preferential alignment of the particles close to the wall. In order to facilitate the investigation in the bulk, we first need to characterize such an influence quantitatively.

Following the ideas in Galanis et al. (2006), we calculate the wall-rod angular correlation function , where is the shortest distance from the rod center to the container wall, the angle quantifies the tangential direction of the corresponding point on the wall (see inset in Fig. 4), and denotes an average over all the particles at a distance . Either homeotropic or planar alignment of the particles with respect to the wall results in . In Fig. 4, is presented as a function of the rescaled distance to the wall with a binning width of . For all aspect ratios investigated, decays exponentially with . To minimize the influence of the wall, we consider only those particles with to be in the ROI. In this region, is always smaller than and remains in a range comparable to the experimental uncertainties. We characterize the state of the system by measuring the area fraction and in circular regions with radius inscribed in the ROI. Subsequently, we calculate from accumulated over all the regions that share the same .

Figure 5: (Color online) Tetratic and uniaxial order parameters in the ROI as a function of the area fraction in a system of rods with aspect ratio : (a) For three ranges of at  Hz and (b) for four with accumulated over all . The threshold is obtained through fits to the data (straight lines) accumulated over all for  Hz (see text for details). The order parameters are not exactly zero in the isotropic state due to the finite size of the system Cuesta and Frenkel (1990).

Figure 5 shows the order parameters as a function of for short rods with . It indicates an area fraction above which the tetratic order parameter grows from its initial low value, while the uniaxial order parameter remains low. Such a combination of and suggests a gradual isotropic-tetratic (I-T) transition. As shown in (a), the behavior of and does not depend on the peak vibration acceleration. This is further confirmed through a comparison among data obtained for all in the range of and also for all aspect ratios investigated. As shown in (b), a variation of the vibration frequency from  Hz to  Hz for also yields the same behavior of .

Such agreements indicate that the details of how the rods are effectively ‘thermalized’ in our nonequilibrium system are not essential in determining the ordering transitions, providing us the opportunity to draw connections to the corresponding equilibrium systems. Accordingly, we accumulate the data over all at  Hz for a more accurate characterization of the transition threshold . By fitting with a constant value in the isotropic region and with a straight line in the ordered state, we obtain as the intersection point which minimizes the standard error. Only data with sufficient statistics (i.e., error bar ) and are chosen for the fits.

Moreover, the height of the container is found to play a minor role in determining the ordering transition: A variation of from to leads to the same behavior of . Experiments with for and give rise to slightly lower transition thresholds . More specifically, we find a decrease of for short rods and of for long rods, which is in both cases within the uncertainty of the fit. In addition, for a specific aspect ratio of , the same experiments have been performed for two different rod diameters. The results agree with each other within the error bar, suggesting that the mass of the rods does not play a dominating role in the ordering transition.

Figure 6: (Color online) Tetratic and uniaxial order parameters as a function of the area fraction for (a) and (b) . The insets show the corresponding results from MC simulations. The experimental data is an accumulation over all at  Hz. Straight lines are linear fits to determine the threshold .

iii.3 Experiments vs. simulations and DFT

Based on the above characterizations of the boundary influence, we compare the ordering transitions of granular rods in the ROI to the corresponding thermal system. Figure 6 shows the averaged order parameters obtained in both experiments and MC simulations (insets) for rods with (a) and (b). As discussed above, tetratic ordering occurs in the system of short rods. For long rods, both order parameters start to grow above , suggesting a gradual I-U transition. Qualitatively, the agreement between experiments and MC simulations on the behavior of both tetratic and uniaxial order parameters is remarkable for both aspect ratios. Such similarities indicate that the ordering of granular rods is governed by the geometric constrain of non-overlapping rods, which is the only interaction considered in the simulations. Quantitatively, the threshold obtained experimentally for rods with agrees with the one obtained from MC simulations within the error. However, the experimentally obtained threshold for rods with is larger than the one obtained for the corresponding thermal system, .

As and are the key parameters determining the state of the system, we compare the experimental (nonequilibrium) results with the MC (equilibrium) simulations in a state diagram shown in Fig. 7. In both systems short rods form a tetratic state and long rods a uniaxial nematic state at sufficiently high area fractions. The aspect ratio at which the ordered state changes from tetratic to uniaxial nematic agrees quantitatively. It is found to be in both experiments and simulations 111This value represents the average between (uniaxial) and (tetratic).. This result agrees with previous simulations in which a tetratic phase was found for and some evidence of uniaxial ordering for Martínez-Ratón et al. (2006). The quantitative agreement of across systems in and out of thermodynamic equilibrium illustrates the universal aspects of the ordering transitions.

On the other hand, the threshold for agitated rods differs from that in MC simulations, indicating the non-universal aspects of the ordering transitions. First, the experimentally determined exhibit a peak around . In contrast, MC simulations show a monotonic decay with . Second, measured in experiments deviates systematically from that obtained via MC simulations as grows. For the largest aspect ratio investigated experimentally, , much higher area fraction is required for the uniaxial state to develop. This difference might be attributed to the following mechanisms. (i) The strong fluctuations in the nonequilibrium steady states of granular rods may lead to temporal disorder in a system that could in principle relax into an ordered state. (ii) Due to the dissipative rod-rod interactions, the tendency of clustering for granular rods is larger than in MC simulations, especially for large [compare panels (d) and (e) of Fig. 3]. (iii) Finally, the container wall may frustrate the orientational order of the agitated rods in the entire cavity. Further experiments using containers with different sizes and shapes might shed light on such a discrepancy.

Concerning the fluctuations, it is known that the velocity distributions of agitated granular spheres are non-gaussian and exhibit exponential tails, no matter whether the particles form clusters Olafsen and Urbach (1998) or not Kudrolli and Henry (2000). As the dissipative nature does not depend on the shape of the particles, we expect a similar behavior in our system. This feature sets agitated granular rods apart from thermally driven liquid crystals, and triggers the question of how to define an effective ‘thermal’ energy scale for a nonequilibrium system. Monitoring the mobility of individual granular rods with high speed photography could help to shed light on the difference between equilibrium and nonequilibrium systems found here.

Figure 7: (Color online) State diagram in the plane of aspect ratio , and area fraction obtained via experiments of agitated granular rods (circles) and MC simulations of thermally driven hard rectangles (triangles). The labels denote the states: isotropic (I), uniaxial nematic (U), and tetratic (T). They are colored according to the experimental data. The diameter of the rods used in the experiments is  mm for (gray symbols) and  mm for (black symbols). Closed and open symbols indicate the I-T and I-U transitions, respectively. The inset shows the state diagram of equilibrium hard rods according to DFT in comparison to MC simulations in an extended region of . Dashed lines are continuous transitions and solid lines denote first order transitions.

In the inset of Fig. 7 we show the state diagram according to DFT together with the thresholds obtained from MC simulations in an extended region of . It is similar to the one predicted by the scaled particle theory Martínez-Ratón et al. (2005). DFT also predicts I-T transitions for small and I-U transitions for large . However, the tetratic state is stable only for , most likely because only two-body correlations are considered in the theory Martínez-Ratón and Velasco (2009); Martínez-Ratón et al. (2006). Concerning the ordering transition threshold , there is a good agreement between DFT and MC simulations for . For low aspect ratios, the deviations between both approaches are due to the mean field character of the theory. For , DFT predicts a T-U transition at very high area fractions. Due to the limitations in both experiments and MC simulations, the region of very high area fractions, where the T-U transition may arise, has not been explored.

Iv Conclusions

To conclude, the ordering in agitated granular rod monolayers is found to be determined predominately by the aspect ratio of the rods and the area fraction, while the frequency and the strength of the agitation are not essential. It suggests that the detailed nature of energy injection into such a nonequilibrium system is not important, analogous to the role that temperature plays in equilibrium hard rod models. In comparison to previous experimental investigations on monolayer systems, we have focused on the bulk region of the container and found both tetratic and uniaxial nematic order for cylindrical rods. This enables a direct comparison to the state diagram of the corresponding equilibrium system. We find that, depending on whether the aspect ratio is smaller or larger than , a gradual isotropic-tetratic or an isotropic-uniaxial nematic transition arises as the area fraction grows, in both experiments and simulations. This agreement suggests some degree of universality for the ordering of rod shaped particles across systems in and out of thermodynamic equilibrium. Nevertheless, we have also found a qualitative difference between both systems, namely the trend of the area fraction threshold at the ordering transitions.

Further investigations will focus on characterizing the area fraction and velocity fluctuations of the system, in order to find an effective ‘thermal’ energy scale for such a nonequilibrium system. Moreover, a comparison to molecular dynamics simulations Volfson et al. (2004) with tunable rod-rod dissipation energy could help to elucidate how fluctuations influence the ordering transition threshold.

The authors would like to thank Wilhelm August for the preliminary work on the experimental set-up. Inspiring discussions with M. Schmidt, D. van der Meer, and C. Krülle are greatly acknowledged. TM and KH acknowledge the support from the DFG through Grant No. HU1939/2-1.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description