text:
# Shell Model of Two-dimensional Turbulence in Polymer Solutions

## Abstract

We address the effect of polymer additives on two dimensional turbulence, an issue that was studied recently in experiments and direct numerical simulations. We show that the same simple shell model that reproduced drag reduction in three-dimensional turbulence reproduces all the reported effects in the two-dimensional case. The simplicity of the model offers a straightforward understanding of the all the major effects under consideration.

The effect of soluble polymers on two-dimensional turbulent flows appear quite surprising 02AK ; 03BCM . Instead of reducing the drag as in three-dimensional turbulence 75Vir ; 00SW  (thus increasing the large scale velocity components 03ACBP ; 03BDGP ; 03BCHP ), in two dimensions polymers appear to suppress the large scale velocity components. In addition, the polymers affect the probability distribution functions of the velocity field, changing them from sub-gaussian to super-gaussian, with approximately exponential tails. These phenomena appear generic, and were observed both in experiments in fast flowing soap films 02AK  and in direct numerical simulation of the two dimensional viscoelastic equations 03BCM . In this Letter we argue that these effects can be understood using a simple shell model of viscoelastic turbulence.

The shell model used to describe the effects of polymers on the turbulent velocity field had been derived 03BDGP ; 03BCHP  on the basis of the FENE-P model of visco-elastic flows 87BCAH . It was employed recently to understand successfully drag reduction in homogeneous three dimensional visco-elastic turbulence 03BDGP ; 03BCHP . It reads

$`{\displaystyle \frac{du_n}{dt}}`$ $`=`$ $`{\displaystyle \frac{i}{3}}\mathrm{\Phi }_n(u,u){\displaystyle \frac{i}{3}}{\displaystyle \frac{\nu _p}{\tau }}P(B)\mathrm{\Phi }_n(B,B)\gamma _nu_n+F_n,`$
$`{\displaystyle \frac{dB_n}{dt}}`$ $`=`$ $`{\displaystyle \frac{i}{3}}\mathrm{\Phi }_n(u,B){\displaystyle \frac{i}{3}}\mathrm{\Phi }_n(B,u){\displaystyle \frac{1}{\tau }}P(B)B_n\nu _Bk_n^2B_n,`$
$`P(B)`$ $`=`$ $`{\displaystyle \frac{1}{1_nB_n^{}B_n}}.`$
$`\gamma _n`$ $`=`$ $`\nu k_n^4+\mu k_n^4.`$ (1)
In these equations $`u_n`$ and $`B_n`$ stand for the Fourier amplitudes $`u(k_n)`$ and $`B(k_n)`$ of the two respective fields, but as usual in shell models we take $`n=0,1,2,\mathrm{}`$ and the wavevectors are limited to the set $`k_n=2^n`$. The nonlinear interaction terms take the explicit form

$`\mathrm{\Phi }_n(u,B)=k_n\left[(1b)u_{n+2}B_{n+1}^{}+(2+b)u_{n+1}^{}B_{n+2}\right]`$
$`+k_{n1}\left[(2b+1)u_{n1}^{}B_{n+1}(1b)u_{n+1}B_{n1}^{}\right]`$
$`+k_{n2}\left[(2+b)u_{n1}B_{n2}+(2b+1)u_{n2}B_{n1}\right],`$ (2)
with an obvious simplification for $`\mathrm{\Phi }_n(u,u)`$ and $`\mathrm{\Phi }_n(B,B)`$. Here $`b`$ is a parameter which can be used to distinguish between three-dimensional and two-dimensional behavior (and see below for details). In accordance with the generalized energy of the FENE-P model, the non linear terms in our shell model conserves the total energy:

$$E\frac{1}{2}\underset{n}{}|u_n|^2\frac{1}{2}\frac{\nu _p}{\tau }\mathrm{ln}\left(1\underset{n}{}|B_n|^2\right).$$
(3)
The dissipative term $`\gamma u_n`$ contains both a hyper viscosity at the small scales and a friction term at the large scales, to mimic wall friction in 2-dimensional turbulence. With $`\nu _p=0`$ the first of Eqs. 1 reduces to the well-studied Sabra model of Newtonian turbulence 98LPPPV . As in the FENE-P case we consider $`\nu _p/\nu `$ to be the concentration of the polymer $`c`$. For $`\nu _p0`$ we refer to the model as the SabraP shell model. The forcing in Eqs. (1) is chosen to input a fixed amount of energy per unit time, i.e.

$$F_n=\delta _{n,f}\stackrel{~}{F}/u_n^{},\stackrel{~}{F}=10^3\sqrt{2}\text{ for }f=10,11.$$
(4)
The Sabra model with the parameter $`b`$ in the range $`1b0`$ agrees with the 3-dimensional Navier-Stokes equations in the sense that the energy cascade is normal, from the large scales to the smaller scales. For $`b<b_c=12^{2/3}`$ the situation changes qualitatively, since the energy flux changes direction, going from smaller to larger scales, as in 2-dimensional turbulence 02GLPP . The role of the direct energy flux is taken by an enstrophy-like conserved variable $`H`$ which cascades from larger to smaller scales:

$$H\underset{n}{}\left(\frac{1}{1+b}\right)^n|u_n|^2.$$
(5)
Thus, if the forcing is applied to an intermediate level $`k_f`$, there are two scaling regimes with an inverse energy flux supported on scales $`k_n<k_f`$ and a direct enstrophy-like flux for $`k>k_f`$. The spectrum $`S_2(k_n)u_nu_n^{}`$ has a Kolmogorov exponent$`S_2(k_n)k_n^{2/3}`$ in the inverse cascade regime, whereas in the direct cascade regime the spectrum reads (up to intermittency corrections) 02GLPP

$$S_2(k_n)k_n^{2[1+\mathrm{log}_2(1/1+b)]/3}.$$
(6)
When the parameter $`b`$ is in the range $`1b0`$ the coupling to the $`B`$-field results in drag reduction, as had been discussed in refs. 03BDGP ; 03BCHP . We demonstrate here that for $`b<12^{2/3}`$ the coupling to the $`B`$-field suppresses the large scale velocity components. In Fig. 1 we show the respective spectra of the Sabra and the SabraP models, in which the forcing is identical. Note that we force at shells 10 and 11, and therefore we observe only the direct cascade, in agreement with the experimental and simulational works 02AK ; 03BCM . The inverse cascade regime is destroyed by the large scale friction and is not observable here, as in 02AK ; 03BCM . Eq. (6) predicts a spectral slope of -0.88 for the present value of $`b`$, while the observed slope is -0.96$`\pm 0.05`$; the difference is attributed to the usual intermittency correction. The main point is that the spectrum with the coupling to the $`B`$ field exhibits a strong suppression of the large scales, and the slope reduces to about -0.66$`\pm 0.05`$.

We also show, in the same plot, the spectrum of the $`B_n`$ field, $`|B_n|^2`$. This spectrum has a typical scale $`k_c`$ which is known as the Lumley scale 69Lu ,

$$\sqrt{S_2(k_c)}k_c\tau ^1.$$
(7)
For wavevector smaller than $`k_c`$ the $`B_n`$-spectrum is strongly suppressed, whereas for $`k_n>k_c`$ the spectrum is a power law.

The advantage of the present simple model is that the drag enhancement observed in the spectra lends itself to a straightforward interpretation. The standard newtonian energy balance equation is now changed due to the second term in the generalized energy (3) which contributes to the dissipation a positive definite term of the form $`(\nu _p/\tau ^2)P^2(B)_n|B_n|^2`$. Thus the energy balance equation reads

$$\stackrel{~}{F}=ϵ+(\nu _p/\tau ^2)P^2(B)\underset{n}{}|B_n|^2+\mu \underset{n}{}k_n^4|u_n|^2,$$
(8)
where

$$ϵ=\nu \underset{n}{}k_n^4|u_n|^2$$
(9)
is the dissipation due to hyper viscosity. In 3-dimensional turbulence the dissipation is independent of the viscosity $`\nu `$ (This fact is known in the jargon as “the viscous anomaly”). In the present case where the energy flux reverses direction $`ϵ`$ vanishes in the limit $`\nu 0`$, and for large Reynolds numbers can be neglected 02BCMV . Observing the fact that the last sum in (8) is strongly dominated by the forcing scale $`k_f`$ we can safely estimate

$$|u_f|^2\frac{Ck_f^4}{\mu }\left[\stackrel{~}{F}(\nu _p/\tau ^2)P^2(B)\underset{n}{}|B_n|^2\right],$$
(10)
where $`C`$ is a constant of the order of unity. The first term in the parenthesis represents the Sabra value of $`|u_f|^2`$. The subtraction of the positive definite second term is the explanation of the observed spectra in Fig. 1.

Another phenomenon that was discovered in the experiments and the direct numerical simulations is a significant change in the pdf’s of the velocity fluctuations. We demonstrate that the same phenomenon is recaptured in our SabraP model, showing that it is generic to the interaction with the $`B`$ field. In Fig. 2 we present the probability distribution functions (pdf) of $`u_n`$ for $`n=15`$, which is within the bulk of the direct enstrophy cascade. The qualitative change in the pdf is obvious. We quantify the difference by fitting the tails of the pdf’s to stretched exponential forms

$$p[\mathrm{}(u_n)]\mathrm{exp}(\beta [\mathrm{}(u_n)]^{\alpha _n})\text{for}\mathrm{}(u_n)>\sigma _n,$$
(11)
where $`\sigma _n\sqrt{|u_n|^2}`$. We find $`\alpha _{15}3.45\pm 0.15`$ and $`\alpha _{15}1.85\pm 0.2`$ for the Sabra and SabraP models respectively. We note that in 03BCM  the forcing was Gaussian. We repeated our simulations for a Gaussian force and found results in correspondence with 03BCM , see Fig. 3.

The change in pdf’s indicates that the coupling to the $`B`$-field results in increased intermittency. This can be quantified by measuring the scaling exponents of the standard structure functions $`S_q(k_n)`$,

$$S_q(k_n)|u_n|^qk_n^{\zeta _q}.$$
(12)
Indeed, the scaling exponents of the SabraP models are significantly more nonlinear (as a function of $`q`$) then the corresponding exponents of the Sabra model, cf, Fig. 4.

The simplicity of our model again allows us to offer a qualitative understanding of the increased intermittency. Focus on the Lumley scale $`k_c`$, Eq. 7. We denote $`u_c=u(k_c)`$ and $`B_c=B(k_c)`$. The total “energy” of the B field will be denoted by $`Q`$ where $`Q=_n|B_n|^2`$. In light of the spectrum shown in Fig. 1 we can estimate $`Q|B_c|^2`$ up to coefficients of the order of unity. The flux of energy from the flow to the polymer is proportional to $`k_cQu_c`$. cf. Eq. 2. Therefore we can write an approximate equation, in which we disregard the difference between complex numbers and their moduli:

