$1/f$ noise on the brink of wet granular melting

The collective behavior of a two-dimensional wet granular cluster under horizontal swirling motions is investigated experimentally. Depending on the balance between the energy injection and dissipation, the cluster evolves into various nonequilibrium stationary states with strong internal structure fluctuations with time. Quantitative characterizations of the fluctuations with the bond orientational order parameter $q_{\rm 6}$ reveal power spectra of the form $f^{\alpha}$ with the exponent $\alpha$ closely related to the stationary states of the system. In particular, $1/f$ type of noise with $\alpha\approx-1$ emerges as melting starts from the free surface of the cluster, suggesting the possibility of using $1/f$ noise as an indicator for phase transitions in systems driven far from thermodynamic equilibrium.

1. Introduction f 1 noise (also called flicker or pink noise) is a class of signal exhibiting power spectrum a f with exponent a~-1, typically in the range a -< < -1.4 0.8 [1]. Due to its ubiquity in nature, f 1 noise has drawn much more attentions from various communities for almost a century, in comparison to white (a = 0), Brownian (a = -2) as well as diffusion noises [2][3][4][5]. The first and most extensively investigated noise of this class is electric noise in conducting or semiconducting materials [1,[6][7][8]. Later on, the ubiquity of f 1 noise is unveiled as a surprisingly large number of systems are found to fall into this class, ranging from blinking of stars and sunspot activity [2] in astrophysics, earthquake triggering [9] and undersea ocean currents in geophysics [10], to the loudness fluctuations of music [11], gene expression [12], and human cognition [13,14].
Because of this ubiquity, it is important to understand why such a common feature persists and to explore whether there is a universal mechanism behind or not. Following the widely used the Bernamont-Surdin-McWhorter (BSM) model [15][16][17][18], f 1 noise can be considered as a superposition of independent events with Lorentzian spectra and a broad distribution of relaxation time. Although it was originally introduced to explain electric noises, this model has been successfully applied to a large variety of systems exhibiting f 1 noise [1,19]. Nevertheless, recent investigations also reveal that the assumption of elementary events with Lorentzian spectra may not be necessary, since models based on random or point process [20,21] have also been reported to produce f 1 noise. In order to provide a general understanding of f 1 noise, Bak et al proposed the concept of self-organized criticality (SOC) [22,23], claiming that in spatially extended dissipative systems f 1 noise can be viewed as an indication of self-organized critical state. In contrast to the critical state in equilibrium thermodynamics, the system is self-organized, i.e. no external tuning is needed. Using a cellular automata (CA) model describing the avalanches in a pile of sand, they demonstrated that the scale invariance in time is associated with fractal structures, i.e. scale invariance in space. Triggered by this concept, a tremendous amount of investigations have been conducted to test the universality of SOC [19,24,25]. However, the CA model proposed by BTW was soon found to produce f 1 2 noise in both one-and two-dimensions [26,27] and thus cannot directly be used to explain f 1 noise. Moreover, contradictory results were found in model experiments on the avalanches of a sandpile [24,28,29]. Those investigations led to the conclusion that SOC is not as universal as it was claimed. Instead, its validity relies on the dissipation mechanisms of specific systems [30], and also on the way of analysis [31,32]. Thus from the SOC perspective, the ubiquitous f 1 noise is still not completely understood and whether a general theory for f 1 noise exists or not is still unclear. Nevertheless, the advance on the relation between the invariance of time and space gives us the hint that noise may provide useful insights into the selforganized stationary states of systems driven far from thermodynamic equilibrium [19].
Wet granular matter, such as wet sand used to build sand sculptures, has been drawing more and more interest from both physical and engineering communities in the past decades [33][34][35]. This is not only because of its ubiquity in nature, industries and our daily lives, but also due to the fact that it can be treated as a nonequilibrium model system with granular particles replacing molecules and liquid bridges formed between adjacent particles replacing molecular bonds. Because of the strong energy dissipation associated with wet impacts [36], continuous energy injection is necessary to maintain a certain stationary state of such a nonequilibrium system. Former investigations have revealed that the dissipative energy scale arising from the rupture of liquid bridges between adjacent particles plays an important role in determining the collective behavior of wet granular matter, such as melting [37,38], clustering [39], and phase separation [40,41]. More recently, an analogue of surface melting was found in agitated wet granular matter [42,43], suggesting that existing knowledge on phase transitions in equilibrium systems can be extended to explain the wide-spreading nonequilibrium systems in nature. Here, focusing on the noise spectra extracted the nonequilibrium stationary states of a wet granular cluster, I demonstrate the possibility of using f 1 type noise as an indicator for phase transitions far from thermodynamic equilibrium.

Methods
The experimental setup is sketched in figure 1. Various number N of polished soda lime glass beads (SiLibeads P) with a density of r = -2.58 g cm p 3 and a diameter of =  d 4 0.02mm are used as granular sample. After being mixed with a certain volume V l of purified water (specific resistance 18.2 MΩcm, LaborStar TWF), the particles are added into a cylindrical polytetrafluoroethylene container with an inner diameter of D = 102 mm and a height of 6 mm. The latter confinement ensures a strict mono-layer of particles. The liquid content = W V V l s with V s the total volume of the spheres is fixed at 0.02, which ensures short ranged attractive interactions between adjacent particles via the formation of liquid bridges [44]. The container is closed with a glass lid to keep the liquid content constant. During the experiments, the glass lid is slightly heated to avoid water vapor condensed there. The wet granular sample is illuminated from below with a LED array, and viewed from top with a CCD camera (Lumenera LU125M) mounted in the comoving frame (see figure 1(a) for a typical snapshot). The swirling table is leveled within´-5.7 10 3 degrees to avoid the influence from gravity. The horizontal swirling motion, i.e. a superposition of two perpendicular sinusoidal vibrations with identical amplitude A and a phase difference of p 2, is provided by a mechanical orbit shaker (Thermolyne, AROS160). The amount of energy injection into the system is adjusted by varying the driving frequency f d at a fixed amplitude A = 31.8 mm. With a computer controlled precision resistance decay box (Burster 1424), the swirling frequency is controlled with an accuracy of ∼10 −4 Hz. f d is measured by tracing a fixed particle on the swirling table with another CCD camera (Lumenera LU075M) mounted in the lab frame. A computer program is developed to adjust the swirling frequency in steps and capture images with a fixed rate of 5 frames per second.
Using an image processing procedure based on the Hough transformation [45], I determine the positions of all particles in each frame captured. For each particle, the bond orientational order parameters (BOOP) are calculated with [46][47][48] is an average of the local order parameter over all bonds connecting this particle to its nearest neighbors, which are identified with a critical bond length = r d Here corresponds to the spherical harmonics of a bond located at  r . Finally, I determine the local structure of each particle by comparing the order parameter q 6 with the standard values for perfectly hexagonal, square, line structures, as well as for isolated particles. Here, q 6 is chosen as the order parameter because of its sensitivity to the hexagonal order. As shown in figure 1(b),the particles are colored by their local structures after the above analysis in order to visualize the local structure changes.
The collective behavior of a wet granular cluster is obtained through initializing the system at a relatively large = t 700 s. The long initial time ensures a homogeneous wetting liquid distribution and a stable temperature gradient to avoid the condensation of water on the lid. Immediately after the lowest driving frequency is reached, f d is swept up with the same Df d and Dt . The whole cycle is repeated seven times for each particle number N. Figure 2 shows the collective behavior of a wet granular cluster at four different driving frequencies and corresponding internal structure fluctuations. The internal structure is characterized by P 6 , the percentage of particles in a local hexagonal structure.

Nonequilibrium stationary states
Hz, the cluster is composed of a few 'crystals' loosely connected with each other. Here, 'crystals' refer to clusters composed of particles in a local hexagonal packing. As the positions of the cluster center (see the cloud of dots in figure 2(a)) indicate, the mobility of the cluster is low. This suggests that the external driving force is not sufficient for the 'crystals' to overcome the friction from the ground and move freely. The P 6 order parameter stays at a high value close to 0.9, indicating that almost all particles are kept in a hexagonal structure. Thus, this state is called initial crystalline state. The small fluctuations of P 6 mainly arise from the interactions between the loosely connected 'crystals' and the frequent attaching and detaching of individual particles from them.
As the driving frequency f d increases to 1.285 Hz, the enhanced mobility of the particles leads to the merging of the 'crystals' into a single one with dislocation defects. The defect lines, along which local structures of particles deviate from hexagons to squares, suggest a certain amount of potential energy being stored there. This energy is obtained from the frequent collisions between the 'crystal' and the container wall, because of the higher mobility of the whole cluster in comparison to the initial crystalline state. Note that the state with defects is unstable: a small perturbation will lead to a healing of the defects due to the cohesion between adjacent particles. As the ground state with a perfect hexagonal packing is not easily achieved, the healing process typically ends up with a metastable state corresponding to a local potential energy minimum. The corresponding fluctuations of P 6 are much stronger than the initially crystalline state, because the defects associated structure changes occur at the length scale of the whole 'crystal'. Moreover, the fluctuations occur at time scales much longer than the swirling period, i.e. f 1 d . Such type of fluctuations persist as the driving frequency increases, as long as the cluster stays in the amorphous state. As f d increases further to 1.409 Hz, an abrupt transition into the surface melting state arises, owing to a balance between the energy injection and the rupture energy of a liquid bridge [43]. This transition is manifested by the collapse of the voids inside the 'crystal' together with the emerging of a liquid like layer covering the crystalline core, as shown in figures 2(c) and (g). Consequently, the whole cluster moves in a circular trajectory with permanent contact to the container wall. (Note that the shortest distance from the circle to to the container wall matches the mean radius of the cluster.) This feature indicates that the hard impacts between the cluster and the container wall with a finite effective coefficient of restitution (COR) is now replaced with soft ones (i.e. zero COR), because the liquid like layer acts as an effective damper.
It is remarkable to see that the fluctuations of P 6 also behave dramatically different in the surface melting state in comparison to the amorphous one: the large fluctuations at long time scales in the amorphous state are not present any more. This can be understood from the fact that the energy injection and dissipation in this state is localized in the molten surface layer of the cluster, because of the much more frequent collisions there. Hence, the metastability associated with the change of the crystalline structures at a length scale of the whole cluster is suppressed, leading to a white noise type fluctuations. Further increasing f d leads to a growth of the molten layer thickness, until the whole cluster behaves like a liquid at = f 1.719 d Hz. Along with the growth of the molten layer thickness, the center of the more deformable cluster is driven closer to the rim of the container, as a comparison between figures 2(c) and (d) reveals. Moreover, the transition into the liquidlike state leads to a decrease of the average P 6 , although the behavior of the fluctuations persists.
The various stationary states and the associated structure fluctuations are robust: variations of Df d , Dt by at least one order of magnitude and a change the sweeping direction of f d yield the same behavior. The lower limit of f d is selected such that the wet granular assemblies are immobile, i.e. the driving force cannot overcome the total frictional force from the ground.

Power spectrum of internal structure fluctuations
From the perspective of f 1 noise, it is intuitive to analyze the power spectrum of the internal structure fluctuations in various nonequilibrium stationary states of the wet granular cluster. As shown in figure 3, the power spectra of P 6 exhibit predominately power law behavior for a wide range of f d . In the initial crystalline or amorphous state (  f 1.285 d Hz), all the power spectra decay with an exponent of a » -1.5. The peak at = f f d indicates the influence from the driving frequency, as the energy injection is associated with periodic collisions between the cluster and the container wall. The peak becomes less and less pronounced as f d grows, which can be attributed to the 'softening' of the impact: as the 'crystal' becomes more amorphous, the effective COR for its impact with the container wall decreases, and the internal structure change of the cluster becomes more susceptible. This trend will eventually lead to more random collisions with the container wall and consequently a vanishing influence of the container wall. As a = -1.5 is typically associated with noise generated by diffusion mechanisms [3][4][5], it is interesting to explore possible diffusion mechanisms of driven wet granular particles in an amorphous state. For a monolayer of noncohesive dry granular particles, former investigations (see e.g. Reis et al [49]) reveal a non-trivial subdiffusive behavior in the vicinity of melting transition, which arises from the caging dynamics. However, it is still unclear what the diffusion mechanisms of wet granular particles in various nonequilibrium stationary states are and how it is related to the power spectra observed here. Further clarifications of these questions require a detailed characterization of the particle mobility and will be a focus of further investigations.
The exponent α decreases to »-1.2 as f d increases to 1.347 Hz, suggesting a f 1 type noise in the vicinity of melting transition. In the surface melting regime with even higher f d , α grows continuously until it saturates at »-0.3, suggesting a trend toward a white noise type of fluctuations without any long time correlations.
Meanwhile, a cut-off frequency of »0.1 Hz, below which the power spectrum decays much slower, arises.
As the power law exponent α is closely related to the stationary states of the system, it is intuitive to use it as an order parameter for the melting transition. Figure 4(a) (lower triangle) shows the fitted power law exponent as a function of the driving frequency for the spectrum shown in figure 3. The transition from an amorphous state to a melted state can be clearly distinguished from the deviation of α from its initial low value. To avoid the influence from the peak at f d , I fit the spectrum with an upper limit f lim , which is determined by decreasing f lim in steps from f d until the standard error of the fit converges. Because of the logarithmic scale of the power spectra, the data below the cut-off frequency play a minor role in determining α. Thus, a lower limit of f is not included in the fitting algorithm. A comparison to the results obtained by increasing f d (upper triangle) shows a good agreement, except for the region where the exponent starts to deviate from its initial low value. In this region, the slightly smaller α for the case of increasing f d indicates hysteresis, which is in agreement with a former investigation [43].
To check possible influence of the finite system size, I vary the particle number from 180 to 406, which corresponds to a range of area fraction from 0.277 to 0.624. As shown in figure 4(b), similar behavior of the power law exponent is found for all the area fractions explored. This demonstrates the robustness of using f 1 noise as an order parameter for the melting transition. More specifically, the trend of a rapid deviation from the initially low a » -1.5 into a saturated region persists. Hence, I fit the data with a constant value followed by a parabolic growth, and determine the threshold f th as the intersection point that minimizes the standard error. Quantitatively, the threshold frequency stays at »1.3 Hz, with larger error bars for small cluster sizes. Note that the smallest cluster with N = 180 corresponds to a radius of roughly seven particle diameters. For such a small 'crystal', the energy obtained from the container wall presumably leads to relatively strong fluctuations of the internal structures, and consequently large data scattering for the initial crystalline and amorphous states. For even smaller cluster size, the reduced probability for all particles to stay in one cluster hinders the statistical measurements performed here.
In the following part of this section, I try to provide a clue to the emerging f 1 noise and its connections to the dynamics of the wet granular clusters in the vicinity of the melting transition. As a starting point, we need to identify the essential ingredients leading to the f 1 noise. As already discussed in the above analysis, each assembly of the wet granular particles can be considered as one metastable state. This is reminiscent to the packing of equal sized spherical particles in three dimensions: The face-centered cube packing (fcc) with a density of 0.74 is rare [50]. From time to time, we obtain a distribution of metastable states with a packing density range from random close to loose packings [51,52]. The impacts with the container wall lead to frequent switches of the cluster from one metastable state to another one, which can be considered as individual activation-relaxation processes. For example, the collision of a cluster in the amorphous state with the container wall leads to a certain energy injection  µ - where ò is the effective COR. On one hand, E inj leads to an activation of defects inside the cluster, a corresponding reduction of P 6 , and an increase of the total potential energy through breaking of liquid bridges. On the other hand, it effectively 'heats' the cluster up, i.e. enhances the kinetic energy of granular particles therein. The latter effect leads to a healing of defects generated, because the structure with four fold rotational symmetry along the dislocation defects is locally unstable (see figure 2 and corresponding discussions). Moreover, the additional kinetic energy obtained from the healing process of a large defect may in turn generate small defects. A cascade of such energy transfer processes toward smaller and smaller length scales eventually leads to a relaxation of the injected energy through the dissipative interactions between individual particles. Following the BSM model [15][16][17][18], we may assume that the relaxation process after an impact has the following autocorrelation function where A and t 0 denote the activation energy and the relaxation time, respectively. Note that the two parameters may differ dramatically from one impact to another, because the relaxed state after one impact is not necessarily the one with the lowest potential energy. Instead, a distribution of various metastable states arises with each state corresponding to a certain configuration of the wet granular particles. As t ( ) C is an even function, one can then calculate the corresponding power spectrum for this process with ] , 0 1 2 , will lead to a limit on the power law decay: The spectrum flattens below t 1 2 and steepens above t 1 1 . Therefore, one can determine the intrinsic time scale of the relaxation process from the power spectra.
Note that the relaxation process through a cascade of defects nucleation and healing events relies on a certain rigidity of the cluster. As melting starts at relatively large f d , the spatial as well the associated time correlation vanishes, and the energy injection and dissipation occur locally in the molten layers close to the impact point. As a consequence, a white noise type fluctuations with vanishing time correlation arise. As the internal structure fluctuations can be considered as a superposition of individual activation-relaxation processes, the dependency of the power law exponent on the stationary states of the nonequilibrium system represents the change of the corresponding correlation time.
So far, we have a qualitative understanding on why the power law exponent can be used to determine the transitions between various stationary states in the nonequilibrium model system. However, one has to face the following obstacles toward a quantitative understanding: There exists possible coupling between individual activation-relaxation events. A hard impact with the container wall will change the structure of a cluster substantially, which in turn will influence the COR for the next impact. Hence the assumption of independent subsequent events may not be justified. Moreover, such an influence will also lead to a certain distribution of the activation factor A, i.e. the energy injection, and a coupling between A and the relaxation time t 0 . Therefore it is necessary to find the relation between A and t 0 for more quantitative comparisons to the experiments.

