Time-delayed feedback control of shear-driven micellar systems
Benjamin von Lospichl Sabine H. L. Klapp
physics.flu-dyn cond-mat.soft








 blospichl@itp.tu-berlin.de
[present address: ]Stranski-Laboratorium für Physikalische und Theoretische Chemie, Institut für Chemie, Sekretariat TC 7, Technische Universität Berlin, Straße des 17. Juni 124, D-10623 Berlin, Germany
Institut für Theoretische Physik, Sekretariat EW 7-1, Technische Universität Berlin, Hardenbergstraße 36, D-10623 Berlin, Germany


 klapp@physik.tu-berlin.de
Institut für Theoretische Physik, Sekretariat EW 7-1, Technische Universität Berlin, Hardenbergstraße 36, D-10623 Berlin, Germany



             


Suspensions of elongated micelles under shear display complex non-linear behaviour including shear banding, spatio-temporal oscillatory patterns and chaotic response. Based on a suitable rheological model [S. M. Fielding and P. D. Olmsted, Phys. Rev. Lett. 92, 084502 (2004)], we here explore possibilities to manipulate the dynamical behaviour via closed-loop (feedback) control involving a time delay τ. The model considered relates the viscoelastic stress of the system to a structural variable, that is, the length of the micelles, yielding two time- and space-dependent dynamical variables ξ_1, ξ_2. As a starting point we perform a systematic linear stability analysis of the uncontrolled system for (i) an externally imposed average shear rate and (ii) an imposed total stress, and compare the results to those from extensive numerical simulations. We then apply the so-called Pyragas feedback scheme where the equations of motion are supplemented by a control term of the form K(a(t)-a(t-τ)) with a being a measurable quantity depending on the rheological protocol. For the choice of an imposed shear rate, the Pyragas scheme for the total stress reduces to a non-diagonal scheme concentrating on the viscoelastic stress. Focusing on parameters close to a Hopf bifurcation, where the uncontrolled system displays oscillatory states as well as hysteresis in the shear rate controlled protocol, we demonstrate that (local) Pyragas control leads to a full stabilization to the steady state solution of the total stress, while a global control scheme does not work. In contrast, for the case of imposed total stress, global Pyragas control fully stabilizes the system. In both cases, the control does not change the space of solutions, rather it selects the steady state solutions out of the existing solutions. This underlines the non-invasive character of the Pyragas scheme.


47.57.Qk, 05.45.Gg, 47.52.+j, 02.30.Yy
		
		
		

Time-delayed feedback control of shear-driven micellar systems
    Sabine H. L. Klapp
    December 30, 2023
==============================================================







§ INTRODUCTION


Exposing soft matter systems to shear flow often results in a non-linear behavior of measurable rheological quantities (such as viscosity, stress) and to instabilities. The most prominent example is the shear banding instability, where the system splits into two or multiple spatial regions (“bands”) of different local shear rate and local viscosity <cit.>. This instability is indicated by a non-monotonicity of the flow curve, that is, the relationship between the total stress  in the steady state and the imposed shear rate γ̇. Shear banding occurs in many complex fluids like polymer solutions, emulsions or liquid crystals. To capture this phenomenon theoretically, several non-linear models have been proposed revealing multiple branches of the underlying constitutive curve (γ̇). A classical example is the diffusive Johnson-Segalman model <cit.>.

In the present article we concentrate on one of the most prominent examples of a shear-banding system, namely solutions of elongated or wormlike micelles. Within the past decades the self-assembly processes as well as the rheological properties of such systems have attracted a lot of interest <cit.>. Due to their flexibility, on the one hand, and their ability to break and reform, on the other hand, wormlike micelles often called “living polymers”. Further, wormlike micellar solutions exhibit an unique viscoelastic behaviour which can be exploited for many applications such as daily care products <cit.> or as fracturing fluids in oil recovery <cit.>. Indeed, already the first rheological experiments done by Rehage and co-workers demonstrated that the viscosity of these elongated micelles reveals a highly non-linear flow behaviour <cit.>: Unlike covalently bound polymers, wormlike micelles can not only break but also recombine after exposing them to high shear forces. The first theory to account for the resulting viscoelastic spectra of wormlike micelles has been proposed by Turner and Cates <cit.>.

Beyond stationary shear bands and nonlinear viscosity, wormlike micellar solutions can moreover display oscillatory states and irregular or even chaotic spatio-temporal patterns, as indicated based on a bunch of experimental data based on particle image velocimetry and bulk-rheology <cit.>. Theoretical approaches to account for this complex behavior in micellar (and other complex) fluids may be subdivided into three classes. The first are spatially homogeneous models for the rheological (mechanical) dynamical variables such as  (or the viscoelastic stress σ) and γ̇; this approach has been applied to shear-thickening fluids <cit.>. The same variables are also at the center of several inhomogeneous or non-local models for shear-thinning <cit.> and shear-thickening fluids <cit.>. These models additionally take into account that the macroscopic rheological behaviour is coupled to a intrinsic (microscopic) variable such as the polymer concentration <cit.>, the micellar length <cit.>, or the director of orientational ordering. The third class focuses directly on the (spatio-temporal) dynamics of liquid-crystalline substances, involving evolution equations for the (tensorial) nematic order parameter and its backcoupling to the flow <cit.>.

The present study starts from a (non-local) model for micellar solutions proposed by Fielding and Olmsted <cit.>, where the viscoelastic stress σ(y,t) is coupled to the length n(y,t) of the elongated micelles. This coupling of a mechanical and a structural variables destabilizes the high shear rate branch of the flow curve, giving rise to spatio-temporal oscillatory and chaotic solutions beyond a Hopf bifurcation <cit.>. The goal of the present paper is twofold: First, we systematically examine the system's dynamics by a linear stability analysis for two rheological protocols (fixed average γ̇ or ), thereby extending the results presented in Ref. <cit.>. Second, we explore methods to manipulate the occurring dynamical states by using time-delayed feedback control.

Indeed, in many technical applications of wormlike micelles the occurrence of oscillatory or chaotic states is not desired. Motivated by this, we here seek to stabilize stationary states inside the parameter range where the uncontrolled system displays oscillations. Specifically, we apply the concept of time-delayed feedback control (TDFC) as proposed by Pyragas in 1992 <cit.>. Within the Pyragas control scheme, the dynamical equations for the relevant variables a (which, ideally, should be measurable quantities) are supplemented by control terms involving the difference a(t)-a(t-τ), with τ being the delay time. This type of control is noninvasive as the control forces vanish when the steady state (or a periodic state with period T_osc=mτ with m =1,2,…) is reached. The Pyragas scheme has been applied to a broad variety of non-linear systems from various fields (see <cit.> for overviews) including semiconductor nanostructures, lasers, neural systems and general reaction-diffusion systems <cit.>. There are also recent applications to the flow of soft-matter systems, both from the experimental <cit.> and from the theoretical side <cit.>, an example being the stabilization of steady shear-aligned states in sheared liquid crystals <cit.>. Within the present study, we focus on a non-diagonal control scheme which targets the average total stress or shear rate, respectively, yielding a Pyragas term involving only the viscoelastic stress. We apply this scheme to states close to a Hopf bifurcation, where the system displays hysteresis and the relation between shear rate and total stress is not unique. The control then selects the stationary solution.

