Lina Hu

Claude-Michel Brauner
email: cmbrauner@gmail.com

AND Jie Shen
email: shen@math.purdue.edu

Gregory I Sivashinsky

MODELING AND SIMULATION OF FINGERING PATTERN FORMATION IN A COMBUSTION MODEL

Keywords: Gas-solid combustion, Smouldering, Free interface, Spectral method, Asymptotics, Fingering, Pattern formation

We consider a model of gas-solid combustion with free interface proposed by L. Kagan and G.I. Sivashinsky. Our approach is twofold: (I) we eliminate the front and get to a fully nonlinear system with boundary conditions; (II) we use a fourth-order pseudo-differential equation for the front to achieve asymptotic regimes in rescaled variables. In both cases, we implement a numerical algorithm based on spectral method and represent numerically the evolution of the char. Fingering pattern formation occurs when the planar front is unstable. A series of simulations is presented to demonstrate the evolution of sparse fingers (I) and chaotic fingering (II).

1. Introduction 1.1. Physical background. Combustion is basically a process of fast oxidation accompanied by substantial heat release. Combustion phenomena are very diverse in nature and may occur in different regimes. One of them of a great practical importance, is when combustion assumes a form of a self-sustained or driven reaction wave spreading subsonically or supersonically at a well-defined speed.

Over the past forty years the theory of combustion waves, in both homogeneous and heterogeneous systems, reached rather a high level of conceptual coherence, and to date it is perhaps one of the most elegant areas of classical phenomenology, presenting a graphic example of how much natural phenomena could be deduced from a few fundamental principles. Some of the main achievements are summarized in recent surveys [START_REF] Buckmaster | Combustion theory and modeling[END_REF][START_REF] Williams | The role of theory in combustion science[END_REF].

A major leap forward in development of the theory occurred in the late sixties with renewed attention to the high activation energy limit introduced earlier by Zeldovich and Frank-Kamenetsky in their classical theory of flame propagation [START_REF] Zeldovich | A theory of thermal propagation of flame[END_REF]. It was then first realized that the flame structure consists basically of two layers: the narrow reaction zone where the role of convective transport is relatively minor, and the preheat zone dominated by convective and molecular transport but where the chemical reaction is comparatively weak and may be ignored. This crucial observation allowed, on one hand, effective evaluation of the flame velocity, and, on the other, reduction of the pertinent reaction-diffusion-advection system to a free interface system where the reaction rate term is replaced by a localized source of variable intensity. The free-interface formulation greatly facilitated theoretical exploration of the curved and unsteady flames. Moreover, having been derived for elementary one-step kinetics the resulting model appears to be relevant to a much wider range of reactions.

Since deflagrative combustion is normally strongly subsonic one may often utilize a quasiisobaric limit with the instantaneous equalizing of the pressure throughout the volume. The flow-field obviously has a profound effect both on the flame structure and its dynamics. However, for understanding many aspects of flame behavior the inverse effect, that is the influence of the flame on the background flow-field, may often be ignored and the flow-field may be regarded as prescribed. Formally this may be achieved by ignoring the burning gas thermal expansion, i.e. assuming the gas density to be independent of temperature. Although in real life systems the thermal expansion is never small, the constant-density models (or models where the density variation is accounted for but treated as a small perturbation [START_REF] Sivashinsky | Instabilities, pattern formation and turbulence in flames[END_REF] ) proved to be highly successful in interpreting even rather subtle and complex effects and far beyond the models' nominal range of validity. Note that at small density variations one may sometimes ignore the comparatively weaker effect of the flame-generated vorticity, and enjoy all the technical advantages this entails. In many physically relevant situations the structure of a curved flame may be regarded as quasiplanar and quasi-steady. In these circumstances the problem sometimes allows for separation of the longitudinal and transversal variables and thereby reduction of the effective dimensionality of the system. As a result, instead of dealing with an often difficult free-interface problem one ends up with an explicit equation for the flame interface advantageous both for physical analysis and numerical simulations [START_REF] Sivashinsky | Some developments in premixed combustion modeling[END_REF].

The aforementioned are the so-called rational approximations extractable from the firstprinciple conservation laws by pushing certain parameters to their limits. Yet, simplifications obtained by this strategy are sometimes still far below the tractability threshold. To get to the core of the phenomenon one therefore has no choice but to turn to more risky ad hoc approximations, which cannot be obtained as distinguished limits, but which are nevertheless believed to have a good enough contact with the first-principle formulation, and thus capable of capturing the essence of the phenomenon being studied. Here one may mention the so-called sandwich type models often allowing description of the essentially multi-dimensional processes by means of the appropriately coupled one-dimensional models [START_REF] Brailovsky | On stationary and travelling flame balls[END_REF][START_REF] Buckmaster | Edge-flames and their stability[END_REF][START_REF] Kagan | Pattern formation in flame spread over thin solid fuels[END_REF][START_REF] Weber | Modeling in Combustion Science[END_REF]. An impressive fact about combustion dynamics is that many of its fundamental aspects are quite robust against chemical and fluid-dynamical uncertainties and complexities and in this sense are perfectly suited for mathematical modeling. This is one of the miracles of combustion dynamics and certainly a blessing to the theorists whose analytical means are, as always, rather scanty. Apart from their technological relevance, combustion waves constitute a truly fascinating dynamical system, displaying an amazingly rich variety of phenomena such as non-uniqueness of possible propagation regimes, their birth (ignition) and destruction (extinction), chaotic self-motion and fractal-like growth, and various hysteretic transitions.

1.2. The smouldering model. In this paper, we are interested in the model of gas-solid combustion proposed in [START_REF] Kagan | Pattern formation in flame spread over thin solid fuels[END_REF]. It was motivated by experimental studies of Zik and Moses (see [START_REF] Zik | Fingering instability in solid fuel combustion: the characteristic scales of the developed state[END_REF][START_REF] Zik | Fingering instability in combustion: an extended view[END_REF][START_REF] Zik | Fingering instability in combustion[END_REF]), where the authors observed a striking fingering pattern in flame spread over thin solid fuels (although the physical process is actually flameless, we use the traditional term "flame" or "flame front" for the moving reaction zone). Among other relevant references, see [START_REF] Ijiomaa | Pattern formation in reverse smouldering combustion: a homogenisation approach[END_REF][START_REF] Ikeda | Mathematical treatment of a model for smoldering combustion[END_REF][START_REF] Lu | Pattern formation in reverse filtration combustion[END_REF]. Oxygen flows to the right and the combustion front to the left (Fig. 2 in [START_REF] Zik | Fingering instability in combustion[END_REF], courtesy of E. Moses).

The authors of [START_REF] Kagan | Pattern formation in flame spread over thin solid fuels[END_REF] interpreted the fingering phenomenon in terms of the diffusive instability similar to the one occurring in laminar flames of low-Lewis-number premixtures see, e.g. [START_REF] Buckmaster | Theory of Laminar Flames[END_REF]. They derived a free-interface system for the scaled temperature θ, the excess enthalpy S, the surface mass fraction of the solid product (char) Σ, and the moving front x = ξ(t, y):

U θ x = ∆θ, x < ξ(t, y), θ = 1, x ≥ ξ(t, y), (1.1)
θ t + U S x = ∆S -α∆θ, x ̸ = ξ(t, y), (1.2)
Σ = 0, x ≤ ξ(t, y), Σ t = 0, x > ξ(t, y), (1.3)
where α is a real number, see [19, p. 274]. Note that System (1.1)-(1.3) is rather unusual since it is the temperature's time derivative which appears in the enthalpy equation (1.2), in contrast to the classical NEF system (see [START_REF] Matkowsky | An asymptotic derivation of two models in flame theory associated with the constant density approximation[END_REF][START_REF] Sivashinsky | On flame propagation under conditions of stoichiometry[END_REF]). At the corrugated front x = ξ(t, y), [θ ] = [S ] = 0, whereas the following jump conditions occur:

[ θ n ] = -exp(S), [ S n ] = α [ θ n ] , V n [ Σ ] = exp(S), (1.4)
where

V n = - ξ t √ 1 + (ξ y ) 2
is the normal velocity.

Eventually, as x → ±∞,

