text:
# Quantum Melting of the Charge Density Wave State in 1T-TiSe2

## Abstract

We report a Raman scattering study of low-temperature, pressure-induced melting of the CDW phase of 1T-TiSe2. Our Raman scattering measurements reveal that the collapse of the CDW state occurs in three stages: (i) For P$`<`$5 kbar, the pressure dependence of the CDW amplitude mode energies and intensities are indicative of a “crystalline” CDW regime; (ii) for 5 $`<`$ P $`<`$ 25 kbar, there is a decrease in the CDW amplitude mode energies and intensities with increasing pressure that suggests a regime in which the CDW softens, and may decouple from the lattice; and (iii) for P$`>`$25 kbar, the absence of amplitude modes reveals a melted CDW regime.

There has been a great deal of interest in the relationship between various diverse and exotic low temperature phases of strongly correlated systems, including the antiferromagnetic insulating and unconventional superconducting phases of the high Tc cuprates,Ginsberg  the charge-ordered insulating and ferromagnetic metal phases of the manganese perovskites,Dagottobook  the orbital-ordered and ferromagnetic metal phases of the ruthenates,Cao2002 ; Nakatsuji ; SnowRuthenate  and the charge-density-wave (CDW) and superconducting phases of layered dichalcogenides such as 2H-NbSe2.Moncton  Of particular interest is the exotic phase behavior that is expected between fully-ordered (crystalline) and disordered (isotropic) phases as one tunes the interactions in these systems using some control parameter other than temperature. These include, for example, electronically phase-separated regimes,Dagottobook  and “quantum liquid crystal” smectic and nematic phases, which are expected to be observed between charge-ordered insulating and ‘disordered’ metallic or superconducting phases as one increases the interactions between the charge stripes.Kivelson  Clearly, therefore, it is of great interest to carefully explore the manner in which quantum ordered phases collapse, or ‘melt’, into quantum disordered phases as a function of some control parameter that tunes the competing interactions in the material at low temperatures.

In this paper, we report a pressure-dependent low temperature Raman scattering study of the CDW system 1T-TiSe2, in which we are able to explore the manner in which the CDW state ‘melts’ with increasing pressure near T$``$0 K. Because of its layered structure and simple commensurate CDW phase, 1T-TiSe2 is an ideal system for such an investigation. 1T-TiSe2 is also of interest because the CDW transition is not driven by conventional Fermi surface nesting, but rather by an unconventional mechanism involving electron-hole coupling and an “indirect” Jahn-Teller effect.Kidd  Our low-temperature, pressure-dependent light scattering approach allows us to explore unique details associated with quantum mechanical melting of the CDW in 1T-TiSe2. In particular, this study reveals that the CDW state evolves with increasing pressure in a manner reminiscent of classical 2D melting, with ‘crystalline’ and ‘disordered’ regimes, as well as an intermediate ‘soft’ CDW regime in which the CDW may be incommensurate with the lattice. Furthermore, this study demonstrates the efficacy of pressure-tuned light scattering for examining a much broader class of quantum phase transitions in strongly correlated systems.

The 1T-TiSe2 samples used in this study were grown by iodine vapor transport with a temperature gradient of 570-640 oC.Kidd  The sample stoichiometry was verified by x-ray and resistivity measurements. The Raman spectra were taken in a true backscattering geometry with 647.1 nm incident photons. Variable low temperature, high-pressure measurements were obtained with a modified SiC-anvil cell inserted into a flow-through helium cryostat, allowing continuous adjustment of both the temperature (3.5-300 K) and pressure (0-100 kbar).SnowRuthenate  Argon was used as the pressure transmitting medium, and the pressure inside the cell was determined from the shift of the ruby fluorescence line; argon is quasi-hydrostatic in the temperature and pressure range of interest.Venkateswaran  From the linewidth variation of the ruby line as a function of pressure and temperature, we estimate the degree of non-hydrostaticity in the pressure cell to be 1$`\%`$ or less.

Figure 1 shows the temperature-dependent Raman scattering spectrum below the CDW transition temperature at TC $``$ 200 K. Several new modes develop in the CDW phase. Of particular interest are an Eg mode near 75 cm-1 and an A1g mode near 115 cm-1. These two modes are CDW-coupled “amplitude” modes associated with the zone-boundary transverse acoustic phonons from the $`L`$-point in the Brillouin zone, which are folded to the zone-center due to the formation of the CDW superlattice.Holytise2 ; Sugai  That these two modes are indeed coupled to the CDW mode is confirmed by their temperature dependence in Fig. 2: Both the 75 cm-1 Eg and 115 cm-1 A1g amplitude mode energies soften as the temperature is increased toward the CDW transition temperature. By contrast, the Eg zone-center optical phonon near 134 cm-1 is nearly temperature-independent, indicating that it is not strongly influenced by the development of the CDW state.

A few comments concerning the two CDW amplitude modes observed in Fig. 1 will help in the interpretation of the pressure-dependent light scattering data. First, the CDW amplitude modes are excitations of the CDW ground state, involving fluctuations of the CDW state that modulate the amplitude of the charge density wave. In particular, the A1g amplitude mode near 115 cm-1 is a breathing mode distortion, and therefore involves fluctuations of the CDW amplitude that preserve the symmetry of the CDW ground state. Consequently, the A1g amplitude mode displacement pattern is the same as the static CDW displacement pattern shown in the inset of Fig. 1. On the other hand, the Eg amplitude mode involves out-of-phase fluctuations of the CDW amplitude away from the ground state symmetry. In combination, these two amplitude modes serve as ideal “probes” with which to study the stability and elasticity of the CDW state as it evolves, and eventually melts, as a function of increasing pressure.

Figure 3 shows the pressure dependent Raman spectra of TiSe2 at 3.5 K. The energy and intensity of the 115 cm-1 A1g and 75 cm-1 Eg CDW amplitude modes, along with the 134 cm-1 Eg optical phonon mode, are summarized in Fig. 4. One of the chief effects of pressure on the low temperature Raman spectrum of 1T-TiSe2 is the gradual suppression of the 75 cm-1 and 115 cm-1 CDW amplitude mode intensities with increasing pressure, and the complete collapse of the CDW state near a T$``$0 critical pressure of P$`{}_{}{}^{}`$25 kbar. This value of the critical pressure is similar to that observed in pressure-induced CDW-to-metal transitions measured using transport measurements, including NbSe3 (P$`{}_{}{}^{}`$24 kbar)Yasuzuka  and Lu5Ir4Si10 (P$`{}_{}{}^{}`$21 kbar).Lue

Even more interesting than the observed pressure-induced collapse of the CDW phase, is the manner in which this collapse occurs. This process can be carefully studied by examining the energies and integrated intensities of the 75 cm-1 Eg and 115 cm-1 A1g amplitude modes as a function of pressure, which are summarized in Fig. 4(b) and (c), respectively. For comparison, the pressure-dependent energy of the 134 cm-1 optical phonon mode is illustrated in Fig. 4(a).

Notably, Fig. 4 reveals several regimes of behavior associated with the pressure-induced (T$``$0) quantum melting of the CDW state in 1T-TiSe2: (i) Crystalline CDW regime \- From P=0 to 5 kbar, the A1g CDW amplitude mode’s intensity increases, and its energy increases at a rate of approximately $`d\omega _o/dP`$ \+ 1 cm-1/kbar. This behavior is consistent with increased stiffening of the CDW state with increased pressure, characteristic of a typical “crystalline” phase. Significantly, the rate at which the A1g amplitude mode energy increases with pressure is similar to that observed for the 134 cm-1 optical phonon mode (Fig. 4(a)), providing evidence that the CDW remains commensurate with the lattice in this regime. (ii) Soft CDW regime \- Between P = 5 to 25 kbar, the A1g amplitude mode exhibits a number of interesting and anomalous changes as a function of increasing pressure: its energy softens - revealing an anomalous, slightly negative Grüneisen mode parameter d$`\omega _o`$/d$`P`$ in this regime - its intensity decreases rapidly, and its linewidth increases substantially. All of this behavior is consistent with a distinct softening of the CDW state to breathing-mode fluctuations of the CDW amplitude in this regime. Equally interesting is the fact that the energies of the A1g amplitude and Eg optical phonon modes exhibit distinctly different pressure-dependencies in this regime, suggesting that the CDW decouples from the lattice (i.e., becomes incommensurate) in this regime.commensurate  Evidence for a pressure-induced transition to an incommensurate CDW phase has also been observed in 2H-TaSe2 by high-pressure x-ray scattering.McWhan  Interestingly, the Eg amplitude mode also exhibits a particularly anomalous pressure-dependence in this regime, including a decrease in energy and a rapid reduction in intensity with increasing pressure. The rapid disappearance of the Eg mode, in particular, indicates that out-of-phase fluctuations of the CDW amplitude are not well-defined excitations above roughly P$``$5 kbar, even though there is clearly some vestige of the CDW state at these pressures, as evidenced by the presence of the A1g amplitude mode. This may indicate a breakdown of long-range CDW order in this phase regime; indeed, calculations of 1D CDW systems show that a softening of the CDW to shear deformations, and a consequent breakdown of long-range translational order, occurs when the coupling between CDW stripes reaches a critical value.ZaitsevZotov  In sum, the observed pressure-dependencies of the CDW and lattice modes are characteristic of a regime in which the CDW has not only begun to soften, or ‘melt’, but in which the CDW may have an increased propensity to move (‘slide’) in response to an applied field due to decreased pinning to the lattice. (iii) Disordered regime \- Finally, above P$`{}_{}{}^{}`$ 25 kbar, there is no evidence in the spectra for long- or short-range CDW order, suggesting that the CDW state has melted completely into a disordered phase.

It is interesting to compare the pressure-induced ‘melting’ process described above for 1T-TiSe2 to melting in other systems.HalperinNelson ; Dai  For example, calculations of classical melting in 2D solids suggest the presence of a “hexatic” phase - in which long-range orientational order is preserved, but long-range translational order is lost - intermediate between the crystalline and disordered phases. This topological phase arises because dislocations cause translational order to decay exponentially, but cause a much weaker suppression of orientational order. Such a hexatic-like phase has indeed been observed as a function of increasing disorder ($`x`$) in the NbxTa1-xS2 system using scanning-tunnelling-microscopy (STM): these measurements reveal a system that is, in turn, crystalline (0 $`<x<`$ 0.04), hexatic (0.04 $`x`$ 0.07), and amorphous ($`x>`$ 0.07), as a function of increasing disorder ($`x`$).

While the above phase behavior in 2D systems is reminiscent of the phase regimes we observe as a function of pressure in 1T-TiSe2, two important points should be made with respect to this comparison. First, in contrast to the examples above, the melting process we observe in 1T-TiSe2 is driven at T$``$0 K by pressure-tuning the competing interactions in this system. To understand the nature of this competition, note first that the zero-pressure CDW state in 1T-TiSe2 is unconventional, arising from an indirect Jahn-Teller-type interaction that splits and lowers the unoccupied conduction band:Kidd  as a result of the electron-hole interaction between the conduction and valence bands, the lowering of the split conduction band “repulses” and flattens the valence band, resulting in a lowering of the system’s energy, and the formation of a small gap CDW state. Upon applying pressure to this CDW state, one expects several regimes of behavior: at low applied pressures, increasing pressure will increase the matrix element associated with the Jahn-Teller interaction; this will result in a further lowering the conduction band, and via the electron-hole coupling, a lowering of the valence band and a consequent stiffening of the CDW state. This behavior is similar to that observed in the “crystalline CDW” regime of Fig. 4. As the pressure is increased beyond some critical pressure, however, the increasing strength of the Jahn-Teller interaction will eventually overwhelm the electron-hole interaction between the conduction and valence bands, leading to a collapse of the CDW gap, and a pressure-induced transition to a metallic phase in which the CDW distortion is completely suppressed. Most significantly, our results indicate that prior to the complete collapse of the CDW gap above P, there is a distinct “soft CDW” phase regime in which the CDW loses its stiffness; this is likely to result from an increase in fluctuations of the CDW near P. In theoretical support of this, Zaitsev-Zotov et al. have shown that increased coupling between CDW stripes leads to an increase in dynamic fluctuations and a decreased CDW stiffness.ZaitsevZotov  Interestingly, such long wavelength lattice fluctuations are expected to destroy long-range translational order but preserve long-range orientational order.Mermin

Second, the degree of interlayer coupling in 1T-TiSe2 should influence the nature of the quantum melting transition.mcmillan  For example, while classical melting of a 2D system occurs via a continuous phase transition with an intermediate topological phase regime, classical 3D melting occurs via an abrupt first-order phase transition. Significantly, the continuous nature of the pressure-induced evolution of the CDW state we observe in 1T-TiSe2, and evidence for an intermediate “soft CDW” regime in which there may only be short-range CDW order, suggests that pressure-induced quantum melting of the CDW in 1T-TiSe2 occurs in a manner akin to classical 2D melting.

In summary, Raman scattering studies of T$``$0 pressure-induced melting of the CDW state in 1T-TiSe2 reveals both a low-pressure (P$`<`$5 kbar) “crystalline” regime, and a high-pressure (P$`>`$25 kbar) “disordered” regime. Most interesting, however, is evidence for an intermediate pressure regime (5$``$P$``$25 kbar) in which the CDW loses its stiffness and may have no long-range translational order.

