Manipulating shear-induced non-equilibrium transitions in colloidal films by feedback control

Abstract

Using Brownian Dynamics (BD) simulations we investigate non-equilibrium transitions of sheared colloidal films under controlled shear stress σ_xz. In our approach the shear rate [small gamma, Greek, dot above] is a dynamical variable, which relaxes on a time scale τ_c such that the instantaneous, configuration-dependent stress σ_xz(t) approaches a pre-imposed value. Investigating the dynamics under this “feedback-control” scheme we find unique behavior in regions where the flow curve σ_xz([small gamma, Greek, dot above]) of the uncontrolled system is monotonic. However, in non-monotonic regions our method allows to select between dynamical states characterized by different in-plane structure and viscosities. Indeed, the final state strongly depends on τ_c relative to an intrinsic relaxation time of the uncontrolled system. The critical values of τ_c are estimated on the basis of a simple model.

---

1. Introduction

Soft matter under shear flow can display rich nonlinear behaviour involving transitions between different dynamical states,^1,2 spontaneous spatial symmetry-breaking,³ shear-banding^4–7 (for recent reviews, see ref. 8–12), rheochaos^13–15 and heterogeneous local dynamics.^16,17 These intriguing phenomena often strongly affect the rheological properties of the system. Understanding shear-induced effects in, e.g., complex surfactant solutions¹⁸ or liquid crystals,^1,2 colloids,^19–21 soft glasses,^16,17,22 and active fluids,²³ is thus a major current topic connecting non-equilibrium physics and soft material science.

A quantity of particular interest is the flow (or constitutive) curve,^7,12 that is, the shear stress σ as function of the shear rate [small gamma, Greek, dot above], both of which can serve as experimental control parameters. In many examples, the curve σ([small gamma, Greek, dot above]) behaves not only nonlinear (indicating shear-thinning,^24,25 shear-thickening,^25,26 sometimes connected irregular (chaotic) rheological response^15,27), but becomes also multivalued, i.e. different shear rates lead to the same stress. In complex fluids of e.g. wormlike micelles, this multivalued property is in fact a universal indicator of a shear-banding instability, specifically, gradient banding, where the (formerly homogeneous) system separates in gradient direction into coexisting bands characterized by a smaller and a larger local shear rate¹² (note that this is different from the more exotic vorticity banding, i.e., the formation of bands along the vorticity direction as discussed e.g., in ref. 11, 12 and 28). In soft (colloidal) glasses, multivalued functions σ([small gamma, Greek, dot above]) occur as transient phenomena after a sudden switch-on of shear stress (Bauschinger effect),²⁹ or in the vicinity of the so-called yield stress;³⁰ in these systems one observes strong dynamical heterogeneities.¹⁷ A further intriguing feature is that close to such nonmonotonicities, the system's behaviour strongly depends on whether one uses the shear stress or the shear rate as a control parameter.^14,31 In fact, both choices are common in rheological experiments.^22,32,33

Here we present BD computer simulation results of yet another system with multivalued stress–shear curve, that is, a thin colloidal film sheared from an (equilibrium) state within in-plane crystalline order in a planar Couette geometry. As shown in previous experimental³⁴ and simulation^21,35 studies, such films can display successive non-equilibrium transitions from square crystalline over molten into hexagonal states. Here we demonstrate that the structural transitions lead again to non-monotonic flow curves, with a shape reminding that of materials which perform a solid-to-liquid transition beyond a critical (yield) stress.³⁶

Based on this nonlinear behaviour, we then investigate the films in presence of controlled shear stress. In fact, so far most simulations of sheared colloids have been done under fixed shear rate [small gamma, Greek, dot above], exceptions being e.g.ref. 17, 30 and 37, where constant σ has been realized by fixing the force acting on the atoms forming the walls. Here we introduce an alternative, easy-to-apply scheme to control σ which has been previously used by us in continuum approaches of sheared liquid crystals.³¹ In that scheme [small gamma, Greek, dot above] (which directly enters the BD equations of motion) becomes a dynamical variable whose time dependence is governed by a relaxation equation involving a time scale τ_c. The relaxation is based on the difference between the instantaneous (configuration-dependent) stress σ(t) and a preimposed value σ₀. The idea of adapting the shear rate [small gamma, Greek, dot above] is inspired by experiments of stress-controlled systems.^32,38 Our scheme differs from earlier schemes^17,30,37 where the desired value σ₀ is imposed instantaneously. Moreover, due to the coupling to the particle positions our method corresponds to a “feedback” (closed-loop) control scheme, which is similar in spirit to e.g. a Berendsen thermostat for temperature control.³⁹ However, here the choice for τ_c is found to be crucial for the observed dynamical behaviour. In particular we demonstrate that, if our scheme is applied within the multivalued regime of σ([small gamma, Greek, dot above]), the final state strongly depends on the magnitude of τ_c relative to important intrinsic time scales of the system. Thus, the stress-control can be used to deliberately select a desired dynamical state.