$$\frac{dQ}{dt}\frac{Q}{\tau }+Qk_cu_c.$$
(13)
From this equation we see that, because $`k_cu_c1/\tau `$, that the fluctuations of Q are expected to be very strong. Moreover, even if $`u_c`$ in Eq. (13) were a Gaussian process, $`Q`$ would turn out to be log-normally distributed, due to the multiplicative role of $`u_c`$ in Eq. (13). Indeed, a direct measurement of the distribution of $`Q`$, see Fig. 5, confirms the existence of the anomalously long tail in this pdf. Finally, due to Eq. (10), we see that an anomalous tail in $`Q`$ would directly induce an anomalous tail in the pdf of $`u_f`$, as seen in the simulations.

Finally, it is interesting to ask whether in 2-dimensional flows we can have enstrophy drag reduction, since this quadratic invariant replaces the energy as the quantity cascaded to smaller scales. In our model this is impossible; our forcing with constant energy flux also produces a constant enstrophy flux, and as can be seen from Eq. (5), the enstrophy is dominated by the larger scales (small $`k_n`$). Thus the reduction in energy is accompanied with the reduction of enstrophy. To switch the role of the scales to get the enstrophy to be dominated by the small scales (large $`k_n`$) requires changing $`b`$ to compensate for the decrease of $`|u_n|^2`$ with increasing $`k_n`$. It is easy to check however that the change of roles occurs precisely at the value of $`b`$ where the inverse cascade of energy disappears, i.e. $`b=12^{2/3}`$. We thus cannot have enstrophy drag reduction in this model, and in our opinion neither in the 2-dimensional Navier-Stokes case.

In summary, we conclude that the results observed for the effects of dilute polymers on 2-dimensional turbulence, as observed in experiments and direct numerical simulations, can be readily understood with a simple model of viscoelastic flows once the cascade of energy is inverted. The same model exhibits drag reduction when the parameters allow a direct energy cascade. These observations lend further support for the use of shell models as simple and transparent tools for understanding a variety of turbulent phenomena that remain more obscure when the full hydrodynamics equations are employed.

###### Acknowledgements.

This work was supported in part by the European Commission under a TMR grant, the Minerva Foundation, Munich, Germany, and the US-Israel Binational Science Foundation.


format: text

text:
# Hybrid atomistic-continuum methods for multiscale hydrodynamics

## 1 Introduction

By limiting the molecular treatment to regions where it is needed, a hybrid method allows the simulation of complex thermo-fluid phenomena which require modeling at the microscale without the prohibitive cost of a fully molecular calculation. In what follows we provide an overview of this rapidly expanding field and discuss recent developments. We also present archetypal hybrid methods for incompressible and compressible flows; the hybrid method for incompressible gas flow is based on the Schwarz alternating coupling method which uses Chapman-Enskog boundary condition imposition, the hybrid method for compressible flow is the recently developed  flux-coupling based, multispecies adaptive mesh and algorithm refinement scheme that extends adaptive mesh refinement by introducing the molecular description at the finest level of refinement.

Over the years a fair number of hybrid simulation frameworks has been proposed leading to some confusion over the relative merits and applicability of each approach. Original hybrid methods focused on dilute gases , which are arguably easier to deal with within a hybrid framework than dense fluids, mainly because boundary condition imposition is significantly easier in gases. The first hybrid methods for dense fluids appeared a few years later . These initial attempts have led to a better understanding of the challenges associated with hybrid methods.

To a large extent, the two major issues in developing a hybrid method is the choice of a coupling method and the imposition of boundary conditions on the molecular simulation. Generally speaking, these two can be viewed as decoupled, in the sense that the coupling technique can be developed on the basis of matching two compatible and equivalent over some region of space hydrodynamic descriptions and can thus be borrowed from the already existing and extensive continuum-based numerical methods literature. The choice of coupling technique is further discussed in section 2.1. Boundary condition imposition can again be considered in a decoupled sense and can be posed as a general problem of imposing “macroscopic” boundary conditions on a molecular simulation. In our opinion, this is a very challenging problem that has not been, in general, resolved to date completely satisfactorily. Boundary condition imposition on the molecular sub-domain is discussed in section 2.2. Boundary condition imposition on the continuum sub-domain is generally well understood, as is the process of extracting macroscopic fields from molecular simulations (typically achieved through averaging).

In section 3 we give a brief description of the direct simulation Monte Carlo (DSMC), the dilute gas simulation method used in this work. In section 4 we demonstrate a hybrid scheme suitable for low speed, incompressible gaseous flows based on the Schwarz alternating method. This scheme introduces Chapman-Enskog boundary condition imposition in incompressible hybrid formulations. Subsequently, in section 5 we discuss a recently developed  multi species compressible formulation for gases which introduces the molecular simulation as the finest level of refinement within a fully adaptive mesh refinement scheme. We finish with some concluding remarks.

## 2 Developing a hybrid method

### 2.1 The choice of coupling method

Coupling a continuum to a molecular description is meaningful in a region where both can be presumed valid. In choosing a coupling method it is therefore convenient to draw upon the wealth of experience and large cadre of coupling methods nearly 50 years of continuum computational fluid dynamics have brought us. Coupling methods for the compressible and incompressible formulations generally differ, since the two correspond to two different physical and mathematical hydrodynamic limits. Faithful to their mathematical formulations, the compressible formulation lends itself naturally to time explicit flux-based (control-volume-type) coupling while incompressible formulations are typically coupled using either state properties (Dirichlet) or gradient information (Neumann).

Given that the two formulations have different limits of applicability and physical regimes in which each is significantly more efficient than the other, care must be exercised when selecting the ingredients of the hybrid method. In other words, the choice of a coupling method and continuum sub-domain formulation needs to be based on the degree on which compressibility effects are important in the problem of interest and not on a preset notion that a particular coupling method is more appropriate than all others. The latter approach was recently pursued in a variety of studies which enforce the use of a control-volume-type approach to steady and essentially incompressible problems to achieve coupling by time explicit flux matching. This approach is not recommended. On the contrary, for an efficient simulation method, similarly to the case of continuum solution methods, it is important to allow the flow physics to dictate the appropriate formulation, while the numerical implementation is chosen to cater to the particular requirements of the latter. Below, we expand on some of the considerations which influence the choice of coupling method under the assumption that the hybrid method is applied to problems of practical interest and therefore the continuum subdomain is appropriately large. Our discussion focuses on timescale considerations that are more complex but equally important to limitations resulting from lengthscale considerations, such as the size of the molecular region(s).

It is well known  that the timestep for explicit integration of the compressible formulation, $`\tau _c`$, scales with the physical timestep of the problem, $`\tau _{\mathrm{\Delta }x}`$($`=\mathrm{\Delta }x/U`$, where $`\mathrm{\Delta }x`$ is the numerical grid spacing and $`U`$ is the characteristic velocity), according to

$$\tau _c\frac{M}{1+M}\tau _{\mathrm{\Delta }x}$$
(1)
where $`M`$ is the Mach number. As the Mach number becomes small, we are faced with the well-known stiffness problem whereby the numerical efficiency of the solution method suffers  because the compressible formulation resolves the acoustic modes when those are not important. For this reason, when the Mach number is small, the incompressible formulation is used which allows integration at the physical timestep $`\tau _{\mathrm{\Delta }x}`$. In the hybrid case matters are complicated by the introduction of the molecular integration timestep, $`\tau _m`$, which is at most of the order of $`\tau _c`$ (in some cases in gases when $`\mathrm{\Delta }x\lambda `$, where $`\lambda `$ is the molecular mean free path) and in most cases significantly smaller. One consequence of equation (1) is that as the global domain of interest grows, the total integration time grows, and transient calculations in which the molecular subdomain is explicitly integrated in time become more computationally expensive and eventually infeasible. The severity of this problem increases with decreasing Mach number and makes unsteady incompressible problems very computationally expensive. New integrative frameworks which coarse grain the time integration of the molecular subdomain are therefore required.

Fortunately, for low speed steady problems implicit (iterative) methods exist which provide solutions without the need for explicit integration of the molecular domain to the global problem steady state. The particular method used here is known as the Schwarz method and is discussed further in Section 4. This method decouples the global evolution timescale from the molecular evolution timescale (and timestep) by achieving convergence to the global problem steady state through an iteration between steady state solutions of the continuum and molecular subdomains. Because the molecular subdomain is small, explicit integration to its steady state is feasible. Although the steady assumption may appear restrictive, it is interesting to note that the vast majority of both compressible and incompressible test problems solved to date and all incompressible practical problems of interest  solved by hybrid methods have been steady. A variety of other iterative methods may be suitable as they provide for timescale decoupling. The choice of the Schwarz coupling method, which uses state variables instead of fluxes to achieve matching, was motivated by the fact (as explained below) that state variables suffer from smaller statistical noise and are thus easier to prescribe on a continuum formulation.

The above observations do not preclude the use of the compressible formulation in the continuum subdomain for low speed flows. In fact, preconditioning techniques which allow the use of the compressible formulation at very low Mach numbers have been developed . Such a formulation can, in principle, be used to solve the continuum sub-problem while this is being coupled to the molecular sub-problem via an implicit (eg. Schwarz) iteration. What should be avoided is a time-explicit control-volume-type coupling procedure for solving essentially incompressible steady state problems.

The issues discussed above have not been very apparent to date because in typical test problems published so far, the continuum and atomistic subdomains are of the same size (and, of course, small). In this case the large cost of the molecular subdomain masks the cost of the continuum subdomain and also typical evolution timescales (or times to steady state) are small. It should not be forgotten, however, that hybrid methods make sense when the continuum subdomain is significantly larger than the molecular subdomain.

Although the continuum subdomain stiffness in the limit $`M0`$ (see eq (1)) may be remedied by implicit timestepping methods  or preconditioning approaches, control-volume-based hybrid approaches suffer from adverse signal to noise ratios in connection with the averaging required for imposition of boundary conditions from the molecular sub-domain to the continuum sub-domain. In the case of an ideal gas (where compressible formulations are typical) it has been shown in  that, for the same number of samples, flux (shear stress, heat flux) averaging exhibits relative noise, $`E_f`$, which scales as

$$E_f\frac{E_{sv}}{Kn}$$
(2)
where $`E_{sv}`$ is the relative noise in the corresponding state variable (velocity, temperature). Here $`Kn=\lambda /L`$ is the Knudsen number based on the characteristic lengthscale of the transport gradients, $`L`$, and $`\lambda `$ is the mean free path which is expected to be much smaller than $`L`$ since, by assumption, a continuum sub-domain is present. It thus appears that flux coupling will be significantly disadvantaged in this case, since the number of samples required to make $`E_fE_{sv}`$ scales as $`1/Kn^2`$ times the number of samples required by state-variable averaging. On the other hand, Schwarz-type iterative methods based on the incompressible physics of the flow require a fair number of iterations for convergence ($`O(10)`$). These iterations require the re-evaluation of the molecular solution. This is an additional computational cost that is not shared by control-volume-type or explicit incompressible approaches. At this time, the choice between a time-explicit formulation or an iterative (Schwarz-type) iteration for incompressible unsteady problems is not clear and may be problem dependent. Despite the fact that as L grows the advantage seems to shift towards iterative methods, we should recall that from equation (1), unless time coarse-graining techniques are developed, large, low-speed, unsteady problems are currently too expensive to be feasible by either approach.