We thank M.V. Klein and E. Fradkin for useful discussions. This work was supported by the National Science Foundation under grant DMR02-44502, and by the Department of Energy under grant DEFG02-91ER45439.

Corresponding author: s\_cooper@mrl.uiuc.edu
+ Present address: Brookhaven National Laboratory


format: text

text:
# Large Scale Evolution of Premixed Flames

## I introduction

Propagation of plane flames in gaseous mixtures is well known to be unstable. An efficient way of investigating this instability is to consider the flame front as a surface of discontinuity, expanding all quantities of interest in powers of $`L_\mathrm{f}/\lambda ,`$ where $`L_\mathrm{f}`$ is the flame front thickness, and $`\lambda `$ characteristic length scale of a flame perturbation. The leading term of the perturbation growth rate expansion has the form

$`\sigma =c_0{\displaystyle \frac{U_\mathrm{f}}{\lambda }}`$ (1)
where $`U_\mathrm{f}`$ denotes an adiabatic velocity of a plane flame front with respect to the fuel, and $`c_0=c_0(\theta )`$ a function of the gas expansion coefficient $`\theta `$ defined as the ratio of the fuel density ($`\rho _\mathrm{u}`$) and the density of burnt matter ($`\rho _\mathrm{b}`$), $`\theta =\rho _\mathrm{u}/\rho _\mathrm{b}>1.`$ According to Refs. landau ; darrieus , $`c_0(\theta )`$ has positive values for all $`\theta ,`$ implying an unconditional instability of zero-thickness flames, the Landau-Darrieus (LD) instability. In the next order in $`L_\mathrm{f}/\lambda ,`$ account of the transport processes inside the flame front modifies Eq. (1) to

$`\sigma =c_0{\displaystyle \frac{U_\mathrm{f}}{\lambda }}\left(1c_1{\displaystyle \frac{L_\mathrm{f}}{\lambda }}\right),`$ (2)
where $`c_1`$ depends on $`\theta `$ as well as on the ratio of the heat and mass diffusivities (the Lewis number) markstein ; pelce ; matalon . The product $`c_1L_\mathrm{f}\lambda _\mathrm{c},`$ the so-called cut-off wavelength, is the short wavelength limit of unstable perturbations. By the order of magnitude, $`\lambda _\mathrm{c}`$ represents also the characteristic length of the so-called cellular structure of the flame front, which is formed eventually as a result of the nonlinear flame stabilization. For many flames of practical interest, $`c_1=1520.`$ It is the fact that $`L_\mathrm{f}`$ is relatively small in comparison with $`\lambda _\mathrm{c}`$ which underlies the above point of view on flame dynamics.

Consider an arbitrary initially smooth front configuration. As a result of the rapid growth of the unstable flame perturbations with wavelengths $`\lambda _\mathrm{c},`$ the flame front becomes corrugated within the time interval

$`\mathrm{\Delta }t{\displaystyle \frac{\lambda _\mathrm{c}}{U_\mathrm{f}}}.`$ (3)
Dynamics of the short wavelength modes are mainly determined by the transport processes inside the flame front, and are affected only slightly by the large scale flow. One can say that the small scale cellular structure of the flame front develops on the “background” of its smooth large scale configuration (see Fig. 1). It follows all developments of the background, Eq. (3) playing the role of the characteristic time of cell adaptation to the large scale front evolution.

In practice, it is the large scale evolution of the flame front, rather than its exact local structure, which is often of the main concern. In this respect, an important question arises about the reverse influence of the flame cellular structure on the front evolution at scales much larger than $`\lambda _\mathrm{c}.`$ More precisely, one can state the problem as follows. Imagine that we have smeared the small scale rapid variations of all relevant quantities by averaging them over many flame cells. Then the question is what equations governing dynamics of the averages are.

In connection with the above statement of the problem it should be noted that the exact cellular structure of flames is actually unknown. This is because the process of cell formation is essentially nonlinear, in the sense that it cannot be treated perturbatively in principle, which is the main reason of lack of its theoretical description. Only in the case $`\theta 1,`$ which is practically irrelevant, can this structure be determined analytically henon ; siv1 ; sivclav ; kazakov1 ; kazakov2 . The question of principle, therefore, is to what extent the large scale dynamics of averages depend on the exact local flame structure in the regime of fully developed LD-instability.

The main purpose of the present paper is to show that to the leading order in the ratio of $`\lambda _\mathrm{c}`$ to the characteristic length of the problem ($`L_0`$), dynamics of the averaged quantities are actually independent of particularities of the local flame structure. The latter determines essentially only one parameter characterizing the large scale evolution – the effective normal velocity of the flame front.

Perhaps, it is worth to explain the essence of the problem in a little bit more detail. The above point of view on the flame propagation is based on the possibility to separate the local cellular dynamics from the large scale evolution of the background. This possibility is underlined by the following common property of the transport processes. From the mathematical point of view, all these processes are of higher differential order than those governing dynamics of an ideal fluid. Therefore, their relative role increases at smaller scales. In particular, in the limit $`L_\mathrm{f}0,`$ cell formation is completely determined by the transport processes inside the flame front. On the contrary, the role of these processes at scales $`LL_\mathrm{f}`$ is relatively small. It should be fully realized, however, that this reasoning is inherently linear. It tacitly assumes that if every quantity of interest, say $`A,`$ is represented as a sum of its averaged value $`AA_0`$ and the small scale fluctuation $`A_1,`$ then the dynamics of $`A_0`$’s can be determined solely in terms of $`A_0`$’s themselves. Because of the high nonlinearity of basic equations governing the flame propagation, this assumption is far from being self-evident. For instance, averaging of a cubic combination of $`A`$’s gives rise to a term $`A_1^2A_0`$ comparable with $`A_0^3,`$ since the small and large scale parts of the flow variables are generally of the same order of magnitude. Clearly, equations for $`A_0`$’s involving such terms would not be of great value, since the local flame dynamics, and therefore, the coefficients $`A_1^2,`$ are unknown. The main result of the present work is the proof that such terms actually do not arise in the leading order with respect to $`\lambda _\mathrm{c}/L_0.`$ One can say that the governing equations for the quantities $`A_0`$ and $`A_1`$ decouple from each other. The proof consists of two parts corresponding to decoupling of the flow equations in the bulk, and decoupling of the jump conditions at the flame front, given in Secs. II.2 and II.3, respectively. As an application of this result, the problem of nonlinear front stabilization in a gravitational field will be considered in Sec. III. The results obtained are discussed in Sec. IV.

## II The decoupling theorem

### II.1 The averaging procedure

Let us begin with the precise formulation of the averaging procedure. Denote $`L_0`$ the characteristic length of the problem in question. For instance, $`L_0`$ can be the tube width, in the case of a flame propagating in a tube, or be related to an external field acting on the system. In practice, this length largely exceeds the flame cell size,

$$L_0\lambda _\mathrm{c}.$$

Assuming this, let us choose a length $`L`$ satisfying

$`\lambda _\mathrm{c}LL_0.`$
Analogously, denoting the characteristic time interval by $`T_0,`$ and noting that

$$T_0\frac{L_0}{U_\mathrm{f}},$$

we can choose $`TL/U_\mathrm{f}`$ such that

$$\frac{\lambda _\mathrm{c}}{U_\mathrm{f}}TT_0.$$

Given a function $`A(𝒙,t),`$ we define its space-time average over $`\{𝒙,t:𝒙(𝒙_0,𝒙_0+\mathrm{\Delta }𝒙),`$ $`\mathrm{\Delta }x_i=L,`$ $`t(t_0,t_0+T)\}`$

$`A={\displaystyle \frac{1}{L^3T}}{\displaystyle \underset{T}{}}𝑑t{\displaystyle \underset{V}{}}d^3𝒙A(𝒙,t)A_0(𝒙_0,t_0).`$ (4)
By the definition, $`A`$ varies noticeably over space distances $`\mathrm{\Delta }𝒙_0L_0,`$ and time intervals $`\mathrm{\Delta }t_0T_0.`$ The function $`A`$ thus turns out to be decomposed into two parts corresponding to the two scales, $`L_0`$ and $`\lambda _\mathrm{c}:`$

$$A=A_0+A_1,A_1=0.$$

As was mentioned in the Introduction, flame dynamics can be analyzed in the framework of the power expansion with respect to the small ratio $`L_\mathrm{f}/L_0\epsilon `$ (or equivalently, with respect to $`\lambda _\mathrm{c}/L_0,`$ since $`\lambda _\mathrm{c}=O(L_\mathrm{f})`$). Thus, we write $`A_0`$ and $`A_1`$ as follows

$$A_0=A_0^{(0)}+\epsilon A_0^{(1)}+\mathrm{},A_1=A_1^{(0)}+\epsilon A_1^{(1)}+\mathrm{},$$

dots denoting terms of higher order in $`\epsilon .`$ In this notation, the large scale flame dynamics in zero order approximation with respect to $`\epsilon `$ are described by the quantities $`A_0^{(0)}.`$ Our main purpose below will be to investigate coupling between $`A_0^{(0)}`$ and $`A_1^{(0)},`$ $`A_1^{(1)},`$ etc., and to obtain effective equations governing dynamics of $`A_0^{(0)}.`$ Accordingly, all quantities will be measured in units relevant to the large scale dynamics. Namely, space coordinates $`𝒙`$ and time $`t`$ are assumed to be normalized on $`L_0`$ and $`L_0/U_\mathrm{f},`$ respectively. Furthermore, $`U_\mathrm{f}`$ will be taken as the unit of gas velocity $`𝒗,`$ while $`\rho _\mathrm{u}U_\mathrm{f}^2`$ as the unit of gas pressure $`p.`$ For future reference, let us write down their expansions explicitly

$`𝒗=𝒗_0+𝒗_1,p=p_0+p_1,`$ (5)
$`𝒗_0=𝒗_0^{(0)}+\epsilon 𝒗_0^{(1)}+\mathrm{},𝒗_1=𝒗_1^{(0)}+\epsilon 𝒗_1^{(1)}+\mathrm{},`$ (6)
$`p_0=p_0^{(0)}+\epsilon p_0^{(1)}+\mathrm{},p_1=p_1^{(0)}+\epsilon p_1^{(1)}+\mathrm{}.`$ (7)
Within our choice of units, we have the following estimates

$`𝒗_0`$ $`=`$ $`O(1),𝒗_1=O(1),p_0=O(1),p_1=O(1),`$ (8)
$`{\displaystyle \frac{v_{0i}}{x_k}}`$ $`=`$ $`O(1),{\displaystyle \frac{v_{0i}}{t}}=O(1),{\displaystyle \frac{p_0}{x_k}}=O(1),`$ (9)
$`{\displaystyle \frac{v_{1i}}{x_k}}`$ $`=`$ $`O\left({\displaystyle \frac{1}{\epsilon }}\right),{\displaystyle \frac{v_{1i}}{t}}=O\left({\displaystyle \frac{1}{\epsilon }}\right),{\displaystyle \frac{p_1}{x_k}}=O\left({\displaystyle \frac{1}{\epsilon }}\right),`$ (10)
$`{\displaystyle \frac{^2v_{0i}}{x_kx_l}}`$ $`=`$ $`O(1),{\displaystyle \frac{^2v_{1i}}{x_kx_l}}=O\left({\displaystyle \frac{1}{\epsilon ^2}}\right),i,k,l=1,2,3.`$ (11)
Let us now proceed to the examination of flame dynamics in terms of $`𝒗_{0,1}^{(0,1)},p_{0,1}^{(0,1)}.`$

The main result concerning the large scale flame dynamics which will be proved below can be expressed in the form of the following

Decoupling theorem: The large scale dynamics of a flame are unaffected by its local cellular structure up to a rescaling. More precisely, the form of dynamical equations for $`𝒗_0^{(0)},`$ $`p_0^{(0)},`$ as well as of jump conditions for these quantities at the flame front, is the same as in the case of zero-thickness flame propagating in an ideal fluid at constant speed $`𝔘U_\mathrm{f}`$ with respect to the fuel, the number $`𝔘>1`$ describing the flame front length increase due to the local front wrinkling.

The proof consists of two parts corresponding to decoupling of the flow equations in the bulk, and decoupling of the jump conditions at the flame front, presented in Secs. II.2 and II.3, respectively.

### II.2 Decoupling of dynamical equations

For definiteness, we will assume in what follows that external field acting on the system is the gravitational field, denoting its strength by $`𝒈.`$ Accordingly, $`L_0`$ will be identified with the characteristic length associated with this field:111If the gravitational field is not homogeneous, it is assumed to vary noticeably over distances larger than $`L_0.`$

$`L_0={\displaystyle \frac{U_\mathrm{f}^2}{𝒈}}.`$ (12)
Then the dimensionless velocity and pressure fields obey the following equations in the bulk

$`\mathrm{div}𝒗`$ $`=`$ $`0,`$ (13)
$`{\displaystyle \frac{𝒗}{t}}+(𝒗\mathbf{})𝒗`$ $`=`$ $`{\displaystyle \frac{1}{\varrho }}\mathbf{}p+𝑮+\epsilon Pr\mathrm{}𝒗,`$ (14)
where