2. Model and simulation details

We consider a model colloidal suspension where two particles with distance r_ij interact via a combination of a repulsive Yukawa potential, u_Yuk(r_ij) = W exp[−κr_ij]/r_ij and a repulsive soft-sphere potential u_SS(r_ij) = 4ε_SS(d/r_ij)^1/2 involving the particle diameter d.²¹ The interaction parameters are set in accordance to a real suspension of charged silica macroions [for details, see ref. 40 and 41]. Specifically, at the density considered (see below), the interaction strength W/(k_BTd) ≈ 123 (where k_BT the thermal energy) and the inverse screening length κ ≈ 3.2d. Spatial confinement is modeled by two plane parallel, smooth, uncharged surfaces separated by a distance L_z along the z direction and of infinite extent in the x–y plane. We employ a purely repulsive fluid-wall decaying as z⁻⁹ [see ref. 21]. This is motivated by previous investigations of the equilibrium layer structure, where we found a very good agreement with experiments.^40,41

Our investigations are based on standard BD simulations in three dimensions, where the position of particle i is advanced according to⁴²

where F_i is the total conservative force acting on particle i. Within the framework of BD, the influence of the solvent on each colloidal particle is mimicked by a friction constant and a random Gaussian displacement δW_i. The friction constant is set to (D₀/k_BT)⁻¹, where D₀ is the short-time diffusion coefficient. The value δW_i has zero mean and variance 2D₀δt for each Cartesian component. The time scale of the system was set to τ = σ²/D₀, which defines the so-called Brownian time. We impose a linear shear profile [see last term in eqn (1)] representing flow in x- and gradient in z-direction, characterized by uniform shear rate [small gamma, Greek, dot above]. We note that, despite the application of a linear shear profile, the real, steady-state flow profile can be nonlinear.⁴³ This approximation has also been employed in other recent simulation studies of sheared colloids;^19,44,45 the same holds for the fact that we neglect hydrodynamic interactions.

One quantity of prime interest in our study is the x–z component of the stress tensor,

Thus, we consider only the configuration-dependent contribution to σ_xz; the kinematic contribution (which involves the velocity components in x- and z-direction) is negligible under the highly confined conditions here. We note that, apart from the kinematic contribution, eqn (2) also neglects higher-order contributions involving gradients.^46,47 In continuum approaches, one typically includes non-local terms which are essential for the description of interfaces between shear-bands.^7,14 The importance of such terms in our highly confined system, which is characterized by pronounced layer formation, remains to be investigated.

Based on the shear stress, we introduce a feedback scheme as follows. Starting from an initial value for [small gamma, Greek, dot above] we calculate, in each time step, the configuration-dependent stress σ_xz from eqn (2) and adjust [small gamma, Greek, dot above] via the relaxation equation

where σ₀ is a pre-imposed value of σ_xz, and τ_c determines the time scale of relaxation. Also, η₀ is the shear viscosity obtained for [small gamma, Greek, dot above] → 0 (Newtonian regime). This control scheme is inspired by experiments under fixed stress [see, e.g.ref. 32], where the adaptation of the shear rate to a new stress value always takes a finite time.

From a more formal point of view, we note that through eqn (3), [small gamma, Greek, dot above] becomes an additional dynamical variable. Therefore, and since σ_xy(t) depends on the instantaneous configuration {r_i(t)} of the particles, simultaneous solution of the N + 1 equations of motion (1) and (3) forms a closed feedback loop with global coupling. Interestingly, this feedback scheme is in accordance with the common view that, in a stable system, the shear stress σ_xz should increase with the applied shear rate. This can be shown (at least for a homogeneous system) by a linear stability analysis of eqn (3) as outlined in Appendix A.

In our numerical calculations, we focus on systems at high density, specifically ρd³ = 0.85, and strong confinement, L_z = 2.2d. The corresponding equilibrium system forms a colloidal bilayer with crystalline in-plane structure.⁴¹ We also show some data with L_z = 3.2d, yielding three layers. The values L_z = 2.2 and 3.2 have been chosen because the equilibrium lattice structure is square⁴¹ (other values would lead to hexagonal equilibrium structures which do not show the shear-induced transitions discussed here). The number of particles was set to N = 1058 and 1587, the width of the shear cell to L ≈ 23.8d and 24.2d for L_z = 2.2d and 3.2d, respectively. Periodic boundary conditions were applied in flow (x) and vorticity (y) direction. The time step was set to δt = 10⁻⁵τ where τ = d²/D₀ the time unit.

The system was equilibrated for 10 × 10⁶ steps (i.e. 100τ). Then the shear force was switched on. After the shearing was started the simulation was carried out for an additional period of 100τ. During this time the steady state was reached. Only after these preparations we started to analyze the considered systems.

3. Shear-induced transitions

