The simplest microscopic model of a complex fluid: flow phenomena and constitutive relation

# The simplest microscopic model of a complex fluid: flow phenomena and constitutive relation

R. M. L. Evans111corresponding author    Craig A. Hall School of Mathematics, University of Leeds, LS2 9JT, United Kingdom    R. Aditi Simha IIT Madras    Tom Welsh
August 2016
###### Abstract

It was shown in [PRL 114, 138301 (2015)] that a remarkably simple dynamical model exhibits many of the complex flow regimes and non-equilibrium phase transitions characteristic of complex fluids. By removing extraneous detail, this simplest microscopic model of non-Newtonian flow can reveal the universal physics relevant to all complex fluids. Here we present more detailed results and a full derivation of the model’s compact mean-field constitutive relation, with great potential scope for insights into universality and tractable mathematics. By enforcing local conservation of angular momentum, the one-dimensional (1D) XY-model (originally used for equilibrium magnetic systems) can be driven into various flow regimes, including simple Newtonian behaviour, shear banding, solid-liquid coexistence and slip-plane motion. The model demonstrates that the phenomenon of shear banding does not rely on details of tensorial stress fields, but can exist in 1D.

###### pacs:
75.10.Hk, 83.60.Rs, 83.10.Gr, 82.70.-y, 64.60.Ht

## I Introduction

A detailed and elaborate model of a particular physical systems cannot be used to investigate the physics underlying universal phenomena, because it would not be clear which of the model’s features are crucial to the results. Only by removing a particular feature from a model can one ascertain whether it was important. Appreciating this, the pioneers of condensed matter theory investigated universality in equilibrium phase transitions using some archetypal models, including the Ising model and XY model XY0 ().

Experiments and simulations also hint at some universal phenomenology far from equilibrium, in the flow-induced transitions of non-Newtonian fluids. Transitions between Newtonian flow, shear-banding and slip-planes/fractures can be observed in systems as diverse as foams foamCoex (); Debregeas01 (), dense colloids colloidBanding1 (); colloidBanding2 (); Haw04 (), surfactant solutions Schmitt94 (); surfBanding1 (); surfBanding2 (); Jones95 () and polymers polyBanding1 (); polyBanding2 (); Tapadia03 (); Pouyan09 (). In Ref. PRL (), we demonstrated that those same phenomena exist in the one-dimensional (1D) classical XY model (also known as the classical rotor model or model) when it is given the simplest of angular-momentum-conserving Langevin dynamics. Here we provide more details of the model’s steady-state behaviour, and the full derivation of its remarkably compact and useful constitutive relation.

As depicted in Fig. 1 the 1D XY model consists of a chain of lattice sites, each with one degree of freedom, characterized by the orientation of a two-component unit vector , which was originally invented to represent the constant-magnitude magnetisation of an atom with a continuously orientable magnetic moment constrained to lie in the -plane.

Each site couples only to its two nearest neighbours, via an interaction potential that favours parallel alignment, so that the model’s Hamiltonian (in units of the uniform coupling constant) is given by

 H=N∑j=1[−sj⋅sj−1+12˙θ2j] (1)

where the rotor’s moment of inertia (scaling the final, kinetic energy term) has been set to unity without loss of generality. We may impose periodic boundary conditions by defining , or free boundaries by , or other possibilities, discussed further below. As with any 1D system with short-range interactions, the 1D XY model has trivial equilibrium phase behaviour, with a single transition to an ordered state at zero temperature Mattis84 ().

Away from equilibrium, we have found that this simple one-dimensional model has an elaborate steady-state phase diagram with transitions between different flowing phases. By choosing to represent the -plane perpendicular to the direction of the 1D lattice, we can compare the model’s local angular velocity with the rotation flow profile of a fluid in a parallel-plate rheometer rheometry1 (); rheometry2 (). Of course, our model is not truly three-dimensional like the fluid in a rheometric experiment, but it has many of the same features, including microscopic internal interactions, a macroscopic steady-state flow profile and shear rate, energy exchange, conserved angular momentum and torque propagation.

Equilibrium phase behaviour is independent of the model’s dynamics, so long as a Hamiltonian is defined and detailed balance is respected, but non-equilibrium states depend on the specific dynamics. We use the simplest physically motivated dynamics for a thermal system: Langevin dynamics, with all forces (torques) respecting Newton’s (rotational) third law of motion. Hence, all forces (conservative, dissipative and stochastic) are equal and opposite on nearest neighbours, and and angular momentum is conserved. The equations of motion are

 ¨θj = τj−τj−1, (2a) τj = sinΔθj+μΔ˙θj+ηj(t) (2b)

where is the torque between rotors and , which have relative angle . The coefficient of friction is , and are the usual delta-correlated stochastic functions at temperature with

 ⟨ηi(t)ηj(t′)⟩=2μTδ(t−t′)δij (3)

and zero mean. Hamilton’s equations are recovered in the conservative case . Note that no centrifugal or corriolis forces are in effect, because there is no radial component of motion, so the system is independent of absolute angular velocity and respects an angular version of Galilean relativity.

We can nevertheless drive the system into a non-equilibrium steady state by counter-rotating its boundaries to induce shear flow. Thus, non-zero mean work is done by the boundary, but no body-forces are applied, so the Hamiltonian remains that of the XY model.

We have conducted both numerical and mathematical investigations of the XY model’s steady-state behaviour in shear flow. In the next section, we define the numerical protocol, including the physically appropriate boundary conditions for imposing torsional shear flow in a simple manner. The results, including the steady-state phase diagram of the XY model, are presented in section III. A theoretical analysis is given in section IV, and conclusions and open questions are discussed in section V.

## Ii Numerical method

### ii.1 Algorithm

The Dissipative Particle Dynamics (DPD) algorithm DPD () was used to evolve the state of the system. It shares the stability of the velocity Verlet algorithm, but can be applied to situations with velocity-dependent forces, such as ours, and includes the random thermal forces of Langevin dynamics.

