Shear banding in nematogenic fluids with oscillating orientational dynamics

We investigate the occurrence of shear banding in nematogenic fluids under planar Couette flow, based on mesoscopic dynamical equations for the orientational order parameter and the shear stress. We focus on parameter values where the sheared homogeneous system exhibits regular oscillatory orientational dynamics, whereas the equilibrium system is either isotropic (albeit close to the isotropic--nematic transition) or deep in its nematic phase. The numerical calculations are restricted to spatial variations in shear gradient direction. We find several new types of shear banded states characterized by regions with regular oscillatory orientational dynamics. In all cases shear banding is accompanied by a non--monotonicity of the flow curve of the homogeneous system; however, only in the case of the initially isotropic system this curve has the typical $S$--like shape. We also analyze the influence of different orientational boundary conditions and of the spatial correlation length.


Introduction
The emergence of banded structures in complex fluids under shear flow is a paradigmatic example of an instability in a correlated soft matter system far from equilibrium. Above a critical value of the applied shear rate (or shear stress), the formerly homogeneous system becomes unstable and separates into macroscopic bands with different local shear rates (stresses), see Refs. [1,2] for recent reviews. Typical systems where shear band formation has been observed experimentally are wormlike micelles [3], liquid-crystalline polymers [4], colloidal suspensions [5], but also non-ergodic soft systems such as glasses [6]. In all cases, the flow leads to reorganization of the fluid's microstructure which then feeds back into the flow field. This eventually leads to a non-monotonicity of the flow curve, that is, the relation between shear stress and shear rate. In that sense, non-monotonic flow curves are signatures of shear banding. Theoretically, shear banding (and the related vorticity banding [7]) has been studied mainly via continuum models. An important example is the diffusive (non-local) Johnson-Segelman (DJS) model [8,9] for shear thinning systems, i.e., systems in which the viscosity decreases with the stress, which form shear bands along the gradient direction. Moreover, particle resolved simulations [10,11] have revealed insight into microscopic mechanisms accompanying shear thinning and shear thickening systems.
In the present paper we focus on shear banding in nematogenic fluids whose anisotropic constituents can ar-range into orientationally ordered, yet transitionally disordered states. Prominent examples are wormlike micelles and colloidal suspensions of rod-like particles, both of which display flow-induced spatial instabilities in experiments. From the theoretical side, the shear-induced behavior of nematogenic fluids has been intensely studied on the basis of nonlinear equations for the dynamics of the orientational order parameter the so-called Q-tensor, with coupling to the concentration [12] or for dense systems deep in the nematic phase [13,14,15,16]. Contrary to the DJS model, the Q-tensor models allow to investigate directly the impact of shear on the structure (on a coarse-grained, order parameter level), from which the shear stress can then be derived by additional relations.
Already for homogeneous systems, these Q-tensor models predict many interesting effects such as the shearinduced shift of the isotropic-nematic transition [17], and the occurrence of dynamical states with regular or even chaotic oscillatory motion of the nematic director [13,15]. Indications of such time-dependent dynamical states under steady shear flow have also been observed in manyparticle simulations [10,11] and in experiments [18,19]. In addition, Q-tensor models have been used to explore spatial inhomogeneities, yielding shear banding between differently steady (aligned) states close to the isotropicnematic transition [12] and in parameter regimes where the dynamics of the homogeneous system is chaotic [20].
In the present article we consider (as in [20]) systems at constant concentration where deviations from the applied flow profile are taken into account in the Stokesian limit. Our purpose is to extend the earlier studies [21,22] on shear banding in nematogenic fluids along the gradient direction in several ways. First, we focus on parameters where the homogeneous systems exhibits regular oscillatory dynamics. These are, on the one hand, initially (i.e., atγ = 0) isotropic systems at low tumbling parameters and high shear rates, which display wagging motion [16] and, on the other hand, initially nematic systems where dynamical modes such as tumbling and kayaking are well established (see e.g. ref. [15]). Second, we consider both, homogeneous and inhomogeneous systems (with inhomogeneities in gradient direction) in order to identify the signatures of shear banding in the homogeneous flow curves. In fact, only one of the considered systems is characterized by a "van-der-Waals" like flow curve (which is familiar, e.g., within the DJS model [9]); the other ones rather exhibit discontinuities. Third, we explore systematically the impact of different orientational boundary conditions and different correlation lengths. We show, in particular, that appropriate boundaries can induce banding in otherwise homogeneous states. Further, we discuss the consequences for the stress of the system in the banded state.
The article is organized as follows. In sect. 2 we introduce the set of dynamical equations for spatially inhomogeneous, anisotropic fluids at constant density (with inhomogeneities along the direction of the shear gradient). Numerical results are presented in sect. 3. There, we first discuss (sect. 3.1) spatially homogeneous systems sheared from isotropic or nematic states and obtain the corresponding flow curves. Section 3.2 is then devoted to the appearance of shear bands and the impact of boundary conditions. This work finishes with concluding remarks and an outlook in sect. 4.