### 2.2 Boundary condition imposition

Consider the molecular region $`\mathrm{\Omega }`$ on the boundary of which, $`\mathrm{\Omega }`$, we wish to impose a set of hydrodynamic (macroscopic) boundary conditions. Typical implementations require the use of particle reservoirs $``$ (see Fig 1) in which particle dynamics may be altered in such a way that the desired boundary conditions appear on $`\mathrm{\Omega }`$; the hope is that the influence of the perturbed dynamics in the reservoir regions decays sufficiently fast and does not propagate into the region of interest, that is, the relaxation distance both for the velocity distribution function and the fluid structure is small compared to the characteristic size of $`\mathrm{\Omega }`$.

In a dilute gas, the non-equilibrium distribution function in the continuum limit has been characterized  and is known as the Chapman-Enskog distribution. Use of this distribution to impose boundary conditions on molecular simulations of dilute gases results in a robust, accurate and theoretically elegant approach. Typical implementations  require the use of particle generation and initialization within $``$. Particles that move into $`\mathrm{\Omega }`$ within the simulation timestep are added to the simulation whereas particles remaining in $``$ are discarded. See chapter 4 for details.

Unfortunately, for dense fluids where not only the particle velocities but also the fluid structure is important and needs to be imposed, no theoretical results for their distributions exist. A related issue is that of domain termination; due to particle interactions, $`\mathrm{\Omega }`$, or in the presence of a reservoir, $``$, needs to be terminated in a way that does not have a big effect on the fluid state inside of $`\mathrm{\Omega }`$.

As a result, researchers have experimented with possible methods to impose boundary conditions. It is now known that similarly to a dilute gas, use of a Maxwell-Boltzmann distribution for the velocities leads to slip . Li et al.  used a Chapman-Enskog distribution to impose boundary conditions to generate a dense-fluid shear flow. In this approach, particles crossing $`\mathrm{\Omega }`$ acquire velocities that are drawn from a Chapman-Enskog distribution parametrized by the local values of the required velocity and stress boundary condition. Although this approach was only tested for a Couette flow, it appears to give reasonable results (within molecular fluctuations). Because in Couette flow no flow normal to $`\mathrm{\Omega }`$ exists, $`\mathrm{\Omega }`$ can be used as symmetry boundary separating two back-to-back shear flows; this sidesteps the issue of domain termination.

In a different approach, Flekkoy et al.  use external forces to impose boundary conditions. More specifically, in the reservoir region they apply an external field of such magnitude that the total force on the fluid particles in the reservoir region is the one required by momentum conservation. They then terminate their reservoir region by using an ad-hoc weighing factor for the distribution of this force on particles within $``$ that prevents particles from leaving the reservoir region. In particular, they chose a weighing factor that diverges as particles approach the edge of $``$ such that particles do not escape the reservoir region while particles introduced there move towards $`\mathrm{\Omega }`$. Particles introduced into the reservoir are given velocities drawn from a Maxwell-Boltzmann distribution, while a Langevin thermostat keeps the temperature constant. The method appears to be successful although the non-unique (ad-hoc) choice of force fields and Maxwell-Boltzmann distribution makes it not very theoretically pleasing. It is also not clear what the effect of these forces are on the local fluid state (it is well known that even in a dilute gas  gravity driven flow exhibits significant deviations from Navier-Stokes behavior) but this effect is probably negligible since force fields are only acting in the reservoir region. Delgado-Buscalioni and Coveney  refined the above approach by using an Usher algorithm to insert particles in the energy landscape such that they have the desired specific energy, which is beneficial to imposing a desired energy current while eliminating the risk of particle overlap at some computational cost. This approach uses a Maxwell–Boltzmann distribution however for the initial velocities of the inserted particles. Temperature gradients are imposed by a small number of thermostats placed in the direction of the gradient. Although no proof exists that the disturbance to the particle dynamics is small, it appears that this technique is successful at imposing boundary conditions with moderate error. Boundary conditions on MD simulations can also be imposed through the method of constraint dynamics . Although the approach in  did not allow hydrodynamic fluxes across the matching interface, this feature can be integrated into this approach with a suitable domain termination.

A method for terminating molecular dynamics simulations with small effect on particle dynamics has been suggested and used in . This simply involves making the reservoir region fully periodic. In this manner, the boundary conditions on $`\mathrm{\Omega }`$ also impose a boundary value problem on $``$, where the inflow to $`\mathrm{\Omega }`$ is the outflow from $``$. As $``$ becomes bigger, the gradients in $``$ become smaller and thus the flowfield in $``$ will have a small effect on the solution in $`\mathrm{\Omega }`$. The disadvantage of this method is the number of particles that are needed to fill $``$ as this grows, especially in high dimensions.

We believe that significant contributions can still be made by developing methods to impose boundary conditions in hydrodynamically consistent and, most importantly, rigorous approaches.

## 3 The direct simulation Monte Carlo

The DSMC method was proposed by Bird  in the 1960’s and has been used extensively to model rarefied gas flows. A comprehensive discussion of DSMC can be found in the review article by Alexander et. al. . The DSMC algorithm is based on the assumption that a small number of representative “computational particles” can accurately capture the hydrodynamics of a dilute gas as given by the Boltzmann equation. Air under standard conditions narrowly meets the dilute gas criterion. Empirical results  show that a small number ($`20`$) of computational particles per cubic molecular mean free path is sufficient to capture the relevant physics. This is approximately 2 orders of magnitude smaller than the actual number of gas atoms/molecules contained in the same volume. This is one source of DSMC’s significant computational advantage over a fully molecular simulation.

DSMC solves the Boltzmann equation using a splitting approach: the time evolution of the system is approximated by a sequence of discrete timesteps, $`\mathrm{\Delta }t`$, in which particles undergo successively collisionless advection and collisions. Collisions are performed between randomly chosen particle pairs within small cells of linear size $`\mathrm{\Delta }x`$. The flow solution is determined by averaging the individual particle properties over space and time. This approach has been shown to produce correct solutions of the Boltzmann equation in the limit $`\mathrm{\Delta }x`$, $`\mathrm{\Delta }t0`$ . The splitting approach eliminates the computational cost associated with integrating the equations of motion of all particles, but most importantly allows the timestep to be significantly larger (see also below) than a typical timestep in a hard sphere molecular dynamics simulation. This is another reason why DSMC is significantly more computationally efficient than “brute force” molecular dynamics.

Recent studies  have shown that for steady flows, or flows which are evolving at timescales that are long compared to the molecular relaxation times, a finite timestep leads to a truncation error that manifests itself in the form of timestep-dependent transport coefficients; this error has been shown to be of the order of 5% when the timestep is of the order of a mean free time and goes to zero as $`\mathrm{\Delta }t^2`$. Quadratic dependence of transport coefficients on the collision cell size $`\mathrm{\Delta }x`$ was shown in .

## 4 The Schwarz method for incompressible formulations

Although in some cases compressibility may be important, a large number of applications are typically characterized by flows where use of the incompressible formulation results in a significantly more efficient approach . As explained in the introduction section, our definition of incompressible formulation is based on the flow physics and not on the numerical method used. Although we have used here a finite element discretization based on the incompressible formulation, we believe that a preconditioned compressible formulation could also be used to solve the continuum subdomain problem if it could be successfully matched to the molecular solution through a coupling method which takes into account the elliptic nature of the (low speed) problem to provide solution matching consistent with the flow physics.

Here, matching is achieved through an iterative procedure based on the Schwarz alternating method for the treatment of steady-state problems. The Schwarz method was originally proposed for molecular dynamics-continuum methods in  and extended in , but it is equally applicable to DSMC-continuum hybrid methods . This approach was chosen because of its ability to couple different descriptions through Dirichlet boundary conditions (easier to impose on dense-system molecular simulations compared to flux conditions, because fluxes are non-local in dense systems), and its ability to reach the solution steady state in an implicit manner. The importance of the latter characteristic cannot be overemphasized; the implicit convergence in time guarantees timescale decoupling that is necessary for the solution of macroscopic problems; the integration of molecular trajectories at the molecular timestep for total times corresponding to macroscopic evolution times is, and will for a long time be, infeasible.

Within the Schwarz coupling framework, an overlap region facilitates information exchange between the continuum and atomistic subdomains in the form of Dirichlet boundary conditions. A steady state continuum solution is first obtained using boundary conditions taken from the atomistic subdomain solution. For the first iteration this latter solution can be a guess. A steady state atomistic solution is then found using boundary conditions taken from the continuum subdomain. This exchange of boundary conditions corresponds to a single Schwarz iteration. Successive Schwarz iterations are repeated until convergence, i.e. until the solutions in the two subdomains are identical in the overlap region.

The Schwarz method was recently applied  to the simulation of flow through micromachined filters. These filters have passages that are sufficiently small that require a molecular description for the simulation of the flow through them. Depending on the geometry and number of filter stages the authors have reported computational savings ranging from 2 to 100. The approach in  used a Maxwellian velocity distribution and a “control mechanism” to impose the flowfield on the molecular simulation. This approach, although succesful in quasi one-dimensional flows, is not very general; additionally, it is well known that using a Maxwellian distribution to impose hydrodynamic boundary conditions, in general, if uncorrected will lead to slip (discrepancy between the imposed and observed boundary conditions). General boundary condition imposition on dilute-gas molecular simulations can be performed using the Chapman-Enskog velocity distribution . This approach eliminates the need for a feedback correction since supplying the correct local distribution function eliminates slip. A Chapman-Enskog procedure for the Schwarz method is described below. Extensions of the Schwarz method to time-dependent problems is currently under investigation , although, as shown by equation (1), when the Mach number is low, the disparity between the molecular and hydrodynamic timescales makes this a very stiff problem.

### 4.1 Driven cavity test problem

In this section we discuss the Schwarz alternating method in the context of the solution of the driven cavity problem. We pay particular attention to the imposition of boundary conditions on the DSMC domain using a Chapman-Enskog distribution which is arguably the most rigorous and general approach. For illustration and verification purposes we solve the steady driven cavity problem (see Figure 2), in which the continuum subdomain is described by the Navier-Stokes equations solved by finite element discretization. The hybrid solution is expected to recover the fully continuum solution since the atomistic subdomain is far from solid boundaries and from regions of large velocity gradients. This test therefore provides a consistency check for the scheme.

Standard Dirichlet velocity boundary conditions for a driven cavity problem were applied on the continuum subdomain; the $`u`$ velocity component on the left, right and lower walls were held at zero while the upper wall $`u`$ velocity was set to $`50`$ m/s. The $`v`$ velocity component on all boundaries was set to zero. Despite the high velocity, the flow is essentially incompressible and isothermal. The pressure is scaled by setting the middle node on the lower boundary at atmospheric pressure ($`1.013\times 10^5`$ Pa).