In order to eliminate edge effects, all rotors were treated equally by the simulation algorithm, through the application of periodic boundary conditions. It is perhaps surprising that boundary conditions exist that are periodic but also apply a relative torque across the system. The appropriate condition is an angular version of the Lees-Edwards “sliding brick” boundary condition LeesEdwards (), whereby the rotors at opposite ends of the chain are designated as neighbours, but see each other through an angular offset that increases at a constant angular velocity. This can be done without altering Eqs. 2b simply by defining

 τ0≡τN (4)

and

 ΔθN≡θ0−θN+N˙γt (5)

so that

 Δ˙θN=˙θ0−˙θN+N˙γ. (6)

Hence the -rotor system occupies a periodic space which has a topological twist that increases at a rate . That twist need not be localized between rotors and , but may be distributed across the system, while the boundary rotors may choose their velocities such that their interaction torque fluctuates in a statistically identical way to the other inter-rotor torques. Thus all rotors are governed by the same equations of motion, while the system as a whole experiences a shear rate (angular velocity difference per rotor) of .

The distribution of the random torques was given a Gaussian form with variance

 σ2=2μTdt (7)

and zero mean and a discrete time step .

### ii.2 Calibration

The time step was chosen adaptively to ensure that the interaction potential was explored in sufficient detail, irrespective of rotor speed. On each time step, the maximum increment in angular difference between any interacting pair of rotors was limited to a maximum magnitude . Thus, no pair of rotors could pass over the potential barrier without feeling the appropriate torque. If a pair of interacting rotors had a sufficiently high relative angular velocity that their relative angular increment exceeded the specified limit, all data for that time-step was discarded and the step was repeated with a reduced time interval . To maintain simulation speed, the time-step for the next iteration was increased whenever the maximum relative angular increment between all neighbours was found to be substantially below the threshold. For the majority of states, results were found to be reliable for a threshold as high as 1 radian. This becomes less surprising when one considers that the vast majority of relative angular increments within the system are much smaller than the threshold value. Nevertheless, we chose never to set above , and used lower values in some cases (often ).

To ensure that no discretization artefacts, start-up transients or finite-size effects were present in the final results, trial simulations were run in a variety of different regions of the parameter space, and repeated for different values of and system size , and different initial conditions. A system of rotors was found to be sufficiently large to eliminate finite-size effects from all states except for those with very widely separated slip-planes (see Results, section III). Hence a safe value of was used for most simulations, occasionally increasing to 1024 or 2048 for extreme cases.

In experiments and simulations, both in and out of equilibrium, it is never straightforward to determine the decay time of all start-up transients in order categorically to establish a steady state. For instance, equilibrium investigations can be misled by long-lived metastable states that are impossible to distinguish from the true steady state except by infinite patience. (Without evidence of the existence of graphite, diamond would be mistaken for the ground-state of carbon.) We must therefore accept that our study of the XY model’s steady states cannot be definitive. Subject to that proviso, we judge that initial transients have decayed when simulations with very different initial conditions converge to a statistically similar state of motion. The initial conditions tried included a state of uniform shear rate () and states in which plateaux of uniform velocity (where ) are separated by sharp discontinuities (“slip planes”). The convergence time varies greatly with the values of the parameters (), becoming very large close to phase transitions. In most cases the initially uniform system arrived most quickly at its final steady state.

Having established our protocols, we performed simulations at a large number of points in parameter space in order to chart the system’s non-equilibrium steady-state phase behaviour.

## Iii Results

Four distinct types of flow behaviour where found, separated by phase transitions or narrow cross-over regions: a state of uniform time-averaged shear flow, a shear-banding state, a coexistence between solid and fluid regions, and a state of solid regions separated by localized slip planes. Examples of each are shown in Fig. 2, in steady states that were found to be well established after . Animations of some of the steady states are available online animations ().

On increasing temperature or friction coefficient or shear rate, the general trend is to progress from the slip-plane regime, through solid-fluid coexistence and then the shear-banding regime, to a state of statistically uniform shear flow. However, for some parameter values, one or both of the two intermediate states (solid-fluid coexistence and shear-banding) are bypassed.

While a wide variety of relative angles, velocities and accelerations exist at any instant within the uniform phase, it is statistically uniform, in the sense that is independent of , where indicates a time-average. The uniform limit at large parameter values is easily explained; all structure due to interactions is washed out when either thermal energy or shear-induced kinetic energy of relative motion swamps the scale of the interaction potential, or when frictional forces dominate over the potential gradient. This argument, for high shear rate, will be made more rigorous in section IV.2, and a full theoretical investigation of the phase behaviour will be presented in section IV.

### iii.1 Slip planes

The existence of slip planes takes more consideration. Some velocity profiles in the slip-plane state are shown in Fig. 3 for parameter values and and various shear rates .

In the slip-plane regime, the majority of rotors have a time-averaged angular velocity equal to that of their neighbours. So the inter-rotor angle exhibits only small-amplitude excursions about a constant value, within the potential well. These rotor-pairs are “locked” in the terminology of Büttiker et al. LABEL:Buttiker. Hence, most of the system behaves as an elastic solid, while all of the shear flux is concentrated in a small population of isolated inter-rotor gaps (the slip planes), where the relative angle is “running” in the terminology of Büttiker et al. Buttiker (). Crucially, the position of each slip plane remains fixed for a long duration, far exceeding its internal rotation time , so that the pair of rotors straddling it execute many full turns of relative motion before the slip plane stochastically wanders to a new site. Indeed in the majority of cases studied, after initial transients, slip planes remained fixed for the entire simulation, whilst traversing their potential barrier many thousands of times.

The stability of each solid region or “plateau” within the steady-state requires the overall time-averaged torque within the system (applied by the boundary condition) to be smaller than the steepest gradient in the interaction potential (unity), which is the threshold torque at which an elastic solid region would yield. (This follows from time-averaging Eq. 2b given that, within this region, remains bounded, by definition of “solid”, so that its averaged time-derivative vanishes.) Although the time-averaged torque is uniform at steady state (from time-averaging Eq. 2a), and smaller than the yield threshold, it is nonetheless sufficient to maintain flow within the slip planes, where the potential barrier is repeatedly overcome by virtue of the rotors’ momenta. The traversing of the periodic potential creates a periodic variation in torque at the slip plane, which propagates outwards as damped torsion waves in the elastic solid regions.

