Hydrodynamic bidimensional stability of detonation wave solutions for reactive mixtures

The structure of a planar detonation wave is analyzed for an Eulerian mixture of ideal gases undergoing the symmetric reversible explosive reaction . The chemical rate law is derived from the reactive Boltzmann equation, showing a detailed chemical kinetics in terms of a second-order reaction rate. The hydrodynamic bidimensional stability of the detonation wave is also investigated using a normal mode approach, when small time-space transverse disturbances affect the shock wave location. A suitable numerical technique is here proposed in order to solve the stability problem and numerical results are provided illustrating the detonation wave structure and its instability spectrum.


Introduction
The steady detonation for chemically reacting gases has been investigated in the framework of the macroscopic equations deduced in the hydrodynamic limit of kinetic models. This approach has been considered in several works. See, for example, papers [1][2][3]. Furthermore, the linear stability of the steady detonation, when unidimensional disturbances perturb the steady wave, has been analised in the kinetic frame in papers [4][5][6].
On the other hand, the detonation stability in the presence of bidimensional perturbations has been considered in the kinetic frame for the first time in paper [7]. In that paper, the stability equations have been constructed adopting a normal mode approach and the necessary closure condition has been derived as a radiation condition at the end of the reaction zone, according to the classical approach, see for example paper [8].
The stability analysis presented in paper [7] provides a consistent mathematical modelling with a rather detailed chemical kinetics, so that the analysis is highly https://doi.org/10.1088/1742-5468/ab3341 J. Stat. Mech. (2019) 083217 sensitive to the eects of the chemical reaction. Such result is mainly due to the fact that the dispersion relation includes the contribution of several chemical kinetic parameters, since the chemical reaction rate is built in the hydrodynamic limit of a molecular dynamics within kinetic theory of reacting gases. At the same time, the model is suitable to develop the computational treatment of the detonation stability.
Therefore, starting from the mathematical modelling of paper [7], for what concerns the detonation dynamics and its stability properties, a pertinent solution technique is proposed here in order to solve numerically the stability problem and to characterize the instability spectrum. The technique extends the one proposed in paper [6], where unidimensional perturbations were considered, to the case of bidimensional transverse disturbances. The numerical strategy combines the shooting method by Lee and Stewart (see [8]) with the Erpenbeck's idea of counting the zeros of a complex function (see [9]). However, the fact that bidimensional perturbations are taken into account renders the stability problem much more complete but, at the same time, rather dicult to be treated from both theoretical and computational sides. Therefore, additional eorts are required to set up the numerical technique and to provide satisfactory results at the end of characterizing the instability spectrum of the solution in an appropriate parametric space.
In the present paper, a numerical application is treated for a chosen set of input data, such as specific heat ratio, overdrive degree, reaction heat and activation energy, for which at least one instability mode does exist, as shown by the contour plot of the residual function. In particular, the influence of the disturbance wave number in the instability spectrum is studied, representing the disturbance growth rate and disturbance frequency of the instability mode for dierent values of the wave number.
The results presented in this paper can be viewed as the first computational treatment of the stability problem with bidimensional perturbations formulated in paper [7]. At the same time, the numerical application developed in this paper validates the extension of the numerical technique to the more complex problem of bidimensional stability.
The content of this paper is organized as follows. After this introduction, the principal aspects of the modelling are summarized in section 2. The mathematical description of the steady detonation wave is revisited in section 3 with reference to paper [7]. The problem of the hydrodynamic bidimensional stability of the steady detonation wave is formulated in section 4, with emphasis on its closure radiation condition and on the resolubility strategy in section 4.1. In particular, the numerical technique adopted to solve the stability problem is presented in section 4.2. At last, a numerical application in an appropriate parametric space is treated in section 5. The procedure presented in the previous sections is applied in view of producing some representative results for what concerns the steady detonation wave and its related stability problem.

