# Mixing-induced activity in open flows

## Abstract

We develop a theory describing how a convectively unstable active field in an open flow is transformed into absolutely unstable by local mixing. Presenting the mixing region as one with a locally enhanced effective diffusion allows us to find the linear transition point to an unstable global mode analytically. We derive the critical exponent that characterizes weakly nonlinear regimes beyond the instability threshold and compare it with numerical simulations of a full two-dimensional flow problem. The obtained scaling law turns out to be universal as it depends neither on geometry nor on the nature of the mixing region.

###### pacs:

47.54.-r; 47.70.-n; 89.75.Kd## 1 Introduction

In many geophysical and laboratory flows active chemical and biological processes occur. These active processes are crucially dependent on the nature of the flow. Especially important from the theoretical point of view becomes understanding the role of mixing. As we demonstrate below, even localized region of strong mixing introduced in the laminar flow is able to significantly change the overall picture of activity. Because the active processes accompanied by mixing are widespread, the outcomes of our analysis are highly relevant to practically important problems varying from chemical reactions in micro-mixers to plankton growth in ocean, see [1] and references therein.

What turns out to be generic, the active processes quite often occur in an open rather than in a closed geometry. Here the main issue is whether the throughflow is stronger or weaker than the activity. One has to compare the velocity of the throughflow with the velocity of the activity spreading due to diffusion. If the throughflow is stronger, the activity is blown away like a candle flame in a strong wind, in the opposite case a sustained activity can be observed [2, 3, 4]. This simple picture is valid, however, only for homogeneous media. Often additional vortexes are superimposed on a constant throughflow, due to, e.g., wakes behind islands in ocean currents or mixing enforced by revolving fan blades in laboratory experiments. We want to study under which conditions such an additional kinematic mixing in a strong open flow can lead to a transition to a sustained activity, and to characterize this transition quantitatively.

Our main model is a reaction-advection-diffusion equation for the dimensionless concentration of an active scalar field

(1) |

Here is a constant throughflow in -direction, is a molecular diffusion of the scalar field. Activity is assumed to be of the simplest form: a linear growth with rate with a saturation at . The nonlinearity index is typically integer (1 or 2) for chemical reactions, while for biological populations a wide range of values of has been recently reported [5]. Mixing is described by a spatially localized incompressible velocity field , its intensity is denoted as . Note that in the absence of fluid flow Eq. (1) is reduced to the famous Kolmogorov-Petrovsky-Piskunov-Fisher (KPPF) model of an active medium with diffusion (see, e.g., [6] for original references, analysis, and applications of KPPF), while for Eq. (1) describes a linear evolution of a passive scalar in a flow. Model (1) can be used for the description of biological activity, where is, e.g., concentration of a growing plankton, advected by oceanic currents [7]; for a possible laboratory realization see recent experiments with an autocatalytic reaction in a Hele-Shaw cell with a throughflow [8].

In the absence of flow, the diffusion causes the active state to spread forming eventually a front with velocity [9]. Thus, for vanishing mixing , the activity is blown away provided . For this parameter range the instability in Eq. (1) is convective and in the absence of external sources, no activity is observed in the medium. A nontrivial state is, however, possible if there is a localized source of the field : then a growing tail stretches from this source in the downstream direction, where it eventually saturates at . The linearized problem with a point (-function) source of intensity can be readily solved, yielding

(2) |

in one-dimensional setups and

(3) |

in two dimensions (where is the modified Bessel function, ). Note that in both solutions as , where .

In this paper we demonstrate, that beyond some critical intensity , a localized mixing turns the convective instability locally into the absolute one, which results in a stationary (in statistical sense) profile of (see Fig. 1 for a sketch of the profile and Figs. 4 and 7 below for numerical examples). Beyond criticality , the mixing region acts as an effective source of the field, in Fig. 1 this region is denoted as a “source.” A “tail” where the field grows exponentially as in (2), (3) extends downstream of the source. Further downstream stretches the “plateau” domain where . Our main quantitative result, obtained by matching solutions in these three domains, is the critical exponent relating the intensity of the effective source with the mixing intensity: .

