Shear-density coupling for a compressible single-component yield-stress fluid

Flow behavior of a single-component yield stress fluid is addressed on the hydrodynamic level. A basic ingredient of the model is a coupling between fluctuations of density and velocity gradient via a Herschel-Bulkley-type constitutive model. Focusing on the limit of low shear rates and high densities, the model approximates well---but is not limited to---gently sheared hard sphere colloidal glasses, where solvent effects are negligible. A detailed analysis of the linearized hydrodynamic equations for fluctuations and the resulting cubic dispersion relation reveals the existence of a range of densities and shear rates with growing flow heterogeneity. In this regime, after an initial transient, the velocity and density fields monotonically reach a spatially inhomogeneous stationary profile, where regions of high shear rate and low density coexist with regions of low shear rate and high density. The steady state is thus maintained by a competition between shear-induced enhancement of density inhomogeneities and relaxation via overdamped sound waves. An analysis of the mechanical equilibrium condition provides a criterion for the existence of steady state solutions. The dynamical evolution of the system is discussed in detail for various boundary conditions, imposing either a constant velocity, shear rate, or stress at the walls.


I. INTRODUCTION
Heterogeneous flow and shear banding are ubiquitous phenomena, commonly occurring in a variety of complex fluids such as polymer solutions and worm-like micelles [1][2][3], colloidal gels [4], hard sphere colloidal glasses [5,6] and granular media [7].In line with this diversity of the physical systems, one encounters different underlying mechanisms as being responsible for localized flow.Classically, shear banding occurs in systems with a strongly shear-thinning flow curve (stress versus imposed shear rate) [8].Alternatively, banding can result for a non-monotonic flow curve stemming from a shear-induced phase transformation.In this case, an instability occurs if the globally imposed shear rate lies between the two solutions corresponding to homogeneous steady flow.The system divides into two regions, each flowing with one of the stable shear rates [1].For colloidal gels, on the other hand, the mechanism of shear localization is attributed to a competition between formation and growth of fractal-like clusters and its shear-induced breakage [4].
An interesting case occurs in dense suspensions of hard sphere colloidal particles and granular materials, where the underlying flow curve is monotonic, yet the flow can develop spatio-temporal heterogeneities [5,7].In these "soft glassy materials" [9], shear-induced rejuvenation competes with the sluggish relaxation (aging) kinetics and may lead to a heterogeneous flow in the glassy state [10][11][12].
Flow localization in dense hard-sphere suspensions has been recently rationalized in terms of the so-called shearconcentration coupling (SCC) [5,6], a hydrodynamic model, first proposed in Ref. [13], which couples the local flow to the concentration field.This coupling is encoded in a non-Newtonian stress and in a shear-rate dependent osmotic pressure.
Within SCC, one considers a background fluid which transports-and is influenced by-a concentration field.While this picture emerges naturally in the case of polymer solutions, the role of the background fluid is less obvious in hard sphere colloidal glasses.Indeed, there is a common consensus that the effect of hydrodynamic interactions can be neglected in colloidal hard-sphere systems close to the glass transition in the low shear rate limit, which is of primary interest to the present study [14,15].Accepting this standpoint, it is tempting to fully neglect the background fluid and investigate the issue of flow heterogeneity within hydrodynamic equations of a single-component non-Newtonian fluid.This paper presents such a study.
In Refs.[5,13,16], the instability of a sheared colloidal suspension has been investigated based on an advectiondiffusion equation for the colloid concentration ρ, embedded in a solvent of velocity u, where j denotes the total particle flux, ζ is a friction coefficient, and µ is a (shear-rate dependent) generalized chemical potential [17].Equation (1) asserts that the total flow velocity j/ρ of the colloidal particles consists of an imposed "background" flow u, onto which a contribution −(1/ζ)∇µ due to the diffusive motion of the particles is superimposed.The flow velocity u is assumed to be governed by the Stokes equation, where σ is the viscous stress tensor, which is typically given in terms of an expansion in gradients of u.The Greek symbols stand for spatial directions (α, β ∈ {x, y} in the present 2D study) and Einstein's sum rule over repeated indices is used.In Ref. [6], the possibility of a coupling between shear and concentration has been investigated in a system of hard spheres.A constant kinetic temperature has been imposed by continuously rescaling the particle velocity during the simulations.Notably, there is no background fluid in the system investigated in Ref. [6].Thus, it can be considered as an isothermal compressible single-component fluid, described by a continuity equation for the particle density ρ and a transport equation for the fluid momentum ρu: As in Ref. [6], Π and σ denote the reversible and the irreversible (viscous) stress tensors.In close analogy to shear concentration coupling, one postulates a coupling between fluid density and local shear rate, which we shall call "shear-density coupling" (SDC) in the following.As shown in section II, this coupling is generated by reversible and viscous stresses being functions of the shear rate and density, respectively.In equilibrium, the divergence of the reversible stress tensor can be related to a chemical potential via ∂ β Π αβ = ρ∂ α µ.Beyond equilibrium, this relation serves as a definition of a shear-rate dependent chemical potential.
Before proceeding further with our analysis, a comment on the above equations is at order here.The advectiondiffusion equation ( 1) is central to dynamic density functional theory and widely used for the description of driven colloidal suspensions [18][19][20].In these approaches, u represents the velocity of the background fluid, which consists of an externally imposed component (e.g., shear flow) and a contribution arising from the hydrodynamic inter-particle interactions.Notably, the dynamics of a subset of tagged particles in a single-component fluid flowing with velocity u is formally also described by eqs.( 1) and (2) [21,22].In this case, the chemical potential and the viscosity would react only to the fluctuations of the tagged particles.However, the viscosity and the pressure are actually sensitive to the total density, since this quantity describes the caging and trapping responsible for the dynamic slowing down near the glass transition.
In view of these arguments on the single-component fluid nature of the problem, it appears more appropriate to analyze the hard-sphere system of Ref. [6] in terms of the isothermal compressible fluid equations in eq. ( 3), rather than an advection-diffusion equation.In passing, we remark that the different nature of the two models is also crucial in the case of critical phenomena: here, the advection-diffusion and momentum transport equations define the universality class of "model H", which primarily describes a binary fluid mixture [23].The isothermal single-component fluid, instead, is described in terms of a continuity equation and a momentum equation, giving rise to a critical behavior distinct from model H [24].

