A theory for the flow of chemically-responsive polymer solutions: equilibrium and shear-induced phase separation

Chemically-responsive polymers are macromolecules that respond to local variations of the chemical composition of the solution by changing their conformation, with notable examples including polyelectrolytes, proteins and DNA. The polymer conformation changes can occur in response to changes to the pH, the ionic strength or to the concentration of a generic solute that interacts with the polymer. These chemical stimuli can lead to drastic variations of the polymer flexibility and even trigger a transition from a coil to a globule polymer conformation. In many situations the spatial distribution of the chemical stimuli can be highly inhomogeneous, which can lead to large spatial variations of polymer conformation and of the rheological properties of the mixture. In this paper, we develop a theory for the flow of a mixture of a solute and chemically-responsive polymers. The approach is valid for generic flows and inhomogeneous distributions of polymers and solutes. To model the polymer conformation changes introduced by the interactions with the solute, we consider the polymers as linear elastic dumbbells whose spring stiffness depends on the solute concentration. We use the Onsager's variational formalism to derive the equations governing the evolution of the variables, which unveils novel couplings between the distribution of dumbbells and that of the solute. Finally, we use a linear stability analysis to show that the governing equations predict an equilibrium phase separation and a distinct shear-induced phase separation whereby a homogeneous distribution of solute and dumbbells spontaneously demix. Similar phase transitions have been observed in previous experiments using stimuli-responsive polymers and may play an important role in living systems.

Chemically-responsive polymers are macromolecules that respond to local variations of the chemical composition of the solution by changing their conformation, with notable examples including polyelectrolytes, proteins and DNA. The polymer conformation changes can occur in response to changes to the pH, the ionic strength or to the concentration of a generic solute that interacts with the polymer. These chemical stimuli can lead to drastic variations of the polymer flexibility and even trigger a transition from a coil to a globule polymer conformation. In many situations, the spatial distribution of the chemical stimuli can be highly inhomogeneous, which can lead to large spatial variations of polymer conformation and of the rheological properties of the mixture. In this paper, we develop a theory for the flow of a mixture of a solute and chemically-responsive polymers. The approach is valid for generic flows and inhomogeneous distributions of polymers and solutes. To model the polymer conformation changes introduced by the interactions with the solute, we consider the polymers as linear elastic dumbbells whose spring stiffness depends on the solute concentration. We use the Onsager's variational formalism to derive the equations governing the evolution of the variables, which unveils novel couplings between the distribution of dumbbells and that of the solute. Finally, we use a linear stability analysis to show that the governing equations predict an equilibrium phase separation and a dis-

I. INTRODUCTION
Chemically-responsive polymers constitute a class of macromolecules that can change their conformation in response to chemical stimuli such as changing the chemical composition of a suspension 1 . These types of polymers is commonly encountered in several biomedical and industrial applications 2,3 . For instance, polyelectrolytes can change their conformation in response to changes of the pH or of the ionic strength of the solution. Change of local salt concentration screens the charges along the backbone of the polyelectrolytes thus making the polymer chain more flexible and promoting a transition from an extended semiflexible chain configuration to a flexible and coiled one [4][5][6][7] . The gradual change in the polymer conformation has been observed in many experiments that measured the persistence length of the chain at different salt concentrations using optical and magnetic tweezers [8][9][10] . A more extreme type of chemically-responsive behavior is displayed by linear poly(N-isopropylacrylamide)(PNIPAM) polymers. Experiments and simulations measuring the extension of PNIPAM chains at different salt and solute concentrations [11][12][13] found a transition between a coiled conformation and a globular one. Such coil to globule transition is often accompanied by a phase separation of the polymeric phase.
Not surprisingly, the conformation of many bio-polymers that are crucial for the correct functioning of living organisms is sensitive to small changes to the chemical composition of its environment. This feature can be exploited by eukaryotic cells and bacteria to actively control the function of these important macromolecules 14,15 . For instance, DNA and RNA can change their conformation depending on the local concentration of binding proteins 16,17 . Similarly, proteins can change their conformation in response to chemical stimuli such as changes of pH or ionic strength 18 . In some cases, the conformational changes of these macromolecules are associated to phase separation and organization of these biopolymers into membrane-less compartments, which carry out critical tasks for cell survival [19][20][21][22] . Finally, a last class of chemically-responsive polymers are those that weakly bind in the presence of molecules that bridge different groups along the polymer backbone, thus changing its conformation 23,24 .
In general, the conformational changes of chemically-responsive polymers affect the rheological properties of the suspension 25 . For instance, in the case of polyelectrolytes, increasing the salt concentration leads to more compact configurations that are not easily stretched by straining flows. Indeed, experiments in shear flow show that the viscosity decreases as the salt concentration is increased, which is a signature of the smaller dissipation generated by the macromolecules in their compact configuration [26][27][28][29] . Besides changing the viscosity of the solution, conformational changes can also impact its viscoelastic properties because more compact macromolecules relax faster to their equilibrium shape than those that are in an extended conformation 25 . This qualitative picture is confirmed by small-amplitude oscillatory shear and microfluidics experiments showing a change of the relaxation time as a function of the concentration of a solute [30][31][32] .
Standard theories for polymer solutions such as the Flory-Huggins theory 25 capture the interaction between a polymer and a solvent whose concentration is homogeneous and has been used to predict the phase separation, the coil-globule transition and the dependence of the elastic modulus of polymer solutions depending on the solvent quality. However, this theory cannot describe polymers in a solvent interacting with solute molecules since in general solute concentration can be heterogeneous as a consequence of chemical reactions, external gradients, solute or polymer advection or other nonequilibrium processes that are ubiquitous in living systems. Even in the simplest case of a shear flow with a homogeneous polymer concentration, flow instabilities can lead to gradients of concentration transverse to the flow direction [33][34][35][36][37][38][39][40] .
Another relevant question arises in the simulation of inhomogeneous flows of chemicallyresponsive polymers where the solute concentration changes in space leading to inhomogeneous rheological properties. How should one modify the viscoelastic constitutive equations 41 to include the solute-dependent rheological properties? A naive answer would be to simply use spatiallyvariable viscosity and relaxation time. However, this procedure might ignore important couplings between the gradients of solute and polymer density. Furthermore, since the polymer conformation depends on the local chemical composition, it is interesting to investigate if unexpected couplings between solute and polymer conformation could occur.
In this paper, we address these points by proposing a thermodynamically self-consistent model for the isothermal flow of a dilute mixture of solvent, solute and chemically-responsive polymers.
Our model is applicable to general flows with spatially-variable strain rates, not necessarily unidirectional, and inhomogeneous solute and polymer concentrations. To develop a general model of chemically-responsive polymers, we assume that conformation changes driven by the presence of a solute can be modeled as a change of the polymer stiffness. By using the Onsager's variational formalism, we derive systematically the set of equations that govern the flow of a mixture of solute and chemically-responsive polymers. These equations reveal couplings between the gradients of solute, the gradients of polymer density and the gradients of polymer conformation, which are not present in previous viscoelastic models. Finally, we show that the solute-polymer coupling leads to distinct mechanisms of demixing of a homogeneous solution, equilibrium and shear-induced.