$$𝑮=\frac{𝒈L_0}{U_\mathrm{f}^2},𝑮=1,$$

$`\varrho `$ is the fluid density normalized on the fuel density $`\rho _\mathrm{u},`$ and $`Pr`$ the Prandtl number representing the ratio of viscous and thermal diffusivities, $`Pr=\nu /\chi .`$

Substituting expansions (5)–(7) into Eqs. (13), (14), taking into account the estimates (8)–(11), and extracting $`O(1/\epsilon )`$ terms yields

$`\mathrm{div}𝒗_1^{(0)}`$ $`=`$ $`0,`$ (15)
$`{\displaystyle \frac{𝒗_1^{(0)}}{t}}+\left(\left[𝒗_0^{(0)}+𝒗_1^{(0)}\right]\mathbf{}\right)𝒗_1^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{\varrho }}\mathbf{}p_1^{(0)}+\epsilon Pr\mathrm{}𝒗_1^{(0)}.`$ (16)
Next, collecting $`O(1)`$ terms gives

$`\mathrm{div}𝒗_0^{(0)}`$ $`=`$ $`0,`$ (17)
$`{\displaystyle \frac{𝒗_0^{(0)}}{t}}+\epsilon {\displaystyle \frac{𝒗_1^{(1)}}{t}}+\left(\left[𝒗_0^{(0)}+𝒗_1^{(0)}\right]\mathbf{}\right)𝒗_0^{(0)}`$ $`+`$ $`\epsilon \left(\left[𝒗_0^{(0)}+𝒗_1^{(0)}\right]\mathbf{}\right)𝒗_1^{(1)}`$
$`={\displaystyle \frac{1}{\varrho }}\mathbf{}p_0^{(0)}`$ $``$ $`{\displaystyle \frac{\epsilon }{\varrho }}\mathbf{}p_1^{(1)}+𝑮+\epsilon ^2Pr\mathrm{}𝒗_1^{(1)}.`$ (18)
Equation (18) involves both slowly and rapidly varying terms. The slowly varying part of this equation, determining dynamics of the fields $`𝒗_0^{(0)},`$ $`p_0^{(0)},`$ can be separated out by averaging it according to Eq. (4). Under this operation, all terms linear in $`𝒗_1,p_1`$ give rise to $`o(1)`$ contribution. For instance,

$`\epsilon {\displaystyle \frac{𝒗_1^{(1)}}{t}}={\displaystyle \frac{\epsilon }{L^3T}}{\displaystyle \underset{V}{}}d^3𝒙𝒗_1^{(1)}|_{t_0}^{t_0+T}={\displaystyle \frac{L_\mathrm{f}}{L_0}}O\left({\displaystyle \frac{T_0}{T}}\right)=O\left({\displaystyle \frac{L_\mathrm{f}}{L}}\right)=o(1),`$ (19)
in view of the estimates (8), and the choice of $`L,T.`$ The same argument applies to $`(𝒗_0^{(0)}\mathbf{})𝒗_1^{(1)},`$ as well as to the second and fourth terms in the right hand side of Eq. (18). Furthermore,

$$\left(𝒗_1^{(0)}\mathbf{}\right)𝒗_0^{(0)}0$$

according to the definition of $`𝒗_1.`$ Finally, contribution of the last term in the left hand side of (18) also is $`o(1).`$ Indeed, integrating by parts and taking into account Eq. (15), we have

$`\left(𝒗_1^{(0)}\mathbf{}\right)𝒗_1^{(1)}={\displaystyle \frac{1}{L^3T}}{\displaystyle \underset{T}{}}𝑑t{\displaystyle \underset{V}{}}d^3𝒙\left(𝒗_1^{(0)}\mathbf{}\right)𝒗_1^{(1)}={\displaystyle \frac{1}{L^3T}}{\displaystyle \underset{T}{}}𝑑t{\displaystyle \underset{S}{}}\left(d𝒔𝒗_1^{(0)}\right)𝒗_1^{(1)},`$ (20)
where $`S`$ is the surface of the cube $`V=\{𝒙:𝒙(𝒙_0,𝒙_0+\mathrm{\Delta }𝒙),`$ $`\mathrm{\Delta }x_i=L\},`$ $`d𝒔`$ being its element. Using Eqs. (8), the right hand side of Eq. (20) is estimated as $`O(L_0/L).`$ Hence,

$$\epsilon \left(𝒗_1^{(0)}\mathbf{}\right)𝒗_1^{(1)}=\frac{L_\mathrm{f}}{L_0}O\left(\frac{L_0}{L}\right)=o(1).$$

Thus, Eq. (18) reduces upon averaging to the ordinary Euler equation for the functions $`𝒗_0^{(0)},`$ $`p_0^{(0)}`$

$`{\displaystyle \frac{𝒗_0^{(0)}}{t}}+\left(𝒗_0^{(0)}\mathbf{}\right)𝒗_0^{(0)}={\displaystyle \frac{1}{\varrho }}\mathbf{}p_0^{(0)}+𝑮+o(1),`$ (21)
which proves the first part of the decoupling theorem. It is worth of mentioning that the large scale flow dynamics in the bulk turn out to be ideal at zeroth order in $`\epsilon .`$

In connection with Eqs. (15), (16) the following circumstance should be emphasized. These equations describe bulk dynamics of the small scale parts of the flow variables at zeroth order in $`L_\mathrm{f}/L_0,`$ i.e., when the influence of the large scale flows on the flame cellular structure is completely neglected. This might seem to be in contradiction with the structure of Eq. (16), because it involves $`𝒗_0^{(0)}`$ explicitly. However, the functions $`𝒗_1^{(0)}(𝒙,t),`$ $`p_1^{(0)}(𝒙,t),`$ satisfying Eq. (16) in a given space-time region $`\{𝒙,t:𝒙(𝒙_0,𝒙_0+\mathrm{\Delta }𝒙),`$ $`\mathrm{\Delta }x_i=L,`$ $`t(t_0,t_0+T)\},`$ can be written as

$$𝒗_1^{(0)}(𝒙,t)=\stackrel{~}{𝒗}_1^{(0)}(𝒙𝒗_0^{(0)}(𝒙_0,t_0)t,t),p_1^{(0)}(𝒙,t)=\stackrel{~}{p}_1^{(0)}(𝒙𝒗_0^{(0)}(𝒙_0,t_0)t,t),$$

where $`\stackrel{~}{𝒗}_1^{(0)},`$ $`\stackrel{~}{p}_1^{(0)}`$ satisfy

$`{\displaystyle \frac{\stackrel{~}{𝒗}_1^{(0)}}{t}}+\left(\stackrel{~}{𝒗}_1^{(0)}\mathbf{}\right)\stackrel{~}{𝒗}_1^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{\varrho }}\mathbf{}\stackrel{~}{p}_1^{(0)}+\epsilon Pr\mathrm{}\stackrel{~}{𝒗}_1^{(0)}.`$
In other words, the role of $`𝒗_0^{(0)}`$ in Eq. (16) is purely kinematical: it describes the large scale “drift” of the flame cellular structure.

### II.3 Decoupling of jump conditions

The proof of decoupling of the jump conditions is more complicated, since this is the place where the transport processes inside the flame front come into play. These conditions express the conservation of energy and momentum across the flame front. For freely propagating flames, and within the accuracy of the first order in the small front thickness, they were derived in the most general form in Ref. matalon . To take into account the influence of gravity, it is sufficient to note that the bulk equations (14) can be rendered formally free by substituting $`p=\stackrel{~}{p}+\varrho (𝑮𝒙).`$ However, gravity reappears through the jump conditions at the flame front. On the other hand, the influence of gravity on gas dynamics inside the flame front is small in comparison with the transport effects; their relative value is known to be given by the inverse Froude number $`Fr^1=L_\mathrm{f}𝒈/U_\mathrm{f}^2=\epsilon .`$ To the leading order in $`\epsilon ,`$ therefore, contribution of the gravitational field to the jump conditions is the same as in the case of a zero-thickness flame.

For simplicity, we will consider two-dimensional (2D) case, assuming also that the Lewis number is equal to unity. No assumption is made concerning the incoming flow, except that its characteristic length $`\stackrel{~}{L}L_0.`$ Let the flame front position be described by an equation $`z=f(x,t),`$ where the Cartesian coordinates $`(x,z)`$ are scaled on $`L_0,`$ and chosen so that $`z`$-axis is parallel to $`𝑮.`$ The $`x`$\- and $`z`$-components of the flow velocity will be denoted by $`w`$ and $`u,`$ $`𝒗=(w,u).`$ We also introduce the unit vector $`𝝉`$ tangential to the flame front, and $`𝒏`$ orthogonal to it (pointing to the burnt matter). In components,

$`𝝉=({\displaystyle \frac{1}{N}},{\displaystyle \frac{f/x}{N}}),𝒏=({\displaystyle \frac{f/x}{N}},{\displaystyle \frac{1}{N}}),N\sqrt{1+\left({\displaystyle \frac{f}{x}}\right)^2}.`$ (22)
Rewriting Eqs. (5.32)–(5.43) of matalon  in this notation for the 2D case, and taking into account the contribution of the archimedean force to the pressure jump yields

$`(𝒗_+𝒏)(𝒗_{}𝒏)`$ $`=`$ $`(\theta 1),`$ (23)
$`(𝒗_+𝝉)(𝒗_{}𝝉)`$ $`=`$ $`\epsilon \left(\mathrm{ln}\theta +(\theta 1)Pr\right){\displaystyle \frac{1}{N}}\left(\widehat{D}w_{}+{\displaystyle \frac{f}{x}}\widehat{D}u_{}+{\displaystyle \frac{1}{N}}\widehat{D}{\displaystyle \frac{f}{x}}\right),`$ (24)
$`p_+p_{}`$ $`=`$ $`\stackrel{~}{p}_+\stackrel{~}{p}_{}+\left(\varrho G_zz\right)_+\left(\varrho G_zz\right)_{}`$ (25)
$`=`$ $`(\theta 1){\displaystyle \frac{\theta 1}{\theta }}Gf+\epsilon (\theta 1){\displaystyle \frac{}{x}}\left({\displaystyle \frac{1}{N}}{\displaystyle \frac{f}{x}}\right)`$
$`+`$ $`{\displaystyle \frac{\epsilon \mathrm{ln}\theta }{N}}\left({\displaystyle \frac{^2f}{t^2}}+2w_{}{\displaystyle \frac{^2f}{tx}}+w_{}^2{\displaystyle \frac{^2f}{x^2}}+2\widehat{D}N{\displaystyle \frac{1}{N}}{\displaystyle \frac{f}{x}}{\displaystyle \frac{N}{x}}\right),`$
where

$$\widehat{D}\frac{}{t}+\left(w_{}+\frac{1}{N}\frac{f}{x}\right)\frac{}{x},GG_z,$$

and the subscripts “$``$” and “$`+`$” mean that the corresponding quantity is calculated for $`z=f(x,t)0`$ and $`z=f(x,t)+0,`$ respectively.

Finally, to complete the system of hydrodynamic equations and jump conditions, one needs an expression for the local burning rate. This expression, the so-called evolution equation, has the following form (Cf. Eq. (6.1) in Ref. matalon )

$`(𝒗_{}𝒏){\displaystyle \frac{1}{N}}{\displaystyle \frac{f}{t}}=1{\displaystyle \frac{\epsilon \theta \mathrm{ln}\theta }{N(\theta 1)}}\left({\displaystyle \frac{N}{t}}+{\displaystyle \frac{}{x}}(Nw_{})+{\displaystyle \frac{^2f}{x^2}}\right).`$ (26)
At this point, it is worth to make the following comment on the meaning of the above asymptotic relations. Equations (23)–(26) were derived in Ref. matalon  under assumption that the terms in the right hand sides of these equations, proportional to $`\epsilon ,`$ are small, which is only true if the gas flow is characterized by a length scale much larger than the flame front thickness. However, the rapidly developing LD-instability makes any smooth flame configuration highly corrugated within the time interval of the order (3). As a result, the $`\epsilon `$-terms turn out to be of the order $`L_\mathrm{f}/\lambda _\mathrm{c}=O(1).`$ Similarly, account of the $`\epsilon ^2`$-corrections in the above equations would give rise to terms of the order $`L_\mathrm{f}^2/\lambda _\mathrm{c}^2=O(1)`$ etc., questioning thereby validity of the small $`\epsilon `$-expansion. However, it was mentioned in the Introduction that in practice, the cut-off wavelength $`\lambda _\mathrm{c}`$ is noticeably larger than the flame front thickness $`L_\mathrm{f}.`$ Thus, in the regime of fully developed LD-instability, the right hand sides of Eqs. (23)–(26) are to be considered the leading order terms of the asymptotic expansion in powers of $`L_\mathrm{f}/\lambda _\mathrm{c}=1/c_1,`$ rather than $`L_\mathrm{f}/L_0=\epsilon .`$ On the contrary, $`\epsilon `$ is the true parameter of the power expansions (6), (7), which determines the relative order of successive terms in these expansions.

