Emergent Bistability and Switching in a Nonequilibrium Crystal

TL;DR

This paper demonstrates that global bistability in nonequilibrium systems can occur without bistable elements, through experimental and simulation evidence of switching between crystalline and gas-like states driven by conservative forces, damping, and noise.

Contribution

It reveals that system-wide bistability can emerge without local bistable components, challenging previous assumptions about multistability origins.

Findings

Experimental observation of state switching in microparticle system

Simulations show forces, damping, and noise cause bistability

Bistability occurs over seconds to hours

Abstract

Multistability is an inseparable feature of many physical, chemical and biological systems which are driven far from equilibrium. In these nonequilibrium systems, stochastic dynamics often induces switching between distinct states on emergent timescales, for example, bistable switching is a natural feature of noisy, spatially-extended systems that consist of bistable elements. Nevertheless, here we present experimental evidence that bistable elements are not required for the global bistability of a system. We observe temporal switching between a crystalline, condensed state and a gas-like, excited state in a spatially-extended, quasi-two-dimensional system of charged microparticles. Accompanying numerical simulations show that conservative forces, damping, and stochastic noise are sufficient to prevent steady-state equilibrium, leading to switching between the two states over a range of time scales, from seconds to hours.

Full text

Systems driven far from equilibrium, either through external fields or from active constituents, exhibit a diverse range of collective phenomena. Although no universal framework exists for describing nonequilibrium systems, in many cases a steady state is reached where some aspects of equilibrium theories can be applied Mukamel (2000). However, nonequilibrium systems often never reach a steady state, and constantly switch behavior Hinrichsen (2000). In biological systems, switching is a common theme in the regulation of the cell cycle Pomerening et al. (2003), gene expression Maamar et al. (2007), and neural activity Pigorini et al. (2015). Certain system features can affect switching behavior over time, for example, stochastic noise can both facilitate and hinder switching Lindner et al. (2004). The role of noise is often discussed in terms of stochastic resonance or stochastic facilitation, where noise serves a constructive role to amplify switching between two states Gammaitoni et al. (1998); McNamara et al. (1988); McDonnell and Ward (2011).

For spatially-extended systems such as a neural network, the existence of collective switching between states is supported by the bistable nature of the individual neurons. Other well-known examples include arrays of diode resonators Löcher et al. (1996) and nematic liquid crystals Sharpe et al. (2003). More interestingly, bistable switching can emerge from collective system behavior, even if the system constituents are not bistable. Bistability has been demonstrated in several numerical studies on collective systems composed of monostable elements Vilar and Rubi (1997); Cubero et al. (2009); Miloshevich et al. (2009). However, observations of similar behavior in experimental systems are lacking. Using an electrostatically-levitated, quasi-two-dimensional collection of charged microparticles, here we show how dynamical switching between distinct states can emerge as a collective system behavior rather than a consequence of bistable elements.

Experiments — Figure 1a shows the experimental conditions for particle levitation. The micron-scale particles are levitated in a weakly-ionized argon plasma Shukla and Eliasson (2009); Merlino and Goree (2004). Direct visualization of the individual particles allows using such “dusty” or complex plasmas as model systems for studying various condensed matter phenomena Joyce et al. (2002); Sütterlin et al. (2009). Gas pressure and plasma density control the residual charge on the particles, and thus their ability to form a crystalline lattice with little or no Brownian motion (Fig. 1b), or a kinetic gas-like state (Fig. 1c). The particles interact via a repulsive, screened Coulomb potential, where the typical screening length is $\lambda_{D}\approx$ 1 mm. The particles levitate at a vertical position defined by a balance between gravity and electrostatic repulsion from the negatively-biased aluminum electrode at voltage $V_{bias}$ (Fig. 1d). Weak horizontal confinement is provided by a gradient in ion density which keeps the negatively-charged particles near the center of the electrode.

Our experiments use melamine-formaldehyde particles obtained from microParticles GmbH. The mean diameter of the particles was 9.46 $\mu$m, with a coefficient of variation, $c_{v}\approx 1.1\%$. This inherent polydispersity led to a variation in the equilibrium height of the particles and their corresponding vertical oscillation frequency. To remove “dimer” particles consisting of two particles stuck together, we pulsated the plasma using the duty cycle on the RF power supply (Seren model R301, 13.56 MHz). The duty cycle was adjusted so that heavier particles impacted the electrode, while smaller particles were levitated once the plasma was re-initiated. The bias voltage on the electrode was adjusted by varying the RF power ($\sim$ 2-10 W) delivered to the plasma. Lower voltages were achieved by detuning the matching network to reduce the power. A 300 mW diode laser was used as a source for the horizontal and vertical laser sheets (Laserglow Technologies). The thickness of the laser sheet was $\approx$ 100 $\mu$m. The particle motion was imaged from above and from the side using USB 3.0 digital video cameras (Point Grey). We tracked the particle positions using a python-based software package Allan et al. (2016). Particle velocities were calculated by subsequent differences in position in each frame.