Our publication is structured as follows. In Sec. <ref> we review the model equations proposed in Ref. <cit.>. In Sec. <ref> we describe the linear stability analysis for the two rheological protocols and compare the results with those from extensive numerical calculations of the full system. In Sec. <ref> we introduce a time-delayed feedback control scheme to both rheological protocols. We then use the results of linear stability analysis derived in Sec. <ref> to find the neutral stability curves determining the parameter ranges where the control scheme is successful. The article closes with a brief summary and outlook.



§ MODEL EQUATIONS AND BACKGROUND




 §.§ Flow geometry and force balance equation


We assume our micellar system to be subject to a planar Couette geometry consisting of two plates separated by a distance L along the y-direction. The lower plate is fixed at the position y=0, while the upper plate located at the distance y=L moves with a constant velocity in x-direction. This is equivalent to an externally imposed velocity field ^0=v_x^0(y,t), with  being the unit vector in x-direction and v_x^0(y,t)=y γ̇(t). Here, γ̇(t) is the external shear rate corresponding to ^0. The local shear rate can be defined as spatial gradient of the actual velocity profile v_x(y,t) in the system (which may differ from v_x^0(y,t)), that is, γ̇(y,t)=∂_y v_x(y,t). In general, the total stress tensor  of a viscoelastic fluid is defined by <cit.>

    =+2η -p  ,


where  is the viscoelastic contribution which may be related to a structual quantity. The second term on the r.h.s. of Eq. (<ref>) denotes the Newtonian part, with η being the Newtonian viscosity and 2=∇ +(∇ )^T the deformation tensor. The third term in Eq. (<ref>) is an isotropic contribution, where p is the hydrostatic pressure and  the identity matrix. 

The dynamics of a viscoelastic fluid can be described through the generalized Navier-Stokes equation (or momentum balance equation)

    ρ (∂_t+ ∇)     =∇ 
        =∇ (+2η -p ) ,


where ρ is the fluid density. We assume our system to be incompressible, i.e. ∇·=0. Additionally, the system is assumed to be overdamped (limit of low Reynolds numbers). In this case the momentum balance equation (<ref>) reduces to the force balance equation

    ∇ (+2η -p )=0 .


For the present imposed linear velocity field, Eq. (<ref>) reduces to

    ∂_y (σ(y,t)+η γ̇(y,t))=0 ,


where σ(y,t) is the local value of the xy-component of the viscoelastic stress tensor , and γ̇(y,t) is the local shear rate. The restriction to one spatial dimension is motivated by our flow geometry. The force balance equation (<ref>) demands the xy-component (y,t) of the stress tensor  to be spatially constant. Thus

    (y,t)=σ(y,t)+η γ̇(y,t)≡(t)


for all positions y and times t. Note that the constancy in space of (t) does not necessarily hold for the individual term σ(y,t) and γ̇(y,t). Averaging the second and third part of Eq. (<ref>) over the y-coordinate yields

    σ(t)+η γ̇(t)=(t) .
 

where A(t)=⟨ A(y,t)⟩_y=1/L ∫_0^L A(y,t) dy defines the spatial average over the y-direction of a quantity A(y,t), which can be either the viscoelastic stress or the shear rate.

In the present study we consider two rheological protocols: First, we impose a temporarily constant average shear rate γ̇(t)=γ̇=const. to the system, yielding the local viscoelastic stress σ(y,t) as the main dynamical variable. Equating the right hand sides of Eqs. (<ref>) and (<ref>) yields an expression for the local shear rate γ̇(y,t), that is

    γ̇(y,t)=γ̇+1η [σ(t)-σ(y,t)] .


The total stress in this shear rate controlled protocol is determined from Eq. (<ref>) or (<ref>). Secondly, we impose a temporarily constant stress to the system, i.e., (t)=^fix=const. From Eq. (<ref>) it then follows that

    γ̇(y,t)=1η [^fix-σ(y,t)] .




 §.§ Model equations in the shear rate controlled protocol


In the present work we employ the rheological model proposed in Ref. <cit.>, which couples the viscoelastic stress σ(y,t) to a variable characterizing the microscopic structure. The latter is given by the length distribution n(y,t) of elongated micelles.

The dynamics of n(y,t) in the shear rate controlled protocol is governed by the non-linear relaxation equation

    ∂_t n(y,t)=-n(y,t)τ_n+1τ_n (n_01+[τ_n γ̇(y,t)]^β) ,


where τ_n is the relaxation time of the underlying rates of scission and recombination of the micelles, n_0>0 is the (average) length at zero shear and β≥0 is a power law exponent. While the first term on the right hand side of Eq. (<ref>) is a pure relaxation term, the second term accounts for shear-induced changes of scission or recombination processes. The local shear rate γ̇(y,t) is determined through Eq. (<ref>).

The homogeneous steady state solution (∂_t n(y,t)=0) of Eq. (<ref>) reads

    n_s(γ̇)=n_01+(τ_n γ̇)^β .


In the low shear rate limit, γ̇→0, the steady state solution converges to n_0, while in the limit of (infinitely) high shear rates (γ̇→∞) we find n_s→0. The rate of convergence in this high shear rate limit strongly depends on the power law exponent β as illustrated in Fig. <ref> (a). In the case β=0 the steady state solution will be n_s=n_0/2=const. for all shear rates. The dependence of the function n_s(γ̇) on the relaxation time τ_n is shown in Fig. <ref> (b). We further note that there is a critical shear rate γ̇_c=τ_n^-1, at which the zero-shear length n_0 is halved upon increase of γ̇. Generally, the decrease of the micellar length with increasing shear rates reflects that scission processes become dominant. A similar behaviour has been observed in particle-based (Molecular Dynamics) computer simulations <cit.>. Experimentally, the micellar length is hardly accessible.



The spatio-temporal dynamics of the viscoelastic stress σ(y,t) is determined by the partial differential equation

    ∂_t σ(y,t)=     -σ(y,t)τ(n)+γ̇(y,t)1+[τ(n) γ̇(y,t)]^2
        +𝒟 ∂_y^2 σ(y,t) ,


with 𝒟 being the (stress-)diffusion constant and τ(n) a length dependent relaxation time [not to be confused with the time τ_n appearing in Eq. (<ref>)]. The quantity τ(n) is defined as

    τ(n)=τ_0 (n(y,t)n_0)^α .


Here, τ_0>0 is the zero-shear relaxation time of the micelles, and α≥0 is another power law exponent. The system of Eqs. (<ref>), (<ref>) and (<ref>) is closed by Eq. (<ref>) for γ̇(y,t).