II. MODEL
We consider a fluid described by eq. ( 3), bounded by walls at y = 0 and y = L (see fig. 1).The flow is assumed to be homogeneous along the vorticity direction (z) as well as along the flow direction (x), such that generally In this and the following section, we focus on bulk dynamics, such that specification of the boundary conditions at the walls is not necessary.We shall therefore merely assume the presence of a constant steady background shear rate γ0 .(We return to the effect of boundary conditions in section IV, where we numerically solve the Navier-Stokes equations in a finite domain.)The pressure tensor Π is isotropic, Π αβ = Πδ αβ , where Π denotes the scalar pressure.Consequently, eq. ( 3) reduces to Figure 1.Slit geometry considered in the present study.The fluid is subjected to a steady shear flow in lateral direction (x) with a spatially constant (background) shear rate γ0, but can develop arbitrarily large deviations described by a local shear rate γ(y, t).We assume the fluid to be homogeneous both in the lateral and the vorticity direction (z, pointing normal to the figure plane).
Analogously to Refs.[16,25,26], we take the following form for the viscous stress tensor σ: which corresponds to an expansion in gradients of the flow field u respecting certain symmetry properties of the stress [27].Here, η and ζ denote the shear and bulk viscosity, respectively, which are generally functions of the density and the shear rate (see below).The parameter κ denotes the shear-curvature viscosity, and the stress contribution associated with it serves to stabilize the flow field against large gradients.Analogously, the parameter κ controls the corresponding contribution stabilizing the bulk viscous stress.The yield stress σ yield is independent of γ and is nonzero only in the glassy phase (ρ > ρ g ).In contrast to the shear viscosity η (see below), detailed data for the bulk viscosity ζ and the curvature viscosities κ and κ in a hard-sphere fluid near the glass transition are not available.Following Ref. [16], we shall therefore assume these viscosities to have the same functional form as η, i.e., Here, η 0 and κ 0 are the shear (curvature) viscosities in the zero-shear rate (Newtonian) limit [see eq.(13b)], and b is a free dimensionless parameter.Typically, we set b = 1 and κ 0 /η 0 (10 − 100)a 2 , where a is a microscopic length scale, e.g., the average particle diameter in a colloidal glass.This choice gives rise to an effective interface width, ∼ κ 0 /η 0 , of the shear band of a few particle diameters a [16].Using eq. ( 4), the relevant components of the viscous stress tensor follow as with b ≡ b + 4/3 = 7/3.In order to track the influence of the bulk viscosity, we shall carry along the parameter b in our calculations.Summarizing, eq. ( 5) reduces to where we defined which are generally functions of ρ and γ.
It seems reasonable to assume where σ yield is a common yield stress function.In the liquid phase (ρ < ρ g ), the yield stress vanishes and the shear viscosity is well described by a Krieger-Dougherty relationship (cf.Ref. [16]): Here and in the following, Φ ≡ ρ/ρ m , where ρ m = 0.67 (in appropriate units, see below) is the packing fraction corresponding to random close packing of (polydisperse) hard spheres.
In the glassy phase (ρ > ρ g ), instead, MD simulations of a hard-sphere system indicate [6] where the parameters σ 0 0.0119 k B T /a 3 , A = 34.5 (η 0 a 3 /(k B T )) n , p 2.355, n 0.4 result from a fit.k B T denotes the thermal energy and ρ g = 0.585 is the density of the glass transition.The pressure is given, for any ρ, by [6] with The shear-rate dependence of Π is a manifestation of the flow-induced distortion of the pair-correlation function.We remark that the parameters in eqs.( 14) and ( 15) have been obtained in Ref. [6] from a fit to the global flow curves, taking γ ≡ γ0 , but are assumed here to apply also locally in the system.We shall henceforth fix the units of mass, length and time by setting k B T = a = η 0 = 1.With these choices, the fundamental "microscopic" time scale t 0 ≡ η 0 a 3 /k B T = 1.Using the fact that η 0 is the fluid viscosity in the dilute limit [see eq. ( 13b))] and invoking the Stokes-Einstein relation, one obtains t 0 ∼ a 2 /D with the self-diffusion coefficient D. In other words, t 0 is the time needed for a particle to explore, in the dilute limit, a distance comparable to its own size.Noteworthy, this is also a measure of the structural relaxation time.Accordingly, the microscopic time scale t 0 determines, together with thermal energy and particle size, the viscosity and stress scale.In the context of macroscopic fluid dynamics, however, a more natural dimensionless measure of time, which we shall use in the discussion of our results, is instead given by the strain t γ.