We start by considering flow curves for systems at constant [small gamma, Greek, dot above]. The functions σ_xz([small gamma, Greek, dot above]) for both, two- and three layer systems are plotted in Fig. 1, where we have included data for the viscosities η = σ_xz/[small gamma, Greek, dot above]. Note that we have defined η via the externally applied rate [small gamma, Greek, dot above] rather than via the effective shear rate within the system (characterized by the average velocity of the layers²¹), which can show significant deviations from [small gamma, Greek, dot above]. As seen from Fig. 1, both the bi- and the trilayer systems are characterized by a nonmonotonic flow curve, accompanied by pronounced shear-thinning. At small shear rates, the systems display linear stress, related to an Newtonian response of the square in-plane lattice structure. An exemplary simulation snapshot for the three-layer system is shown in Fig. 2a [see ref. 21 for corresponding results for the bilayer]. In fact, within the Newtonian regime (square state) the layer velocities are zero, i.e., the particles are “locked”.²¹ In this state the lattice structure persists and the system reacts to the displacement of the particles in elastical manner. Further increase of the shear rate then destroys the square order. In the bilayer, the system then enters a “shear-molten” state characterized by the absence of translational order within the layers (as indicated, e.g., by a non-zero long-time diffusion constant in lateral directions²¹). At the same time, the entire system starts to flow, that is, the layer velocities increase from zero to non-zero values.²¹ From the function σ_xz([small gamma, Greek, dot above]) plotted in Fig. 1a, the appearance of the shear-molten state is indicated by a sudden decrease of σ_xz, implying the onset of shear-thinning. In fact, with the shear-molten regime the slope of the flow curve is negative. For bulk systems, such a behavior implies that the system is mechanically unstable.^48,49 Here we are considering a strongly confined system, where the macroscopic arguments cannot be immediately applied. Nevertheless, it is an interesting question to which extent the flow curve pertains to a true steady state in the regime where the shear rate has values corresponding to shear-molten configurations. Investigating the shear stress as function of strain (see Appendix B) we find that σ_xz assumes indeed a constant value on the time scale of our simulations; however, the relaxation time is extremely long (see next section). We also mention that our observation of shear-molten (long-time) configurations in the regime, where σ_xz decreases with [small gamma, Greek, dot above], is consistent with findings in an earlier theoretical study of a colloidal suspension under shear.⁵⁰

Fig. 1 (Color online) steady state shear stress and shear viscosity (insets) for bi- and trilayer systems as function of the applied shear rate. Regions indicated as I, II, III are discussed in the main text.

Fig. 2 (Color online) snapshots of a colloidal trilayer at three different shear rates corresponding to (a) square, (b) laned, and (c) hexagonal state. In the main figures, the three colours correspond to particles in the three layers. In the insets, the two colours indicate the particles of the middle layer which are positioned at z > 0 (red) and z < 0 (yellow), respectively (where z = 0 is the middle plane).

Somewhat different behaviour is found in the trilayer system which displays, before melting, an intermediate state [see Fig. 2b]: this state is characterized by a non-zero layer velocity. In addition, the middle layer separates into two sublayers, in which the particles are ordered in “lanes” [see inset of Fig. 2b] and move with the velocity of the corresponding outer layer (a more detailed discussion of this “laned” state will be given elsewhere⁵¹). Only further increase of [small gamma, Greek, dot above] then yields a shear-molten state characterized by a decreasing flow curve (in analogy to the bilayer). Finally, both systems transform into a state with in-plane hexagonal ordering [see Fig. 2c] and low viscosity. In this hexagonal state the layer velocity is non-zero,²¹ that is, the systems flows. As demonstrated earlier by us²¹ the mechanism of relative motion involves collective oscillations of the particles around lattice sites, consistent with recent experiments of 3D sheared colloidal crystals.²⁰ Regarding the stress, we see that the hexagonal regime is (in both systems) characterized by a slight increase of σ_xz with [small gamma, Greek, dot above]. As a consequence, there is a parameter range (indicated as region “II” in Fig. 1) where the flow curve is multivalued, that is, different [small gamma, Greek, dot above] lead to the same σ_xz. In many contexts (such as for worm-like micelles), multivalued flow curves are associated with the phenomenon of shear-banding, that is, the separation of the system into spatial regions with different shear rates. In our case, where the system consists of two or three layers such a separation can not occur. Instead, the non-monotonic stress curve is a consequence of the structural transitions of the system induced by the shear.

4. Intrinsic time scales

Before exploring the impact of shear–stress control, which involves a time scale itself through the parameter τ_c [see eqn (3)], we take a closer look at the intrinsic time scales characterizing the uncontrolled systems. We focus on the bilayer (the same findings apply qualitatively on the trilayer) and consider the response of the unsheared equilibrium system, which is in a square state, to a sudden switch-on at time τ_on of shear with rate [small gamma, Greek, dot above]_new. The resulting time dependence of the instantaneous stress is plotted in Fig. 3.