Model equations
The model equations are those for an Eulerian reacting mixture undergoing the reversible chemical reaction A 1 + A 1 A 2 + A 2 . Such equations express the balance laws for the mass density ρ, flow velocity u = (u 1 , u 2 ) and pressure p of the mixture, and for the J. Stat. Mech. (2019) 083217 progress variable of the chemical reaction defined as the mass concentration z of the products. They are given by where e = ε + u 2 /2 is the total specific energy of the mixture, ε = ε(ρ, p, z) is the specific internal energy of the mixture, and r = r(ρ, p, z) is the reaction rate. The closure of the system (1)-(5) is assured by the thermal and caloric equations of state, where k is the Boltzmann constant, m the constituent's molecular mass, T the mixture temperature and γ the specific heat ratio, and by the reaction rate law where a is the molecular diameter, f and r are the forward and backward activation energies, and ∆ = 1 − 2 is the binding energy dierence. Such reaction rate has been obtained in the kinetic theory framework for chemically reactive mixtures (see book [10]), using an approximate solution of the microscopic system of Boltzmann equations (BE) describing the reactive mixture in terms of a collisional molecular dynamics. Adopting reactive dierential cross sections of hard spheres with activation energy, and assuming a chemical flow regime close to the thermodynamical equilibrium, an approximate solution of the BE has been obtained as Maxwellian velocity distribution functions which does not assure the chemical equilibrium condition, that is where n α is the particle number density and C α the modulus of the peculiar velocity of each constituent. The chemical equilibrium condition is https://doi.org/10.1088/1742-5468/ab3341 where is the reaction heat.

Detonation wave structure
The structure of the detonation wave characterized here is consistent with the classical ZND model, proposed by Zeldovich, von Neumann and Doering, for a one-dimensional detonation wave. The configuration of the ZND wave consists of a planar, non-reactive shock wave propagating with constant velocity in the positive x-direction, followed by a finite reaction zone where the chemical reaction takes place, see [11]. The reaction zone connects the von Neumann state N just behind the shock wave, where the chemical reaction is triggered, to the final state S, where the reaction reaches the chemical equilibrium. The ZND configuration is steady in the shock-attached frame, so that the steady variable ξ is introduced in the form with ξ < 0 ahead of the shock and ξ > 0 behind the shock. The governing equations are the reactive Euler equations (1)-(5) in one-dimensional form, where u represents the flow velocity of the mixture in the x-direction. The ZND structure of the steady detonation wave solution is specified when the von Neumann state N, the equilibrium final state S and the continuous reacting flow in the reaction zone are determined for each value of the detonation velocity s. The N and S states are determined resorting to a Rankine-Hugoniot analysis, so that the existence of a minimum acceptable value of the wave speed s can be proved. Such minimum defines the Chapman-Jouguet velocity, s j , and the corresponding equilibrium final state is the CJ state. Since the present detonation model considers a reversible reaction, the Hugoniot diagram must be completed with the so called equilibrium Hugoniot curve, locus of all equilibrium final states (see, for instance [11,12]).
In paper [7], the equations of such curve have been explicitly derived, and a solution procedure has been proposed in order to determine the steady detonation wave structure for the present model. Accordingly, the wave structure is characterized as follows.

Final states
For each value of s, the equilibrium final states in the v-p plane, with v = 1/ρ being the specific volume of the mixture, are characterized by where z is obtained by the chemical equilibrium condition with Above the upper sign identifies a state compatible with a strong detonation wave solution, whereas the lower sign identifies a state compatible with a weak detonation wave solution. In this paper, the analysis is addressed to strong detonation solutions.

CJ state
The CJ state is obtained when the wave speed s reaches the minimum value s j for which the corresponding Rayleigh line is tangent to the equilibrium Hugoniot curve. Therefore, β(z) = 0 and s j and z CJ eq are determined by solving conditions and 1 + s 2 γ + 1 The CJ state is then characterized by