A. Linearization of the dynamics
We consider small fluctuations of the density and the shear-rate, i.e., ρ(y, t) = ρ 0 + δρ(y, t), γ(y, t) = γ0 + δ γ(y, t), where ρ 0 and γ0 denote the uniform background values.In linear order in the fluctuations and derivatives, eq. ( 10) becomes where now the coefficients σ γ,ρ , Π γ,ρ , η, and κ are understood to be evaluated for the background values ρ 0 and γ0 .In order to develop a basic understanding of the transport mechanisms in the compressible fluid, note that, inserting eq.(16a) into eq.(16b), the latter becomes a generalized diffusion equation for the shear rate fluctuation δ γ, While the last term on the r.h.s. is typically negligible, the first and the second term induce a smoothing of shear rate inhomogeneities.However, due to the third term, which is not present in a Newtonian fluid, a positive density fluctuation can effectively lower the local shear rate.Such a negative shear rate fluctuation drives, via the first term on the r.h.s. of eq.(16c), a flow which [via eq. ( 16a)] further enhances the density in that region.This gives rise to a feedback mechanism, which is further analyzed in section III B. In passing, we note that eqs.(16a) and (16c) can be combined into a generalized "sound-wave" equation The dynamics induced by the above compressible fluid equations is further discussed and contrasted to a diffusive transport model in appendix A.
In order to investigate the linear stability, we solve eq. ( 16) via the ansatz where ω and k represent the growth rate and wavenumber of a fluctuation, respectively, and the bared quantities denote the fluctuation amplitudes.This ansatz transforms eq. ( 16) into which can be written in matrix form as with the abbreviations and A nontrivial solution of eq. ( 21) requires the coefficient matrix to be singular and, correspondingly, the determinant to vanish: Note that, in order to obtain a purely real solution, the ansatz in eq. ( 19) must be linearly combined with an expression of the same form but where k is replaced by −k.The three roots w 1,2,3 of the cubic equation ( 24) are independent of the sign of ±k.Accordingly, we can write the general solution to the linearized hydrodynamic equations ( 16) as The coefficient vectors A, B, . . .are obtained by inserting each root ω j into eq.( 21) and determining the null-space of the resulting linear mapping.