This motion can be seen in on-line animations animations () and is also visible in Fig. 4, which is a space-time diagram of the system’s motion. A grey scale represents the value of the inter-rotor potential at each of the 512 sites along the chain (set out horizontally in the figure). Time increases vertically in the figure.

In Fig. 4, all of the shear is concentrated in a single slip plane between two rotors, at a location (selected by spontaneous symmetry breaking in the system) near the right-hand end of the chain. The slip plane is clearly distinguished in the space-time diagram by the repeated set of black marks where the maximum of the potential is repeatedly traversed. High-frequency waves (visible as short slanting lines) propagate outwards only a short distance to either side of the slip plane, before decaying away. Meanwhile very long wavelength torsional waves of lower amplitude travel across the entire system length, visible as long, faint slanting marks in Fig. 4. The gradient of the marks indicates a wave speed of approximately in this case.

### iii.2 Solid-fluid coexistence

Whereas the space-time diagram for slip planes (Fig. 4) shows linear propagation of small-amplitude waves within the solid regions, large-amplitude non-linear waves of yielding events propagate through the fluid region of a solid-fluid coexistence. This can be seen in Fig. 5, the space-time diagram for a steady state of solid-fluid coexistence, where the black marks indicate events where the maximum of the potential is overcome.

It is clear, in Fig. 5, that each yielding event locally displaces stress sufficient to trigger another yielding event nearby, leading to a propagating sequence of rotor-pairs crossing their potential barrier within a finite region of the system. The stress excitations from these events do not propagate across the whole system. The (uniform) time-averaged torque remains just below the yield-point value of unity (the maximum of in Eq. 2b), so a portion of the system is able to exist in a solid state. In fact, due to the finite temperature, some relatively rare yield events are visible in the “solid” region also, hence, it is not strictly solid, but a high-viscosity fluid, so that this coexistence may not be qualitatively distinct from a shear-banding state.

It is surprising that this one-dimensional system has sufficient complexity to self-organise into a stable coexisting state. Some light is shed on the mechanisms responsible for this phenomenon by the mean-field theory in section IV.3.1.

### iii.3 Shear banding and critical states

In the shear-banding state, the X-Y model self-organises into macroscopic regions with different effective viscosities. As in solid-fluid coexistence (which we have seen is a special case of shear banding with a large ratio of effective viscosities), the stress in the shear-banding state is only a little below its maximum (yield-point) value of unity so that a small stress perturbation can trigger a large-amplitude rearrangement.

It appears that a continuous transition separates the shear-banded and uniform states. However, we have not established whether this is an isolated critical point, a line of continuous transitions, or a critical point that terminates a line of first-order phase transitions. The space-time diagram in Fig. 6 is for a system close to criticality, at the transition between uniform and shear-banded states.

A number of critical phenomena are apparent here. One characteristic feature of a nearby critical point is a large correlation length, indicating highly collective/cooperative motion. Although diverging correlation lengths are not possible in one-dimensional, locally-interacting systems at equilibrium, our non-equilibrium steady state exhibits this feature. It can only be recognized by observing the space-time domain, as in Fig. 6, where regions of correlated configurations (in which ) extend across large distances, but not at constant time.

Secondly, the very slow decay of initial transients is typical of critical slowing-down. A third critical phenomenon is the divergence of response functions (discussed more in section III.4 and Fig. 11b), due to the fact that the system appears to be poised in a highly susceptible state. Fig. 7 shows the value of the interaction potential for each rotor pair. The yield point on this potential, at which a stable elastically-bound interaction fails because , is at , where , i.e. halfway up the potential barrier. Note that the system forms a distribution that is clustered around this yield point. By comparison, the inset of Fig. 7 shows the values of in a system at equilibrium (), which are Boltzmann-distributed close to the potential minimum at .

### iii.4 Phase diagrams and measurements of stress and energy

Representative results are shown on some mutually orthogonal slices through the parameter space in Fig. 8. In most cases, the qualitative type of flow behaviour was identified, from time-averaged velocity profiles such as those in Fig. 2. The identification was unambiguous in most cases, but the distinction between shear-banding states and uniform or solid-liquid-coexisting states was sometimes unclear, particularly near criticality.

It is informative to plot data from all of our steady-state simulations on axes versus , allowing us to present all of our simulated phase data in one figure, while Fig. 8 contains only a subset. This is shown in Fig. 9. The simulations were performed at values of ranging from to and from to . Clearly, there is not a full data collapse, but some structure is visible in Fig. 9 — in particular, the uniform region for values greater than on either axis, and general trends along the major and minor diagonals. Further understanding, and some justification for a partial data collapse on axes of versus are provided by the upper bounds in section IV.2 and the approximate scaling analysis in appendix .1.

Empirically, we find a marginally better collapse of the data into fairly clearly-defined phase regions on axes of and , as shown in Fig. 10, although we have no explanation for this.

Measurements of the mean torque in the system, , as a function of shear rate are shown in Fig. 11. These are reminiscent of stress-versus-strain-rate measurements of non-Newtonian fluids measured in rheometric studies of their constitutive relations (e.g. Ref. Schmitt94 ()). A simple Newtonian fluid, on the other hand, would produce a straight line through the origin with gradient equal to its viscosity.

At and , the systems crosses a phase transition directly from the uniform state at high shear rate to the slip-plane state at low shear rate. Notice that, in the slip-plane state, the mean torque remains very constant, independent of . It achieves this by varying the number of slip planes (running pairs) as shear rate varies, as shown in the measurements in Fig. 12, and also in the velocity profiles of Fig. 3. At and , the data in Fig. 11a cross a phase transition from a uniform phase at high shear rate to a solid-fluid coexistence at low shear rate, while, at (and possibly also at ), the system remains uniform.