θ(t, -∞, y) = S(t, -∞, y) = S x (t, +∞, y) = 0. (1.5)
Here U is the prescribed flow intensity with 0 < U < 1. As far as the variable y is concerned, it belongs to some fixed interval (-ℓ/2, ℓ/2) with periodic boundary conditions. Hereafter in Part II, we will take ℓ asymptotically large as in [START_REF] Brauner | Asymptotic analysis in a gas-solid combustion model with pattern formation[END_REF]. Thus, the parameters of the problem are α and ℓ. System (1.1)-(1.5) admits a planar solution traveling at constant velocity -V , where V = -U ln U . In the coordinate x ′ = x + V t and omitting the prime, the Traveling Wave reads:

θ = exp(U x), x ≤ 0, θ = 1, x ≥ 0, (1.6)
S = (α -ln U )U x exp(U x) + (ln U ) exp(U x), x ≤ 0, S = ln U, x ≥ 0, (1.7) Σ = 0, x < 0, Σ = -(ln U ) -1 , x > 0. (1.8)
It transpires that the fingering pattern formation is related to the instability of the planar front, whenever the parameter α is larger than the critical value

α c (ℓ, U ) = 1 + ln U + 16π 2 ℓ 2 U 2 .
(1.9)

As usual in Free Boundary Problems one fixes the moving front: we set

ξ(t, y) = -V t + φ(t, y),
where φ is the perturbation of the planar front. In the coordinate x ′ = x -ξ(t, y) and again omitting the prime, System (1.1)-(1.3) reads:

U θ x = ∆ φ θ, x < 0, θ = 1, x ≥ 0, (1.10)
θ t + (V -φ t )θ x + U S x = ∆ φ S -α∆ φ θ, x ̸ = 0, (1.11) Σ = 0, x < 0, Σ t + (V -φ t )Σ x = 0, x > 0, (1.12)
where

∆ φ = (1 + (φ y ) 2 )D xx + D yy -φ yy D x -2φ y D xy .
The front is now fixed at x = 0. The jump conditions at x = 0 read: First, in Section 2, we use the method of [START_REF] Brauner | Instabilities in a two-dimensional combustion model with free boundary[END_REF][START_REF] Lorenzi | A free boundary problem stemmed from combustion theory. I. Existence, uniqueness and regularity results[END_REF][START_REF] Lorenzi | A free boundary problem stemmed from combustion theory. II. Stability, instability and bifurcation results[END_REF] to eliminate the front in System (1.10)-(1.13) and get to a fully nonlinear problem (see [START_REF] Lunardi | Analytic Semigroups and Optimal Regularity in Parabolic Problems[END_REF]): the price to pay is the nonlocal and second-order nonlinear terms f i , see below (2.7). The free-interface conditions are replaced by fixed boundary conditions at x = 0. Although the new system looks more involved because of the nonlocal second order terms, this method has appeared to be quite efficient also from a numerical viewpoint, see e.g. [START_REF] Baconneau | Computation of bifurcated branches in afree boundary problem arising in combustion theory[END_REF] for the computation of two-dimensional bifurcated branches in a DDT model.

√ 1 + (φ y ) 2 [ θ x ] = -exp(S), [ S x ] = α [ θ x ] , (V -φ t ) [ Σ ] = √ 1 + (φ y )
The linearized problem around the null solution is studied in Section 3 where we perform a linear stability analysis (we refer to the forthcoming article [6] for the nonlinear stability). Theorem 3.1 states that the threshold of linear stability occurs at α c , see (1.9). The main feature of the proof is writing the linearized problem as the coupling of two subsystems, respectively for v and (w, h). The elliptic equation for v presents abstract issues that we examine in the spirit of [START_REF] Brauner | Bifurcation of nonplanar travelling waves in a free boundary problem[END_REF]. We investigate the spectral properties of the linear operator in weighted spaces.

Section 4 is devoted to the numerical method for solving the linearized system. We implement a backward-Euler scheme in time and discrete Fourier transform in the y direction to reduce the system to a sequence of one-dimensional problems, then solve the one-dimensional system by a Chebyshev collocation method. Beside the validation of the algorithm, an important issue is the determination of a numerical approximation of the critical α c with high accuracy. The (perturbation of the) char, a by-product of v, w and h, is computed via a transport equation. The method extends nicely to the nonlinear system in Section 5. Section 6 is devoted to the numerical results. We observe for α > α c the formation of a gamut of sparse fingers (as in Fig. 1.1(e)) which are pictured for ℓ = 2π and 0 < t < 100.

Moreover, we have been interested in identifying more intricate fingering regimes, including chaotic fingers, as already observed asymptotically in [START_REF] Kagan | Pattern formation in flame spread over thin solid fuels[END_REF], for very large ℓ and time. This is the goal of Part II of the paper. Whereas Part I's main feature was the elimination of the (perturbation of the front) φ, the second part is based on the opposite idea of deriving a selfconsistent equation for φ which, in the spirit of the derivation of the Kuramoto-Sivashinsky equation in [START_REF] Sivashinsky | On flame propagation under conditions of stoichiometry[END_REF] (see also [START_REF] Brauner | On the κ-θ model of cellular flames: existence in the large and asymptotics[END_REF][START_REF] Brauner | Weakly nonlinear asymptotics of the κ-θ model of cellular flames: the Q-S equation[END_REF][START_REF] Brauner | A fully nonlinear equation for the flame front in a quasi-steady combustion model[END_REF][START_REF] Brauner | On a strongly damped wave equation for the flame front[END_REF]), may capture most of the dynamics and yields a reduction of the effective dimensionality of the system.

In Section 7, we revisit the fourth-order pseudo-differential equation derived in [START_REF] Brauner | Asymptotic analysis in a gas-solid combustion model with pattern formation[END_REF] within the framework of a simplified, quasi-steady version of (1.10)-(1.13):

d dt U √ U 2 I -4D yy (φ) + 4φ yyyy + (γ -1)U 2 φ yy = F ((φ y ) 2 ), (1.14)
with the notation γ = α -ln U , which eventually enables a simple computation of the char. Clearly (1.14) is a generalization of the Kuramoto-Sivashinsky equation (see, e.g., [START_REF] Temam | Infinite-Dimensional Dynamical Systems in Mechanics and Physics[END_REF]), however it is fully nonlinear since the nonlinearity is also of the fourth-order. An asymptotic model is formulated by rescaling the dependent and independent variables, in line with [START_REF] Sivashinsky | On flame propagation under conditions of stoichiometry[END_REF]. As in [START_REF] Brauner | Asymptotic analysis in a gas-solid combustion model with pattern formation[END_REF], we use a Fourier spectral method to solve numerically the latter problem, see Section 8. Finally, a series of fingering pattern formation is presented in Section 9. A richer dynamics appears for large time t and spatial coordinate y, up to a chaotic fingering. Our numerical computations confirm the observations of [START_REF] Kagan | Pattern formation in flame spread over thin solid fuels[END_REF] and the relevance of the model.

PART I 2. Elimination of the moving front

To decouple the system we argue as in [START_REF] Brauner | Instabilities in a two-dimensional combustion model with free boundary[END_REF], introducing the new unknowns v, w and σ defined by: θ(t, x, y) = θ(x) + φ(t, y)θ x (x) + v(t, x, y), x < 0, -ℓ/2 < y < ℓ/2, (2.1a)

S(t, x, y) = S(x) + φ(t, y)S x (x) + w(t, x, y), x < 0, -ℓ/2 < y < ℓ/2, (2.1b) Σ(t, x, y) = Σ(x) + σ(t, x, y), x > 0, -ℓ/2 < y < ℓ/2. (2.1c)
Thanks to the jump conditions

[ θ ] = [ θ ] = 0, [ θ x ]
= -U and since θ(0) = 1, the main feature is that we can eliminate the front φ(t, y), getting:

φ(t, y) = -U -1 v(t, 0, y). (2.2)
It is convenient to work with the system of independent variables (v, w, h, σ), where

h(t, x, y) = w(t, -x, y), σ(t, x, y) = σ(t, -x, y), x < 0, -ℓ/2 < y < ℓ/2.
Then, the interface conditions become new boundary conditions at x = 0, which read:

    
αv(t, 0, y) + h(t, 0, y) -w(t, 0, y) = 0, (2.3a) αU v(t, 0, y) -U (ln U )v(t, 0, y) -h x (t, 0, y) -w x (t, 0, y) + αv x (t, 0, y) = 0, (2.3b) U v(t, 0, y) -v x (t, 0, y) = U -(1 + (U -1 v y (t, 0, y)) 2 ) -1 2 exp(S(t, 0, y)).

(2.3c)

We may write (2.3c) as:

v x (t, 0, y) -U v(t, 0, y) = U h(t, 0, y) + g,
where g = g(v y (t, 0, y), h(t, 0, y)) contains the higher order terms:

g = (1 + (U -1 v y (t, 0, y)) 2 ) -1 2 exp(S(t, 0, y)) -U -U h(t, 0, y).
Finally, the condition at x = 0 for the (perturbation of the) char reads:

σ(t, 0, y) = (1 + (φ y ) 2 ) 1 2 U e h(t,0,y) -U -φt ln U V -φ t , (2.4)
where

φ y (t, y) = -U -1 v y (t, 0, y) and φ t (t, y) = -U -1 v t (t, 0, y).
Eventually, we get to a system for the unknowns v, w, h and σ only. We denote:

f 1 (v) =(U -1 v y (0)) 2 (θ x -U -1 v(0)θ xx + v x ) x + 2U -1 v y (0)(-U -1 v y (0)θ x + v y ) x + U -1 v yy (0)(-U -1 v(0)θ x + v) x , (2.5)
f 2 (v, w) =(U -1 v y (0)) 2 ( Sx -U -1 v(0) Sxx + w x ) x + 2U -1 v y (0)(-U -1 v y (0) Sx + w y ) x + U -1 v yy (0)(-U -1 v(0) Sx + w) x -αf 1 (v), (2.6)
f 3 (v, h) =(U -1 v y (0)) 2 h xx -2U -1 v y (0)h xy -U -1 v yy (0)h x .
(2.7)

The new system reads, for x < 0, -ℓ/2 < y < ℓ/2:

             ∆v -U v x = -f 1 (v), ∆w -U w x = α∆v + V v x + v t -f 2 (v, w) + U -1 v t (t, 0, y)(-U -1 v(t, 0, y) θx + v) x , ∆h + U h x = -f 3 (v, h), (2.8)
       v x (0) -U v(0) = U h(0) + g(v y (0), h(0)), h(0) -w(0) = -αv(0), h x (0) + w x (0) = αv x (0) + αU v(0) -U (ln U )v(0), (2.9)
together with the transport equation for σ:

       σ t = (V + U -1 v t (0)) σ x , σ(0) = 1 V ( 1 + v t (0) V U ) -1 { ( 1 + 1 U 2 (v y (0)) 2 ) 1 2 U e h(0) -U + v t (0) U ln U } .
(2.10)

All the independent variables are periodic in y variable, and, as x → -∞:

v(t, -∞, y) = w(t, -∞, y) = h x (t, -∞, y) = 0.
(2.11)

Linear stability analysis

In this section, we may omit the (perturbation of the) char which is a by-product of (v, w, h). Clearly, the null solution is an equilibrium of System (2.8)-(2.9). We consider the linearized system around the null solution which reads, whenever x < 0, -ℓ/2 < y < ℓ/2:

     ∆v -U v x =0, ∆w -U w x =α∆v + V v x + v t , ∆h + U h x =0, (3.1)
together with the linearized boundary conditions at x = 0:

     v x (0) -U v(0) =U h(0), h(0) -w(0) = -αv(0), h x (0) + w x (0) =αv x (0) + αU v(0) -U ln U v(0). (3.2)
We will prove the following result of linear stability:

Theorem 3.1. Let: α c (ℓ, U ) = 1 + ln U + 16π 2 ℓ 2 U 2 . (3.
3)

It holds: (i) if α < α c ,

then the null solution is linearly asymptotically stable;

(ii) if α > α c , then it is linearly unstable.

We rewrite (3.1)(3.2) as the coupled system of the elliptic problem for v:

{ ∆v -U v x = 0, v x (0) -U v(0) = U h(0), (3.4)
and of the subsystem for w and h:

           ∆w -U w x = α∆v + V v x + v t , ∆h + U h x = 0, h(0) -w(0) = -αv(0), h x (0) + w x (0) = αv x (0) + αU v(0) -U (ln U )v(0). (3.5)
The functions v, w and h and their derivatives, whenever they are defined and continuous, verify the periodic boundary condition u(•, -ℓ/2) = u(•, ℓ/2).

3.1.

Abstracts results for the elliptic equation. In a first step, we will focus on the subsystem (3.4) and omit the time, which appears only as a parameter in this equation. In the domain Ω = {(x, y) ∈ R 2 : x < 0, y ∈ (-ℓ/2, ℓ/2)}, we define the linear operators:

L u(x, y) = ∆u(x, y) -U u x (x, y), x < 0, -ℓ/2 < y < ℓ/2, Bu(y) = u x (0, y) -U u(0, y), -ℓ/2 < y < ℓ/2.
Since we are dealing with a stability problem, we define a weighed space, with weight q(x) = exp(-U x/2):

for 0 < β < 1, let C β q (Ω) = {u : (x, y) → q(x)u(x, y) ∈ L ∞ (Ω) ∩ C β (Ω)}, and X β = { u ∈ C β q (Ω), u(•, -ℓ/2) = u(•, ℓ/2) } . (3.6)
In the space X β , the realization L of the operator L with homogeneous boundary conditions has nice spectral properties. The results of [START_REF] Brauner | Bifurcation of nonplanar travelling waves in a free boundary problem[END_REF][START_REF] Brauner | Multi-dimensional stability analysis of planar travelling waves[END_REF] hold with minor changes, see [6] for a thorough investigation:

Lemma 3.2. For 0 < β < 1, let:

X 2+β = { u ∈ C 2+β q (Ω), D (j) y u(•, -ℓ/2) = D (j) y u(•, ℓ/2), j ≤ 2 } , D(L) = { u ∈ X 2+β , Bu = 0 at x = 0 } , L : D(L) → X β , Lu = L u = ∆u -U u x . Then: (i) L is a sectorial operator; (ii) σ(L) = (-∞, -U 2 /4] ∪ {0}
. Moreover, 0 is a single eigenvalue of L, the kernel of L being spanned by the function θ = exp(U x);

(iii) the spectral projection P associated to the eigenvalue 0 is given by:

P u = U ℓ ∫ Ω u(x, y)dxdy θ, u ∈ X β . (3.7)
Since System (3.4) is non-homogenous, we need a lifting operator which can be constructed as in [START_REF] Brauner | Multi-dimensional stability analysis of planar travelling waves[END_REF]Lemma 4.1], see also [23, p. 10]: there exists a linear continuous operator N such that

N : {g ∈ C 1+β ([-ℓ/2, ℓ/2]), g(-ℓ/2) = g(ℓ/2)} → X 2+β , (BN g)(y) = g(y), y ∈ [-ℓ/2, ℓ/2].
Let us consider a slightly more general version of (3.4):

L u = f, Bu = g, (3.8)
where f ∈ X β and g ∈ C 1+β ([-ℓ/2, ℓ/2]), g(-ℓ/2) = g(ℓ/2). We look for a solution of (3.8) in the form u = z + N g. It is easy to see that

z ∈ D(L), Lz = f -L N g, (3.9)
and, since 0 is a simple eigenvalue of operator L, (3.9) is uniquely solvable if and only if P (f -L N g) = 0. An elementary computation based on the explicit expression of the lifting operator yields that the solvability condition simply reads:

∫ Ω f (x, y) dxdy = ∫ ℓ/2 -ℓ/2
g(y) dy := Π(g).

(3.10)

Returning to (3.4), we have proved that the solvability condition for (3.4) is:

Π(h(0)) = 0. (3.11)
Based on this observation, it is easy to see that Π(v) verifies for x < 0:

L 0 Π(v) = 0, B 0 Π(v) = 0, (3.12)
where

L 0 u(x) = u xx (x) -U u x (x), x < 0, and B 0 u = u x (0) -U u(0).
In the space X 0 of functions u such that x → u(x)q(x) is uniformly continuous and bounded on R -, the realization L 0 of L 0 with domain D(L 0 ) = {u ∈ X 0 , u x , u xx ∈ X 0 } is a sectorial operator, see [START_REF] Brauner | Nonlinear stability analysis of singular travelling waves in combustion: a one-phase problem[END_REF]. Moreover, 0 is a simple eigenvalue of L 0 , the kernel of L 0 is generated by θ and the associated spectral projection is

P 0 f = U ∫ 0 -∞ f (x) dx.
Therefore, we infer from (3.12) that there exists a constant c such that Π(v) = c θ, which is likely a function of the time t. However the model does not provide any information about this constant. We argue as in [10, p. 258], observing that the fact that 0 is an eigenvalue is related to the translation invariance of the original problem. The translation invariance provides us with an additional degree of freedom, therefore we look for v such that Π(v) ∈ (I -P 0 )X 0 , in other words we fix c = 0 and:

Π(v(x, •)) = 0, x < 0. (3.13)
Finally, we observe that (3.13) yields P (v) = 0, hence v ∈ (I -P )X β . We will take advantage of (3.13) in the numerical scheme, see Section 4.1.

3.2.

Proof of Theorem 3.1. We label the eigenvalues of the realization of the operator D yy in L 2 (-ℓ/2, ℓ/2) with periodic boundary conditions: 0

= -λ 0 (ℓ) > -λ 1 (ℓ) = -λ 2 (ℓ) > . . . , with λ 1 (ℓ) = 4π 2 /ℓ 2 , etc.
First, we apply the discrete Fourier transform to (3.4) and, since we look for v such that Π(v) = 0, hereafter k ≥ 1. It comes:

       v xx (t, x, k) -λ k v(t, x, k) -U v x (t, x, k) = 0, x < 0, v x (•, 0, k) -U v(•, 0, k) = U h(•, 0, k), v(•, -∞, k) = 0, (3.14)
whose solution is clearly:

v(t, x, k) = U ν k -U h(t, 0, k)e ν k x , (3.15)
for any k = 1, 2, . . ., with the notation

ν k = U/2 + (1/2) √ U 2 + 4λ k .
Second, we solve System (3.5) via Fourier, assuming that (w, h) is a sufficiently smooth solution. We end up with the infinitely many equations for k ≥ 1:

                                         w xx (t, x, k) -λ k w(t, x, k) -U w x (t, x, k) = α U ν k -U (ν 2 k -λ k ) h(t, 0, k)e ν k x +V U ν k ν k -U h(t, 0, k)e ν k x + U ν k -U h t (t, 0, k)e ν k x , (3.16a
)

h xx (t, x, k) -λ k h(t, x, k) + U h x (t, x, k) = 0, (3.16b) h(t, 0, k) -w(t, 0, k) = -α U ν k -U h(t, 0, k), (3.16c
)

h x (t, 0, k) + w x (t, 0, k) = α U ν k ν k -U h(t, 0, k) + αU U ν k -U h(t, 0, k) -U ln U U ν k -U h(t, 0, k). (3.16d) From (3.16b), it is easy to show that h(t, x, k) = c(t, k)e (ν k -U )x , that we plug in (3.16a), whose solution is w(t, x, k) = b(t, k)e ν k x + Axe ν k x , where A = 1 2ν k -U ( α U ν k -U (ν 2 k -λ k ) h(t, 0, k) + U ν k -U h t (t, 0, k) + V U ν k ν k -U h t (t, 0, k)
) .

Finally we use the boundary conditions (3.16c)-(3.16d) for w(t, x, k) and h(t, x, k), to infer the ODE for h(t, 0, k):

h t (t, 0, k) = ω k (α) h(t, 0, k), (3.17)
with

U ω k (α) = -(2ν k -U ) 2 (ν k -U ) -αU (ν 2 k -λ k ) + α(2ν k -U )U 2 -U 2 (ln U )(ν k -U ). Lemma 3.3. ω k (α) = 0 at the critical value α c,k = 1 + ln U + 4λ k /U 2 , k ≥ 1.
Proof. ω k (α) = 0 is clearly equivalent to:

α(2ν k U 2 -U 3 -ν 2 k U + λ k U ) = 4ν 3 k -8U ν 2 k + 5U 2 ν k -U 3 + U 2 (ln U )ν k -U 3 ln U. Since 0 < U < 1, 2ν k U 2 -U 3 -ν 2 k U + λ k U = 1 2 (U 2 √ U 2 + 4λ k -U 3 ) > 0, hence: α c,k = 4λ k √ U 2 + 4λ k + U 2 √ U 2 + 4λ k + U 2 (ln U ) √ U 2 + 4λ k -4U λ k -U 3 (1 + ln U ) U 2 √ U 2 + 4λ k -U 3 =1 + ln U + 4λ k U 2 .
□ To complete the proof of Theorem 3.1, we observe that the sequence α c,k is non-decreasing and converges to zero as k → +∞. Therefore, ω k (α) < 0 if and only if α < α c,1 := α c , and

ω k (α) > 0 if and only if α > α c,1 := α c . Since λ 1 = 4π 2 /ℓ 2 , α c = 1 + ln U + 16π 2 /ℓ 2 U 2 .
4. Numerical method for solving the linearized system 4.1. Description of the method. In view of (2.11), we may restrict the domain Ω to Ω L := {(x, y) ∈ R 2 : x ∈ (-L, 0), y ∈ (-ℓ/2, ℓ/2)}, with L sufficiently large. Next, we transform the domain Ω L to Ω L := {(x ′ , y ′ ) ∈ R 2 : x ′ ∈ (-1, 1), y ′ ∈ (0, 2π)}, by setting x ′ = 2x/L + 1, y ′ = 2πy/ℓ + π, and v(t, x, y) = v(t, (x ′ -1)L/2, (y ′ -π)ℓ/2π) := ṽ(t, x ′ , y ′ ). Omitting the primes and the tilde (the same applies to w and h), System (3.1)-(3.2) becomes:

             ( 2 L ) 2 v xx + ( 2π ℓ ) 2 v yy -U 2 L v x =0, ( 2
L ) 2 w xx + ( 2π ℓ ) 2 w yy -U 2 L w x =α( 2
L ) 2 v xx + α( 2π ℓ ) 2 v yy + V 2 L v x + v t , ( 2
L ) 2 h xx + ( 2π ℓ ) 2 h yy + U 2 L h x =0, (4.1)
with the set of boundary conditions:

           2 L v x (t, 1, y) -U v(t, 1, y) =U h(t, 1, y), h(t, 1, y) -w(t, 1, y) = -αv(t, 1, y), 2 L h x (t, 1, y) + 2 L w x (t, 1, y) = 2 L αv x (t, 1, y) + αU v(t, 1, y) -U ln U v(t, 1, y), (4.2)
and

v(t, -1, y) = w(t, -1, y) = 2 L h x (t, -1, y) = 0. (4.3)
First, we discretize the above system with a backward-Euler scheme in time:

                     ( 2 L ) 2 v n+1 xx + ( 2π ℓ ) 2 v n+1 yy -U 2 L v n+1 x =0, ( 2
L ) 2 w n+1 xx + ( 2π ℓ ) 2 w n+1 yy -U 2 L w n+1 x =α( 2 L ) 2 v n+1 xx + α( 2π ℓ ) 2 v n+1 yy + V 2 L v n+1 x + v n+1 -v n ∆t , ( 2
L ) 2 h n+1 xx + ( 2π ℓ ) 2 h n+1 yy + U 2 L h n+1 x =0, (4.4)
           2 L v n+1 x (1, y) -U v n+1 (1, y) =U h n+1 (1, y), h n+1 (1, y) -w n+1 (1, y) = -αv n+1 (1, y), 2 L h n+1 x (1, y) + 2 L w n+1 x (1, y) = 2 L αv n+1 x (1, y) + αU v n+1 (1, y) -U ln U v n+1 (1, y), (4.5
)

v n+1 (-1, y) = w n+1 (-1, y) = 2 L h n+1 x (-1, y) = 0. (4.6)
Second, we use the discrete Fourier transform in the y direction to reduce the above system to a sequence of one-dimensional problems. More precisely, we plug in (4.4)-(4.6):

v n+1 = k= N 2 ∑ k=-N 2 vn+1 k (x)e iky , w n+1 = k= N 2 ∑ k=-N 2 ŵn+1 k (x)e iky , h n+1 = k= N 2 ∑ k=-N 2 ĥn+1 k (x)e iky ,
see [START_REF] Shen | Spectral Methods: Algorithms, Analysis and Applications[END_REF]Chapt. 2] for the slight difference in ordering the modes with respect to Section 3.2.