Overdrive degree
For the admissible values of the detonation wave speed, s s j , the overdrive degree f is defined by f = (s/s j ) 2 . Thus, f = 1 corresponds to the Chapman-Jouguet detonation whereas f > 1 to the overdriven detonation.

Von Neumann state
The von Neumann state N is the non trivial solution of the Rankine-Hugoniot conditions of the model referred to the upper state of the reaction zone where the chemical reaction is not yet initiated. Such post-shock state is characterized by where w = s − u is the steady waveframe velocity, and by the condition

Continuos reacting flow
The continuous flow within the reaction zone is determined by solving (for ξ > 0) the following dimensionless system which is deduced by referring the governing equations (12)- (15) to the steady variable (11) and by considering p , T, w and z as variables, In the above equations, p , T, w, ∆ , ξ and r are given, respectively, in units of p 0 , T 0 , is a reference time of order of mean free time. The dimensionless form of the reaction rate r reads https://doi.org/10.1088/1742-5468/ab3341 with f given in units of kT 0 . The initial conditions at ξ = 0 for the system of equations (25)-(28) are provided by the post-shock conditions (23).

Graphic representation of the solution
The structure of the detonation solution can be represented in the Hugoniot diagram including Rayleigh lines, partial Hugoniot curves, and also the equilibrium Hugoniot curve defined in parametric form by the equations (16) and (17), where z is the parameter obeying the equilibrium condition (18).
As an illustrative example of the detonation wave structure, figure 1 shows the Hugoniot diagram constructed using the present model, for the following choice of the material properties For this choice, the CJ velocity s j and the corresponding equilibrium concentration z j are The CJ Rayleigh line is drawn in figure 1 with another line for s = 9.007 24. The figure contains the equilibrium Hugoniot curve, gathering all final states for dierent values of the wave velocity s s j . Such curve contains the final states proper of the strong solution and, for completeness, also those states in the weak branch. The Hugoniot curve for z = 0 represents all von Neumann states in the strong branch and other post-shock states in the weak branch.
For completeness, figure 1 also contains the reference Hugoniot curve for z = 1 obtained as a particular hyperbola of the family of the Hugoniot curves of fixed product concentration z. The points on this curve do not belong to the wave solution, because the reversibility of the chemical reaction is not compatible with a state for which all reactants of the forward reaction are transformed into products, namely z = 1.
Further aspects of the detonation wave structure will be investigated in section 5, where some numerical simulations are performed regarding both the steady detonation solution and its bidimensional hydrodynamical stability.

Hydrodynamic stability
The hydrodynamic linear stability of the steady wave structure is studied in presence of bidimensional disturbances. A small rear boundary perturbation is instantaneously assigned and a distortion on the shock wave location, ψ(y, t), is induced by small transverse disturbances. Assuming that the instability of the detonation solution results uniquely from the interaction between the perturbed shock and the reaction zone, the https: In what follows, the stability problem is formulated assuming a normal mode representation of the disturbances, and a numerically technique is proposed to solve the problem.

Mathematical formulation of the stability problem
The analysis is based on the bidimensional reactive Euler equations (1)-(5), re-written in the form and then transformed to the perturbed shock-attached frame (35) A restricted class of asymptotic solutions to the transformed system, which deviate by a small amount from the known steady detonation solution, is seeked. Accordingly, a linear approach is appropriate to analyze this problem and equations (31)-(34) are linearized through a normal mode expansion around the known solution, represent the state vector, the known steady-state vector and the complex perturbation amplitude vector, respectively. Moreover, Re(α) is the disturbance growth rate, Im(α) the disturbance frequency, k the disturbance wave number, and ψ the spatial perturbation amplitude of the shock location. The linearized equations for the perturbation amplitudes q results then in the following homogeneous linear system constituting the stability equations of the model, Above, α is given in units of τ −1 0 , k in units of (sτ 0 ) −1 and ξ in units of sτ 0 , where τ 0 is the reference time defined by expression (29). The linearization of the reaction rate leads to the form where r * v , r * p and r * z denote the partial derivatives of r. The stability equation (38) constitute a set of ten first-order linear dierential equations for the real and imaginary parts of the complex perturbations, with spatially varying coecients depending on the steady solution q * (ξ). https://doi.org/10.1088/1742-5468/ab3341