The imposition of boundary conditions on the atomistic subdomain is facilitated by a particle reservoir as shown in Figure 3. Particles are created at locations $`x,y`$ within the reservoir with velocities $`C_x,C_y`$ drawn from a Chapman-Enskog velocity distribution. The Chapman Enskog distribution is generated  by using the mean and gradient of velocities from the continuum solution, that is, the number and spatial distribution of particles in the reservoir are chosen according to the overlying continuum cell mean density and density gradients. After particles are created in the reservoir they move for a single DSMC timestep. Particles that enter DSMC cells are incorporated into the standard convection/collision routines of the DSMC algorithm. Particles that remain in the reservoir are discarded. Particles that leave the DSMC domain are also deleted from the computation.

The rapid convergence of the Schwarz approach is demonstrated in Figure 4. The continuum numerical solution is reached to within $`\pm 10\%`$ at the 3rd Schwarz iteration and to within $`\pm 2\%`$ at the 10th Schwarz iteration. Our error estimate which includes the effects of statistical noise  and discretization error due to finite timestep and cell size is approximately 2.5%. Similar convergence of the $`v`$ velocity field is also observed.

The close agreement with the fully continuum results indicates that the Chapman-Enskog procedure is not only theoretically appropriate but also robust. Despite a Reynolds number of $`Re1`$, the Schwarz method (originally only shown to converge for elliptic problems ) converges with negligible error. This is in agreement with the findings of Liu  who has recently shown that the Schwarz method is expected to converge for $`ReO(1)`$.

## 5 Adaptive mesh and algorithm refinement for compressible formulations

As discussed above, consideration of the compressible equations of motion leads to hybrid methods which differ significantly from their incompressible counterparts. The hyperbolic nature of compressible flows means that steady state formulations typically do not offer a significant computational advantage, and as a result, explicit time integration is the preferred solution method and flux matching is the preferred coupling method. Given that the characteristic evolution time, $`\tau _h`$, scales with the system size, the largest problem that can be captured by a hybrid method is limited by the separation of scales between the molecular integration time and $`\tau _h`$. Local mesh refinement techniques  minimize the regions of space that need to be integrated at small CFL timesteps (due to a fine mesh), such as the regions adjoining the molecular subdomain. Implicit timestepping methods  can also be used to speed up the time integration of the continuum subdomain. Unfortunately, although both approaches enhance the computational efficiency of the continuum sub-problem, they do not alleviate the issues arising from the disparity between the molecular timestep and the total integration time.

As explained in the introduction, overwhelming computational costs can be incurred when using a control-volume-type approach to capture steady phenomena where compressibility effects are negligible as is in most cases in dense fluids. In this case the integration timestep of the continuum subdomain also becomes of the order of the molecular timescale, while the continuum subdomain is, presumably, much larger than the molecular subdomain and evolves at a much longer timescale. For an example, see in , where a compressible formulation for liquid flow results in a CFL timestep of 0.17 $`\tau _{LJ}`$ where $`\tau _{LJ}O(10^{12}\mathrm{s})`$ is the Lennard-Jones timescale. This appears to not have been fully appreciated by various groups which have attempted to develop dense-fluid hybrid methods based on the compressible continuum formulation and control-volume-type matching procedures to solve steady and essentially incompressible problems.

In hybrid continuum-DSMC methods, locally refining the continuum solution cells to the size of DSMC cells leads to a particularly seamless hybrid formulation in which DSMC cells differ from the neighboring continuum cells only by the fact that they are inherently fluctuating. (The DSMC timestep required for accurate solutions–see –is very similar to the CFL timestep of a compressible formulation.) Coupled to a mesh refinement method, this approach can tackle true multiscale phenomena as shown below.

Another characteristic inherent to compressible formulations is the possibility of describing parts of the domain by the Euler equations of motion . In that case, consistent coupling to the molecular formulation can be performed using a Maxwell-Boltzmann distribution .

In a recent paper , Alexander et al. have shown that explicit time-dependent flux-based formulations preserve the fluctuating nature of the molecular description within the molecular regions but the fluctuation amplitude decays rapidly within the continuum regions; correct fluctuation spectra can be obtained in the entire domain by solving a fluctuating hydrodynamics formulation  in the continuum sub-domain.

### 5.1 Fully adaptive mesh and algorithm refinement for a dilute gas

The compressible formulation of Garcia et al. , referred to as AMAR (Adaptive Mesh and Algorithm Refinement), pioneered the use of mesh refinement as a natural framework for the introduction of the molecular description in a hybrid formulation. In AMAR the typical continuum mesh refinement capabilities are supplemented by an algorithmic refinement (continuum to atomistic) based on continuum breakdown criteria. This seamless transition is both theoretically and practically very appealing.

In what follows we briefly discuss a recently developed  fully adaptive AMAR method. In this method DSMC provides an atomistic description of the flow while the compressible two–fluid Euler equations serve as the continuum–scale model. The continuum and atomistic representations are coupled by matching fluxes at the continuum–atomistic interfaces and by proper averaging and interpolation of data between scales. This is performed in three steps; a) the continuum solution values are interpolated to create DSMC particles in the reservoir region, here called buffer cells, b) the conserved quantities in each continuum cell overlaying the DSMC region are replaced by averages over particles in the same region and c) fluxes recorded when particles cross the DSMC interface are used to correct the continuum solution in cells adjacent to the DSMC region. This coupling procedure makes the DSMC region appear as any other level in an AMR grid hierarchy. Similarly to the overlap region described for the Schwarz method above, the Euler solution information is passed to the particles via buffer cells surrounding the DSMC region. At the beginning of each DSMC integration step, particles are created in the buffer cells using the continuum hydrodynamic values.

The above algorithm allows grid and algorithm refinement based on any combination of flow variables and their gradients. Density gradient based refinement has has been found to be generally robust and reliable. Concentration gradients or concentration values within some interval are also effective refinement criteria especially for multi-species flows involving concentration interfaces. In this particular implementation, refinement is triggered by spatial gradients exceeding user defined tolerances. This approach follows from the continuum breakdown parameter method proposed by Bird .

Using the AMR capabilities provided by the Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI) developed at the Lawrence Livermore National Laboratory , the above adaptive framework has been implemented in a fully three-dimensional, massively parallel form in which, multiple molecular (DSMC) patches can be introduced or removed as needed.

Figure 5 shows the adaptive tracking of a shockwave of Mach number 10 used as a validation test for this method. Density gradient based mesh refinement ensures the DSMC region tracks the shock front accurately. Furthermore, as shown in Figure 6 the density profile of the shock wave remains smooth and is devoid of oscillations that are known to plague traditional shock capturing schemes .

## 6 Discussion

One of the most important messages of this paper is that boundary condition imposition on molecular domains is quite independent of the choice of solution coupling approach. As a example, consider the Schwarz method which provides a recipe for making solutions in various subdomains globally consistent subject to exchange of Dirichlet conditions. The imposition of these boundary conditions can be achieved through any method and no certain method is favored by the coupling approach. Flexibility in adopting appropriate elements from previous approaches, and the importance of chosing the coupling method according to the flow physics are key steps to the development of more sophisticated, next-generation hybrid methods.

Although hybrid methods provide significant savings by limiting molecular solutions only to the regions where they are needed, solution of time-evolving problems which span a large range of timescales is still not possible if the molecular domain, however small, needs to be integrated for the total time of interest. New frameworks are therefore required which allow timescale decoupling or coarse grained time evolution of molecular simulations.

Significant computational savings can be obtained by using the incompressible formulation when appropriate for steady problems. Neglect of these simplifications can lead to a problem that is simply intractable when the continuum subdomain is appropriately large. It is interesting to note that, when a hybrid method was used to solve a problem of practical interest  while providing computational savings, the Schwarz method was preferred because it provides a steady solution framework with timescale decoupling.

For dilute gases the Chapman-Enskog distribution provides a robust and accurate method for imposing boundary conditions. Further work is required for the development of similar frameworks for dense liquids.

## 7 Acknowledgements

The authors wish to thank R. Hornung and A. L. Garcia for help with the computations and valuable comments and discussions, and A. T. Patera and B. J. Alder for helpful comments and discussions. This work was supported in part by the Center for Computational Engineering, and the Center for Advanced Scientific Computing, Lawrence Livermore National Laboratory, US Department of Energy, W-7405-ENG-48. The authors also acknowledge the financial support from the University of Singapore through the Singapore-MIT alliance.


text:
# Enhancement of the stability of genetic switches by overlapping upstream regulatory domains

(October 31, 2003)

## Abstract

We study genetic switches formed from pairs of mutually repressing operons. The switch stability is characterised by a well defined lifetime which grows sub-exponentially with the number of copies of the most-expressed transcription factor, in the regime accessible by our numerical simulations. The stability can be markedly enhanced by a suitable choice of overlap between the upstream regulatory domains. Our results suggest that robustness against biochemical noise can provide a selection pressure that drives operons, that regulate each other, together in the course of evolution.

Biochemical networks are the analog computers of life. They allow living cells, to detect, transmit, and amplify environmental signals, as well as integrate different signals in order to recognize patterns in, say, the food supply. Indeed, biochemical networks can perform a variety of computational tasks analogous to electronic circuits. However, their design principles are markedly different. In a biochemical network, computations are performed by molecules that chemically and physically interact with each other. These interactions are stochastic in nature. This becomes particularly important when the concentrations are low. In gene regulatory networks, this is generally the case: not only the DNA, but also the proteins that regulate gene expression are often present in very small numbers, which can be as low as ten, or even fewer. Hence, one would expect that gene regulatory networks, in contrast to electronic circuits, are highly stochastic and error prone McAdams and Arkin (1997); Elowitz and Leibler (2000); Barkai and Leibler (2000); Ozbudak et al. (2002); Elowitz et al. (2002). An important question, therefore, is how the ability to resist biochemical noise constrains the design of the network Elowitz and Leibler (2000); Barkai and Leibler (2000).

In prokaryotes, the expression of operons—groups of contiguous genes that are transcribed into single mRNA molecules—is regulated by the binding of transcription factors (TFs) to upstream regulatory domains on the DNA. A spatial arrangement in which two operons are transcribed in diverging directions allows the upstream regulatory domains to interfere with each other. This affords additional regulatory control. In particular, biochemical noise in the expression of operons can become correlated or anticorrelated. Just as the existence of operons provides for correlated *gene* expression, interference between the regulatory domains of two diverging operons allows for correlated or *anticorrelated* expression of *operons*. Here, we show that this can have a dramatic influence on the stability of gene regulatory networks.

Recently, we have performed a statistical analysis of the spatial distribution of operons on the genome of Escherichia coli col . The analysis identified a large number of motifs in which the operator regions of two operons overlap and interfere col . Among them are well-known examples such as the lysA-lysR and the araBAD-araC operon pairs Wagner (2000). But perhaps the best-known and arguably the most studied example of such a motif is provided by the $`\lambda `$-phage switch, which consists of two adjacent operons that mutually repress each other Ptashne (1992). Here, we study a minimalist model of such a switch, as shown in Fig. 1. In particular, we compare the stability of an ‘exclusive’ (XOR) switch, for which the simultaneous binding of the repressive TFs for both operons is inhibited, to that of a general switch. We find that the exclusive switch is much more stable than the general switch. The basis reason is that one cannot simultaneously turn off both operons in the exclusive switch. This demonstrates the potential importance of such motifs in making gene regulatory networks robust against biochemical noise. Although these constructions seem peculiar to prokaryotes, evidence for correlated and anticorrelated gene expression has also been reported for eukaryotes Willy et al. (2000); Szallasi (2001).