Fig. 3 (Color online) response of σ_xz(t) to a sudden switch-on (at time τ_on) of shear with different rates [small gamma, Greek, dot above]_newτ. The simulations were started from the equilibrium (square) state in a bilayer system. The inset shows the fit of the relaxation times τ₁ according to eqn (4).

If [small gamma, Greek, dot above]_newτ has a value pertaining to the square state, the shear stress jumps at τ_on to a non-zero value but then settles quickly to its steady-state value [see Fig. 1]. At shear rates corresponding to the shear-molten state we can also observe a relaxation at some non-zero value. However, it should be emphasized that this value is transient in character. The true, steady state value is only achieved at much longer simulation times (see the stress–strain relations presented in Appendix B). Finally, for shear rates related to the hexagonal steady state ([small gamma, Greek, dot above]_newτ > [small gamma, Greek, dot above]_hexτ ≈ 257, [see Fig. 1]), we observe a well-pronounced stress overshoot, similar to what is observed e.g. in soft glassy systems,²⁹ wormlike micelles⁵² and polymer melts.⁵³ Closer inspection shows that the actual value of τ₁ as well as the functional behavior of σ_xz(t) strongly depends on the distance between [small gamma, Greek, dot above]_newτ and the threshold between shear-molten and hexagonal state, [small gamma, Greek, dot above]_hexτ: the smaller this distance is, the larger becomes τ₁, and the more sensitive it is against small changes of the shear rate. Moreover, a sudden quench deep into the hexagonal state leads to an oscillatory relaxation of the stress σ_xz(t) [see curves [small gamma, Greek, dot above]_newτ = 400, 500], with τ₁ (which now corresponds to the relaxation time of the envelope) being still quite large. Taken together, for [small gamma, Greek, dot above]_newτ > [small gamma, Greek, dot above]_hexτ, τ₁ can be fitted according to (see inset in Fig. 3)

where we find a₁/τ = 0.21, b₁ = 0.52 (setting [small gamma, Greek, dot above]₁τ = [small gamma, Greek, dot above]_hexτ). The oscillations occurring at large [small gamma, Greek, dot above]_new induce yet a different time scale τ_os, which is smaller than τ₁. Specifically, we find τ_os/τ = [scr O, script letter O](10⁻³). The physical reason for these oscillations is the “zig-zag” motion of particles in adjacent layers.²¹ This motion is accompanied by periodic variations of nearest-neighbor distances, and thus, pair forces, which eventually leads to oscillations of σ_xz(t).

Furthermore it is interesting to relate the relaxation times emerging from Fig. 3 to the structural transition from square to hexagonal state in Fig. 1. To this end we consider the Peclet number Pe ≡ [small gamma, Greek, dot above]τ_Pe where τ_Pe is a “typical” relaxation time.⁵⁴ Identifying τ_Pe with τ₁ and considering shear rates [small gamma, Greek, dot above] close to the transition from the quadratic into the shear-molten state, we find Pe = [scr O, script letter O](10⁰). In other words, the shear-induced structural transitions happen at Pe ≧ 1, consistent with our expectations.⁵⁴

For comparison we have also investigated the reverse situation, where the system is initially in a hexagonal steady state at shear rate [small gamma, Greek, dot above]_init = 400. We then suddenly change the shear rate towards a much smaller value and consider the relaxation towards the square equilibrium state. The corresponding behaviour of σ_xz(t) is shown in Fig. 4. Again we find that, the smaller the difference [small gamma, Greek, dot above]_newτ − [small gamma, Greek, dot above]₂τ is, the larger τ₂ becomes (and the more pronounced is the sensitivity to small changes in [small gamma, Greek, dot above]_new). The resulting relaxation times can be also fitted viaeqn (4) with a₂/τ = 22.58, b₂ = 1.73 and [small gamma, Greek, dot above]₂τ = 215, whereby [small gamma, Greek, dot above]₂τ = [small gamma, Greek, dot above]_sqτ denotes the threshold between the square and the molten states. The result for this fit is visualized in the inset of the Fig. 4.

Fig. 4 (Color online) response of σ_xz(t) to a sudden change (at time τ_on) of shear. The new shear rates [small gamma, Greek, dot above]_newτ result in a relaxation in the square state. The simulations were started from the hexagonal steady state at [small gamma, Greek, dot above]τ = 400 in a bilayer system. The inset shows the fit of the relaxation times τ₂ according to eqn (4).

5. Impact of feedback control

We now discuss the impact of our shear stress control scheme defined in eqn (3). The latter involves the zero-shear viscosity, η₀, which is estimated from Fig. 1 to η₀ = 0.086/dτ and 0.323/dτ for the bilayer and trilayer, respectively.

The overall dynamical behaviour under feedback control strongly depends on the value of σ₀ (imposed shear stress) relative to the flow curve of the original system [see Fig. 1]. We can differentiate between regimes I, II, and III, which are indicated in Fig. 1.