The homogeneous steady state solution of Eq. (<ref>) (i.e., ∂_t σ(y,t)=0) in the non-diffusive limit (𝒟=0) takes the following form:

    σ_s(γ̇)=τ(n_s) γ̇1+(τ(n_s) γ̇)^2 .


Here we have to distinguish several limiting cases: For α=0, the relaxation time reduces to a constant (τ(n_s)=τ_0), and therefore the coupling to the structural variable n(y,t) vanishes. This case has been discussed, e.g., in Ref. <cit.>. The steady state solution σ_s is then equivalent to the steady state solution found for the xy-component of the viscoelastic stress tensor in the Johnson-Segalman model <cit.>. For 0<α<1, we find n_0^α τ(n_s)/τ_0>n_s^α for all shear rates. On the other hand, n_0^α τ(n_s)/τ_0<n_s^α for α>1 for all shear rates. Finally, in the special case α=1, the relaxation time follows instantaneously the micellar length, that is τ(n_s)∝ n_s. The plausibility of the ansatz made in Eq. (<ref>) is underlined by results from Molecular Dynamics simulations showing that the breaking time, which is similar to the relaxation time defined in Eq. (<ref>), follows a power law like behaviour <cit.>.

Using Eqs. (<ref>) and (<ref>) we can state an explicit expression for the total stress in the homogeneous stationary limit,

    _s(γ̇)=σ_s(γ̇)+η γ̇=τ(n_s) γ̇1+(τ(n_s) γ̇)^2+η γ̇ .




From the Johnson-Segalman model it is known that σ_s as function of γ̇ displays a maximum at low shear rates and then converges towards zero for increasing shear rates. For the present model, the behaviour of σ_s as function of γ̇ strongly depends on the choice of the power law exponent α: For α<0.5, we find the same behaviour as in the Johnson-Segalman model. However, for α≥0.5 this behaviour changes such that the maximum at low shear rate is followed by two other extrema, that is one minimum and one maximum, before tending to zero at high shear rates. This non-monotonic behaviour reflects the highly non-linear character of the model equations (<ref>) and (<ref>). To get a more detailed insight into the impact of changing the relaxation time τ_n, which is attributed to scission and recombination processes, we plot in Fig. <ref> the steady state solution of the total stress, _s(γ̇), for various τ_n. Upon increase of τ_n the maximum of the constitutive curve stays at the same position, while the minimum of the curve shifts to smaller shear rates. A further effect of increasing τ_n is a change of the slope of the high shear rate branch. This is again an indication for a stronger coupling to the (micro-)structure.



 §.§ Model equations in the stress controlled protocol


The dynamical equations corresponding to the stress controlled protocol ((t)=^fix=const.) follow from Eqs. (<ref>) and (<ref>) by replacing the shear rate via Eq. (<ref>), that is, γ̇(y,t)=1/η [^fix-σ(y,t)]. This leads to

    ∂_t n(y,t)=     -n(y,t)τ_n
        +1τ_n (n_01+[τ_n η^-1 (^fix-σ(y,t))]^β)


and

    ∂_t σ(y,t)=     -σ(y,t)τ(n)+(^fix-σ(y,t))η+[τ(n)·(^fix-σ(y,t))]^2
        +𝒟 ∂_y^2 σ(y,t) ,


with

    τ(n)=τ_0 (n(y,t)n_0)^α .


The system of equations (<ref>) to (<ref>) is closed, i.e., there is no further equation required to calculate the dynamics of the system. 

The stress controlled protocol has already been discussed in Ref. <cit.> for the non-diffusive limit (𝒟=0), where the dynamical variables are only time-dependent. To evaluate the impact of a finite diffusion constant (𝒟>0) within the stress controlled protocol, we have performed numerical test calculations of Eqs. (<ref>) – (<ref>) using the procedure described in Sec. <ref> (which allows for inhomogeneous solutions). We found the system to be spatially constant in all cases considered. Therefore, we set 𝒟=0 in all following calculations related to the stress controlled protocol, yielding a spatially homogeneous dynamics.



§ LINEAR STABILITY ANALYSIS AND STATE DIAGRAMS
 

In order to get a first insight into the (complex) dynamics of our model system, we perform a linear stability analysis. We follow the strategy outlined in Ref. <cit.>. We then compare the results with those from a full numerical solution. As we will demonstrate, the overall dynamics strongly depends on the rheological protocol. We therefore discuss the two protocols separately.



 §.§ Shear rate controlled protocol
 



  §.§.§ Linear stability analysis
 

We start by introducing a two-dimensional vector containing the dynamical variables n(y,t) and σ(y,t),

    ξ(y,t)=[n(y,t),σ(y,t)]^T .


Within this vector notation we can rewrite Eqs. (<ref>) and (<ref>) as

    ∂_t ξ(y,t)=Q(ξ,γ̇)+∂_y^2 (𝒟 ξ(y,t)) ,


where Q(ξ,γ̇) contains the non-linear terms on the r.h.s. of Eqs. (<ref>) and (<ref>). Further, we have introduced the diffusion matrix 𝒟 which is given by

    𝒟=[ 0 0; 0 𝒟 ] .


The homogeneous steady state solutions of the dynamical system Eq. (<ref>) are given in Eqs. (<ref>) and (<ref>), respectively. We summarize these solutions in a vector ξ_s=[n_s,σ_s]^T. To determine whether ξ_s is stable or unstable, we add small heterogeneous perturbations δξ(y,t), yielding

    ξ(y,t)=ξ_s+δξ(y,t) .


Due to the coupling between shear rate and viscoelastic stress given in Eq. (<ref>), we additionally have to consider heterogeneous perturbations δγ̇(y,t) in the shear rate

    γ̇(y,t)=γ̇+δγ̇(y,t) .


Substituting Eqs. (<ref>) and (<ref>) into Eq. (<ref>) and linearising up to first order in the perturbations leads to

    ∂_t δξ(y,t)=    ·δξ(y,t)+q·δγ̇(y,t)
        +∂_y^2 (𝒟·δξ(y,t)) ,


where =.∂_ξ Q|_ξ_s,γ̇ and q=.∂_γ̇ Q|_ξ_s,γ̇. To eliminate the perturbations in the shear rate we use the force balance equation (<ref>), which in its linearised form is given by

    0=p·δξ(y,t)+η δγ̇(y,t) .


Here p=(0,1) is a projection vector. Combining Eqs. (<ref>) and (<ref>) then yields

    ∂_t δξ(y,t)=(+∂_y^2 𝒟) δξ(y,t) ,


where  is the Jacobian matrix given by

    =-1η q·p .


Explicit expressions for the elements of  are given in Appendix <ref>. We now assume the heterogeneous perturbations to be of exponential form such that

    δξ(y,t)=C e^λ t+iky ,


where C=(C_1,C_2)^T, C_1,C_2∈ℝ, is an arbitrary constant vector, k is the wave number and λ=μ+iω is a complex number with μ,ω∈ℝ. The homogeneous steady state ξ_s is said to be stable, if μ<0, i.e., the perturbation dies out with time. Otherwise it is unstable. Inserting Eq. (<ref>) into Eq. (<ref>) leads to the eigenvalue problem

    det[-k^2 𝒟-λ ]=0 .