The starting point of our analysis is a set of chemical reactions that constitute the switches shown in Fig. 1. As chemical species, we introduce a pair of TFs which can exist as monomers, $`\mathrm{A}`$ and $`\mathrm{B}`$, or multimers, $`\mathrm{A}_n`$ and $`\mathrm{B}_m`$. The state of the genome is represented by $`\mathrm{O}`$, $`\mathrm{OA}_n`$, etc, depending on the binding of the TF multimers. Adopting a condensed notation in which ‘$`|`$’ indicates alternative sets of reactants and ‘$``$’ indicates that the reactants are not destroyed by the reaction, the set of chemical reactions are

$`n\mathrm{A}\mathrm{A}_n,m\mathrm{B}\mathrm{B}_m,`$ $`(k_\mathrm{f},k_\mathrm{b})`$ (1a)
$`\mathrm{O}+\mathrm{A}_n\mathrm{OA}_n,\mathrm{O}+\mathrm{B}_m\mathrm{OB}_m,`$ $`(k_{\mathrm{on}},k_{\mathrm{off}})`$ (1b)
$`\mathrm{OA}_n+\mathrm{B}_m|\mathrm{OB}_n+\mathrm{A}_m\mathrm{OA}_n\mathrm{B}_m,`$ $`(k_{\mathrm{on}},k_{\mathrm{off}})`$ (1c)
$`\mathrm{O}|\mathrm{OA}_n\mathrm{A},\mathrm{O}|\mathrm{OB}_m\mathrm{B},`$ $`(k_\mathrm{A}),(k_\mathrm{B})`$ (1d)
$`\mathrm{A}|\mathrm{B}\mathrm{}.`$ $`(\mu _\mathrm{A}),(\mu _\mathrm{B})`$ (1e)
These reactions account for, respectively, the formation of multimers, the binding of TF multimers to the genome (Eqs. (1b) and (1c)), the expression of TF monomers, and the degradation of TF monomers. Repression of gene expression is implicit in Eqs. (1d), thus $`\mathrm{A}`$ is expressed if and only if $`\mathrm{B}_m`$ is not bound, etc. Reaction rates are as indicated, and we define equilibrium constants for multimerisation, $`K_\mathrm{d}=k_\mathrm{f}/k_\mathrm{b}`$, and binding to the genome, $`K_\mathrm{b}=k_{\mathrm{on}}/k_{\mathrm{off}}`$.

Whilst detailed and biologically faithful models can be constructed as has been done for the $`\lambda `$-phage switch Arkin et al. (1998); Aurell et al. (2002), the above model is intentionally ‘as simple as possible’. We believe such an approach is as important as detailed biological modelling in elucidating the basic physics behind switches. Thus, for example, we have condensed the details of transcription and translation into a single reaction step in Eqs. (1d), governed by rate coefficients $`k_\mathrm{A}`$ and $`k_\mathrm{B}`$. On the other hand, as Cherry and Adler have shown Cherry and Adler (2000), the TF binding isotherms must satisfy certain criteria in order to make a working switch. In the present model this is effected by introducing co-operativity through the binding of TF multimers rather than monomers.

In our model the genome is in one of four states $`\{\mathrm{O},\mathrm{OA}_n,\mathrm{OB}_m,\mathrm{OA}_n\mathrm{B}_m\}`$. We now include the effect of interference between the upstream regulatory domains by disallowing some of these states (this is in the spirit of simplicity, strictly speaking the effect is to modify the probabilities of the states). Since the empty genome is always a possibility and both $`\mathrm{A}_n`$ and $`\mathrm{B}_m`$ should be allowed to bind otherwise they would not be TFs, it turns out that there are only five possibilities, two of which are related by symmetry. The four distinct cases are shown in Table 1, and are implemented by excluding some of the reactions in Eqs. (1b) and (1c). For example, the exclusive switch is obtained by discarding the reactions in Eqs. (1c) thereby removing the state $`\mathrm{OA}_n\mathrm{B}_m`$.

We first use mean-field theory to analyse the behaviour of Eqs. (1). Switching behaviour corresponds to the appearance of two distinct stable states in the space of TF molecule numbers. Previously, general switches were studied in detail by Cherry and Adler Cherry and Adler (2000), and a specific example of the exclusive switch was studied by Kepler and Elston Kepler and Elston (2001). We extend the analysis of Cherry and Adler to determine where switching behaviour can occur, for all the cases in Table 1. Firstly, for $`n=m=1`$, no switching behaviour can be found for *any* case. This confirms that some form of co-operative binding is required. For the totally co-operative switch though, switching behaviour cannot be found for *any* values of $`n`$ and $`m`$. For the remaining cases, we have analysed in detail the situation for $`n=m=2`$ where both TFs bind as dimers. Fig. 2 shows the regions in the $`(\mu _\mathrm{A}/k_\mathrm{A},\mu _\mathrm{B}/k_\mathrm{B})`$ plane where switching behaviour is found. Clearly, switching behaviour is more extensive for the exclusive switch compared to the general switch, and is strongly suppressed for the partially co-operative switch. Thus we conclude that, at least in mean field theory, the structure of the switch has a strong influence on the extent of switching behaviour.

To go beyond mean field theory, we have simulated the reactions in Eqs. (1) using Gillespie’s kinetic Monte-Carlo scheme which generates sample trajectories appropriate to the chemical master equation Gillespie (1977). We focus on dimerising ($`n=m=2`$) general and exclusive switches, and on the symmetry line $`k_\mathrm{A}=k_\mathrm{B}=k`$ and $`\mu _\mathrm{A}=\mu _\mathrm{B}=\mu `$. We will use the expression rate $`k0.1`$$`1\mathrm{s}^1`$ Aurell et al. (2002) as a unit of (inverse) time and the degradation rate $`\mu `$ as the main control parameter. The choice of the rate constants is biologically motivated, in particular we expect expression to be a slow step and the binding equilibrium to be biased in favour of bound states Aurell et al. (2002). For a baseline set we use $`k_\mathrm{f}/V=k_\mathrm{b}=k_{\mathrm{on}}/V=5k_{\mathrm{off}}=5k`$ ($`K_\mathrm{d}/V=1`$ and $`K_\mathrm{b}/V=5`$), where $`V2\mu \mathrm{m}^3`$ is the cell volume Aurell et al. (2002); we assume one copy of the genome is present.

We monitor the total numbers of the TFs, $`N_\mathrm{A}`$ and $`N_\mathrm{B}`$, including those in dimers and those bound to the genome. If the system is behaving as a switch then we typically see that one of the TFs is strongly repressed compared to the other one. A switching event occurs when the roles of the two TFs flip spontaneously, as shown in Fig. 3.

We can obtain more insight into the switching behaviour by sampling the probability distribution $`P(N_\mathrm{A},N_\mathrm{B})`$ for states in the $`(N_\mathrm{A},N_\mathrm{B})`$ plane, as shown in Fig. 4. Switching behaviour appears as a double maximum in probability in this representation, and the transition state is seen to lie at low numbers of both TFs. Three points are worthy of note. First of all, it is seen that the positions of the two stable steady states do not depend much on the architecture of the switch. This is not surprising, because if one species dominates, both switches will behave similarly. What is perhaps more surprising, is that the pathways for switching are different. The transition paths of the exclusive switch cross the transition state surface at higher values of $`N_\mathrm{A}=N_\mathrm{B}`$, as compared to the general switch. The reason is that in the general switch both genes can be repressed simultaneously, while in the exclusive switch only one gene can be turned off at the time. More importantly, however, the barrier for flipping the switch is higher for the exclusive switch than for the general switch, as can be seen in the insets in Fig. 4. This is because for a switch to flip, two events have to happen. First of all, the system has to wait for a rare fluctuation by which the concentration of the dominant species decreases; this allows for the synthesis of the other component. Subsequently, the latter component has to bind to its operator site in order to toggle the switch. In the general switch, the latter event is more probable, because the minor component can bind to its site as soon as it is synthesized, while in the exclusive switch the dominant species first has to dissociate from the DNA. This is the main reason why the exclusive switch is more stable than the general switch.

We have also characterised the switching dynamics, by constructing the cumulative distribution function $`F(\mathrm{\Delta }t)`$ for the time intervals $`\mathrm{\Delta }t`$ between zero crossings of the order parameter $`N_\mathrm{A}N_\mathrm{B}`$. About 50% of $`F`$ arises from noise on a time scale $`\mathrm{\Delta }tk^1`$ as the system jitters around the transition state, but for $`\mathrm{\Delta }tk^1`$, and provided we are well into the switching regime, we invariably see Poisson statistics with $`F1\mathrm{exp}[\mathrm{\Delta }t/\tau ]`$ (see Fig. 3(b)). This firstly confirms that the switch states *have* a well defined lifetime $`\tau `$, and secondly allows us to extract an accurate estimate of the value of $`\tau `$.

Bialek has suggested that the switch lifetime may grow exponentially with the number of molecules involved in switching between states Bia . Motivated by this, we monitor the mean number $`\overline{N}`$ of the most-expressed TF, defined to be the time average of $`\mathrm{max}(N_\mathrm{A},N_\mathrm{B})`$. We can also calculate $`\overline{N}`$ from mean field theory, and we find good agreement between this and the value measured in the simulations, as $`\mu `$ varies.

Qualitative support for Bialek’s conjecture comes from Fig. 5. It shows that $`\tau `$ grows very rapidly with $`\overline{N}`$, which is the basic reason why extremely stable switches can be built with at most a few hundred expressed proteins. In contrast to Bialek’s conjecture, however, $`\tau (\overline{N})`$ appears to be sub-exponential in $`\overline{N}`$. Interestingly, we can fit $`\tau (\overline{N})`$ to a form that is suggested by an analysis of a related problem, namely that of switching between broken-symmetry phases in a driven diffusive model Evans et al. (1995); Godrèche et al. (1995). This suggests that the ultimate scaling is $`\tau \overline{N}{}_{}{}^{\alpha }\mathrm{exp}[b\overline{N}]`$, where $`\alpha `$ and $`b`$ are constants. Note that this corresponds to Bialek’s conjecture, but with a logarithmic correction in $`\overline{N}`$.

More importantly, however, Fig. 5 clearly demonstrates that the switch construction does indeed have a marked influence on the stability of the switch. It shows that the lifetime of the exclusive switch grows much more rapidly with mean copy number than that of the general switch not . Our simulations cover $`10\overline{N}30`$, but if we extrapolate our results to $`\overline{N}100`$, then $`k\tau 10^4`$$`10^6`$ for the general switch but $`k\tau 10^8`$$`10^{10}`$ for the exclusive switch. In the latter case, this corresponds to lifetimes measured in tens of years. Such extremely long lifetimes have been reported for $`\lambda `$-phage Aurell et al. (2002).

