Persistent localized states for a chaotically mixed bistable reaction

We describe the evolution of a bistable chemical reaction in a closed two-dimensional chaotic laminar flow, from a localized initial disturbance. When the fluid mixing is sufficiently slow, the disturbance may spread and eventually occupy the entire fluid domain. By contrast, rapid mixing tends to dilute the initial state and so extinguish the disturbance. Such a dichotomy is well known. However, we report here a hitherto apparently unremarked intermediate case, a persistent highly localized disturbance. Such a localized state arises when the Damkoehler number is great enough to sustain a hotspot, but not so great as to lead to global spread. We show that such a disturbance is located in the neighborhood of an unstable periodic orbit of the flow, and we describe some limited aspects of its behavior using a reduced, lamellar model.


I. INTRODUCTION
The progress of chemical or biological processes mixed by fluid flows is of great environmental and industrial significance, in a wide variety of contexts (see [1,2,3,4] for some broad applications; in this paper we shall, for brevity, mostly use the language of chemical processes in our discussions). The underlying processes often possess multiple possible long-time states, characterized by whether some reaction proceeds or dies out. The ultimate fate of a localized initial disturbance can depend subtly on the combined influences of advection by the fluid, diffusion of the reactants, and the reaction kinetics.
This paper concerns the influence of chaotic, laminar fluid mixing on the course of a bistable chemical reaction in a closed flow domain. In a spatially homogeneous environment, such a reaction tends to one of two steady states (which we describe below as the "quenched" and "excited" states) according to whether or not the initial reactant concentrations exceed some threshold. However, the large-time behavior of the reaction from a spatially localized disturbance subject to fluid mixing is a more subtle question, since the mixing has two competing effects on the course of the reaction. On the one hand, the mixing promotes the chemical reaction, since, according to the kinematics of the fluid motion, interfaces are exponentially stretched. However, the mixing also tends to suppress the reaction, since the repeated stretching and folding that is characteristic of chaotic flows generates a complicated filamentary pattern; if the thinning of the filaments is sufficiently rapid, it may ultimately quench the reaction, through diffusive dilution of the reactants. The loss of a stable excited state as the mixing grows in strength (or as the reaction rate is reduced) may be * Electronic address: stephen.cox@nottingham.ac.uk; URL: http://www.maths.nottingham.ac.uk/personal/smc viewed as a saddle-node bifurcation, and has been identified in a number of reaction schemes [5,6,7,8,9,10,11]. For example, Neufeld et al. [6] note that for both bistable and autocatalytic reactions there are two regimes: when mixing is slow, localized excited perturbations propagate throughout a closed domain; by contrast, under rapid mixing, localized perturbations decay to the uniform quenched state. In the propagative case, perturbations spread in the form of filaments. For an excitable medium, Neufeld et al. [7] describe three scenarios for a closed flow: in the first, there is eventually a global excitation, which is spatially incoherent, so different parts of the domain are excited at different times; in the second, there is again a global excitation, but this time it is spatially coherent, and ultimately the system homogeneously decays to the unexcited state; in the third, the excitation enjoys no significant propagation before its ultimate decay. Kiss et al. [12] investigate a chemical reaction scheme capable of supporting both steady and oscillatory homogeneous states, subject to chaotic mixing. They find that ultimately the system reaches a spatially homogeneous state, which may be of either steady or oscillatory type, depending on the parameters.
We shall show here the existence of a third regimebeyond the global quenched and excited states -comprising a persistent localized disturbance. For definiteness, we consider the model bistable reaction scheme in which a chemical concentration A(x, t) satisfies, in dimensionless form, where 0 < α < 1/2 [6,7,13]. Here Pe = U L/D is the Péclet number, Da = kC 2 /U L is the Damköhler number and u(x, t) is the dimensionless fluid velocity (here U , L and C are scales for the fluid velocity, lengths and chemical concentrations, respectively; D is the diffusion coefficient and k is the reaction constant). In view of the need to carry out well resolved simulations at moderately high Péclet number, we shall consider a two-dimensional problem, in which u is the so-called "sine flow" [14]. This choice of flow permits a particularly efficient and accurate spectral numerical solution of the problem [15,16]. Let us now discuss the behavior of (1). First consider the reaction in the spatially homogeneous case. Then if the initial value of A exceeds the threshold α, A → 1 at large times. By contrast, if A is initially less than the threshold, then A → 0 at large times. We expect a similar scenario to hold in (1) if the initial state is nonhomogeneous, but the mixing is sufficiently rapid: an initial state is first rapidly averaged in space, then the long-time behavior of the system depends on whether or not, in this averaged (nearly spatially homogeneous) state, A exceeds α. In the opposite limit of no fluid motion, an initial localized disturbance spreads or becomes extinct depending on the initial state, but if the disturbance takes the form A = A 0 (a constant) in some finite region Ω, and A = 0 outside Ω, then the condition for the disturbance to spread is roughly that A 0 should exceed α. We expect a similar scenario for slow fluid mixing. These considerations, and the results of Refs. 6, 7, 9, 12, suggest that, in a closed flow, the fate of a localized disturbance in (1) is as follows: the system tends to a uniform state at large time, with A ≡ 0 (for sufficiently rapid mixing or slow reaction or weak initial perturbation) or A ≡ 1 (for sufficiently slow mixing or rapid reaction). Furthermore, given an initial disturbance that leads to the latter globally excited state for given values of Da and Pe, one would expect that a sufficient reduction in Da would cause the reaction to proceed instead to the quenched state. However, we show here that there is an additional, intermediate scenario, in which the perturbation appears to remain localized for all time, spread along a segment of the unstable manifold of an unstable periodic orbit of the flow. We investigate this localized state through numerical simulation of (1), in particular exploring its dependence on Damköhler and Péclet numbers, and on the initial state. Finally, we shed some limited analytical light by consideration of a reduced, lamellar model that represents the balance between spreading of the hotspot due to reaction and diffusion, and compressive effects of the flow.