For a σ₀ chosen in region I, the response of the system is unique, that is, the final state is independent of the control time scale τ_c, as well as of the initial shear rate [small gamma, Greek, dot above]_init and the initial microstructure. Thus, when starting from a square state, with a corresponding initial shear rate [small gamma, Greek, dot above]_init, the system immediately settles at this state. As a more critical test of the injectivity of the flow curve in region I, we plot in Fig. 5a and b the functions [small gamma, Greek, dot above](t) and σ_xz(t) for the bilayer system at σ₀dτ² = 6 and various τ_c, starting from a hexagonal configuration (and [small gamma, Greek, dot above]_initτ = 400). In all cases, the shear rate decreases towards the value [small gamma, Greek, dot above]τ ≈ 70 and the structure relaxes into the square state pertaining to the value σ_xzdτ² = 6 in the uncontrolled system. This indicates that the square state in region I is indeed the only fixed point of the dynamics. We also see from Fig. 5a that the relaxation time into this steady state increases with τ_c. Fig. 5b additionally shows that σ_xz(t) displays a pronounced peak. The peak indicates the time window in which the initial hexagonal ordering transforms into a square one. In fact, the high values of σ_xz at the peak reflect the large friction characterizing the intermediate molten state. Similar behaviour occurs in region I of the trilayer system [see Fig. 5c and d] where, however, fluctuations of σ_xz(t) are generally larger.

Fig. 5 (Color online) time dependence of the instantaneous shear rate and shear stress for a bilayer- [a and b] and a trilayer system [c and d] in presence of feedback control within region I. The imposed stress was set to σ₀dτ² = 6(2) for the bilayer (trilayer) system. Various values of τ_c/τ are considered. The initial configuration is hexagonal.

We now choose σ₀ within region II of the flow curve, where there are three different shear rates (and thus, three fixed points) pertaining to the same stress [see Fig. 1]. We focus on systems which are initially in a square configurations, whereas the initial shear rate [small gamma, Greek, dot above]_init has a value pertaining to the hexagonal state (other initial conditions will be discussed below). The impact of τ_c on the time dependence of [small gamma, Greek, dot above](t) and σ_xz(t) is shown in Fig. 6. For small values of the control time scale the systems stays in the initial lattice configuration, i.e., [small gamma, Greek, dot above] relaxes towards the value pertaining to the square state ([small gamma, Greek, dot above]τ ≈ 90). Different behaviour occurs at larger values of τ_c/τ: although the initial structure is square, the final state is hexagonal, and the shear rate essentially remains at its high initial value. We stress that these findings crucially depend on the choice of [small gamma, Greek, dot above]_init. In particular, the dependency of the long-time behaviour on τ_c/τ only arises for large values of [small gamma, Greek, dot above]_init; for small values the system remains in the square state irrespective of τ_c. An overview of the final dynamical states in the feedback-controlled bilayer at σ₀dτ² = 8 and various combinations of [small gamma, Greek, dot above]_init and τ_c/τ (assuming a square initial structure) is given in Fig. 7. The colour code indicates the ratio of local bond-order parameters 〈Ψ₆/Ψ₄〉 [for a definition of the Ψ_n see e.g., ref. 21]. The restriction to values 〈Ψ₆/Ψ₄〉 ≤ 6 is related to the actual values observed in the simulations. From Fig. 7 one clearly sees that for initial shear rates [small gamma, Greek, dot above]_initτ > [small gamma, Greek, dot above]_hexτ ≈ 257, the final state of the feedback-controlled system depends on τ_c/τ. This is in contrast to the uncontrolled system which becomes hexagonal for all [small gamma, Greek, dot above]_init > [small gamma, Greek, dot above]_hex. For a hexagonal initial configuration the diagram (not shown here) looks similar from a qualitative point of view; however, the range of control times where the system retains a hexagonal state despite of [small gamma, Greek, dot above]_init < [small gamma, Greek, dot above]_quad (with [small gamma, Greek, dot above]_quad being the threshold between square/molten states) is much smaller.

Fig. 6 (Color online) same as Fig. 5, but for σ₀dτ² = 8(2.7) for the bilayer (trilayer) system (region II). The initial configuration is square.

Fig. 7 (Color online) state diagram indicating long-time lattice structures. All simulations were started from a square initial structure and the imposed shear stress was set to σ₀dτ² = 8. The line shows the result from eqn (8).

We conclude that, by varying τ_c and the initial structure, we can “switch” between the two stable, steady-state configurations arising in the multivalued region of the uncontrolled system. That these states are stable also under feedback (stress) control is supported by the linear stability analysis presented in Appendix A. Indeed, the dynamics under feedback control never evolves towards the intermediate, shear-molten states, consistent with the view that these states are mechanically unstable. This holds also in region III of the flow curve of the uncontrolled system, e.g., for σ₀dτ² = 16(5) for the bilayer (trilayer): here, a small value of τ_c yields relaxation towards the square state, whereas for large τ_c, the system evolves into a hexagonal state. Finally, we note that completely analogous behaviour is found in the trilayer system [see Fig. 6b and c] for a σ₀ pertaining to the regime where square, molten and hexagonal states exist.

