format: text

text:
# Flow phase diagrams for concentration-coupled shear banding

## I Introduction

For many complex fluids, the intrinsic constitutive curve of shear stress $`\mathrm{\Sigma }`$ as a function of shear rate $`\dot{\gamma }`$ is non-monotonic, admitting multiple values of shear rate at common stress. For semi-dilute wormlike micelles, theory cates90 ; SpenCate94 ; SCM93  predicts the form ACEG of Fig. 1. In the range $`\dot{\gamma }_{\mathrm{c1}}<\dot{\gamma }<\dot{\gamma }_{\mathrm{c2}}`$ where the stress is decreasing, steady homogeneous flow (Fig. 2a) is unstable Yerushalmi70 . For an applied shear rate $`\overline{\dot{\gamma }}`$ in this unstable range, Spenley, Cates and McLeish SCM93  proposed that the system separates into high and low shear rate bands ($`\dot{\gamma }_\mathrm{h}`$ and $`\dot{\gamma }_{\mathrm{}}`$; Fig. 2b) and that any change in the applied shear rate then merely adjusts the relative fraction of the bands, while the stress $`\mathrm{\Sigma }_{\mathrm{sel}}`$ (which is common to both) remains constant. The steady state flow curve then has the form ABFG. Several constitutive models augmented with interfacial gradient terms have captured this behaviour olmsted99a ; lu99 ; olmstedlu97 ; spenley96 .

Experimentally, this scenario has been widely observed in semidilute wormlike micelles berret94b ; Call+96 ; grand97 . The steady state flow curve (which is often attained only after very long transients grand97 ) has a well defined, reproducible stress plateau $`\mathrm{\Sigma }_{\mathrm{sel}}`$. Coexistence of high and low viscosity bands has been observed by NMR spectroscopy Call+96 ; MairCall96 ; MairCall96c ; BritCall97 . Further evidence comes from small angle neutron scattering (SANS) berret94a ; schmitt94 ; Capp+97 ; rehage91 ; berret94b ; and from flow birefringence (FB) Decr+95 ; Makh+95 ; DCC97 ; BPD97 , which reveals a (quasi) nematic birefringence band coexisting with an isotropic one (but see FisCal01 ; FisCal00 ).