II. MODEL FOR POLYMER-SOLUTE INTERACTIONS
To describe chemically-responsive polymer solutions in general flow conditions and heterogeneous concentrations, we make the simple choice of modeling polymers as dumbbells. i.e. two beads connected by a spring, see Figure 1. The vector joining the two spheres of the dumbbell represents the end-to-end vector of the real polymer. The dumbbell model is a good approximation of polymer suspensions in dilute and semidilute regimes where the polymer chains are not entangled with each other 42 . We assume that the stiffness of the spring connecting the two beads of a dumbbell depends on the concentration of solute. By doing so, we provide a simple, yet general, model for the conformation change of polymers introduced by the presence of a solute. Since the extension of the dumbbells represents the end-to-end vector of the polymer, by considering a stiffness that depends on the solute concentration we can model chemically-induced conformational changes. For instance, the coil to globule transition of polymers caused by the interactions with a solute can be modeled through a steep change of the spring stiffness with the concentration of solute, see Figure 1.
We assume that the dumbbells and the solute are suspended in a solvent of viscosity η and that the solution is incompressible and isothermal. We are interested in spatially-extended systems in which the concentration of polymers and of solute can change in space. Therefore, we characterize the distribution of dumbbells inside a volume, Ω, through the function ψ(x, r,t) , where x denotes the position of the center of a dumbbell in space, r denotes its end-to-end vector and t is the time, see Figures 1(a-b). The local number density of dumbbells at a given point in space is given by the integral over the end-to-end distance space: The distribution of solute is simply characterized by its number density c(x,t). We assume that the elastic energy stored by each dumbbell molecule is given by a linear spring: The elastic energy, given by Eq. (2), has an entropic origin that reflects the resistance of the polymer chains to be stretched away from their equilibrium conformation. In the equation above, k(c) represents a solute-dependent stiffness, which models the effects of the solute on the conformation of a polymer molecule. This function could be obtained from experiments measuring the persistence length of polymers at different solute concentrations or could be estimated from scaling laws. In the case of polyelectrolytes, the persistence length of a polymer depends on the local salt concentration and the function k(c) can be estimated directly from theoretical considerations 7 .
The persistence length, l p , of polyelectrolytes is proportional to the Debye-length, which scales as c −1/2 with c the number density of salt. For a flexible worm-like chain, the Kuhn length, b, is twice the persistence length 25 . It follows that the entropic elasticity of a polymer, which is equal to 3k B T /Nb 2 with N the number of monomers, is a function of the local salt concentration. Using the scaling l p ∝ c −1/2 and b = 2 l p , we can estimate the function k(c) for the case of a polyelectrolyte in an ionic suspension as k(c) ∝ c, where c is the local salt concentration. We emphasize that the elastic energy given by Eq. (2) is a simple but rather crude approximation of the true elastic energy stored by a stretched polymer and it is expected to be valid only for small deformations. It is known that the choice of a linear spring leads to the divergence of the extensional viscosity at large extension rates 41 . Nevertheless, in this work, we focus on the effect of a solute-dependent stiff-ness and we leave the extension of Equation (2) to include finite spring extensibility or excluded volume effects 42 to future works. Throughout the rest of the paper, we drop the dependence of the dumbbell distribution, ψ, on x, t and r and the dependence of the solute distribution, c, on x and t for clarity.

III. DERIVATION OF THE GOVERNING EQUATIONS
We begin the derivation of the governing equations by identifying the relevant conservation equations for the dumbbell and solute distribution. The conservation of the number density of each species implies that the distribution of dumbbells must satisfy 43 In Eq. (3),ṙ represents the rate of change of the end-to-end vector of the dumbbells and w d is the velocity of the center of mass of the dumbbells. The conservation of the number density of the where w is the velocity of the solute. Mass conservation for an incompressible solvent implies: where v is the velocity of the solvent. We assume that the boundaries of the volume, δ Ω, are impermeable to the dumbbells, to the solute and to the solvent: Where n is the vector normal to the volume boundaries. Note that the framework developed here can be extended to consider adsorption and desorption of solute and dumbbells from the surfaces.
Finally, the distribution of dumbbells goes to zero ψ → 0 as the end-to-end vector goes to infinity r → ∞: there are no dumbbells that are infinitely stretched.
To derive the balance of momentum and the expressions forṙ, w d and w, we neglect inertial effects and we employ the Onsager's variational formalism 44,45 . This framework allows us to identify the relevant couplings between the fluxes of solute, of dumbbells and of momentum directly from the principle of least dissipation. Briefly, the governing equations are obtained by minimizing the Rayleighian functional, constrained by the incompressibility condition, with respect to the process variables. In our case, the process variables are given by the velocity of the dumbbells, w d , by the velocity of the solute, w, by the velocity of the solvent, v, and by the rate of change of the end-to-end vectorṙ. To constrain the Rayleighian to fulfill the incompressibility condition ∂ ∂ x · v = 0, we include the pressure p as Lagrange multiplier and we define the Lagrangian functional as In Eq. (7),Ḟ and D are the rate of change of the free energy and the energy dissipation of the system, respectively. By minimizing L with respect to the process variables and with respect to the pressure we obtain the governing equations. Recently, Zhou and Doi 46 showed that the Onsager's variational principle can be used to derive the Oldroyd-B constitutive model in the case of a homogeneous suspension of dumbbells. Our approach bears some similarities to the GENERIC and bracket frameworks of nonequilibrium thermodynamics 47,48 , which have been used successfully to derive constitutive equations of complex fluids that couple different physics [49][50][51][52][53][54] . While it is more transparent than the GENERIC formalism, our approach implicitly assumes that the momentum of all the fields are fast variables and relax quickly compared to the evolution of the system 55 .

