*2.3. Mean-Field Approach*

The explanatory potential of the model is first examined by means of a mean-field approximation which mainly relies on the replacement of fluctuating quantities by space-averaged values and the neglect of correlations. The observables involved in the mean-field expressions are the ensemble averages of the microscopic states *S*(*i*, *j*). Their values at *t* + 1 are obtained by taking averages of the governing Equation (2) using the expression of the configurational probabilities given in (3). By assumption/definition *S* is the mean outcome of *pCx* averaged over all the possible configurations, where space dependence (*i*, *j*) is temporarily kept: *S* (*i*, *j*) = *p*[*S*1,*S*2,*S*¯ 4,*S*5]. This gives a set of two equations:

$$
\langle B'(i,j)\rangle \quad = 1 - \left\langle (1 - p\_1 B(i,j))(1 - p\_2 B(i,j+1))(1 - p'\_4 B(i-1,j))(1 - p\_5 B(i-1,j)) \right\rangle, \tag{4}
$$

$$\left< R'(i,j) \right> \quad = \quad 1 - \left< (1 - p\_1 R(i,j))(1 - p\_2 R(i+1,j))(1 - p'\_4 B(i,j-1))(1 - p\_5 R(i,j-1)) \right>. \tag{5}$$

The approximation now enters the evaluation of the products on the right hand side of the equation. Each variable is replaced by its average and the spatial dependence is dropped: *B*(*i*, *j*) → *B* and *R*(*i*, *j*) → *R*. Further, correlations are neglected so that the average of a product is just the product of averages. The expansions of (4) and (5) in powers of *B* and *R* are then readily obtained. Forgetting for a moment the intricacy linked to transversal splitting/collisions, the general expression for the dummy variables *S* and *S* reads:

$$
\langle S' \rangle = \sum\_{k} p\_k \langle S\_k \rangle - \sum\_{k\_1, k\_2} p\_{k\_1} p\_{k\_2} \langle S\_{k\_1} \rangle \langle S\_{k\_2} \rangle + \text{h.o.t.} \tag{6}
$$

with *pk* ∈ {*p*1, *p*2, *p* <sup>4</sup>, *p*5} and where h.o.t. stands for the higher order terms, formally cubic, quartic, etc. The first sum in (6) corresponds to the contribution of the elementary configurations introduced in Figure 6, and the second sum to binary configurations, in particular the nontrivial ones corresponding to transversal splittings and collisions examined in Figure 8 (right). Orders of magnitude among the *pk*, further support neglecting the contribution of configurations populated with three or more active sites, involving products of three or more probabilities *pk*, and among contributions of a given degree, those not containing *p*<sup>5</sup> when compared to those that do, recalling the assumption *p*<sup>5</sup> 1 and {*p*1, *p*2} 1 implied by the nearly deterministic propagation of states in position 5 of Figure 6. A number of terms can, therefore, be neglected in the expanded forms of (4) and (5), which after simplification read:

$$
\langle B' \rangle\_{\phantom{\cdot}} = \left( p\_1 + p\_2 + p\_5 \right) \langle B \rangle + p\_4' \langle R \rangle - p\_5 (p\_1 + p\_2) \langle B \rangle^2 - p\_5^2 \langle B \rangle \langle R \rangle \,\tag{7}
$$

$$
\langle R' \rangle \quad = \left( p\_1 + p\_2 + p\_5 \right) \langle R \rangle + p\_4' \langle B \rangle - p\_5 (p\_1 + p\_2) \langle R \rangle^2 - p\_5^2 \langle R \rangle \langle B \rangle \,. \tag{8}
$$

This system presents itself as the discrete time counterpart of the differential system introduced in [14] to interpret the symmetry-breaking bifurcation observed at decreasing Re in the simulations. As a matter of fact, subtracting *B* and *R* on both sides of (7) and (8) respectively, one gets:

$$<\langle B' \rangle - \langle B \rangle \approx \frac{\mathbf{d} \langle B \rangle}{(\mathbf{d}t \equiv 1)} = (p\_1 + p\_2 + p\_5 - 1) \langle B \rangle + \dots \tag{9}$$

