A bistable reaction--diffusion system in a stretching flow

We examine the evolution of a bistable reaction in a one-dimensional stretching flow, as a model for chaotic advection. We derive two reduced systems of ordinary differential equations (ODE's) for the dynamics of the governing advection--reaction--diffusion partial differential equation (PDE), for pulse-like and for plateau-like solutions, based on a non-perturbative approach. This reduction allows us to study the dynamics in two cases: first, close to a saddle--node bifurcation at which a pair of nontrivial steady states are born as the dimensionless reaction rate (Damkohler number) is increased, and, second, for large Damkohler number, far away from the bifurcation. The main aim is to investigate the initial-value problem and to determine when an initial condition subject to chaotic stirring will decay to zero and when it will give rise to a nonzero final state. Comparisons with full PDE simulations show that the reduced pulse model accurately predicts the threshold amplitude for a pulse initial condition to give rise to a nontrivial final steady state, and that the reduced plateau model gives an accurate picture of the dynamics of the system at large Damkohler number.


Introduction
In natural and industrial environments active processes such as chemical or biological reactions are often embedded in a fluid flow. Examples include mixing of reactants within continuously fed or batch reactors [1,2], the development of plankton blooms and occurrence of plankton patchiness [3][4][5][6], increased depletion of ozone caused by chlorine filaments [7] and flame filamental structures in combustion systems [8]. The fluid flow is very often time-dependent and stirring. This gives rise to interesting effects, and the presence of compression and stirring in chaotic flows typically leads to filamentation of the chemical or biological components. This evidently changes the reaction dynamics. The filamentation caused by the chaotic advection leads to an increased surface area of the reacting components and an increase of the reaction output. However, if the stirring rate is too strong, the filaments become thinner and thinner and the reaction may stop. The corresponding flow-mediated saddle-node bifurcation has been reported in [9][10][11][12][13][14][15][16].
A natural question in such stirred reactions is: Given a chaotic flow environment and an initial distribution of chemical or biological reactants, is it possible to predict whether the reaction will take place and develop or whether it will die out? This question is of paramount importance in an environmental context where the reactants may be ozone in the atmosphere or pollutants in the ocean, or in an industrial context where one is interested in maximizing the reaction output. This problem is well-studied in the non-stirred case [17]. However, so far the complex nature of chaotic flows has prohibited an analytical treatment of this problem in the stirred case.
In principle, these phenomena can be studied directly in a two-or three-dimensional reactionadvection-diffusion system, albeit with huge computational effort. An analytical treatment of the full system is prohibited by the complicated nature of the underlying equations, which involve multiple-scale processes. Simplified models are needed to capture essential features of the influence of the stirring on the reaction kinetics. In filamental or lamellar models [18], the two-dimensional problem of reacting tracers is replaced by a one-dimensional problem of the for n reacting tracers (u i , i = 1, . . . , n) with diffusion coefficients D i , reaction rates k i and stir-ring rate λ. Such models have been applied by several groups, to autocatalytic, bistable and excitable media in several physical, chemical and biological contexts [4,5,[9][10][11][12][13][14]. Phenomenological filamental models such as (1) can be justified by the following consideration: The chaotic advection causes filaments to be stretched in one direction and compressed in another. In the stretched direction, the concentration is homogenized and gradients along the filaments can be neglected. This motivates a one-dimensional reduction for the concentration in the direction transverse to the filament, subject to the effect of stirring and compression. The parameter λ can be thought of as the Lagrangian mean strain in the contracting direction, and may be argued to be given by the absolute value of the negative Lyapunov exponent or the (slightly larger) topological entropy. For a different approach to this problem see [19,20]. The validity of such simplified models has been numerically investigated in [21].
The solutions of this model are found to be generically either of a plateau-like or a pulse-like character. Plateau solutions are found stably when the stirring rate is low compared to the reaction rate. Pulse-like solutions are found close to the saddle-node bifurcation, where the stirring rate is strong enough to suppress the reaction; moreover the unstable solution for weak stirring rates is also of a pulse-like character. The unstable solution will be of interest when we consider the initial-value problem and ask when an initial perturbation will develop into a plateau-like solution. We extend a non-perturbative approach developed in [26] to derive from (1) a corresponding low-dimensional system of ordinary differential equations which describes the time evolution of plateau-like and pulse-like solutions. We determine the equilibrium solutions, accurately describe the saddle-node behaviour and derive an analytical expression for the asymptotic width of a stationary front. By studying the phase-portrait of the reduced system we are able to determine the fate of various initial concentration profiles in (1); these predictions are then verified numerically by comparison with simulations of the full partial differential equation (1).
In Section 2 we introduce the bistable model and discuss its solutions and their bifurcations. In Section 3 we describe the non-perturbative variational approach. This method will be used in Section 4 to describe pulse-like solutions and in Section 5 to describe plateau-like solutions. In Section 6 we use the reduced system derived in Sections 4 and 5 to characterise the initial-value problem. The paper concludes with a discussion and an outlook for further research.

The model
In this paper we study a bistable model where 0 < α < 1. The Damköhler number Da = k/λ measures the ratio of the time scales of fluid motion and reaction: small Damköhler numbers correspond to fast stirring/slow reaction, while for large Damköhler numbers the system behaves asymptotically like an unstirred system.
The ODE corresponding to (2) in the absence of spatial effects has two stable fixed points, u = 0 and u = 1, which are separated by an unstable fixed point at u = α. For the unstirred case of the PDE (2), the system is well known and well described in textbooks such as [22][23][24][25]: an initial perturbation which is larger than α over a finite range will spread over the whole domain if 0 < α < 0.5; by contrast, if 0.5 < α < 1, an initial perturbation will decay to the stable state u = 0. The stirred case is much less well understood, and the initialvalue problem has to our knowledge never been analytically treated. The stirred case was investigated numerically in [5,10], and semi-analytically in [15]: stationary fronts between the u = 0 and u = 1 states exist as a balance between the x-dependent stirring and the diffusionmediated counterpropagating fronts, for large enough values of the Damköhler number. An initial sufficiently large perturbation seeded at x = 0 spreads as a pair of fronts, driven by its reaction kinetics and diffusion, until the fronts reach the locations x = ±ν where their velocities equal the ambient spatially dependent velocity of the chaotic stirring.
It has been observed in [10,15] that there is a critical Damköhler number such that no stationary pulses exist for Da < Da c , i.e., when the time scale of the chaotic advection, τ f = 1/λ, becomes too short in comparison with the time scale of the reaction, t r = 1/k. For large Damköhler numbers, an asymptotic expression for the scaling of the total concentration was developed in [10]. However, the techniques used in [10] cannot describe the behaviour close to the bifurcation point Da = Da c . Moreover, these techniques are stationary and cannot answer questions about the evolution of initial perturbations of reactants.
We now describe in more detail the asymptotic form of the two nontrivial steady states for large Da. u(x) x Figure 1: The two nontrivial steady states of (2) for α = 0.2 and Da = 100. The 'pulse' solution is unstable; the 'plateau' solution is stable.

Asymptotic steady states for large Da
We aim to find steady states u(x) that satisfy when Da 1. Numerical simulations show that for Da 1 the stable solution is plateau-like and the unstable solution is bell-shaped (see Figure 1).
For the unstable bell-shaped solution u = V (x), we introduce with an error of order Da −1 . We may now readily determine the appropriate leading-order solution, which gives rise to We note that this gives which we shall later compare with a corresponding result obtained from a reduced test function ansatz using a bell-shaped test function, derived in Section 4.
In the large-Da limit, the plateau solution u = U(x) is, to leading order in Da, given by U = 1 in a region −ν < x < ν and U = 0 outside this region, where ν = √ Da ω, for some ω = O(1).
The locations of the interfaces are determined by a balance between the tendency of the pulse to spread (with constant velocity) and the compressive effects of the imposed velocity field (for which the velocity is proportional to x) [10,15]; we return to this point in more detail in Then from (3) we obtain with an error of order 1/ √ Da. Correspondingly, with ω = D/2(1 − 2α). There is a similar transition region near x = −ν. We shall in Section 5 use tanh profiles such as these as an ansatz for our second reduced model; since the ansatz captures exactly the large-Da form of the fronts, we shall find correspondingly that the reduced model provides an excellent approximation to the dynamics of the full PDE (2). Figure 1 shows the solutions of (2) obtained by a shooting algorithm for high Damköhler number (Da = 100). It clearly shows that the stable solution is plateau-like and the unstable solution is pulse-like. Note that the stable solution for smaller Da, close to the bifurcation (not shown here), is also pulse-like.

Asymptotic solutions for small Da
We now seek an appropriate ansatz for the solution behaviour close to the saddle-node bifurcation. Some insight into the form of the solution may be gained by examining the small Damköhler number limit of the full time-dependent problem (2). Beyond the saddle-node bifurcation point, i.e., for Da < Da c , no steady-state solutions exist, and any initial condition decays to zero. In fact, if we expand the solution as a series in the Damköhler number according , the equation is satisfied at leading order by where w 0 → 1/ √ 2D and f 0 → 0 as t → ∞; the latter limit reflects the nonexistence of nontrivial steady states for Da < Da c . Although Da c = O(10), which is clearly not small, it seems that the leading-order term (7) provides an accurate representation of the solution to (2) close to the saddle-node point Da = Da c -see Figure 2. So in the first instance we shall approximate both stable and unstable solutions close to the saddle-node point by a Gaussian; we shall later consider, in somewhat less detail, other possible profiles of 'bell-shaped' form.
Using a non-perturbative method developed in [26] we now derive a set of ordinary differential equations which describe the dynamics of these plateau-and pulse-like solutions. In the next section we briefly explain this variational approach.

Nonperturbative, variational method
A method was developed in [26] to study critical wave propagation of single pulses and pulse trains in excitable media in one and two dimensions. It was based on the observation that close to the bifurcation point the pulse shape is approximately a bell-shaped function. Numerical simulations and the asymptotic analysis in Section 2.1 and 2.2 show that this is also the case for the bistable model (2) close to the saddle-node bifurcation at which the stable and unstable solutions are born, and on the unstable branch at large Damköhler numbers. A testfunction approximation that optimises the two free parameters of a bell-shaped function, i.e., its amplitude and its width, allows us to find the actual bifurcation point, Da c , and determine the pulse shape for close-to-critical pulses at Damköhler numbers near Da c . We note that the framework of asymptotic techniques, such as inner and outer expansions where the solution is separated into a steep narrow front and a flat plateau, are bound to fail close to the bifurcation point, because the pulse is clearly bell-shaped, and such a separation is no longer possible. We shall make explicit use of the shape of the pulse close to the critical point and parameterise the pulse appropriately.
We choose u(x, t) of the general form where U(η) is a symmetric, bell-shaped function of unit width and height, and f (t) is the amplitude of the pulse. Close to the saddle-node bifurcation and for the unstable solution at large Damköhler numbers, the solution is indeed a bell-shaped pulse, which we approximate by a Gaussian (see Figure 2). However, our results do not depend on the specific choice of the test function, and the numerical values differ only marginally when sech-functions are used instead, for instance (see Figure 5). We restrict the solutions to a subspace of a bell-shaped function f U(η), which is parameterised by the amplitude, f (t), and the inverse pulse width, w(t). The evolution of these parameters is then determined by minimizing the error made by the restriction to the subspace defined by (8); this is achieved by projecting (2) onto the tangent space of the restricted subspace, which is spanned by ∂u/∂f = U and ∂u/∂w = f ηU η /w. We set to zero the integral of the product of (2) with the basis functions of the tangent space (over the entire η-domain). This leads to ordinary differential equations for f (t) and w(t); it also yields an approximation to the critical Damköhler number Da c .
One may in fact choose any plausible test function and optimize the restriction to the particular space of solutions by varying the parameters. This variational idea is not restricted to pulselike solutions. Far away from the saddle-node bifurcation, at large Damköhler numbers, where the stable asymptotic solution comprises a pair of fronts (see Section 2.1) we may use as a test function a combination of tanh-functions where a(t) = w(t)ν(t) and the fronts are located around x = ±ν. In a manner analogous to that described above for a pulse ansatz, the variational approach allows us in this case to determine the time evolution of the solution amplitude f (t), the inverse interface width w(t) and the front locations ±ν(t) by projecting (2) onto the tangent space of the restricted solution subspace, which is now spanned by ∂u/∂f , ∂u/∂w and ∂u/∂ν.

Stable and unstable pulse solutions
Near the turning point at Da = Da c (α), both nontrivial steady-state solutions take the form of bell-shaped pulses near x = 0. We note that in the literature a pulse commonly refers to a stable homoclinic solution; however, here we use the term merely to mean a bell-shaped function, regardless of whether it is stable or unstable. This terminology allows us to distinguish between the two test function ansätze (8) and (9).
To determine a reduced model for such pulse solutions, we therefore introduce the ansatz (8) and to determine the evolution of f (t) and w(t) we enforce the vanishing of the two inner products For a Gaussian pulse, with U(η) = e −η 2 , we note the requisite integrals Thus after some algebra we find that f (t) and w(t) evolve according tȯ In order to examine the dynamics of this system in the phase plane it is useful for brevity of notation to write it in the forṁ and note that the µ i and λ i are all positive; also µ 2 = λ 2 (this equality does not necessarily hold for other profiles U(η)).
We now examine the steady states (w, f ) = (w, f) of (12), (13) and their stability. Note that the steady-state problem has been considered previously [15], although only the solutions identified below as 'P 4 ' were considered there. Although not all steady states of (12), (13) correspond to distinct physical solutions, they help us to map out the structure of the phase plane, and hence determine the dynamics of (12), (13).
This steady state exists for all parameter values, and a linearisation about P 1 yields the eigenvalues λ 1 and −µ 1 . Therefore P 1 is a saddle point, unstable to perturbations in w and stable to perturbations in f .
This steady state exists for all parameter values. A linearisation about this point yields the eigenvalues −µ 1 − µ 2 λ 1 /λ 2 < 0 and −2λ 1 < 0. Thus P 2 is a stable node. Here This quadratic has two real, positive solutions for all parameter values; indeed f − and f + are independent of Da. Note that, since P 3 corresponds to infinitely wide pulses (w = 0), the roots of the quadratic would be α and 1 if there were no approximation due to the Gaussian pulse ansatz. A linearisation about either steady state yields the eigenvalues  Here f satisfies We note the coefficients (cf. [15]) The quadratic (15) has two real, positive roots provided and Da > Da pulse Note that α m is slightly smaller than 0.5, the threshold for existence of non-zero solutions in the unstirred case.
In this case, we denote the two solutions of (15) by f * Since the right-hand side of (17) must be non-negative, we deduce, by comparison with (14), A linearisation about either of the P 4 steady states yields, for f = f + δf (t) and w = w + δw(t), Hence the growth rate σ satisfies Note that for f * +   we show that good agreement is also given for different ansätze such as U(η) = sech n η. We thus conclude that the actual form of the pulse shape is not particularly important (and hence we cannot expect to deduce this pulse shape from any asymptotic theory).
By contrast, the asymptotic approaches which are usually employed, and which assume plateaulike solutions with well-separated constant plateaus, fail to describe the bifurcation behaviour.
We note that the agreement between the pulse model and the full problem (2) is excellent on the lower branch. Figure 6 shows the excellent agreement between u(0) = max u(x) computed according to the pulse result (15) in the limit Da → ∞ and a direct large-Da asymptotic analysis of (2), given by (5).

Plateau solutions
Away from the saddle-node bifurcation, the stable solution takes the form of a nearly uniform region around x = 0, surrounded by a pair of fronts. This solution may be captured by writing where φ(η) = tanh(η + a(t)) − tanh(η − a(t)).
Then the evolution of f (t), w(t) and ν(t) = a/w (or, equivalently, f , w and a) may be determined by forcing the inner products of (2) with each of u f , u w and u ν to be zero, where with ψ(η) = sech 2 (η + a) + sech 2 (η − a).
The corresponding steady-state calculation has been considered in [15]; here we extend their analysis to the time-dependent problem. This permits us to determine which initial conditions lead to excitation and which to extinction. Details of our derivation, including analytical expressions for all the requisite inner products, are given in the Appendix: we find that the equations for the time evolution of f , w and ν take the form The various coefficients λ i , µ i and β i are complicated functions of a and α, and are listed in the Appendix.

Steady states
We shall not attempt here a full analysis of the three-dimensional phase space for f , w and ν; instead we consider only nontrivial steady-state solutions for which f , w and ν are all nonzero.
We expect two such solutions when Da exceeds some threshold, one stable and one unstable.
The equations for the (nontrivial) steady states may be written in the form These equations are readily solved without recourse to nonlinear root finding, for example in the following way. First we fix values for α and a. Then we use (23) to write Then, upon substitution of this expression for w 2 in (24), we have Likewise Now by eliminating Da between (27) and (28), we find the following equation for f : which is of the form The quadratic equation (30) needs careful interpretation. First we note that the coefficients The asymptotic results are indistinguishable from the numerical results for f and a.
(26) the corresponding values of w, we find numerically that the root of (29) corresponding to positive Da seems always to give rise to w 2 > 0, and so is indeed physically reasonable. Thus we emphasise that both branches of solutions evident in Figures 7 and 8 arise from one of the roots of (30) for fixed α, as a is increased monotonically.

Large Damköhler number, Da 1, upper branch
We now construct the plateau solutions on the upper branch, in the limit Da 1. In that limit, we find Recall that a = O(Da). The result for β 5 indicates that this quantity is exponentially small in Da. Thus at leading order in Da we have, from (23)- (25), 2 ) = 0.
The first of these equations givesf = 1 orf = α; the latter solution must be rejected since the last equation would then giveā = −(2 − α)α/4 < 0. The choicef = 1 is consistent with our numerical calculations, which indicate that this is the solution relevant to the upper branch.
Withf = 1, the remaining equations give Dw 2 = 1 8 andā = 1 4 (1 − 2α). Thus in the large-Da limit, Note, as a consequence, that for large Da we have in agreement with the phenomenological argument given in [10,15]. This formula can be understood if we note that the location of a stationary front is given approximately through a balance between the front velocity v 0 = √ 2DDa 1 2 − α of the unstirred problem [10] and the velocity of the chaotic stirring (i.e., x in the scaling of (2)). Thus in the stirred problem, the front has zero velocity when v 0 = x, which implies v 0 = ν, and (33) follows.
Finally, the solution takes the form A comparison between this expression and a numerical calculation of the stable steady state of (2) is shown in Figure 9; to graphical resolution, the two graphs are indistinguishable. Figure 10 shows how well the ansatz (19) captures the width of the steady-state solution.
Plotted is 2(ν approx − ν exact ), where ν approx and ν exact are, respectively, the half-widths of the solutions to the plateau model and the full problem (2). In each case only the stable solution on the upper branch is considered, and the 'half-width' is defined as ν approx , ν exact = x + > 0, where u(x + ) = 0.5. The difference 2(ν approx −ν exact ) clearly decreases with increasing Damköhler number, which is expected, because the tanh-front ansatz gets closer to the actual solution, which asymptotically is a pair of tanh-fronts (see Section 2.1).
Solutions of the plateau model on the lower branch are of less interest because they do not correspond closely to steady-state solutions of the full PDE, where the unstable solution is pulse-like.

Stability
Suppose that (f, w, ν) is a nontrivial steady state of (20)-(22) (i.e., f, w, ν = 0). Since the various coefficients in (20)- (22) are functions of a, it is useful to consider the linearised evolution of δf , δw and δa, where in the linearisation δa = wδν + νδw. Thus we have where ∂λ i ≡ ∂λ i ∂a and where all nonperturbation quantities are evaluated at the fixed point. There are two further equations, with λ i → µ i and λ i → β i . The three may be written together as where the constant-coefficient matrices A and B are readily determined from the three perturbation evolution equations.

Large Damköhler number, Da 1, upper branch
On the upper solution branch, in the limit Da 1, we have where the zero entries indicate terms that are exponentially small in Da. Thus, in view of the scalings (31), It then follows (after some algebra) that We require the dominant eigenmodes of A −1 B. There are three eigenvalues: The eigenvector corresponding to Λ 1 has δf = δw = 0 and δa = 1, say; thus disturbances in a decay like e −t . We find Λ 2 ∼ −(1 − α)Da; recalling the scalings f = O(1), w = O( √ Da) and , we find that the associated eigenvector has (δf, δw, δa) ∼ (F, W Da 1/2 , ADa), The final eigenvalue is Note that this calculation confirms the stability of the upper branch solution, since Λ 1,2,3 < 0. Since 0 < α < 1 2 and since 3 2 (π 2 − 6) −1 ≈ 0.3876, we have 6 Comparison between dynamics of the reduced models and the full PDE Our models extend the steady-state analysis of [15] by providing approximations to the dynamics of the PDE (2). This allows us to examine the evolution of initial conditions which are either plateau-like or pulse-like. For instance, using the Gaussian pulse ansatz, the dynamics is approximately reduced to that of the two-dimensional system (12)- (13). Given an  Figure 11. The agreement is excellent. There is similar agreement at higher Da; for example, at Da = 100, and w(0) = 1, the Gaussian ansatz gives a threshold f (0) = 0.2537, while from the PDE we find 0.2525 ± 0.0005.
Assessing the usefulness of the plateau model is rather harder. This is because if we start with an initial profile of the form (19), with well separated fronts, and u(x, 0) approximately uniform between those fronts, then in the uniform region the initial evolution of u is well approximated The agreement between the plateau model and the full PDE is remarkable because, as shown in Section 2.1, the unstable solution is pulse-like and not plateau-like. Furthermore, as shown in Figure 4, the unstable solution predicted by the plateau model is quantitatively incorrect.
However, some insight into the agreement may be obtained by noting that the unstable pulselike steady-state has eigenvalues σ = O(Da), which are large in the limit Da → ∞, and which indicate that the associated instability is reaction-dominated. Our simulations of the PDE (2) confirm that the initial condition illustrated in Figure 12 first rapidly undergoes an adjustment of the central concentration to u ≈ 1, consistent with a time scale O(1/Da), then a slower relaxation of the front width and location.

Conclusions
We have developed reduced systems of ordinary differential equations describing the time evolution of both plateau-like and pulse-like solutions of a one-dimensional bistable reactiondiffusion partial differential equation for a stirred flow environment. The PDE possesses a trivial steady state in which all the reactant is used up, and, when the dimensionless reaction rate Da is great enough, a pair of nontrivial steady states, one stable and one unstable, born in a saddle-node bifurcation at Da = Da c . The pulse model was shown to provide a good approximation close to the bifurcation point and for the unstable steady solution at large Da.
The plateau model applies for large Da, and captures well the form of the stable steady state.
By studying equilibrium pulse solutions, we were able to accurately describe the bifurcation behaviour of the system near Da c . For large Damköhler numbers our reduced plateau model allowed us to describe the stationary fronts, and derive a commonly used phenomenological formula for the width of a stationary front.
The time-dependent ODE models allowed us to study the initial-value problem for (2). We found that although the unstable solution is actually of a pulse-like shape, the dynamics of (2) is well captured by the plateau model for high Damköhler numbers. This remarkable and nonintuitive result is due to the fact that the growth rate of the instability of the unstable solution is O(Da), whereas the fastest growth rate for the shape-change of a front is O(1). Our numerical solutions of (2) seem rapidly to approach the plateau form, and thus the dynamics of the PDE (2) is described accurately by the plateau model.
Besides being instructive in themselves, one-dimensional models such as (1) aim to provide information about the evolution of chemical reactions in a two-or three-dimensional chaotic flow. It is now known [21] that simple variants of reduced models perform poorly in predicting the yield of multi-stage reactions, and so a more sophisticated analysis, as presented here for the bistable case, is therefore warranted. One particularly promising extension [27] of onedimensional models such as (1) involves two stages. In the first stage, a point x in the twoor three-dimensional problem and a time t are selected, at which the chemical composition is required. This point and its associated stable manifold are then advected backwards in time, to t = 0, and the initial condition for the system is imagined to be sampled by the stable manifold (which is now exponentially stretched and highly contorted). The second stage involves carrying out a one-dimensional simulation of the advection, reaction and diffusion along the evolving stable manifold (now forwards in time), using a model such as (1), but with a nonconstant compression rate [28].