A. Rate of change of the free energy
We begin by defining the free energy functional of the system, F , as the sum of the contribution of the dumbbells and the solute where k B is the Boltzmann's constant, T is the absolute temperature and U is the elastic energy, defined in Eq. (2). The first two terms represent the solute and dumbbell free energy of mixing and the third term represents the elastic energy stored by the springs. Since the stiffness of the spring, k(c), depends on the local solute concentration, the energy stored by the dumbbells also depends on c. The rate of change of the free energyḞ is obtained by taking the time derivative of Eq. (8), see the Appendix A for details,

B. Energy dissipation
We assume that the energy dissipation functional is comprised of three different contributions: (i) the friction due to the viscosity of the solvent, which is proportional to the solvent shear viscosity η and to the rate of deformation tensor; (ii) the friction due to the relative motion of the solute molecules and the solvent, which is proportional to the friction coefficient of the solute ξ ; and (iii) the friction between the dumbbells and the solvent, which is proportional to the friction coefficient of the dumbbell beads, ξ d . For simplicity, we assume ξ d to be a constant and not to depend on the dumbbell extension 56 or on the conformation of the nearby dumbbells as in the Giesekus model 57 .
This is a simplification since we expect that changes in polymer conformation do impact the dissipation when moving relative to the solvent. Under these assumptions, the dissipation functional is given by: The first term in Eq. (10) represents the dissipation due to the solvent rate of deformation, the second term represents the dissipation due to the relative motion of solute and solvent molecules. The last two terms of Eq. (10) are the dissipation generated by the relative motion between the dumbbells and the solvent and can be derived from a Taylor expansion of the flow field around the dumbbell center of mass, as shown in the Appendix B.

C. Minimization of the Lagrangian
To derive the governing equations, we substitute the rate of change of the free energy, given by Eq. (9), and the dissipation functional, given by Eq. (10), in the definition of the Lagrangian, given by Eq. (7). The minimization of L with respect to the solute velocity, w, yields an expression for the solute flux: The minimization of L with respect to the translational velocity of the dumbbells, w d , yields an expression for the spatial flux of dumbbells: The minimization of L with respect to the rate of change of the end-to-end vector of the dumbbells,ṙ, yields: Finally, minimizing L with respect to the solvent velocity, v, and using integration by parts yields the momentum balance: where we substituted the definition of elastic energy, U, as given by Eq. (2). Besides the classic Fickian diffusion, with diffusion coefficient, k B T /ξ , and the advection by the solvent flow, the right-hand side of Eq. (16) displays one more contribution to the solute flux that is proportional to dk(c) dc . It represents a coupling between the distribution of dumbbells and the distribution of solute, which is introduced by the solute-dependent stiffness of the dumbbells. In the case of constant spring stiffness, the transport of solute reduces to the classic advection-diffusion equation. Since the integral of ψ |r| 2 over r represents the local average extension of the dumbbells, the new coupling term predicts that a gradient of dumbbell extension drives a flux of solute. As a consequence of this coupling, in the case of inhomogeneous flows where the stretching of polymers changes in space, we expect that the solute concentration is also inhomogeneous. The solute will migrate towards regions of higher or lower dumbbell extensions, depending on the sign of dk(c) dc . From a physical standpoint, the migration of solute is driven by the tendency of the solute to reduce the elastic energy stored by highly stretched dumbells. Therefore, if dk(c) dc < 0, the solute migrates towards regions of large dumbbell stretch to make their spring softer and reduce the elastic energy.
Conversely, if dk(c) dc > 0, the solute migrates away from regions of large dumbbell extensions so to decrease the stiffness of the dumbbells in these regions.
By substituting Eqs. (13)- (14) into the Eq. (3) we obtain an equation for transport of dumbbells: The first bracket on the right-hand side of Eq. (17) represents the flux in the end-to-end space, while the second term on the right-hand side represents the spatial flux of dumbbells. The flux of dumbbells in the end-to-end vector space, r, is composed of three terms: the first one is a Brownian relaxation term, the second one is the effect of the restoring force due to the spring and the third is the stretching due to straining solvent flows. This flux is similar to that derived in other dumbbell kinetic models 42 , with the exception that here the spring stiffness depends on the local solute concentration. As a consequence, the relaxation time of the end-to-end vector, which is proportional to ξ d /k(c), can change in space. The spatial flux of dumbbells has three contributions: a Fickian diffusion term with an effective diffusion coefficient given by k B T /2 ξ d , a contribution originating from the coupling with the solute field, and a contribution due to the advection by the solvent flow. The second contribution, proportional to dk(c) dc , is novel and it predicts that a gradient of solute drives a spatial flux of dumbbells. The flux arising from gradients of solute is proportional to ψ|r| 2 , which represents the local extension of the dumbbells. The more stretched are the dumbbells the larger is their spatial flux resulting from a gradient of solute. The direction of the dumbbell flux depends on the sign of dk(c) dc . The physical mechanism driving the migration of dumbbells is similar to that of the solute. In the case of an inhomogeneous distribution of solute, the dumbbells migrate towards regions of higher or lower concentrations of solute to reduce the elastic energy stored by the spring. We emphasize that this mechanism of polymer migration is fundamentally different from the well-known mechanisms driven by curved streamlines 58,59 or by stress gradients 36,60-62 .
By substituting Eqs. (13)- (14) into the momentum balance, given by Eq. (3), we obtain: where rr denotes the dyadic product and I denotes the identity tensor. In Eq. (18), we have defined a modified pressure, P, which is the sum of the hydrodynamic pressure and of the osmotic pressure: where n d is the local number density of dumbbells defined in Eq. (1). The stress generated by the elasticity of the dumbells is proportional to the spring k(c), which can vary in space. Finally, by maximizing the Lagrangian, given by Eq. (7), with respect to the pressure, p, we obtain the By neglecting the dependence of the spring stiffness on the solute concentration the transport of solute decouples from the transport of dumbbells. If we also neglect the spatial variations of the dumbbell concentration, then Eqs. (17) and Eq. (18)