Conclusions
In conclusion, I investigate the self-organization of a two-dimensional wet granular cluster under horizontal swirling motion, and use the percentage of particles in a local hexagonal packing to characterize the internal structure fluctuations with time. In the frequency domain, the power spectra of all fluctuations exhibit predominately a power law behavior with decay exponent ranging from a » -1.5 to »-0.3. More specifically, the exponent is closely related to the nonequilibrium stationary states of the system: for the initial crystalline and amorphous states at low driving frequencies, it stays low at »-1.5, suggesting a diffusion noise. For the molten state at large driving frequencies, it saturates at »-0.3, suggesting a trend toward white noise. In the intermediate states associated with the melting transition, f 1 type of noise with a » -1 arises. A variation of the cluster size yields qualitatively the same power law behavior and quantitatively the same dependency of the decay exponent with f d . Such a robust behavior suggests an alternative way of detecting phase transitions in systems driven far from thermodynamic equilibrium. Finally, I show that the f 1 noise can be qualitatively understood from the impact mechanisms between the wet granular cluster and the container wall.
This investigation demonstrates the importance of noise in characterizing systems driven out of thermodynamic equilibrium. Concerning the robustness of such a characterization, future questions to ask are: Will any type of fluctuation related to the energy flux in and out of the nonequilibrium system, e.g. the sound energy arising from the inelastic collisions, exhibits a similar behavior? Will the way of energy injection or dimensionality of the system influence the behavior? Concerning a more quantitative understanding of the emerging f 1 noise, it is also important to characterize the coupling between the activation energy and the relaxation time experimentally through a measurement of the effective coefficient of restitution for wet granular clusters.