In some systems, the coexistence plateau is not perfectly flat, but slopes upward slightly with increasing shear rate (B’F’ in Fig. 1). See, for example, Ref. LerDecBer00  for CTAB(0.3M)/$`\mathrm{NaNO}_3(1.79\mathrm{M})`$/H2O at micellar volume fraction $`\varphi =11\%`$. This effect is much more pronounced in other, more concentrated systems that are near an underlying (zero-shear) isotropic-nematic (I-N) transition ($`\varphi 30\%`$schmitt94 ; BRL98 ; berret94a .

In a cylindrical Couette geometry, this upward slope is qualitatively consistent with the inhomogeneous stress arising from the cell curvature: as the high-shear band at the inner cylinder expands outward with increasing applied shear rate, the applied torque must increase to ensure that the interface between the bands stays at the selected stress $`\mathrm{\Sigma }_{\mathrm{sel}}`$ olmsted99a . However a more generic explanation, independent of geometry, is that the shear banding transition is coupled to concentration schmitt95 ; olmstedlu97 . In this case, the properties of each phase must change as the applied shear rate is tracked through the coexistence regime, because material is redistributed between the bands as the high shear band grows to fill the gap.

Generically, one expects flow to be coupled to concentration in viscoelastic solutions where the different constituents (polymer and solvent) have widely separated relaxation timescales brochdgen77 ; HelfFred89 ; DoiOnuk92 ; milner93 ; WPD91 ; beris94 ; SBJ97 ; tanak96 . This was explained by Helfand and Fredrickson (HF) HelfFred89  as follows. In a sheared solution, the parts of an extended polymer molecule (micelle for our purposes) in regions of lower viscosity will, upon relaxing to equilibrium, move more than the parts mired in a region of high viscosity and concentration. A relaxing molecule therefore on average moves towards the higher concentration region. This provides a positive feedback mechanism whereby micelles can move up their own concentration gradient, and leads to flow-enhanced concentration fluctuations perpendicular to the shear compression axis. This was observed in steadily sheared polymer solutions in the early 1990’s WPD91 . In a remarkable paper, Schmitt et al. schmitt95  discussed the implications of this feedback mechanism for the onset of flow instabilities. Strongly enhanced concentration fluctuations were subsequently observed in the early time kinetics of the shear banding instability in Ref. DecLerBer01 .

Recently, therefore, we introduced a model of concentration-coupled shear banding FieOlm02 ; FieOlm02b  by combining the diffusive Johnson Segalman (d-JS) model johnson77 ; olmsted99a  with a two-fluid approach brochdgen77 ; milner91 ; deGen76 ; Brochard83  to concentration fluctuations. This “d-JS-$`\varphi `$” model does not address the microscopics of any particular viscoelastic system, but instead should be regarded as a minimal model that combines (i) a constitutive curve like that of semidilute wormlike micelles (Fig. 1) with (ii) the non-local (interfacial) terms required for selection of a unique banded state lu99  and (iii) a simple approach to concentration coupling.

In Refs. FieOlm02 ; FieOlm02b , we examined the linear stability of initially homogeneous shear states in this d-JS-$`\varphi `$ model with respect to coupled fluctuations in shear rate $`\dot{\gamma }`$, micellar strain $`\underset{¯}{\underset{¯}{W}}`$ and concentration $`\varphi `$. We thereby calculated the “spinodal”, inside which such homogeneous states are unstable. We also calculated the selected length scale at which inhomogeneity first emerges during startup flows in the unstable region. In the limit of zero concentration coupling, the unstable region coincides with that of negative slope in the homogeneous constitutive curve, as expected; but no length scale is selected during startup. Concentration coupling enhances this instability at short length scales. It thereby broadens the region of instability, and selects a length scale at which inhomogeneity must emerge.

In the present paper, we compute the corresponding steady-state flow phase diagram (the “binodals” and their tie-lines). As far as we are aware, this is the first concrete calculation aimed at qualitatively describing concentration-coupled shear banding for systems such as semi-dilute wormlike micelles. We start in Sec. II by describing the experimental background in more detail. We also compare our present calculation with the only other existing one for concentration-coupled shear banded states, in concentrated solutions of rigid rods olmsted99c . In Sec. III we summarize our d-JS-$`\varphi `$ model. We then review our results for the spinodal onset of instability in Sec. IV. In Sec. V we describe our numerical procedure for computing the banded steady states, with brief discussion of our careful study of mesh and finite size effects. We then (Sec. VI) present our results for the flow phase diagrams and shear-banded profiles. We conclude in Sec. VII.

## II Experimental background; theoretical context

In this section, we discuss in more detail the experimental evidence for concentration coupling in the shear banding of wormlike micelles. We survey both (i) concentrated systems near the zero-shear I-N phase transition and (ii) semidilute systems, in which underlying nematic interactions are likely to be less important. Correspondingly, we compare the present calculation (aimed at the semi-dilute systems) with an earlier calculation of flow phase diagrams in concentrated solutions of rigid rods (near the I-N transition) olmsted99c .

The earliest observations of an upwardly sloping stress plateau in wormlike micelles were made by Schmitt et al. schmitt94  and Berret et al. berret94a ; BRL98 . Schmitt et al. schmitt94  studied CpClO3/NaClO3(0.05M)/H2O at the high micellar volume fraction $`\varphi 31\%`$, just below the onset of the I-N transition at $`\varphi _{\mathrm{}}=34\%`$. In the steady-state flow curve, the stress increased smoothly up to the critical shear rate $`\dot{\gamma }_{\mathrm{}}`$, where it showed a pronounced downward kink before curving upward again for $`\dot{\gamma }>\dot{\gamma }_{\mathrm{}}`$ (qualitatively like B’F’ in Fig. 1). SANS measurements confirmed a superposition of nematic and isotropic contributions in this regime $`\dot{\gamma }>\dot{\gamma }_{\mathrm{}}`$, with the nematic contribution rising linearly from zero at $`\dot{\gamma }=\dot{\gamma }_{\mathrm{}}`$.

Berret et al. berret94a ; BRL98  studied CpCl/hexanol/NaCl for several micellar volume fractions, again at a volume fraction just below the onset of the zero-shear I-N transition ($`\varphi _{\mathrm{}}32\%`$). The overall height of the coexistence plateau (which again sloped upwards in $`\dot{\gamma }`$) was found to fall with increasing surfactant concentration $`\varphi \varphi _{\mathrm{}}`$, extrapolating to zero at $`\varphi \stackrel{>}{}\varphi _{\mathrm{}}`$, which is already biphasic in zero shear. They also found an increasing nematic contribution to SANS patterns for increasing shear rates above $`\dot{\gamma }_{\mathrm{}}`$. They further used the SANS data to show that the nematic (high shear) band was more concentrated than the low shear band.

As noted above, the majority of existing calculations of shear-banded states have assumed uniform concentration. An important exception is the calculation of Olmsted and Lu olmsted99c . Although this model was aimed at concentrated solutions of rigid rods, it broadly captured some of the experimental phenomenology for the concentrated ($`\varphi 30\%`$) wormlike micelles schmitt94 ; berret94a ; BRL98 . For example, the overall height of the coexistence plateau $`\mathrm{\Sigma }_{\mathrm{sel}}`$ increased from zero as the concentration was reduced below the threshold $`\varphi _{\mathrm{}}`$ of the zero-shear I-N biphasic regime. The coexistence plateau sloped upward markedly in shear rate. In further agreement with experiment, the high shear (nematic) phase had a higher volume fraction of rods. It should be noted that model of Ref. olmsted99c  was explicitly aimed at concentrated systems, which in zero shear are already close to the I-N transition: hence, the dynamics of the relevant order parameter $`\underset{¯}{\underset{¯}{Q}}`$ was driven by a free energy that already contained a phase transition. In contrast, the simple free energy $`F_\mathrm{e}(\underset{¯}{\underset{¯}{W}})`$ we consider below has no underlying phase transition, and flow-induced efects are driven by convective, rather than dissipative (relaxational) dynamics.

There have also been several experimental studies of concentration dependence in the shear banding of more dilute wormlike micellar solutions. For example, Berret et al. berret94b  investigated the non-linear rheology of CpCl/NaSal/brine in the concentration range $`5\%20\%`$, well below the I-Nc transition at $`\varphi 36\%`$. (Nc is “nematic calamitic”.) In contrast to the more concentrated “prenematic” systems, the plateau height $`\mathrm{\Sigma }_{\mathrm{sel}}`$ decreased with decreasing concentration: see the left 5 points in Fig. 3. The width of the plateau also decreased so that the difference $`\dot{\gamma }_\mathrm{h}\dot{\gamma }_{\mathrm{}}`$ fell to zero at a critical point $`\overline{\varphi }_\mathrm{c},\overline{\dot{\gamma }}_\mathrm{c},\mathrm{\Sigma }_\mathrm{c}`$ (leftmost point in Fig. 3). (These trends are the same as those in Fig. 6(d) below.) In contrast the scaled plateau height $`\mathrm{\Sigma }_{\mathrm{sel}}/G(\varphi )`$ (where $`G`$ is the plateau modulus) decreased with increasing concentration and extrapolated to zero in the zero-shear biphasic (I-N) regime at $`\varphi \stackrel{>}{}\varphi _{\mathrm{}}33\%`$. According to this extrapolation (which is actually well beyond the final data point at $`\varphi 22\%`$), and in the absence of a divergence in $`G(\varphi )`$, the unscaled plateau height $`\mathrm{\Sigma }_{\mathrm{sel}}`$ must itself fall to zero at $`\varphi 33\%`$ (rightmost data point in Fig. 3), consistent with the behaviour of the concentrated systems discussed above. To summarise, the plateau height $`\mathrm{\Sigma }_{\mathrm{sel}}`$ appears to be a non-monotonic function of concentration, increasing with $`\varphi `$ through the studied regime $`5\%<\varphi <20\%`$ before (probably) falling to zero in the zero-shear biphasic $`IN`$ regime $`\varphi 33\%`$.

Although this experiment showed that shear banding depends on the overall concentration of the solution, there is relatively little evidence for concentration-coupling (i.e. concentration differences between the bands) in such dilute systems, far from the I-N transition. Indeed, the experiment just described revealed no discernible upward slope in the coexistence plateau. We are not aware of any measurements of concentration differences between the coexisting shear bands in such systems. Nonetheless, recent experiments on CTAB(0.3M)/$`\mathrm{NaNO}_3(1.79\mathrm{M})`$/H2O at $`\varphi =11\%`$ LerDecBer00  did reveal a stress plateau with slight upward slope. Along with the generic expectation that flow should be coupled to concentration in these viscoelastic solutions, this suggests that an explicit calculation of concentration difference in the shear bands of systems far from an I-N transition might be worthwhile.

In this paper, therefore, we present the first such calculation, using our d-JS-$`\varphi `$ model FieOlm02 ; FieOlm02b . In contrast to the work of Olmsted et al. olmsted99c  for concentrated rigid rods, shear banding in the d-JS-$`\varphi `$ model is not due to any underlying nematic feature of the elastic free energy $`F^\mathrm{e}(\underset{¯}{\underset{¯}{W}})`$. Instead the instability results mainly from the non-linear effects of shear (the intrinsic constitutive curve has a region of negative slope), though it can be strongly enhanced by concentration coupling in systems close to an underlying Cahn-Hilliard (CH) demixing instability (governed by the osmotic free energy $`F^\mathrm{o}(\varphi )`$). Indeed, the d-JS-$`\varphi `$ model captures a broad crossover between (i) instabilities that are mainly mechanical (governed by the negative slope of the flow curve) and (ii) instabilities that are essentially CH demixing (governed by $`F^\mathrm{o}(\varphi )`$), but now triggered by shear. \[Likewise, in practice there should be no sharp distinction between concentrated micelles with an underlying nematic feature in $`F^\mathrm{e}(\underset{¯}{\underset{¯}{W}})`$ on the one hand and “non-nematic” (more dilute) systems on the other: any more refined model should allow a smooth crossover between the two cases. This will be the focus of a future publication FieOlm02f .\]

## III Model

In this section we outline the d-JS-$`\varphi `$ model, which couples shear banding instabilities to concentration in a simple way by combining the non-local Johnson-Segalman (d-JS) model johnson77  with a 2-fluid framework brochdgen77 ; milner93  for concentration fluctuations. While this description is self-contained, readers are referred to Ref. FieOlm02b  for fuller details.

### III.1 Free energy

In a sheared fluid, one cannot strictly define a free energy. Nonetheless, for realistic shear rates, many internal degrees of freedom of a polymeric solution relax quickly on the timescale of the moving constraints and are therefore essentially equilibrated. Integrating over these fast variables, one obtains a free energy for a given fixed configuration of the slow variables. For our purposes, the relevant slow variables are the fluid momentum and micellar concentration $`\varphi `$ (which are both conserved and therefore truly slow in the hydrodynamic sense), and the micellar strain $`\underset{¯}{\underset{¯}{W}}`$ that would have to be reversed in order to relax the micellar stress (which is slow for all practical purposes):

$$W_{\alpha \beta }=\frac{R_\alpha ^{}}{\underset{¯}{R}}\frac{R_\beta ^{}}{\underset{¯}{R}}\delta _{\alpha \beta }$$
(1)
where $`\delta \underset{¯}{R}^{}`$ is the deformed vector corresponding to the undeformed vector $`\delta \underset{¯}{R}`$.

The resulting free energy is assumed to comprise separate osmotic and elastic components,

$$F=F^\mathrm{o}(\varphi )+F^\mathrm{e}(\underset{¯}{\underset{¯}{W}},\varphi ).$$
(2)
The osmotic component is

$`F^\mathrm{o}(\varphi )`$ $`=`$ $`{\displaystyle d^3x\left[f(\varphi )+\frac{g}{2}(\underset{¯}{}\varphi )^2\right]}`$ (3)
$``$ $`\frac{1}{2}{\displaystyle d^3q(1+\xi ^2q^2)f^{\prime \prime }|\varphi (q)|^2},`$
where $`f^{\prime \prime }`$ is the osmotic susceptibility and $`\xi `$ is the equilibrium correlation length for concentration fluctuations. The elastic component is

$$F^\mathrm{e}(\underset{¯}{\underset{¯}{W}},\varphi )=\frac{1}{2}d^3xG(\varphi )\mathrm{tr}\left[\underset{¯}{\underset{¯}{W}}\mathrm{log}(\underset{¯}{\underset{¯}{\delta }}+\underset{¯}{\underset{¯}{W}})\right]$$
(4)
in which $`G(\varphi )`$ is the micellar stretching modulus.

### III.2 Dynamics

The basic assumption of the two-fluid model is a separate force balance for the micelles (velocity $`\underset{¯}{v}_\mathrm{m}`$; volume fraction $`\varphi `$) and the solvent (velocity $`\underset{¯}{v}_\mathrm{s}`$) within any element of solution.

The micellar force balance equation is:

$$\rho _\mathrm{m}\varphi (_t+\underset{¯}{v}_\mathrm{m}.\underset{¯}{})\underset{¯}{v}_\mathrm{m}=\underset{¯}{}.G(\varphi )\underset{¯}{\underset{¯}{W}}\varphi \underset{¯}{}\frac{\delta F(\varphi )}{\delta \varphi }+2\underset{¯}{}.\varphi \eta _\mathrm{m}\underset{¯}{\underset{¯}{D}}_\mathrm{m}^0\zeta (\varphi )\underset{¯}{v}_{\mathrm{rel}}\varphi \underset{¯}{}p.$$
(5)
In this equation, $`G(\varphi )\underset{¯}{\underset{¯}{W}}2G(\varphi )\underset{¯}{\underset{¯}{W}}.\frac{\delta F}{\delta \underset{¯}{\underset{¯}{W}}}`$ is the viscoelastic micellar backbone stress due to deformation of the local molecular strain, while the osmotic stress $`\frac{\delta F^\mathrm{o}}{\delta \varphi }`$ results from direct monomeric interaction. The Newtonian stress $`2\varphi \eta _\mathrm{m}\underset{¯}{\underset{¯}{D}}_\mathrm{m}^0`$ describes fast micellar processes (e.g. Rouse modes) with $`\underset{¯}{\underset{¯}{D}}_\mathrm{m}^0`$ the traceless symmetric micellar strain rate tensor. The force $`\zeta \underset{¯}{v}_{\mathrm{rel}}`$, where $`\underset{¯}{v}_{\mathrm{rel}}=\underset{¯}{v}_\mathrm{m}\underset{¯}{v}_\mathrm{s}`$, impedes relative motion; $`\zeta `$ is the drag coefficient (Eq. 8). Incompressibility determines the pressure $`p`$.

Likewise, the solvent force balance comprises the Newtonian viscous stress, the drag force (equal and opposite to the drag on the micelles), and the hydrostatic pressure:

$$\rho _\mathrm{s}(1\varphi )(_t+\underset{¯}{v}_\mathrm{s}.\underset{¯}{})\underset{¯}{v}_\mathrm{s}=2\underset{¯}{}.(1\varphi )\eta _\mathrm{s}\underset{¯}{\underset{¯}{D}}_\mathrm{s}^0+\zeta (\varphi )\underset{¯}{v}_{\mathrm{rel}}(1\varphi )\underset{¯}{}p.$$
(6)
Equations 5 and 6 contain the basic assumption of “dynamical asymmetry”, i.e. that the viscoelastic stress acts only on the micelles and not on the solvent. Adding them, and assuming equal mass densities $`\rho _\mathrm{m}=\rho _\mathrm{s}\rho `$, we obtain the overall force balance equation for the centre of mass velocity, $`\underset{¯}{v}=\varphi \underset{¯}{v}_\mathrm{m}+(1\varphi )\underset{¯}{v}_\mathrm{s}`$:

$$\rho (_t+\underset{¯}{v}.\underset{¯}{})\underset{¯}{v}D_t\underset{¯}{v}=\underset{¯}{}.G(\varphi )\underset{¯}{\underset{¯}{W}}\varphi \underset{¯}{}\frac{\delta F(\varphi )}{\delta \varphi }+2\underset{¯}{}.\varphi \eta _\mathrm{m}\underset{¯}{\underset{¯}{D}}_\mathrm{m}^0+2\underset{¯}{}.(1\varphi )\eta _\mathrm{s}\underset{¯}{\underset{¯}{D}}_\mathrm{s}^0\underset{¯}{}p.$$
(7)
Subtracting them (with each predivided by its own volume fraction), we obtain an expression for the relative motion $`\underset{¯}{v}_{\mathrm{rel}}=\underset{¯}{v}_\mathrm{m}\underset{¯}{v}_\mathrm{s}`$, which in turn specifies the concentration fluctuations:

$$D_t\varphi =\underset{¯}{}\varphi (1\varphi )\underset{¯}{v}_{\mathrm{rel}}=\underset{¯}{}\frac{\varphi ^2(1\varphi )^2}{\zeta (\varphi )}\left[\frac{\underset{¯}{}G(\varphi )\underset{¯}{\underset{¯}{W}}}{\varphi }\underset{¯}{}\frac{\delta F}{\delta \varphi }+\frac{2\underset{¯}{}\varphi \eta _\mathrm{m}\underset{¯}{\underset{¯}{D}}_\mathrm{m}^0}{\varphi }\frac{2\underset{¯}{}(1\varphi )\eta _\mathrm{s}\underset{¯}{\underset{¯}{D}}_\mathrm{s}^0}{1\varphi }\right]$$
(8)
which defines the micellar diffusion coefficient $`Df^{\prime \prime }(\varphi )\varphi ^2(1\varphi )^2/\zeta `$. We have omitted negligible inertial corrections to Eqs. (7) and (8) FieOlm02b .

The essence of the 2-fluid model is that the physically distinct elastic and osmotic stresses appear together in the force-balance equation (7) and also in the generalised CH equation (8). This allows micellar diffusion in response to gradients in concentration and in the viscoelastic stress. We will see below that this gives rise to a positive HF feedback between concentration and flow HelfFred89 , allowing micelles to diffuse up their own concentration gradient.

For the dynamics of the viscoelastic micellar backbone strain we use the phenomenological d-JS model johnson77 ; olmsted99a :

$$(_t+\underset{¯}{v}_\mathrm{m}.\underset{¯}{})\underset{¯}{\underset{¯}{W}}=a(\underset{¯}{\underset{¯}{D}}_\mathrm{m}.\underset{¯}{\underset{¯}{W}}+\underset{¯}{\underset{¯}{W}}.\underset{¯}{\underset{¯}{D}}_\mathrm{m})+(\underset{¯}{\underset{¯}{W}}.\underset{¯}{\underset{¯}{\mathrm{\Omega }}}_\mathrm{m}\underset{¯}{\underset{¯}{\mathrm{\Omega }}}_\mathrm{m}.\underset{¯}{\underset{¯}{W}})+2\underset{¯}{\underset{¯}{D}}_\mathrm{m}\frac{\underset{¯}{\underset{¯}{W}}}{\tau (\varphi )}+\frac{l^2}{\tau (\varphi )}\underset{¯}{}^2\underset{¯}{\underset{¯}{W}},$$
(9)
where $`2\underset{¯}{\underset{¯}{\mathrm{\Omega }}}_\mathrm{m}=\underset{¯}{}\underset{¯}{v}_\mathrm{m}(\underset{¯}{}\underset{¯}{v}_\mathrm{m})^T`$ with $`(\underset{¯}{}\underset{¯}{v}_\mathrm{m})_{\alpha \beta }_\alpha (v_\mathrm{m})_\beta `$. $`\tau (\varphi )`$ is the Maxwell time and $`l`$ is a length scale discussed in Sec. III.4 below. The slip parameter $`a`$ measures the non-affinity of the molecular deformation, i.e. the fractional stretch of the polymeric material with respect to that of the flow field. For $`|a|<1`$ (slip) the intrinsic constitutive curve in planar shear is capable of the non-monotonicity of Fig. 1.

We use Eqns. 78 and 9, together with the incompressibility condition, $`\underset{¯}{}.\underset{¯}{v}=0`$, as our model for the remainder of the paper.

### III.3 Flow geometry. Boundary conditions

We consider idealised planar shear bounded by infinite plates at $`y=\{0,L\}`$ with $`(\underset{¯}{v},\underset{¯}{}v,\underset{¯}{}\underset{¯}{v})`$ in the $`(\underset{¯}{\overset{^}{x}},\underset{¯}{\overset{^}{y}},\underset{¯}{\overset{^}{z}})`$ directions. We allow variations only in the flow-gradient direction, and therefore set all other derivatives to zero: $`_x\mathrm{}=0`$, $`_z\mathrm{}=0`$. In appendix A we give all the relevant components of the model equations 78 and 9 in this coordinate system.

The boundary conditions at the plates are as follows. For the velocity we assume there is no slip. For the concentration we assume

$$_y\varphi =_y^3\varphi =0,$$
(10)
which ensures (in zero shear at least) zero flux of concentration at the boundaries. Following Ref. olmsted99a , for the micellar strain we assume

$$_yW_{\alpha \beta }=0\alpha ,\beta .$$
(11)
Conditions 10 and 11 together ensure zero concentration flux at the boundary even in shear. For the controlled shear rate conditions assumed throughout,

$$\overline{\dot{\gamma }}=_0^L𝑑y\dot{\gamma }(y)=\mathrm{constant}.$$
(12)
### III.4 The interfacial terms

The model contains two different interfacial terms. The first is the gradient term on the RHS of Eqn. (9). The length $`l`$ in this term could, for example, be set by the mesh size or by the equilibrium correlation length for concentration fluctuations. Here we assume the former, since the dynamics of the micellar conformation are more likely to depend on gradients in molecular conformation than in concentration. Physically, one can interpret the gradient term in equation 9 as resulting dynamically, from the diffusion of stretched molecules across the interface elkareh89 , or statically, from nematic interactions between the micelles, or both. There is, at present, no accepted theory for these gradient terms in semi-dilute solutions. The equilibrium correlation length $`\xi `$ of course still enters our analysis through our second interfacial term, in the osmotic free energy of Eqn. 3.

Together, $`l`$ and $`\xi `$ set the length scale of any interfaces. Throughout this paper, we study the physical limit in which $`l`$ and $`\xi `$ are small compared to the system size so that we have a sharp interface connecting two bulk homogeneous phases. In this case, the solution to Eqns. 78 and 9 naturally fits the zero-gradient boundary conditions, and is invariant under $`yy/2`$, $`ll/2`$ and $`\xi \xi /2`$. Therefore, a simultaneous reduction in $`l`$ and $`\xi `$ by the same factor only changes the overall length of the interface, and not the values of the order parameters in each phase (which determines the phase diagram). However the phase diagram does depend slightly on the ratio $`r=l/\xi `$: below we will give results for $`r=0,r=\mathrm{}`$ and $`r=O(1)`$. This provides a concrete example of the early insight of Lu and co-workers lu99 , that the banded state must depend on the nature of the interfacial terms. This contrasts notably with equilibrium phase coexistence, in which the dynamical equations are integrable and therefore insensitive to interfaces.

### III.5 Model parameters

The d-JS-$`\varphi `$ model (Eqns. 789) has the following parameters: the solvent viscosity $`\eta _\mathrm{s}`$ and density $`\rho `$; the plateau modulus $`G`$; the Maxwell time $`\tau `$; the Rouse viscosity $`\eta _\mathrm{m}`$; the mesh size $`l`$; the osmotic modulus $`f^{\prime \prime }(\varphi )`$ and the equilibrium correlation length $`\xi `$ (recall Eqn. 3); the drag coefficient $`\zeta `$ and the slip parameter $`a`$. We also need to know the typical rheometer gap, $`L`$. A reference set of parameter values at $`\varphi =0.11`$ is summarised in table 1. These values were taken from experiment or calculated using scaling arguments: see Ref. FieOlm02b  for details. Note that explicit data is not available for $`f^{\prime \prime }(\varphi )`$; however dynamic light scattering gives the diffusion coefficient

$$D\frac{f^{\prime \prime }(\varphi )\varphi ^2(1\varphi )^2}{\zeta (\varphi )}.$$
(13)
In this paper we will be guided by these parameter values, but subject to the following considerations.

First, we are only interested in steady states so for convenience can take the limit of zero Reynolds number ($`\rho =0`$) and rescale the kinetic coefficient $`1/\zeta `$ so that the diffusion time $`L^2/D`$ is of order the Maxwell time. These choices have no effect on the steady state, but make our numerical calculation of it much more efficient (by evolving the dynamical equations 78 and 9). Second, realistic interfaces are much narrower than the typical rheometer gap, with $`l`$ and $`\xi `$ both of $`O(10^4L)`$. To resolve such interfaces (allowing a minimal 10 numerical mesh points per interface) would therefore require $`O(10^5)`$ grid points, while in practice we are limited to $`O(10^2)`$. We will therefore use artificially large values of $`l`$ and $`\xi `$. However this does not affect the phase diagram, provided the interface is still small compared with the gap size: see Sec. V for details, and Fig. 5 for a typical banded profile. Finally, we artificially increase the Rouse viscosity $`\eta _\mathrm{m}`$ by a factor $`50`$ to ensure, again for numerical convenience, that the shear rate of the high shear phase is not too large. This does quantitatively change the phase diagram, but we checked that the qualitative trends are not affected.

Exploring this large parameter space is a daunting prospect so we shall not, in general, vary the parameters independently of each other. Instead we simply tune the concentration $`\varphi `$, relying on known semi-dilute scaling laws for the $`\varphi `$-dependence of the other parameters (column 4 of table 1). However we will, in separate $`\varphi `$sweeps, vary the degree of concentration coupling, which is dictated by ratio of the elastic term $`\underset{¯}{}.G(\varphi )\underset{¯}{\underset{¯}{W}}`$ to the osmotic term $`\underset{¯}{}\frac{\delta F}{\delta \varphi }`$ and which we encode in the parameter

$$\alpha \frac{G^{}(\varphi =0.11)}{2f^{\prime \prime }(\varphi =0.11)}$$
(14)
(where a prime denotes a derivative). In other $`\varphi `$sweeps we will vary the characteristic interface widths $`l(\varphi =0.11)`$ and $`\xi (\varphi =0.11)`$, to investigate any dependence of the phase diagram on the ratio $`l/\xi `$ in the double limit $`l/L0`$, $`\xi /L0`$. In what follows, we adopt the convenient shorthand of $`l`$ for $`l(\varphi =0.11)`$ with the understanding that $`l`$ does actually varies with $`\varphi `$ according to the scaling given in table 1. We do likewise for $`\xi `$.

Throughout we rescale stress, time and length so that $`G(\varphi =0.11)=1`$, $`\tau (\varphi =0.11)=1`$, and $`L=1`$.

## IV Intrinsic flow curves; spinodals

The homogeneous intrinsic steady state flow curves $`\mathrm{\Sigma }(\overline{\dot{\gamma }},\overline{\varphi })=G(\overline{\varphi })W_{xy}+\eta (\overline{\varphi })\overline{\dot{\gamma }}`$ that satisfy $`_t\underset{¯}{v}=_t\varphi =_t\underset{¯}{\underset{¯}{W}}=0`$ are shown as dotted lines in Fig. 4. (The average viscosity $`\eta (\varphi )\varphi \eta _\mathrm{m}+(1\varphi )\eta _\mathrm{s}`$.) The region of negative slope ends at a “critical” point $`\overline{\varphi }_\mathrm{c}0.015`$. CPCl/NaSal in brine berret94b  shows the same trend. For completeness, in App. B we give analytical results for the steady state conditions in homogeneous shear flow.

In Ref. FieOlm02 ; FieOlm02b  we linearised in fluctuations about these homogeneous states to find the spinodal region in which the homogeneous states are unstable. The spinodals are shown in Fig. 4 for different levels concentration coupling, $`\alpha `$.

In the limit of zero concentration coupling $`\alpha 0`$, fluctuations in the “mechanical variables”, $`\underset{¯}{\underset{¯}{W}}`$ and $`\dot{\gamma }`$ decouple from those in concentration, and are unstable in the region of negative constitutive slope, as expected. Separately, the concentration could have its own Cahn-Hilliard demixing instability, when the diffusion coefficient $`D<0`$; however we are interested only in flow-induced instabilities and set $`D>0`$ throughout. For finite $`\alpha >0`$, the region of mechanical instability is broadened by coupling to the concentration fluctuations, as seen in Fig. 4. This can be understood as follows. Consider the first term in the square brackets of Eqn.(8). This causes micelles to move up gradients in the viscoelastic stress $`\underset{¯}{\underset{¯}{W}}`$, thereby increasing the concentration in stressed regions. If $`G^{}(\varphi )>0`$ (assumed here), the increased concentration causes the stress to increase further, closing a positive HF HelfFred89  feedback loop whereby the micelles can diffuse up their own concentration gradient.

## V Numerical details

In this section, we outline our numerical procedure for solving the dynamical equations 78 and 9 and discuss our careful study of time-step, mesh size and finite size effects. Readers who are not interested in these issues can skip this section.

We consider variations only in the flow gradient direction, in which we discretise $`y0,1`$ on an algebraic grid $`y_n=n/N_y`$ for $`n=0,1\mathrm{}N_y`$. We stored $`\varphi `$ and $`\underset{¯}{\underset{¯}{W}}`$ on these grid points. The velocities $`\underset{¯}{v}_\mathrm{m}`$ and $`\underset{¯}{v}_\mathrm{s}`$ were stored on half grid points $`y_{n+1/2}`$, and we used linear interpolation between the half and full grid points. Likewise we discretized time such that $`t_n=n\mathrm{\Delta }t`$. We evolved the discretized equations 78 and 9 using the Crank-Nicholson algorithm which is semi-implicit in time, with centred space derivatives

For each run, we seeded an initial profile that was either homogeneous up to a small random contribution, or inhomogeneous according to $`\varphi =\overline{\varphi }[1+\mathrm{\Delta }\mathrm{cos}(\pi x)]`$ (with $`\mathrm{\Delta }0.1`$). We then evolved the discretized equations 78 and 9 under an imposed wall velocity until a steady banded state was reached. We checked that the homogeneous phases between the interfaces were insensitive to the initial conditions. However, for the random initial condition several bands could form (and did not coarsen over any accessible timescale). Therefore in most runs we used the co-sinusoidal initial profile, to conveniently obtain just two bands (as in Fig. 12).

For the dynamics to be independent of time-step, a very small time-step has to be used. However the numerically-attained steady is much less sensitive hence allowing much larger time-steps. A typical steady state changes by less than $`10^3\%`$ for a factor-two reduction in timestep 5. For the special case of $`\xi =0`$, timesteps $`\mathrm{\Delta }tN_y^2`$ can be used, since the highest spatial derivative is second order. For $`\xi 0`$, we have a fourth order derivative in Eqn. 8 and much smaller timesteps $`\mathrm{\Delta }tN_y^4`$ must be used.

In all our calculations, we are interested in the physical limit where the interface width is much smaller than the rheometer gap. This creates a delicate balance, since narrow interfaces require a very fine grid. Therefore we adopted the following procedure. For any fixed value of the interfacial lengthscales $`l`$ and $`\xi `$, we performed several runs with progressively finer meshes (but always with a small enough time-step) until the shear banded profile and selected stress didn’t depend on the mesh. This is quite easy to achieve: a typical steady state presented below changes by less than $`0.1\%`$ upon doubling the number of mesh points. We then reduced $`l`$ and $`\xi `$ (in fixed ratio) until the order parameters in the homogeneous phases changed by less than $`0.5\%`$ upon further halving of $`l`$ and $`\xi `$ (but always ensuring convergence with respect to the number of grid points). A sample study of these issues is presented in Fig. 5 for the special case $`\xi =0`$.

## VI Results

We now present our results for the steady-state flow phase diagrams, flow curves and shear banded profiles. Because one of our aims is to show that the shear-banded state depends on the nature of the interfacial terms, we consider three separate cases: (A) interfacial terms only in the viscoelastic constitutive equation 9 ($`l0,\xi =0`$, $`rl/\xi =\mathrm{}`$); (B) interfacial terms in both the constitutive and concentration equations, 9 and 8 ($`l0,\xi 0`$, $`r=O(1)`$); and (C) interfacial terms only in the concentration equation 8 ($`l=0,\xi 0,r=0`$).

### VI.1 Interfacial terms only in the viscoelastic constitutive equation: $`l0`$, $`\xi =0`$.

In this section, we set the correlation length for concentration fluctuations, $`\xi `$, to zero and consider small but non-zero values the interfacial lengthscale $`l`$ in the constitutive equation 9.

#### VI.1.1 Flow phase diagrams

For any given shear banded profile, the values of the order parameters in each of the two homogeneous phases specify the two ends of one tie line in the phase diagram. Analogously to equilibrium tie lines, the concentrations and strain rates of the coexisting states are related to the mean strain rate $`\overline{\dot{\gamma }}`$ and mean concentration $`\overline{\varphi }`$ by the lever rule,

$`\overline{\varphi }`$ $`=\beta \varphi _1+(1\beta )\varphi _2`$ (15)
$`\overline{\dot{\gamma }}`$ $`=\beta \dot{\gamma }_1+(1\beta )\dot{\gamma }_2,`$ (16)
where $`\beta `$ is the volume fraction of material in state $`(\varphi _1,\dot{\gamma }_1)`$. For each of several values of the concentration coupling, $`\alpha `$, we calculated the full phase diagram via a succession of shear startup runs, all at the critical shear rate $`\overline{\dot{\gamma }}_\mathrm{c}(\alpha )`$ (determined from Fig. 4), for average concentrations ranging from $`\overline{\varphi }=0.15`$ down to the critical value $`\overline{\varphi }_\mathrm{c}(\alpha )`$. For concentrations below the critical point the response of the system is smooth as a function of stress. In our model, this arises because decreasing concentration reduces the viscosity of the low shear rate branch faster than it reduces the viscosity of high shear rate branch. Hence the stress maximum decreases with decreasing concentration, disappearing when the stress maximum vanishes. Alternatively, in a more dilute system the plateau modulus and Maxwell time are both smaller, and one expects a smaller stress and higher strain rate at the onset of instability.

The results for $`\alpha =10^2`$, which gives rather strong concentration coupling, are shown in Fig. 6a,b,c. Because the width, $`\delta `$, of the interface in the banded state is set by $`l`$, but with a prefactor that diverges at the critical point, in each successive run we rescaled $`l`$ so that $`\delta `$ remained (approximately) equal to its value ($`L`$) in the first run at $`\overline{\varphi }=0.15`$. We return below to study the divergence of $`\delta /l`$ at the critical point.

To illustrate the finite size considerations of Sec. V (above), in Fig. 6(a) we show the tie lines obtained for two different (starting) values of $`l`$. All the results are converged with respect to mesh fineness and timestep (not explicitly shown), but the tie lines differ slightly between the two values of $`l`$. However all seem to be consistent with one given binodal line: we do not have any explanation for this apparent consistency.

To investigate the effect of reducing the coupling to concentration, we repeat the phase diagram for $`\alpha =10^2`$ alongside that for $`\alpha =10^3`$ and $`\alpha =10^4`$ in Fig. 7. As expected, the concentration difference between the bands tends to zero as $`\alpha 0`$.

#### VI.1.2 Flow curves

So far, we have discussed the flow phase diagrams. Measurement of these diagrams still presents an open challenge to experimentalists, due to the difficulty in measuring the concentration of micelles in each band (although SANS data has been used to estimate the bands’ concentrations in systems near the I-N transition BRL98 ). In this section we discuss the macroscopic flow curves, which are relativly easily measured using conventional bulk rheology. However it is important to realise that a set of flow curves $`\mathrm{\Sigma }(\overline{\dot{\gamma }},\overline{\varphi })`$ measured for several values of $`\overline{\varphi }`$ actually contains the same information as the phase diagram: reconstruction of the latter from the former is described in Fig. 10 A full set of flow curves could therefore be used to check measurements of concentration differences.

In this work, we take the opposite approach for convenience, and reconstruct the steady-state flow curves from the tie lines of the phase diagram. The results are shown in Fig. 6(d). The inset shows the same data on a log-log plot, to enable comparison with Ref.berret94a  in which the coexistence plateau in a log-log representation was a reasonably straight line (over the shear-rate range investigated) with slope $`0.3`$. Note that the results shown in Fig. 6(d) are in units of $`G(\varphi =0.11)`$ and $`\tau (\varphi =0.11)`$. In Ref. BRP94 , Berret replotted the flow curves in units of $`G(\overline{\varphi })`$ and $`\tau (\overline{\varphi })`$, finding scaling collapse of the family $`\mathrm{\Sigma }(\overline{\dot{\gamma }},\overline{\varphi })/G(\overline{\varphi })`$ vs. $`\overline{\dot{\gamma }}\tau (\overline{\varphi })`$ in the low shear regime $`\overline{\dot{\gamma }}0`$. We do not find this scaling collapse (Fig. 11) because we have used an artificially large high-shear Newtonian contribution $`\eta \dot{\gamma }`$ for numerical convenience (recall Sec. III.5): the overall zero shear viscosity, $`G(\varphi )\tau (\varphi )+\eta (\varphi )`$ therefore does not scale as $`G(\varphi )\tau (\varphi )`$, even approximately.

To check the reconstruction of flow curves from the phase diagram, we also explicitly calculated the flow curve at a single $`\overline{\varphi }=0.11`$. To do this, we first performed a shear startup at a given $`\overline{\dot{\gamma }}`$ in the unstable region. We then (without reinitialising the system) decreased $`\overline{\dot{\gamma }}`$ in steps to the edge of the coexistence regime, ensuring that a steady state was reached before measuring the total stress. We then reinitialised the system and repeated the entire procedure, but now with increasing $`\overline{\dot{\gamma }}`$jumps. The results are shown in 8 for two different values of $`l`$. The slight discrepancy between the directly measured flow curve “plateaus” (i.e. the inhomogeneouse part of the flow curve) and those reconstructed from the tie lines is due to the finite size of the interface $`\delta `$ relative to the cell $`L`$, and so is smaller for the smaller value of $`\delta /L`$. The construction described in Fig. 10 implicitly assumes that $`\delta /L=0`$.

As expected for this value of $`\alpha `$ (which gives a large concentration difference between the bands; Fig. 6(a)), the steady state flow curve “plateau” slopes strongly upwards in $`\overline{\dot{\gamma }}`$. In Fig. 9 we compare the (directly measured) macroscopic flow curve for the three levels of concentration coupling shown in Fig. 7: as expected, the slope of the flow curve tends to zero with the degree of concentration coupling $`\alpha `$.

The upturn in the measured flow curve at the edge of the coexistence plateau (apparent at the lower binodal for $`\alpha =10^2`$ in Fig. 9) results again from the finite value of $`\delta /L`$: the interface bumps into the edge of the rheometer when one of the bands gets very narrow. We expect this (steady-state) effect to be much less pronounced in experimental systems, since realistic interfaces are much smaller than those used in our numerical study. Only near a critical point, where the interface becomes very broad (for fixed $`l`$), would we expect to see a true steady-state bump at the edge of the plateau. Nonetheless, pronounced bumps are often apparent in data obtained via upward strain-rate sweeps. However in most cases this is likely to be a metastable effect, so that the bump could be eliminated (or at least reduced) by reducing the rate of the sweep grand97 .

As noted in Sec. I, in a curved Couette geometry the “plateau” (B’F’ of Fig. 1) in the flow curve will slope upwards due to the inhomogeneity of the stress field, even without concentration coupling. It should be noted that all calculations in this paper are for a planar shear geometry, and the slope of our flow curves in the coexistence regime results solely from concentration coupling. In fact, the slope in Fig. 6(d) is far greater than one would typically expect from curvature effects: for a Couette cell with radius $`R`$ and gap $`\delta R`$, the stress measured at the inner Couette wall would change by $`\delta \mathrm{\Sigma }/\mathrm{\Sigma }=2\delta R/R`$ over the coexistence regime, and so too would the relative change in torque through the coexistence “plateau”. The slope of Fig. 6(d) would therefore require an atypically large curvature of $`\delta R/R0.5`$.

#### VI.1.3 Interfacial profiles; divergence of interface width at the critical point

We now turn to the interfacial profiles and widths. A full steady state banded profile for $`\alpha =10^2`$ (corresponding to the rightmost/uppermost tie line in fig. 6a,b,c) is shown by the thick lines in Fig. 12. As required, the interface is smooth on the scale of the mesh, but sharp on the scale of the gap size, i.e. $`L/N_y\delta L1`$ where $`\delta `$ is the width of the interface. Note that the shear rate is negative across the gap since we have chosen to move the wall at $`y=0`$; accordingly we have plotted $`W_{xy}`$, since $`W_{xy}`$ is antisymmetric in shear rate. $`W_{wy}`$ is rather small in the high shear band, as expected from the underlying constitutive non-monotonicity. Meanwhile $`W_{xx}`$ is very large, while $`W_{yy}0.5`$ (recall that $`\underset{¯}{\underset{¯}{W}}`$ measures deformation relative to the unit tensor $`\underset{¯}{\underset{¯}{\delta }}`$): this corresponds to the micelles being highly stretched along the flow direction and is consistent with the experimental observation that the first normal stress difference progressively increases throughout the banding regime rehage91 . The concentration is lower in the high shear band, where $`W_{yy}`$ is smaller (more negative): this is a direct result of the tendency of mielles to move up gradients in $`W_{yy}`$, as determined by Eqn. 8 above.

In fact the interface width, $`\delta `$, is slightly different for each order parameter: we define it to be the distance between the two points where the change in that order parameter between the two homogeneous phases is $`25\%`$ and $`75\%`$ complete. For a fixed value of $`l`$ (which sets the overall scale of the interface width), $`\delta `$ diverges at the critical point (for each order parameter). In tracking $`\overline{\varphi }`$ down towards the critical point, therefore, we continually rescaled $`l`$ to ensure that the interface width remained approximately constant. In each case, we measured $`\delta /l`$, for each of $`W_{xy}`$, $`W_{xx}`$, $`W_{yy}`$ and $`\varphi `$: see Fig. 13. According to mean field theory, the divergence should be of the form $`\delta /l(\mathrm{\Sigma }\mathrm{\Sigma }_\mathrm{c})^{1/2}`$. The power $`1/2`$ is accordingly shown in Fig. 13 as a guide for the eye.

### VI.2 Interfacial terms in both the viscoelastic constitutive equation, and in the concentration equation: $`l0`$, $`\xi 0`$.

We now study the effect of including interfacial gradient terms in the concentration equation 8 (so that now $`\xi 0`$) as well as in the viscoelastic equation 9, $`l0`$. Hence, while in the previous section we considered $`rl/\xi =\mathrm{}`$, then, we now consider $`r=O(1)`$. In Fig. 14(a), we give the phase diagram for $`r=0.4`$. Comparing it with our results for $`r=\mathrm{}`$ (also shown in Fig. 14(a)), we see that the slopes of the tie lines and the overall binodal both depend quantitatively on $`r`$. \[The difference between the results for $`r=\mathrm{}`$ and $`r=0.4`$ is far greater than any “error” associated with the fact that we are not quite in the limit $`\mathrm{\Delta }t0`$, $`lN_y\mathrm{}`$, $`\xi N_y\mathrm{}`$, $`l/L0`$ and $`\xi /L0`$.\] This provides a concrete example of the fact that shear-banding coexistence is determined by, and non-universal with respect to, the interfacial terms lu99 . As noted above, this contrasts sharply with the equilibrium case, in which the equations of motion are integrable and so the phase diagram is independent of the interfacial terms. Although conceptually important, this dependence is in practice rather weak: the overall features of the phase diagram are unchanged. The critical point is unaffected.

In Fig. 14(b) we show the corresponding macroscopic flow curves, reconstructed using the tie lines of Fig. 14(a). Because we only calculated a few tie lines in this case, the recontruction is rather sparse. Nonetheless, the slight difference between $`r=\mathrm{}`$ and $`r=0.4`$ is apparent.

In Fig. 12, we compare a full banded profiles for $`r=\mathrm{}`$ and $`r=0.4`$. The slight dependence on $`r`$ is again apparent.

### VI.3 Interfacial gradient terms only in the concentration equation: $`l=0`$, $`\xi 0`$.

Finally we set the interfacial length $`l`$ in the constitutive equation equal to zero. The constitutive equation is now local, and the only source of spatial gradients is the equilibrium correlation length for concentration fluctuations (Eqns. 8 and 3): $`rl/\xi =0`$. In the absence of concentration coupling, $`\alpha =0`$, it is known that there is no uniquely selected, smoothly shear banded state when $`l=0`$ lu99 . Here we investigate whether a smoothly banded state is selected for $`\alpha 0`$, by virtue of the interfacial terms in the concentration equation.

Our numerics only gave a smoothly banded profile for stresses near the critical point, even for the largest accessible values of $`\xi `$ and $`N_y`$. The profiles shown from left to right in Fig. 15(b) are progressively further above the critical point. The tie lines correponding to the smooth profiles near the critical point are shown in Fig. 15(a), alongside the correponding results at $`r=0.4`$ for comparison. Consistent with the discussion of non-universality in the previous section, the phase diagram for $`r=0.0`$ is slightly different from that for $`r=0.4`$ (and is different again from the case $`r=\mathrm{}`$; not shown).

For the spiky profiles, further from the critical point, the binodal of the associated tie lines is irregular (not shown in Fig. 15(a)), suggesting that the steady state is not uniquely selected. In view of this, a natural question is whether selection could occur in principle (but is inaccessible with any realistic mesh due to the pronounced non-monotonicity in $`W_{xy}(y)`$), or whether selection cannot occur, even in principle. In Fig. 15(c) we show that the steady state depends on the initial condition; so state selection appears to be lost when $`l=0`$. This numerical observation is backed up by the following analytical argument.

In steady state, the system must obey:

* The force-balance equation,

$$S(\dot{\gamma },\varphi )G(\varphi )W_{xy}[\dot{\gamma }\tau (\varphi )]+\overline{\eta }(\varphi )\dot{\gamma }=\mathrm{\Sigma }=\mathrm{const}.$$
(17)
* The (now local) constitutive equation, equation (9),

$$W_{\alpha \beta }=W_{\alpha \beta }[\dot{\gamma }\tau (\varphi )]\text{for }\alpha \beta =xx,xy,yy.$$
(18)
* The steady-state of equation 8. For the purposes of this analytical argument we use a simplified version of this equation, which we believe still captures the essential physics:

$$0=_y^2\left\{\mathrm{\Lambda }_y^2\varphi \right\}$$
(19)
with

$$\mathrm{\Lambda }=f^{}(\varphi )\frac{G(\varphi )W_{yy}[\dot{\gamma }\tau (\varphi )]}{\varphi }.$$
(20)
Integrating Eqn. 19 twice, and using the boundary conditions $`_y\varphi =0,_y^3\varphi =0`$ for $`y=0,L`$, we obtain

$$\mu =\mathrm{const}.=\mathrm{\Lambda }_y^2\varphi ,$$
(21)
where $`\mu `$ is an integration constant.

We now show that a solution satisfying Eqns. 17 and 18 cannot in general simultaneously satisfy Eqn. 21.

Consider firstly Eqns. 17 and 18. Sustituting $`W_{xy}`$ from Eqn. 18 into Eqn. 17, we obtain an expression for $`S(\dot{\gamma },\varphi )`$: this is just the family of homogeneous constitutive curves, as plotted in Fig. 4 above. Because the constitutive equation is local, the solution at all points across the rheometer cell must lie on one of these intrinsic constitutive curves. Indeed, as the shear rate changes across the interface, the system must pass through constitutive curves of differing concentrations to maintain a uniform stress $`\mathrm{\Sigma }`$. In other words, a relation $`\varphi =\varphi (\dot{\gamma },\mathrm{\Sigma })`$ must be obeyed. The family of these curves is shown as dotted lines in Fig. 16a. For the range of stresses at which $`\varphi (\dot{\gamma },\mathrm{\Sigma })`$ is non-monotonic, $`\varphi `$ must have the form shown in Fig. 17b in which the derivative $`_y^2\varphi `$ changes sign three times across the interface, as in Fig. 17c. \[Actually, the forms of Fig. 17b,c assume that the profile in $`\dot{\gamma }`$ increases monotonically through the interface (Fig. 17a). However this monotonicity will emerge self consistently from our argument below.\]

However we know from Eqn. 21 that $`_y^2\varphi =\mathrm{\Lambda }\mu `$. $`\mathrm{\Lambda }`$ is plotted in Fig. 16b,c using Eqn. 20 together with the constraint $`\varphi =\varphi (\dot{\gamma },\mathrm{\Sigma })`$ (imposed from Eqns. 17 and 18, as discussed above). From this plot we see that, for any $`\mu `$, a solution that starts and ends in homogeneous phases (for which $`_y^2\varphi =\mathrm{\Lambda }\mu =0`$111While this boundary condition is intuitively clear (the homogeneous phases are in general large compared with the interface), in our numerics we only actually imposed the boundary condtion $`_y\varphi =0`$. We therefore did not, a priori, “overspecify” the differential equation 21. However we shall show below that $`_y^2\varphi `$ must automatically equal zero at each interface. can only involve at most one sign change of $`_y^2\varphi `$ between the boundaries. This inconsistency with Fig. 17c means that a steady banded solution cannot exist for these stress values for which $`\varphi (\dot{\gamma },\mathrm{\Sigma })`$ is non-monotonic. To summarize: for stresses far enough above the critical point that $`\varphi (\dot{\gamma },\mathrm{\Sigma })`$ is non-monotonic, a steady state solution cannot simultaneously satisfy Eqns. 17 and 18 (which imply three sign changes of $`_y^2\varphi `$) at the same time as Eqn. 21 (which only allows one sign change). Therefore there a steady, smoothly banded profile cannot exist for such stresses.

This argument is consistent with the sharp numerical profiles of Fig. 15(b), which are replotted in Fig. 16a (solid lines): each solution should have followed a local (dotted) curve $`\varphi (\dot{\gamma })`$, but instead has jumped across the region in which this curve is non-monotonic. (The stresses used to generate the homogeneous solutions $`\varphi (\dot{\gamma },\mathrm{\Sigma })`$ in Fig. 16a were slightly different from those of the numerical profiles: however the trend is still clear.) Although Eqn. (19) is highly oversimplified, we believe that the failure to negotiate the interface due to the conflict described above is the reason for non-selection in the full, numerically solved model.

We return finally to justify our assumption that the shear rate must increase monotonically through the interface, and to discuss in more detail the nature of the banded solution when it can exist (i.e. for stresses near crticial point where $`\varphi (\dot{\gamma })`$ is monotonic). Multiplying Eqn. 21 across by $`d\varphi /dy`$, integrating on $`\varphi `$, and imposing $`_y\varphi =0`$ at each boundary, we find, for the simplified model of Eqn. 19,

$$_{\varphi _l}^{\varphi _r}𝑑\varphi \left[\mu \mathrm{\Lambda }\right]=0,$$
(22)
which is an “equal areas” construction. ($`\varphi _l`$ and $`\varphi _r`$ denote the boundary values at $`y=0,L`$.)

If, in addition, we were to impose that $`_y^2\varphi =0`$ at each boundary, then the construction must automatically be as shown in Fig. 18a. However we did not actually impose this condition in our numerics, so the construction of Fig. 18b is also possible. This in fact corresponds to a finite system, where the true homogeneous state $`_y^n\varphi =0n`$ is not quite reached at the boundaries. Any other equal areas construction (Fig. 18c) is not possible, for the following reason. Consider starting at point C with $`_y\varphi =0`$ (which we do impose at the boundary in our numerics). Eqn. 21 then tells us that $`_y^2\varphi <0`$ at this point, so the function $`\varphi (y)`$ must curve downwards from its starting point of zero slope. Therefore $`\varphi `$ locally decreases, and the system moves to point $`C^{}`$. Repeating this argument, we find that the system can never cross to the point $`D`$. By similar reasoning, the shear rate must rise monotonically through the interface since any initial fall (from the side of the low shear band) would be similarly unstable to point $`C`$ in Fig. 18c above.

Of course the concentration equation (19) is highly oversimplified. For instance, a more realistic model (such as the one of Eqn. 8) would have $`\varphi `$ dependent prefactors to the $`_y^2\varphi `$ term. The equal areas result of Eqn. (22) is therefore specific to our oversimplified Eqn. (19), and does not hold in general. Nonetheless we believe that Eqn. 19 correctly predicts the absence of a uniquely banded solution for stresses far above the critical point, via the basic conflict between the number of sign changes of $`_y^2\varphi `$ across the interface, described above.

## VII Conclusions

In this paper, we have studied the role of concentration coupling in the shear banding of complex fluids using the two-fluid, non-local Johnson-Segalman model. We have calculated phase diagrams for different degrees of coupling between concentration and mechanical degrees of freedom (molecular strain), and found a phase diagram qualitatively consistent with experiments on micellar solutions at dilutions well below the equilibrium isotropic-to-nematic transition berret94b . Specific points to note are as follows.

1. The coexistence plateau in the steady-state flow curve slopes upward with shear rate, because of the concentration difference between the coexisting bands. The overall plateau height and width decrease with average concentration, terminating in a non-equilibrium critical point. CPCl/NaSal in brine berret94b  shows the same trend.
2. Of the two coexisting bands, the high shear band has a smaller concentration due to the fact that concentration tends to move up gradients in the normal micellar strain component $`W_{yy}`$ (where $`y`$ is the flow-gradient direction). ($`\underset{¯}{\underset{¯}{W}}`$ describes deformation relative to the unit tensor $`\underset{¯}{\underset{¯}{\delta }}`$, and $`W_{yy}`$ is more negative in the high-shear phase than in the lower shear phase.) Tie lines of the phase diagram in the $`\dot{\gamma },\varphi `$ plane therefore have negative slope.
3. The concentration gap is smaller for smaller values of concentration-coupling $`\alpha G^{}(\varphi )/f^{\prime \prime }(\varphi )`$, and tends to zero in the limit $`\alpha 0`$. Accordingly, the coexistence region of the steady-state flow curve becomes flat in this limit.
4. We have described the way in which the flow phase diagram can be reconstructed from the family of flow curves $`\mathrm{\Sigma }(\overline{\dot{\gamma }},\overline{\varphi })`$, measured for several average concentrations $`\overline{\varphi }`$ (Fig. 10).
5. The phase diagram and flow curves depend slightly on the relative size of the interfacial term in the viscoelastic constitutive equation to that in the equation that specifies the concentration dynamics. This is a concrete demonstration of how stress selection and the coexistence conditions of driven systems depend on the nature of the interface, in contrast to equilibrium coexistence.
6. We find no unique state selection when there are no gradient terms in the viscoelastic constitutive equation, except for stresses that are close to the critical point. This implies that, for a model to reproduce a uniquely selected stress, it is not enough to simply have gradient terms only in, for example, the concentration dynamics. The dynamical equations of motion for each degree of freedom must possess inhomogeneous terms to attain selection in all situations. Conversely, in situations where such terms are physically absent, one can expect, under certain conditions, no selection and hence a range of control parameters (shear stress or strain rate) for which the steady states are intrinsically history-dependent.
7. The interface width diverges at the critical point as a power law $`(\mathrm{\Sigma }\mathrm{\Sigma }_\mathrm{c})^n`$ with $`n0.5`$, although $`n`$ differs slightly across the different order parameters.

Although our d-JS-$`\varphi `$ model is highly oversimplified, we believe that it contains the basic ingredients required for a first description of wormlike micellar surfactant solutions at concentrations well below the isotropic-nematic (I-N) transition. In particular, it incorporates the minimal set of realistic degrees of freedom (tensorial order parameter for the micellar strain together with concentration), and unifies a non-monotonic flow curve with the Helfand-Fredrickson coupling between concentration and flow. Similar techniques could be applied to more involved Cates non-linear theory for wormlike micelles cates90 ; SpenCate94 .

We recall a previous calculations by Olmsted et al. was aimed at systems of rigid rods near the I-N transition olmsted99c . In future work we hope to unify these two approaches into a description of wormlike micelles that is valid over the entire concentration range. This should provide a first step towards understanding the crossover regime in the data of Fig. 3, in which the coexistence plateau stress is a non monotonic function of the micellar concentration.

###### Acknowledgements.

We thank J-F Berret and Paul Callaghan for useful discussions, and EPSRC GR/N11735 for financial support.

## Appendix A  d-JS-$`\varphi `$ Equations in Cartesian Coordinates

In this appendix, we give the components of the d-JS-$`\varphi `$ model’s equations for planar shear flow along the $`x`$ direction, allowing gradients only in the flow-gradient direction, $`y`$, as described in Sec. III.3, above. The $`x`$ component of force-balance is (in the zero-Reynolds limit considered in this paper)

$$0=_y\left[G\left(\varphi \right)W_{xy}\right]+\eta _\mathrm{m}_y\left[\varphi _yv_{mx}\right]+\eta _\mathrm{s}_y\left[\left(1\varphi \right)_yv_{sx}\right].$$
(23)
The $`y`$ component of force-balance is fixed by incompressibility, $`\underset{¯}{}.\underset{¯}{v}=0`$, along with the boundary condition $`v_y=0`$:

$$0=\varphi v_{my}+(1\varphi )v_{sy}.$$
(24)
The relative velocity between the micelles and solvent (again ignoring inertial terms) is

$`v_{my}v_{sy}`$ $`={\displaystyle \frac{\varphi (1\varphi )}{\zeta }}\left\{{\displaystyle \frac{1}{\varphi }}_y\left[G\left(\varphi \right)W_{yy}\right]+2{\displaystyle \frac{1}{\varphi }}\eta _\mathrm{m}_y\left[\varphi _yv_{my}\right]2{\displaystyle \frac{1}{1\varphi }}\eta _\mathrm{s}_y\left[\left(1\varphi \right)_yv_{sy}\right]_y𝔉\right\}`$ (25)
$`v_{mx}v_{sx}`$ $`={\displaystyle \frac{\varphi (1\varphi )}{\zeta }}\left\{{\displaystyle \frac{1}{\varphi }}_y\left[G\left(\varphi \right)W_{xy}\right]+{\displaystyle \frac{1}{\varphi }}\eta _\mathrm{m}_y\left[\varphi _yv_{mx}\right]{\displaystyle \frac{1}{1\varphi }}\eta _\mathrm{s}_y\left[\left(1\varphi \right)_yv_{sx}\right]\right\}`$ (26)
$`𝔉`$ $`=f^{}\left(\varphi \right)g_y^2\varphi +\frac{1}{2}G^{}\left(\varphi \right)\left[W_{yy}+W_{xx}\mathrm{ln}\left(W_{yy}W_{xx}+W_{yy}+W_{xx}+1W_{xy}^{}{}_{}{}^{2}\right)\right].`$ (27)
The evolution of the micellar strain tensor is given by

$`_tW_{xy}+v_{my}_yW_{xy}`$ $`=\frac{1}{2}(a1)W_{xx}_yv_{mx}+\frac{1}{2}(1+a)W_{yy}_yv_{mx}+aW_{xy}_yv_{my}+_yv_{mx}{\displaystyle \frac{W_{xy}}{\tau \left(\varphi \right)}}+{\displaystyle \frac{l\left(\varphi \right)^2_y^2W_{xy}}{\tau \left(\varphi \right)}},`$ (28)
$`_tW_{yy}+v_{my}_yW_{yy}`$ $`=(a1)W_{xy}_yv_{mx}+2aW_{yy}_yv_{my}+2_yv_{my}{\displaystyle \frac{W_{yy}}{\tau \left(\varphi \right)}}+{\displaystyle \frac{l\left(\varphi \right)^2_y^2W_{yy}}{\tau \left(\varphi \right)}},`$ (29)
$`_tW_{xx}+v_{my}_yW_{xx}`$ $`=(1+a)W_{xy}_yv_{mx}{\displaystyle \frac{W_{xx}}{\tau \left(\varphi \right)}}+{\displaystyle \frac{l\left(\varphi \right)^2_y^2W_{xx}}{\tau \left(\varphi \right)}}.`$ (30)
Finally, the concentration dynamics are

$`_t\varphi `$ $`=`$ $`_y\left\{{\displaystyle \frac{\varphi ^2\left(1\varphi \right)^2}{\varphi \zeta }}\left\{_y\left[G\left(\varphi \right)W_{yy}\right]+2\eta _\mathrm{m}_y\left[\varphi _yv_{my}\right]{\displaystyle \frac{2\varphi }{1\varphi }}\eta _\mathrm{s}_y\left[\left(1\varphi \right)_yv_{sy}\right]_y𝔉\right\}\right\}.`$ (31)
## Appendix B Stationary homogeneous solutions of the d-JS-$`\varphi `$ model

In planar shear, the stationary homogeneous solutions to Eqns. (7-9) for given $`\dot{\gamma }`$ and $`\varphi `$ are $`\underset{¯}{v}_{\mathrm{rel}}\underset{¯}{v}_\mathrm{m}\underset{¯}{v}_\mathrm{s}=0`$ and

$`W_{xy}`$ $`={\displaystyle \frac{\dot{\gamma }\tau (\varphi )}{1+b\dot{\gamma }^2\tau ^2(\varphi )}},`$ (32a)
$`W_{yy}`$ $`={\displaystyle \frac{a1}{1+a}}W_{xx}={\displaystyle \frac{1}{(1+a)}}{\displaystyle \frac{b\dot{\gamma }^2}{1+b\dot{\gamma }^2}}`$ (32b)
$`W_{zz}`$ $`=W_{xz}=W_{yz}=0,`$ (32c)
where $`b=1a^2`$. The steady state shear stress is given by

$$\mathrm{\Sigma }_{xy}=G(\overline{\varphi })W_{xy}+\varphi \eta _\mathrm{m}+(1\varphi )\eta _\mathrm{s}\dot{\gamma }=\text{constant}.$$
(33)


text:
# Location of H2O maser in the double-nuclei system of NGC 6240

## 1 Introduction

### 1.1 Background of extragalactic H2O masers

A number of attempts to find more extragalactic H2O masers (616–523) have been made in recent times, stimulated by the discovery of the highly Doppler-shifted maser emission symmetrically straddling the systemic velocity of the galaxy in NGC 4258 (Nakai et al. (1993)). VLBI observations of the maser emission in NGC 4258 revealed the presence of a sub-parsec-scale H2O maser disk rotating around a central massive object (Miyoshi et al. (1995)). After the discovery of the masing torus surrounding a ’black hole’ in NGC 4258, very few other such sources suitable for studying the central sub-parsecs of Active Galactic Nuclei (AGN) have been discovered, but astronomers need more candidates to probe the full range of the kinematics and dynamical structures of sub-parsec-scale circumnuclear regions in AGN. This kind of study is possible through the technique of radio interferometry. The extragalactic H2O masers observed and studied to date can be categorized into two types in terms of the origin of the emission. Most of $`lowluminosity`$ masers (L$`_{\mathrm{H}_2\mathrm{O}}`$ $`<`$ $``$ 1 L) have been observed outside galactic nuclear regions, and are associated with star forming regions and (compact-)HII regions, suggesting that the masers have no direct connection to AGN-activity. $`Highluminosity`$ masers (L$`_{\mathrm{H}_2\mathrm{O}}`$ $`>>`$ 10 L) are apt to be located in or nearby an active nucleus, or to be associated with jet-activity. Such a water maser which associates with the AGN-activity is a $`nuclearmaser`$. In either case they trace the nuclear activity within the central few parsecs of the AGN. There are also a number of masers with luminosity on the order of 10 L, for instance, those in Mrk 1 and NGC 5506 (Braatz et al. (1994)).

The origins of H2O masers with luminosities in the range 10 L $`<`$ L$`_{\mathrm{H}_2\mathrm{O}}`$ $`<`$ 100 L have not been well studied in comparison with those of higher luminosity because the H2O emission is generally not sufficiently bright for radio interferometric observations with reasonable sensitivities. Intensity variabilities on time scales of weeks to months are commonly observed in these masers, which has posed difficulties in determining the optimum timing of interferometric experiments. Currently, about 25–30 extragalactic H2O masers, excluding the low luminosity masers, are known. However, only about ten of them have been so far imaged at milliarcsecond resolution using VLBI (Moran et al. (1999); Greenhill (2000); Hagiwara et al. 2001a ; Ishihara et al. (2001)). Radio interferometric studies of this sort of H2O maser are quite challenging, but are necessary to investigate the overall properties of the H2O masers in active galaxies regardless of their maser intensity. The apparent luminosity of the maser is not always relevant to its intrinsic properties.
In a previous paper, we reported on the single-dish detection of the H2O maser in NGC 6240 (Hagiwara et al. 2002a ), and here present the results of follow-up observations at subarcsecond angular resolution aimed at investigating the spatial distribution of the maser.

### 1.2 NGC 6240

NGC 6240 (D=97 Mpc, assuming that H0 is 75 km s-1Mpc-1; 0.1”=48.5 pc) is a well-studied galaxy at various wavelengths. It is a prototypical ultraluminous FIR galaxy with a complex optical structure, hosting a LINER nucleus (Heckman et al. (1987)). The huge FIR luminosity ($`>10^{11}`$ L) is considered to originate from re-emission of ultraviolet radiation processed by a dust shell surrounding a hot core, suggesting ongoing star formation in the galaxy (e.g., Sanders et al. (1988)). X-ray satellite observations detected 6.4 keV Fe K$`\alpha `$ line emission in the galaxy, showing a definite presence of an AGN together with a large X-ray luminosity of $``$ 1043 \- 1044 erg s-1 in the 2 - 10 keV bands (e.g., Iwasawa & Comastri (1998)). At sub-arcsecond resolution a number of hard X-ray blobs were imaged by Chandra, and some of them extend several kiloparsec from the center of the galaxy (Lira et al. (2002)). 5 GHz high-resolution continuum MERLIN maps made by Beswick et al. (2001) showed two compact radio sources (double nuclei) with non-flat spectral indices ($``$ 0.7 – 0.8 at 4.9 –15 GHz (Colbert et al. (1994), Carral et al. (1990))), and they are separated by 1”.5 at the highest resolution so far attained ($``$ 0.1 arcsecond) in the nuclear region. The velocity fields of the broad HI absorption towards these components are seen running along the position angle joining the two radio nuclei. Beswick et al. (2001) found that the radio continuum flux in the galaxy is mostly dominated by a starburst component, with only a small portion arising in a low luminosity AGN. Tacconi et al. (1999) and Bryant et al. (1999) found a significant peak of CO molecular gas between the double radio nuclei. Their speculation is that the CO gas settling down between the two stellar nuclei will eventually form a central thin disk. With lower angular resolutions but higher sensitivities than MERLIN, the VLA revealed HI and OH absorption structures, which agree with those of H2, CO, and Hα (Baan & Hofner (2002)). There has been an open discussion on the need to introduce an AGN in NGC 6240 as an energy supplier sustaining the large infrared luminosity, however, no one has yet succeeded in determining the location of the predicted AGN in the double nuclei system at any wavelengths.

### 1.3 H2O maser in NGC 6240

Tentative detections of an H2O maser were made by Henkel et al. (1984) and Braatz et al. (1994), although never confirmed. Hagiwara et al. (2002a) reported on the first clear detection of the water maser in the galaxy in May 2002. They also detected intensity variability on timescales of weeks to months by means of single-dish monitoring of the two distinct components, one of which lies at VLSR= 7565 km s$`^1`$and the other at VLSR= 7609 km s-1. Given the fact that the apparent maser luminosity of the two components is about 40 L, well above that of the typical starforming masers, and the emission is highly redshifted by some 400 km s$`^1`$with respect to the adopted systemic velocity of VLSR= 7131 km s$`^1`$(Hagiwara et al. 2002a ), the maser is likely to arise from the AGN activity and can be used to probe the kinematics in the galaxy.

## 2 Observations and data analysis

### 2.1 The 100-m radio telescope in Effelsberg

After the first single-dish detection of the H2O maser emission in the rotational 616 – 523 transition (rest frequency = 22.23508 GHz) towards NGC 6240 in May 2001 (Hagiwara et al. 2002a ), we observed the maser at several epochs until September 2002. All the observations were performed with the Max-Planck-Institut für Radio astronomie (MPIfR) 100-m radio telescope in Effelsberg using a two channel K-band HEMT receiver with dual polarizations and a digital autocorrelator (AK90) as a backend. The backend correlator provided a total bandwidth of 4 $`\times `$ 40 MHz for each polarization (velocity coverage $``$ 2100 km s-1) with a spectral resolution of 78 kHz (1.1 km s-1). The beam size of the telescope is so small ($``$ 40″(HPBW) at 22 GHz) that the frequent calibration of telescope pointing was critical. It was determined by observing a nearby quasar 1655+077 every 1.5 – 2 hrs during the observations, yielding the pointing accuracy of 4″– 6″. More details of the observing parameters were remarked in Hagiwara et al. (2002a). Data analysis consisted of averaging the spectra and subtracting a linear baseline in general, but occasionally a non-linear baseline.

### 2.2 Very Large Array

NRAO111The National Radio Astronomy Observatory (NRAO) is operated by Associated Universities, Inc., under a cooperative agreement with the National Science Foundation. Very Large Array (VLA) observations were carried out in the B configuration on 16 July 2002. Prior to that epoch, VLA or VLBA observations had been attempted several times in 2001 with no clear detections; a list of observations is displayed in Table 1. The current observations (this paper) employed a single IF band of 12.5 MHz divided into 64 channels of width 2.64 km s$`^1`$ each. The band was centered at an LSR velocity (VLSR) of 7585 km s-1, which is about 450 km s$`^1`$redshifted from the systemic velocity of the galaxy. After removal of several channels at the band edges, the usable velocity range was 7510 – 7665 km s-1. To detect maser lines with weak flux densities of $`<`$ 50 mJy, we conducted phase-referencing observations using a nearby quasar, 16582+074 as a phase-reference source. The phase-referencing observations were executed in a sequence of 3 min scans with cycling interval of 140 sec for NGC 6240 and 40 sec for 16582+074. The total observing time was three hours, and the resultant synthesized beam has a FWHM size of $`0\text{.}{}_{}{}^{\prime \prime }29\times 0\text{.}{}_{}{}^{\prime \prime }26`$ with uniform weighting and a position angle of $``$34. The amplitude and bandpass calibrations were made by observing 3C 286, assuming $`S(3\mathrm{C}\mathrm{\hspace{0.17em}286})=`$ 2.59 Jy at 21.7 GHz. The data were calibrated and mapped in the standard way using AIPS. The continuum emission was subtracted from the spectral line visibility cube using UVLSF, and then velocity channel maps of the H2O emission were made for each velocity channel. The resultant 64 channel maps each had a typical rms noise level of 1.3 mJy beam-1. The continuum emission was imaged with an rms noise of 0.26 mJy beam-1.

## 3 Results

### 3.1 Effelsberg

Fig. 1 shows H2O maser spectra obtained from the MPIfR 100-m telescope for 4 epochs between 10 May 2001 and 14 September 2002. Each spectrum contains the two distinct redshifted components at VLSR= 7565 km s$`^1`$and 7609 km s$`^1`$and other weaker components that require confirmation. We searched for maser emission primary from VLSR$``$ 6800 – 7800 km s-1. Occasionally, the highly blueshifted velocities ranging to VLSR= 5900 km s$`^1`$were searched (Table 1). No emission was distinctly detected except for the two components with 3 $`\sigma `$ limits of 15 – 60 mJy. The spectra of the two epochs in 2001 have better signal-to-noise ratios than those in 2002 largely because the integration times are much longer. Given the fact that the velocities of the two components were not changed over 8 months, an upper limit of drift rate of each line is 1.3 km s-1 year-1. The flux density of the VLSR= 7565 km s$`^1`$component was measured to be $``$ 50 mJy in March 2002, the highest value recorded during our single-dish monitoring. The maser was not detected in September 2002 with an rms of $``$ 20 mJy. Neither component was detected by the VLA snapshot made in early Oct 2001 with a 3 $`\sigma `$ level $``$ 15 mJy per 2.6 km s$`^1`$spectral channel. This implies that the maser once faded after July 17 2001.

### 3.2 VLA

A naturally weighted VLA 22 GHz continuum map is shown in Fig. 2 at 0.3” resolution; it displays double-peaked radio emission in the nuclear region in radio robes around each peak. The radio morphology is consistent with the 8.4 GHz and 15 GHz continuum images observed with the VLA in A configuration (Colbert et al. (1994)). The components ’S’ and ’N’ in Fig.2 correspond to N1 and N2 in Fig.2 or Fig. 3 of Colbert et al. (1994). Spectra of the 22 GHz maser emission towards S and N are also shown in Fig.2, both of which were determined from the areas containing the continuum peaks.

The maser line spectrum towards S in Fig.2 shows a single line profile with a flux density of about 35 mJy (L$`_{\mathrm{H}_2\mathrm{O}}`$ $``$ 20 L) , and a centroid velocity of VLSR= 7611 km s-1, while the spectrum towards N does not show any emission within a 3$`\sigma `$ level of $``$ 4 mJy per velocity channel (L$`_{\mathrm{H}_2\mathrm{O}}`$ $``$ 2 L), suggesting that N houses no high luminosity maser at this time. Since the velocity resolutions of the 100-m and VLA spectra are different, there is the ambiguity in estimating the centroid velocity of each line in each spectrum. The uncertainties of velocity are not better than the VLA channel spacings of 2.64 km s-1. Therefore, we conclude that the VLA component lying at VLSR= 7611 km s$`^1`$is identical to the 7609 km s$`^1`$component detected by the 100-m and shown in Fig.1. The position of the maser is $`(\alpha ,\delta )_{\mathrm{J2000}}`$ = $`16^\mathrm{h}52^\mathrm{m}58\text{.}^\mathrm{s}88\pm `$ 0.s05, $`+02^{}24^{}03\text{.}\mathrm{}3\pm 0\mathrm{}.1`$. The relative positional error, dominated by statistical noise, between the maser and the 22 GHz continuum peak S is estimated to be $``$ 0.007 arcsec. Accordingly, the maser coincides with the continuum peak within an astrometric uncertainty of $``$ 0.007 arcsec, or $``$ 3 pc. The maser emission remains unresolved at our angular resolution of 0.26”, or 126 pc. Assuming a Gaussian distribution of the H2O emission, the source size deconvolved by the clean beam is 0.12 $`\times `$ 0.08 arcsec (58 pc $`\times `$ 39 pc). H2O line emission was not detected with a 4 $`\sigma `$ limit of 5.2 mJy from any other observed points in the galaxy, including towards N. We reported earlier on the marginal interferometric detection of the VLSR= 7565 km s$`^1`$component obtained by the VLA (C configuration) in Aug 2001, the position of which was significantly displaced from both S and N (Hagiwara et al. 2002b ). The sensitivity of that map is, however, quite poor, and the detection level of the emission was $``$ 3 $`\sigma `$. We will not use that result in the later discussion.

Time variability of the maser intensity is that the two distinct components show changes in flux density from 10 – 50 mJy over weeks to months. There seems to be a weak anti-correlation of intensities between the two components.

## 4 Discussion

### 4.1 Location of the H2O maser in NGC 6240

The most significant result from our high-resolution imaging of NGC 6240 is that the spatially unresolved maser is exactly coincident with the center of the southern nucleus (S) and there is no detectable H2O emission from the northern nucleus (N). There is compelling evidence that luminous H2O maser emission marks the radio continuum nucleus of the galaxy as in NGC 4258 (e.g., Miyoshi et al. (1995)). The detected maser in NGC 6240 may well probe the position of one of the double-nuclei, S. Given the detection of features redshifted by $`<`$ 400 km s$`^1`$from the systemic velocity and the narrowness of each line, the maser might originate from a well-defined region, for example, a tangential point in the receding side of a rotating disk. In addition such extremely high velocity shifts of each line have never been measured in stellar jets or outflows. Alternatively, it is possible that the redshifted maser emission is associated with the jet-activity of AGN. We cannot support the latter at this moment because there has not been detection of any spatially resolved core-jet structure in S and N at the highest resolution of 0.055 arcsec (Beswick et al. (2001)). The interferometric observing results at different radio frequencies leave the possibility that these two radio sources are stellar, i.e. non-galactic, components as they were resolved by milliarcsecond VLBI observations at 1.4 GHz resulting in non-detection of ultra-compact radio sources on scales of 0.5 pc – 5 pc in N and S (unpublished results; Hagiwara et al.), and neither N nor S shows the flat radio spectral indices (Table 2) that are expected for such resolved nuclei (Colbert et al. (1994)). The apparent H2O maser luminosity of 20 L is nevertheless brighter by two orders of magnitude than that of starforming masers such as NGC 253 and M 82; their typical H2O maser luminosity is $``$ 0.1 L or less, and neither of them lies in a nucleus, i.e. the dynamical center of the galaxy (Ho et al. (1987); Baudry & Brouillet (1996)). If S is one of the real nuclei of the galaxy, the interpretation of the maser will be straightforward - a nuclear maser with the emission physically confined within 1 pc of S. It is plausible that the other maser component will be detected towards N. The absence of compact radio emission between N and S is not consistent with the detection of the CO (1-0) between them, but a radio-quiet AGN is not ruled out. The dust shell foreground to the double peak system which is responsible for the strong far-infrared emission may hide a true nucleus in the background, or obscure these nuclei. The dusty environment itself is favorable to maser emission, and the maser could be enhanced through the thick dust lane in our line of sight, amplifying the background continuum source. It is interesting that the dust in the galaxy causes highly variable extinction with AV = 2 – 8 at optical bands (Scoville et al. (2000)).

### 4.2 Physical conditions for the luminous H2O maser

The observation that the presence of extragalactic OH and H2O maser emission is mutually exclusive appears to be due to the different properties of host galaxies: the OH megamaser galaxies show large FIR luminosities and gas-rich properties in the circumnuclear region but do not harbor a distinct active nucleus, while the H2O maser galaxies on the whole contain an AGN which provides the conditions necessary for maser amplification. NGC 6240 however exhibits both OH and H2O molecular gas towards the double nuclei system. OH molecules in the galaxy are not seen to be producing maser emission but exist as absorbing gas that is extended on kiloparsec scales, suggesting that the water masing cloud and OH absorbers do not coexist. In general, do OH megamasers coexist with H2O megamasers under similar physical environments? We note one theoretical study by Neufeld et al. (1994) on the dense circumnuclear gas in a thin gaseous layer, illuminated by X-rays from AGN, which produces the conditions that are necessary for H2O maser excitation. When the temperature of a molecular gas phase heated by nuclear X-ray radiation from a central engine rises sufficiently ($``$ 400 – 1000 K), a pair of chemical reactions of O $`+`$ H2 = OH $`+`$ H and OH $`+`$ H2 = H2O $`+`$ H will occur (e.g., Elitzur (1992)). The density of OH, hence, will be relatively smaller, and the abundance ratio of water to hydrogen will be greater than 10-4 as the reaction progresses. Note that the temperature required to enhance the reactions is $`>`$ 400 K, but under such conditions the OH molecule is less likely to emit as a maser. With such H2O abundance, the apparent luminosity of the saturated H2O maser can range from about 30 - 300 L pc-2 over a wide range of parameters of X-ray fluxes, gas pressures, densities, gas temperatures, and effective scales of X-ray heating from the surface of a slab of material, all of which are derived from reasonable observational values. Even with the introduction of a torus model instead of the slab, the luminosity will be of similar value (Wallin & Watson (1997)). Such X-ray heating models producing luminous H2O masers will not give rise to OH masers in the warm gas.
The high-luminosity maser must be located in a smaller area than the VLA synthesized beam (e.g., Claussen & Lo (1986)). The location of the unresolved maser is well defined within a radius of the beam-deconvolved size of $``$ 30 pc. If we adopt the total maser luminosity of 20 L, the luminosity per square parsec is L$`_{\mathrm{H}_2\mathrm{O}}`$ $``$ 0.01 L pc-2. The relative location of the maser in NGC 6240 is measured to be within 3 parsecs of its continuum peak, yielding L$`_{\mathrm{H}_2\mathrm{O}}`$ $``$ 0.7 L pc-2. Recent milliarcsecond VLBI studies have demonstrated that the distribution of H2O megamaser emission is, in general, confined to $`1`$ pc of the center of galaxies (Moran et al. (1999)). With this in mind, the lower limit of the NGC 6240 maser luminosity will be comparable to 20 L pc-2, which agrees with the minimum value of the theoretically estimated luminosity. The H2O maser in NGC 6240 can, therefore, be explained by the AGN-activity.

### 4.3 Gas kinematics in NGC 6240

If the maser is a nuclear maser, it should lie in the dynamical center of the galaxy, where the observed dense molecular gas will eventually settle down between the two radio sources and form a disk (Bryant & Scoville (1999); Tacconi et al. (1999)). However, our data show that the maser lies in one of the nuclei displaced by some 100 pc from the CO emission peak lying between the twin nuclei. Is this the real dynamical center of the galaxy? It is understood that there is no well-identified active nucleus in the galaxy, and the position of the detected 6.4 keV Fe line emission is uncertain (Ikebe et al. (2000)). Beswick et al. (2001) found different centroid velocities of the HI absorption towards N and S; they peaked at VLSR= 7260 km s$`^1`$and 7087 km s-1, respectively. The velocity ranges of the absorption do not overlap with any of the known maser components. The radial HI velocity towards S is blueshifted by about 50 km s$`^1`$from the galaxy’s systemic velocity (VLSR=7131 km s-1), while the velocity towards N is redshifted by about 130 km s$`^1`$to the systemic velocity (Table 2). The nucleus S clearly shows different kinematics from N. This situation is similar to the radial velocity difference between the nuclei of 150 km s$`^1`$in H2 $`\nu `$ = 10 S(1) and 100 km s$`^1`$in the CO (J = 2 $``$ 1) (Tecza et al. (2000)).

If the maser were to exhibit a symmetric velocity distribution, we might assume that the components we detect at VLSR= 7565 km s$`^1`$and 7611 km s-1are red-shifted components and we might therefore search for emission around the ’systemic’ velocity of S (VLSR$``$ 7080 km s-1), and the blueshifted counterparts of the red-shifted emission around VLSR= 6500 km s$`^1`$and 6550 km s-1. However, in the latest Effelsberg observations of the maser carried out on 14 September 2002, no new components were detected within $`\pm `$ 800 km s$`^1`$centered on the systemic velocity to an rms of $``$ 20 mJy. On the other hand, the velocities of CO emission (J = 1 $``$ 0, 2 $``$ 1) span a velocity width of $`\pm `$ 500 km s$`^1`$centered on the systemic velocity, reaching to the redshifted maser velocity. The CO emission that peaks between the twin nuclei seems to trace a circumnuclear region on larger galactic disk scales. The orientation and position of the galactic scale disks such as those imaged by HI, CO or HCN gas have little correlation in general with those of a maser disk on sub-parsec scales (e.g., Hagiwara et al. 2001b ). If the two radio nuclei are highly obscured by the thick dust lane and scattered by the foreground materials in the galaxy, they will be observed as less intense due to absorption and are likely to be resolved with milliarcsecond resolution. Note that the radio spectral properties of N and S are not significantly different (Colbert et al. (1994)), suggesting that both of the nuclei could house AGN.

In conclusion, we propose that one of the nuclei lies in S, and the AGN activity which is more dominant at S than any other locations in the galaxy gives rise to the maser emission. However, the starburst or composite starburst plus AGN components at the nuclei cannot be ruled out, despite the fact that such a luminous ($`>`$ 10 L) water maser has not ever been observed in starburst galaxies (e.g., Claussen & Lo (1986)). It is essential to measure the locations of all visible maser components with VLBI. Our interpretation would have to be reconstructed if any of the other velocity components are found to be located outside the nuclei in future observations.

###### Acknowledgements.

YH thanks R. Beswick and W. Baan for useful comments and discussion. The authors would like to thank the referee, Jim Braatz, for his critical and careful reading of the manuscript.