IV. EQUILIBRIUM PHASE SEPARATION OF A HOMOGENEOUS SUSPENSION
The well-established Flory-Huggins theory shows that changing the solvent quality can lead to phase separation 25 . Since changing the solute concentration can be interpreted as a change in solvent quality, it is interesting to ask if the simple model that we propose also predicts an equilibrium phase separation. To do so, we study the stability of a homogeneous suspension of solute and dumbbells at equilibrium using linear stability analysis. We assume small perturbations to the equilibrium concentration of solute, c eq and to the equilibrium dumbbell distribution ψ eq : The equilibrium distribution of dumbbells is given by the Boltzmann distribution: where n d,eq is the number density of dumbbells and k eq denotes the stiffness of the dumbbell spring at equilibrium: it is a shorthand notation for k(c eq ). It is straightforward to show that the integral of ψ eq over the end-to-end vector space gives n d,eq .
By substituting the perturbations, given by Eqs. (21), into the solute transport equation, given by Eq. (16), and neglecting the nonlinear terms, we obtain: where we have Taylor expanded dk(c) dc around c = c eq . We can further simplify Eq. (23) and carry out the last integral in the right-hand side: and we define the perturbation to the equilibrium stretching of the dumbbells as S : We then rewrite Eq. (23) as: Next, we substitute Eq. (21) into Eq. (17) and we neglect the nonlinear terms: The transport of solute, given by Eq. (26), shows that the evolution of the perturbations to the solute concentration depends only on the dumbbell stretch S and not on the higher moments of the dumbbell distribution. To obtain an equation for S , we multiply Eq. (27) by r 2 and we integrate over the end-to-end vector space using integration by parts: Where n d is the perturbation to the local number density of dumbbells defines as n d = r ψ dr.
The last integral in the right hand side of Eq. (28) can be computed to give: which substituted in Eq. (28) yields: To obtain an equation for n d , we integrate Eq. (27) over the end-to-end vector space using integration by parts: Finally, it can be shown that the distribution of the end-to-end vector remains isotropic during the evolution of a perturbation. Eq. (18) shows that an isotropic end-to-end vector distribution only modifies the pressure but it does not generate a velocity field. As a consequence, we do not need to consider perturbations of the velocity field.
Eqs. (26), (30) and (31) represent a linear system of partial differential equations that govern the evolution of small perturbations to the homogeneous state, which can be solved using standard methods. Here we assume an infinite domain and we consider the following ansatz for the perturbations c , n d and S : Where i is the imaginary unit, k is the wavevector of the perturbation, λ is the growth rate and (26), (30) and (31) we obtain an eigenvalue problem for the growth rate λ , which has to be solved numerically. Once the parameters are specified, the numerical solution yields the eigenvalues as a function of the wavevector.
For all the parameters tested we find three real eigenvalues. The homogeneous state is stable when all the three eigenvalues are negative but it is unstable when one or more eigenvalues are positive. We find that one eigenvalue can be positive, depending on the parameters chosen. Since in this section we examine thermodynamic phase separation, the stability of the homogeneous state depends on three dimensionless numbers that do not contain any transport coefficient: the number density of dumbbells relative to that of the solute α = n d,eq /c eq , the dimensionless slope of the dumbbell spring stiffness evaluated at the homogeneous state β = (c eq /k eq )dk(c)/dc c=c eq , and the dimensionless curvature of the dumbbell spring stiffness evaluated at the homogeneous . To discriminate between stable and unstable homogeneous states, we look for the set of parameters for which at least one eigenvalue has a positive real part at any wavevector.
In appendix C we show that the stability of the homogeneous phase can also be studied considering perturbations of the free energy, given by Eq. (8). Briefly, we look at how the total free energy of the system changes when the fields c and ψ are perturbed with respect to their equilibrium values. If the perturbations reduce the total free energy, then the homogeneous system is unstable and phase separates. Using this approach, we find that the equilibrium homogeneous It follows that, the separation between stable and unstable regions in the parameter space is determined by the equation: This analytical expression can be used to verify the results obtained from the numerical solution of the eigenvalue problem. Interestingly, Eq. (36) also shows that the homogeneous equilibrium state is always stable if γ > 15/6 β 2 .
In Figure 2, we report the stability diagram of a homogeneous distribution of solute and dumb- In Figure 3, we report the effects of a nonzero second derivative of the spring stiffness with respect to the dimensionless solute concentration γ. In Figure 3(a), we find that a positive curvature γ > 0 tends to stabilize the homogeneous distribution. Compared to the case shown in Figure 2, the boundary of stability is pushed towards larger values of β . This means that a steeper slope of k(c) is required to make the homogeneous distribution unstable. Conversely, Figure 3(b) shows that a negative curvature γ < 0 tends to destabilize the homogeneous distribution. In this case, the homogeneous distribution can become unstable even in the case of β = 0. In the same plots, we show that the numerical solution of the eigenvalue problem matches perfectly the stability boundaries predicted by Eq. (36), which is shown as a red dashed line. In the present work, we do not consider any mechanism that stabilizes the phase boundaries at longer stages of the phase separation. If one is interested in the long-time and nonlinear behavior of phase separating systems, a term that penalizes the gradients of dumbbell concentration should be included in the formulation of the free energy, given by Eq. (8).
The spontaneous phase separation of a homogeneous distribution of solute and dumbbells can be understood, qualitatively, as an entropy-driven phase separation. The number of microstates that each dumbbell has access to depends on the stiffness of its spring. Dumbbells with softer springs access more microstates than dumbbells with a stiffer spring, therefore their entropy is larger. In the case of stiffness that changes significantly with the solute concentration, the dumbbells aggregate in regions where their spring is softer to maximize the total entropy. A phase-separated state, which apparently has a larger degree of spatial order, has, in fact, a larger total entropy than a homogeneous state. Finally, the results displayed in