For the higher value of the friction coefficient shown in Fig. 11b, the system remains in the uniform phase for the temperatures and shear rates shown. Note that the axes are both linear, so very low shear rates are not interrogated by this graph. Despite remaining statistically uniform, the system exhibits large-scale fluctuations, particularly for the data at , as it is close to a continuous phase transition, which we believe to be the same critical point that influenced the data in Fig. 6, discussed in section III.3. The precise location of this critical point (or line) is difficult to establish, as its critical phenomena are observed across a large region of the parameter space.

Close to criticality, response functions diverge. In this case, we can identify the response of the strain rate to changes in stress, as the reciprocal of the gradient in Fig. 11b, which diverges as a continuous transition (possibly to a shear-banded state) is approached from above.

Figure 13 shows the mean potential energy density as a function of shear rate. We define as the value of averaged with respect to position and time. Potential energy density is a more suitable thermodynamic quantity to study than the total internal energy density, because the kinetic energy of a continuously sheared system is non-extensive. Features, due to phase transitions being crossed or narrowly avoided, are again visible in these data.

## Iv Theoretical analysis

Equation. (2a) yields an equation of motion for the relative angle between neighbouring rotors,

 d2Δθjdt2=τj+1+τj−1−2τj (8)

which, together with Eq. (2b), forms a closed set of equations in the relative angles only, independent of any absolute angular positions , as would be hoped for this rotationally symmetric system.

Time-averaging both sides of Eq. (8), and imposing the steady-state condition that, while rotors may have a finite time-averaged angular velocity, their time-averaged angular acceleration must vanish, yields a balance in mean torque-differences . Applying the periodic boundary condition, that the torques are equal at opposite ends of the chain, sets the constant to zero, so that the time-averaged torque is uniform (independent of index ) throughout the steady-state system,

 ⟨τj⟩=¯¯¯τ∀j. (9)

This resembles the uniform-stress condition for steady shear flow of a fluid.

Let us define the local shear rate to be that part of the neighbours’ relative velocity that is persistent rather than fluctuating:

 sj≡⟨Δ˙θj⟩ (10)

which is time-independent, due to the temporal averaging . Note that the local shear rate should not be confused with the global shear rate . The two are equal in the uniform phase. In other phases, is equal to the positional average (with respect to ) of .

Now time-averaging the torque Eq. (2b) and applying Eq. (9) and the fact that the noise has zero mean yields a relationship between the local shear rate and the overall mean torque,

 ¯¯¯τ=μsj+⟨sinΔθj⟩. (11)

Equation (11) can be regarded as a rheological stress-strain-rate relation. The final term in Eq. (11) would vanish for constant relative rotational motion, yielding a Newtonian fluid of viscosity , the constant of proportionality between “stress” (i.e. torque in this rotational case) and strain rate. However, quasi-periodic velocity fluctuations synchronised with (i.e. having a frequency component equal to) the local shear rate lead to a non-zero average of due to a non-uniform angular distribution of , and hence to non-Newtonian rheology.

### iv.2 Upper bounds

Clearly the magnitude of the final term in Eq. (11) cannot exceed unity. Furthermore, it is reasonable to assume that all terms in Eq. (11) are positive for a positive overall shear rate , since the second law of thermodynamics forbids a negative response of torque to imposed shear rate, and local shear rate is expected to share the direction of the global shear rate, as is the conservative part of the torque. Hence .

In any part of the chain that is behaving as a solid, by definition, the local shear rate (the first term on the RHS of Eq. (11)) vanishes, yielding a condition on the overall torque for the existence of any solid regions. Beyond this threshold, no bound structures, even as small as a single pair of rotors, can survive. Hence we anticipate a uniform state when the condition is violated.

The fluid regions (portions of the chain in which the local shear rate is non-zero) also put conditions on . Since the positionally-averaged local shear rate is constrained to equal , the maximum (with respect to ) local shear rate cannot be less than . Hence, at that fastest-shearing location, Eq. (11) gives . Combining this with the other inequality on , and recalling that is independent of position , implies that the steady state is uniform if .

This upper bound on non-uniform states is consistent with the data in Fig. 9. The highest value of for which non-uniform behaviour was identified in any of our steady-state simulations was for a state of closely spaced slip planes at . However, our numerical study was not exhaustive.

A more stringent condition on the parameters applies in the slip-plane regime. In this regime, by definition, no macroscopic fluid regions exist. If is the average length of a solid region within the slip-plane steady-state (i.e.  is the ratio of the system size to the number of slip planes) then the highest local shear rate at any of the slip planes is so that, from Eq. (11), . Hence, the condition for non-yielding of the solid regions gives , which is a more stringent condition than that for survival of any micro-solid structures. The definition of the slip-plane regime requires , the equality occurring in the limit where half of the inter-rotor spaces are slip planes, so that each rotor is solidly bonded to only one neighbour. Hence the condition for the slip-plane phase is . As noted above, this upper bound was realised in one of our simulations (see Fig. 9), but never exceeded.

### iv.3 Linearization

The full equations of motion (8 and 2b) are non-linear in the unknown functions . To find solutions, we shall linearize them. This is not simply a case of replacing the sinusoid in Eq. 2b by its small-angle limit, as one might at equilibrium, since, in the boundary-driven system, (some of) the relative angles are ever-increasing functions of time.

At steady state, can be written in terms of a part that grows linearly at rate (the local shear rate), a constant offset , and a bounded fluctuating part with zero mean, thus:

 Δθj(t)=sjt+cj+εj(t). (12)

Then the (still exact) equations of motion become

 ¨εj(t)=τj+1+τj−1−2τj\vspace−4mm (13a) τj=sin(sjt+cj+εj(t))+μsj+μ˙εj(t)+ηj(t). (13b)

Time averaging Eq. 13b and using the boundedness of allows Eq. 11 to be recast as

 ¯¯¯τ=μs+⟨sin(sjt+cj+εj(t))⟩ (14)

Next, we expand in terms of Fourier modes with coefficients thus:

 εj(t)=∑ω^εj(ω)eiωt (15)