From Eq. (<ref>) it is seen that the complex eigenvalues λ depends on the wave number k and also on the imposed shear rate γ̇ (via ). To see the impact of k we plot in Fig. <ref> the largest real part μ_max=μ_max(k) (out of the two real parts of the two eigenvalues) at a fixed structural relaxation time τ_n and a fixed Newtonian viscosity η for three imposed shear rates γ̇. In all cases, μ_max(k) reaches its largest value at k=0. For k>0 the real part is a monotonic decaying function. This behaviour is a consequence of the asymmetric structure of Eq. (<ref>) regarding the diffusion terms. As a result, the eigenvalue equation (<ref>) contains only quadratic powers of k. The behaviour of μ_max(k) indicates that the most unstable mode corresponds to k=0. In the following we therefore restrict our attention to this mode alone. Physically, this restriction implies to consider the system as homogeneous, in other words, 𝒟=0.







The various stability scenarios occurring in this case can be classified according to table <ref>. Investigating the complex eigenvalues for a range of (average) shear rates γ̇ and relaxation times τ_n we obtain a stability diagram, which is presented in Fig. <ref>.

Below a critical point (CP) <cit.> the system shows essentially the same behaviour as the Johnson-Segalman model: There are two stable regions (stable fix point and stable focus), corresponding to homogeneous steady states, and one unstable region (unstable saddle). Within the Johnson-Segalman model the occurrence of the unstable saddle is associated with the regime, where the slope of the constitutive curve _s(γ̇) is negative, respectively where shear banding can occur <cit.>. Above the critical point the present system displays additionally an unstable fix point and an unstable focus. The latter is separated from the stable focus by a Hopf bifurcation (corresponding to μ_1,2=0 and ω_1>0, ω_2<0). The Hopf bifurcation is absent in the Johnson-Segalman model and may thus be considered as a characteristic feature of the present model.



  §.§.§ Numerical results
 

Our stability analysis at k=0 clearly does not provide complete information on the system. To explore the full spatio-temporal dynamics at 𝒟0 we solved Eqs. (<ref>) and (<ref>) numerically by using a Crank-Nicolson algorithm <cit.>. We sample over N_t=30,000 temporal and N_y=150 spatial discrete steps with a corresponding step size dt=0.005 and dy=0.007, yielding a total integration time of T=N_t· dt=150, and a gap size of L=1.0. In analogy to Ref. <cit.>, we fixed the two power law exponents at α=1.2 and β=1.5, and the diffusion constant at 𝒟=0.0016. The equilibrium parameters, namely the zero-shear length and relaxation time, are set to n_0=τ_0=1.0, respectively. To ensure numerical stability we use the Courant-Friedrichs-Lewy criterion for forward time stepping given through √(2 𝒟 dt)<dy <cit.>. To realize unique initial conditions for all numerical calculations, we prepared the stress to be homogeneous in space and the micellar length to be in a shear banded state given by n(y,t=0)=N·[1+Δ·cos(y π)]. In order to check the sensitivity of the system with respect to the initial conditions, we varied the constants N and Δ. The results presented here have been obtained with N=0.5 and Δ=0.5 corresponding to n(t=0)=n_0/2. At the confining walls we employ Dirichlet boundary conditions through .∂_y σ(y,t)|_y=0=0 and .∂_y σ(y,t)|_y=L=0.



In our numerical calculations we performed “continuous shear ramps”. That is, starting from a given configuration of the dynamical variables, we solved the equations at constant γ̇ up to the steady state limit, which is typically achieved at T=150. To proceed towards the next (larger or smaller) value of the shear rate, we took the steady state configuration from the previous one as an initial configuration. By this strategy we are able to avoid stress “overshoots” in the start up regime of each shear step, leading to more stable numerical results. Further, this procedure is more realistic from an experimental point of view.

In Fig. <ref> we present the state diagram for a `down' ramp, that is starting from high shear rates and decreasing the shear rate step by step. It is seen that the system displays a variety of dynamical states, which differ from each other by the spatio-temporal behaviour of the viscoelastic stress σ(y,t) and the other dynamical variables, that is, n(y,t) and γ̇(y,t). We can distinguish six types of patterns (illustrated in Fig. <ref>): spatially homogeneous (h), (static) shear banding (sb), travelling pulses (tp), spatio-temporal oscillatory (o), mixture of spatio-temporal oscillatory and chaotic (oc) and purely chaotic in space and time (c). Some of the results have already been reported in Ref. <cit.> and in fact, our results are consistent with this earlier study. Here we have substantially extended the range of relaxation times τ_n and imposed shear rates γ̇. In order to verify our classification made for the dynamical states we also determined the largest Lyapunov exponent.

Inspecting the numerical results in Fig. <ref> for small relaxation times τ_n≤0.112, we see that the system behaves as predicted by linear stability analysis (and consistent with the Johnson-Segalman model). The regimes at low and high shear rates corresponding to homogeneous steady states are separated by a shear banded region associated with the unstable saddle. By increasing the relaxation time towards values above the critical point (τ_n>0.112), where the Hopf bifurcation arises, we enter a region rich of dynamical states, which have already been observed by Fielding and Olmsted <cit.>. In the subsequent discussion we especially focus on the purely oscillatory and chaotic states (as identified by the Lyapunov exponents), respectively.



Within the oscillatory state the system consists of a well defined high and low shear rate band, which change position periodically in time. Similar states have been observed experimentally <cit.>. In the present model the resulting oscillations of γ̇(y,t) are associated to oscillations of the micellar length n(y,t). The time evolution of these quantities over one period as plotted in Fig. <ref> (a) can be subdivided into four stages:


	
  (i) At very low local shear rate, the micellar length is maximal. From a physical point of view this corresponds to a situation, where the system consists of isotropically distributed and entangled micelles.
	
  (ii) The local shear rate starts to increase, while the micellar length decreases: This can be interpreted such that the micellar network starts to break down (due to increasing shear forces), while the micelles align in the shear flow.
	
  (iii) The local shear rate reaches its maximum, whereas the micellar length is found to be minimal: In this case the shear forces are so large that scission processes are dominating, yielding particularly small values of the micellar length.
	
  (iv) At the end of the cycle the local shear rate is still quite high but so is the degree of alignment of the micelles. Therefore, end-end-recombination processes are more likely to happen, and the micelles start to grow again.


To quantify the oscillatory dynamics inside the unstable regime we extract the fundamental frequency f^∗ from a time series of the total stress (t) [see Eq. (<ref>)] by creating a power spectrum. We then compare f^∗ to the frequency f=ω/(2π) obtained from linear stability analysis by introducing the ratio θ=f/f^∗. Numerical results for θ are shown in Fig. <ref> (b).



Inspecting the frequency ratio θ for the oscillatory and chaotic states, we find that for low and moderate relaxation times, i.e. τ_n=0.140 and τ_n=0.160, no clear trend is visible. For high relaxation times, here τ_n=0.180, the system reveals only oscillatory states. Generally, we find the linear stability results underestimates the frequency for the total investigated shear rate range. Especially at the lower bound of the unstable focus the numerically determined frequency is about a factor 2.5 higher than the estimated one. We recall, however, that the linear stability analysis is restricted to homogeneous states, therefore, deviations are expected. Interestingly, close to the Hopf bifurcation the numerical and analytical results match quite well. These results already give hints where it is possible to manipulate the dynamics of the system using time-delayed feedback control.

So far we have been focusing on the `down' ramp. The corresponding state diagram for an `up' ramp (i.e. increasing γ̇ from low values) is presented in Fig. <ref> (a). Comparing the two state diagrams in Figs. <ref> and <ref> (a) at low relaxation times, i.e. below the critical point (where only shear banded states exist) we conclude that the system reveals the same hysteric behaviour, which was already described for the Johnson-Segalman model elsewhere <cit.>. However, comparing the state diagrams above the critical point, we find for the up ramp the oscillatory region to be extended towards higher shear rates. To clarify this point, we plot in Fig. <ref> (b) numerical steady state results for the total stress for the two ramp protocols. While for the `down' ramp the numerical data leave the analytical steady state solution [see Eq. (<ref>)] exactly at the Hopf bifurcation, the corresponding data for the `up' protocol remains below the analytical one even after crossing the predicted Hopf bifurcation. This indicates that the system reveals a second instability, which is not covered by our linear stability analysis.



 §.§ Stress controlled protocol