B. Boundary of stability and growth dynamics
Before turning to the discussion of the cubic equation (24) in the full parameter space, we first focus on the region near the boundary of stability, where the analysis is simplified by the fact that the real part of at least one ω j must be small.We proceed by discussing the two possible cases admitted by the solutions to a cubic equation with real coefficients, like eq. ( 24).
Case 1: All the three roots are purely real (but not necessarily distinct).The general solution given in eq. ( 25) consists in this case only of exponentially growing, decaying or constant contributions.In the stable region, ω j ≤ 0 for all j.Directly at the boundary to the unstable region, one must have ω j = 0 for at least one mode index, say j = 1.Setting ω 1 = 0 in eq. ( 24) readily yields As is shown below, the quantity B 1 defined here determines the boundary of stability.Inserting eq. ( 26) into eq.( 24), the other two decay rates result as For typical systems, one has In fact, for the constitutive relations reported in eqs.( 14) and ( 15), this inequality is violated only for unrealistically small shear rates γ 10 −12 and extreme densities ρ ρ m , where the hydrodynamic model considered here is doubtful.Since generally σ γ ≥ 0 and θ(k) > 0, it follows that ω 2,3 ≤ 0 -still assuming purely real ω j .Accordingly, provided that eq. ( 28) holds, none of the frequencies ω 2 and ω 3 vanishes and, consequently, the boundary of stability is solely defined by the condition ω 1 = 0 in this case.Close to the boundary of stability, nonlinear terms in ω 1 can be neglected in eq. ( 24), such that one readily obtains the growth rate (case 1, all frequencies real) ( We thus infer that, under the condition in eq. ( 28), the system is linearly stable if This inequality is consistent with the stability of the Navier-Stokes equations for a purely Newtonian fluid, since σ ρ = Π γ = 0 and thus B 1 > 0 in that case.As discussed below, eq. ( 30) in fact describes the boundary of stability of the whole relevant parameter space for the compressible single-component fluid.
In order for eq.( 25) to be real, one must have Re Â = ReA, Im Â = −ImA, with analogous conditions applying for B and C.These conditions are indeed fulfilled by the solution in eq. ( 25), which can be seen by writing eq. ( 21) as where M and M denote the real and imaginary parts of the matrix in eq. ( 21).Now let A = A + iA be a solution to eq. (31).Comparison of the real and imaginary parts of the resulting expression in eq.(31) shows that Case 2: One root is real and the other two are complex conjugates.Let ω 1 denote the purely real and ω 2,3 = Ω ±iΩ the complex conjugate solutions to eq. (24).For the imaginary part of eq. ( 25) to vanish, B and C must be complex conjugates of one another, while A must be purely real.Taking B = C * = B + iB allows one to write the general solution as Analogously to case 1, at least either ω 1 or Ω must vanish at the boundary of stability.If ω 1 = 0, we recover eq. ( 26) as a necessary consequence and eq.( 27) shows that Ω ≤ 0. Thus, in this case, the growing mode will be a monotonic function as in case 1 with a growth rate given by eq. ( 29).In contrast, the oscillatory modes will in general be decaying functions of time and will not give rise to any linear instability.
In order to analyze the case Ω = 0, we consider Vieta's formulas [28] for the solutions to the cubic equation (24) in case 2: If Ω = 0, eq.(33a) immediately implies ω 1 < 0, i.e., the purely real mode is stable.Moreover, combining the relations in eqs.(33a) to (33c) results in In order to determine the stability boundary for the complex conjugate pair of solution, we consider in eq. ( 24) small variations around Ω = 0. Accordingly, we insert ω = δΩ ± iΩ into eq.( 24), where Ω is determined by eq.(33b).Neglecting terms of O(δΩ 2 ) and higher in eq. ( 24) (keeping, however, all orders in Ω , as this quantity is not necessarily small), yields with Under the condition (28), the denominator on the r.h.s. in eq. ( 35) is positive for all k, allowing one to conclude that the system is linearly stable for As expected, the condition B 2 = 0 coincides with eq. ( 34).Note furthermore that B 2 > 0 for k → ∞ and, generally, A numerical analysis reveals that the condition B 2 | k=0 > 0 and thus B 2 (k) > 0 is fulfilled for all physically relevant ρ 0 and γ0 of the present model.The main result of the above analysis is that the cubic equation eq. ( 24) admits instability only through a single monotonically growing mode, the two other modes being decaying functions of time, either in a monotonic (case 1) or an oscillatory (case 2) fashion.

C. Stability diagram and discussion
As shown above, within the present linear stability analysis of the hydrodynamic equations for a compressible single component yield-stress fluid, the instability occurs uniquely via a monotonic (and thus non-oscillatory) mode, which grows exponentially with a rate given by eq. ( 29).In fact, an extensive numerical evaluation of the solutions of the dispersion relation in eq.(24) indicates that, over the whole relevant parameter space, all unstable modes have a non-oscillatory character.Owing to eq. ( 28) and the fact that B 1 ∼ k 2 > 0 in the limit k → ∞, the growth rate in eq. ( 29) becomes negative for sufficiently large k, as is necessary for a physically reasonable model.In particular, asymptotically for k → ∞ one obtains Note, however, that the continuum model in eq. ( 3) is not expected to be valid at arbitrarily small scales.We remark that, in the absence of the shear-and bulk-curvature viscosities κ and κ , one has ω(k → ∞) −ρ 0 B 1 /(bησ γ ).On the other hand, neglecting all contributions related to bulk viscosity, yields . From these results one infers that the stability of the system for large k is indeed due to the shear-curvature viscosity.Since σ γ is a growing function of k, the boundary of stability of the whole phase diagram of a bulk system is determined by eq. ( 30) for k = 0. Specifically, if B 1 (k = 0) < 0, the system is unstable for all wavenumbers, satisfying B 1 (k) = B 1 (k = 0) + κΠ ρ k 2 < 0 (cf.eq. ( 29)).In other words, all wavenumbers 0 < k < k c are unstable, where the critical wavenumber k c is defined via  Remarkably, this expression as well as the condition for stability threshold, B 1 (k) = 0, are identical to the corresponding expressions obtained in Ref. [16] for the advection-diffusion model described by eqs.( 1) and ( 2) [29].However, as can be inferred from eq. ( 29), owing to the presence of bulk viscosity, the fastest growing mode behaves differently as a function of k for the compressible fluid.Such a finite bulk viscosity is to be expected, since colloidal suspensions exhibit a certain degree of local compressibility even in the highly concentrated regime.Figure 2 shows the stability diagram obtained for k = 0.For illustrative purposes, it is more convenient to consider instead of B 1 [eq.(26)] the (dimensionless) stability parameter according to which the system is unstable for values S > 1.As seen in fig.2, the instability occurs only in the glassy phase (ρ > ρ g ).
The wavenumber k m and the growth rate ω m of the fastest growing mode has to be determined numerically from eq. ( 24) in the general case.Figure 3 shows k m and ω m as functions of the background density ρ 0 and shear rate γ0 .When expressed in terms of the fundamental time and length scales t 0 and a (which are unity for our choice of units), ω m and k m reach a maximum for intermediate shear rates and generally grow upon increasing the density.When taking instead the inverse shear rate as the fundamental time scale, ω m / γ0 grows with increasing distance from the boundary of stability [see main plot of fig.3(c)].At intermediate shear rates, an effective algebraic behavior ω m / γ0 ∼ γ−0.5 0 can be inferred from the numerics.As illustrated in figs.3(c) and 3(d), changing the value of the shear-curvature parameter κ 0 [eq.( 7)] has only a moderate effect on k m and ω m .
Close to the stability boundary, B 1 [eq.( 26)] and therefore k c are small, such that a Taylor expansion of the growth rate in eq.(29) to O(k 4 ) is sufficient to determine k m .(The leading term of the expansion of ω is of O(k 2 ).)Within this approximation, the wavenumber of the fastest growing mode follows by evaluating the condition dω/dk = 0 as with In the last expression, the fact that B (0) 1 0 close to the stability boundary has been used.Notably, due to bulk viscous effects, eq. ( 41) is generally different from the corresponding result obtained in Ref. [16].In fig.4, the typical behavior of the growth rate ω as a function of the wavenumber k is illustrated.We have chosen here values of the parameters ρ 0 and γ0 near the stability boundary, where eq. ( 41) provides an accurate approximation to the actual maximum wavenumber.While ω ∝ k 2 for small k [see eq. ( 29)], ω eventually becomes negative for sufficiently large k [see eq. ( 38)], as required for reasons of stability.Figure 5 illustrates the direction of growth and the magnitude of the most unstable mode (having wavenumber k m ).In order to obtain the amplitude vector v(ρ 0 , γ0 ), the nullspace solution v 0 = (ρ, γ, ūy ) of eq. ( 21) is determined and normalized, v 0 /||v 0 ||, and then projected onto the space spanned by ρ 0 and γ0 , additionally normalizing the components by ρ 0 and γ0 , respectively.As illustrated in fig.5, in general, ρ and γ have opposite signs in the unstable region, as expected for the SDC instability.A similar anti-correlation has been reported in molecular dynamics studies of heterogeneous flow in a hard sphere glass [6].Note that, with v, also −v is a valid solution of eq. ( 21); in the plot, we have chosen the positive sign of ρ.One notes that the development of the instability is dominated by a strong relative change of the shear rate, while the growth of the density is rather weak.This feature of the linear regime will also prevail in the nonlinear case discussed below.Comparing with figs.24)] as a function of the wavenumber k in the unstable region of the parameter space.For wavenumbers k with 0 < k < kc, the system is unstable.The dotted line indicates the location of the critical wavenumber kc [eq.( 39)].Near the boundary of stability, the wavenumber km of the fastest growth mode is estimated by eq. ( 41) (dash-dotted line).Near the boundary of stability and for small k, the growth rate ω is well approximated by eq. ( 29), implying ω ∝ k 2 .The values ρ0 = 0.91ρm, γ0 3.5 × 10 −4 , and κ0 = 100 are used for the calculation.
log ||v||  19).The coloring indicates the magnitude of v in a logarithmic scale, while the arrows indicate the growth direction (non-logarithmic scale along both axes), which is determined up to a sign.Accordingly, in the unstable region, fluctuations grow indefinitely by reducing the shear rate and increasing the density (or vice versa), providing a nonlinear feedback mechanism for the SDC-instability.