Thanks to (3.13), the zeroth order modes are null. Then, for each k ̸ = 0 we obtain:

                     ∂ xx ( vn+1 k (x) ) ( 2 L ) 2 -k 2 vn+1 k (x)( 2π ℓ ) 2 -U ∂ x ( vn+1 k (x) ) 2 L =0, ∂ xx ( ŵn+1 k (x) ) ( 2 L ) 2 -k 2 ŵn+1 k (x)( 2π ℓ ) 2 -U ∂ x ( ŵn+1 k (x) ) 2 L =∂ xx ( vn+1 k (x) ) ( 2
L ) 2 α -k 2 vn+1 k (x)( 2π ℓ ) 2 α + V ∂ x ( vn+1 k (x) ) 2 L + vn+1 k (x) -vn k (x) ∆t , ∂ xx ( ĥn+1 k (x) ) ( 2
L ) 2 -k 2 ĥn+1 k (x)( 2π ℓ ) 2 + U ∂ x ( ĥn+1 k (x) ) 2 L =0, (4.7)
with

                 2 L ∂ x ( vn+1 k (x) ) (1) -U vn+1 k (1) =U ĥn+1 k (1), ĥn+1 k (1) -ŵn+1 k (1) = -αv n+1 k (1), 2 L ∂ x ( ĥn+1 k (x) ) (1) + 2 L ∂ x ( ŵn+1 k (x) ) (1) = 2 L α∂ x ( vn+1 k (x) ) (1)
+ αU vn+1 k (1) -U ln U vn+1 k (1), (4.8)
and

vn+1 k (-1) = ŵn+1 k (-1) = 2 L ∂ x ( ĥn+1 k (x) ) (-1) = 0. (4.9)
We will solve the above one-dimensional system by mean of a Chebyshev collocation method (see [START_REF] Shen | Spectral and High-Order Method with Application[END_REF][START_REF] Shen | Spectral Methods: Algorithms, Analysis and Applications[END_REF]). Let {l j (x)} N j=0 be the Lagrange polynomials based on the Chebyshev-Gauss-Lobatto points {x j } N j=0 = {cos( πj N )} N j=0 . We denote the differentiation matrix of order m associated to the {x j } N j=0 by D m = (d

(m) ij ) i,j=0,••• ,N , where d (m) ij = l
(m) j (x i ) and l j (x i ) = δ ij , here δ ij is the Kronecker Delta symbol. Setting:

vn+1 k (x) = N ∑ j=0 ( vn+1 k ) j l j (x), ŵn+1 k (x) = N ∑ j=0 ( ŵn+1 k ) j l j (x), ĥn+1 k (x) = N ∑ j=0 ( ĥn+1 k ) j l j (x),
and inserting the above quantities in (4.7)-(4.9), we find

                                     N ∑ j=0 ( ( 2
L ) 2 d (2) ij -k 2 ( 2π ℓ ) 2 δ ij -U 2 L d (1) ij ) ( vn+1
k ) j =0, ∆t N ∑ j=0 ( ( 2
L ) 2 d (2) ij -k 2 ( 2π ℓ ) 2 δ ij -U 2 L d (1) ij ) ( ŵn+1 k ) j =∆t N ∑ j=0 ( ( 2
L ) 2 αd (2) ij -k 2 ( 2π ℓ ) 2 αδ ij + V 2 L d (1) ij + δ ij ∆t ) ( vn+1 k ) j - N ∑ j=0 δ ij ( vn k ) j , N ∑ j=0 ( ( 2
L ) 2 d (2) ij -k 2 ( 2π ℓ ) 2 δ ij + U 2 L d (1) ij ) ( ĥn+1
k ) j =0,
with the boundary conditions

                         N ∑ j=0 2 L d (1) 0j ( vn+1 k ) j -U ( vn+1 k ) 0 =U ( ĥn+1 k ) 0 , ( ĥn+1
k ) 0 - ( ŵn+1
k ) 0 = -α ( vn+1 k ) 0 , 2 L N ∑ j=0 d (1) 0j ( ĥn+1 k ) j + 2 L N ∑ j=0 d (1) 0j ( ŵn+1 k ) j = 2 L α N ∑ j=0 d (1) 0j ( vn+1 k ) j + αU ( vn+1 k ) 0 -U ln U ( vn+1 k ) 0 ,
and

( vn+1 k ) N = ( ŵn+1 k ) N = 2 L N ∑ j=0 d (1) N j ( ĥn+1 k ) j = 0.
We may write the above system in the following compact form:

A k X k = B k where B k 3(N +1)×1 = ( 0, 0, 0, 0, -(v n k ) 1 , 0, 0, -(v n k ) 2 , 0, • • • , 0, -(v n k ) N -1 , 0, 0, 0, 0 ) T X k 3(N +1)×1 = ( (v n+1 k ) 0 , ( ŵn+1 k ) 0 , ( ĥn+1 k ) 0 , • • • , (v n+1 k ) N , ( ŵn+1 k ) N , ( ĥn+1 k ) N ) T
and A k is 3(N + 1) × 3(N + 1) matrix. Note in particular that A k is the same at each time step, so it can be pre-factorized and stored in LU form.

( ( U L 2 - 1 4 )x 2 + x 2 + 3 4 - U L 2
) ; w = cos(t) cos(y) ( x 2 4 (αU L -LU ln U -α -1) + x 2 (1 + α) + 1 4 (3 + 3α + LU ln U -αLU ) ) ; h = cos(t) cos(y) ( αU L 4 (x 2 + 2x -3) + 1
) .

Since the exact solution are low-order polynomial in x, the error is dominated by time discretization as soon as N ≥ 8 (see Fig. 4.1 (a)).

4.3.

Numerical determination of the critical α c (linear case). While the critical α c is determined explicitly in Theorem 3.1, it is useful, particularly when α c is not known explicitly, to show that our numerical scheme is capable of approximating α c to high accuracy. To this end, we use a bisection method to find a numerical approximation α cN of the critical value α c , based on the following remark: if the solution goes to zero in time or the energy decreases, then α c is on the right side of the parameter interval; if the solution blows up in time or the energy increases, then α c is on the left side of the interval. α cN error = |α c -α cN | 0.8 7.0268564487e+000 7.0268563829e+000 6.5800000000e-008 0.5 1.6306852819e+001 1.6306852639e+001 1.8000000000e-007 0.2 9.9390562088e+001 9.9390561745e+001 3.4300000000e-007 e -1 2.9556224396e+001 2.9556224135e+001 2.6100000000e-007

Table 4.2. L = 10, U = 1, N = 32, dt = 0.1 ℓ α c α cN error = |α c -α cN | 2π
5.0000000000e+000 4.9999999706e+000 2.9400000000e-008 3π 2.7777777778e+000 2.7777711575e+000 6.6203000000e-006 4π 2.0000000000e+000 1.9999195019e+000 8.0498100000e-005 5π 1.6400000000e+000 1.6396664907e+000 3.3350930000e-004

We infer from the above two tables that our numerical method is able to approximate α c to a high accuracy. Moreover, Tables 4.1-4.2 show that α cN depends upon the quantities U and L in the same way as α c , see (3.3).

Next, we present some solution profiles in the three cases α < α c , α = α c and α > α c , to exhibit the dynamics of v(t, 0, y) and the energy of v(t, x, y).

Solving the nonlinear system

We are concerned with the nonlinear systems for (v, w, h, σ

) in Ω L = {-L < x < 0, -ℓ/2 < y < ℓ/2}:              ∆v -U v x = -f 1 (v), ∆w -U w x = α∆v + V v x + v t -f 2 (v, w) + U -1 v t (t, 0, y)(-U -1 v(t, 0, y) θx + v) x , ∆h + U h x = -f 3 (v, h), (5.1)
together with the nonlinear boundary conditions at x = 0:

       v x (0) -U v(0) = U h(0) + g(v y (0), h(0)), h(0) -w(0) = -αv(0), h x (0) + w x (0) = αv x (0) + αU v(0) -U (ln U )v(0), (5.2)
and the transport equation:

       σ t = (V + U -1 v t (0)) σ x , σ(0) = 1 V ( 1 + v t (0) V U ) -1 { ( 1 + 1 U 2 (v y (0)) 2 ) 1 2 U e h(0) -U + v t (0) U ln U } . (5.3)
By treating all the nonlinear terms explicitly, the above nonlinear system can be solved by using essentially the same algorithm described in the previous section for the linearized system. We look for v such that its zeroth mode is null according to (3.13), whereas the zeroth mode of w and h are free in the nonlinear problem (see [6]).

Set L = 10, U = 1, ℓ = 2π, α = α c . We now validate the algorithm with an exact solution satisfying the boundary conditions (5.2) with g = 0, and in Ω L : v = cos(t) cos(y)