For the linear stability analysis of the stress controlled protocol we follow the same procedure outlined in the previous section. As argued in Sec. <ref>, the stress controlled protocol leads to spatially homogeneous dynamics. We therefore set 𝒟=0 and define the vector ξ(t)=[n(t),σ(t)]^T containing the dynamical variables of the system. Within this vector notation Eqs. (<ref>) and (<ref>) can be written as

    ∂_t ξ(t)=Q(ξ,^fix) ,


where Q(ξ,^fix) summarizes the right hand sides of the two dynamical equations. We recall here that, unlike in the shear rate controlled protocol, the total stress is set to a constant. Further, we introduce the vector ξ_s=[n_s,σ_s]^T containing the steady state solutions of the dynamical system [see Eq. (<ref>)]. To determine the stability of these steady states, we add small time-dependent perturbations δξ(t) to ξ_s

    ξ(t)=ξ_s+δξ(t)=ξ_s+C e^λ t ,


with C=(C_1,C_2)^T being a constant vector with C_1,C_2∈ℝ and λ=μ+iω, μ,ω∈ℝ, being the complex eigenvalue. Inserting this ansatz into Eq. (<ref>) and linearising up to first order in the perturbations leads to the eigenvalue problem

    det[-λ ]=0 .


Here =.∂_ξ Q|_ξ_s,^fix is the Jacobian matrix of the dynamical system, which is explicitly derived in Appendix <ref>. When solving the eigenvalue equation (<ref>) we find a quite simple stability diagram, see Fig. <ref>, (compared to the shear rate controlled protocol). It reveals a stable fix point (homogeneous steady state), an unstable focus (homogeneous oscillatory state) delimited by a Hopf bifurcation, and a stable focus (homogeneous steady state). The transition between the region of stable fix points and unstable focus is characterised by hysteretic behaviour (bistability). This reduced complexity compared to Fig. <ref> is a consequence of the fact that within the stress controlled protocol, there are only time-dependent states, such that the dimensionality of the system is reduced to d=2. This excludes to observe any chaotic states, since chaos requires dimension d≥3 <cit.>.





§ FEEDBACK CONTROL


In this section we discuss the application of time-delayed feedback control (TDFC) to our model system. The primary goal thereby is to suppress oscillatory and chaotic states in favour of steady states. This is particularly motivated by an experimental study, where TDFC was used to control chaotic states appearing in Taylor-Couette flow of a complex fluid <cit.>. Generally, the dynamics of a shear-driven system can be manipulated by feeding back a measurable quantity such as the total stress (t) or the shear rate γ̇(t) through a (closed) feedback loop. The control in our system is realized through a Pyragas scheme <cit.>, where the equations of motion of the original system are supplemented by a control term of the form a(t)-a(t-τ). Here, a is a measurable quantity, which can be either the total stress or the shear rate, at the present time t and an earlier time t-τ with τ being the delay time. 

In the following we explore the impact of TDFC on both the shear rate and the stress controlled protocol. Thereby, we restrict ourselves to states close to the Hopf bifurcation, where the model system exhibits a hysteretic behaviour see, e.g.,  Fig. <ref> for the shear rate controlled protocol. In this situation, the relation between (average) shear rate and stress is not unique. The TDFC proposed here aims at selecting one of the branches, particularly, the stationary branch of the two solutions.



 §.§ Shear rate controlled protocol
 

The measurable quantity within the shear rate controlled protocol is the total stress (t) defined in Eq. (<ref>). Using this definition a Pyragas-like feedback term involving  has the form

    (t)-(t-τ)     =σ(t)+η γ̇(t)-σ(t-τ)-η γ̇(t-τ)
        =σ(t)-σ(t-τ) .


Here, we assumed that the (mean) shear rate γ̇(t)=γ̇(t-τ)=γ̇=const. for all t. This seems somewhat non-physical due to the fact that, under normal circumstances, a control of the total stress should affect the shear rate (and vice versa). However, we demonstrated in Sec. <ref> by numerically solving our (uncontrolled) model system, there are state points where the relationship between the total stress and the shear rate is not unique, rather one observes a hysteresis. In the following we will demonstrate that in such a situation, an application of the Pyragas feedback scheme leads to a “stress selection” in the system, that is for certain sets of control parameters K and τ, the system selects the steady state solution _s(γ̇) given in Eq. (<ref>).

Equation (<ref>) further implies that a feedback of the total stress  is equivalent to a (global) feedback of the spatially averaged viscoelastic stress σ. In other words, the TDFC affects only the dynamics of the viscoelastic stress and not that of the micellar length n(y,t), yielding a non-diagonal control scheme. However, linear stability analysis suggests that the global TDFC cannot stabilize spatially inhomogeneous dynamics (see Appendix <ref>), as confirmed by numerical test calculations.

Therefore, we focus here on a local version of the control scheme, i.e. we use the local difference σ(y,t)-σ(y,t-τ) for all following calculations in Sec. <ref>. We are aware that this type of control scheme is difficult to realize in an experimental set up, but from a theoretical point of view it is the most suitable way to introduce TDFC to the shear rate controlled protocol. However, we note already here that a global feedback scheme is successful within the stress controlled protocol (see Sec. <ref>).

In our vectorial notation introduced in Eq. (<ref>) for the uncontrolled system, the dynamical equation for the system under (local) TDFC is given by

    ∂_t ξ(y,t)=     Q(ξ,γ̇)+∂_y^2 (𝒟 ξ(y,t))
        -K  [ξ(y,t)-ξ(y,t-τ)] ,