In summary, a genetic switch is intrinsically stochastic, because of the molecular character of its components. However, our simulations demonstrate that the stability of a genetic switch can be strongly enhanced by spatially arranging the operons such that *competing* regulatory molecules *mutually exclude* each other at the operator regions. Such a spatial arrangement can be achieved if the two operons lie next to each other on the DNA and are transcribed in diverging directions – a network motif that has been identified by our statistical analysis of the gene regulatory network of E. coli col . Hence, our simulations suggest that robustness against biochemical noise can provide a selection pressure that drives pairs of operons, that either regulate each other or are controlled by a common transcription factor, towards each other in the course of evolution.

We thank Daan Frenkel, Fred MacKintosh and Sander Tans for useful discussions and for a critical reading of the manuscript, and Martin Evans for drawing our attention to the work on driven diffusive systems. This work is supported by the Amsterdam Centre for Computational Science (ACCS). The work is part of the research program of the “Stichting voor Fundamenteel Onderzoek der Materie (FOM)”, which is financially supported by the “Nederlandse organisatie voor Wetenschappelijk Onderzoek (NWO)”.


text:
# X-rays from the Dawn of the Modern Universe. Chandra and XMM-Newton Observations of 𝒛>𝟒 Quasars

## 1. Introduction

Our knowledge of the X-ray properties of quasars at $`z>4`$ has advanced rapidly over the past few years. In particular, the Sloan Digital Sky Survey (SDSS; e.g., York et al. 2000) has generated large and well-defined samples of $`z>4`$ quasars (e.g., Anderson et al. 2001); most of these quasars are suitable for X-ray studies. The X-ray observational strategy has comprised archival studies of high-redshift quasars with ROSAT (Kaspi, Brandt, & Schneider 2000; Vignali et al. 2001, hereafter V01), snapshot ($``$ 4–10 ks) observations with Chandra to define basic quasar X-ray properties such as fluxes and luminosities (e.g., V01; Brandt et al. 2002, 2003; Bechtold et al. 2003; Vignali et al. 2003a,b, hereafter V03a, V03b) and longer observations with XMM-Newton to derive either tight constraints on the X-ray emission (e.g., Brandt et al. 2001) or spectral parameters by direct X-ray fitting (Ferrero & Brinkmann 2003; Grupe et al. 2004). Chandra snapshot observations have also allowed joint spectral fitting of subsamples of quasars drawn from two main samples at $`z>4`$: the optically luminous Palomar Digital Sky Survey (e.g., Djorgovski et al. 1998) and the SDSS. The X-ray spectral results provide no evidence of strong spectral evolution in radio-quiet quasar (RQQ) X-ray emission from local samples up to $`z`$ 5; the spectrum at high redshift is well parameterized by a power law in the $``$ 2–40 keV rest-frame band with $`\mathrm{\Gamma }=1.8`$–2 (V03a; V03b). Furthermore, no evidence for widespread intrinsic X-ray absorption has been found, although it seems likely that a few individual objects may be X-ray absorbed (e.g., V01; V03b). These overall results have been supported recently by direct X-ray spectroscopy of QSO 0000$``$263 at $`z=4.10`$ with XMM-Newton (Ferrero & Brinkmann 2003).

The color selection of the SDSS has been proven to be effective in finding high-redshift optically luminous quasars up to $`z`$ 5.7 (see Fan et al. 2003 for SDSS quasars at higher redshifts). On the other hand, moderately deep Chandra observations and the ultra-deep (2 Ms) survey of the Chandra Deep Field-North (CDF-N; Alexander et al. 2003) can detect Active Galactic Nuclei (AGN) at $`z>4`$ that are typically $`\stackrel{>}{}`$ 10–30 times less luminous than the SDSS quasars (e.g., Barger et al. 2002; Silverman et al. 2002; Vignali et al. 2002, hereafter V02; Castander et al. 2003). These AGN are much more numerous and therefore more representative of the AGN population at high redshift than the rare SDSS quasars; however, their X-ray emission does not appear to contribute significantly to reionization at $`z`$ 6 (Barger et al. 2003). A detailed X-ray spectral analysis of the $`z>4`$ AGN in the CDF-N is presented in V02.

Below we present some new X-ray spectral results obtained by joint spectral fitting of all the RQQs at $`z>4`$ thus far detected by Chandra. A spectral analysis performed on a smaller but more X-ray luminous sample of $`z>4`$ radio-loud quasars is presented by Bassett et al. (in preparation).

## 2. Joint X-ray Spectral Results

To define the overall X-ray properties of $`z>4`$ RQQs, we selected all of the RQQs detected by Chandra with $`>2`$ counts in the observed 0.5–8 keV band. The sample comprises 46 quasars with a median redshift of 4.43; the number of source counts is $``$ 750. Note that these quasars represent a large fraction ($``$ 70%) of the optically selected RQQs at $`z>4`$ with X-ray detections at present.111See http://www.astro.psu.edu/users/niel/papers/highz-xray-detected.dat for a regularly updated listing of X-ray detections and sensitive upper limits at $`z>4`$. Although it is possible that individual objects are characterized by “peculiar” X-ray properties, our approach obtains average spectral parameters for the quasar population at $`z>4`$ using a much larger sample than those presented in V03a and V03b. In the X-ray spectral analysis, the Cash statistic (Cash 1979) has been adopted. Our preliminary analysis shows that a power law fits the X-ray data reasonably well; the photon index in the rest-frame $``$ 2–40 keV band is $`\mathrm{\Gamma }=1.9\pm 0.1`$ (see Fig. 1).

This is consistent with previous results obtained for RQQ samples at high redshift observed with Chandra (V03a; V03b) and XMM-Newton (Ferrero & Brinkmann 2003; Grupe et al. 2004), as well as with quasar X-ray spectral results at low and intermediate redshift (e.g., George et al. 2000; Page et al 2003). Our analyses indicate that the X-ray spectral properties of $`z>4`$ RQQs and local RQQs are similar; the only significant differences have been found in their broad-band properties using the SDSS Early Data Release quasar catalog (Schneider et al. 2002) and are likely due to luminosity effects (Vignali, Brandt, & Schneider 2003; see also Brandt, Schneider, & Vignali, these proceedings). Thus, despite the strong changes in large-scale environment and quasar number density that have occurred from $`z`$ 0–6, individual quasar X-ray emission regions appear to evolve relatively little. From the joint X-ray spectral fitting we also find no significant evidence for absorption above the Galactic value; the upper limit in the source rest frame is $`N_\mathrm{H}\stackrel{<}{}9\times 10^{20}`$ cm-2 (see Vignali et al., in preparation, for detailed discussion).

## 3. The Future

The correlation found between quasar AB magnitude at a rest-frame wavelength of 1450 Å and the observed 0.5–2 keV flux (e.g., V03b) is a powerful tool to select samples of $`z>4`$ quasars suitable for follow-up X-ray observations. The combination of snapshot observations with Chandra and longer exposures with XMM-Newton should continue to be highly effective in allowing the study of the overall X-ray properties of quasars at high redshift. In the coming years, as the SDSS is completed and several thousand Chandra and XMM-Newton archival observations become available to the scientific community, our knowledge of the broad-band properties of quasars at the highest redshifts will significantly increase. However, detailed X-ray spectroscopic analyses of large samples of $`z>4`$ quasars must await the more distant future and X-ray missions such as Constellation-X, XEUS, and Generation-X.

### Acknowledgments.

Support from the Italian Space Agency under contract ASI I/R/073/01 (CV), NASA LTSA grant NAG5-13035 (WNB, DPS), NSF CAREER award AST-9983783 (WNB), and NSF grant AST-9900703 (DPS) is gratefully acknowledged.

## References

Alexander, D.M., et al. 2003, AJ, 126, 539

Anderson, S.F., et al. 2001, AJ, 122, 503

Barger, A.J., et al. 2002, AJ, 124, 1839

Barger, A.J., et al. 2003, AJ, 126, 632

Bechtold, J., et al. 2003, ApJ, 588, 119

Brandt, W.N., Guainazzi, M., Kaspi, S., Fan, X., Schneider, D.P., Strauss, M.A., Clavel, J., & Gunn, J.E. 2001, AJ, 121, 591

Brandt, W.N., et al. 2002, ApJ, 569, L5

Brandt, W.N., et al. 2003, in New X-ray Results from Clusters of Galaxies and Black Holes, ed. C. Done, E.M. Puchnarewicz, & M.J. Ward, in press (astro-ph/0212082)

Cash, W. 1979, ApJ, 228, 939

Castander, F.J., Treister, E., Maccarone, T.J., Coppi, P.S., Maza, J., Zepf, S.E., & Guzman, R. 2003, AJ, 125, 1689

Djorgovski, S.G., et al. 1998, in Wide Field Surveys in Cosmology, ed. S. Colombi, & Y. Mellier (Paris: Editions Frontieres), 89

Fan, X., et al. 2003, AJ, 125, 1649

Ferrero, E., & Brinkmann, W. 2003, A&A, 402, 465

George, I.M., Turner, T.J., Yaqoob, T., Netzer, H., Laor, A., Mushotzky, R.F., Nandra, K., & Takahashi, T. 2000, ApJ, 531, 52

Grupe, D., Mathur, S., Wilkes, B., & Elvis, M. 2004, AJ, in press (astro-ph/0310188)

Kaspi, S., Brandt, W.N., & Schneider, D.P. 2000, AJ, 119, 2031

Page, K.L., Turner, M.J.L., Reeves, J.N., O’Brien, P.T., & Sembay, S. 2003, MNRAS, 338, 1004

Schneider, D.P., et al. 2002, AJ, 123, 567

Silverman, J.D., et al. 2002, ApJ, 569, L1

Vignali, C., Brandt, W.N., Fan, X., Gunn, J.E., Kaspi, S., Schneider, D.P., & Strauss, M.A. 2001, AJ, 122, 2143 (V01)

Vignali, C., Bauer, F.E., Alexander, D.M., Brandt, W.N., Hornschemeier, A.E., Schneider, D.P., & Garmire, G.P. 2002, ApJ, 580, L105 (V02)

Vignali, C., Brandt, W.N., & Schneider, D.P. 2003, AJ, 125, 433

Vignali, C., Brandt, W.N., Schneider, D.P., Garmire, G.P., & Kaspi, S. 2003a, AJ, 125, 418 (V03a)

Vignali, C., et al. 2003b, AJ, 125, 2876 (V03b)

York, D.G., et al. 2000, AJ, 120, 1579


text:
# Spatio-temporal oscillations and rheochaos in a simple model of shear banding

## Abstract

We study a simple model of shear banding in which the flow-induced phase is destabilised by coupling between flow and microstructure (wormlike micellar length). By varying the strength of the instability and the applied shear rate, we find a rich variety of oscillatory and rheochaotic shear banded flows. At low shear and weak instability, the induced phase pulsates in width next to one wall of the flow cell. For stronger instability, single or multiple high shear pulses ricochet across the cell. At high shear rates we observe oscillating bands on either side of a defect. In some cases, multiple such defects exist and propagate across the cell to interact with each other. We discuss our results in the context of recent observations of oscillating and fluctuating shear banded flows.