II. SIMULATIONS
In this section we describe some numerical simulations of (1). Stirring of the fluid is accomplished by the "sine flow", which is time-periodic with period T ; then [14] u = (sin 2πy, 0) for mT ≤ t < (m + 1 2 )T, (0, sin 2πx) for (m + 1 2 )T ≤ t < (m + 1)T, (2) for m = 0, 1, 2, . . .. The simulations take place in a doubly periodic square box 0 ≤ x, y ≤ 1. We consider the "globally chaotic" case T = 1.6, for which the Poincaré section appears to the eye to have no significant regular islands. The initial condition is given by the Gaussian concentration profile for 0 ≤ x, y ≤ 1, so that the initial spatially averaged concentration is A = 0.01 × Q at t = 0. In most simulations, Q = 1, although we also briefly consider other values of Q, in the range 0 < Q ≤ 2. The choice to center the initial disturbance on the point (0.25, 0.25) was arbitrary, but, as we shall see, it proves serendipitous.
We simulate (1) subject to the flow field (2) and the initial state (3) using a pseudospectral code, which is described in Ref. 16 (see also Ref. 15). The Péclet number in our simulations is limited to 10 5 , to ensure well resolved results with our available computational resources. Most simulations up to Pe = 10 4 are carried out using 512 Fourier modes in each of the x and y directions, and a time step of 2 × 10 −3 ; unlike the simulations reported elsewhere [16], which exclusively use the second-order time-stepping scheme ETDRK2, we mostly use here the fourth-order scheme ETDRK4 (see Ref. 17 for details of these two schemes).
In all the numerical simulations below, the sine flow is fixed to have period T = 1.6, in order to avoid the appearance and disappearance of regular islands with varying T . Thus our investigation of the effects of fluid mixing is carried out by varying the Péclet number (smaller and larger values of this parameter corresponding, respectively, to more rapid and slower mixing). We also investigate the influence of the reaction rate and the initial amount of A upon the long-time behavior of the system through varying the Damköhler number Da and Q, respectively.