6. Transition line

The most significant observation from Fig. 7 is that at high values of [small gamma, Greek, dot above]_init, the feedback-controlled system can achieve either the hexagonal or the square configuration, provided that we start from a square configuration and choose τ_c/τ accordingly. We now propose a simple model which allows us to estimate the transition values of the control time, τ^trans_c.

The physical idea behind our model is that, with the initial conditions described above, relaxation into the hexagonal state only occurs if the reorganization time τ_reorg required by the system to transform from a square into a hexagonal configuration, is smaller than the time τ_decay in which [small gamma, Greek, dot above] decays to a value pertaining to the square state. We can estimate τ_decay from eqn (3) if we assume, for simplicity, a linear relationship σ_xz(t) = m[small gamma, Greek, dot above](t) (note that such a relationship is indeed nearly fulfilled within the square and hexagonal states, see Fig. 1). Under this assumption eqn (3) can be easily solved, yielding

[small gamma, Greek, dot above](t) = m⁻¹e^{−mt/η₀τ_c}(m[small gamma, Greek, dot above]_init − σ₀ + e^mt/η₀τ_cσ₀). (5)

From eqn (5) we find that the decay time of [small gamma, Greek, dot above] to the threshold value [small gamma, Greek, dot above]_hexτ ≈ 257 (below which the hexagonal state of the uncontrolled system is unstable) is given by

To estimate the reorganization time τ_reorg (from the initial square into a hexagonal configuration), we assume that its dependence on [small gamma, Greek, dot above]_init is analogous to that of the relaxation time τ₁ introduced for the uncontrolled system [see eqn (4)]. Specifically, we make the ansatz

As stated above, a crucial assumption of our model is that the system can only reach the hexagonal state if τ_reorg does not exceed τ_decay. Note that the latter involves (in fact, is proportional to) the time τ_c. By equating expressions (6) and (7) for τ_decay and τ_reorg, respectively, we can therefore find an expression for the minimal control time, τ^trans_c, above which the system reaches the hexagonal state, that is

Due to the square initial configuration, we set m = η₀ and σ_xz(t) = m[small gamma, Greek, dot above](t) as defined in our ansatz. The remaining parameters a′ and b′ are determined by fitting the numerical results for τ_c/τ at the boundary [see Fig. 7] to expression (8), yielding a′/τ = 54.127 and b′ = 1.503. The resulting line τ^trans_c([small gamma, Greek, dot above]_init) is included in Fig. 7, showing that our estimate describes the transition between square and hexagonal states very well.

Similar considerations are possible, when we use a hexagonal initial lattice structure. Choosing then a small value of γ_initτ we find that we can switch between hexagonal and square state. This is illustrated in Fig. 8. To obtain the corresponding transition values of τ_c, we use the same strategy as before, but take a different ansatz for the stress. Specially, we set σ_xz(t) = n + m[small gamma, Greek, dot above](t) which approximately describe the flow curve in the hexagonal state of the uncontrolled system. From the results plotted in Fig. 1 we find n = 7.0477/dτ² and m = 0.0025. The analog of the eqn (8) then reads

with a′/τ = 0.012 and b′ = 0.237 [see Fig. 8]. The result is visualized in Fig. 8. Comparing the typical control time scales at the transition with those seen in Fig. 7 we find that τ_c/τ which is necessary to switch from square into the hexagonal state [see Fig. 7] is about two orders of magnitude larger than switching from hexagonal into the square state. We suspect that this difference results from the differences of the slope of the shear stress in the square and hexagonal regimes.

Fig. 8 (Color online) state diagram indicating long-time lattice structures. All simulations were started from the hexagonal initial structure and the imposed shear stress was set to σ₀dτ² = 8.

7. Conclusions

Using numerical simulation we have studied the complex dynamical behaviour of sheared colloidal films under a specific type of shear–stress control. Our approach involves relaxation of the shear rate in a finite relaxation time τ_c, until the instantaneous stress matches its desired value. This approach is inspired by rheological experiments³² where the instantaneous shear rate as function of time can be measured. Focusing on systems with multivalued flow curves (resulting from successive non-equilibrium transitions) we have found that, by tuning τ_c and the initial conditions, it is possible to select a specific dynamical state. In the present system these are either a state with square in-plane ordering and high viscosity, or a hexagonal state with low viscosity. Therefore, our study suggests a way to stabilize states with desired rheological properties, particularly shear viscosities. Moreover, we have proposed a model which relates the transition values of τ_c to relevant intrinsic relaxation times under sudden change of [small gamma, Greek, dot above].