In order to extract from Eqs. (23)–(25) jump conditions for the quantities $`𝒗_0^{(0)},`$ $`p_0^{(0)},`$ we need to introduce an auxiliary operation of averaging along the front. Let a quantity $`A`$ be defined on the flame front, i.e., for $`\{x,z,t:z=f(x,t)\}.`$ Given a point $`x_0,`$ choose $`\mathrm{\Delta }x=\mathrm{\Delta }x(x_0)`$ such that the front length $`𝔏(t)`$ between the points $`(x_0,f(x_0,t))`$ and $`(x_0,f(x_0+\mathrm{\Delta }x,t))`$ satisfies

$$\lambda _\mathrm{c}𝔏(t)L_0,𝔏(t)=O(L)$$

for all $`t(t_0,t_0+T).`$ This is always possible since $`𝔏(t)`$ is of the order of distance between the two points. Then the average value of $`A`$ over $`\{x,z,t:x(x_0,x_0+\mathrm{\Delta }x),t(t_0,t_0+T),z=f(x,t)\}`$ is defined as

$`A_l={\displaystyle \frac{1}{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑lA,𝔚={\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t𝔏(t)={\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑l,`$ (27)
$`dl`$ being the front line element, $`dl=Ndx.`$

Using the operation introduced, the quantities $`𝝉,𝒏`$ defined on the front, as well as the flame front position itself, can be decomposed into two parts corresponding to the scales $`L_0`$ and $`\lambda _\mathrm{c},`$ in a way analogous to Eqs. (5)–(7):

$`𝝉=𝝉_0+𝝉_1,(n_x,n_z)=(\tau _z,\tau _x),f=f_0+f_1,𝝉_1_l=0,f_1_l=0,`$ (28)
$`𝝉_{0,1}=𝝉_{0,1}^{(0)}+O(\epsilon ),f_0=f_0^{(0)}+O(\epsilon )f_1=O(\epsilon ).`$ (29)
To obtain jump conditions for the quantities $`𝒗_0^{(0)},`$ $`p_0^{(0)},`$ expansions (5)–(7) and (28) should be inserted into Eqs. (23)–(25), with the subsequent averaging of the latter along the flame front.

Let us begin with the jump of the normal component of the gas velocity. In view of Eq. (13), one can introduce the stream function $`\psi =\psi (x,z,t)`$ according to

$`u={\displaystyle \frac{\psi }{x}},w={\displaystyle \frac{\psi }{z}}.`$ (30)
Using the operation of the bulk averaging (4), the function $`\psi `$ can be decomposed as

$`\psi =\psi _0+\psi _1,\psi _1=0.`$ (31)
It follows from Eqs. (8) that

$`\psi _0=O(1),\psi _1=O(\epsilon ).`$ (32)
Substituting Eq. (31) into Eqs. (30), and averaging gives

$`u_0={\displaystyle \frac{(\psi _0+\psi _1)}{x}}={\displaystyle \frac{\psi _0}{x}}+{\displaystyle \frac{1}{L^2T}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{z_0}{\overset{z_0+L}{}}}𝑑z\psi _1|_{x_0}^{x_0+L}={\displaystyle \frac{\psi _0}{x}}+O\left({\displaystyle \frac{L_\mathrm{f}}{L}}\right),`$
and analogous equation for $`w_0.`$ Thus, up to $`o(1)`$ terms, one has

$`u_0={\displaystyle \frac{\psi _0}{x}},w_0={\displaystyle \frac{\psi _0}{z}}.`$ (33)
Inserting Eqs. (30) into Eq. (23), and using Eq. (22), one finds

$`(𝒗_\pm 𝒏)_l`$ $`=`$ $`{\displaystyle \frac{1}{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑x\left({\displaystyle \frac{\psi }{x}}+{\displaystyle \frac{\psi }{z}}{\displaystyle \frac{f}{x}}\right)_\pm ={\displaystyle \frac{1}{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑x{\displaystyle \frac{d\psi _\pm }{dx}}`$ (34)
$`=`$ $`{\displaystyle \frac{1}{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t\psi _{0\pm }|_{x_0}^{x_0+\mathrm{\Delta }x}+o(1).`$
Since the variation of $`\psi _0`$ over space distances $`L`$ and time intervals $`T`$ is small, taking into account Eqs. (29), (33), and neglecting $`O(\epsilon )`$ terms, one can write

$`{\displaystyle \frac{1}{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t\psi _{0\pm }|_{x_0}^{x_0+\mathrm{\Delta }x}`$ $`=`$ $`{\displaystyle \frac{T}{𝔚}}\left[\psi _{0\pm }(x_0+\mathrm{\Delta }x,f_0(x_0+\mathrm{\Delta }x,t_0),t_0)\psi _{0\pm }(x_0,f_0(x_0,t_0),t_0)\right]`$ (35)
$`=`$ $`{\displaystyle \frac{T}{𝔚}}\left({\displaystyle \frac{\psi _0}{x}}+{\displaystyle \frac{\psi _0}{z}}{\displaystyle \frac{f_0}{x}}\right)_\pm \mathrm{\Delta }x={\displaystyle \frac{T}{𝔚}}\left(u_{0\pm }w_{0\pm }{\displaystyle \frac{f_0}{x}}\right)\mathrm{\Delta }x.`$
Substituting this into Eq. (23) gives

$`\left[u_{0+}u_0(w_{0+}w_0){\displaystyle \frac{f_0}{x}}\right]{\displaystyle \frac{T\mathrm{\Delta }x}{𝔚}}=(\theta 1)+o(1).`$ (36)
Note that

$`n_{x0}={\displaystyle \frac{f/x}{N}}_l={\displaystyle \frac{1}{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑x\left({\displaystyle \frac{f_0}{x}}+{\displaystyle \frac{f_1}{x}}\right)={\displaystyle \frac{T\mathrm{\Delta }x}{𝔚}}{\displaystyle \frac{f_0}{x}}+o(1),`$ (37)
and similarly,

$`n_{z0}={\displaystyle \frac{T\mathrm{\Delta }x}{𝔚}}+o(1).`$ (38)
Hence, Eq. (36) can be rewritten as

$`\left(𝒗_{0+}𝒏_0\right)\left(𝒗_0𝒏_0\right)=(\theta 1)+o(1),`$
or, with the same accuracy,

$`\left(𝒗_{0+}^{(0)}𝒏_0^{(0)}\right)\left(𝒗_0^{(0)}𝒏_0^{(0)}\right)=(\theta 1)+o(1).`$ (39)
The inverse norm of $`𝒏_\mathrm{𝟎}`$

$`𝒏_0^1={\displaystyle \frac{{\displaystyle \frac{1}{T}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t𝔏(t)}{\mathrm{\Delta }x\sqrt{1+\left({\displaystyle \frac{f_0}{x}}\right)^2}}}𝔘`$ (40)
has a clear geometrical meaning. Namely, $`(𝔘1)`$ represents the relative increase of the flame front length due to its small scale wrinkling.

Let us turn to the examination of the remaining jump conditions. Because of the $`\epsilon `$-terms in the right hand sides of Eqs. (24), (25), which involve highly nonlinear combinations of the flow variables and the function $`f(x,t),`$ there seems to be a very little hope that the large scale parts of the flow variables eventually decouple from their small scale parts describing flame cellular structure. Nevertheless, they do, as will be shown presently.

According to Eqs. (8)–(11) and analogous estimates for the space-time derivatives of $`f(x,t),`$

$`{\displaystyle \frac{f}{x}}`$ $`=`$ $`O(1),{\displaystyle \frac{f}{t}}=O(1),`$ (41)
$`{\displaystyle \frac{^2f}{t^2}}`$ $`=`$ $`O\left({\displaystyle \frac{1}{\epsilon }}\right),{\displaystyle \frac{^2f}{xt}}=O\left({\displaystyle \frac{1}{\epsilon }}\right),{\displaystyle \frac{^2f}{x^2}}=O\left({\displaystyle \frac{1}{\epsilon }}\right),`$ (42)
the $`\epsilon `$-terms are $`O(1).`$ Let us show first that the average value of the right hand side of Eq. (24) along the flame front is actually $`o(1).`$ Using the evolution equation (26), expression in the parentheses in the right hand side of Eq. (24) can be rewritten as follows:222It was mentioned after Eq. (26), that the $`\epsilon `$-terms in the jump conditions represent the leading order terms of the asymptotic expansion in powers of $`1/c_1.`$ In transforming these terms, therefore, one can use the evolution equation with the $`\epsilon `$-term omitted.

$`\widehat{D}w_{}+f^{}\widehat{D}u_{}+{\displaystyle \frac{1}{N}}\widehat{D}f^{}`$
$`=\dot{w}_{}+f^{}\dot{u}_{}+{\displaystyle \frac{\dot{f}^{}}{N}}+w_{}w_{}^{}+{\displaystyle \frac{(f^{}w_{})^{}}{N}}+{\displaystyle \frac{\left(f^{}\right)^2u_{}^{}}{N}}+f^{}w_{}u_{}^{}+{\displaystyle \frac{N^{}}{N}}`$
$`=\dot{w}_{}+f^{}\dot{u}_{}+{\displaystyle \frac{\dot{f}^{}}{N}}+w_{}w_{}^{}+{\displaystyle \frac{\left(u_{}N\dot{f}\right)^{}}{N}}+{\displaystyle \frac{N^21}{N}}u_{}^{}`$
$`+u_{}^{}\left(u_{}N\dot{f}\right)+{\displaystyle \frac{N^{}}{N}}`$
$`=\dot{w}_{}+{\displaystyle \frac{(f^{}u_{})}{t}}(u_{}\dot{f})^{}+{\displaystyle \frac{\left(u_{}^2+w_{}^2\right)^{}}{2}},`$ (43)
where the dot and the prime denote differentiation with respect to $`t`$ and $`x,`$ respectively. Hence, taking into account the estimates (8), (41), one has

$$\frac{1}{N}\left(\widehat{D}w_{}+f^{}\widehat{D}u_{}+\frac{1}{N}\widehat{D}f^{}\right)_l=$$

$`={\displaystyle \frac{1}{𝔚}}\left({\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑x[w_{}+f^{}u_{}]_{t_0}^{t_0+T}+{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t\left[u_{}\dot{f}+{\displaystyle \frac{u_{}^2+w_{}^2}{2}}\right]_{x_0}^{x_0+\mathrm{\Delta }x}\right)`$
$`=O\left({\displaystyle \frac{T_0}{T}}\right)+O\left({\displaystyle \frac{L_0}{𝔏}}\right).`$ (44)
In view of the choice of $`T,𝔏,`$ the right hand side of Eq. (24) turns out to be $`O(\epsilon L_0/𝔏)=O(L_\mathrm{f}/𝔏)=o(1).`$ Thus, averaging of Eq. (24) gives

$$(𝒗_+𝝉)(𝒗_{}𝝉)_l=o(1).$$

To the leading order in $`\epsilon ,`$ this equation can be written as

$`\left(\left[𝒗_{0+}^{(0)}𝒗_0^{(0)}\right]𝝉^{(0)}\right)_l+\left(\left[𝒗_{1+}^{(0)}𝒗_1^{(0)}\right]𝝉^{(0)}\right)_l=o(1),𝝉^{(0)}=𝝉_0^{(0)}+𝝉_1^{(0)}.`$ (45)
It was mentioned in the end of Sec. II.2 that in the approximation considered, the large scale flows do not affect the flame cellular structure. In particular, the two directions along the flame front, $`𝝉`$ and $`𝝉,`$ are left equivalent. This implies that up to $`O(\lambda _\mathrm{c}/𝔏)`$ terms, the value of $`(𝒗_{1+}^{(0)}𝒗_1^{(0)})𝝉^{(0)}_l`$ must be invariant under the change $`𝝉^{(0)}𝝉^{(0)}.`$ Hence,

$$\left(\left[𝒗_{1+}^{(0)}𝒗_1^{(0)}\right]𝝉^{(0)}\right)_l=\left(\left[𝒗_{1+}^{(0)}𝒗_1^{(0)}\right]𝝉^{(0)}\right)_l=0.$$

Thus, taking into account definition of $`𝝉_1,`$ we have from Eq. (45)

$`\left(\left[𝒗_{0+}^{(0)}𝒗_0^{(0)}\right]𝝉_0^{(0)}\right)=o(1).`$ (46)
Consider next the pressure jump, Eq. (25). The last term in the right hand side of this equation can be transformed as333See the footnote 2.

$`\ddot{f}+2w_{}\dot{f}^{}+w_{}^2f^{\prime \prime }+2\widehat{D}N{\displaystyle \frac{f^{}N^{}}{N}}`$
$`=\ddot{f}+2w_{}\dot{f}^{}+\left(w_{}^2f^{}\right)^{}2w_{}^{}(u_{}N\dot{f})+2\left(\dot{N}+w_{}N^{}+{\displaystyle \frac{f^{}N^{}}{N}}\right){\displaystyle \frac{f^{}N^{}}{N}}`$
$`=\ddot{f}+2\dot{N}+2(w_{}\dot{f})^{}+2(w_{}N)^{}+\left(w_{}^2f^{}\right)^{}2w_{}^{}u_{}+{\displaystyle \frac{f^{}N^{}}{N}}`$
$`=\ddot{f}+2\dot{N}(w_{}^2f^{})^{}+2w_{}u_{}^{}+{\displaystyle \frac{f^{}N^{}}{N}}.`$ (47)
As before, the first three terms give rise to $`o(1)`$ terms upon averaging. The remaining two terms, however, do not reduce to the full $`x`$\- or $`t`$-derivatives. Their contribution is, therefore, $`O(1),`$ in general. Notice that the fourth term is quadratic in the gas velocity. This fact can be used to show that its contribution is independent of the functions $`𝒗_0,p_0.`$ Indeed, one has, to the leading order in $`\epsilon ,`$