Often condensed matter systems have some intrinsic disorder and fluctuations. In complex plasmas, variations in particle size introduce quenched disorder, leading to non-reciprocal interactions Ivlev et al. (2015) and melting of the crystal under weak confinement Couëdel et al. (2010); Liu et al. (2010). More importantly for our experiments, single particles experience self-induced vertical oscillations at low gas pressures. A number of proposed mechanisms for the source of the vertical oscillations have been suggested, initially, a delayed charging of the particles Nunomura et al. (1999) and subsequently, spatial charge variation Samarian et al. (2001) and plasma sheath inhomogeneities Resendes and Sorasio (2002). Although there is no consensus on the source of the oscillations, they are crucial to the emergence of bistability in our experiments. Figure 1e shows the amplitude of vertical oscillation for a single particle as a function of pressure, $P$, for two different values of $V_{bias}$. The inset shows the $z$-position of a particle in time, where the characteristic oscillation frequency is $\approx$ 20 Hz.

By varying $P$ and $V_{bias}$, which control the dissipative and interaction forces, respectively, we explore the system behavior (Fig. 1f). A similar phase diagram for smaller particles is shown in Fig. S1. For large values of $P$ and $-V_{bias}$, the particles form a stable crystalline monolayer. Initially, reducing $P$ or $-V_{bias}$ destabilizes the monolayer as small fluctuations become more pronounced. However, at low $P$ and $-V_{bias}$, $\lambda_{D}$ increases due to the low plasma density, leading to an increase in the inter-particle distance and the formation of stable crystals. Surprisingly, in this regime we also observe a “bistable” phase, where the particle monolayer continually switches between a crystalline state and a gas-like state on a minutes-long timescale. This timescale is not inherent in the dynamics of individual particles, but rather emerges from the many-body dynamics of a stochastically-driven system with damping.

Figures 2a-f show images from a typical bistable sample over 112 s (see Supplementary Videos 1 and 2 for the full evolution of similar samples). The vertical oscillations play a key role here. As the particles begin to oscillate vertically, they move in-and-out of the laser sheet, as evidenced by the apparent “holes” in the crystalline sample. Most of the particles oscillate in a synchronous fashion (see Supplementary Video 3 for a side view). Since each particle’s oscillation frequency is slightly different due to variations in particle mass, some particles begin to move out-of-phase with the rest of the sample. As the amplitude of vertical oscillation increases, nonlinear interactions can lead to a nucleation of melting in the sample, where energy is rapidly redistributed into the horizontal degrees of freedom. Figures 2g-i show a typical nucleation event in the sample. The redistribution of kinetic energy spreads at the speed of sound in the crystalline lattice, roughly 1 cm/s (Fig. 2j), and is reminiscent of similar phenomena observed in vibrated granular layers Losert et al. (1999).

After a nucleation event, the particles form a highly-energetic, gas-like state. However, hydrodynamic damping causes the kinetic energy to decay and consequently the particles are able to re-crystallize until the vertical oscillation of particles begins again, and the cycle repeats itself for hours without modifying any external parameters (see Fig. S2). Switching between the two states occurs over a wide range of time scales, depending on the environmental conditions and sample size, although it is often quasi-periodic. Figure 2k shows the horizontal kinetic energy per particle, $KE_{xy}$ over a few bistable periods. This switching can also be observed in the temporal evolution of structural order of the sample, as evidenced by the six-fold bond orientational order parameter Brock et al. (1986), $\Psi_{6}$, which is equal to $1$ for a perfect hexagonal crystal and 0 for a completely disordered sample (Fig. 2l).

Numerical Simulations — In order to understand the essential ingredients for the emergence of bistability in our system, we constructed a numerical simulation based on our experiments. The simulation consists of interacting charged particles with vertical and horizontal electrostatic confinement, resulting in a system of conservative forces. We add two other forces on each particle: a hydrodynamic drag force as given by the Epstein law Shukla and Eliasson (2009), which relates the damping experienced by a particle to the ambient pressure, and a spatially-uniform, Langevin force in the vertical direction with relative strength $\beta$. The amplitude of oscillation is determined by a balance between the stochastic forcing, $\beta$, and dissipation due to damping, $\gamma$. As discussed above, the nature of the vertical oscillations is not well-understood. We initially chose to drive the system using periodic forcing in the vertical direction, however, this did not cause bistable switching. This is likely due to the fact that a sinusoidally-driven harmonic oscillator responds at the driving frequency, so that all particles move together, whereas a stochastically-driven harmonic oscillator has a peak response at its fundamental frequency, leading to out-of-phase oscillations between neighboring particles and melting of the crystal (Fig. 2g-i).