Although most of our results pertain to a colloidal bilayer, the fact that we found analogous results for trilayers suggests that the proposed technique can also be applied for systems with larger number of layers. In fact, we think that this method, after some minor adaptations such as consideration of the kinematic (and, possibly, also the non-local) contributions in eqn (2), should also be applicable and fruitful in bulk systems. Indeed, we expect the method to allow for state selection in any shear-driven system with multivalued flow curve. For example, in an earlier study we have used an analogous approach (based, however, on continuum equations) to select states and even suppress chaos in shear-driven nematic liquid crystals.³¹ It therefore seems safe to assume that the capabilities of the present scheme are quite wide. For colloidal layers one may envision, e.g., stabilization of time-dependent structures such as oscillatory density excitation, which may have profound implications for lubrication properties.⁵⁵

Finally, our findings are quite intriguing in the broader context of manipulating nonlinear systems by feedback control. In our case, the feedback character stems from the fact that the stress control involves the configuration-dependent instantaneous stress. Mathematically, this scheme can be viewed as feedback control with exponentially distributed time-delay⁵⁶ (as can be seen by formally integrating eqn (3) and inserting it into eqn (1)). Similar schemes are used to stabilize dynamical patterns in laser networks,⁵⁷ neural systems,⁵⁸ and more generally, coupled oscillator systems.⁵⁹ The implications of these connections are yet to be explored.

Appendix

A. Stability of the feedback controlled system

In this Appendix we investigate the stability of the solutions of eqn (3). Specifically, we consider the impact of small variations of the shear rate from its steady state value [small gamma, Greek, dot above]₀ related to the imposed stress σ₀. Expanding the right side of eqn (3) with respect to the difference [small gamma, Greek, dot above] − [small gamma, Greek, dot above]₀ yields

For long times we expect the first term on the right side of eqn (10) to vanish, since σ_xz([small gamma, Greek, dot above]₀,t) → σ₀. To linear order, eqn (10) then reduces to

Noting that the values of τ_c and η₀ are both positive, we can follow that the feedback controlled shear rate approaches a steady-state value only if

This corresponds to the usual criterion of mechanical stability.⁴⁸

B. Strain–stress relation under constant shear rate

In this Appendix we present results for stress–strain relations at different fixed values of [small gamma, Greek, dot above]. These can be obtained from the data shown in Fig. 3 by rescaling the time axis with the applied shear rate. Numerical results are shown in Fig. 9. Similar to the stress–time relations shown in Fig. 3, one observes simple, monotonic behavior for the case [small gamma, Greek, dot above]_newτ = 200 (quadratic regime), while the curves for larger shear rates display pronounced stress overshoots. The width of these overshoots is largest at [small gamma, Greek, dot above]_newτ = 240, where the system is in the molten state. This is consistent with the appearance of a particularly large intrinsic relaxation time as discussed in Section 4.

Fig. 9 (Color online) stress–strain relations in the colloidal bilayer for different shear rates, starting from the equilibrium (square) configuration.

Acknowledgements

This work was supported by the Deutsche Forschungsgemeinschaft through SFB 910 (project B2).

References