To develop the theory we use the concept of global modes [10, 11, 12]. In this concept a self-sustained non-advected pattern arises due to inhomogeneities of the system. Typically, one considers inhomogeneities of the growth rate : if has a hump where locally the front velocity is large , then a global mode appears, located at this hump and downstream. In this paper we are interested in another, mixing-based mechanism of a global mode appearance. It can be easily understood if the concept of effective diffusion (see, e.g., [13]) is used to describe the mixing term in (1). In this approach, we phenomenologically introduce effective diffusivity that accounts for the coarse grained mixing dynamics, and write instead of (1) an equation with a non-homogeneous diffusion

(4) |

A hump of diffusivity leads to an increase of the local front velocity , and one expects that when the front propagation prevails over the throughflow, a stationary global mode can appear, producing a mixing-induced sustained structure. We focus on a geometry shown in Fig. 1(a): in a constant open flow there is a localized region of strong mixing, which, as we will see below, is not necessarily chaotic or turbulent. The theory below will be developed for a one-dimensional case, which is relevant, e.g., for flows in a micropipe; the results will be supported by numerical simulations of two-dimensional flows. We restrict ourselves to this case because of computational simplicity, and also because two-dimensional flows are relevant for many geophysical and laboratory (especially in microfluidics) experimental situations.

## 2 Linear stability analysis

### 2.1 Effective diffusion model

We start with a linear analysis of a one-dimensional situation, described by the linearized at Eq. (4):

(5) |

We look for an exponentially growing in time solution and with an ansatz obtain

(6) |

As we have a homogeneous medium with , here the solution should tend to values at which the r.h.s. of (6) vanishes. More precisely, as we have and as we have . With these two boundary conditions one easily finds the solution of (6) numerically, matching at integrations starting at large from the values . In a particular analytically solvable case of a piecewise-constant diffusivity: for and for , one can perform the integration analytically and obtain the equation for the growth rate :

(7) |

The value corresponds to the onset of global instability, in this case (7) gives the relation between the critical values and :

(8) |

From (8) it follows that as . In other words, there exists a minimal size of the mixing region, so that for even a very strong mixing, with a very large effective diffusion, cannot create a global mode (the same is true for a smooth profile of ; note also that the size of the mixing region is not limited from above). This is in contrast to the situation when the global mode is induced by a local hump of the growth rate (cf. [15]): here one can obtain instability even when is highly localized (a delta-function), a global mode then looks as in (2), (3).

A similar analysis can be performed for a two-dimensional inhomogeneous linear reaction-advection-diffusion equation

(9) |

where we assume for and for . Presenting the solutions in the inner and outer domains as for and for , we obtain equations for whose solutions can be written down as series in Bessel functions and . Matching these series at leads to rather cumbersome expressions in terms of expansion coefficients. As a result, the eigenvalue problem reduces to a matrix equation which we solved numerically. In Fig. 2 we show the critical line (corresponding to the condition ) on the plane for . One can see, that similar to the one-dimensional case, there exists a minimal radius of the higher diffusivity spot, so that for the global mode cannot arise.

### 2.2 Non-uniform flow

Now we discuss the linear stability not in the framework of the effective diffusion model (4), but in the full reaction-advection-diffusion problem as in Eq. (1). After the linearization we arrive at a linear stability problem

(10) |

which is non-stationary if the velocity field is time dependent. For a time-independent the stability is defined by the growth rate as before. For time dependent fields , the proper way to determine the stability is to calculate the largest Lyapunov exponent (LE) . This can be done numerically, as described in [16]. Noteworthy, in this consideration we are not restricted to a deterministic flow, as the LE can be calculated also for a randomly or chaotically time-dependent field .

We first calculate the growth rate for a linearized quasi-one-dimensional reaction-advection-diffusion equation (10) subject to a constant open flow with a superimposed stationary vortex, described by a stream function

(11) |

We fix parameters , and evaluate the growth rate for different molecular diffusion constants and vortex intensities (Fig. 3). Note that the parameter simply shifts , therefore we plot , where is the growth rate for a non-mixed flow with . We see that the mixing-induced enhancement of field growth is mostly pronounced for small diffusion and saturates at .