Theoretical framework
In the present work we consider systems of uniaxial, rigid rod-like particles (such as suspensions of fd-viruses [23]) whose orientation is characterized by the unit vector u i parallel to the symmetry axes of particle i. Following earlier studies within the Doi-Hess theory, the dynamics of the many-particle system is described by the space-and time-dependent tensorial order parameter Q(r, t) (thus, density variations are neglected) [24]. This second-rank Q-tensor is defined as where f (r, u, t) is the space-and time-dependent orientational distribution function (for a microscopic definition, see, e.g., [16]), and x stands for the symmetric traceless part of the tensor x. Explicitly, one has x µν = (x µν + x νµ )/2 − Tr(x)I µν /3, where µ and ν are Cartesian indices, I µν is the unit matrix, and Tr denotes the trace. Thus, the Q-tensor is, by definition, a symmetric traceless tensor. For the special case of uniaxial nematic phases, the Q-tensor reduces to the form Q = µ 3 3/2 nn , where n is the system-averaged nematic director, i.e., the eigenvector related to the largest eigenvalue µ 3 . In equilibrium, µ 3 is proportional to the well-known Maier-Saupe order parameter, µ 3 = 10/3S where S ≡ P 2 (u · n) and P 2 denotes the second Legendre polynomial [25].

Mesoscopic dynamical equations
In equilibrium, the orientational order of systems of rodlike particles is controlled by the temperature T (typical for molecular fluids) and/or by the number density (concentration) ρ; the latter case is characteristic of colloidal suspensions of fd-viruses [23]. On a mesoscopic (Q-tensor) level, the stability of homogeneous (isotropic or nematic) phases is governed by the Landau-type orientational free energy density where A, B and C are dimensionless coefficients. These can be related to system parameters such as ρ and the molecular aspect ratio κ as shown, e.g., in [26]. For spatially inhomogeneous phases, an additional contribution to the free energy arises due to the energy cost associated with local deformations of the alignment field. For uniaxial nematic order, an expansion of this elastic energy in terms of the director n was derived by Oseen [27] and Frank [28]. Rewriting the expression in terms of the full Q-tensor yields where ξ is the elastic correlation length, which is related to the pitch length of the Frank elastic theory [24,25]. On a microscopical level, ξ is related to the direct correlation function of the fluid [29,30,31]. The presence of shear strongly affects the overall orientational ordering due to the competition between flowinduced effects on individual molecular orientations (such as Jeffery orbits [32]) and the relaxation of the entire system towards equilibrium (governed by the free energy density). Moreover, in inhomogeneous systems (with spacedependent Q-tensor) the orientational dynamics feedback into the flow profile through the system's stress tensor. Within the Doi-Hess theory this interplay between flow (characterized through the velocity profile v(r)) and orientational ordering is described by a set of coupled nonlinear equations for Q and v which can be derived on the basis of non-equilibrium irreversible thermodynamics [33,34]. Disregarding non-convective flow (corresponding to fourth-order derivatives of the alignment), these equations are given by Mathematically, eq. (4) is a parabolic equation with H(Q, v) acting as a source term. From a physical point of view, it describes the dynamical evolution of the order parameter including a diffusive term ∝ ∇ 2 Q due to the elastic energy in eq. (3). The source term is given by [35,36] which describes the interplay between flow, entering via the vorticity Ω = 1/2(∇v T − ∇v) and the deformation rate Γ = 1/2(∇v T + ∇v), and relaxation towards equilibrium entering via the free energy derivative Equation (6) further involves the relaxational times τ qp and τ q [33]. As seen from eqs. (4) and (6), the ratio of these two times [appearing as a prefactor of Γ in eq. (6)] quantifies the perturbation of the system in the absence of orientational order. It is thus convenient to introduce, as a coupling parameter, the so-called tumbling parameter λ = −τ qp /τ q which is related to the molecular aspect (lengthto-breadth) ratio κ viz (see ref. [37]) Since λ depends only on the aspect ratio, it follows from eq. (8) that the MDHT is suitable to study spherical particles (λ = 0), disk-like particles (λ < 0) and rod-like particles (λ > 0). As stated before, our focus is on rod-like particles, thus, we consider positive values of the tumbling parameter.
The second mesoscopic equation (5) is the usual momentum balance equation [38]; it describes how the flow field changes due to spatial variations of the stress tensor T. We here consider a planar Couette flow (see fig. 1) where the fluid is confined between two infinitely extended, parallel plates (separated by a distance 2L along the ydirection) moving in opposite directions. The flow profile (for a Newtonian fluid) is then given by v(r) =γyê x , witḣ γ the shear rate.
The full stress tensor can be written as T tot = −p I + T asy + T, where the first term represents the (isotropic) hydrostatic pressure and the two other terms describe flow-induced effects [33]. Specifically, T asy includes asymmetric contributions, whereas T is a symmetric traceless tensor which can be written as According to eq. (9), T splits into a Newtonian contribution (already present in fluids with vanishing orientational order) and a contribution depending explicitly on the Qtensor, that is, In the present work we neglect the asymmetric part of the stress tensor (i.e., T asy = 0), since it typically relaxes faster to zero than the relevant hydrodynamic processes considered. Further, in Newtonian flow regimes this antisymmetric stress is zero anyway [25,35].