where 𝒟 is the diffusion matrix defined in Eq. (<ref>) and  is a non-diagonal 2×2 matrix defined by

    =[ 0 0; 0 1 ] .


This form ensures that the control is only applied to σ(y,t). We note that a similar scheme has been proposed by Kyrychko and co-workers for the one-species control in the Gray-Scott model <cit.>, which belongs to the class of reaction-diffusion systems.

We start again by deriving analytic expressions for the controlled eigenvalues ν=μ+i ω in the framework of linear stability analysis, following the procedure outlined for the uncontrolled case in Sec. <ref>. Assuming heterogeneous perturbations δξ(y,t) and δγ̇(y,t), substituting these perturbations into Eq. (<ref>) and linearising up to first order leads to

    ∂_t δξ(y,t)=    (+∂_y^2 D-K ) δξ(y,t)
        +K  δξ(y,t-τ) .


The eigenvalue problem then takes the form (with δξ(y,t)=C e^ν t+iky)

    _11-ν   _12
    _21   _22-k^2 𝒟-K (1-e^-ντ)-ν=0 ,


where _ij (i,j=1,2) are the components of the Jacobian matrix given in Eq. (<ref>). In principle, we have to distinguish here again between two cases, namely k=0 and k>0. In analogy to the uncontrolled case it can be shown that the most unstable mode occurs at k=0 (see also Fig. <ref> (a)). Thus we consider here only the latter case, yielding

    ν^2     +ν {-K (e^-ντ-1)-A}
        +{_11 K (e^-ντ-1)+B}=0 ,


with A=_11+_22 and B=_11 _22-_12 _21.

To simplify the analysis we consider a state close to the Hopf bifurcation, where the real μ part of the eigenvalue ν is approximately zero and thus, ν=i ω is purely imaginary. Substituting this assumption into Eq. (<ref>) leads to

    -ω^2     +iω [-K (e^-iωτ-1)-A]
        +[_11 K (e^-iωτ-1)+B]=0 .


Separating Eq. (<ref>) into real and imaginary part yields

    ω^2-K[_11cos(ωτ)-ωsin(ωτ)-_11]-B     =     0 , 
    ω A-K[-ωcos(ωτ)-_11sin(ωτ)+ω]     =     0 .


Solving Eqs. (<ref>) and (<ref>) yields a pair of solutions (τ,K), which depend on the imaginary part ω of the eigenvalue ν. Specifically, the delay time τ=τ(ω) is given by


    τ(ω)=1ω [arctan(2ω (A _11-B+ω^2) (B _11+(A-_11) ω^2))(B _11)^2-C_1 ω^2+C_2 ω^4-ω^6)+2jπ] ,  j∈ℤ



with C_1=B^2+(A _11)^2+2B _11 (_11-2A) and C_2=A^2-2B-4A _11+_11^2. These two constants depend on the system parameters α, β, τ_n, η and γ̇. Further, the argument of the arctangent in Eq. (<ref>) is by definition limited to the interval [-π/2,+π/2]. The control strength K=K(ω) follows as

    K(ω)=(A ω)^2+(ω^2-B)^22 (B _11+ω^2 (A-_11)) .




These results are similar to the ones derived in Ref. <cit.> for the Gray-Scott model. In Fig. <ref> (b) we plot the stability borders corresponding to the expressions for τ and K in Eqs. (<ref>) and (<ref>) for several shear rates at a fixed micellar relaxation time τ_n=0.180. With increasing γ̇ we observe the minima of the stability curves to be shifted towards smaller delay times. Further, the width of the stability domains is increasing. Interestingly, the results in Fig. <ref> (b) are similar to those obtained for a diagonal control scheme (see Refs. <cit.>), where a (local) Pyragas term is applied to both, the viscoelastic stress and the micellar length.

As indicated at the beginning of Sec. <ref>, we apply the feedback control to states close to the (analytical) Hopf bifurcation (red line in Fig. <ref> (a)), where the system reveals a hysteretic behaviour as shown in Fig. <ref> (a). Specifically, for fixed γ̇, we find the numerical solution for the total stress to have two solution branches, one corresponding to the (constant) steady state solution _s(γ̇) determined in Eq. (<ref>), and the other one to a solution where the total stress fluctuates in time (the solution plotted in Fig. <ref> (a) corresponds to the time average).

In order to see the impact of the (non-diagonal) control scheme on the spatio-temporal dynamics, we record the total stress (t) (exemplary shown in Fig. <ref> (b) for a fixed control strength and two delay times). In the uncontrolled case (t<75) the system is in an oscillatory state (I). When the control term is switched on (t≥75) we find two scenarios: The initial state (I) can be either fully stabilized (II), i.e. the total stress reaches a constant value corresponding to the steady state solution _s(γ̇) in the long time limit, or the control scheme is not working (III) and the oscillating character of the initial state is preserved, but the amplitude of the oscillations is damped. 

We stress that, in case of reaching the steady state after applying the feedback control, the solution is not chosen arbitrarily, but is already one of the possible solutions (namely the stationary branch) given by the uncontrolled system (due to the hysteretic behaviour shown in Fig. <ref> (a)). This selection of an existing solution underlines the non-invasive character of the Pyragas scheme.

From the time series (t) we calculate the largest Lyapunov exponent Λ_max by using the Wolf-algorithm <cit.>. The results for Λ_max are presented in Fig. <ref> (c). In case of successful control we find Λ_max<0 (corresponding to region (II) in Fig. <ref> (b)). If the control scheme is not successful, the Lyapunov exponent is positive (Λ_max>0). To visualize the impact of the feedback control on the dynamics of the system we present in Fig. <ref> the spatio-temporal evolution of σ(y,t), which correspond to the time series of the total stress presented in Fig. <ref> (b). The situation shown in Fig. <ref> (a) refers to a successful application of the non-diagonal TDFC scheme, i.e. a spatio-temporal inhomogeneous state is stabilized to a homogeneous steady state, while in Fig. <ref> (b) the control scheme does not work in the sense that the inhomogeneous character of the initial state is preserved.







 §.§ Stress controlled protocol


Similarly to the shear rate controlled protocol, we can also construct a Pyragas-like feedback term within the stress controlled protocol. The corresponding Pyragas-like term to be added in the equation of motion has now the form

    γ̇(t)-γ̇(t-τ)     =1η [^fix-σ(t)-^fix+σ(t-τ)]
        =-1η [σ(t)-σ(t-τ)] ,


where we used the definition of the shear rate in the stress controlled protocol given in Eq. (<ref>). Henceforth, we focus on the non-diffusive limit with 𝒟=0, following the arguments presented in Sec. <ref>. We again restrict ourselves to states close to the Hopf bifurcation, where the system reveals temporal oscillations. Within this regime we find the shear rate oscillating around its steady state values (see Fig. <ref> (a)). We will show that the feedback-controlled system (analogous to the situation considered in the previous section) selects this steady state solution for certain sets of control parameters K and τ.