CASE OF GENERIC FLOWS
One of the drawbacks of Eq. (17) is that it has to be solved over the space domain x and the end-to-end vector space r. This was not an obstacle to performing the linear stability analysis in the previous section but it complicates significantly numerical simulations. In this section, we discuss how this problem can be mitigated. We begin by realizing that the transport of solute, Eq.
which is proportional to the local average stretch of the dumbbells.
In principle, one needs to solve the transport equation of the dumbbells, given by Eq. (17), over both the space domain x and the end-to-end vector space r and then integrate it to compute the conformation tensor and its trace. This requires the solution of a set of partial differential equations over a six-dimensional space, which is notoriously complicated, especially using numerical simulations. Instead, inspired by previous works on polymer suspensions 42,43 , we derive an equation for the conformation tensor, C. By multiplying Eq. (17) with rr and then integrating it over the end-to-end vector space, r, we obtain: We apply integration by parts in the first integral of the right hand side of Eq. (39). After some manipulation, we obtain: In Eq. (40), the left-hand side is the upper convected time derivative of the conformation tensor and n d represents the local number density of dumbbells. In the last term on the right-hand side, the scalar product is carried out between the two spatial derivatives. Eq. (40) shows that the conformation tensor diffuses with a diffusion coefficient proportional to k B T /2ξ d . The last term in the right-hand side of Eq. (40) shows that gradients of solute also contribute to a spatial flux of the conformation tensor. Their contribution is proportional to higher moments of the endto-end vector distribution. To obtain this term, we can multiply Eq. (17) by |r| 2 rr and then integrate over the end-to-end vector space. However, this would involve again higher moments of the distribution. Repeating this procedure iteratively leads to an infinite relation between higher and higher moments of ψ, which is not practical. A similar problem is encountered in microscopic models for finite-extensive springs or rigid rods 42,63 , and has been addressed by relating higherorder moments to products of lower-order ones. Inspired by these early works, we propose the following closure approximation: where the factor 5/3 is included so that the right and the left-hand side are equal at equilibrium.
The proposed closure approximation can be interpreted as pre-averaging the term of Eq. (17) proportional to |r| 2 ψ.
By inserting this approximation in Eq. (40), we obtain: where denotes the upper convected derivative of the tensor C. Note that, in contrast to previous works 42 , the closure approximation is expected to have limited effects on the evolution of the conformation tensor. Indeed, as we discuss more in depth at the end of section VI, the relaxation of the dumbbells is usually much faster than the time required to diffuse over macroscopic lengths. In this case, the last two terms on the right-hand side of Eq. (42) where in the last term of the right-hand side the scalar product is carried out between the two spatial derivatives. Eq. (44) yields the evolution of the local number density of dumbbells where the coupling with the solute field is given by the last term on the right-hand side. We can rewrite the momentum balance using the definition of the conformation tensor as where the term inside the square brackets in the left-hand side can be identified as the divergence of the stress tensor. The equation that governs the evolution of the solute field remains unchanged by the closure approximation and can be written in terms of C as: In summary, Eqs. only. This set of equations can be solved numerically using standard numerical techniques 64 . In the Appendix D, we report Eqs. (42) -(46) using Einstein index notation.

A. Boundary conditions
The governing equations given by Eqs.
where τ N is the tangential traction imposed on the boundary. This boundary condition is complemented by the condition that the wall is impenetrable to the solvent: v · n = 0 .
Alternatively, the velocity of the solvent can be specified at the boundary: where U is the velocity at which the boundary is moving. For a stationary boundary U = 0.
The boundary of the domain is also considered to be impermeable to the solute. It follows that the component of the solute velocity normal to the wall, w, vanishes: By substituting the solute velocity given by Eq. (12) we obtain: Finally, the velocity of the dumbbells normal to the wall must vanish: By substituting the dumbbell velocity given by Eq. (13) we obtain: To obtain the boundary conditions for the local dumbbell number density, n d , and for the conformation tensor, C, we integrate the equation above over the end-to-end vector space. By integrating over r, we obtain By multiplying the boundary condition with rr, integrating and applying the closure approximation, we obtain: In the equation above, the scalar product is carried out between the normal vector and the gradient.

VI. DIMENSIONLESS EQUATIONS
To assess the importance of each term appearing in the governing equations in the case of generic flows, we make Eqs. (42) - (46) and Eq. (20) dimensionless. We assume that the flows occur over a characteristic length L, with a characteristic shear rateγ, so that the characteristic velocity scale is Lγ. We identify a characteristic dumbbell concentration n d,0 and a characteristic solute concentration c 0 , which could be the equilibrium number densities in the absence of flow.
Since the spring stiffness of the dumbbells is a function of the local solute number density, we choose a characteristic spring stiffness k 0 , which we assume to be the stiffness evaluated at c 0 , k 0 = k(c 0 ). The characteristic stress can be chosen in different ways; we choose to scale it using the solvent viscosity and the characteristic shear rate as ηγ. Finally, the components of the conformation tensor are scaled using k B T n d,0 /k 0 , which is proportional to their value at equilibrium in a homogeneous suspension.
By plugging these characteristic scales into the governing equations we obtain their dimensionless version, with the dimensionless quantities denoted with a superscript. The dimensionless momentum balance becomes where the Deborah number, De, compares the characteristic flow timescale with the relaxation time of the dumbbells and it is defined as De =γξ d /4 k 0 . In Eq. (56), we have introduced the nondimensional viscosity η * = η/η p as the ratio of the solvent viscosity and the characteristic polymer viscosity η p = k B T n d,0 ξ d /k 0 . The spring stiffness k * (c * ) and its derivative dk * (c * ) dc * are made dimensionless using the k 0 and c 0 so that k * (1) = 1. The dimensionless continuity equation The dimensionless equation of transport of solute reads: where we have introduced the solute Péclet number as Pe =γ L 2 ξ /k B T and the ratio of the dumbbell and solute characteristic concentrations n * 0 = n d,0 /c 0 . The Péclet number compares the importance of solute transport due to fluid flow to that due to diffusion. The evolution of the dimensionless dumbbell number density is given by where we defined the dumbbell Péclet number as Pe d =γ L 2 ξ d /k B T . Finally, employing the closure approximation described in section V, the equation governing the evolution of the dimensionless conformation tensor reads: Let us now assess the validity of the closure approximation that we introduced in the previous section. Eq. (60) shows that the spatial flux of the conformation tensor is proportional to De/Pe d = k B T /4k 0 L 2 . This ratio is independent of the shear rate and depends on the dumbbell relaxation time, on its friction coefficient and on the characteristic length of the problem only. It can be interpreted as the ratio between the time required to diffuse over a length L and the relaxation time of the dumbbell. Typically, the relaxation time of a polymer is in the order of 1s and its diffusion coefficient is in the order 1µm 2 s −1 . Therefore, unless the characteristic lengthscale of the problem is smaller than about 1µm, dumbbells relax faster than they diffuse, which results in De/Pe d 1.
As a consequence, we expect that in many cases the closure approximation will have a limited impact on the fluid flow.