The numerical simulations used a custom molecular dynamics code which treated particles as point-like spherical capacitors where the charge is proportional to the particle voltage, $V_{p}$. Particle diameters were chosen from a Gaussian distribution with mean = 10 $\mu$m and coefficient of variation = $c_{v}$. The quenched disorder in the system was due to variation in both mass and charge among the particles. The inter-particle potential was $U(r_{ij})=q_{i}q_{j}e^{-r_{ij}/\lambda_{D}}/(4\pi\epsilon_{0}r_{ij})$, where $q$ is the charge on each particle, and $r_{ij}$ is the distance between the centers of particle $i$ and $j$, and $\lambda_{D}=1$ mm. The vertical confining potential was $V_{v}=-\chi_{v}z^{2}/2$ and the horizontal confining potential was $V_{h}=-\chi_{h}(x^{2}+y^{2})/2$, where $\chi_{v}=4.13\times 10^{6}$ V/m2 and $\chi_{h}=2000$ V/m2. These values were chosen to satisfy the typical oscillation frequencies of particles in the experiment (i.e. 20 Hz for vertical oscillations and 0.5 Hz for horizontal oscillations). The hydrodynamic drag force was $\vec{\bf{F}}_{d}=-\gamma m\vec{\bf v}$, where $m$ is the particle mass and $\gamma$ is the damping coefficient. The stochastic force in the vertical direction used to model the self-induced vertical oscillations was given by $\vec{\bf F}_{s}=\eta(\beta,t)m\sqrt{\Delta t_{0}/\Delta t}$ $\hat{\bf z}$, where $\eta(\beta,t)$ is a random number chosen at each time step from a Gaussian distribution with standard deviation $\beta$, $\Delta t$ is the time step in the simulation, and $\Delta t_{0}=0.001$ s. The particle positions and velocities were advanced in time using velocity-Verlet integration. Given the large parameter space, we used $\gamma$ = 0.2 s*-1*, $\beta$ = 0.5 m/s2, $V_{p}$ = -6 V, and $c_{v}$ = 1.25% as default parameters unless otherwise noted.

There are two system properties which facilitate bistable switching. The first is quenched disorder due to variation of particles sizes, given by $c_{v}$. The second is dynamical noise derived from the vertical stochastic forcing, given by $\beta$. Figure 3a-b shows the time evolution of the ratio of horizontal to total kinetic energy, $\Delta_{xy}=KE_{xy}/(KE_{z}+KE_{xy})$, for two values of $c_{v}$. For $c_{v}=1.25\%$, which is close to the experimental value, we observe bistable switching. For larger values of $c_{v}$, the system remains unstable and the energy is nearly equipartitioned between the various degrees of freedom, which corresponds to $\Delta_{xy}=2/3$. For $c_{v}\lesssim 0.75\%$, the particles oscillate in-phase and the sample remains crystalline. The switching timescales can be extracted from the Fourier spectrum of $\Delta_{xy}$, as shown in Fig. 3c. Bistable switching corresponds to a large Fourier amplitude at low frequencies. If the system was self-organized near a critical point in the dynamics, then one would expect a 1/$f$ behavior in the Fourier spectrum Bak et al. (1987); Kadanoff et al. (1989); Manna (1991). The data shows good agreement with this scaling at higher frequencies, but rolls off at low frequencies.

By varying both $\beta$ and $c_{v}$, we find that the intensity of switching is maximized for intermediate values of both quenched disorder and dynamical noise. Figure 3d shows a heat map of the integrated Fourier amplitude below $f=0.02$ Hz, which we use as a proxy for the intensity of switching. The map verifies that too little noise is insufficient to induce switching between the states, whereas too much noise leaves the system stuck in the gas-like state. However, $\beta$ and $c_{v}$ can not be easily controlled in the experiment. In order to better connect the simulation with the experiment, in Fig. 3e we plot the switching intensity as a function of $V_{P}$ and $\gamma$. The switching occurs for a fairly narrow range of $V_{P}$ and $\gamma$, which is consistent with Fig. 1f. The agreement is good considering that in the experiment, changing $V_{bias}$ also changes the vertical and horizontal confinement. Finally, in Figs. 3f-g, we plot the average vertical oscillation amplitude of a single particle as a function of $\beta$ and $\gamma$, which also shows good agreement with the typical amplitudes measured in the experiment (Fig. 1e).

Summary — The experiments and simulations presented here display a broad class of nonequilibrium phenomena in a single system with minimal ingredients and rich dynamics. We have experimentally demonstrated global bistability in a spatially-extended system composed of non-bistable elements. Given the underlying first-order phase transition between the condensed and gas-like states, our experiment may be a realization of self-organized bistability di Santo et al. (2016). The inter-state switching is facilitated by both quenched disorder and dynamical noise. The time scales of individual stable and unstable periods are not symmetric (see Fig. S3). Durations of instability are mostly determined by the damping time, whereas the stability durations can be much longer and depend on the nucleation of an energy-redistribution event. This is a common property in many excitable systems Lindner et al. (2004), where the relaxation path is more deterministic than the excitation path. However, the distribution of switching timescales in the experiment is narrower than in the simulation. This may be due to a weak periodic signal in the experiment which couples with the noise to induce switching. The source of the periodicity, in addition to controlling the vertical oscillations through modulating the electrode voltage, are subjects of current investigation in our lab.

Acknowledgments — This work was supported by the NSF DMR Grant No. 1455086. We thank E. R. Weeks, K. Wiesenfeld, R. Khomeriki, J. Goree, I. Nemenman, G. H. Hentschel, F. Family, D. Hofmann, B. Beal, and C. Carter for fruitful discussions.