A. Numerical results
In all our simulations, we set α = 0.25. Thus for sufficiently rapid mixing the reaction is quenched (since the initial averaged concentration A = 0.01Q < 0.25), and for slow enough mixing the disturbance spreads throughout the domain (since the peak initial concentration  (2), plotted at times T , 2T , . . . , 6T (recall that the line x = 0 is identified with the opposite edge x = 1). Also shown are the contours A(x, y, t) = 0.5 at corresponding times, for a hotspot found at Da = 6 and Pe = 10 5 , indicating that the hotspot is located in the vicinity of the period-6 orbit.
We first set Pe = 10 4 and Q = 1 and consider the effect of varying Da. Fig. 1 shows some illustrative results. For Da = 6, the reaction is too slow to sustain itself and A → 0 everywhere as t → ∞. For Da = 12, the reaction is fast enough to overcome the diluting influence of the fluid mixing; in this case A → 1 everywhere as t → ∞. The choice Da = 8 provides an intermediate case, for which the reaction proceeds to neither homogeneous state as t → ∞, but rather generates a persistent time-dependent state in which A = 0 in most of the flow domain, but A ≈ 1 in a small time-dependent region. This region turns out, upon further examination, to lie in the vicinity of a period-6 orbit of the underlying flow field u. The localized state appears to be time-periodic, with a period of 6T , while the spatial mean concentration A appears to have period 3T . We have tracked this solution for many tens of periods, and it is robust. A further simulation, at Da = 9, at first appears to produce another (somewhat larger) localized "hotspot", but the hotspot does not approach a time-periodic state; instead it gradually grows, eventually filling the entire flow domain, so that A → 1.
The localization on an unstable periodic orbit is readily understood in qualitative terms, as follows. In the vicinity of such an orbit there are directions of local compression and stretching. In the compressive direction, a balance is maintained between the spreading of the reaction-diffusion front and the compressive velocity field [6,7,13], resulting in a finite "width" for the hotspot. The hotspot is also stretched along the unstable manifold of the periodic orbit; growth in its length seems to be moderated by the inability of the reaction to sustain a thin filamentary structure in the vigorous fluid mixing.
We have located this unstable period-6 orbit, using the stabilization algorithm of Pingel, Schmelcher and Diakonos [18] (see  We compute that the Jacobian matrix associated with the period-6 orbit has eigenvalues 18.73 and 0.053. In the next section we shall use this information regarding the stretch rate close to the periodic orbit in the development of a reduced model for the hotspot. Fig. 3 summarizes the results of a large set of simulations, all at Q = 1, for different Damköhler and Péclet numbers. Since the goal of the simulations was to determine the boundaries in parameter space between the three different long-time behaviors, we have not attempted to fill in data where the behavior is (presumably) predictable. (It should be noted that each simulation currently takes on the order of a day or two to run on a desktop PC.) Gaps in the data towards the lower right-hand part of the figure reflect the increased length of runs required there for the long-time behavior to become apparent.
Consider first the parameter dependence to the left of Fig. 3: if we start in the region of hotspots and decrease either Da or Pe then we move into the A → 0 region. This is because either the reaction rate becomes too low to sustain the hotspot, or the initial rate of spreading of the chemical becomes too great for an adequate level of material to be delivered to the period-6 orbit in the first instance. Hence the system is quenched. To the right of the figure, if we start in the hotspot region and increase Da then we move into the A → 1 region, because the reaction rate becomes too great for the hotspot to remain localized. Similarly, if we instead reduce Pe, then again the hotspot cannot remain localized, due to its enhanced rate of spreading through diffusion, and the system is globally excited. The parameter dependence of the system towards the lower right-hand side of the data in Fig. 3 is rather subtle, since a reduction in Pe may lead either to the quenched state or to the excited state. Indeed, we find that our numerical simulations in this region of parameter space often involve a transient, slowly evolving hotspot before eventually selecting the ultimate quenched or excited state. We next discuss the influence of the initial amount of A upon the long-time behavior of the system. To this end, we have carried out simulations with Pe = 10 4 and Da = 8 for a range of values of Q. As Fig. 4 shows, Q ≥ 1.4 leads to the globally excited state, Q ≤ 0.38 leads to the globally quenched state, and intermediate values of Q lead to a hotspot at long times. Although it is not apparent from the figure, which focuses on relatively short-time results, the long-time hotspot solution is the same for all of these simulations with Q in the appropriate range (this is evident in a sample of longer simulations that we have carried out, but which are not reported here). In other words, for these simulations the ultimate hotspot solution depends only on the system parameters (i.e., T , Pe and Da); the role of the initial condition is merely to determine which of the three longtime states is selected.
We now turn to the Péclet number dependence of the hotspot, characterizing the rate of fluid mixing. In Fig. 5, we fix the reaction rate Da = 6 and Q = 1 and vary Pe. Here we find that at Pe = 10 4 , as we have described above, the relatively rapid mixing quenches the reaction, so that A → 0 everywhere. At Pe = 10 5 /3 and at Pe = 10 5 , by contrast, the mixing is slower, and we again find a persistent hotspot. If we accept the qualitative picture above, that the width of the hotspot is determined by a balance between compression and reaction- If the qualitative argument holds, then one should expect that as Pe is increased beyond 10 5 , the hotspot persists, with A continuing to scale with Pe −1/2 . Unfortunately we are unable to simulate reliably beyond Pe = 10 5 , so we are unable to confirm this prediction.

III. LAMELLAR MODEL
There has recently [5,6,7,8,9,10,11,13,19,20] been significant progress in understanding the dynamics of reaction-advection-diffusion problems such as (1) by analysing reduced "filamental" or "lamellar" models [21], which concern the simpler problem of the evolution of the reaction in a one-dimensional setting, in our case Here ξ represents a coordinate in the compressive direction near an unstable periodic orbit of the sine flow. This is the "Lagrangian filament slice model", which has been thoroughly reviewed by Tél et al. [22], largely, but not entirely, in the context of open rather than closed flows. This model might be expected to provide at least a partial picture of the behavior of the hotspot. To see this, we first note that the hotspot is a thin structure, which appears, from the two-dimensional numerical simulations of Sec. II A, to be approximately concentrated around the unstable manifold of an unstable periodic orbit. Thus we approximate the flow in the neighborhood of the orbit by its linearization. The key balance in determining the existence or otherwise of the hotspot is that between the compressive flow in the direction of the stable manifold (represented by the term −λξA ξ in (4)) and the tendency of the hotspot to spread as a reactiondiffusion front (represented by the right-hand side of (4)), and so we reduce consideration to one-dimensional variations along the stable manifold, and ignore variations along the unstable manifold. This remains a complicated model, and so we take the further gross liberty of assuming that the flow field along the stable manifold is timeindependent, with a uniform and constant rate of compression λ. For example, for the period-6 orbit, the value of λ is chosen such that exp(−6λT ) corresponds to the compressive eigenvalue of the Jacobian matrix of the periodic orbit (i.e., λ = (− log 0.0534)/(6T ) = 0.3052). We are thus led to (4) as a simplified one-dimensional model. Since nontrivial solutions to (4) are localized around the origin, we may make the further simplification that the corresponding spatial domain is −∞ < ξ < ∞.
In the presence of the advection term, such an exact characterization of solutions to (4) is no longer possible, although (large-Da) asymptotic and numerical treatments are available [13]. Of particular interest in the present context is the behavior of solutions to (4) on −∞ < ξ < ∞ from localized initial disturbances, which is well documented [6,11,13,22]: if Da/λ is less than some threshold f (α) then all initial states decay to A ≡ 0; however, if Da/λ > f (α) then an above-threshold initial state generally tends to a nontrivial, stable, steady localized state as t → ∞. The transition at Da/λ = f (α) is a saddle-node bifurcation [11,24]. For α = 0.25, as in our two-dimensional simulations, we find this saddlenode bifurcation to be at f (0.25) = 11.13 (by solving (4) numerically and tracking the stable steady state as Da is reduced). Since λ = 0.3052, we thus conclude that no hotspots should be possible in the two-dimensional simulations for Da < 3.4. It is certainly the case that we have observed no hotspots for such small Damköhler numbers; indeed, the lowest Damköhler number for which we have observed hotspots is 6. Given the small sample of parameter space provided by our simulations and the gross assumptions made in arguing for (4), this seems to represent reasonable agreement between this model and the full simulations. Further quantitative comparisons may be made by examining the structure of the hotspot itself. To this end, Fig. 6 shows a cross-section of the concentration through the hotspot, from two-dimensional simulations of (1) and from the model (4). The general concentration profile is captured well by the model, although the hotspot itself is rather narrower than the one-dimensional model prediction. The one-dimensional model also predicts that only sufficiently wide Gaussian initial conditions will generate a hotspot, for supercritical Damköhler numbers [13]. This dependence upon the initial condition is qualitatively consistent with the full system, as illustrated in Fig. 4.
Of course, the one-dimensional model (4) has only limited predictive power for the full two-dimensional system. While the saddle-node bifurcation of the model seems to capture the demise of the hotspot as the Damköhler number is decreased, as illustrated in Fig. 1, there is no mechanism in the model for global excitation (instead the model predicts a hotspot whose width grows like Da 1/2 for large Da [13]). In other words, the model predicts behavior qualitatively consistent with the left-hand boundary of the hotspot region in Fig. 3, but not the right-hand boundary. A further shortcoming of the model is that its dependence upon the Péclet number can be removed by a straightforward scaling of ξ, unlike the full twodimensional problem. As a consequence, the saddle-node bifurcation in (4) is independent of Pe and, for a given Damköhler number, the stable steady-state solution simply becomes thinner as the Péclet number is increased, with width proportional to Pe −1/2 . Such behavior is in contrast to the precipitate loss of the hotspot illustrated in Fig. 5 for Pe = 10 4 , and the left-hand boundary of the hotspot region in Fig. 3.

IV. CONCLUSIONS AND DISCUSSION
For a bistable chemical reaction mixed in a chaotic laminar fluid flow, we have demonstrated numerically localized states that seem to persist indefinitely. In a sweep of parameter space, these localized states provide a case intermediate between global quenching and global excitation of the reaction, as the Damköhler number is varied. They are associated with unstable an periodic orbit, and their existence may partially be understood in terms of a simple lamellar model. The long-time behavior of the system is strongly dependent on the size of the initial disturbance: too large or small and the system approaches, respectively, globally excited or quenched states; only for initial disturbances of intermediate size can a hotspot be found.
The broad influence of the Damköhler number upon the long-time state of the system is straightforward to summarize: if we begin with parameter values that lead to a hotspot, then a sufficient increase or decrease in the Damköhler number takes the system instead to the globally excited or quenched states.
The influence of the Péclet number is slightly more complicated. Again suppose that we start with parameter values leading to a hotspot, and consider the effects of reducing the Péclet number (i.e., enhancing diffusion). For lower reaction rates (i.e., smaller Damköhler numbers), we find that the enhanced initial smearing of a localized initial disturbance results ultimately in a globally quenched state. When the reaction rate is greater, however, the nascent hotspot is able to sustain itself even in the presence of the enhanced diffusive spreading, but is unable to remain localized when the Péclet number is reduced sufficiently; thus the system approaches the globally excited state.
Since we are unable accurately to simulate Péclet numbers greater than approximately 10 5 , any large-Péclet number trends are inaccessible to our numerics. Nevertheless, we conjecture that a hotspot generally persists as the Péclet number is increased (all other parameters being fixed), with the amount of chemical present scaling as Pe −1/2 . Of course, an exception to this suggestion would arise if one chose an initial condition for which diffusion was essential in transporting enough of an initial patch of chemical to the vicinity of the periodic orbit upon which the hotspot would ultimately be attached.
Then it is possible that too large a Péclet number would lead to insufficient material reaching the intended periodic orbit, and hence to a globally quenched state.
Our discovery of these "hotspots" resulted from a lucky choice of localized initial state, which happened to be centered close to an unstable period-6 orbit. However, since the initial disturbance lay in the chaotic region (apart from regular islands too small to see in a Poincaré section of the flow) the initial state must also have covered infinitely many other unstable periodic orbits, and we have no explanation for why one particular orbit should have been selected as the location of the hotspot. It seems likely that other choices of initial conditions and parameter regimes will lead to robust hotspots besides the period-6 orbit evidenced here.
Although we have demonstrated hotspots for a single specific model reaction scheme, we expect similar hotspots to be observable in other stirred multistable chemical reactions (subject to sufficient luck or judgment in choosing the initial conditions and parameter values).
Finally, we note that the existence and indefinite coherence of the hotspots reported here is due to the periodic nature of the mixing flow u. In a randomized sine flow, for example, in which the phase of the underlying shear flows is chosen randomly at each period, we would not expect a similar persistent phenomenon.