Explicit equations of motion
Equations (4) and (5) can be rewritten in terms of scaled variables; this is described in detail in the Appendix (see also [26]). One obtains where the tilde indicates scaled quantities and the parameter β appearing in eq. (12) is defined as This coefficient, or rather the ratio between β and the scaled viscosityη iso defines the Reynolds number of the solvent Experiments of shear-induced instabilities are typically performed at low Reynolds numbers, Re 1 [20,39]. In this limit the momentum balance equation (12) reduces to∇T = 0 .
We note that due to the time dependence ofQ(t), the total stressT evaluated through eqs. (9) and (10) generally also depends on time. However, at each time the total stress has to fulfill eq. (15). The resulting velocity profile [obtained by solving eq. (9) under the condition eq. (15)] can therefore deviate from the linear profile assumed initially. This is clearly essential for the description of spatial symmetry-breaking such as shear banding. In the following we drop the tilde (∼) on all variables; all quantities then appear with the same symbols as originally. We also set σ = 0 since this parameter has minor effect on the dynamics of the system for planar Couette flow geometry (see [40,41]). Regarding the spatial variation of the Q-tensor and the stress, we restrict ourselves to a one-dimensional investigation along the y-axis, i.e., the direction of the shear gradient (see fig. 1). Thus, we here exclude the possibility of banding in vorticity (z-) direction.
The resulting dynamical equations are simplified by expressing the Q-tensor in terms of a standard orthonormal tensor basis, that is, [42], where q 0 , q 1 and q 2 are related to ordering within the x-y-plane (i.e., the shear plane), whereas q 3 and q 4 refer to out-of-plane ordering [42]. From the orthogonality of the basis functions (B i : B j = δ ij ) it follows that q i = Q : B i , which allows to rewrite eqs. (11) into a set of scalar equations. Explicitly, one has In eqs. (16) the quantities Φ i are non-linear functions of the q i (stemming from the free-energy derivatives); they are given by where q 2 = 4 i=0 q 2 i . Finally, the momentum balance equation (12) becomes (in the regime of low Reynolds numbers) Equation (18) indicates that the only non-vanishing component of the stress tensor is T 2 , which corresponds to the in-plane stress T xy .
The gradient terms are discretized by a central finite difference scheme of fourth order. At the boundaries, asymmetric stencils (using only available grid points) are implemented [43]. The calculations are initialized with values of q 0 , · · · , q 4 matching the boundary conditions (see below), to accelerate the calculations we additionally use a small random perturbation. To find steady configurations of the system we monitor the evolution of q 0 , · · · , q 4 and T 2 , disregarding transient behavior. The resulting dynamical states are characterized employing an algorithm that recognizes the periodicity and sign change of the time dependent components q 0 , · · · , q 4 [14,13].
It turns out that the solution of eqs. (16)-(18) is quite sensitive to initial conditions; thus all calculations have been repeated several times. Further, we have checked that the steady-state solutions do not change with decreasing ∆t (however, the numerical stability does depend on the grid spacing).
Regarding the boundary conditions at the plates (y = ±L), we assume "strong anchoring" conditions, that is, the Q-tensor at the plates is constant in time, but may have different symmetries. We further assume that the degree of ordering is given by the corresponding equilibrium value.
Though these assumptions are clearly an idealization, we note that, from an experimental point of view, it is indeed possible to realize different boundary conditions by chemical or mechanical treatment of the plates [44,45,46,47]. Here we focus on the following cases (see fig. 2).
where µ eq 3 is the equilibrium value of the alignment tensor in the nematic phase. Equations (19) and (22) describe disordered states where the rod orientations are distributed, either in all three directions ("isotropic") or within the plane of the plates, i.e., in the x-z plane ("planar degenerate"). The other two boundary conditions correspond to nematic states, with the director lying either in the plane of the plates ("planar nematic") or along the gradient (y-) direction ("vertical nematic"). We note that the latter boundary condition is sometimes referred to as "homeotropic". Regarding the velocity field, we implement no-slip boundary conditions, that is v x (y = ±1) = ±1 (in reduced units).