( ( U L 2 - 1 4 )x 2 + x 2 + 3 4 - U L 2 ) 10 -4 ; w = cos(t) cos(y) ( x 2 4 (αU L -LU ln U -α -1) + x 2 (1 + α) + 1 4 (3 + 3α + LU ln U -αLU ) ) 10 -4 ; h = cos(t) cos(y) ( αU L 4 (x 2 + 2x -3) + 1 ) 10 -4 .
Again, since the solutions are low-order polynomials in x, the errors are dominated by the time discretization as soon as N ≥ 10, see Fig. 5.1. We observe also from Fig. 5.2 that the scheme is first-order accurate in time.

5.1. Determination of the critical α c (nonlinear case). Next, we use the bisection method to determine numerically the critical α cN in the nonlinear case. Numerical results in Tables 5.1 and 5.2 show that, with only N = 32, our algorithm is able to approximate α c with a very good accuracy.

Table 5.1. L = 10, ℓ = 2π, N = 32, dt = 0.1 U α c α cN error = |α c -α cN | 0
.8 7.0268564487e+000 7.0268583875e+000 1.9388000000e-006 0.5 1.6306852819e+001 1.6306852650e+001 1.6900000000e-007 0.2 9.9390562088e+001 9.9392551959e+001 1.9898710000e-003 e -1 2.9556224396e+001 2.9556230949e+001 6.5530000000e-006 Table 5.2. L = 10, U = 1, N = 32, dt = 0.1 ℓ α c α cN error = |α c -α cN | 2π 5.0000000000e+000 5.0000042870e+000 4.2870000000e-006 3π 2.7777777778e+000 2.7777777777e+000 1.0000000000e-010 4π 2.0000000000e+000 1.9999281988e+000 7.1801200000e-005 5π 1.6400000000e+000 1.6396754829e+000 3.2451710000e-004 Next, we plot some solution profiles in the three different cases α < α c , α = α c and α > α c , to exhibit the dynamics of v(t, 0, y) and the energy of v(t, x, y). We first transform the domain Ω L to Ω L := {(x ′ , y ′ ), x ′ ∈ (-1, 1), y ′ ∈ (0, 2π)}. Omitting the primes, and recalling φ(t, y) = -1/U v(t, 1, y) in the rescaled coordinates, the transformed function σ satisfies:

                 σt (t, x, y) -(V + U -1 v t (1))σ x (t, x, y) 2 L = 0, x ∈ [-1, 1), σ(t, 1, y) = 1 V ( 1 + vt(t,1,y) U V ) -1 { ( 1 + 4π 2 ℓ 2 U 2 (v y (t, 1, y)) 2 ) 1 2 U e h(t,1,y) -U + vt(t,1,y) U ln U } , x = 1, σ(0, x, y) = 0, t = 0, (5.4)
where v, w and h are solution of (5.1)-(5.2) in Ω L . Next, we use a similar algorithm as before to solve the system (5.4). We discretize in time as follows:

                 σn+1 (x,y)-σ n (x,y) ∆t -V σn x (x, y) 2 L -1 U 2 L v n+1 (1,y)-v n (1,y) ∆t σn x (x, y) = 0, x ∈ [-1, 1), σn+1 (1, y) = 1 V ( 1 + v n+1 (1,y)-v n (1,y) U V ∆t ) -1 { ( 1 + 4π 2 ℓ 2 U 2 (v n+1 y (1, y)) 2 ) 1 2 U e h n+1 (1,y) -U + v n+1 (1,y)-v n (1,y) ∆tU ln U } , x = 1, σ0 (x, y) = 0.
(5.5) Once σ is computed, we may calculate Σ as follows:

Σ(t, x, y) = Σ + σ(t, -x, y), Σ = -(ln U ) -1 , x > 0.
Then, we are in position to compute the temperature θ by mean of formulas (2.1a), (2.2) and (1.6).

Figures (I)

In view of a better representation of the char Σ(t, x, y), of the values of the char at the front Σ(t, 0, y) (physically it means immediately behind the front) and of the temperature θ, we return to the original "physical" coordinates (x, y) in which the front is moving. However, for the convenience of the reader, we focus on the front zone and therefore fix the Traveling Wave at the origin. In other words, we do not represent the front itself but the perturbation φ(t, y) of the planar front and its immediate neighboring.

For all the figures hereafter, initial data are:

v 0 = cos( 2π ℓ y + π) ( ( 2 L x + 1) 3 -( 2 L x + 1) 2 -( 2 L x + 1) + 1 ) 10 -4 , w 0 = cos( 2π ℓ y + π) ( -0.5 L( 2 L x + 1) 2 + 0.5L ) 10 -4 , h 0 = cos( 2π ℓ y + π) ( 0.25 L( 2 L x + 1) 2 + 0.5 L( 2 L x + 1) -0.75 L ) 10 -4 ,
where ℓ = 2π and L = 10. We have performed numerical experiments with different gamuts of initial data without any significant difference in the fingering evolution after a brief transition period.

We first present the evolution of sparse fingers over three spatial periods in Figs. 6.1-6.2, respectively for U = 0.2 in Fig. 6.1 and U = 0.9 in Fig. 6.2. Second, we illustrate in Fig. 6.3 the relation between the flame front, the temperature and the char over three spatial periods at t = 200, for U = 0.2 and U = 0.9. Clearly, the fingering pattern and the temperature profile depend on the structure of the flame front. The temperature is represented in the "fresh" side x ≤ φ(t, y), whereas the char pattern lies in the "burnt" side x ≥ φ(t, y). Near the flame front, the value of the temperature approaches 1.

Next, Fig. 6.4 focuses on the fingers for different values of U , over two spatial periods. We observe that the width of fingers increases as U becomes larger. For convenience, in this figure the fingers are rotated upwards.

PART II

In contrast to the previous part where we eliminated the (perturbation of) the front, this part is based on a self-consistent equation for the (perturbation of) the front derived in [START_REF] Brauner | Asymptotic analysis in a gas-solid combustion model with pattern formation[END_REF]. In this respect, we use a simplified, quasi-steady version of System (1.10)-(1.13), motivated from the fact that not far from the instability threshold, the time derivatives exhibit a relatively small effect: the dynamics appears to be essentially driven by the front. We perform a rescaling of the dependent and independent variables, in the spirit of [START_REF] Sivashinsky | On flame propagation under conditions of stoichiometry[END_REF]. A rich dynamics appears for large time t and spatial coordinate y, up to observing a chaotic fingering evolution.

7. The asymptotic model 7.1. Derivation of the front equation. We briefly recall some basic results of [START_REF] Brauner | Asymptotic analysis in a gas-solid combustion model with pattern formation[END_REF]: we consider a simplified, quasi-steady version of the model, replacing (1.10)-(1.12) by

U θ x = ∆ φ θ, x < 0, θ = 1, x ≥ 0, (V -φ t )θ x + U S x = ∆ φ S -α∆ φ θ, x ̸ = 0, Σ = 0, x < 0, Σ x = 0, x > 0,
with the interface conditions (1.13), namely:

√ 1 + (φ y ) 2 [ θ x ] = -exp(S), [ S x ] = α [ θ x ] , (V -φ t ) [ Σ ] = √ 1 + (φ y ) 2 exp(S). (7.1)
In this quasi-steady model, the dynamics appears to be essentially driven by the front which is still an unknown of the problem. Next, we define the respective perturbations of temperature u, enthalpy v and char σ:

θ = θ + u, S = S + v, Σ = Σ + σ.
As in [START_REF] Brauner | Asymptotic analysis in a gas-solid combustion model with pattern formation[END_REF], in the equations satisfied by u, v and φ, we keep only linear and second-order terms for φ, and first-order terms for u and v.

The final system for u and v reads:

                 U u x -∆u = (∆ φ -∆)θ, x < 0, V u x -∆(v -αu) + U v x -φ t θ x = (∆ φ -∆)(S -αθ), x ̸ = 0, u(0 -) = [v] = 0, U v(0) -u x (0 -) = 1 2 (φ y ) 2 U, [v x ] = -αu x (0 -). (7.2)
Let us examine the (perturbation of the) char more carefully: σ(t, x, y) = σ(t, 0, y), x > 0, -ℓ/2 < y < ℓ/2, (7.3)

(V -φ t )σ(t, 0, y) = (1 + (φ y ) 2 ) 1
2 U e v(t,0,y) -U + φ t Σ. (7.4) Keeping again only linear and second-order terms for φ, and first-order terms for v, it is not difficult to find the relation:

V σ(0) = U v(0) - φ t ln U + U 2 (φ y ) 2 - 1 ln U v(0)φ t + (φ t ) 2 U (ln U ) 2 . (7.5)
To compute σ(0), hence σ, thanks to formula (7.5), the next step is to link v(0) with φ.

Using the convenient notation γ = α -ln U , it follows from [5, Section 3.2] that φ verifies the pseudo-differential fourth-order equation:

d dt U √ U 2 I -4D yy (φ) = -4φ yyyy -(γ -1)U 2 φ yy + F ((φ y ) 2 ), (7.6)
with

F = 1 4 (U 2 I -4D yy ) 3 2 - 3 4 U (U 2 I -4D yy ) -γU 2 ( √ U 2 I -4D yy -U ) ,
and periodic boundary conditions at y = ±ℓ/2. Actually, (7.6) is a fully nonlinear equation since the nonlinear part F ((φ y ) 2 ) is of the same order as the linear operator. We refer to [START_REF] Brauner | Asymptotic analysis in a gas-solid combustion model with pattern formation[END_REF]Cor. 3.1] for the nonlinear stability analysis of the null solution of (7.6). The threshold of stability occurs at γ c (ℓ) = 1 + 16π 2 /ℓ 2 U 2 , namely at α c , see (3.3).

Let us use the discrete Fourier transform as in Section 3.2. We set

X k = √ U 2 + 4λ k . According to [5, Section 3.1]: v(t, 0, k) = 2U (X k + U )X k φ t (t, k) + (U -X k )U 2 γ 2X 2 k φ(t, k) + 4γU 2 X k -4γU 3 2X 2 k (X k + U ) (φ y ) 2 (t, k), which is equivalent to 1 2 X 2 k (X k + U ) v(t, 0, k) =(X k U ) φ t (t, k) -U 2 γλ k φ(t, k) + 1 4 (4γU 2 X k -4γU 3 ) (φ y ) 2 (t, k), (7.7)
and in physical variable:

1 2 (U 2 I -4D yy ) ( √ U 2 I -4D yy + U ) v(0) = U √ U 2 I -4D yy φ t + γU 2 φ yy + γU 2 ( √ U 2 I -4D yy -U ) (φ 2 y ), (7.8)
which determines v(0).

7.2.

Rescaling. We define a small perturbation parameter ε > 0 by γ = 1 + ε, and rescale dependent and independent variables accordingly:

t = τ ε 2 U 2 , y = η √ εU , φ = ε U ψ.
One links ε and the period ℓ, which will blow up as ε → 0: for ℓ 0 > 4π fixed, we take ℓ = ℓ 0 / √ εU , hence γ c (ℓ) → 1 as ε → 0. After division by ε 3 and U 3 we get:

∂ ∂τ

( √ I -4εD ηη ) ψ = -4D ηηηη ψ -D ηη ψ + 1 4
{ (I -4εD ηη ) 3 2 -3(I -4εD ηη ) -4(1 + ε) ( √ I -4εD ηη -I )} (D η ψ) 2 . (7.9)
In the same way, we rescale v(t, 0, y) and σ(t, 0, y) as follows:

v(t, 0, y) = ε 2 v 1 (τ, 0, η), σ(t, 0, y) = ε 2 U σ 1 (τ, 0, η), therefore V σ 1 (τ, 0, η) = v 1 (τ, 0, η) + ε ( - ψ τ ln U + 1 2 (ψ η ) 2 ) + o(ε) (7.10)
(we will forget the lower order terms hereafter). According to (7.8):

1 2

{ (I -4εD ηη ) ( √ I -4εD ηη + I )} v 1 (0) = ε ∂ ∂τ ( √ I -4εD ηη ) ψ + γD ηη ψ + ε 1 4 { 4(1 + ε) ( √ I -4εD ηη -I )} (D η ψ) 2 . (7.11)

Numerical method for the asymptotic model

Let us briefly recall the method in [START_REF] Brauner | Asymptotic analysis in a gas-solid combustion model with pattern formation[END_REF] for solving (7.9). First, we set ζ = η/2 l0 + π (where l0 = ℓ 0 /4π) to reformulate (7.9) on the interval [0, 2π). Next, we define the bifurcation parameter β = 4 l0 (see also [START_REF] Hyman | The Kuramoto-Sivashinsky equation: a bridge between PDEs and dynamical systems[END_REF]) and rescale the time by setting t ′ = τ /β 2 . For simplicity we denote ε/β by ϵ. Then, omitting the prime, (7.9) becomes:

∂ ∂t ( √ I -4ϵD ζζ ) ψ = -4D ζζζζ ψ -βD ζζ ψ + β 4
{ (I -4ϵD ζζ ) 3 2 -3(I -4ϵD ζζ ) -4(1 + ε) ( √ I -4ϵD ζζ -I )} (D ζ ψ) 2 . ( 8.1)
Similarly, we have

1 2 β 2 { (I -4ϵD ζζ ) ( √ I -4ϵD ζζ + I )} v 1 (0) =ε ∂ ∂t ( √ I -4ϵD ζζ ) ψ + γβD ζζ ψ + ε β 4 { 4(1 + ε) ( √ I -4ϵD ζζ -I )} (D ζ ψ) 2 , (8.2)
and (7.10) yields:

V σ 1 (t, 0, ζ) = v 1 (t, 0, ζ) + ε 2β (ψ ζ ) 2 - ε β 2 ψ t ln U . (8.3)
Now we are in position to implement a numerical method to solve (8.1), (8.2) and (8.3). As in [START_REF] Brauner | Asymptotic analysis in a gas-solid combustion model with pattern formation[END_REF], we use a Fourier spectral method to solve these periodic problems. After discretizing (8.1) in time by a semi-implicit first order scheme which treats implicitly the linear terms and explicitly the nonlinear terms, then taking the discrete Fourier transform, the discrete Fourier coefficient of ψ n+1 is determined by:

ψ n+1 k = ( (1 + 4ϵk 2 ) 1 2 + 4k 4 ∆t -βk 2 ∆t ) -1 × ( (1 + 4ϵk 2 ) 1 2 ψ n k (8.4) + β∆t 4 { (1 + 4ϵk 2 ) 3 2 -3(1 + 4ϵk 2 ) -4(1 + ε)[(1 + 4ϵk 2 ) 1 2 -1] } { [(ψ x ) n ] 2 } k ) ,
together with the initial condition given by ψ(0, •) = ψ 0 . Similarly, the discrete Fourier coefficient of v 1 at x = 0 is determined by:

1 2 β 2 { (1 + 4ϵk 2 ) ( √ 1 + 4ϵk 2 + 1 )} {v 1 (0)} n+1 k = ε ( √ 1 + 4ϵk 2 ) ψ n+1 k -ψ n k ∆t -(1 + ε)βk 2 ψ n+1 k + ε β 4 { 4(1 + ε) ( √ 1 + 4ϵk 2 -1 )} {[(ψ x ) n+1 ] 2 } k , (8.5)
where {(ψ x ) 2 } k represents the k-th Fourier coefficient of (ψ x ) 2 . It yields:

{v 1 (0)} n+1 k = 2 β 2 ( (1 + 4ϵk 2 )( √ 1 + 4ϵk 2 + 1) ) -1 × ( ε( √ 1 + 4ϵk 2 ) ψ n+1 k -ψ n k ∆t -(1 + ε)βk 2 ψ n+1 k + ε β 4 { 4(1 + ε) ( √ 1 + 4ϵk 2 -1 )} {[(ψ x ) n+1 ] 2 } k ) , (8.6)
with the initial {v 1 (0)} 0 k given by

{v 1 (0)} 0 k = 2 β 2 ( (1 + 4ϵk 2 )( √ 1 + 4ϵk 2 + 1) ) -1 × ( -(1 + ε)βk 2 ψ 0 k + ε β 4 { 4(1 + ε) ( √ 1 + 4ϵk 2 -1 )} {[(ψ x ) 0 ] 2 } k ) .
Once ψ and v 1 (0) are computed and plugged into (8.3), we get:

σ n+1 1 (0) = 1 V v n+1 1 (0) + ε 2βV (ψ n+1 ζ ) 2 - ε β 2 V ln U ψ n+1 -ψ n ∆t .
Finally, we compute

Σ = Σ + σ = - 1 ln U + σ(0) = - 1 ln U + ε 2 U σ 1 (0).
As far as the temperature is concerned in this part, we recall the results in [5, Section 3.1], namely:

u(t, x, 0) = -U (φ y ) 2 (t, 0)xe U x , x ≤ 0, u(t, x, k) = U (λ k ) -1 ( U (φ y ) 2 (t, k) + λ k φ(t, k) ) ( e U x -e ν k x ) , x ≤ 0, k ≥ 1.
According to the relation θ = θ + u, we may recover θ immediately.

Figures (II)

As in Section 6, in view of a better representation of the char Σ(t, x, y), of the char at the front Σ(t, 0, y) and of the temperature θ, we return to the original "physical" coordinates (x, y) in which the front is moving. Also we only represent the perturbation φ(t, y) of the planar front and its immediate neighboring.

The initial data are as follows: The temperature is represented in the "fresh" side x ≤ φ(t, y), whereas the char pattern lies in the "burnt" side x ≥ φ(t, y).

(i) cos ( U ( ε β ) 1

Figure 1 . 1 .

 11 Figure 1.1. Instability of combustion front when V O2 decreases from (a) to (e).Oxygen flows to the right and the combustion front to the left (Fig.2in[START_REF] Zik | Fingering instability in combustion[END_REF], courtesy of E. Moses).

4. 2 .

 2 Algorithm validation. In this section, we let L = 10, U = 0.5, ℓ = 2π, α = α c , and validate the algorithm with the exact solution of (4.1)-(4.2)-(4.3): v = cos(t) cos(y)

Figure 4 . 1 .

 41 Figure 4.1. (a) Error vs. polynomial degree with dt = 0.0001, T = 1, and (b) error vs. time step with N = 16, T = 1 On the other hand, we observe clearly from Fig. 4.1 (b) that our algorithm is first order accuracy in time.

Table 4 . 1 .

 41 L = 10, ℓ = 2π, N = 32, dt = 0.1 U α c

Figure 4 . 2 .Figure 4 . 3 .

 4243 Figure 4.2. v(t, 0, y) with ℓ = 2π, U = 1 in the cases α < α c , α = α c , α > α c .

Figure 5 . 1 .

 51 Figure 5.1. Error vs. polynomial degree with dt = 0.0001, T = 1.

Figure 5 . 2 .

 52 Figure 5.2. Error vs. time step with N = 16, T = 1.

Figure 5 . 3 .Figure 5 . 4 . 5 . 2 .

 535452 Figure 5.3. v(t, 0, y) with ℓ = 2π, U = 1 in the cases α < α c , α = α c , α > α c .

  (a) t = 0 (b) t = 10 (c) t = 20 (d) t = 30 (e) t = 40 (f) t = 80

Figure 6 . 1 .

 61 Figure 6.1. Two-dimensional evolution of the char Σ at different times t=0, 10, 20, 30, 40, 80 with ℓ = 2π, U = 0.2 and α = α c + 0.5. At t = 80, the char pattern develops into sparse elongated fingers without tip-splitting. Propagation from right to left. The color grading enlightens the levels of the solid product distribution, the median level 0.6213 corresponding to the planar front.

Figure 6 . 2 .

 62 Figure 6.2. As in Fig. 6.1, two-dimensional evolution of the char Σ at different times t=0, 10, 20, 30, 40, with ℓ = 2π, U = 0.9 and α = α c + 0.5. At t = 50, the char develops into sparse elongated fingers without tip-splitting. The median level 9.4912 corresponds to the planar front.

Figure 6 . 3 .

 63 Figure 6.3. Two-dimensional profile of the flame structure at t = 200. (a) and (d): temperature, (b) and (e): flame front, (c) and (f ): char pattern, respectively for U = 0.2 (top) and 0.9 (bottom).Here ℓ = 2π, α = α c +0.5. The temperature is represented in the "fresh" side x ≤ φ(t, y), whereas the char pattern lies in the "burnt" side x ≥ φ(t, y). Near the flame front, the value of the temperature approaches 1.

Figure 6 . 4 .

 64 Figure 6.4. Two-dimensional profiles of the char pattern for different values of U within two periods at t = 200. Here ℓ = 2π, α = α c + 0.5. The darkest part of the charred area is viewed as the width of the finger: it increases as the scaled flow velocity U becomes larger.Finally, we represent the evolution of the char Σ(t, 0, y) at the flame front over two spatial periods, for U = 0.2.

Figure 6 . 5 .

 65 Figure 6.5. Evolution of the char pattern immediately behind the flame front over two periods, represented in different time scales. Here ℓ = 2π, U = 0.2, α = α c + 0.5. Propagation from bottom to top.

. 9 .

 9 4 and in Figs. 9.6-9.10. We illustrate in Figs. 9.1-9.5 the sparse finger evolution in time for different values of β: 10, 30, 62 and 100 respectively in Figs. 9.1-9.3 and 9.5, and β = 62 in Fig. 9.4. Chaotic fingering evolution in time is observed in Fig. 9.6 for β = 400. Inspection of these figures reveals that the number of fingers depends upon the values of the bifurcation parameter β: larger β's correspond to a richer dynamics. Hereafter we use the notation: T 0 = β 2 ε -2 U -2 .

  (a) t = 0 (b) t = 1.0e-02 T0 (c) t = 1.0e-01 T0 (d) t = 2.0e-01 T0 (e) t = 3.0e-01 T0 (f) t = 4.0e-01 T0

Figure 9 . 1 .

 91 Figure 9.1. Two-dimensional evolution of the char Σ at different times with initial datum (i). Here β = 10, U = 0.9, ε = 0.01. Following a transition period (a) to (e), the char pattern develops into a single finger without tip-splitting (f). Propagation from right to left. The median level 9.4912 of the color grading corresponds to the planar front.

Figure 9 . 2 .

 92 Figure 9.2. As Fig. 9.1, with larger β = 30. After a transition period (a) to (e), the char pattern develops into two elongated fingers without tip-splitting (f).

Figure 9 . 3 .

 93 Figure 9.3. As in Fig. 9.1, with larger β = 62. After a transition period (a) to (e), the char pattern develops into three elongated fingers without tip-splitting (f).

Figure 9 . 4 .

 94 Figure 9.4. As in Fig. 9.3, with initial datum (ii). After a transition period (a) to (e), the char pattern develops into three elongated fingers without tip-splitting (f).

Figure 9 . 5 .

 95 Figure 9.5. Again as in Fig. 9.1, with larger β = 100. After a transition period (a) to (e), the char pattern develops into four elongated fingers without tip-splitting (f).

Figure 9 . 6 .

 96 Figure 9.6. Two-dimensional evolution of the char Σ at different times. Here β = 400, U = 0.9, ε = 0.01. After a transition period (a) to (e), the char pattern develops into chaotic fingering with tip-splitting (f). Compared with previous figures, collage of chaotic fingering patterns shown in (g), (h) and (i) exhibits a richer dynamics.

Figure 9 . 7 .

 97 Figure 9.7. Two-dimensional profile of the flame structure at t = 0.2 T 0 . (a) and (d): temperature, (b) and (e): flame front, (c) and (f ): char pattern, respectively for U = 0.2 (top) and 0.9 (bottom). Here β = 256 and ε = 0.01.The temperature is represented in the "fresh" side x ≤ φ(t, y), whereas the char pattern lies in the "burnt" side x ≥ φ(t, y).

Figure 9 . 8 .

 98 Figure 9.8. As in Fig. 9.7, two-dimensional profile of the flame structure at t = 0.5 T 0 . (a) and (d): temperature, (b) and (e): flame front, (c) and (f ): char pattern, respectively for U = 0.5 (top) and 0.9 (bottom).

Figure 9 . 9 .

 99 Figure 9.9. Char pattern at the flame front, for different time scales. Here β = 400, U = 0.9, ε = 0.01. Propagation from bottom to top. Lighter colors correspond to more advanced parts of the char at the flame front.

  Organization of the paper. The paper is divided into two parts, corresponding to two supplementary approaches and different fingering regimes. Part I goes from Section 2 to Section 6, and Part II from Section 7 to Section 9.

	2 exp(S).	(1.13)
	1.3.

Acknowledgment

The authors would like to thank Elisha Moses for providing us with Figure 1.1.

J. S. is supported in part by NSF of China grants 91130002 and 11371298. G.I. S. is supported in part by the US-Israel Binational Science Foundation (grant 2012-057) and the Israel Science Foundation (grant 335/13). C.-M. B. wishes to thank the State Administration of Foreign Experts Affairs, PRC.