1. M. Ripoll, P. Holmqvist, R. G. Winkler, G. Gompper, J. K. G. Dhont and M. P. Lettinga, Phys. Rev. Lett., 2008, 101, 168302.
2. D. Strehober, H. Engel and S. H. L. Klapp, Phys. Rev. Lett., 2013, 88, 012505.
3. M. Das, B. Chakrabarti, C. Dasgupta, S. Ramaswamy and A. Sood, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 021707.
4. R. Moorcroft and S. M. Fielding, Phys. Rev. Lett., 2013, 110, 086001.
5. R. Besseling, L. Isa, P. Ballesta, G. Petekidis, M. E. Cates and W. C. K. Poon, Phys. Rev. Lett., 2010, 105, 268301.
6. R. L. Moorcroft, M. E. Cates and S. M. Fielding, Phys. Rev. Lett., 2011, 106, 055502.
7. C.-Y. D. Lu, P. D. Olmsted and R. C. Ball, Phys. Rev. Lett., 2000, 84, 642.
8. P. T. Callaghan, Rheol. Acta, 2008, 47, 243–255.
9. J. K. G. Dhont and W. J. Briels, Rheol. Acta, 2008, 47, 257–281.
10. S. Manneville, Rheol. Acta, 2008, 47, 301–318.
11. P. D. Olmsted, Rheol. Acta, 2008, 47, 283–300.
12. S. M. Fielding, Soft Matter, 2007, 3, 1262.
13. A. Aradian and M. E. Cates, Europhys. Lett., 2005, 70, 397.
14. A. Aradian and M. E. Cates, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 041508.
15. D. Lootens, H. V. Damme and P. Hebraud, Phys. Rev. Lett., 2003, 90, 178301.
16. J. Zausch and J. Horbach, Europhys. Lett., 2009, 88, 60001.
17. P. Chaudhuri and J. Horbach, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 040301.
18. M. E. Cates and S. M. Fielding, Adv. Phys., 2006, 55, 799.
19. T. H. Besseling, M. Hermes, A. Fortini, M. Dijkstra, A. Imhof and A. van Blaaderen, Soft Matter, 2012, 8, 6931.
20. D. Derks, Y. L. Wu, A. v. Blaaderen and A. Imhof, Soft Matter, 2009, 5, 1060.
21. T. A. Vezirov and S. H. L. Klapp, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 052307.
22. J. S. Raynaud, P. Moucheront, J. C. Baudez, F. Bertrand, J. P. Guilbaud and P. Coussot, J. Rheol., 2002, 46, 709.
23. M. E. Cates, S. M. Fielding, D. Marenduzzo, E. Orlandini and J. M. Yeomans, Phys. Rev. Lett., 2008, 101, 068102.
24. G. L. Burrell, N. F. Dunlop and F. Separovic, Soft Matter, 2010, 6, 2080.
25. X. Cheng, J. H. McCoy, J. N. Israelachvili and I. Cohen, Science, 2011, 333, 6047.
26. A. Fall, A. Lematre, F. Bertrand, D. Bonn and G. Ovarlez, Phys. Rev. Lett., 2010, 105, 268303.
27. R. Bandyopadhyay and A. K. Sood, Europhys. Lett., 2001, 56, 447.
28. K. Kang, M. P. Lettinga and J. K. G. Dhont, Rheol. Acta, 2008, 47, 499–508.
29. F. Frahsa, A. K. Bhattacharjee, J. Horbach, M. Fuchs and T. Voigtmann, J. Chem. Phys., 2013, 138, 12A513.
30. F. Varnik, L. Bocquet and J. L. Barrat, J. Chem. Phys., 2004, 120, 2788.
31. S. H. L. Klapp and S. Hess, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 051711.
32. Y. T. Hu, P. Boltenhagen and D. J. Pine, J. Rheol., 1998, 42, 1185.
33. V. S. Rudraraju and C. M. Wyandt, Int. J. Pharm., 2005, 292, 53–61.
34. R. Biehl and T. Palberg, Europhys. Lett., 2004, 66, 291–295.
35. R. Messina and H. Löwen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 011405.
36. P. H. T. Uhlherra, J. Guoa, C. Tiua, X.-M. Zhanga, J. Z.-Q. Zhoua and T.-N. Fang, J. Non-Newtonian Fluid Mech., 2005, 125, 101–119.
37. V. Mansard, A. Colin, P. Chaudhuri and L. Bocquet, Soft Matter, 2013, 9, 7489.
38. J. Paredes, N. Shahidzadeh-Bonn and D. Bonn, J. Phys.: Condens. Matter, 2011, 23, 284116.
39. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684.
40. S. H. L. Klapp, Y. Zeng, D. Qu and R. von Klitzing, Phys. Rev. Lett., 2008, 100, 118303.
41. S. Grandner and S. H. L. Klapp, J. Chem. Phys., 2008, 129, 244703.
42. D. L. Ermak, J. Chem. Phys., 1975, 62, 4189.
43. J. Delhommelle, J. Chem. Phys., 2003, 119, 21.
44. B. Lander, U. Seifert and T. Speck, J. Chem. Phys., 2013, 138, 224907.
45. J. Rottler and D. J. Srolovitz, Phys. Rev. Lett., 2007, 98, 175503.
46. J. H. Irving and J. G. Kirkwood, J. Chem. Phys., 1950, 18, 817.
47. B. D. Todd, D. J. Evans and R. J. Daivis, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 1627.
48. R. I. Tanner, Engineering Rheology, Oxford Engineering Science Series, 2002, vol. 52, p. 493.
49. S. M. Fielding and P. D. Olmsted, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 036313.
50. P. Harrowell, Phys. Rev. A, 1990, 42, 3421–3431.
51. T. Vezirov, S. Gerloff and S. H. L. Klapp, to be published.
52. J. T. Padding, E. S. Boek and W. J. Briels, J. Phys.: Condens. Matter, 2005, 17, S3347–S3353.
53. J. M. Delay and R. G. Larson, Structure and Rheology of Molten Polymers, Hanser, Munich, 2006.
54. M. J. Stevens and M. O. Robbins, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 48, 5.
55. A. Vanossi, N. Manini and E. Tosatti, PNAS, 2012, 109, 16429.
56. P. Hoevel and E. Schöll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 046203.
57. T. Juengling, A. Gjurchinovski and V. Urumov, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 046213.
58. Y. Song, Y. Han and Y. Peng, Neurocomputing, 2013, 121, 442.
59. Y. N. Kyrychko, K. B. Blyuss and E. Schöll, Eur. Phys. J. B, 2011, 84, 307.

---