IV. NONLINEAR DYNAMICS AND STEADY STATES
A. Dynamics The one-dimensional Navier-Stokes equations for a compressible fluid given in eq. ( 5) are numerically solved in a slit geometry (see fig. 1) in the following way: the flux j = ρu is introduced and the partial differential equations are converted to ordinary ones by spatial discretization on a grid of L/∆h nodes [30,31].Specifically, we use second-order accurate central differences for the approximation of the spatial derivatives.The grid spacing is taken as ∆h = a, which is thus unity in our choice of units.A vanishing normal flux j y is assumed at the boundaries.The lateral flux j x at the boundaries is determined by imposing, at both walls, either a constant wall velocity u w , a constant wall shear rate γ w , or a constant wall stress σ w .Values of j x exterior to the computational domain ("ghost nodes") are calculated via linear interpolation from the adjacent bulk nodes [32].Exterior values of ρ are determined by assuming a vanishing gradient of ρ at the boundary.We have checked that the total density, L 0 dy ρ(y), remains practically constant during the time evolution.
As initial configuration we use a density and shear rate profile with a weak sinusoidal modulation (barely visible in the plots) in order to trigger the SDC instability.In the case of fixed wall stress, we initialize the shear rate with the constant value γ0 calculated from eq. ( 43) below.The dynamical evolution is, however, not significantly altered if instead a different value for the initial shear rate is used, except for a short transient at early times.After this transient, the evolution of the shear rate is found to be essentially enslaved to the density dynamics.In all cases, the wavelength of the maximally unstable mode predicted by the dispersion relation [eq.(24), see also fig.3(d)] is somewhat larger than the system size.Accordingly, the instability is realized here with the largest wavelength that fits in the system (cf.fig.4), provided that π/k c L, i.e., the system size exceeds half of the critical wavelength [see eq. ( 39)].This condition constrains, inter alia, also the value of the curvature parameter κ 0 [see eq. ( 7)], which sets the width of the shear band interface.Simulations with π/k m L are typically found to be unstable at late times since the nonlinear feedback mechanism leads to a singularity in the integration of the Navier-Stokes equations (see the discussion below eq.( 44)).This singularity manifests itself in a diverging viscosity and vanishing shear rate.For sufficiently large interface widths, global mass conservation stabilizes the stationary state before the singularity is reached.
Figures 6 to 8 illustrate the time evolution of the density ρ, flow velocity u x , and local shear rate γ = ∂ y u x across the slit in the unstable region for various boundary conditions.One observes that, in all cases, the system evolves from an essentially homogeneous initial state towards a steady state with inhomogeneous density and shear-rate profiles.The steady state of u x (or, correspondingly, the shear rate) is typically reached within a time scale 1/ γ0 determined by the average shear rate γ0 .The latter is given by (u x (L) − u x (0))/L = γav in the case where a fixed wall velocity is imposed and by γw in the case where a fixed wall shear rate is used.
At late times, the evolution slows down due to the slow transport of mass towards the boundaries.This effect is particularly pronounced in the case of a fixed wall velocity (fig.6), where the shear rate profile is essentially fully developed at times t γav 1, while the density at the left wall reaches the steady state only for times t γav O( 103 ).One observes that the time evolution is fastest if a fixed wall-shear rate is imposed.The broken left-right symmetry with respect to the walls in figs.6 and 8 is a direct consequence of the asymmetry of the initial configuration.In fact, using an initial sinusoidal density profile with a maximum in the right half of the system leads to spatially mirrored evolution.
The density dynamics is generally overdamped, which is expected based on an analysis of the linear equations in eq. ( 16): for typical values ρ 0 , γ0 of the density and shear rate in the unstable regime and for wavenumbers k ∼ O(10 −2 ) (in units of a), one finds for the present constitutive model [eq.( 14)] that the viscous term η∂ 2 y dominates over the restoring force Π ρ ∂ 2 y in the sound-wave equation [eq.(18)].An interesting question here regards the existence of multiple shear bands.In order to obtain such a structure, we decrease the shear-curvature viscosity κ, whereby, as stated above [see eq. ( 7)], the width of the interface between regions of low and high shear rates is reduced.As shown in Fig. 9, where we used a value of κ 0 = 2 (in dimensionless units), multiple shear bands are indeed observed in our model.The figure also highlights the important role of the initial perturbations for the formation of a multi-banded structure, since the number of nodes of the band directly depends on the period of the initial sinusoidal profile.The shear-curvature viscosity κ, in contrast, plays a subordinate role for detailed band structure.We emphasize that the profiles in fig. 9 are not in the steady state but correspond to times t γav ∼ O(0.1).In fact, the presently used constitutive model does not allow us to reach the steady state in these cases due to the occurrence of an intrinsic singularity [see eq. ( 44) below].