where the summation is over a discrete but infinite set of frequencies that are not necessarily equally (or even finitely) spaced, since the motion may be aperiodic. A continuous Fourier transform was not used in Eq. 15 because has a discrete spectrum of delta functions, excited by the continual rotation of the interacting pairs of rotors. Similarly, the torque is expanded as

 τj(t)=∑ω^τj(ω)eiωt. (16)

Real-valued and require the positive and negative-frequency spectral components to respect the symmetry and where an asterisk () indicates the complex-conjugate.

The inverse transformations (Fourier analysis) can be performed by applying temporal averaging on the infinite time domain

 ⟨e−iωtεj(t)⟩=∑ω′^εj(ω′)⟨ei(ω′−ω)t⟩=^εj(ω) (17)

since where the Krönecker delta has its usual definition,

 δ(ω′,ω)={1 if ω′=ω0 otherwise,

despite the fact that it arguments are non-integer. Note that and and Eq. 13a becomes

 −ω2^εj(ω)=^τj+1(ω)+^τj−1(ω)−2^τj(ω). (18a) Analysis in this section is thus far exact. Henceforth, we shall assume, as observed in many of the steady states simulated, that fluctuations about constant motion (with respect to time) are small, and so linearize Eq. 13b in εj(t), yielding the torque-spectrum ^τj(ω)=μsjδ(ω,0)+12ib∗jδ(ω,−sj)−12ibjδ(ω,sj)+^ηj(ω)+iωμ^εj(ω)+12bj^εj(ω−sj)+12b∗j^εj(ω+sj) (18b)

where and . The zero-frequency component of the torque-spectrum is the mean torque

 ¯¯¯τ=12i(b∗j−bj)δ(sj,0)+μsj+12bj^ε(−sj)+12b∗j^ε(sj). (19)

The first term on the right-hand side of Eq. 19 represents the elastic torque in a solid region of the chain (where the local shear rate vanishes), which depends (via ) only on the mean angle by which neighbours are twisted. This is the only non-zero term for , in which case Eq. 19 evaluates to . The second term represents the Newtonian contribution to the viscosity. The non-Newtonian contribution of the final two terms arises only (at linear order) from fluctuations at the same frequency as the local shear rate, which synchronise with the periodic traversing of the potential, thus causing non-uniform occupancy of the potential, due to the rotor lingering for longer at some relative angles than others.

In principal, one can solve Eqs 18a and 18b for , then substitute and its complex conjugate into Eq. 19 to find the mean torque for any given configuration of local shear rates and offsets , yielding a spatial constitutive relation for the XY model. The uniformity (independence of ) of the torque then restricts the allowed configurations of and to those exhibited by the steady-state phases.

In practice, solution of Eqs 18a and 18b remains difficult in general, and requires some further approximation. We present two alternative methods in sections IV.3.1 and IV.3.2 below. However, for the particular case of a solid region of the chain (one for which ), the Fourier modes become decoupled, so Eq. 18b simplifies to

 ^τj(ω)=¯¯¯τ+^ηj(ω)+(iωμ+√1−¯¯¯τ2)^εj(ω) for sj=0. (20)

In this case, Eq. 18a becomes a discretized noisy wave equation. At zero temperature, is has a simple travelling-wave solution of the form

 ^εj+1=α(ω)e±iω/v(ω)^εj (21)

with a frequence-dependent (i.e. dispersive) wave speed due to dissipation and discretization. In the low-frequency, low-dissipation limit, the wave speed varies with overall torque as , vanishing in the limit of stability . Closed-form expressions for the real values and are straightforwardly obtained by substitution of Eqs. 20 and 21 into 18a, but are not reproduced here, as they are large. We shall use this travelling-wave solution for the elastic solid in section IV.3.2 to find an approximate solution for slip planes. Firstly, we shall develop an approximation for a general steady state.

#### iv.3.1 Effective Medium Theory

Here, we shall develop a mean-field theory, in which the fluctuations of neighbours and of a given site are uncorrelated with it. In Eq. 18a we average (with respect to complex phase) over the fluctuations of the neighbours, so that only their zero-frequency contributions remain, yielding

 ω2^εj(ω)=2^τj(ω)−2¯¯¯τδ(ω,0). (22)

This is equivalent to replacing by in Eq. 8, by placing a single rotor-pair in a medium of constant torque. In a fluid region, where , the one remaining phase constant becomes arbitrary, so we may choose implying . On the other hand, if , we must retain as a parameter.

Now that we are considering only a single shear rate , Eq. 18b and 22 will generate harmonic spectra, i.e. only at integer multiples of . So, to simplify the notation, for integer , let us define new coefficients and . Hence and .

So Eq. 19 yields

 ¯¯¯τ=sin(cj)δ(s,0)+μs+R(a1) (23)

where if . For non-zero , Eqs. 22 and 18b yield

 ns(ns−2iμ)an=2(μs−¯¯¯τ)δ(n,0)+iδ(n,−1)−iδ(n,1)+2ξn+an−1+an+1. (24)

To make further progress, in the mean-field spirit, we drop (average over) the noise term , equivalent to setting zero temperature. For , Eq. 24 yields

 a1=ia2a1−s(s−2iμ). (25)

For , Eq. 24 yields a simple recurrence relation for the ratio of successive harmonic amplitudes for ,

 1Rn=ns(ns−2iμ)−Rn+1. (26)

Using Eq. 26 to define the as-yet undefined variable allows us to combine Eqs. 23, 25 and 26 into a neat expression of the fluid’s consititutive relation, in terms of the imaginary part of a continued fraction,

 ¯¯¯τ = sin(cj)δ(s,0)+μs+I(R1) R1 = 1s(s−2iμ)−12s(2s−2iμ)−13s(3s−2iμ)−1…. (27)

Equation 27 is plotted in Fig. 14 for various values of the friction coefficient .

Note that we find values of mean torque that depend on , not only on the combination and as approximated in appendix .1. This is not unexpected, since the scaling approximation in the appendix requires the correlated motion of several nearby rotors to behave as a continuum, whereas the mean-field theory yielding Eq. 27 analyses only a single rotor-pair. It is worth remarking that the scaling collapse of Fig. 9 is only approximate, so that the qualitative effect of varying in Fig. 14 is not inconsistent.

A local constitutive relation, or “flow curve” (a fluid’s stress as a function of local shear rate) of the kind plotted in Fig. 27 has a standard interpretation Olmsted (). A negative gradient of the curve indicates an unstable regime. No fluid can remain uniform at a shear rate where the gradient is negative, as the unbalanced stresses on a region of the fluid that is infinitesimally perturbed from uniformity will act to enhance the perturbation. Furthermore, because the time-averaged torque must be uniform at steady state, any inhomogeneous fluid must exhibit a coexistence of local shear rates at equal values of . Hence, a shear-banding state, or a solid-fluid coexistence, can only exist at values of the parameters and for which a horizontal line has multiple intersections with the local flow curve.

The above criteria are not sufficient to uniquely define the extent of the coexistence behaviour from the curves in Fig. 14, but do put bounds on it. Specifically, since the overall shear rate is a volume-weighted average of the coexisting local shear rates, a uniform steady state cannot exist at an overall shear rate where the flow curve has negative slope. And, inverting the flow curve, a uniform steady state must exist wherever the local shear rate is single-valued as a function of .

For , Eq. 27 describes a solid region capable of supporting any torque in the range (depending on the value of ), corresponding to the vertical line segment from to in Fig. 14. The horizontal dotted line is the upper limit on the solid regions’ torque. So curves that cross the dotted line represent possible solid-fluid coexistences.

Hence, Fig. 14 predicts (for ) only uniform states for . On decreasing , a critical point is predicted just below and , where a negative gradient first appears. Below this critical value of , shear-banding takes place for a finite range of shear rates. Only at lower values (below ) is solid-fluid coexistence possible, according to the approximate solution.

Comparing with the simulation results, and extrapolating the simulated trends to zero temperature where simulation results are difficult to obtain due to long-lived transients, we see that the approximate theory, captures some qualitative features. In particular, in simulations, a uniform state exists at high and, on decreasing the shear-banded state is entered via a continuous transition, whereas the “freezing” transitions (to solid-fluid coexistence or slip planes) appear to be first-order, and occur at lower .

Furthermore, quantitative comparison can be made with values of the mean torque measured in simulation. For the torques shown in Fig. 11a, from simulations at , temperature has a significant effect. Only the lowest-temperature data (shown by squares) should be compared with the mean-field theory, which is effectively a zero-temperature theory. In agreement with the theoretical curve, the data (Fig. 15a) are almost indistinguishable from a straight line (a Newtonian fluid) for all high shear rates.

A phase transition to a coexisting state is characterized, in the simulation data, by a transition to a region of vanishing gradient. Note that the data are plotted against the global shear rate, averaged over the coexisting regions, whereas the mean-field curve is a function of local shear rate, and so can only put bounds on the volume-weighted mean of coexisting phases. (Recall that coexistence is necessary wherever the local constitutive relation has a negative gradient, and forbidden wherever local shear rate is single-valued as a function of .) Hence, Fig. 11a demonstrates that the mean-field theory is completely consistent with the measured phase transition(s) and quantitatively fairly accurate at .

For any single-phase regime, the theory makes a unique prediction. A single phase exists for the data in Fig. 11b at , for which temperature has little effect. They are compared with the mean-field curve in Fig. 15b. Although imperfect, the mean-field theory is not far from the data.

#### iv.3.2 An isolated slip plane

We can solve for the motion of the fully interacting linearized system (i.e. no longer neglecting the correlated fluctuations and offsets of neighbouring rotor-pairs) for the specific case of the slip-plane phase at , subject to the approximation, to be introduced below, that the motion is dominated by a single frequency (the rotation rate of the local slip plane). We consider an isolated running pair (where increases with time as Eq. 12 with finite ) in a chain of otherwise locked () rotors. Let us assume that the slip planes are sufficiently far apart that we may consider only one isolated slip plane.

The running pair (slip plane), which we locate at without loss of generality, sends out small-amplitude waves into the surrounding solid chain of rotors. The motion of the locked rotor pairs, for all is given in terms of by Eq. 21, and similarly the variables at negative are eliminated in favour of . Finally, are eliminated in favour of by applying Eq. 18a for with the torque given by Eq. 20 for and Eq. 18b for .

Given that the pair at is the only running pair, the absolute phase of its motion is irrelevant, so we may set hence , so that Eq. 18b for the torque at simplifies a little. Nevertheless it couples Fourier coefficients at all harmonic frequencies, thus eventually yielding a recurrence relation for these harmonics, analogous to the one found for effective-medium theory (Eq. 24), but much less neat (hence we do not show it here). Unfortunately, we have not found a full solution to this recurrence relation, as we did in section IV.3.1. So we further approximate that the motion is dominated by the frequency , and drop all higher harmonics (setting for ).

The resulting relationship between mean torque and rotation rate of the running pair is plotted in Fig. 16 for various .

Recall (as in section IV.2) that, if the (effectively isolated) slip planes are separated by an average distance of rotors, then the overall shear rate is related to the local shear rate at the slip plane by . Since is confined to a finite (and, in some cases small) range of values (see Fig 16) — indeed only to the subset of stable values for which the gradient of is positive — while is fixed by the boundary conditions, it follows that is confined to vary approximately with . Hence the number of slip planes varies approximately linearly with , as observed in the simulation results of Fig 12. The gap (in values), between the stable regimes at zero shear rate and at finite , is responsible for the existence of the slip-plane phase, which is unable to re-distribute its localized points of high shear rate into more diffuse regions of lower shear rate.

While different in detail from the curves in Fig. 14, we note that the solutions in Fig. 16 are close in value to those in Fig. 14 which were derived (at mean-field level) for local shear rate in an arbitrary phase. This agreement therefore lends further support to the validity of the earlier mean-field result.

Example trajectories for the running pair of rotors and some of its neighbours in the solid region are shown in Fig. 17. The further into the solid region, the smaller the amplitude of the fluctuations about the average value. The fluctuations are successively out of phase by the same amount, indicating a travelling torsional wave, like those visible in the simulation data of Fig. 4 as diagonal lines emanating a short distance from the slip plane. Hence, we see that each slip plane disturbs its environment via small-amplitude waves that decay as they propagate. Two nearby slip planes will interact appreciably only if their separation is small compared with the decay length of the torsional waves.

## V Conclusions and prospects for further work

It is remarkable how closely the shear flow of the angular-momentum-conserving 1D XY model mimics the rich behaviour of non-Newtonian complex fluids. We have observed that, independent of its initial conditions, the three-parameter model exhibits reproducible steady-states with four main types of velocity profile, as shown in Fig. 2: uniform flow, shear-banding, solid-fluid coexistence and slip-plane states. There may be yet other, more exotic states in small regions of parameter-space.

All four of the model’s main flow regimes are routinely seen in diverse types of soft matter, which have hitherto each been analysed using different ad hoc phenomenological models that are often quite complex.

The macroscopic phenomenology of non-Newtonian fluids is typically captured theoretically using constitutive models, normally tensorial models that define the stress tensor in terms of the mesoscopic structure, for example the (diffusive) Johnson-Segalman model Johnson77 (); Olmsted00 (). The dynamics of the structure is specified and then the model can be solved for the flow properties either in or out of the steady state, usually numerically. While tensorial models are needed to describe real fluids, simpler non-tensorial models have previously been defined phenomenologically, without microscopic degrees of freedom (e.g. Fielding04 (); Dubbeldam09 ()), and can reproduce phenomena such as the instability in the constitutive curve described above.

The reproduction of so many non-Newtonian macroscopic phenomena in a single, simple, microscopic model was quite unexpected, and without precedent. It offers further insights into the microscopic and macroscopic features of flow-induced phase transitions.

By numerically charting some of the phase diagram and explaining some of the behaviour by linearization and mean-field analysis, we have demonstrated the XY model’s great richness and its usefulness as an idealized model of non-equilibrium fluids. Furthermore, we have shown that there is much scope for further research on this system, which we would encourage others to investigate.

In particular, more numerical work is required, to explore the phase diagram more comprehensively and in greater detail, and also to fully characterize the nature of the various phase transitions, which will contribute to a wider understanding of non-equilibrium condensed matter in general. Once those phase transitions have been characterized, it will be interesting to make comparative experimental studies of complex fluids that exhibit equivalent behaviour.

Theoretical progress is possible in a number of directions. It may be feasible to solve the linearized version of the model (i.e. to first order in deviations from constant but non-uniform shear rate ) without applying the effective medium approximation that removes spatial features. Alternatively, a perturbative expansion in , perhaps to second or third order, would reveal localized interactions that govern the curvature of velocity profiles such as Fig. 2b and create an effective interfacial energy between coexisting phases, leading to their macroscopic separation.

It would be appealing to find an approximate continuum version of the model, expressed in terms of differentiable fields that are functions of continuous position and time only. The difficulty is that, although the local torque can, in some regimes, be assumed smoothly varying (as in Appendix .1), the relative angle cannot, and must remain indexed by a discrete index (although its time-derivative may be made continuous). The reason is that is periodic modulo , so that the torque (which depends on ) is invariant under a full turn of one rotor relative to its neighbour. A continuous field would acquire a topological defect with increasing winding number each time a full turn was introduced, which would be appropriate for modelling an elastic solid with memory of total acquired strain (i.e. a twisting rubber belt), but not a fluid under continuous shear flow. Overcoming this difficulty, to derive a Ginsburg-Landau-like model for the steady-state profile of local shear rate would require approximately “integrating out” (i.e. solving for) the fluctuations .

Some other theoretical approaches with potential for further progress, including some ideas for renormalization group analysis, can be found in Ref. WelshThesis (), along with some further numerical data on this elegant model.

## Vi Acknowledgements

We are grateful for the advice of Adrian Baule. Parts of this work were funded by a Royal Society University Research Fellowship, an EPSRC grant GR/T24593/01 and an EPSRC doctoral training grant.

APPENDIX

### .1 Approximate scaling

To investigate effects of re-scaling the model’s parameters, let us non-dimensionalize time by measuring it in units of reciprocal shear rate. We define the reduced-time , and note that the noise function scales non-trivially with time and with friction coefficient and temperature. Writing it explicitly as a function, , we note, from Eq. 3 the following invariance with respect to re-scaling of its arguments:

 ηj(at,bμ,abT)↔ηj(t,μ,T) (28)

where the equivalence means that statistical properties (specifically, all moments) of the functions are equal. Hence we may write

 ηj(t,μ,T)=ηj(^t˙γ,μ,T)↔ηj(^t,μ˙γ,T). (29)

So Eq. 2b becomes

 τj=sinΔθj+(μ˙γ)dΔθjd^t+ηj(^t,μ˙γ,T). (30)

So we see that the formula for the torque retains the form of Eq. 2b under the temporal re-scaling if the system parameters are re-scaled according to

 (T,μ,˙γ)→(^T,^μ,^˙γ)=(T,μ˙γ,1) (31)

where the transformation of follows from the fact that it is the positional average of .

Although Eq. 2b is invariant under the above transformation, Eq. 8 is not. To find a symmetry of this pair of equations, we must first approximate Eq. 8, by appealing to the fact that its right-hand side is a discrete version of a second derivative with respect to position. If we assume that torque varies smoothly with position, then we may replace the discrete index by a continuous position and approximate the right-hand side of Eq. 8 by a second derivative. Under that assumption, if we define our unit of length to be the distance between rotors, then Eq. 8 becomes

 d2Δθjdt2≈∂2τ(x,t)∂x2. (32)

Our observations of the steady-state dynamics lead us to believe that the assumption of smoothly-varying torque is a reasonably realistic approximation in some but not all regimes.

As explained in section V, although can (in some regimes) be assumed smoothly varying, the relative angle cannot. Hence the left-hand side of Eq. 32 remains a function of , indexed by , and the right-hand side is evaluated at .

Under the temporal re-scaling , Eq. 32 transforms to

 d2Δθjd^t2=1˙γ2d2Δθjd t2≈1˙γ2∂2τ(x,t)∂x2 (33)

which can be written as

 d2Δθjd^t2≈∂2τ(x,t)∂^x2 (34)

if we define a re-scaled position .

Finally, then comparing Eqs. 34 and 30 with 32 and 2b, we observe that, for states with smoothly varying torque, a system with parameters respects approximately the same equations of motion as a system with parameters and re-scaled length and time. Hence, these systems will exhibit the same steady-state phase behaviour. Therefore, we expect some partial data-collapse when plotting the phase diagram on axes of and only. As seen in Fig. 9, the resulting data-collapse is not complete, indicating that local details of the torque variation control a significant amount of the macroscopic behaviour. Nevertheless, the scaling analysis in this appendix — particularly the exact scaling in Eq. 30 — may prove useful for future studies.

## References

• (1) Kosterlitz and Thouless, Ordering, stability and phase transitions in two-dimensional systems. J Phys C 6, 1181 (1973).
• (2) K. Krishan and M. Dennin, Viscous shear banding in foam Phys. Rev. E 78, 051504 (2008).
• (3) G. Debrégeas, H. Tabuteau, and J.-M. di Meglio, Deformation and flow of a two-dimensional foam under continuous shear, Phys. Rev. Lett. 87 (2001), 178305.
• (4) R. Besseling, L. Isa, P. Ballesta, G. Petekidis, M. E. Cates, and W. C. K. Poon, Shear Banding and Flow-Concentration Coupling in Colloidal Glasses. Phys. Rev. Lett. 105, 268301 (2010).
• (5) F. Ianni, R. Di Leonardo, S. Gentilini, and G. Ruocco, Shear-banding phenomena and dynamical behavior in a Laponite suspension. Phys. Rev. E 77, 031406 (2008).
• (6) M.D. Haw, Jamming, two-fluid behaviour and âselffiltrationâ in concentrated colloidal suspensions, Phys. Rev. Lett. 92 (2004), 185506.
• (7) V. Schmitt, F. Lequeux, A. Pousse, and D. Roux. Flow behaviour and shear induced transition near an isotropic/nematic transition in equilibrium polymers. Langmuir 10, 955 (1994).
• (8) M. A. Fardin, D. Lopez, J. Croso, G. Grégoire, O. Cardoso, G. H. McKinley, and S. Lerouge, Elastic Turbulence in Shear Banding Wormlike Micelles. Phys. Rev. Lett. 104, 178303 (2010).
• (9) J.-B. Salmon, S. Manneville, and A. Colin, Shear banding in a lyotropic lamellar phase. I. Time-averaged velocity profiles. Phys. Rev. E 68, 051503 (2003).
• (10) J. L. Jones and T. C. B. McLeish, Rheological Response of Surfactant Cubic Phases. Langmuir 11, 785 (1995).
• (11) J. Cao and A. E. Likhtman, Shear Banding in Molecular Dynamics of Polymer Melts. Phys. Rev. Lett. 108, 028302 (2012).
• (12) I. Kunita, K. Sato, Y. Tanaka, Y. Takikawa, H. Orihara, and T. Nakagaki, Shear Banding in an F-Actin Solution. Phys. Rev. Lett. 109, 248303 (2012).
• (13) P. Tapadia and S.-Q. Wang, Yieldlike constitutive transition in shear flow of entangled polymeric fluids. Phys. Rev. Lett. 91, 198301 (2003).
• (14) E. Pouyan Boukany, S.-Q. Wang, and X. Wang, Step shear of entangled linear polymer melts: new experimental evidence for elastic yielding, Macromolecules 42, 6261 (2009).
• (15) R. M. L. Evans, Craig A. Hall, R. Aditi Simha, and Tom S. Welsh, Classical XY Model with Conserved Angular Momentum is an Archetypal Non-Newtonian Fluid. PPhys. Rev. Lett. 114, 138301 (2015).
• (16) D. C. Mattis, Transfer matrix in plane-rotator model. Phys. Lett. A 104, 357 (1984).
• (17) M. J. Unterberger, H. Weisbecker, G. A. Holzapfel, Mechanical modeling of rheometer experiments: Applications to rubber and actin networks. Int. J. Non-Linear Mech. 67, 300 (2014).
• (18) M. Harvey and T. A. Waigh, Optical coherence tomography velocimetry in controlled shear flow. Phys. Rev. E 83, 031502 (2011).
• (19) R. D. Groot and P. B. Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys. 107:4423 (1997).
• (20) A. W. Lees and S. F. Edwards, J. Phys. C: Solid State Phys. 5, 1921 (1972).
• (21) Animations of rotor dynamics in the uniform and slip-plane states are available in the supplemental material to Ref. PRL () at http://link.aps.org/ supplemental/10.1103/PhysRevLett.114.138301 and also at https://sites.google.com/site/rmlevans/xy- model-animations.
• (22) M. Büttiker, E. P. Harris, and R. Landauer, Thermal activation in extremely underdamped josephson-junction circuits. Phys. Rev. B 28, 1268 (1983).
• (23) P.D. Olmsted, Two-state shear diagrams for complex fluids in shear flow, Europhys. Lett. 48, 339 (1999).
• (24) M. Johnson and D. Segalman. A model for viscoelastic behaviour which allows non-affine deformation. Journal of Non-Newtonian Fluid Mechanics, 2, 255 (1977).
• (25) P. D. Olmsted, O. Radulescu, and C.-Y. D. Lu. Johnson Segalman model with a diffusion term in cylindrical Couette flow. Journal of Rheology, 44, 257 (2000).
• (26) S. M. Fielding and P. D. Olmsted. Spatiotemporal oscillations and rheochaos in a simple model of shear banding. Phys. Rev. Lett. 92, 084502, (2004).
• (27) J. L. A. Dubbeldam and P. D. Olmsted. Two-dimensional perturbations in a scalar model for shear banding. Euro. Phys. J. E, 29, 363 (2009).
• (28) T. S. Welsh, Statistical mechanics of boundary driven systems, University of Leeds thesis (2012).

## References

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