The role of a mixing vortex can be understood in the following way. In the case of a passive advection-diffusion process and for a relatively small molecular diffusivity, the field is trapped by the vortex. In the limit of vanishing diffusivity the field outside the vortex is blown away, whereas everything inside the vortex is trapped and cannot escape for an infinitely long time, as e.g. in [17]. As a result, a pattern in the form of a cloud arises. Because of diffusion, such pattern is no longer stable and its concentration gradually decays. Hence, under an advection-diffusion process, the cloud exists only for a limited time. However, if we now take into consideration activity, this time is spent by the active field to grow, which is enough to compensate the loss caused by diffusion. Thus, a self-sustained pattern is born, Fig. 4.

Next we calculate the LE for a linearized two-dimensional reaction-advection-diffusion equation (10) with a constant open flow and a superimposed oscillating vortex, described by the stream function

(12) |

Here we fix , and and calculate the LE , see Fig. 5. As before, we plot , where is the LE for a non-mixed flow with . Again, the mixing-induced enhancement of field growth is mostly pronounced for small diffusion. In contrast to the previous case the dependence on is non-monotonic and has a maximum at . This is the mixing strength at which a chaotic saddle [1] in the Lagrangean particle trajectories appears. A further increase of the vortex intensity does not lead, however, to a significant growth of the LE.

## 3 Nonlinear analysis. A universal scaling law

Now we develop a nonlinear theory of the global mode. It is clear that the nonlinear saturation stops the exponential growth of a slightly supercritical linear mode and leads to a nonlinear solution with a finite amplitude. Our aim is to describe the dependence of this amplitude on the deviation from criticality. First we notice, that the very notion of the amplitude is here nontrivial. Indeed, the nonlinear solution looks as in Fig. 1(b) (cf. Fig. 7 below); it saturates to in the downstream direction. However, outside the mixing region the field looks like a solution caused by a localized field source. Thus, we can take the effective intensity of this source , which is proportional to the characteristic field amplitude in the mixing region [see relations (2), (3)], as the order parameter of the transition. The deviation from the criticality we will measure with the growth rate , for which holds in the full model (1) or in model (4).

We will consider the simplest possible setup, namely the nonlinear modification of one-dimensional Eq. (5):

(13) |

We look for a stationary global mode , and rewrite this equation as the system

(14) | |||||

(15) |

We consider this system separately in two spatial domains. The first, linear region, includes the inflow and the mixing domains (“source” in Fig. 1): , where the field remains small. In the second, outflow region , the field further grows (“tail”) and nonlinearly saturates (“plateau”). In the linear region, because of smallness of the field, we can neglect in (14), thus we obtain an equation similar to (6). The only difference is that because we look for a stationary solution, in (14) the term is absent. Near the criticality, where is small, we can consider this term as a perturbation, therefore the solution of (14) in the linear region is close to the solution of Eq. (6); it has the asymptotic as . Due to the perturbation term , at the right border of the linear region deviates from : the deviation is proportional to , and, thus, to . At the field is small and .

Next we consider full equations (14), (15) in the nonlinear region . Here the solution should tend as to the saddle fixed point , . Thus, starting integration from large values of in the negative direction, we have to follow the stable manifold of this saddle and match this solution at with the obtained above. Because the value to be matched is very close to , in the region where the solution approaches we can write to obtain

(16) |

Here, since is small, we have approximated the solution of (15) as . Because linear inhomogeneous Eq. (16) is solved in the negative in direction, the solution follows the slowest exponent: , where .

At the criticality, the region of validity of the exponential solution becomes very large. Thus it is dominant for small deviations from criticality , therefore we can estimate the coordinate at which the field saturates (i.e., we reach the state and ) from the relations above: from it follows . Substituting this in the expression for , we obtain . Now we take into account that , and, because the evolution of in the interval only weakly depends on the criticality, . The final expression for the scaling law of the amplitude of the global mode thus reads

(17) |