B. Steady states
In the steady state, the wall normal velocity must vanish, i.e., u y = 0.It thus follows from eqs. ( 5) and ( 8) that the density ρ and the shear rate γ in the steady state are determined by the equations Upon discretizing this boundary value problem using finite differences, the resulting system of nonlinear equations can be solved via Newton's method.In order for this scheme to converge, good initial guesses for ρ(y) and γ(y) are required, which can be obtained from the dynamical equation ( 5), as discussed above.Alternatively, the steady state profiles may also be directly obtained by integrating the PDEs in eq. ( 5) over a sufficiently large time, as is done in figs.6 to 8. According to eq. ( 42), in the steady state, the effective pressure Π − σ yield as well as the viscous stress σ xy must be constant throughout the system.The resulting profiles realizing these constraints are illustrated in figs.6 to 8 (thick black curves).For a fixed wall velocity u w or a fixed wall stress σ w , we obtain here spatially asymmetric steady state profiles for which the maximum density and minimum shear rate is attained close to the walls.These profiles are qualitatively similar to the ones observed in Ref. [16], where a fixed wall stress was considered (see also appendix B).The spatial symmetry of the steady profile in fig.7 is a consequence of the fact that the same shear rate γw is imposed at both walls.Profile shapes similar to the ones in figs.6 and 8 result when the values of γw at each wall are set accordingly (data not shown).
In the unstable region of the phase diagram, a necessary condition (which, however, is not sufficient; see below) for the development of an inhomogeneous steady state profile is the presence an initial perturbation in the system and a system size large enough such that at least one unstable mode can be accommodated.However, eq. ( 42) also admits constant solutions, i.e., ρ = ρ 0 and γ = γ0 .In this case, eq.(42b) with ∂ y γ = 0 readily yields a relation between ρ 0 , γ0 , and the stress in the system σ w (which arises as an integration constant and typically corresponds to the wall stress): where Φ = ρ 0 /ρ m and we used eq.( 14).Alternatively, eq.(42a) provides a relation between ρ 0 , γ0 , and the system (wall) pressure.If, instead of σ w , the wall shear rate γw = γ0 or the wall velocity u w (implying γ0 = u w /L) are prescribed, eq. ( 43) represents a family of solutions for ρ 0 with the integration constant σ w as adjustable parameter.Equation ( 43) is illustrated in fig.10.A solution of eq. ( 43) exists provided that the numerator on the r.h.s. is  5) for a fixed wall stress σw.The profiles shown are obtained at times t γ0 = 0.07, 0.42, 0.48, ∞, where t = ∞ corresponds to the steady state (thick black curve) reached for t γ0 10.An effective shear rate γ0 10 −3 is obtained from eq. ( 43) based in the initial density ρ0.In the steady state, the local shear rate at the boundary at y = 0 is obtained as γ 3.3 × 10 −5 .The initial growth rate (corresponding to the above γ0) of the maximally unstable mode is given by ωm 0.012 [see eq. ( 24 positive, i.e., for or, equivalently, if the yield stress remains below the system stress, For σ yield > σ w , instead, the system becomes increasingly more rigid and thus ceases to flow.This singular behavior is indeed reflected in the solution of the fluid dynamical equations [eq.( 5)] in the SDC-unstable region.In order to gain a heuristic understanding of this singularity, let us assume that, close to the inhomogeneous steady state, the density and shear rate profiles in the fluid consists of large nearly flat portions.In each of these regions then eq. ( 43) and thus eq. ( 44) approximately hold, with σ w being the corresponding local stress.In order to trigger the SDC instability, here typical values of σ w O(10) are required [33], for which eq. ( 44) implies Φ ∼ O(0.94) as an upper limit for the density (see fig. 10).As indicated in fig.5, once the system is unstable, some part of it evolves towards smaller shear rates and larger densities.In small systems, this growth is eventually limited by global mass conservation and the fact that the interface between low and high density regions must keep a certain width [see discussion after eq. ( 7)].In contrast, for systems much larger than the interface width, the limit for Φ implied by eq. ( 44) can be easily exceeded by means of mass transport, leading to σ yield > σ w and thus causing the solution of eq. ( 5) to become singular.The .Relation between shear rate γ0 and density ρ0 in a homogeneous system, as provided by eq. ( 43), for various values of the system stress σw (increasing in the direction of the arrow from σw = 8 to 12 in steps of 1).Equation ( 43) is well-defined only if condition (44) holds, i.e., for sufficiently small densities.The maximum possible density corresponds in the plot to the (σw-dependent) location where γ0 → 0. The dashed black curve represents the threshold of the SDC instability, B1 = 0 [see eq. ( 26) as well as fig.2].
singularity, in particular, makes an observation of stationary inhomogeneous profiles with multiple shear bands rather difficult.Nevertheless, as shown in fig.9, we do observe a multi-band structure up to the instant of the numerical singularity.Based on these findings and the above detailed analysis, we anticipate that a slightly modified version of the present constitutive model with a tamed singularity would exhibit stable multiple shear bands.Physically, the singularity can be understood as "freezing"-a behavior which is a hallmark of yield stress fluids upon increasing the density or decreasing the shear rate.Accordingly, a possibility to avoid this singularity would be to limit the growth of the yield stress and the viscosities in the constitutive equations ( 14) and ( 15).An adequate study of this issue is an interesting topic for future work.

V. SUMMARY
This study addresses the issue of flow heterogeneity, often observed in the glassy state of matter under externally imposed shear.We focus on the limit of gently sheared dense single-component fluids, such as colloidal hard sphere glasses, where hydrodynamic interactions are negligible.Therefore, the effect of solvent is ignored in this study.Instead of a coupling to the concentration field as in the original theory of shear-concentration coupling [13], in the present case, fluctuations of the velocity field are coupled to fluctuations of the fluid density.Analogously to a standard density-dependent thermodynamic pressure, here, a shear-rate dependent pressure drives a transverse flow away from regions of high shear rate, thus further lowering the viscosity in that region.Together with a strongly non-Newtonian viscous stress, this gives rise to a feed-back mechanism which ultimately determines the borders of flow stability.A detailed analysis of the resulting cubic dispersion relation reveals that the instability can occur only via a monotonic growth of fluctuations, thus excluding the possibility of an oscillatory growth mode.Notably, after an initial transient, the velocity and density fields reach a stationary profile.In this stationary state, regions of high shear rate exhibit low density and vise versa.For systems much larger than the characteristic width of the shear-band interface, the fluid model considered here generally develops a singularity in the unstable regime, accompanied by a vanishing shear rate and thus a divergent viscosity.The steady states obtained here in fact all occur for system sizes comparable to the interface width, in which case they are stabilized by means of global mass conservation.
Interestingly, the expression for the stability threshold and the range of unstable wavenumbers are identical to those obtained from an analysis of the advection-diffusion equation based on the original theory of shear-concentration coupling [16].The difference between the present compressible single-component fluid model and the one incorporating the coupling of flow to a concentration field is exhibited in the specific dynamics and underlying timescales, such as the expression for the fastest growing mode.In the compressible fluid, the density relaxes via overdamped sound waves-a transport mechanism which, in contrast to diffusion, leads to wavenumber-independent exponential relaxation in the limit of small frequencies and large wavenumbers.The use of a compressible fluid model is supported by molecular dynamics simulations [6], which show that variations of the density are a typically observed response to fluctuations of the shear rate in single-component hard-sphere colloidal glasses.
The existence of a stationary solution seems, at first sight, to be in conflict with the time dependent behavior of the shear band observed in molecular dynamics simulations [6].A plausible interpretation here would be to invoke the coupling between velocity fluctuations and structural heterogeneity in the glassy state.As also discussed in Ref. [5], this may lead to the formation of a locally depleted zone with a density in the stable regime and a denser packing in the remaining part of the system with enhanced instability and a corresponding temporal evolution.This is also in-line with recent reports on the strong influence of structural heterogeneity on plastic deformation in the amorphous solid state [34][35][36].Notably, the time scale of the shear band dynamics in MD simulations is of the order of the inverse shear rate [6].This is also the time associated with structural fluctuations, since during this time a particle moves a distance comparable to its size and the cage of nearest neighbors around it relaxes to a large extent.This stochastic effect is not included in the present deterministic model.A way to account for this would be to add a noise term into hydrodynamic equations, which is left for future work.

VI. ACKNOWLEDGMENTS
This work is supported by the German Research Foundation (DFG) under the project number VA 205/18-1.ICAMS acknowledges funding from its industrial sponsors, the state of North-Rhine Westphalia and the European Commission in the framework of the European Regional Development Fund (ERDF).where σ xy is given by eq. ( 8).In writing eq.(B1a), we used the fact that, for the constitutive model in eq. ( 15), the effective diffusivity D eff and the shear-gradient coefficient ξ defined in Ref. [16] derive from the pressure Π(ρ, γ) via D eff ≡ ∂ ρ Π and ξ ≡ ∂ γ Π, such that ∂ 2 y Π = ∂ y [D eff ∂ y ρ + ξ∂ y γ].Instead of the Couette geometry considered in Ref. [16], we study here the time evolution of eq.(B1) for a planar shear flow (see fig. 1).As in the main text, we impose either a fixed wall velocity u w , a fixed wall shear rate γw , or a fixed wall stress σ w at the boundaries.For the first two cases, instead of eq.(B1b), we solve the full time-dependent equation eq.(5b) for the velocity u x , with γ = ∂ y u x .We generally impose a vanishing pressure gradient ∂ y Π = 0 at the boundaries, which ensures global mass conservation for the dynamics described by eq.(B1a).Figures 11 to 13 illustrate the time evolution of the density and the shear rate in an unstable system for various boundary conditions.In all cases, the density profile is initialized with a sinusoidal modulation in order to trigger the initial instability.In general, the steady state is reached significantly faster for diffusive dynamics than for the compressible fluid model described by eq. ( 5).Both for fixed u w and σ w , the steady state is typically reached for strains t γav ∼ O(1) and ∼ O(0.1), respectively, while, for fixed γw , instead, the growth rate of the maximally unstable mode ω m provides a better estimate of the dynamical time scale than the strain.(The value of the time scale inferred from simulation depends somewhat on the chosen initial configuration.)The larger steady-state time scale in the compressible fluid model is predominantly caused by the slow transport of mass towards the wall at the late stages of the evolution.In fact, apart from this difference, the spatio-temporal evolution in both the compressible and the diffusive transport model are qualitatively similar (cf.figs.6 to 8).
Since we impose a vanishing pressure gradient at the boundaries when solving eq.(B1a), the steady states resulting from eqs. (42) and (B1) are characterized by Π and σ xy being constant throughout the system.Thus the only difference between the steady states of two models stems from the presence of the yield stress in eq.(42a).Accordingly, the steady state-profiles obtained from eq. (42) and eq.(B1) are very similar, which is illustrated in fig.14 for the case of fixed wall stress boundary conditions.It is therefore not surprising that the time evolution, when starting from the same initial conditions, is similar in the two models.

Figure 2 .
Figure 2. Stability diagram for a liquid (a) below (ρ < ρg) and (b) above (ρ > ρg) the glass transition.The glass transition occurs at a density of ρg/ρm 0.873, where ρm = 0.67 denotes the density of random close packing (in dimensionless units, see section II).The SDC instability occurs for values of the parameter S > 1 [eq.(40)] or, equivalently, for B1 < 0 [eq.(30)].The boundary of stability is indicated by the solid curve in (b), corresponding to S = 1.For comparison, the dashed curve in (b) represents the boundary of stability computed with σ yield yy = 0 in Πρ [eq.(11a)].
3(a) and 3(b), one infers that the growth amplitude ||v|| is largest in those those regions of the phase diagram where the growth rate ω m and the wavenumber k m are relatively small.

Figure 4 .
Figure 4. (c) Typical behavior of the growth rate ω [eq.(24)] as a function of the wavenumber k in the unstable region of the parameter space.For wavenumbers k with 0 < k < kc, the system is unstable.The dotted line indicates the location of the critical wavenumber kc [eq.(39)].Near the boundary of stability, the wavenumber km of the fastest growth mode is estimated by eq.(41) (dash-dotted line).Near the boundary of stability and for small k, the growth rate ω is well approximated by eq.(29), implying ω ∝ k 2 .The values ρ0 = 0.91ρm, γ0 3.5 × 10 −4 , and κ0 = 100 are used for the calculation.

Figure 5 .
Figure 5. Growth direction and magnitude of the fluctuation amplitude v ∝ (ρ/ρ0, γ/ γ0) (up to a normalization factor, see text) of the most unstable mode km, as determined by the nullspace solution of eq.(19).The coloring indicates the magnitude of v in a logarithmic scale, while the arrows indicate the growth direction (non-logarithmic scale along both axes), which is determined up to a sign.Accordingly, in the unstable region, fluctuations grow indefinitely by reducing the shear rate and increasing the density (or vice versa), providing a nonlinear feedback mechanism for the SDC-instability.

Figure 8 .
Figure 8.Time evolution of (a) the density, (b) the velocity, and (c) the shear rate resulting from eq. (5) for a fixed wall stress σw.The profiles shown are obtained at times t γ0 = 0.07, 0.42, 0.48, ∞, where t = ∞ corresponds to the steady state (thick black curve) reached for t γ010.An effective shear rate γ0 10 −3 is obtained from eq. (43) based in the initial density ρ0.In the steady state, the local shear rate at the boundary at y = 0 is obtained as γ 3.3 × 10 −5 .The initial growth rate (corresponding to the above γ0) of the maximally unstable mode is given by ωm 0.012 [see eq.(24)].Parameters L = 200∆h, κ0 = 100, ρ0 = 0.93ρm, and σw = 11 are used.

Figure 9 .
Figure 8.Time evolution of (a) the density, (b) the velocity, and (c) the shear rate resulting from eq. (5) for a fixed wall stress σw.The profiles shown are obtained at times t γ0 = 0.07, 0.42, 0.48, ∞, where t = ∞ corresponds to the steady state (thick black curve) reached for t γ010.An effective shear rate γ0 10 −3 is obtained from eq. (43) based in the initial density ρ0.In the steady state, the local shear rate at the boundary at y = 0 is obtained as γ 3.3 × 10 −5 .The initial growth rate (corresponding to the above γ0) of the maximally unstable mode is given by ωm 0.012 [see eq.(24)].Parameters L = 200∆h, κ0 = 100, ρ0 = 0.93ρm, and σw = 11 are used.
Figure 10.Relation between shear rate γ0 and density ρ0 in a homogeneous system, as provided by eq.(43), for various values of the system stress σw (increasing in the direction of the arrow from σw = 8 to 12 in steps of 1).Equation (43) is well-defined only if condition (44) holds, i.e., for sufficiently small densities.The maximum possible density corresponds in the plot to the (σw-dependent) location where γ0 → 0. The dashed black curve represents the threshold of the SDC instability, B1 = 0 [see eq.(26) as well as fig.2].