J. Stat. Mech. (2019) 083217
Linearization of the bidimensional Rankine-Hugoniot relations at the shock, ξ = 0, leads to the following boundary conditions to be joined to the stability equation (38), The stability equations and their boundary conditions, besides the unknown state variables perturbations, involve the shock spatial perturbation amplitude ψ , as well as the key parameter α which determines the stability character of the steady wave. The amplitude ψ will be normalized with respect to a reference amplitude as detailed in section 5. Therefore a closure condition is needed in order to determine the dynamics of α, and a perturbation acoustic analysis is developed at the end of the reaction zone, where the state variables reach their equilibrium values and remain constant.
The closure condition is obtained following the procedure proposed in [13] and then adopted in several papers, see, for example, [8,14,15].
More specifically, first one assumes that the spatial amplitudes q are represented as a Fourier mode in terms of the wave number k, that is q = qe ikξ , and the dierential system of stability equation (38) reduces to an algebraic system of homogeneous linear equations whose solutions describe dierent modes of propagation, namely vorticity and entropic waves, backward and forward acoustic waves and a mode associated to the chemical reaction. This implies that the spatial perturbation amplitudes q can be expressed as a superposition of acoustic waves propagating at the characteristic speeds associated to such modes Then, the crucial point in this derivation is a consequence of a further condition proper of the stability analysis [14], usually referred to as radiation condition in the acoustic nomenclature. It is assumed that the instability depends only on the interaction between the perturbed shock and the reaction zone so that q does not depend on the forward family of acoustic waves.
Mathematically, the procedure leads to the explicit derivation of the closure condition for α, which represents the dispersion relation of the normal modes (36) and has been deduced in [7] in the form https://doi.org/10.1088/1742-5468/ab3341 The dispersion relation (47) shows a relevant influence of the chemical kinetics, thanks to the presence of the partial derivatives of the reaction rate. These contributions are due to the detailed form of the chemical reaction rate which is built in the hydrodynamic limit of the considered kinetic model.
On the other hand, one can observe that in several well known papers in the classical detonation literature see, for example, papers [14,16] and related bibliography, where a first-order reaction rate of Arrhenius form is usually adopted, the radiation condition has a reduced form with respect to condition (47). In fact, such form is recovered as a particular case of condition (47), when the derivatives r * v and r * p vanish at the end of the reaction zone.
In conclusion, the stability problem is here formulated in terms of the complex perturbation amplitudes q and complex growth rate α, by means of the stability equation (38) for ξ ∈ ]0, ξ eq [, with boundary conditions (42)-(46) at ξ = 0, and closure condition specified by the dispersion relation (47) at ξ = ξ eq . Conversely, the disturbance wave number k remains as a real parameter.
This problem presents a practical diculty for the following reasons: the integration of the dierential equation (38) with boundary conditions (42)-(46) needs the knowledge of α; on the other hand, α must be specified by the dispersion relation (47) which, in turn, involves the unknown amplitudes q specified at ξ = ξ eq . Thus, the integration of the stability equations requires the knowledge of α but, at the same time, the condition that determines α requires the solution of the stability equations. To overcome such diculty, a suitable solution technique should be used. A further diculty is the fact that the coecients of the stability equations are not constant. Therefore, a numerical approach seems to be the most convenient strategy to determine the stability solution and to investigate the influence of the model parameters on the instability spectrum.

Solution technique
The numerical technique adopted in this paper follows the procedure proposed in [6], for one-dimensional linear stability and dierent chemical rate. The technique combines the iterative shooting method first proposed by Lee and Stewart in paper [8], hereinafter referred to as the LS procedure, with the Cauchy's argument principle used by Erpenbeck in paper [9], and is addressed to determine the instability modes which correspond to a positive growth rate Re(α). Since these modes occur in conjugate pairs,

J. Stat. Mech. (2019) 083217
a region R in the upper-right quarter of the complex plane is considered in order to search appropriate values of the perturbation parameter α.
The shooting technique starts with a trial value for α and the stability equations are solved for that α. Then one tests if the solution to the stability equations and corresponding trial α verify the dispersion relation. If this is the case, then α and q constitute a solution to the stability problem. Since, in general, this is not the case, one has to consider a new trial value for α and iterate the procedure until the dispersion relation is verified.
A straightforward application of the LS procedure in the eigenvalue domain of the complex plan requires a discretization of that domain using a very thin grid. This procedure requires a large number of trial values, see [17].
A guide for the search of appropriate trial values for α, with lower computational eort, is provided by the residual function H(α), defined by the expression on the left hand side of equation (47), whose zeros are the solutions to the dispersion relation.
The numerical technique consists then of the following steps.
(i) Consider a tentative region R in the upper-right complex plane, and select a great number of points a j , j = 1, 2, . . . , n, in its contour. For each a j , j = 1, 2, . . . , n, fix another point b j such that Re(b j ) = Re(a j ) + 10 −6 and Im(b j ) = Im(a j ). (iv) Use the argument principle to determine the number Z of zeros of H inside the region R. Since H has no poles in R, such number is given by where ζ : [k, ] → C is a path smooth by parts, describing the contour of R in the positive direction. The derivative H (ζ(t)) is estimated by the ratio If Z = 0, the stability problem does not admit any solution in the region R, and the procedure must be iterated starting from a dierent region.
If Z > 0, the stability problem admits at least one solution in the region R, and a contour plot of |H| can be drawn in order to approximate the location of the instability solutions. After that, a refinement of R is needed around each zero and the procedure https://doi.org/10.1088/1742-5468/ab3341 J. Stat. Mech. (2019) 083217 must be iterated, starting from the corresponding refinement, until the solution satisfies the required precision.

Numerical results
Some representative results for the structure of the steady detonation solution and its instability spectrum are presented in sections 5.1 and 5.2, with reference to the data specified below. Then, some computational considerations are presented and discussed in section 5.3.
The reference parameters in this analysis are the specific heat ratio γ, overdrive degree f , forward activation energy f , reaction heat Q and wave number k. In view of equation (10), the binding energy dierence ∆ can be considered in place of Q. The numerical computations presented in this paper do not cover the variation of all these parameters, but only some of them. In particular, the structure of the steady detonation solution is obtained for two values of the overdrive degree f whereas the stability spectrum is characterized varying the wave number k and the forward activation energy f .
For the numerical application, the selected input data are the dimensionless preshock state of the reactive flow, and the material properties

Structure of the steady detonation solution
The numerical simulations refer to two dierent values of the shock wave velocity, namely s = s j and s = 3.184 54 > s j , that is f = 1 (Chapman-Jouguet detonation) and f = √ 2 (overdriven detonation). In this analysis, the forward activation energy is held fixed at f = 2. The detonation wave structure for the considered reacting gas mixture is represented in figure 1 in terms of the Hugoniot diagram and in figure 2 in terms of state variables profiles. These profiles show the behaviour in the reaction zone of the pressure, temperature, waveframe velocity and progress variable, in dependence of the distance ξ from the shock front. Also observe that the pressure profiles in the upper-left pictures of figure 2 exhibits a rarefaction in the reaction zone, reproducing the typical ZND behaviour.

Instability spectrum
The instability spectrum of the steady detonation solution is investigated solving the stability problem formulated in section 4.1.
The analysis is performed for the steady solution obtained when f = √ 2 and the instability spectrum is investigated in the region R = [0, 0.5] × [0, 1] of the upper-right quadrant of the complex plane.
In particular, for the material properties (50), overdrive degree f = √ 2 and wave number k = 0.18, there exist two instability modes in the region R, which are given by α = 0.307 + 0.24i and α = 0.354 + 0.128i. The contour plot of the residual function |H (α)| represented in figure 3 in a subregion of R shows the domain of attraction of such instability modes. The mode with lowest frequency, i.e. the one with smaller imaginary part, represents the fundamental mode, which has a particular interest in the stability analysis [8,15,17]. Observe that, in presence of bidimensional perturbations, the behaviour of this mode in dependence of the wave number of the disturbances applied to the steady detonation front plays a relevant role in the stability analysis.

5.2.1.
Influence of the wave number on the fundamental mode. The first aspect to be investigated is how the fundamental mode is aected when the activation energy is maintained fixed at f = 2 and the wave number is varying in the interval [0, 2]. The migration of the fundamental mode is followed in the complex region R, and the diagrams of figure 4 show the disturbance growth rate Re(α) and disturbance frequency Im(α) versus the wave number k. The pictures reveal that for increasing values of k, the fundamental mode becomes less unstable, since its growth rate (left frame) tends to zero. One may also notice that the frequency of the instability mode (right frame) is zero for k = 0, grows to 0.12 for k = 0.18 and rapidly decreases for larger values of k.  Conversely, it is purely imaginary for f = 10 when k = 0.3, For each value of the wave number k < 0.5, the fundamental mode has larger real part for the larger value of the forward activation energy f . This implies that an increasing of the activation energy, when the wave number is fixed at k < 0.5, intensifies the instability character of the fundamental mode. For f = 10 and k < 0.3, the fundamental mode becomes stable, since its real part becomes negative. On the other hand, the fundamental mode no longer exists when k > 1.5. These features are in agreement with other results presented in literature, see [15], for example.

Computational considerations
Some aspects on the computational eorts required to determine the instability spectrum are discussed. In particular, some comparisons between the numerical technique proposed in section 4.2 and the LS procedure advanced by Lee and Stewart in [8] are here analysed.
The numerical study of the linear stability problem is computationally demanding, essentially because, for each set of chosen parameters, it requires a large number of trial values. In fact, equations (38) with their boundary conditions (42)-(46) have   to be integrated a huge number of times until the dispersion relation (47) is satisfied or, equivalently, until the zeros of the residual function H(α) are found. In addition, for each set of chosen parameters, a root finding routine has to be implemented in an appropriate refinement of R in order to approximate the location of the corresponding instability modes. A further complexity is that the stability problem is very sensitive to small changes in the parameter space, so that in order to study the influence of one selected parameter on the instability spectrum, the procedure must be iterated for very small variations of that parameter.
• With reference to the contour plot of figure 3 and the corresponding set of parameters, 220 trial values were needed to determine the number of instability modes presented in the region R. Then, four sequential refinements of the region R were considered and a further 480 trial values were added to determine a small subregion of R containing the two modes. To guarantee a precision of order of 10 −2 , other 500 trial values were required in the resolutive procedure. This means that the computational procedure needed a total of 1200 trial values to obtain the contour plot represented in figure 3.
• In order to obtain the migration of the fundamental mode represented in figure 4, when the wave number k varies in the interval [0,2] with increments of order of 10 −2 , the stability problem should be numerically solved for 200 dierent sets of parameters, each one requiring 1200 trial values for α.
• The results shown in table 1 have been obtained solving the stability problem for 27 dierent sets of parameters, with increments of k of order of 10 −2 .
In view of a comparison between the numerical technique proposed here and the LS procedure, at least 5000 trial values in place of 1200 are needed for the latter method when applied to the same region R with a precision of the same order. In fact, without any previous knowledge on the location of the modes, the LS procedure uses a grid of trial values to draw a contour plot in the whole region, corresponding to the carpet search [8]. The computations for each trial value require about three seconds using the resources employed in this paper, namely a computer with a 2.4 GHz intel Core I5 processor, 8 GB RAM memory and 64 bits system operator. Therefore, the LS procedure requires a greater computational time, which can be estimated in an additional 3 h to obtain the contour plot of figure 3. Although the computational time depends on the specific values of the considered parameters, one could expect that the LS procedure requires an additional 600 h for figure 4 and an additional 81 h for table 1.
Observe that the dierence between the computational times would be even more significant when two instability modes exist in R at a distance smaller than 10 −2 . In this case, the precision of order 10 −2 would not be enough to distinguish the two modes and more trial values would be needed.
A further aspect to be underlined is that the numerical technique proposed in this paper starts by counting the number of modes in region R. This is not the case if the LS procedure is straightforward applied in R, since in this procedure the modes are not counted. Therefore, it is not possible to decide a priori if a thinner grid will be needed to distinguish two modes. https://doi.org/10.1088/1742-5468/ab3341 J. Stat. Mech. (2019) 083217

Conclusions
In this work some effort is done to inspect the complex problem of the instability spectrum of the steady detonation wave for an Eulerian binary gas mixture with a reversible chemical reaction and bidimensional transverse perturbations. The modelling arises from a detailed kinetic description of the chemical mechanism, and the treatment is to be considered in the hydrodynamical limit and in flow conditions close to equilibrium. At this scope, the main arguments discussed by some of the authors in two quite recent papers (see [6] and [7]) are revisited for reader convenience and widened in view of the numerical application proposed in this paper. The methodology presented in section 4 is applied for the actual calculation of the instability spectrum. In particular, a renewed use of the equilibrium properties investigated in paper [7] is here employed to characterize the steady solution for an overdriven one dimensional detonation wave. Moreover, a new numerical technique is here proposed to explore the instability spectrum through a two dimensional linear stability analysis, on the basis of the ideas presented in [6] for the one-dimensional linear stability. The results of section 5 seem to provide an accurate and detailed consistent picture of the detonation wave and its instability, specifically due to the kinetic features of the chemical device. Nevertheless, it is convenient to underline that the results obtained through a linear stability analysis cannot be directly compared with physical experiments where the non linear eects play a dominant role. Thus, on the basis of this first attempt to the analysis of the instability spectrum, some future issues could be addressed to enlarge and improve the modelling here proposed, in view of better understanding the detonation instability and the physical mechanism behind it.
A particularly interesting extension of the present study could be the one of considering a dierent reaction scheme, where direct and reverse chemical reactions have a dierent order. Such scheme would induce a significant dependence of the reaction kinetics on the local thermodynamic state in the reaction zone. The starting point would be a dierent chemical mechanism such that the number of reactants diers from that of products, for instance a dissociation-recombination process. The corresponding kinetic model would lead to a reaction rate formally similar to the one in equation (7) but with dierent order non-linearities in the reactant's and product's concentrations.
Another possible extension could be that of going beyond the Euler level up to the Navier-Stokes regime [18,19] and construct a more detailed modelling able to capture significant properties of detonation and its bidimensional stability. More in particular, the Navier-Stokes equations can be used to investigate and understand the eects of boundary layers on the detonation structure, cellular detonation patterns and their formation when a large number of two-dimensional unstable modes appear, pathological detonation and its stability [20][21][22][23][24]. On the other hand, the macroscopic equations derived in the hydrodynamic limit of a kinetic model contain explicit and detailed expressions of the reaction rate and transport coecients which are deduced using the approximate solution to the reactive Boltzmann equation, see, for instance, the book [10] and related bibliography. Therefore, the hydrodynamic equations obtained in this way show the capability to provide a good description, at least from the theoretical and numerical point of view, of detonation processes in which the coupling between chemical kinetics and the flow evolution represents a dicult task [25,26].

J. Stat. Mech. (2019) 083217
At last, another issue could be the extension of the present stability analysis to three dimensional linear perturbations, using the normal mode approach employed in the present paper. Three dimensional perturbations are relevant, for example, in the stability of spinning detonation, since the spinning regime can describe, also experimentally, a cellular detonation, as documented in paper [16].