The critical index depends only on the nonlinearity index and on the dimensionless velocity :

(18) |

This main result of our letter can be physically interpreted as follows. The exponent is determined solely by the nonlinearity index if the throughflow velocity is much larger than the front velocity ( large). Here the field in the plateau domain (see Fig. 1) is effectively uncoupled from the source, and the saturation of the instability is due to the local nonlinearity at the source. For a small throughflow velocity ( close to one) the plateau state interacts with the source via the tail. Due to this “remote control,” the field at the source is saturated more efficiently than due to nonlinearity, here the exponent is determined solely by the form of the tail, which depends on the velocities ratio .

Below we check formula (18) with direct numerical simulations of model (1). A stationary vortex (12) with and was imposed on a constant flow with . Keeping the diffusion constant fixed , for different field growth rates we have found, from the linearized equations, the critical vortex intensities at which the global mode becomes first unstable. Then we solved full nonlinear equations close to criticality and found the exponent according to (17). The stationary problem was solved with a finite difference method in a domain , with periodic boundary conditions in and conditions , . The results are presented in Fig. 6, they are in very good agreement with the theoretical prediction (17), (18). Figure 7 shows the example of the stationary mode appearing beyond the instability threshold. A similar analysis performed for the quasi-one-dimensional flow (11) also provides very good agreement with formula (18).

## 4 Conclusions

We have described the mixing-induced transition from a convectively unstable active field in an open flow to a persistent global mode. Our theoretical approach bases on the representation of the mixing region as a domain with locally enhanced effective diffusion. For a nonlinear regime close to criticality we have derived the critical exponent (18) that depends only on two parameters of the system: the dimensionless flow velocity normalized by that of the front, and the nonlinearity index . For large velocities the critical exponent depends only on the system’s nonlinearity, which means a local in space saturation of the instability. For small velocities the exponent is a function of velocity, here the growing downstream tail of the active field imposes the saturation. Notably, this prediction of the one-dimensional theory is in a good accordance with two-dimensional calculations. The obtained results are independent of the geometry and the nature of the mixing region and for these reasons are expected in different systems and at different scales. The generic nature of our findings indicates that turbulent mixing can play a key role in open flows that involve active chemical and biological processes.

## References

### Footnotes

- Paper published in Phys. Scr. T132, 014035 (2008);

see also Phys. Rev. Lett. 99, 184503 (2007)

### References

- Tél T, de Moura A, Grebogi C and Károlyi G 2005 Phys. Rep. 413, 91–196
- Kuznetsov S P, Mosekilde E, Dewel G and Borckmans P 1997 J Chem. Phys. 106, 7609–7616
- McGraw P N and Menzinger M 2003 Phys. Rev. E 68, 066122
- Sandulescu M, Hernández-García E, López C and Feudel U 2006 Tellus A 58, 605–615
- Sibly R M, Barker D, Denham M C, Hone J and Pagel M 2005 Science 309, 607–610
- Murray J D 1993 Mathematical Biology (Berlin: Springer)
- Huisman J, Pham Thi N N, Karl D M and Sommeijer B 2006 Nature 439, 322–325
- Leconte M, Martin J, Rakotomalala N and Salin D 2003 Phys. Rev. Lett. 90, 128302
- van Saarloos 2003 Phys. Rep. 386, 29–222
- Chomaz J M, Huerre P and Redekopp L G 1988 Phys. Rev. Lett. 60, 25–28
- Couairon A and Chomaz J-M 1997 Physica D 108, 236–276
- Tobias S M, Proctor M R E and Knobloch E 1998 Physica D 113, 43–72
- Dimotakis P E 2005 Annu. Rev. Fluid Mech. 37, 329–356
- Nelson D R and Shnerb N M 1998 Phys. Rev. E 58, 1383–1403
- Dahmen K A, Nelson D R and Shnerb N M 2000 J. Math. Biol. 41, 1–23
- Pikovsky A and Popovych O 2003 Europhys. Lett. 61, 625–631
- Lyubimov D V, Straube A V and Lyubimova T P 2005 Phys. Fluids 17, 063302