VII. SHEAR-INDUCED PHASE SEPARATION
Motivated by the experimental evidence of shear-induced phase separation 42 , we study the stability of a homogeneously sheared suspension of solute and dumbbells to small perturbations.
We focus on homogeneous suspensions of solute, dumbbells and solvent that are thermodynamically stable and we investigate the possibility of a phase separation that is promoted by the shear flow. We address this question using linear stability analysis and a simplified theory. As shown schematically in Figure 4, we consider perturbations around a homogeneous suspension of solute and dumbbells sheared between two parallel walls. We identify the characteristic dumbbell and solute number densities as those in the homogeneous state n d,0 and c 0 in the homogeneous state.
The characteristic dumbbell spring stiffness is given by k(c 0 ). The walls are separated by a distance H and the top wall is moved along the x axis at a speed U. In this case, we identify H as characteristic length, and as characteristic shear rate,γ = U/H. We consider small perturbations to the homogeneous shear state and we assume that the perturbed variables only change along the y axis. We denote the perturbations with a prime: In the homogeneous base state the shear rate, the number densities of solute and dumbbells and the conformation tensor are constant along in the y direction. Under these conditions, the governing equations reduce to the Oldroyd-B model, which predicts the components of the conformation tensor, C * 0 , to be given by C * xx,0 = 1 + 2 De 2 , C * yy,0 = C * zz,0 = 1 and C * xy,0 = C * yx,0 = De. The linear stability analysis proceeds in the standard way. First, we substitute the perturbation in the governing equations and in the boundary conditions retaining the linear terms only. The resulting set of equations is a system of linear partial differential equations, which can be solved using eigenfunctions appropriate for the geometry considered here.
By plugging the perturbation in the dimensionless governing equations, given by Eqs (56) where we used the fact that the perturbations to the solute concentration only vary along the y axis.
The evolution of the perturbations to the dumbbell number density, n * d , is given by Since the velocity is directed along the x axis and the perturbations vary along the y axis only, there is no contribution to the transport of solute or of dumbbells due to advection. By substituting wherex andŷ are the unit vectors along the x and y axis, respectively. Since the conformation tensor is symmetric, there are four independent components that are nonzero: the three diagonal components and the xy component. Substitution of the perturbed quantities in the momentum balance yields for the x component in which only the xy component of the conformation tensor appears. The y component of the momentum balance simply states that the pressure balances variations of C * yy and of the nonlocal solute interactions along the y axis. Finally, since the perturbed velocity field v * x varies along the y axis only, the continuity equation is identically satisfied.
We assume that the walls, positioned at y * = ±1/2, are impermeable to the solute, to the dumbbells and to the solvent. The tangential velocity of the solvent on the top wall, y * = 1/2, is fixed and it is given by U * =x. The boundary conditions for the perturbations are obtained by making Eq. (51) and Eqs. (54)-(55) dimensionless and plugging in the perturbations given by Eqs. (61). The boundary condition at y * = ±1/2 for the perturbation to the solute, c * reads Since the tangential velocity is fixed at the top and bottom walls, the perturbations must vanish The perturbation to the dumbbell number density satisfies and the perturbation to the conformation tensor satisfies at y * = ±1/2.
The linear stability analysis now proceeds in the standard way 65 . We assume that the spatial and temporal dependence of the perturbed quantities can be separated. The spatial variations of the perturbations can be expanded in Fourier series and the temporal dependence is simply given by an exponential. To satisfy the no flux boundary conditions, given by Eqs (66)-(69) the perturbation of the concentration, dumbbell density and conformation tensor must be given by cosine series and where c * n , n * d,n and C * n is the initial value of the perturbation with mode n and λ n is its dimensionless growth rate. The perturbation to the velocity, which satisfies no-slip boundary conditions on the upper and bottom walls, is given by a sine series where v * x,n is the initial value of a velocity perturbation of mode n. It is straightforward to verify that the ansatz, given by Eqs (70)-(73), satisfies the boundary conditions. In the Appendix E it is shown that, by inserting Eqs (70)-(73) into Eqs (62)-(65), we obtain a set of eigenvalue problems for each growth rate λ n . Since there are six independent variables, for each Fourier mode we find six eigenvalues. The homogeneous shear state is unstable if any of the eigenvalues has a positive real part. As shown in the Appendix E, the eigenvalue problem for each mode n can be written in a compact form as: with the vector X being the vector of the initial values of the perturbations: and E being a six by six matrix. The values of λ n that satisfy Eq. (74) are given by the eigenvalues of the matrix E. In our case, there are six eigenvalues for each mode n, one for each field that is being perturbed. To investigate the stability of the homogeneous state for a given set of parameters, we have to solve the eigenvalue problem numerically for each mode n. We verified that, by solving the eigenvalue problem given by Eq. (74) for very small values of De, we correctly recover the results obtained in Section IV for a solution at equilibrium. This is expected because for De 1 the solution is essentially at equilibrium and we expect to recover the behavior studied in Section IV.
Nevertheless, analytical progress can be made by realizing that in typical situations some dimensionless number is very large or very small. In typical experimental conditions the gap between the walls is H ≈ 1mm, the shear rate of the order ofγ ≈ 1s −1 and the dumbbell diffusion coefficient k B T /ξ d ≈ 1µm 2 s −1 , which makes the dumbbell Péclet number, Pe d ≈ 10 6 , much larger than one. The diffusion coefficient of the solute molecules, k B T /ξ ≈ 10 3 µm 2 s −1 , is much larger than that of the dumbbells, therefore it is Pe/Pe d 1. As a result, there is a separation of timescales between the evolution of the perturbations of the solute number density and that of the dumbbells. It follows that the number density of dumbbells remains essentially constant over the evolution of the perturbations. Similarly, there is also a separation of timescales between the relaxation time of the dumbbells and the time required by the solute and by the dumbbells to migrate over the length of the gap: De/Pe d 1 and De/Pe 1. It follows that the conformation of the dumbbells relaxes quasi-statically to the value determined by the local shear rate and by the value of the spring stiffness. Finally, we consider as a simplification that the relative viscosity is very small η * 1, which is typically the case of polymer solutions. Under these assumptions, we consider that the dumbbell concentration remains homogeneous n * d = 0, that the time derivatives of the conformation tensor are zero and that De/Pe d = 0. By doing so (see Appendix F), we find that there exists a critical Deborah number above which the homogeneous shear is unstable There are some similarities with the equilibrium phase separation studied in Section IV. The critical Deborah number depends on the dimensionless spring stiffness slope squared, (dk * (c * )/dc * In this case, the dimensionless derivative, (dk * (c * )/dc * ) 2 = N 2 and the dimensionless curvature, shear-banding instability and the mechanism described here could play a role in the transition from the homogeneous to the shear-banded regimes.
In Figure 5, we compare the critical Deborah number predicted by the simplified expression, given by Eq. (76), and that obtained from the numerical solution of the eigenvalue problem. To In Figure 6, we show the comparison between the theory and the simulations for the case of > 0, the dumbbells contract where the solute concentration has increased and they stretch where the solute concentration has decreased (Figure 7(b)).
This results in a gradient of dumbbell stretch along the gap, which drives solute from areas where the dumbbells are more stretched to areas where they are contracted (Figure 7(c)). The solute flux reinforces the initial perturbation, ultimately leading to an instability. It is straightforward to show that the growth of the initial perturbations is driven by a decrease in the total elastic energy. In dimensionless form, the elastic energy is written as y * =1/2 By expanding the elastic energy up to second order in the perturbations, we obtain y * =1/2 Since the dumbbells relax very quickly to equilibrium and their concentration evolves slowly, by using Eq. (64) we find that Tr(C * ) ≈ −Tr( c * . By plugging this into the expansion of the elastic energy, given by Eq. (78), we obtain y * =1/2 Since the trace of the conformation tensor is a positive quantity representing the extension of the dumbbells and the solute perturbations, c * , are given by a Fourier series, it is straightforward to show that the integral is always negative. It follows that any perturbation to the solute concentration reduces the total elastic energy of the system. The very same mechanism described above takes place in the case dk * (c * )/dc * Previous works showed that instabilities in shear flow can be produced mainly through two mechanisms: (i) a nonmonotonic stress-strain rate relation 35,39 and (ii) polymer migration driven by inhomogeneous elastic stress 36,37,60,61,[66][67][68] . In contrast, our work identifies a new mechanism for shear instability that couples the accumulation of polymers and solute molecules. This mechanism could be the dominant one if the conformation of a polymer is highly sensitive to the concentration of solute molecules.

VIII. CONCLUSIONS
In this paper, we present a theory for the flow of a suspension of chemically-responsive polymers that can change conformation depending on the local concentration of a solute as is the case, for instance, of polyelectrolytes and salt. To capture the solute-dependent conformation changes while keeping the model as simple as possible, we assume that the polymers can be idealized as dumbbells linked by a linear spring whose stiffness depends on the local number density of a solute. The choice of a linear dumbbell model is expected to be a good description of dilute and We investigate the stability of a mixture of solute and dumbbells that is sheared between two parallel walls. We find that the homogeneous state undergoes spontaneous demixing in the case of sufficiently large shear rates. Using a linear stability analysis, we derive an approximated expression of the critical Deborah number above which the suspension of solute and dumbbell phase separates, which agrees with the numerical solution of the full eigenvalue problem. Our results show that the newly-identified couplings between the solute and the dumbbell distribution drive the phase separation. This mechanism is radically different from that driving the shear banding of micelles 35,39 and of entangled polymer solutions 36,37,[66][67][68] and could act in parallel with the other two mechanisms.
In summary, the model presented here is general and can be applied to study the flow of polymers that respond to different chemical stimuli such as local changes in pH, ionic strength or solute concentration changes. In future works, we plan to extend the framework to consider the finite extensibility of the spring, mean-field hydrodynamic interactions and conformation-dependent dumbbell friction coefficients. Finally, a suggestive application of the theory presented here concerns living systems whereby active processes, chemical reactions and inhomogeneous straining flows could be used to control the conformation and the spatial arrangements of macromolecules and proteins. Appendix A: Calculation of the rate of change of the free energy

ACKNOWLEDGMENTS
Here we report the details on the derivation of the rate of change of the free energy, which is obtained by taking the time derivative of Eq. (8): where we used the fact that the integration domain does not change in time and we used the chain rule ∂U/∂t = (∂U/∂ c)(∂ c/∂t). We now substitute the time derivatives of ψ and of c from the continuity equations, given by Eqs. (3)-(4), to obtain: We apply integration by parts to each of the integrals appearing in Eq. (A2). Since the velocity of dumbbells and of the solute normal to the boundaries vanishes and ψ → 0 as r → ∞, the boundary terms coming from integration by parts disappear: Finally, we use the property of the derivatives of the logarithms, which yieldṡ we rewrite Eq. (B3) in terms of w d andṙ: which is the dumbbell contribution to the dissipation used in the Eq. (10) of the main text.
Appendix C: perturbation of the free energy around equilibrium In this Appendix, we analyze the stability of a homogeneous solution of solute and dumbbells at equilibrium. We consider perturbations to the solute and dumbbell distributions and we check if these perturbations result in a reduction of the total free energy of the system. If this is the case, then small perturbations will lead to phase separation. The free energy of the system is defined in the main text as Where the elastic energy is given by U = 1 2 k(c)|r| 2 . We consider perturbations around equilibrium of the type: where c eq is a constant and ψ eq is constant in space but it depends on the magnitude of the end to end distance vector: As we showed in the main text, it is: By inserting the perturbations, given by Eqs. (C2), into the free energy and retaining only the second-order terms, we obtain the perturbation to the free energy: 4 Ω r ψ eq d 2 k(c) dc 2 c eq δ c 2 |r| 2 dΩ dr+ + c eq 2 Ω r ψ eq dk(c) dc c eq δ ψ δ c|r| 2 dΩ dr . (C5) The system is unstable if a perturbation δ c or δ ψ leads to a negative free energy change δ F < 0 To check this, we proceed by expand the perturbations fields in eigenfunctions. The expansion for δ c is given by In the equation above q is the wavevector defines as q = 2 π L {n x , n y , n z } with n x , n y and n z ∈ Z and L is the size of the domain. The amplitude of the perturbation at a given wavevector is set by c q . Since the perturbations are real-valued functions, the coefficients c q must fulfill the condition c −q = c * q where c * q is the complex conjugate of c q . We insert the decomposition of δ c, given by Eq. (C6), into the first integral on the right-hand side of Eq. (C5): where L 3 is the volume of the spatial domain and δ qq is the Kronecker delta symbol. Similarly, by inserting the decomposition of δ c given by Eq. (C6) into the third integral on the right-hand side of the free energy perturbation, Eq. (C5), we find: where we used the fact that r ψ eq |r| 2 dr = 3n d,eq k B T k eq . (C9) Now we need to assume an expansion for the perturbation of the dumbbell distribution, δ ψ, which depends not only on space but also on the end-to-end space r. We assume an expansion of the type: where L 1 2 j is the associated Laguerre polynomial of order j. The value Ψ q, j gives the amplitude of the perturbation with wavelength q and mode j. Basically, we are assuming that δ ψ can be separated in the product of two functions, one that only depends on space and one that only depends on the end-to-end vector r, we then choose the appropriate basis for these functions. For the end-to-end vector space we choose the Laguerre polynomials because they have the following orthogonality property Since the perturbations are real valued, the coefficients Ψ q, j must fulfill the condition Ψ −q, j = Ψ * q, j . Finally, the expansion in Laguerre polynomials ensures that the perturbations, ψ eq δ ψ, decay to zero as r → ∞ because ψ eq ∝ exp (−r 2 ). Therefore, they fulfill the boundary conditions at r → ∞.
By inserting the expansion given by Eq. (C11) into the second integral on the right-hand side of Eq. (C5), we obtain: where we have rewritten the integral over r in spherical coordinates and we have carried out the integration over the theta and phi coordinates. We substitute the expression for ψ eq in the righthand side to obtain k B T 2 Ω r ψ eq δ ψ 2 dΩ dr = 2 πk B T n d,eq These integrals can be carried out to give: Now, the last integral on the right hand side can be evaluate using the same expansion for of the perturbations c eq 2 Ω r ψ eq dk(c) dc c eq δ ψ δ c|r| 2 dΩ dr = 2 π c eq n d,eq 2 √ 2 k eq π k B T By using again the change of variable z = r 2 k eq 2k B T , the integrals on the right-hand side can be rewritten as c eq 2 Ω r ψ eq dk(c) dc c eq δ ψ δ c|r| 2 dΩ dr = 2 c eq n d,eq √ π k B T k eq dk(c) dc c eq ∑ q ∑ j Ψ q, j c q Ω e iq·x e iq ·x dΩ The integrals on the right-hand side can be carried out to obtain: c eq 2 Ω r ψ eq dk(c) dc c eq δ ψ δ c|r| 2 dΩ dr = 3 L 3 c eq n d,eq 4 k B T k eq dk(c) dc c eq Ψ * q, j c q + Ψ q , j c * q δ qq δ j0 − δ j1 . (C18) The perturbations to the free energy can be expressed in matrix form as where the vectors p and p * are given by and the matrix A is given by 3n eq 4c eq c 2 eq k eq d 2 k(c) dc 2 c eq 3n d,eq 4c eq c eq k eq dk(c) dc c eq − 3n d,eq 4c eq c eq k eq dk(c) dc c eq 0 0 3n d,eq 4c eq c eq k eq dk(c) dc c eq n d,eq 2c eq 0 0 0 − 3n d,eq 4c eq c eq k eq dk(c) dc c eq 0 3n d,eq 4c eq 0 0 The homogeneous equilibrium system is unstable if the determinant of the A is negative det(A) < 0.
The condition of stability can be rewritten − α 2 15 4 where, following the nomenclature used in the main text, we introduced the dimensionless groups α = n eq c eq , β = c eq k eq dk(c) dc c eq and γ = c 2 eq k eq d 2 k(c) dc 2 c eq . The critical stability is given when the determinant of the matrix is exactly zero The condition is satisfied when which leads to the relation This relation describes a curve in the parameter space that discriminates regions that are stable from regions that are unstable. In the main text we show that this expression matches the numerical results obtained from the linear stability analysis.
Appendix D: Governing equations using the closure approximation and the Einstein notation In this appendix we rewrite Eqs. (42) -(46) using Einstein index notation where repeated index are summed. We assume that a cartesian coordinate system is employed. The equation for the local conformation tensor, C i j , reads: denotes the upper convected derivative of the tensor C i j and δ i j represents the Kronecker delta.
The equation for the local dumbbell density, n d , rewritten using the Einstein notation reads: The momentum balance reads: The equation that governs the evolution of the solute using the Einstein notation: λ n + ity analysis: (i) we neglect the perturbations of dumbbell density n * d = 0, since they evolve over a very slow timescale; (ii) we assume that the conformation tensor relaxes faster than any other timescale in the system; (iii) we take the limit η * → 0. Given these assumptions, the components of the conformation tensor are given by: We substitute this expression in the equation for the evolution of the solute perturbations to obtain λ n + 2 π 2 n 2 Pe 2 + 16Λ * (n π) 2 + n * 0 (3 + 2 De 2 ) We simplify c * n on the two sides of the equation to obtain a dispersion relation between the growth rate λ n and the dimensionless parameters of the problem.
We define a critical Deborah number, De crit , as the value that leads to a marginally stable case, λ n = 0 2 + n * 0 (3 + 2 De 2 crit ) Any homogeneous state with De > De crit leads to an instability and a shear-induced phase separation. We can rearrange the equation above to write an explicit expression for De crit