Complex fluids commonly undergo flow instabilities and flow-induced phase transitions that result in spatially heterogeneous “shear banded” states. Classically studied systems include wormlike micellar surfactants Britton and Callaghan (1997); dense lamellar onion Diat et al. (1993) or micellar cubic Eiser et al. (2000) phases; and polymer solutions Hilliou and Vlassopoulos (2002). Fluidity banding has also been reported in soft glassy materials Sollich et al. (1997) such as colloidal suspensions Raynaud et al. (2002) and simulated Lennard Jones particles Varnik et al. (2003). Experimentally, the basic observation is of two coexisting shear bands with differing viscosity and microstructure (or fluidity). Theoretically, this is captured by invoking multiple flow branches in the constitutive relation of shear stress vs. shear rate, $`\mathrm{\Sigma }(\dot{\gamma })`$ Olmsted and Goldbart (1990); Spenley et al. (1993); Lu et al. (2000a). The system then separates into a steady state comprising two shear bands, each on its own flow branch (Fig. 1a).

However an accumulating body of data shows this basic picture to be oversimplified: Many shear banding systems display oscillations or irregular fluctuations suggesting chaos in their bulk rheology, rheo-optics or velocimetry. Example systems include onion phases (at the layering transition sal  or showing size oscillations sal (b)); wormlike micelles (WMs) with a stress plateau in the flow curve Bandyopadhyay et al. (2000); shear thinning WMs showing common stress banding (band normals in the flow-gradient direction) cal (1, 2); shear thickening WMs showing common stress Hu et al. (1998) or common strain rate (vorticity) Fischer et al. (2002) banding; and polymer solutions Hilliou and Vlassopoulos (2002). In contrast, present models predict steady banded states Spenley et al. (1993); Olmsted and Goldbart (1990); Lu et al. (2000a).

In the negligible Reynolds limit relevant to these viscoelastic materials, the non linearity underlying this erratic response must arise in the rheological constitutive behaviour of the system Cates et al. (2002). Temporal “rheochaos” has recently been studied in homogeneous models of director dynamics in sheared nematics Grosso et al. (2001); and of shear thickening systems with a single-branched flow curve Cates et al. (2002). In many systems, however, spatial heterogeneity is likely to be a crucial ingredient of rheochaos. In this work, motivated by the above experiments, we introduce the first model of spatio-temporal rheochaos in shear banding systems with multi-branched constitutive curves.

Homogeneous flow is unstable in any region of negative constitutive slope, $`d\mathrm{\Sigma }/d\dot{\gamma }<0`$. This is easily seen in models that take the “mechanical variables” ($`\mathrm{\Sigma }`$, $`\dot{\gamma }`$) as the relevant dynamic variables Yerushalmi et al. (1970). The system can then separate into two shear bands, each on its own stable flow branch. (See Fig. 1a, for shear thinning gradient-banding systems.) In more realistic models, the mechanical variables are coupled to microstructural quantities such as director orientation Olmsted and Goldbart (1990), polymeric concentration Fielding and Olmsted (2003a) or micellar length Cates and Turner (1992). This coupling can destabilize the rising high shear branch (causing, e.g., tumbling and wagging in nematics Doi (1981)). Here we construct a simple model with an unstable high shear branch (Fig. 1b,c) and show that it has oscillatory and erratic (chaotic) shear banded states at imposed global shear rate. Although we have not yet proved the erratic states to be strictly chaotic, we refer to them as chaotic hereafter. Our model resembles other globally coupled reaction diffusion systems, with similar spatio-temporal behaviour nonlin . Aradian and Cates are currently studying spatio-temporal flows of similar models with a single unstable constitutive branch cat .

We use only the minimal ingredients needed to capture the observed phenomena, so do not address the microscopics of any given system. However for concreteness we use the language of shear thinning WMs. There exist numerous reports of apparently unattainable homogeneous high shear rate branches in such systems, in which the sample flows erratically or is ejected from the cell (e.g. Britton and Callaghan (1997)). This is seldom discussed in detail, but sometimes attributed to surface instability. Nonetheless, the possibility of bulk instability remains. Indeed, in some WMs the high shear band breaks into sub-bands Lerouge et al. (2000).

For simplicity we consider just one microstructural variable, the mean micellar length $`n`$, and define our model by

$`\mathrm{\Sigma }`$ $`=`$ $`\sigma +\eta \dot{\gamma },`$ (1)
$`_t\sigma `$ $`=`$ $`{\displaystyle \frac{\sigma }{\tau (n)}}+{\displaystyle \frac{g[\dot{\gamma }\tau (n)]}{\tau (n)}}+D_y^2\sigma ,`$ (2)
$`_tn`$ $`=`$ $`{\displaystyle \frac{n}{\tau _n}}+{\displaystyle \frac{N(\dot{\gamma }\tau _n)}{\tau _n}}.`$ (3)
Eqn. 1 gives the uniform shear stress $`\mathrm{\Sigma }(t)`$ in the non-inertial limit. It comprises a (generally non-uniform) viscoelastic micellar contribution $`\sigma (y,t)`$ and a solvent contribution with viscosity $`\eta `$. The dynamics of $`\sigma `$ (Eqn. 2) has a length-dependent relaxation time Soltero

$$\tau (n)=\tau _0\left(\frac{n}{n_0}\right)^\alpha ,$$
(4)
and a steady homogeneous state $`\sigma =g[\dot{\gamma }\tau (n)]`$ set by

$$g(x)=\frac{x}{1+x^2},$$
(5)
admitting a constitutive curve of negative slope. The spatial gradients $`D_y^2\sigma `$ allow for a smooth interface of width $`l\sqrt{D}`$ between the shear bands.

The micellar length $`n`$ has its own relaxation time $`\tau _n`$ (Eqn. 3), presumably related to the (unknown) underlying rates of micellar scission and recombination. For example, scission could be enhanced by the agitation of shearing; or recombination aided by shear-induced alignment of micellar ends Cates and Turner (1992). We assume the former, taking

$$N(x)=\frac{n_0}{1+x^\beta }.$$
(6)
Because this decrease in length only becomes important for typical shear rates $`\dot{\gamma }>1/\tau _n`$, the strength of coupling between the mechanical quantities $`(\mathrm{\Sigma },\sigma ,\dot{\gamma })`$ and microstructure $`n`$ can be tuned by varying $`\tau _n`$.

Using this model we study flow between two parallel plates at $`y=0,L`$ with boundary condition $`_y\mathrm{\Sigma }=0`$, using units in which $`n_0=1,\tau _0=1`$ and $`L=1`$. We set $`\alpha =1.2`$, $`\beta =1.5`$, though our results are qualitatively robust to reasonable variations in these values.

Homogeneous dynamics— The constitutive curves $`\mathrm{\Sigma }(\dot{\gamma })`$ and domains of instability for homogeneous states are shown in Fig. 1. For small $`\tau _n0.115`$ we find pure mechanical instability (Fig. 1a) that does not involve $`n`$. The constitutive curve is then (at any fixed $`\mathrm{\Sigma }`$) an unstable saddle Strogatz (1994) (one unstable eigenvector). For larger $`\tau _n`$, coupling to micellar length broadens this instability into the rising high shear branch, which is now an unstable focus Strogatz (1994) (two unstable eigenvectors): Fig. 1b,c. This instability terminates in a Hopf bifurcation Strogatz (1994).

We then studied the non-linear dynamics of the local ($`D=0`$) model for fixed $`\mathrm{\Sigma }`$, solving Eqns. (1-3) via a fourth order Runge-Kutte method Press et al. (1992). This confirmed the stability properties of Fig. 1: states near an unstable (stable) segment of the constitutive curve flow away from (towards) that segment (Fig. 2). They also reveal a periodic orbit about the unstable high shear branch for stresses just below the Hopf bifurcation. We also used AUTO Doedel (1997) to trace the amplitude of the periodic orbits, Fig. 1b,c. Periodic orbits are the most complex behaviour possible for the local model since it has only two degrees of freedom, $`d=2`$. Chaos requires $`d3`$ Strogatz (1994).

Spatially heterogeneous dynamics— We now turn to the non-local model, $`D0`$, focusing on the implications of an unstable high shear branch (Fig. 1b,c). The dimensionality $`d`$ is now effectively infinite, since each spatial point has its own value of $`n`$ and $`\dot{\gamma }`$. We solved the non-local equations using a Crank-Nicholson algorithm (checking our results with the Rosenbrock method) Press et al. (1992), with the constraint of fixed global strain rate, $`\overline{\dot{\gamma }}=_0^1𝑑y\dot{\gamma }(y,t)`$. For small $`\tau _n0.115`$ we find stable shear bands (Fig. 1a). In contrast for $`\tau _n=0.145`$ (unstable high shear branch) we find spatio-temporal oscillations and chaotic banding flows (Fig. 3).

Several regimes are evident. At low applied shear, $`\overline{\dot{\gamma }}=1.5`$, a thin pulse of high shear ricochets back and forth across the cell (Fig. 3, top). A thin fluctuating high shear band, away from the rheometer wall, was seen experimentally in Ref. cal (1). At larger shear rates, we find two or more such pulses. For two pulses (not shown), we typically find a periodically repeating movie with the pulses alternately bouncing off each other (mid-cell) and the cell walls. Once three pulses are present, $`\overline{\dot{\gamma }}=7.0`$, periodicity gives way to chaotic behaviour (Fig. 3).

At still higher shear, $`\overline{\dot{\gamma }}=19.2`$, we find regular oscillations of spatially extended bands pinned at a stationary defect. The local shear rates span both the low and high shear constitutive branches. Oscillating (vorticity) bands were seen experimentally in Ref. Fischer et al. (2002); Hilliou and Vlassopoulos (2002). For the intermediate value $`\overline{\dot{\gamma }}=11.35`$ we find intermittency between patterns resembling those for $`\overline{\dot{\gamma }}=7.0`$ and $`\overline{\dot{\gamma }}=19.2`$.

Finally for $`\overline{\dot{\gamma }}=23.0,31.0`$ we find oscillating bands separated by moving defects (Fig. 3, bottom right) with the flow now governed only by the high shear constitutive branch. In each band, the shear rate cycles round the periodic orbit of the local model (Fig. 2 and top data set of Fig. 4).

For different values of $`\tau _n`$ we find a host of other interesting phenomena lon . For example, for weaker instability ($`\tau _n=0.13`$) at low applied shear rates we see a high shear band that can pulsate in width while adhering to the rheometer wall (Fig. 5a), or meander about the cell (Fig. 5b). The former behaviour resembles interfacial motion in seen WMs cal (2); Hu et al. (1998) and onion phases sal .

Finally we discuss in more detail the high shear pulse of Fig. 3 (top right). At times when the pulse is far from the wall, the stress $`\mathrm{\Sigma }`$ is constant (flat regions in Fig. 3, top left). In this regime we transform to the pulse’s comoving frame $`\widehat{y}=yct`$ and eliminate $`\sigma `$ to get