$$
\langle \langle R' \rangle - \langle R \rangle \approx \frac{\mathbf{d} \langle R \rangle}{(\mathbf{d}t \equiv 1)} = (p\_1 + p\_2 + p\_5 - 1) \langle R \rangle + \dots \tag{10}
$$

to be compared with system (1,2) in [14], reproduced here for convenience:

$$\frac{d\mathbf{X}\_{+}}{dt} = a\mathbf{X}\_{+} + c\mathbf{X}\_{-} - b\mathbf{X}\_{+}^{2} - dX\_{+}\mathbf{X}\_{-},\tag{11}$$

$$\frac{d\mathbf{X}\_{-}}{dt} = a\mathbf{X}\_{-} + c\mathbf{X}\_{+} - b\mathbf{X}\_{-}^{2} - dX\_{-}\mathbf{X}\_{+},\tag{12}$$

where *X*<sup>±</sup> represents what are now the densities *B* and *R*. The coefficients in (11) and (12) are then related to the probabilities introduced in the model as *a* ∝ *p*<sup>1</sup> + *p*<sup>2</sup> + *p*<sup>5</sup> − 1, *b* ∝ *p*5(*p*<sup>1</sup> + *p*2), *c* ∝ *p* 4, and *d* ∝ *p*<sup>2</sup> <sup>5</sup>. By omitting the common proportionality constant that accounts for the time-stepping

inherent in the discrete time reduction (featured by the denominator of left-hand sides in (9) and (10) as "(d*t* ≡ 1))," constants *a*, *b*, *c*, and *d* will serve as short-hand notation for the corresponding full expressions in terms of the probabilities *pk*.

Since fixed points given by the condition *S* = *S* is strictly equivalent to d*X*±/d*t* = 0, we can next take advantage of the analysis performed in [14] and predict a supercritical symmetry-breaking bifurcation for an order parameter |*B*−*R*| (denoted "*A*" in [14]) at a threshold given by *<sup>c</sup>*cr = *<sup>a</sup>*(*<sup>d</sup>* − *<sup>b</sup>*)/(*<sup>d</sup>* + <sup>3</sup>*b*). This symmetry-breaking bifurcation takes place for *<sup>p</sup>* <sup>4</sup> = *c* > 0, but the model can deal with the regime below event A at Re ≈ 800 for which *p* <sup>4</sup> ≡ 0. In that case the bifurcation corresponding to global decay at Reg takes the form of two coupled equations generalizing (1) for DP. Using the abridged notation, these equations read:

$$
\langle B' \rangle = (a+1)\langle B \rangle - b \langle B \rangle^2 - d \langle B \rangle \langle R \rangle \, , \qquad \langle R' \rangle = (a+1)\langle R \rangle - b \langle R \rangle^2 - d \langle R \rangle \langle R \rangle \, . \tag{13}
$$

In addition to the trivial solution *R*<sup>0</sup> = *B*<sup>0</sup> = 0 corresponding to laminar flow, we have two kinds of non-trivial solutions, either single-sided (∗) with *R* = 0 and *B* = 0 or *B* = 0 and *R* = 0, the non-vanishing solution being *S*∗ = *a*/*b*, with *S* = *R* or *B*, or double-sided (∗∗) with *B*∗∗ = *R*∗∗ = *a*/(*b* + *d*). A straightforward stability analysis of the fixed points of iterations (13) shows that the one-sided solution is stable when *b* < *d* and unstable otherwise whereas the reversed situation holds for the two-sided solution. Returning to probabilities, the global stability threshold is thus given for *a* = 0; hence, (*p*<sup>1</sup> + *p*<sup>2</sup> + *p*5)cr = 1 and the one-sided solution is expected when *b* < *d*; i.e., *p*<sup>1</sup> + *p*<sup>2</sup> < *p*5. Results of the mean-field approach adapted from [14] to the present formulation will be illustrated in Figure 14 below.