Results and discussion
In this section we present numerical results for sheardriven systems whose equilibrium configuration (γ = 0) is either isotropic [characterized by Θ > 9/8 in the scaled free energy (27)] or nematic [Θ < 0]. We divide the discussion into two parts. In the first part (sect. 3.1) we focus on the spatially homogeneous system. Here we explore how the non-zero component of the stress tensor, T 2 , is affected by the temporal evolution of Q(r, t) and use this information to predict the formation of spatial instabilities. The second part (sect. 3.2) is devoted to the spatial-temporal behavior for initially isotropic and nematic systems, as well as to the impact of boundary conditions.

Homogeneous solutions
Here we consider spatially homogeneous systems where the boundaries do not play a role, corresponding to the limit of infinite plate separation, i.e., L → ∞. The scaled correlation length appearing in eq. (11) becomes zero and thus, all gradient terms in eqs. (16) and (18) vanish. In particular, the stress T 2 takes the form Here we set η iso = 1.0.

Homogeneous systems sheared from the isotropic state
We start by considering systems whose equilibrium state is isotropic (Θ = 1.20). Increasing the shear rate from  zero the system first develops a paranematic (PN), steady state characterized by weak (yet non-zero) nematic order; this is illustrated in fig. 3(a). The behavior upon further increase ofγ then depends on the tumbling parameter λ (recall that the latter is a measure of the aspect ratio). For λ 0.62 one observes a transition from paranematic to shear-aligned (A) state; the latter is also characterized by a time-independent director (as is the PN state), but the degree of ordering which can be quantified, e.g., by the norm ||Q|| = 4 i=0 q 2 i , is larger [34,42]. The PN-A transition is accompanied by a narrow region of bistability (not visible in fig. 3); in this regime the system's degree of ordering depends on the initial condition. This feature is reminiscent of the first-order isotropic-nematic transition in equilibrium. For smaller values of the tumbling parameter (λ 0.62) the A state is unstable; here the system develops wagging (W) motion characterized by regular oscillations of the nematic director within the shear plane. Overall, the behavior found in the present calculations agrees qualitatively with that reported in ref. [16]; the quantitative data for the boundary lines somewhat differ due to the different scaling of the free energy (see Appendix).
For each of the parameter sets (γ, λ) we have calculated the stress T 2 (in the oscillatory W state, we have averaged T 2 (t) over one period of time). Importantly, it turns out that different parameter sets can lead to the same value of T 2 . To illustrate this point, fig. 3(a) includes dashed (blue) lines and solid (green) lines corresponding to two constant values of T 2 . Moreover, there are several regions of λ where T 2 assumes the same value for different shear rates. For example, at λ = 0.55 there are three values oḟ γ with T 2 = 0.6, and for λ ≥ 0.7 one finds two solutions with T 2 = 0.4.
Given this multivalued behavior, it is interesting to consider corresponding flow curves T 2 (γ). Results for λ = 0.55 and λ = 1.25 are plotted in figs. (3(b) and (3(c), respectively. The flow curve for λ = 0.55 [see fig. 3(b)] displays a region with a negative slope (dT 2 /dγ < 0) betweenγ ≈ 0.357 andγ ≈ 0.365. Within this region the homogeneous flow is mechanically unstable, and as one might expect (and will be explicitly shown in sect. 3.2.1), the system forms a spatially inhomogeneous, shear banded state. We also note that the shear rateγ ≈ 0.359 characterized by T 2 = 0.6 and dT 2 /dγ < 0 agrees roughly with the corresponding point on the boundary line PN-W in fig. 3(a). This indicates that the orientational transition from the (steady) PN state to the (oscillatory) W state, on the one side, and the shear banding instability, on the other side, are closely intercorrelated.
In contrast, for λ = 1.25 [see fig. 3(c)] the flow curve does not display a region with negative slope. Rather one observes a discontinuity and, associated with this, hysteretic behavior. Upon increase ofγ from the small values, i.e., from the paranematic (PN) state, the systems discontinuously "jumps" to the aligned (A) state only aṫ γ ≈ 0.14 [which is above the upper blue dashed line in fig. 3(a)]. However, decreasingγ starting from the large shear rates characterizing the A state, the jump occurs at the much smaller shear rateγ ≈ 0.06. As we will show in sect. 3.2.2, the formation of shear bands in this case (λ = 1.25) strongly depends on the boundary conditions.

Homogeneous systems sheared from the nematic state
We now turn to systems which, in thermal equilibrium, are deep within the nematic phase. Specifically, we set Θ = −0.25 in eq. (17). Similar to previous work [15,16] we find that shear can induce a variety of time-dependent dynamical states [in addition to the W motion already appearing in initially isotropic systems, see fig. 3(a)], as well as a shear-aligned (A) steady state. An overview is given in fig. 4(a).
In the wagging and tumbling (T) state, the nematic director performs regular, oscillatory motion within the shear plane (q 4 (t) = q 5 (t) = 0 ∀t), whereas it displays outof-plane (yet regular) oscillations in the kayak-tumbling (KT) and kayak-wagging (KW) state (q i (t) = 0 ∀ i = 0, . . . , 4). Only in the A state the director stays constant in time. Note that, contrary to the case considered before (see fig. 3), there is no paranematic (PN) state at Θ = −0.25 since the system is deep in the nematic regime. We also note that earlier studies [15,16] investigating similar values of Θ have reported the occurrence of a region characterized by irregular and even chaotic motion of the director. This region is located around the point where the KT, KW and A states meet. Here we did not detect such a region because our algorithm does not resolve Lyapunov exponents.  We now turn to the resulting shear stress. The dashed (blue) and solid (green) lines in fig. 4(a) indicate parameter sets at which the orientational dynamics yield the constant stress-values T 2 = 9.0 and T 2 = 7.0, respectively. In the first case, the line provides a unique relation in the sense that an increase ofγ at fixed λ yields only one crossing with this line. This is different for the case T 2 = 7.0 where, depending on λ, one or two crossings can be observed. Exemplary flow curves for two values of the tumbling parameter are shown in the bottom parts of fig. 4. At λ ≈ 0.55 [see fig. 4(b)], where each of the constant-pressure lines is crossed only once, one observes a monotonic increase of T 2 withγ. In particular, there is no discontinuity or cusp even atγ ≈ 6.5, where the underlying orientational dynamics changes from out-of-plane kayaking-tumbling to in-plane tumbling. We note, however, that a systematic bifurcation analysis (such as the one in ref. [16], where a very similar system was considered) would presumably reveal a bistable region characterized by the presence of (at least) two attractors between the pure KT and the pure T state. In such a situation, the pressure T 2 would not be uniquely defined. This aspect certainly deserves more attention in a future study.
We now consider the case λ ≈ 1.25 in fig. 4(a), where an increase ofγ from small values yields two crossings with the constant-pressure curve T 2 = 7.0. As seen from fig. 4(c), the flow curve exhibits a pronounced discontinuity related to the transformation of the (out-of-plane) oscillating KT state into the shear-aligned (A) steady state atγ ≈ 3.67. One also recognizes a strong dependence on initial conditions (hysteresis), similar to the situation discussed in fig. 3(c). As we will later discuss in sect. 3.2.2, the initially nematic system at λ ≈ 1.25 indeed displays shear banding.

Spatiotemporal behavior and shear banding
In the preceding discussion we have found indications of the formation of inhomogeneous states in both, systems sheared from the isotropic and systems sheared from the nematic phase. We now analyze the corresponding systems (characterized by certain values of the tumbling parameter) further by calculating the full, spatiotemporal behavior of the Q-tensor and the resulting shear stress T 2 . To this end we have solved numerically eqs. (16)-(18) using the methodology described at the end of sect. 2.
Our discussion in this section is divided into two parts, covering the role of the two key factors impacting the spatial structure of the inhomogeneous systems. These are, first, the correlation length ξ, which appears as a prefactor of the gradient term in the orientational free energy density [see eq. (3)], and second, the boundary condition for Q at the plates [see eqs. (19)- (22)]. The impact of ξ is discussed in sect. 3.2.1, where we fix the boundary conditions according to the equilibrium configuration of the system. In sect. 3.2.2 we then explore the role of different boundary conditions.

Impact of the correlation length
Initially isotropic system We first consider the system at Θ = 1.20 and λ = 0.55, where we have observed a region of negative slope in the corresponding flow curve, T 2 (γ) [see fig. 3(b)]. Here we focus on a shear rate within this regime,γ = 0.365. In figs. 5(a)-(c) we show the spacetime evolution of the norm of the Q-tensor at three values of the correlation length. Because the equilibrium state is isotropic, we choose the boundary conditions according to eq. (19), that is, the boundaries do not support any orientational ordering.
Still, as seen from fig. 5(a), the system forms spatiotemporal structures with locally large values of ||Q|| already at the smallest correlation length considered, ξ = 10 −5 . Here, the observed pattern is rather "loose" with its width changing in time. We also find that, within the inhomogeneous regions, ||Q|| is oscillating in time. A closer analysis reveals that the oscillations are consistent with a wagging (W) state. Outside the inhomogeneous regions, ||Q|| takes values typical of a paranematic (PN) state. This behavior is, to some extent, expected since the value ofγ considered in fig. 5 is very close to the boundary between the PN and W state (see fig. 3).
Upon increasing the correlation length, ξ, we observe from figs. 5(b) and 5(c) that the regions characterized by W motion become more defined, both in terms of the shape of the emerging shear band, and in terms of the (increasingly regular) oscillatory motion of the order parameter. At the same time the interface between the W and PN region becomes broader. As a consequence, the W oscillations are transferred to some extent into the outer region, however, with a very small amplitude.
A further illustration of the presence of shear bands is plotted in fig. 6(a), where we present the local shear rate,γ(y), for the system at ξ 2 = 10 −3 . It is seen that the band in the middle of the system, where the orientational dynamics is of W type, is characterized by a significantly higher shear rate than the PN state close to the boundaries. To complete the picture, fig. 6(b) shows flow curves obtained for the inhomogeneous (initially isotropic) systems at different correlation length. Following previous studies [48] we have obtained these curves by calculating the mean value of T 2 increasing gradually from lower to larger values ofγ. As a reference, the corresponding homogeneous flow curve [see fig. 3(b)] is included. Interestingly, the value of T 2 corresponding to the banded state is essentially independent of ξ; in other words, the value of T 2 appears to be unique. This observation is consistent with previous studies on the basis of both, a Q-tensor model [12] and the DJS model [9]. We further observe from fig. 6(b) that there is a slight dependence on the value of T 2 on ξ at high shear rates beyond the banded state. This is an effect stemming from the inhomogeneities induced by the confining walls: the larger ξ, the larger is the extent of these inhomogeneities into the bulk-like region between the plates. In fact, by excluding these regions from the calculations, the results for T 2 completely agree for different ξ.
So far we have focused on an initially isotropic system at λ = 0.55. At the larger tumbling parameter λ = 1.25, where the homogeneous calculations [see fig. 3(c)] yield a discontinuous flow curve T 2 (γ) (rather than one with negative slope), the results from the spatially-resolved calculations are more complex. For isotropic boundary conditions we did not find shear banding behavior, regardless of the value of ξ. However, using boundary conditions which support nematic ordering (such as the planar alignment in eq. (21) or the vertical alignment in eq. (22)) we do find a well-defined shear band. This point will be further discussed in sect. 3.2.2. In all cases, one observes a clear spatial separation of the system into an inner band, where the orientational behavior corresponds to the kayak-tumbling (KT) state, and an outer region where the system is in a shear-aligned (A) state. Upon increase of ξ the width of the KT band widens, while the oscillations within the band become more and more regular.
Corresponding results for the local shear rate and the inhomogeneous flow curves are given in fig. 8. Compared to the initially isotropic system, we see from fig. 8(a) that the oscillatory (KT) band is characterized by a larger shear rate than the regions close to the boundaries. A further difference comes up when we consider in fig. 8(b) the values of the stress plateau in the inhomogeneous flow curve. Here we find a dependence on the correlation length; that is, the value of T 2 at the plateau increases with ξ. This contrasts with our corresponding results for the initially isotropic system (see the discussion of fig. 6). To which extent this dependency is subject to the initial conditions of the numerical calculations is a point which has remained elusive so far. equilibrium system. The correlation length is set to a constant value of ξ 2 = 10 −5 (higher values yield very similar results). Numerical results for the resulting flow curve of the inhomogeneous systems already discussed in sect. 3.2.1 are presented in fig. 9. We first consider systems sheared from the isotropic state at λ = 0.55, where we found a clear shear banding instability for fully isotropic boundary conditions [see figs. 5 and 6]. As indicated by the flow curves in fig. 9(a), similar behavior occurs for other (including nematic) boundary conditions. Indeed, as an analysis of the Q-tensor reveals, all of the systems form bands (within a range of shear ratesγ ≈ 0.355 -γ = 0.385) with wagginglike oscillations within paranematic regimes at the plates. Outside the banding region, the systems are characterized by the same value of T 2 .

Role of the boundary conditions
Moreover, the value of the "selected" stress within the banding region, T sel 2 ≈ 0.59, is essentially independent of the boundary conditions. The latter only affect the onset of shear banding upon increasingγ from low shear rates. Specifically, for the two types of isotropic boundary conditions [eqs. (19) and (22)], as well as for nematic ordering within the plane of the plates [eq. (20)], the system stays in the homogeneous paranematic state for allγ up to the maximum of the flow curve. In contrast, the system with vertical alignment at the plates, i.e., alignment in the shear gradient (y-) direction [see eq. (21)], forms bands once the value of T sel 2 is reached (atγ ≈ 0.34). In this sense, the vertical nematic ordering favors the occurrence of the W state characterized by oscillations in the flow-gradient plane. For completeness we also note that, upon decreasingγ from high values, the system goes without hysteresis into the banded state (with the same stress), irrespective of the boundary conditions.
We now turn to the initially nematic case [see fig. 9(b)]. Upon increasingγ from lower values all systems, irrespective of boundary conditions, display shear band formation at shear rates in the rangeγ ≈ 3.0 -γ ≈ 4.5; here they break up into a band with kayaking-tumbling dynamics (i.e., oscillations out of the shear plane) surrounded by regions of shear-alignment at the boundaries [see fig. 7 for results with planar nematic boundaries]. Further, the stress T 2 characterizing the banded state seems to be unique, and the boundary conditions only affect the onset of shear banding (upon starting from the low-shear rate branch, where the system is in the KT state). Specifically, we see from fig. 9(b) that the onset of banding is "delayed" when we use planar degenerate boundaries [see eq. (22)]. These boundary conditions seem to support the KT state, which is understandable as the KT oscillations are out of the shear plane and thus, do involve the plane of the plates. Interestingly, the behavior upon decreasing the shear rate from the aligned state is different: in that case, all systems stay in the aligned state until the lower end of the high-shear rate branch is reached; then they directly jump into the KT state without an intermediate shear banded state.
To summarize, both systems considered in fig. 9 display shear banding irrespective of the detailed nature of the boundary conditions, if the shear rate is increased from low values. The boundary conditions then only influence the "critical" shear rate at which the homogeneous state observed at smallγ breaks up into bands. On the contrary, shear banding upon decreasingγ from high values is only seen in the initially isotropic system. The behavior in the initially nematic system thus depends on the initial conditions, similar to what has been found in the DJS model [49].
Finally, we come back to a system which we already considered briefly in sect. 3.1.1, that is, the case Θ = 1.20 and λ = 1.25. This initially isotropic system is characterized by a discontinuous homogeneous flow curve [related to the transition from paranematic to shear-aligned state, see fig. 3(c)], but does not form clear shear bands when using disordered boundary conditions [see eqs. (19) or (22)]. Interestingly, this changes when we use "nematic" boundary conditions with alignment either in the plane of the plates [eq. (20)] or along the vorticity direction [eq. (21)].
As an illustration, we present in fig. 10(a) the spatiotemporal behavior of the Q-tensor at ξ 2 = 10 −5 , with planar boundaries and an imposed shear rateγ = 0.115. The plot reveals a band with paranematic ordering within the outer regions where the system is shear-aligned. Figure 10

Concluding remarks
In this paper we have investigated the occurrence of shear bands in nematogenic fluids based on the mesoscopic Doi-Hess theory for the orientational order parameter tensor Q coupled to the shear stress T. We have focused on instabilities along the gradient direction (gradient banding). Studying sheared systems with different (isotropic or nematic) equilibrium states, different tumbling parameters (i.e., aspect ratios), and different orientational boundary conditions, our results reveal complex shear banding behavior whose characteristics can be significantly different from that seen in other models, such as the DJS model (where the dynamical variable is the shear stress alone). In most cases considered, the shear bands involve oscillatory (but not chaotic) orientational motion in certain regions of space. In terms of parameters, our study extents earlier investigations focusing on the isotropic-nematic transition [12] or the rheochaotic regime [21,22]. "Classical" behavior characterized by the S-shaped flow curve (obtained from the homogeneous solutions) and a unique value of the stress in the banded regime (such as in the DJS model) is found only in one case, namely an initially isotropic system with relatively small tumbling parameter. This system displays bands with wagging-like motion in the inner part and steady alignment close to the boundaries. In all other cases, the homogeneous flow curves display a discontinuity (rather than the S-shape). The observed band formation then depends on the orientational boundary conditions in the sense that certain boundary conditions can support or hinder the formation of shear bands. Moreover, in one case we found a strong dependence on the pathway of the shear protocol. These observations suggest that the mechanisms of shear band-ing and stress selection are more complex than in simpler (such as the DJS) models [9,50].
An important question arising from the present work concerns the relation between the parameter sets (θ, λ,γ) considered here and the system's "phase" diagram under shear. In many nematogenic systems the relevant variable for the isotropic-nematic transition is the concentration, which does not occur explicitly in our approach. However, there is an implicit concentration dependence through the free energy functional (specifically, the prefactor θ of the quadratic term). Based on that dependence, we have proposed in an earlier study [16] a phase diagram in the concentration-shear rate plane, which qualitatively resembles earlier results [2,12]. In particular, the diagram reproduces the shear-induced shift of the isotropicnematic transition towards smaller concentrations, as well as a critical shear rate beyond which the transition becomes continuous. Based on ref. [16] we find that the initially isotropic state with small tumbling parameter (with θ = 1.20, λ = 0.55) considered here is very close to the critical point (see fig. 7 in [16]). The appearance of shear banding (in gradient direction) in this situation is indeed consistent with earlier predictions for systems of colloidal rods, such as suspensions of f d-viruses [2,23]. The other parameter sets considered in the present work lie deep in the nematic phase of the concentration-shear rate phase diagram.
A further interesting point concerns the occurrence of shear bands in vorticity direction, a scenario which we have implicitly ruled out by focusing on inhomogeneities along the y-direction alone. Indeed, vorticity banding has been predicted to occur in colloidal rod (f d-virus) suspensions at conditions within the isotropic-nematic spinodal under shear [2]. Conceptually, vorticity-banded states are characterized by the same shear rate (rather than same shear stress as gradient-banded states) [51]. With this background, the shape of the (homogeneous) flow curves obtained here (namely those characterized by discontinuities and large hysteresis) may be taken as an indication for vorticity banding. This point certainly warrants further investigation. Another important direction is the investigation of (binary) mixture system which have even more complex dynamics already in the homogeneous case [26]. Work in these directions is in progress.