$`c\eta \dot{\gamma }^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Sigma }}{\tau (n)}}+{\displaystyle \frac{g[\dot{\gamma }\tau (n)]+\eta \dot{\gamma }}{\tau (n)}}D\eta \dot{\gamma }^{\prime \prime }`$ (7)
$`cn^{}`$ $`=`$ $`{\displaystyle \frac{n}{\tau _n}}+{\displaystyle \frac{N(\dot{\gamma }\tau _n)}{\tau _n}}`$ (8)
in the $`3d`$ space $`(n,\dot{\gamma },\dot{\gamma }^{})`$, with parameters $`\mathrm{\Sigma }`$ and $`c`$. (Prime denotes $`/\widehat{y}`$.) The low shear regions in Fig. 3 (top right) lie on the low shear branch of the constitutive curve. Although this is a stable stationary flow branch (Fig. 1), in the space of $`(n,\dot{\gamma },\dot{\gamma }^{})`$ it is actually a fixed point ($`P`$) with one unstable eigenvector and two stable ones. As $`\widehat{y}`$ increases through the high shear pulse, the solution’s trajectory moves away from $`P`$ along the unstable direction and then back to $`P`$ in the stable plane, completing a “homoclinic connection” Strogatz (1994). The requirement for this $`1d`$ trajectory to correctly rejoin the $`2d`$ stable manifold selects (for a given $`\mathrm{\Sigma }`$) the observed propagation speed, $`c`$. For other values of $`c`$, the trajectory abruptly crosses the $`2d`$ manifold, missing $`P`$. (This follows by simple dimension counting Lu et al. (2000a).) Once selected, the solution’s profile $`\dot{\gamma }(\widehat{y})`$ can (for a given value of $`D`$) be integrated over $`\widehat{y}`$ to find the globally applied shear rate $`\overline{\dot{\gamma }}`$. By the reverse reasoning, we find a selected $`\mathrm{\Sigma }`$ and $`c`$ for any $`\overline{\dot{\gamma }}`$ (and $`D`$). We used AUTO Doedel (1997) to check that the homoclinic connections of Eqn. 8 indeed coincide with our numerical results. The locus of the maxima $`\dot{\gamma }_\mathrm{h}(\mathrm{\Sigma })`$ of these connections is marked in Fig. 4. The space-time points for $`\dot{\gamma }>\dot{\gamma }_\mathrm{h}`$ arise during wall collisions; correspondingly, their density is smaller than for $`\dot{\gamma }<\dot{\gamma }_\mathrm{h}`$.

To summarise, we have constructed a simple model of shear banding in which the high shear branch of the underlying constitutive curve is rendered unstable by a coupling between flow and microstructure. Within this model, we have found a rich variety of spatio-temporal oscillatory and rheochaotic flows, many resembling experimental observations in shear banding systems. It remains an open challenge to delineate more fully the spectrum of mechanisms governing rheochaotic banded states. Extension to higher dimensions, allowing fluctuations along the plane of the interface would be interesting. Within the present one-dimensional approach, generalisation to shear thickening systems cat , curved Couette geometries and vorticity-banding is underway lon .

We thank A. A. Aradian, P. T. Callaghan, M. E. Cates, E. Knobloch, A. M. Rucklidge and J.-B. Salmon for interesting discussions; and EPSRC GR/N11735 for funding.


text:
# Determining the Metallicity of Cepheid Stars in the SMC, LMC and the Galaxy

## 1. Introduction

The question if (and by how much) metallicity has any effect on the Cepheid Period-Luminosity (PL) relation is far from being settled. Recent theoretical studies give different results: some authors claim that metallicity has negligible effects on the PL relation (Baraffe & Alibert 2001, Alibert et al. 1999; Sandage et al. 1999); others find that there is a significant dependence of the PL relation on metallicity (Bono et al. 1999; Caputo et al. 2000; Fiorentino et al. 2002).

Although most of the observational efforts aiming at deriving the chemical composition of Cepheids and its effects on the PL relation have used indirect means, like the measurement of O II in HII regions (e.g. Beaulieu et al. 1997; Sasselov et al. 1997; Kennicutt et al. 1998; Macri et al. 2001), some of them have focused on a full and direct chemical analysis (Fry & Carney 1997; Andrievsky et al. 2002a,b; Luck et al. 2003). The goal of our project is to follow this second approach, i.e. to make a direct determination of the metallicity in a newly observed sample of Cepheid stars.

## 2. The Data Sample

The data-set has been assembled according to the following selection criteria: a large coverage of metallicities (hence, we selected stars from three different galaxies: the SMC, the LMC, and our own); the availability of well-determined intensity-mean magnitudes in the B, V, J, H, and K bands (cf Berdnikov & Caldwell 2001; Groenewegen 1999; Laney & Stobie 1994, Laney & Stobie 1992); a wide range of periods in order to cover the entire PL relation (hence we chose objects from 2 to 98 days of period); the availability of Hipparcos parallaxes for all Galactic Cepheids.

Our sample includes a total of 76 stars: 40 Galactic, 22 LMC, and 14 SMC Cepheids. The Magellanic Cepheids have been observed with UVES (R=40,000, spectral coverage: 480–680 nm) at the ESO-VLT and have typical S/N of the order of 50 on the EEV CCD (i.e. for 480$`\lambda `$580 nm) and $``$70 on the MIT CCD (i.e. for 580$`\lambda `$ 680 nm). The spectra of the Galactic sample, instead, have been collected with FEROS (R=45,000, spectral coverage: 380–880 nm) at the ESO 1.5m telescope in La Silla. The typical S/N of these spectra is $``$70 up to $`\lambda `$=500 nm, but larger than 150 in the red.

## 3. The Analysis

After the spectra were reduced with the help of the FEROS and UVES Data Reduction Software pipelines, our main goal has been to develop a reliable analytical procedure in order to accurately determine the iron abundance of these stars.

First, we carefully assembled a reliable list of iron lines (both neutral and ionized). We started from a preliminary selection from the VALD database (Kupka et al. 1999; Piskunov et al. 1995) choosing lines between 450 and 650 nm, for effective temperatures typical of Cepheids (4500-6500 K). Every line profile was then visually inspected on the observed spectra in order to detect (and eliminate) those lines affected by other elemental blends.

Secondly we measured the equivalent widths (EWs) of all the lines assembled as described above. To do this, we tested two methods: we used the gaussian fitting option offered by the task splot in IRAF; we then developed an IDL routine (based on the line\_width widget available in the FUSE IDL suite of tools) in order to speed up the process, but keeping the whole procedure interactive.

With a final list of reliable equivalent widths in hand, we then made use of the code WIDTH9 (Kurucz 1993) for the computation of our final iron abundances. For this purpose, we chose the Kurucz grid of LTE model atmospheres (Kurucz 1993). However, before being able to derive \[Fe/H\]111where \[Fe/H\]=log(Fe/H) \- log(Fe/H) values from our measured EWs, one clearly needs to know the stellar parameters of the stars, i.e. their effective temperature, gravity, and microtubulent velocity. We derived a first guess set of stellar parameters from the photometry available in the literature (see references given in Section 2) and determine spectroscopically the final set to be used in the analysis. This is an interactive process that requires running the code WIDTH9, with our final list of EWs and different model atmospheres, in order to find which model is able to minimize simultaneously the slope of log$`ϵ`$(Fe) vs EW and vs the excitation potential (these constrain the microturbulence and the temperature respectively) while keeping the ionization balance (i.e. the abundance derived from Fe I and from Fe II must agree)to constrain the gravity. Once such model was identified for each individual star, we were finally able to compute the iron abundance of all our stars.

## 4. Preliminary results and Future Developments

We have determined the iron abundances for a sub-sample of stars: 14 Galactic, 14 LMC and 22 SMC Cepheids. Figure 1 of Romaniello et al. (this volume) summarises this preliminary set of results.

Despite the confidence we have gained in the analytical steps tested sofar, we still have a lot of work in front of us. In the very near future, we want to optimize our line-list (both for the number of lines to be used and for a careful check of their gf-values), to complete the analysis for the entire sample, and to determine iron abundances also via spectrum synthesis. A little bit more into the future, we plan to perform the full chemical analysis of the entire data-set (with special attention devoted to the $`\alpha `$ and to the heavy elements).

## References

Alibert, Y., Baraffe, I., Hautschildt, P.& Allard, F. 1999, A&A, 344, 551

Andrievsky, S.M., Bersier, D., Kovtyukh, V.V., Luck, R.E., Maciel, W.J., Lépine, J.R.D. & Beletsky, Yu.V. 2002a, A&A, 384, 140

Andrievsky, S.M., Kovtyukh, V.V., Luck, R.E., Lépine, J.R.D., Maciel, W.J.& Beletsky, Yu.V. 2002b, A&A, 392, 491

Baraffe, I. & Alibert, Y. 2001, A&A, 371, 592

Beaulieu, J.P., Sasselov, D.D., Renault, C., Grison, P., Ferlet, R. et al. 1997, A&A, 318, L47

Berdnikov, L.N. & Caldwell, J.A.R. 2001, JAD, 7, 3

Bono, G., Caputo, F., Castellani, V. & Marconi, M. 1999, ApJ, 512, 711

Caputo, F., Marconi, M. & Musella, I. 2000, A&A, 354, 610

Fiorentino, G.,Caputo, F., Marconi, M. & Musella, I. 2002, ApJ, 576, 402

Fry, A.M. & Carney, B.W. 1997, AJ, 113, 1073

Groenewegen, M.A.T. 1999, A&ASS, 139, 245

Groenewegen, M.A.T. & Oudmaijer, R.D. 2000, A&A, 356, 849

Kennicutt, R.C., Stetson P.B., Saha, A., Kelson, D., Rawson, D.M, et al. 1998, ApJ, 498, 181

Kupka, F., Piskunov, N.E., Ryabchikova, T.A., Stempels, H.C. & Weiss W.W 1999, A&AS, 138, 119

Kurucz, R.L. 1993, CD-ROMs #1, 13, 18

Laney, C.D. & Stobie, R.S. 1992, A&ASS, 93, 93

Laney, C.D. & Stobie, R.S. 1994, MNRAS, 266, 441

Luck, R.E., Gieren, W.P., Andrievsky, S.M., Kovtyukh, Fouqeé, P., Pont, F.& Kienzle, F. 2003, A&A, 401, 939

Macri, L.M., Bothun, G.D., Freedman, W.L., Garnavich, P.M., Jha, S., Madore, B.F. & Richmond, M.W. 2001, ApJ, 559, 243

Piskunov, N.E., Kupka, F.,Ryabchikova, T.A., Weiss W.W & Jeffrey C.S. 1995, A&AS, 112, 525

Sandage, A., Bell, R.A. & Tripicco, M.J. 1999, ApJ, 522, 250

Sasselov, D.D., Beaulieu, J.P., Renault, C., Grison, P., Ferlet, R. et al. 1997, A&A, 324, 471