The form of Eq. (<ref>) already suggests that the feedback control should again be applied as a non-diagonal scheme, i.e. the TDFC is only acting on the viscoelastic stress. Using the vector notation we introduced in Eq. (<ref>) for the stress controlled protocol, we can formulate a (generalized) feedback equation such that

    ∂_t ξ(t)=Q(ξ,^fix)-K^∗  [ξ(t)-ξ(t-τ)] ,


where  is the matrix defined in Eq. (<ref>) ensuring that the TDFC is only acting on the viscoelastic stress. Taking Eq. (<ref>) into account yields the specific feedback equation

    ∂_t ξ(t)=Q(ξ,^fix)+Kη  [ξ(t)-ξ(t-τ)] .


Comparing Eqs. (<ref>) and (<ref>) we identify K^∗=-K/η as generalized feedback strength. Due to reasons of comparability with the previous section, we focus in the following on the generalized delay differential equation (<ref>). To determine the stability domains of Eq. (<ref>) we follow the same procedure as outlined in Sec. <ref> for the shear rate controlled protocol. Since we consider the non-diffusive limit, the perturbations δξ(t)=C e^ν t with C=(C_1,C_2)^T are only time-dependent. Further, the Jacobian matrix now depends on the fixed total stress ^fix (see also Appendix <ref>).

Interestingly, we find the shape of the neutral stability curves to be unchanged in comparison to the ones of the shear rate controlled protocol depicted in Fig. <ref> (b). The reason for this behaviour can be found in the coupling between the shear rate γ̇ and the total stress  given in Eq. (<ref>), i.e. the shear rates chosen in Fig. <ref> (b) can be converted to the corresponding stresses. 

More interesting is the impact of the TDFC scheme on the dynamics of stress controlled protocol, which is depicted in Fig. <ref> (b). Measuring the shear rate as function of time we find the uncontrolled system to reveal oscillations, which are damped out if the control is activated with (τ,K^∗)=(0.2,3.0) for times t>75 (red lines in Fig. <ref> (b)). The resulting constant value is the same as the one determined from the steady state solution (see Fig. <ref> (a)). This again emphasizes the non-invasive character of the Pyragas scheme, i.e. the system takes an already existing solution when applying TDFC. On the other hand for the parameter set (τ,K^∗)=(0.4,3.0) the oscillations are preserved, but the amplitude is changed (grey lines in Fig. <ref> (b)).





§ CONCLUSIONS


In this publication we have applied time-delayed feedback control (specifically, the Pyragas scheme) to manipulate the complex dynamics of a solution of wormlike micelles based on a non-local rheological model proposed by Fielding and Olmsted <cit.>. As a background we have performed a linear stability analysis (at k=0) for uncontrolled systems at constant (average) shear rate, as well as systems at constant total stress. Further, we solved the dynamical system numerically to get an insight into the full spatio-temporal behaviour for a broad range of parameters. In doing this, we considered both an “up” ramp, i.e. increasing the shear from low values, and a “down” ramp, that is, decreasing the shear rate from high values. Comparing the corresponding numerical results with the analytical predictions from linear stability analysis (at k=0), we find a good agreement of the data for the “down” ramp. For the “up” ramp, the system displays hysteric behaviour not captured by the linear stability analysis. This hysteresis effect has already been observed in the Johnson-Segalman model <cit.> within the bistable region related to shear banding. In the model investigated here, the hysteresis occurs near the Hopf bifurcation located at the high shear rate branch of the constitutive curve.

We have then investigated the changes of the dynamics upon application of a global and local Pyragas control scheme <cit.>, where the control variable is associated with a measurable quantity of the corresponding rheological protocol. In case of the shear rate controlled protocol, we have shown that a local control scheme for the stress can stabilize the spatio-temporal oscillatory dynamics close to the Hopf bifurcation, where the system exhibits hysteresis. Specifically, by applying feedback control of suitable strength and delay time, the system selects the steady state solution branch. In contrast, a global scheme acting only on the spatially averaged stress does not work. However, a global scheme is indeed successful within the stress controlled protocol, where the dynamics is only time dependent. Again we only focused on states close to the analytically determined Hopf bifurcation, where the system reveals temporal oscillations. Here, the steady state solution is again part of the full solution space of the system. In analogy to the shear rate controlled protocol, we have shown that the feedback control can select the steady state solution branch. This again underlines the non-invasive character of the Pyragas scheme.

We are aware that the control scheme presented here is not generally applicable to all shear-induced states. The reason is that, outside the hysteretic regions considered here, the relation between  and γ̇ is unique, that is, one cannot control one of these variables without affecting the other one. Therefore, it is reasonable to think of alternative ways of implementing feedback control to sheared systems. On possibility might be a external (mechanical) control, where the stress  or the shear rate γ̇ is serving as a measurable quantity, which is fed back to the system using a transducer manipulating the (externally) imposed variable. To implement such kind of feedback the (full) Navier-Stokes equation is required to get access to the velocity profile and therewith define a coupling to a transducer. Another possible scenario would be the implementation of an internal control, acting on the structural variables: Hereby, the externally imposed rheological variable (γ̇ or ^fix) is kept constant and the structural properties of the system are changed by optical, electric or magnetic fields, that is, the conformation or the stiffness of aggregates and molecules will be tuned through interactions with an external field. In this case, a further dynamical equation for the external field has to be defined. Clearly, it would be very interesting to test the implications of such feedback control (or related schemes) experimentally. Especially the internal control might be relevant in view of the increasing importance of stimuli responsive materials <cit.>.

We also mention recent progress in accessing the spatio-temporal structure of micellar solutions under shear using appropriate scattering techniques (Rheo-SANS) <cit.>. Given these developments, an application of TDFC to avoid or stabilize certain shear-induced behaviours seems very promising. Further, it seems to be interesting to systematically vary the stress-diffusion constant 𝒟 and the power law exponents α and β to determine the impact on the dynamics of the system and to compare them to experimental data.


We gratefully acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG) through project B2 within the Collaborative Research Center 910. The authors thank M. Gradzielski for fruitful discussions and the insight from the experimental point of view.






§ JACOBIAN MATRIX FOR THE SHEAR RATE CONTROLLED PROTOCOL
 

We start from the compact notation for the dynamical system given in Eq. (<ref>) in Sec. <ref>

    ∂_t ξ(y,t)=Q(ξ,γ̇)+∂_y^2 (𝒟 ξ(y,t)) ,


where 𝒟 is the diffusion matrix and the vector Q(ξ,γ̇) follows from Eq. (<ref>) and (<ref>)

    Q(ξ,γ̇)     =[ -n(y,t)/τ_n+1/τ_n (n_0/1+(τ_n γ̇(y,t))^β);   -σ(y,t)/τ(n)+γ̇(y,t)/1+(τ(n) γ̇(y,t))^2 ]
        =[ Q_1(n,σ,γ̇); Q_2(n,σ,γ̇) ] .