$`\epsilon {\displaystyle \frac{w_{}u_{}^{}}{N}}_l`$ $`=`$ $`{\displaystyle \frac{\epsilon }{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑x\left(w_0^{(0)}+w_1^{(0)}\right)\left(u_0^{(0)}+u_1^{(0)}\right)^{}`$
$`=`$ $`{\displaystyle \frac{\epsilon }{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑x\left(w_0^{(0)}+w_1^{(0)}\right)u_1^{(0)}+O(\epsilon )`$
$`=`$ $`{\displaystyle \frac{\epsilon }{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t\left[w_0^{(0)}u_1^{(0)}\right]_{x_0}^{x_0+\mathrm{\Delta }x}+{\displaystyle \frac{\epsilon }{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑xw_1^{(0)}u_1^{(0)}+O(\epsilon )`$
$`=`$ $`{\displaystyle \frac{\epsilon }{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑xw_1^{(0)}u_1^{(0)}+O\left({\displaystyle \frac{L_\mathrm{f}}{𝔏}}\right).`$
Thus, up to $`o(1)`$ terms, the average value of $`\epsilon w_{}u_{}^{}/N`$ turns out to be independent of the functions $`𝒗_0^{(0)}.`$ The quantities $`w_1^{(0)},`$ $`u_1^{(0)}`$ describe variations of the fuel velocity along the front cell in zero order approximation with respect to $`L_\mathrm{f}/L_0,`$ i.e., when the influence of the large scale flow on the local flame structure is completely neglected. In particular, the value of $`\epsilon w_{}u_{}^{}/N_l`$ is independent of the coordinate $`x_0`$ as well as of the time instant $`t_0.`$ Denote this constant by $`\alpha _1/2.`$444 At first sight, $`\epsilon w_{}u_{}^{}/N_l`$ depends on the choice of orientation of the coordinate axes, while the scalar pressure jump must be independent of this choice. It is easy to see, however, that within the accuracy of the above calculations, $`\epsilon w_{}u_{}^{}/N_l`$ is actually coordinate-invariant. In fact, under rotations of the coordinate system, $`w_{},`$ $`u_{}`$ transform as    $`w_{}`$ $``$ $`\mathrm{cos}\phi w_{}+\mathrm{sin}\phi u_{},`$    $`u_{}`$ $``$ $`\mathrm{sin}\phi w_{}+\mathrm{cos}\phi u_{},`$   where $`\phi `$ is the rotation angle. This transformation leaves $`\epsilon w_{}u_{}^{}/N_l`$ unchanged:    $`\epsilon w_{}u_{}^{}/N_l`$ $``$ $`{\displaystyle \frac{\epsilon }{𝔚}}{\displaystyle 𝑑t(\mathrm{cos}\phi w_{}+\mathrm{sin}\phi u_{})d(\mathrm{sin}\phi w_{}+\mathrm{cos}\phi u_{})}`$     $`=`$ $`\epsilon w_{}u_{}^{}/N_l+O\left({\displaystyle \frac{L_\mathrm{f}}{𝔏}}\right).`$  Definition of $`\alpha _1`$ can be written also in an explicitly invariant form:
$$\alpha _1=\epsilon \frac{w_{}u_{}^{}u_{}w_{}^{}}{N}_l.$$

It remains to find the contribution of the last term in Eq. (II.3). We have

$`\epsilon {\displaystyle \frac{f^{}N^{}}{N^2}}_l={\displaystyle \frac{\epsilon }{𝔚}}{\displaystyle \underset{t_0}{\overset{t_0+T}{}}}𝑑t{\displaystyle \underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}}𝑑x{\displaystyle \frac{f^{}N^{}}{N}}=\epsilon \left({\displaystyle \frac{f^{}}{N}}\right)^{}_l+O\left({\displaystyle \frac{L_\mathrm{f}}{𝔏}}\right).`$ (48)
The quantity $`\left(f^{}/N\right)^{}`$ is nothing but the flame front curvature $`k,`$

$$k=\frac{f^{\prime \prime }}{N^3}.$$

Using the definitions (22), (28), one can write

$$k_l=\left(\frac{f^{}}{N}\right)^{}_l=\tau _z^{}_l=\tau _{z0}^{}_l+\tau _{z1}^{}_l=\tau _{z0}^{}+\tau _{z1}^{}_l.$$

According to Eqs. (41), (42), $`\tau _{z0}^{}=O(1),`$ $`\tau _{z1}^{}=O(1/\epsilon ).`$ Thus, to the leading order in $`\epsilon ,`$

$$\epsilon k_l=\epsilon \tau _{z1}^{(0)}_l=O(1).$$

Similarly to $`𝒗_1^{(0)},`$ $`\tau _{z1}^{(0)}`$ describes geometry of a front cell neglecting the influence of the large scale flame structure on it. Hence, the value of $`\tau _{z1}^{(0)}_l`$ is independent of the particular choice of the point $`x_0`$ on the flame front (and of the time instant $`t_0`$). Denoting this constant by $`k_1,`$ we thus obtain from Eq. (25) the following expression for the jump of $`p_0^{(0)}`$ at the flame front

$`p_{0+}^{(0)}p_0^{(0)}={\displaystyle \frac{\theta 1}{\theta }}Gf_0^{(0)}(\theta 1)\pi _1+\epsilon [(\theta 1)\mathrm{ln}\theta ]k_1+\epsilon \mathrm{ln}\theta \alpha _1+o(1),`$ (49)
where $`\pi _1`$ is another constant defined by

$$\pi _1=p_{1+}^{(0)}p_1^{(0)}_l.$$

It is independent of $`x_0,t_0`$ on the same grounds as $`\alpha _1,k_1.`$

Finally, we have to average the evolution equation (26). As before, contribution of the $`\epsilon `$ term on the right hand side of this equation is $`o(1).`$ Furthermore,

$$\frac{1}{N}\frac{f}{t}_l=\frac{1}{𝔚}\underset{x_0}{\overset{x_0+\mathrm{\Delta }x}{}}𝑑xf_0|_{t_0}^{t_0+T}+o(1)=\frac{T\mathrm{\Delta }x}{𝔚}\frac{f_0}{t}+o(1).$$

Using Eqs. (34)–(38), we find

$`(𝒗_0𝒏_0)n_{z0}\dot{f}_0=1+o(1).`$ (50)
Equations (39), (46), (49), and (50) constitute the proof of the second part of the decoupling theorem. To make this more transparent, let us rewrite these equations in a more convenient notation. The jump conditions for the large scale parts of the flow variables read555The vector quantities denoted by Gothic letters are designated with arrows.

$`\left(\stackrel{}{𝔳}_+\stackrel{}{𝔫}\right)\left(\stackrel{}{𝔳}_{}\stackrel{}{𝔫}\right)=(\theta 1)𝔘+o(1),`$ (51)
$`\left(\stackrel{}{𝔳}_+\stackrel{}{𝔱}\right)\left(\stackrel{}{𝔳}_{}\stackrel{}{𝔱}\right)=o(1),`$ (52)
$`𝔭_+𝔭_{}={\displaystyle \frac{\theta 1}{\theta }}G𝔣+\mathrm{\Pi }+o(1),`$ (53)
where

$`\stackrel{}{𝔫}`$ $`=`$ $`{\displaystyle \frac{𝒏_0^{(0)}}{𝒏_0^{(0)}}}=({\displaystyle \frac{𝔣^{}}{𝒩}},{\displaystyle \frac{1}{𝒩}}),𝔣=f_0^{(0)},𝒩=\sqrt{1+\left(𝔣^{}\right)^2},\stackrel{}{𝔳}=𝒗_0^{(0)},`$
$`\stackrel{}{𝔱}`$ $`=`$ $`{\displaystyle \frac{𝝉_0^{(0)}}{𝝉_0^{(0)}}}=({\displaystyle \frac{1}{𝒩}},{\displaystyle \frac{𝔣^{}}{𝒩}}),𝒏_0^{(0)}=𝝉_0^{(0)}=𝔘^1,`$
$`𝔭`$ $`=`$ $`p_0^{(0)},\mathrm{\Pi }=(\theta 1)\pi _1+\epsilon [(\theta 1)\mathrm{ln}\theta ]k_1+\epsilon \mathrm{ln}\theta \alpha _1.`$
The evolution equation takes the form

$`\left(\stackrel{}{𝔳}_{}\stackrel{}{𝔫}\right){\displaystyle \frac{\dot{𝔣}}{𝒩}}=𝔘+o(1).`$ (54)
As we have seen, in zero order approximation with respect to $`\epsilon ,`$ the quantities $`k_1,`$ $`\alpha _1,`$ $`\pi _1,`$ and $`𝔘`$ entering these equations are independent of the coordinate $`x_0`$ and the time instant $`t_0.`$ These constants can in principle be calculated provided that the exact small scale flame structure is known. However, as it follows from the above equations, this information is actually unnecessary, because dynamics of the large scale fields are independent of the specific values of these constants. Indeed, since the gas pressure enters dynamical equation only through its gradient, it is determined up to a constant. Therefore, the constants $`k_1,`$ $`\alpha _1,`$ $`\pi _1,`$ are irrelevant. Unlike these, however, the constant $`𝔘`$ has a direct physical meaning. As is seen from Eq. (54), $`𝔘`$ plays the role of the effective dimensionless velocity of the curved flame propagation. According to our choice of units, the gas velocity is scaled on the adiabatic velocity of a plane flame front $`U_\mathrm{f}`$ which is also used to define the units of pressure, time, and length \[see the definition (12) of $`L_0`$\]. In analyzing the large scale flame dynamics, it is more natural to choose $`𝔘U_\mathrm{f},`$ rather than $`U_\mathrm{f},`$ as the velocity unit. Then $`𝔘`$ disappears from the jump conditions at the flame front, as well as from the flow equations (17), (21) in the bulk, which thus take the form of the equations governing propagation of a zero-thickness flame in an ideal fluid at constant speed $`𝔘U_\mathrm{f}`$ with respect to the fuel, $`𝔣`$ being the flame front position, while $`\stackrel{}{𝔫}`$ and $`\stackrel{}{𝔱}`$ the normal and tangential unit vectors to $`𝔣,`$ respectively. The decoupling theorem is proved.

## III Nonlinear flame stabilization in gravitational field

The proved theorem considerably widens the scope of issues in flame dynamics accessible for analytical investigation. One of the most important consequences of the decoupling theorem is that unlike the local cellular dynamics, the large scale flame dynamics can be investigated in the framework of perturbation expansion with respect to the flame front slope, provided that the external field exerts a stabilizing influence on the flame. This fact will be illustrated below in the case of a flame propagating in an initially quiescent fluid in the direction of the gravitational field. As is well known (see, e.g., Ref. zel ), in this case gravity plays the stabilizing role at the linear stage of development of the LD-instability. Our aim will be to determine the role of the nonlinear effects, and to explore the possibility of a full stabilization of the curved flame front by the gravitational field.

We will follow the general method of deriving weakly nonlinear equations for the flame front position, developed in Ref. kazakov3 . It consists in bringing the system of hydrodynamic equations together with the jump conditions at the flame front to the so-called transverse representation in which dependence of all flow variables on the coordinate in the direction of flame propagation ($`z`$) is rendered purely parametric, and then reducing this system to a single equation for the front position. The calculation in the presence of gravity is very similar to that in the case of a freely propagating flame. Therefore, derivation of the equation will be only sketched below, referring the reader to the work kazakov3  for more detail.

As was mentioned in the end of the preceding section, it is natural to take $`𝔘U_\mathrm{f}`$ as the velocity unit, redefining accordingly the units of length, time, and pressure to

$$\frac{\left(𝔘U_\mathrm{f}\right)^2}{\stackrel{}{g}}𝔏_0,\frac{𝔏_0}{𝔘U_\mathrm{f}},\left(\rho _\mathrm{u}𝔘U_\mathrm{f}\right)^2,$$

respectively. For simplicity, designation of the flow variables as well as space coordinates and time will be left unchanged. Then the bulk equations read666For brevity, the $`o(1)`$ symbols will be omitted in what follows.