The Jacobian matrix  for the shear rate controlled protocol is defined by

    =-1η q·p ,


with =.∂_ξ Q|_ξ_s,γ̇ and q=.∂_γ̇ Q|_ξ_s,γ̇. We first consider the matrix . The latter can be written in terms of the components of the vector Q

    =.∂_ξ Q|_ξ_s,γ̇=.[ ∂_n Q_1 ∂_σ Q_1; ∂_n Q_2 ∂_σ Q_2 ]|_ξ_s,γ̇ .


Inserting the expressions for Q_1(n,σ,γ̇) and Q_2(n,σ,γ̇) given in Eq. (<ref>) and evaluating the derivatives at the steady state values one obtains

    .∂_n Q_1(n,σ,γ̇)|_ξ_s,γ̇    =     -1τ_n
    .∂_σ Q_1(n,σ,γ̇)|_ξ_s,γ̇    =     0 
    .∂_n Q_2(n,σ,γ̇)|_ξ_s,γ̇    =    α σ_sn_s^α+1-2α γ̇^3 n_s^2α-1[1+(γ̇ n_s^α)^2]^2
    .∂_σ Q_2(n,σ,γ̇)|_ξ_s,γ̇    =     -1n_s^α .


Second, we consider the vector q,

    q=∂_γ̇ .Q|_ξ_s,γ̇=.[ ∂_γ̇ Q_1; ∂_γ̇ Q_2 ]|_ξ_s,γ̇ .


Using Eq. (<ref>) we obtain

    .∂_γ̇ Q_1(n,σ,γ̇)|_ξ_s,γ̇    =     -β (τ_n γ̇)^βτ_n γ̇ [1+(τ_n γ̇)^β]^2 , 
    .∂_γ̇ Q_2(n,σ,γ̇)|_ξ_s,γ̇    =    1-(γ̇ n_s^α)^2[1+(γ̇ n_s^α)^2]^2 .


Inserting the elements of  and q into Eq. (<ref>) for the matrix , one finds the following elements _ij (i,j=1,2)

    _11    =     -1τ_n
    _12    =    β (τ_n γ̇)^βη τ_n γ̇ [1+(τ_n γ̇)^β]^2
    _21    =    α σ_sn_s^α+1-2α γ̇^3 n_s^2α-1[1+(γ̇ n_s^α)^2]^2
    _22    =     -1n_s^α-1-(γ̇ n_s^α)^2η [1+(γ̇ n_s^α)^2]^2 .


For convenience we set n_0=τ_0=1.0 for all calculations.



§ JACOBIAN MATRIX FOR THE STRESS CONTROLLED PROTOCOL
 

In the stress controlled protocol the shear rate is replaced by γ̇(y,t)=1/η [^fix-σ(y,t)], with ^fix being a constant. In vector notation the dynamical system can be written

    ∂_t ξ(y,t)=Q(ξ,^fix) ,


where the vector Q(ξ,^fix) is given by

    Q(ξ,^fix)=    [         -n(y,t)/τ_n+1/τ_n (n_0/1+(τ_n/η [^fix-σ(y,t)])^β); -σ(y,t)/τ(n)+1/η [^fix-σ(y,t)]/1+[τ(n)/η [^fix-σ(y,t)]]^2 ]
        =[ Q_1(n,σ,^fix); Q_2(n,σ,^fix) ] .


The Jacobian matrix  evaluated at the steady state values can be written as

    =.∂_ξ Q|_ξ_s,^fix=.[ ∂_n Q_1 ∂_σ Q_1; ∂_n Q_2 ∂_σ Q_2 ]|_ξ_s,^fix .


Explicitly, the components _ij (i,j=1,2) of the Jacobian matrix are given by

    _11    =     -1τ_n
    _12    =    β [τ_n/η (^fix-σ_s)]^β-1η [1+(τ_n/η (^fix-σ_s))^β]^2
    _21    =    αn_s^α+1 [σ_s+2η (n_s^α (^fix-σ_s))^3[η^2+(n_s^α (^fix-σ_s))^2]^2] 
    _22    =     -1n_s^α+η [(n_s^α (^fix-σ_s))^2-η^2][(n_s^α (^fix-σ_s))^2+η^2]^2 .


For convenience we set n_0=τ_0=1.0 for all calculations.



§ LINEAR STABILITY OF THE GLOBAL FEEDBACK SCHEME IN THE SHEAR RATE CONTROLLED PROTOCOL
 

We start our derivation of the linear stability analysis of the global feedback control scheme in the shear rate controlled protocol with the dynamical equations in vector notation (analogous to Eq. (<ref>)), that is

    ∂_t ξ(y,t)=     Q(ξ,γ̇)+∂_y^2 (𝒟 ξ(y,t))
        -K  [ξ(t)-ξ(t-τ)] ,


where 𝒟 is the diffusion matrix defined in Eq. (<ref>) and  is a non-diagonal 2×2 matrix defined in Eq. (<ref>) ensuring that the feedback control is only acting on the viscoelastic stress σ. Further, ξ represents the spatially averaged version of the dynamical vector ξ(y,t)=[n(y,t),σ(y,t)]^T. Decomposing the dynamical variable ξ(y,t) as well as the spatially averaged variable ξ(t) into a steady state plus small perturbations, substituting this back to Eq. (<ref>) and linearising in the perturbations up to first order gives rise to the following equation

    ∂_t δξ(y,t)=    (+∂_y^2 𝒟) δξ(y,t)
        -K  [δξ(t)-δξ(t-τ)] ,


where  is the Jacobian matrix derived in Appendix <ref>. For the perturbations δξ(y,t) we choose the following ansatz

    δξ(y,t)=C e^ν t+iky ,


with C=(C_1,C_2)^T being an arbitrary constant vector, ν the complex eigenvalues and k the wave number. Averaging the perturbations in Eq. (<ref>) over space yields

    δξ(t)     =1L ∫_0^Ldy δξ(y,t)=C e^ν tL ∫_0^Ldy e^iky
        =C e^ν t δ_k,0 .


Hereby δ_k,0 is the Kronecker-delta. Substituting Eqs. (<ref>) and (<ref>) back in Eq. (<ref>) leads to

    C ν e^ν t+iky=    (+∂_y^2 𝒟) C e^ν t+iky
        -K  C e^ν t δ_k,0 [1-e^-ντ] .


This equation can be transferred to

    (+∂_y^2 𝒟-ν ) e^iky-K   δ_k,0 [1-e^-ντ]=0 .


We have to distinguish two cases: For k=0 Eq. (<ref>) is reduced to the eigenvalue problem we described for the local TDFC in Sec. <ref>. The case of k0 yields a reduction to the eigenvalue problem of the system without feedback control as discussed in Sec. <ref>. This implies that only the long wave limit can be stabilized, i.e. general validity for all wave numbers k is not guaranteed, which is a necessary condition for linear stability analysis. We also confirmed these results by numerical test calculations.


apsrev