$`\mathrm{div}𝔳`$ $`=`$ $`0,`$ (55)
$`\dot{\stackrel{}{𝔳}}+\left(\stackrel{}{𝔳}\stackrel{}{}\right)\stackrel{}{𝔳}`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{}𝔭}{\varrho }}+\stackrel{}{𝔊},\stackrel{}{𝔊}={\displaystyle \frac{𝒈𝔏_0}{𝔘^2}},\stackrel{}{𝔊}=1.`$ (56)
Since the fuel is assumed initially quiescent, the flow is potential upstream, and the general solution of Eqs. (55), (56) can be readily written down. In the reference frame of an initially plane flame front,

$`𝔲`$ $`=`$ $`1+{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}𝑑k𝔲_k\mathrm{exp}(|k|z+ikx),`$ (57)
$`𝔴`$ $`=`$ $`\widehat{\mathrm{H}}(𝔲1),`$ (58)
$`\dot{𝔲}`$ $`+`$ $`\widehat{\mathrm{\Phi }}\left(𝔭+𝔊z\right)+{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{2}}(𝔲^2+𝔴^2)=0,`$ (59)
$`𝔴`$ $``$ $`𝔳_x,𝔲𝔳_z,𝔊𝔊_z.`$
Here $`\widehat{\mathrm{H}}`$ denotes the Hilbert operator defined as

$`\widehat{\mathrm{H}}\mathrm{exp}(ikx)`$ $`=`$ $`i\mathrm{sign}(k)\mathrm{exp}(ikx),k0,`$ (60)
$`\mathrm{sign}(k)`$ $``$ $`{\displaystyle \frac{k}{|k|}}.`$
The LD-operator $`\widehat{\mathrm{\Phi }}`$ is related to $`\widehat{\mathrm{H}}`$ by $`\widehat{\mathrm{\Phi }}=\widehat{\mathrm{H}}/x.`$ Equation (59) is nothing but the Bernoulli equation written in the transverse form.

Because of the vorticity produced by the curved flame front, the flow of products of combustion is not potential. Nevertheless, the following transverse relation between the flow variables downstream can be obtained from Eqs. (55), (56) at the second order of nonlinearity

$`\dot{𝔲}\theta 𝔴^{}\widehat{\mathrm{\Phi }}\left(\theta 𝔭+{\displaystyle \frac{(𝔲\theta )^2+𝔴^2}{2}}\right)+𝔴\left(𝔲^{}+{\displaystyle \frac{\dot{𝔴}}{\theta }}+𝔭^{}\right)=0,`$ (61)
Equations (51)–(54), (58), (59), and (61) constitute the closed system of equations describing flame dynamics in the transverse representation. It can be reduced to the following equation for the function $`𝔣`$

$`(\theta +1)\ddot{𝔣}+2\theta \widehat{\mathrm{\Phi }}\dot{𝔣}+\theta (\theta 1)𝔣^{\prime \prime }+(\theta 1)𝔊\widehat{\mathrm{\Phi }}𝔣+\left(\theta {\displaystyle \frac{(\theta 1)^2}{2}}\right)\widehat{\mathrm{\Phi }}(𝔣^{})^2`$
$`+{\displaystyle \frac{(\theta 1)^2}{\theta }}𝔊(𝔣^{})^2+\left(\theta +{\displaystyle \frac{1}{\theta }}\right)\left(𝔣^{}\dot{𝔣}^{}+\dot{𝔣}^{}\widehat{\mathrm{H}}\dot{𝔣}\right)+{\displaystyle \frac{\theta 1}{2}}\widehat{\mathrm{\Phi }}\left(\dot{𝔣}^2+\left(\widehat{\mathrm{H}}\dot{𝔣}\right)^2\right)`$
$`+(3\theta 1)\widehat{\mathrm{\Phi }}\left(𝔣^{}\widehat{\mathrm{H}}\dot{𝔣}\right)+\left(2\theta 1+{\displaystyle \frac{1}{\theta }}\right)𝔣^{}\widehat{\mathrm{H}}\ddot{𝔣}{\displaystyle \frac{\theta 1}{\theta }}\left(\widehat{\mathrm{H}}\ddot{𝔣}+𝔊𝔣^{}\right)\widehat{\mathrm{H}}\dot{𝔣}=0.`$ (62)
The linear terms in this equation reproduce the well-known equation

$$(\theta +1)\ddot{𝔣}+2\theta \widehat{\mathrm{\Phi }}\dot{𝔣}+\theta (\theta 1)𝔣^{\prime \prime }+(\theta 1)𝔊\widehat{\mathrm{\Phi }}𝔣=0,$$

from which it follows that at the linear stage of development of the LD-instability, the gravitational field plays the stabilizing role in the case of flame propagation in the direction of $`\stackrel{}{𝔊}`$ ($`𝔊=+1`$), and destabilizing in the opposite case ($`𝔊=1`$). To determine the role of the nonlinear terms, let us assume that there exists a stationary regime of flame propagation. It should be stressed that this assumption concerns only the the large scale front structure described by the function $`𝔣.`$ The local cellular structure does not need to be stationary. Then Eq. (III) simplifies to

$`\theta (\theta 1)𝔣^{\prime \prime }+(\theta 1)𝔊\widehat{\mathrm{\Phi }}𝔣+\left(\theta {\displaystyle \frac{(\theta 1)^2}{2}}\right)\widehat{\mathrm{\Phi }}(𝔣^{})^2+{\displaystyle \frac{(\theta 1)^2}{\theta }}𝔊(𝔣^{})^2=0.`$ (63)
It is not difficult to see that the nonlinear term proportional to $`𝔊`$ exerts a stabilizing influence on the flame if $`𝔊=+1.`$ Indeed, if we take

$`𝔣(x,t)e^{\sigma t}\mathrm{sin}(kx),`$ (64)
Eq. (III) can be roughly considered as a “dispersion relation” for the increment $`\sigma .`$ As we see, the nonlinear term decreases $`\sigma `$ if $`𝔊=+1,`$ and vice versa. Therefore, Eq. (63) can only have solutions if $`𝔊=+1,`$ and the question of whether the flame can be stabilized by the gravitational field reduces to the question of existence of nontrivial solutions to this equation.777The last term in the non-stationary equation (III) is also proportional to $`𝔊.`$ Unlike the other two, it has a destabilizing effect on the flame front in the case $`𝔊=+1.`$ Indeed, using the definition of the Hilbert operator (60), and substituting expression (64) into this term, one finds
$$\frac{\theta 1}{\theta }𝔣^{}\widehat{\mathrm{H}}\dot{𝔣}=\frac{\theta 1}{\theta }\sigma |k|\mathrm{cos}^2(kx)<0.$$
This term, however, is irrelevant to the issue of existence of stationary configurations, since it contains $`\dot{𝔣}.`$ Whether the gravitational field has an overall stabilizing effect, or not, depends on solvability of the stationary equation (63).

Equation (63) is a nonlinear integro-differential equation with respect to $`\varphi =𝔣^{}.`$ To solve this equation, we first transform it as follows. Let us rewrite the nonlinear term $`\widehat{\mathrm{\Phi }}(𝔣^{})^2`$ iterating Eq. (63) with respect to $`𝔣^{\prime \prime },`$ i.e., substituting

$$𝔣^{\prime \prime }=\frac{𝔊}{\theta }\widehat{\mathrm{\Phi }}𝔣+O(𝔣^2).$$

One has

$$\widehat{\mathrm{\Phi }}(𝔣^{})^2\widehat{\mathrm{H}}\frac{(𝔣^{})^2}{\eta }=2\widehat{\mathrm{H}}(𝔣^{}𝔣^{\prime \prime })=\frac{2𝔊}{\theta }\widehat{\mathrm{H}}(𝔣^{}\widehat{\mathrm{\Phi }}𝔣)+O(𝔣^3).$$

Using the well-known identity

$`2\widehat{\mathrm{H}}\{\psi \widehat{\mathrm{H}}\psi \}=(\widehat{\mathrm{H}}\psi )^2\psi ^2,`$
we find

$$\widehat{\mathrm{\Phi }}(𝔣^{})^2=\frac{𝔊}{\theta }\{(\widehat{\mathrm{H}}𝔣^{})^2(𝔣^{})^2\}+O(𝔣^3).$$

Hence, within the accuracy of the second order, Eq. (63) takes the form

$`𝔣^{\prime \prime }+\alpha \widehat{\mathrm{\Phi }}𝔣+\beta (𝔣^{})^2\gamma \left(\widehat{\mathrm{\Phi }}𝔣\right)^2=0,`$ (65)
where

$$\alpha =\frac{𝔊}{\theta },\beta =\frac{𝔊}{\theta ^2(\theta 1)}\left(\theta +\frac{(\theta 1)^2}{2}\right),\gamma =\frac{𝔊}{\theta ^2(\theta 1)}\left(\theta \frac{(\theta 1)^2}{2}\right).$$

In connection with the transformation performed, it is worth to emphasize validity of the weak nonlinearity expansion when applied to the investigation of the large scale flame dynamics. As was shown in detail in Refs. kazakov1 ; kazakov2 , this expansion turns out to be self-contradictory in the case of a freely propagating stationary flame treated in the framework of the thin front model. This fact can be seen directly from Eq. (63) with $`𝔊=0,`$ in which case this equation reduces to the equality of two quantities of apparently different orders – $`𝔣^{\prime \prime }`$ and $`\widehat{\mathrm{\Phi }}(𝔣^{})^2.`$888Perhaps, it is worth to stress once more that despite similarity of Eq. (III) with $`𝔊=0`$ to that obtained in Ref. kazakov3 , its meaning is completely different. In the notation of Sec. II.3, the latter equation determines the function $`f_1^{(0)},`$ while Eq. (III) – the function $`f_0^{(0)}.`$ Only if $`\theta 1`$ does the weak nonlinearity expansion of stationary flames make sense, since then $`𝔣^{}=O(\theta 1),`$ $`𝔣^{\prime \prime }=O((\theta 1)^2),`$ so both terms in Eq. (63) are $`O((\theta 1)^3)`$ quantities. On the contrary, in the presence of the gravitational field, $`𝔣^{}`$ can be treated as the first order quantity when $`(\theta 1)`$ is not small, because $`\theta 𝔣^{\prime \prime }`$ and $`𝔊\widehat{\mathrm{\Phi }}𝔣`$ are of the first order in this case, and nothing prevents their sum from being formally a second order quantity. In fact, $`(\theta 1)`$ must be finite in the latter case, since the second order term $`(\theta 1)𝔊\widehat{\mathrm{\Phi }}𝔣`$ would be the only second order term in Eq. (63) otherwise. As we will see below, solutions to this equation turn out to be unbounded for $`\theta 1.`$

Turning back to Eq. (65), let us show first of all that its solutions, if any, are non-periodic. Notice that if $`𝔣(x)`$ is a periodic function, then $`[\widehat{\mathrm{H}}(𝔣\overline{𝔣})](x),`$ where $`\overline{𝔣}`$ is the mean value of $`𝔣,`$ is also periodic with the same period \[substraction of $`\overline{𝔣}`$ is necessary, because $`\widehat{\mathrm{H}}`$ is undefined on constants, see Eq. (60)\]. Let us integrate Eq. (65) over period. The first two terms in this equation give rise to zero:

$`{\displaystyle 𝑑x𝔣^{\prime \prime }}`$ $`=`$ $`𝔣^{}|=0,`$
$`{\displaystyle 𝑑x\widehat{\mathrm{\Phi }}𝔣}`$ $`=`$ $`{\displaystyle 𝑑x\widehat{\mathrm{H}}𝔣^{}}={\displaystyle 𝑑x\widehat{\mathrm{H}}(𝔣\overline{𝔣})^{}}`$
$`=`$ $`{\displaystyle }dx\left\{\widehat{\mathrm{H}}(𝔣\overline{𝔣})\right\}^{}=\widehat{\mathrm{H}}(𝔣\overline{𝔣})|=0.`$
On the other hand, using unitarity of the Hilbert operator, we find

$`{\displaystyle 𝑑x\left\{\beta (𝔣^{})^2\gamma \left(\widehat{\mathrm{H}}𝔣^{}\right)^2\right\}}={\displaystyle 𝑑x(\beta \gamma )(𝔣^{})^2}>0.`$
Of course, the absence of periodic solutions could be inferred already from Eq. (63). It is seen that such solutions are forbidden by the positive definite nonlinear term proportional to $`𝔊.`$

Non-periodic solutions can be found in the form of the pole decomposition for the function $`\varphi =𝔣^{}`$

$`\varphi =a{\displaystyle \underset{k=1}{\overset{2P}{}}}{\displaystyle \frac{1}{xx_k}},`$ (66)
where the value of the amplitude $`a,`$ as well as position of $`P`$ pares of the complex conjugate poles $`x_k,`$ are to be determined by substituting this decomposition into Eq. (65). Using the definition of Hilbert operator, one can show that

$`\widehat{\mathrm{H}}\varphi `$ $`=`$ $`ia{\displaystyle \underset{k=1}{\overset{2P}{}}}{\displaystyle \frac{\mathrm{sign}(\mathrm{Im}x_k)}{xx_k}}.`$
It is not difficult to verify that Eq. (65) is satisfied by (66) provided that

$$a=\frac{1}{\beta +\gamma },$$

and $`x_k`$ satisfy the following system of algebraic equations

$`i\alpha (\beta +\gamma )+2{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{l=1}{lk}}{\overset{2P}{}}}{\displaystyle \frac{\beta \mathrm{sign}(\mathrm{Im}x_k)+\gamma \mathrm{sign}(\mathrm{Im}x_l)}{x_kx_l}}=0,k=1,\mathrm{},2P.`$ (67)
Evidently, this system is only consistent if $`𝔊=+1.`$ Indeed, in the case of a pare of complex conjugate poles $`x_1,`$ $`x_2=x_1^{},`$ one has, assuming $`\mathrm{Im}x_1>0,`$

$$\mathrm{Im}x_1=\frac{\beta \gamma }{\alpha (\beta +\gamma )}=\frac{(\theta 1)^2}{2𝔊},$$

which is consistent with the assumed positivity of $`\mathrm{Im}x_1`$ if $`𝔊=+1.`$ It is not difficult to show that the same is true in the general case of arbitrary number of poles. Consider the imaginary part of Eq. (67) corresponding to the pole uppermost in the complex plane, and take into account that $`\beta >\gamma .`$ The fact that the system (67) is inconsistent for $`𝔊=1`$ does not mean, of course, that the nonlinear stabilization is impossible in this case. Investigation of such a possibility requires account of higher order corrections.

In the case of $`P=1,`$ one has, furthermore,

$$𝔣^{}=\frac{1}{\beta +\gamma }\left(\frac{1}{xx_1}+\frac{1}{xx_1^{}}\right)=\frac{2}{\beta +\gamma }\frac{(x\mathrm{Re}x_1)}{(x\mathrm{Re}x_1)^2+(\mathrm{Im}x_1)^2},$$

and therefore,

$`𝔣={\displaystyle \frac{\theta (\theta 1)}{2}}\mathrm{ln}\left\{(xx_0)^2+{\displaystyle \frac{(\theta 1)^4}{4}}\right\},x_0\mathrm{Re}x_1.`$
The two-pole solutions for the cases $`\theta =5,10`$ and $`x_0=0`$ are shown in Fig. 2.

Equation (65) is derived in the scope of the power expansion with respect to the flame front slope $`\varphi .`$ In the case of the two-pole solution, $`|\varphi |`$ takes its maximal value

$$\varphi _\mathrm{m}=\frac{\theta }{\theta 1}$$

at the points $`x_\mathrm{m}=x_0\pm \mathrm{Im}x_1.`$ We see that the developed weak nonlinearity expansion is valid if $`\theta `$ is not too close to unity. For realistic values of the expansion coefficient ($`\theta =510`$) $`\varphi _\mathrm{m}1.`$

Next, consider the four-pole solution ($`P=2`$). It has the form

$`𝔣={\displaystyle \frac{\theta (\theta 1)}{2}}\mathrm{ln}\left\{\left[(x\mathrm{Re}x_1)^2+\left(\mathrm{Im}x_1\right)^2\right]\left[(x\mathrm{Re}x_2)^2+\left(\mathrm{Im}x_2\right)^2\right]\right\}.`$
Assuming that $`\mathrm{Im}x_{1,2}>0,`$ one has the following equations for the position of poles $`x_1,`$ $`x_2,`$ $`x_3=x_1^{},`$ $`x_4=x_2^{}`$

$`i\alpha (\beta +\gamma )+2\left\{{\displaystyle \frac{\beta +\gamma }{x_1x_2}}+(\beta \gamma )\left({\displaystyle \frac{1}{2i\mathrm{Im}x_1}}+{\displaystyle \frac{1}{x_1x_2^{}}}\right)\right\}`$ $`=`$ $`0,`$
$`i\alpha (\beta +\gamma )+2\left\{{\displaystyle \frac{\beta +\gamma }{x_2x_1}}+(\beta \gamma )\left({\displaystyle \frac{1}{2i\mathrm{Im}x_2}}+{\displaystyle \frac{1}{x_2x_1^{}}}\right)\right\}`$ $`=`$ $`0.`$
Separating the real and imaginary parts, and rearranging yields three equations for the four quantities $`\mathrm{Re}x_{1,2},`$ $`\mathrm{Im}x_{1,2}`$

$`2\alpha {\displaystyle \frac{\beta +\gamma }{\beta \gamma }}\left\{{\displaystyle \frac{1}{\mathrm{Im}x_1}}+{\displaystyle \frac{1}{\mathrm{Im}x_2}}+4{\displaystyle \frac{\mathrm{Im}(x_1+x_2)}{|x_1x_2^{}|^2}}\right\}`$ $`=`$ $`0,`$ (68)
$`4{\displaystyle \frac{\beta +\gamma }{\beta \gamma }}{\displaystyle \frac{\mathrm{Im}(x_2x_1)}{|x_1x_2|^2}}+{\displaystyle \frac{1}{\mathrm{Im}x_2}}{\displaystyle \frac{1}{\mathrm{Im}x_1}}`$ $`=`$ $`0,`$ (69)
$`\mathrm{Re}(x_1x_2)\left({\displaystyle \frac{\beta +\gamma }{|x_1x_2|^2}}+{\displaystyle \frac{\beta \gamma }{|x_1x_2^{}|^2}}\right)`$ $`=`$ $`0.`$ (70)
It follows from Eq. (70) that $`\mathrm{Re}x_1=\mathrm{Re}x_2.`$ This solution describes the “confluence” of poles. Then the remaining Eqs. (68), (69) give

$`\mathrm{Im}x_{1,2}={\displaystyle \frac{1}{\alpha }}\left(1+2{\displaystyle \frac{\beta \gamma }{\beta +\gamma }}\right)\left(1\pm \sqrt{{\displaystyle \frac{\beta +\gamma }{2\beta }}}\right)=\left(\theta ^2\theta +1\right)\left(1\pm \sqrt{{\displaystyle \frac{2\theta }{\theta ^2+1}}}\right).`$
The two-pole and four-pole solutions are compared in Fig. 3 in the case of $`\theta =8`$ and $`x_0=0.`$

The pole confluence is in fact a common property of the solutions (66). To see this, let us take the real part of the equation with $`k`$ corresponding to the rightmost pole in the upper half-plane. We have

$`{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{l=1}{lk}}{\overset{2P}{}}}\left(\mathrm{Re}x_k\mathrm{Re}x_l\right){\displaystyle \frac{\beta \pm \gamma }{|x_kx_l|^2}}=0.`$
In view of the choice of $`k,`$ the left hand side is the sum of non-negative terms. It can be zero only if $`\mathrm{Re}x_k\mathrm{Re}x_l=0`$ for all $`l.`$

The question of which configuration is realized in the given conditions requires carrying out the stability analysis of various pole solutions, and can be solved, of course, only on the basis of the general non-stationary Eq. (III). According to the definition of $`𝔣,`$ such an analysis is to be performed with respect to perturbations with wavelengths $`\lambda L_0.`$

## IV Discussion and conclusions

The large scale flame dynamics are independent of its local cellular structure in zero order approximation with respect to the flame front thickness. This is the main result of the work, proved in Sec. II. The local flame corrugation only affects the value of the normal velocity $`U_\mathrm{f}`$ changing it to $`𝔘U_\mathrm{f},`$ where $`𝔘>1`$ describes increase of the front length due to its wrinkling. In the scope of the thin front model, $`U_\mathrm{f}`$ plays the role of an external parameter specifying the characteristic velocity of the problem under consideration. Thus, the overall effect of the local flame structure on its large scale evolution amounts to a renormalization of this parameter. For flames of practical importance, the $`𝔘`$-factor is about $`1.31.5.`$ In fact, it is $`𝔘U_\mathrm{f},`$ rather than $`U_\mathrm{f},`$ which is more convenient to measure experimentally, since the measurement of $`U_\mathrm{f}`$ requires special facilities to suppress development of the LD-instability, such as those used in Ref. clanet .

The decoupling theorem allows one to avoid the difficult issues arising in investigating flame dynamics at length scales of the order $`\lambda _\mathrm{c},`$ and to go directly to scales characterizing the problem in question. This is particularly important in numerical simulations of the flame dynamics. The computational grid should be chosen so as to well resolve the flame cellular structure, which leaves a little space for investigation of larger scales because of the limitations of computational facilities.

The decoupling theorem also opens the way for analytical investigation of the large scale flame dynamics. As an example, the nonlinear development of the LD-instability in the presence of the gravitational field was considered in Sec. III, where a weakly nonlinear non-stationary equation for the flame front position was obtained \[Eq. (III)\]. This equation admits stationary solutions in the case of flame propagation in the direction of the field, which means that the gravitational field has a stabilizing overall effect in this case. The resulting stationary flame configuration turns out to be essentially non-periodic, and represents a symmetrical “hump” in the direction of the flame propagation, with slowly decreasing logarithmic “tails.” A complete investigation of the non-stationary equation will be given elsewhere.


text:
# Rare 𝐵 Decays into States Containing a 𝐽/𝜓 Meson and a Meson with 𝑠𝑠̄ Quark Content

BABAR-PUB-02/007

SLAC-PUB-9262

hep-ex/0000000

The BABAR Collaboration

## Abstract

We report a study of the $`B`$ meson decays, $`B^+J/\psi \varphi K^+`$, $`B^0J/\psi \varphi K_S^0`$, $`B^0J/\psi \varphi `$, $`B^0J/\psi \eta `$ and $`B^0J/\psi \eta ^{}`$ using 56 million $`B\overline{B}`$ events collected at the $`\mathrm{{\rm Y}}(4S)`$ resonance with the BABAR detector at the PEP-II $`e^+e^{}`$ asymmetric-energy storage ring. We measure the branching fractions $`(B^+J/\psi \varphi K^+)=(4.4\pm 1.4(stat)\pm 0.5(syst))\times 10^5`$ and $`(B^0J/\psi \varphi K_S^0)=(5.1\pm 1.9(stat)\pm 0.5(syst))\times 10^5`$, and set upper limits at 90$`\%`$ confidence level for the branching fractions $`(B^0J/\psi \varphi )<9.2\times 10^6`$, $`(B^0J/\psi \eta )<2.7\times 10^5`$, and $`(B^0J/\psi \eta ^{})<6.3\times 10^5`$.

preprint: BABAR-PUB-02/007preprint: SLAC-PUB-9262

Recent observations of the $`B`$ meson decays $`BJ/\psi \pi `$ babar-charmonium  and $`J/\psi \rho `$ babar-jpsirho  are evidence for the Cabibbo-suppressed transition $`bc\overline{c}d`$ via the color-suppressed diagram shown in Fig. 1 (a). Here we present a search for color-suppressed modes with hidden strangeness, $`s\overline{s}`$, in the final state: $`BJ/\psi \eta ,J/\psi \eta ^{},J/\psi \varphi `$ and $`J/\psi \varphi K`$. The decays $`B^0J/\psi \eta `$ and $`B^0J/\psi \eta ^{}`$ occur via the same diagram, Fig. 1 (a), and should have a rate comparable to $`BJ/\psi \pi `$. If large enough samples can be isolated, these $`CP`$ eigenstates could be used to test $`CP`$ violation beneke . Models based on the heavy quark factorization approximation by A. Deandrea et al. deandrea  are used to predict that the branching fraction for $`B^0J/\psi \eta `$ is a factor of 4 smaller than that for $`B^0`$ $`J/\psi \pi ^0`$. Assuming that the decay $`B^0J/\psi \varphi `$ is a color-suppressed mode with rescattering as shown in Fig. 1 (b), then the absence of a signal would indicate that the rescattering effects are negligible. The decay $`BJ/\psi \varphi K`$ is a Cabibbo-allowed and color-suppressed decay via the transition $`b\overline{q}c\overline{c}s\overline{s}s\overline{q}`$, where the $`s\overline{s}`$ quark pairs are produced from sea quarks or are connected via gluons as shown in Figs. 1 (c) and (d), respectively. This particular three-body decay would be of interest in the search for hybrid charmonium states that decay to the final state $`J/\psi \varphi `$ close . In this paper we report on branching fractions or upper limits for $`J/\psi \eta `$, $`J/\psi \eta ^{}`$, $`J/\psi \varphi `$, $`J/\psi \varphi K^+`$, and $`J/\psi \varphi K_S^0`$.

The data used in this analysis were collected at the PEP-II asymmetric-energy $`e^+e^{}`$ storage ring with the BABAR detector, fully described elsewhere babar-det . The BABAR detector contains a five-layer silicon vertex tracker (SVT) and a forty-layer drift chamber (DCH) in a 1.5-T solenoidal magnetic field. These devices detect charged particles and measure their momentum and energy loss. Photons and neutral hadrons are detected in a CsI(Tl) crystal electromagnetic calorimeter (EMC). The EMC detects photons with energies as low as 20 MeV and identifies electrons by their energy deposition. An internally reflecting ring-imaging Cherenkov detector (DIRC) of quartz bars is dedicated to charged particle identification (PID). Penetrating muons and neutral hadrons are identified by the steel flux return (IFR), which is instrumented with 18-19 layers of resistive plate chambers.

The data correspond to a total integrated luminosity of $`50.9`$ fb-1 taken on the $`\mathrm{{\rm Y}}(4S)`$ resonance and 6.3 fb-1 taken off-resonance at an energy 0.04 GeV below the $`\mathrm{{\rm Y}}(4S)`$ mass and below the threshold for $`B\overline{B}`$ production. In this sample, there are $`55.5\pm 0.6`$ million $`B\overline{B}`$ events ($`N_{B\overline{B}}`$).

In this analysis, all charged track candidates are required to have at least 12 DCH hits and transverse momentum greater than 100 MeV/$`c`$. The track candidates not associated with a $`K_S^0`$ decay must also originate near the nominal beam spot. The muon, electron, and kaon candidates must have a polar angle in radians of $`0.3<\theta _\mu <2.7`$, $`0.410<\theta _e<2.409`$, and $`0.45<\theta _K<2.50`$, respectively. In addition, all charged kaon candidates are required to have a laboratory momentum greater than 250 MeV/$`c`$. These requirements ensure the selection of tracks in the regions where the acceptance is well understood by the PID systems.

Photon candidates are identified from energy deposited in contiguous EMC crystals, summed together to form a cluster with total energy greater than 30 MeV and a shower shape consistent with that expected for electromagnetic showers.

Electron candidates are required to have a good match between the expected and measured energy loss ($`\mathrm{d}E/\mathrm{d}x`$) in the DCH, and between the expected and measured Cherenkov angle in the DIRC. The measurements of the ratio of EMC shower energy to DCH momentum, and the number of EMC crystals associated with the track candidate must be appropriate for an electron.

Muons are selected based on the energy deposited in the EMC, the number and distribution of hits in the IFR, the match between the IFR hits and the extrapolation of the DCH track into the IFR, and the depth of penetration of the track into the IFR.

Charged kaon and pion candidates are selected based on energy loss information from the SVT and DCH and the Cherenkov angle measured by the DIRC.

The intermediate states in the indicated decay modes used in this analysis, $`J/\psi (ee,\mu \mu )`$, $`\varphi \left(K^+K^{}\right)`$, $`\eta (\gamma \gamma ,\pi ^+\pi ^{}`$$`\pi ^0)`$, $`\eta ^{}\left(\eta \left(\gamma \gamma \right)\pi ^+\pi ^{}\right)`$, $`\pi ^0\left(\gamma \gamma \right)`$, and $`K_S^0\left(\pi ^+\pi ^{}\right)`$, are selected with the mass intervals in Table 1. Since $`B^0J/\psi \eta `$ and $`B^0J/\psi \eta ^{}`$ involve decays of a pseudoscalar meson into a vector and a pseudoscalar meson, the angular distribution is proportional to $`\mathrm{sin}^2\theta _{\mathrm{}}`$, where $`\theta _{\mathrm{}}`$ is the helicity angle 111 In the reaction, $`ZX+Y,Xa+b`$, the helicity angle of particle $`a`$ is defined as the angle measured in the particle $`X`$ rest frame between the direction of particle $`a`$ and the direction opposite to particle $`Z`$. of the lepton from the $`J/\psi `$. Hence an additional requirement of $`\left|\mathrm{cos}\theta _{\mathrm{}}\right|<0.8`$ is applied to reject continuum and other backgrounds. The $`\eta `$ candidates are rejected if either of the associated photons, in combination with any other photon in the event, forms a $`\gamma \gamma `$ mass within 20 MeV/$`c^2`$ of the $`\pi ^0`$ mass. For the mode $`B^0J/\psi \eta (\gamma \gamma )`$, the $`\eta `$ candidate is required to have $`\left|\mathrm{cos}\theta _\gamma ^\eta \right|<0.8`$, where $`\theta _\gamma ^\eta `$ is the photon helicity angle in the $`\eta `$ rest frame. This rejects combinatoric background due to random pairs of photons that typically have a photon helicity angle that peaks at 0 or 180 degrees. For the $`\eta ^{}\eta \left(\gamma \gamma \right)\pi ^+\pi ^{}`$ candidates, we use the same $`\eta `$ selection criteria for the $`\eta `$ described above, including the $`\pi ^0`$ veto.

An additional requirement is applied to separate two-jet continuum events from the more spherical $`B`$ meson decays. The angle $`\theta _T`$ between the thrust direction of the $`B`$ meson candidate and the thrust direction of the remaining tracks in the event is calculated. We require $`\left|\mathrm{cos}\theta _T\right|<0.8`$, since these thrust axes are uncorrelated and the distribution in $`\mathrm{cos}\theta _T`$ is flat for $`B\overline{B}`$ events, while the distribution is peaked at $`\mathrm{cos}\theta _T=\pm 1`$ for continuum events.

The intermediate candidates are combined to construct the $`B`$ candidates for the six decay modes under study. The estimation of the signal and the background employs two kinematic variables: the energy difference $`\mathrm{\Delta }E`$ between the energy of the $`B`$ candidate and the beam energy $`E_b^{}`$ in the $`\mathrm{{\rm Y}}(4S)`$ rest frame; and the energy-substituted mass $`m_{\mathrm{ES}}=\sqrt{\left(E_b^{}\right)^2\left(P_B^{}\right)^2}`$, where $`P_B^{}`$ is the reconstructed momentum of the $`B`$ candidate in the $`\mathrm{{\rm Y}}(4S)`$ frame. Typically these two weakly correlated variables form a two-dimensional Gaussian distribution for the $`B`$ meson signal but not for background. The resolutions in $`\mathrm{\Delta }E`$ and $`m_{\mathrm{ES}}`$ are decay mode dependent. A signal region for each mode is defined as a rectangular region in the $`\mathrm{\Delta }E`$ versus $`m_{\mathrm{ES}}`$ plane, listed in Table 2. The $`m_{\mathrm{ES}}`$ range is given in term of $`m_{\mathrm{ES}}m_B`$, where $`m_B`$ is the mass of $`B`$ meson. The number of data events, $`n_0`$, observed in the signal region for each mode is listed in Table 2.

The efficiencies for each mode are determined by Monte Carlo simulation. The simulations of $`J/\psi \varphi K`$ and $`J/\psi \varphi `$ decays assumed three- and two-body phase space, respectively, with unpolarized $`J/\psi `$ and $`\varphi `$ decays. The $`J/\psi \eta `$ and $`J/\psi \eta ^{}`$ simulations used the angular correlations determined by the helicity amplitude.

The backgrounds in the $`m_{\mathrm{ES}}`$ distribution are composed of two components: a combinatoric background, whose shape is described by an ARGUS function argus-function , and a peaking background that peaks in the signal region and is described by a Gaussian function. The sources of combinatoric background are the continuum events and two categories of $`B\overline{B}`$ events: decays with a leptonic $`J/\psi `$ decay, and those without. Monte Carlo simulation studies show that the source of the peaking background is $`B\overline{B}`$ events that contain a leptonic $`J/\psi `$ decay.

The shape of the ARGUS function is determined mode by mode by fitting to the $`m_{\mathrm{ES}}`$ distribution of candidates in an enhanced fake $`J/\psi `$ sample, which is obtained by reversing the normal lepton identification requirements.

The normalization of the combinatoric background for each mode is obtained from a fit to the $`m_{\mathrm{ES}}`$ distributions in the $`\mathrm{\Delta }E`$ signal region of the on-peak data. The integral of the ARGUS function in the signal region is $`n_C`$, the number of combinatoric background events.

The peaking background is determined from a fit to the $`m_{\mathrm{ES}}`$ distribution of Monte Carlo $`B\overline{B}`$ events with leptonic $`J/\psi `$ decays using the sum of a Gaussian and an ARGUS function. The number of peaking background events $`n_P`$ is the integral of the Gaussian function in the signal region.

The total number of background events ($`n_b`$) and the uncertainty on this number ($`\sigma _b`$) are calculated from the fit value of $`n_C`$ and $`n_P`$ and their errors. The values of $`n_b`$ and $`\sigma _b`$ are listed in Table 2 for all modes. The combinatoric background is by far the dominant background in all modes except the $`B^0J/\psi \eta (\pi ^+\pi ^{}\pi ^0)`$ mode, where the peaking component is $`20\%`$ of the total background.

In Table 3, we list the contribution to the systematic error from the uncertainty on each of the following quantities: $`N_{B\overline{B}}`$; secondary branching fractions pdg ; Monte Carlo statistics; PID, tracking, and photon detection efficiencies, which are based on the study of control samples; and background parameterization, which is estimated using $`\mathrm{\Delta }E`$ sideband information.

Additional systematic uncertainties due to the decay model dependence are estimated for the modes $`J/\psi \varphi `$, $`J/\psi \varphi K^+`$, and $`J/\psi \varphi K_S^0`$. Monte Carlo simulations are used to determine how much the efficiency depends on assumptions about intermediate resonances and angular distributions. Two samples are generated for each of the three modes with decay distributions determined by the assumed polarization of the vector daughter mesons, rather than by phase space. One sample is generated with $`100\%`$ transversely polarized $`J/\psi `$ and $`\varphi `$ mesons, and the other with $`100\%`$ longitudinally polarized $`J/\psi `$ and $`\varphi `$ mesons. The resulting relative change in efficiency is entered as a fractional systematic error in Table 3. An additional check based on Monte Carlo samples with an intermediate state gives negligible effect.

The total systematic error for each mode combines all these separate errors in quadrature and is listed (Total) in Table 3.

There is evidence for signals in the $`J/\psi \varphi K^+`$ and $`J/\psi \varphi K_S^0`$ modes. The results are shown in Figs. 2 and 3. The Poisson probability that the background $`n_b`$ fluctuates up to the observed number of events, $`n_0`$, or higher is $`7.7\times 10^6`$ for $`J/\psi \varphi K^+`$ and $`4.2\times 10^5`$ for $`J/\psi \varphi K_S^0`$. The branching fraction for these modes is determined by a simple subtraction of events in the signal region that yields the number of signal events, $`n_s=n_0n_b`$. The calculation of the branching fraction is based on the efficiency, $`n_s`$, $`N_{B\overline{B}}`$, and the secondary branching fractions for the $`J/\psi ,`$ $`\varphi `$, and $`K_S^0`$ from Ref. pdg . The results are summarized in Table 2 where the first error is the statistical error and the second is the systematic error, listed in Table 3. The derived result for $`B^0J/\psi \varphi K^0`$ is also shown in Table 2.

For modes with no signal or limited statistical evidence ($`J/\psi \varphi `$, $`J/\psi \eta `$, $`J/\psi \eta ^{}`$), we determine both a central confidence interval and an upper limit interpretation for the branching fraction. The upper limit method uses $`n_0`$, $`n_b`$, and $`\sigma _b`$, in the signal region, and the total systematic uncertainty $`\sigma _T`$. Assuming the two uncertainties ($`\sigma _b,\sigma _T`$) are uncorrelated and Gaussian, the Bayesian upper limit on the number of events ($`N_{90\%}`$) is obtained by folding the Poisson distribution with two normal distributions for these two uncertainties and integrating it to the $`90\%`$ confidence level (C.L.). In Table 2 we list for each mode the efficiency, the number of observed events, the expected number of background events, the $`90\%`$ C.L. upper limit for observed events, the corresponding branching fraction limit and a central interval for the branching fraction. The upper limit obtained from the combination of the two $`B^0J/\psi \eta `$ modes is shown in Table 2. The mean value of the branching fraction is calculated for $`B^0J/\psi \varphi `$ and $`B^0J/\psi \eta ^{}`$. We also combine the observed numbers of events for the two $`B^0J/\psi \eta `$ modes to calculate a branching fraction of $`(1.6\pm 0.6(stat.)\pm 0.1(syst.))\times 10^5`$. The Poisson probability that the background fluctuates up to the observed number of events or higher is $`2.5\times 10^5`$ for the combined result.

In summary, we determine the branching fraction of $`B`$ $`J/\psi \varphi K`$ in two modes, $``$($`B^+J/\psi \varphi K^+`$)= $`(4.4\pm 1.4\pm 0.5)\times 10^5`$ and $``$($`B^0J/\psi \varphi K_S^0`$)=$`(5.1\pm 1.9\pm 0.5)\times 10^5`$. The branching fraction of $`BJ/\psi \varphi K`$ is consistent with a CLEO cleo  result, $`\left(8.8_{3.0}^{+3.5}\pm 1.3\right)\times 10^5`$. Upper limits have been determined for the modes $`B^0J/\psi \varphi ,`$ $`J/\psi \eta `$, and $`J/\psi \eta ^{}`$. The upper limit on $`B^0J/\psi \eta `$ is a significant improvement over the previous best limit of $`<1.2\times 10^3`$ at $`90\%`$ C.L., from the L3 Collaboration L3 . The combined branching fraction for $`B^0J/\psi \eta `$ is comparable to the $`B^0J/\psi \pi ^0`$ branching fraction babar-charmonium . The search and resulting upper limits on the branching fractions for $`B^0J/\psi \eta ^{}`$ and $`B^0J/\psi \varphi `$ are presented.

We are grateful for the excellent luminosity and machine conditions provided by our PEP-II colleagues, and for the substantial dedicated effort from the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and kind hospitality. This work is supported by DOE and NSF (USA), NSERC (Canada), IHEP (China), CEA and CNRS-IN2P3 (France), BMBF and DFG (Germany), INFN (Italy), NFR (Norway), MIST (Russia), and PPARC (United Kingdom). Individuals have received support from the A. P. Sloan Foundation, Research Corporation, and Alexander von Humboldt Foundation.