*3.2. Decay: 1D vs. 2D*

The model is designed to exemplify a decay according to the DP scenario in a two-dimensional setting, with specificities linked to the anisotropic propagation properties of the LTBs in transitional channel flow, and, in particular, propose an interpretation for the observation of 1D-DP exponents in the absence of transversal splitting (*p* <sup>4</sup> = 0). Accordingly, we examine the role of transverse diffusion (parameter *p*2) modeling the slight upstream shift that may affect LTBs at longitudinal splitting. We focus on a set of experiments with *p*<sup>5</sup> = 0.9, *p*<sup>2</sup> fixed, and control parameter *p*1. When *p*<sup>2</sup> cancels exactly, it is easily understood that transversal expansion is forbidden: An active *B* state at (*i*, *j* + 1) or *R* state at (*i* + 1, *j*) at time *t* cannot give birth to an active state of the same kind at (*i*, *j*) at *t* + 1. The evolution stems from processes associated with configuration C<sup>5</sup> with probability *p*<sup>5</sup> or C<sup>1</sup> with probability *p*1. These processes change occupancy only along direction *i* for *B* states, and *j* for *R* states, precisely in the direction corresponding to the single-sided regime considered (after termination of the transient). The dynamics are, therefore, strictly one-dimensional and decay is expected to follow the 1D-DP scenario. In contrast, introducing some transverse diffusion (*p*<sup>2</sup> = 0) immediately gives some 2D character to the dynamics. This is illustrated in Figures 9–12.

We consider first *p*<sup>2</sup> non-zero and relatively large *p*<sup>2</sup> = 0.1. Figure 9 displays the behavior of the turbulent fraction as a function of *p*1.

**Figure 9.** (**Left**): Time series of the turbulent fraction at different values of *p*1; average over 5 (10) independent simulations (*p*<sup>1</sup> <sup>=</sup> 0.0584 <sup>≈</sup> *<sup>p</sup><sup>c</sup>* <sup>1</sup>, black trace). (**Right**): Mean value of the turbulent fraciton at stationary state as a function of *p*<sup>1</sup> (original data). Inset: once raised to power 1/*β*, with *β* = 0.584 ≈ *β*DP for *D* = 2, the mean turbulent fraction tends to 0 linearly with an extrapolated threshold *p*<sup>c</sup> <sup>1</sup> = 0.05843.

The left panel illustrates the decrease of the turbulent fraction with the number of steps from a uniformly fully turbulent single-sided state (*F*<sup>t</sup> = 1 at *t* = 0) in a domain D = (192 × 192), showing the saturation to a finite value *F*t above threshold, a near power-law decay close to threshold, and an exponential decay below. The right panel presents the mean of *F*t after elimination of an appropriate transient as a function of *p*1, for simulations in domains up to 512 × 512 for the lowest values of *F*t. Once fitted in the range *<sup>p</sup>*<sup>1</sup> ∈ [0.058, 0.064] against the expected power law behavior *F*t = *<sup>a</sup>*(*p*<sup>1</sup> − *<sup>p</sup><sup>c</sup>* 1)*β* one gets *a* = 3.213 (2.936, 3.489), *p<sup>c</sup>* <sup>1</sup> = 0.05844 (0.05842, 0.05845), *β* = 0.5811 (0.566, 0.5962), in very good agreement with the value *β*DP ≈ 0.584 when *D* = 2 [2]. This is confirmed in the inset of Figure 9 (right) showing *F*t1/0.584 as a function of *<sup>p</sup>*<sup>2</sup> for *<sup>F</sup>*<sup>t</sup> small, the linear variation of which extrapolates to zero for *p*<sup>1</sup> ≈ 0.05843.

Having a good estimate of the threshold one can next consider the decay of the turbulent fraction, which is supposed to decrease as a power law at criticality, *p*<sup>1</sup> = *p*<sup>c</sup> <sup>1</sup>: *F*<sup>t</sup> ∼ *t* <sup>−</sup>*α*DP with *<sup>α</sup>*DP <sup>=</sup> *<sup>β</sup>*DP/*ν*DP where *ν*DP ≈ 1.295; hence, *α*DP ≈ 0.451 [2]. Figure 10 (left) shows that this is indeed the case for the compensated turbulent fraction *F*<sup>t</sup> × *t <sup>α</sup>*DP , up to the moment when fluctuations become too important due to size effects and lack of statistics. When *p*<sup>1</sup> is different from *p*<sup>c</sup> <sup>1</sup> but stays sufficiently close to it, the variation of the turbulent fraction keeps trace of the critical situation, except that the number of steps needs to be rescaled by the distance to threshold due to critical slowing down: the time scale *τ* diverging as (*p*<sup>1</sup> − *<sup>p</sup>*<sup>c</sup> 1) <sup>−</sup>*ν* , number of steps is rescaled upon multiplying it by (*p*<sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>c</sup> 1) *<sup>ν</sup>* . Figure 10 (right) indeed shows a good collapse of the compensated curves as a function of the rescaled number of steps when using the exponents corresponding to 2D-DP, *α*DP ≈ 0.451 and *ν*DP ≈ 1.295 [2].

**Figure 10.** (**Left**): Power law decay of the turbulent fraction at *<sup>p</sup>*<sup>1</sup> <sup>=</sup> 0.0584 <sup>≈</sup> *<sup>p</sup>*<sup>c</sup> <sup>1</sup> = 0.05843: compensation with *α*DP confirms the 2D nature of the process. (**Right**): Critical behavior near threshold: compensated turbulent fraction *F*t × *t <sup>α</sup>*DP as a function of the number of steps rescaled by (*p*<sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>c</sup> 1) *<sup>ν</sup>*DP for *<sup>p</sup>*<sup>1</sup> <sup>∈</sup> [4.5:0.2:6.5] <sup>×</sup>10−<sup>2</sup> surrounding the presumed critical value *<sup>p</sup>*<sup>c</sup> <sup>1</sup>, with the exponents corresponding to DP for *D* = 2.

We now consider *p*<sup>2</sup> = 0 which, as argued earlier, should fit the critical behavior of directed percolation when *D* = 1. In that case, when using square or nearly-square rectangular domains, size effects turn out to be particularly embarrassing as will be illustrated quantitatively soon. However, we can take advantage of the fact that, assuming propagation in the one-sided regime, e.g., along the direction for *B* active states, *NB* being the corresponding number of sites involved, the computed turbulent fraction is, in fact, the average of the activity over *NR* independent lines in the complementary direction, while still being sensitive to size effects controlled by *NB*. Accordingly, at given computational load (proportional to *NB* × *NR*), one can freely increase the size artificially in considering a strongly elongated domain D = [(*NB* × *k*) × (*NR*/*k*)], with *k* sufficiently large that the average over *NR*/*k* independent lines still make sense from a statistical point of view, while postponing size effects. With reference to a (192 × 192) domain, we have obtained good results with *k* = 16, i.e., 3072 × 12 up to *k* = 64, i.e., 12288 × 3.

Though this choice is a bit extreme, we present here results about 1D-DP criticality with the 12288 × 3 domain in Figure 11. The left panel displays the variation of the mean turbulent fraction with *<sup>p</sup>*1, which has been fitted against the expected power law, *F*t = *<sup>a</sup>*(*p*<sup>1</sup> − *<sup>p</sup><sup>c</sup>* <sup>1</sup>)*β*. One gets *<sup>a</sup>* = 1.473 (1.446, 1.5), *p<sup>c</sup>* <sup>1</sup> = 0.2682 (0.2682, 0.2683), *β* = 0.2701 (0.2664, 0.2738). This value of *β* is quite compatible with the value *<sup>β</sup>*DP <sup>≈</sup> 0.276 when *<sup>D</sup>* <sup>=</sup> 1 [2]. Furthermore, accepting this value, a linear fit of *F*t1/*<sup>β</sup>* with *p*<sup>1</sup> then provides an extrapolated threshold *p*<sup>c</sup> <sup>1</sup> = 0.26817. As seen in the right panel of Figure 11, in the neighborhood of *p*<sup>c</sup> <sup>1</sup> a good collapse is obtained for the compensated turbulent fraction as a function of the rescaled number of steps when using the exponents *<sup>α</sup>* = 0.159 and *<sup>ν</sup>* = 1.734 corresponding to 1D-DP [2].

**Figure 11.** (**Left**): Mean turbulent fraction at stationary state as a function of *p*1. Once raised at power 1/*β* with *β* = 0.276 ≈ *β*DP for *D* = 1, the mean turbulent fraction tends to 0 linearly with an extrapolated threshold *p*<sup>c</sup> <sup>1</sup> = 0.26817. (**Right**): Critical behavior near threshold: compensated turbulent fraction *F*t × *t <sup>α</sup>*DP as a function of the number of steps rescaled by (*p*<sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>c</sup> 1) *<sup>ν</sup>*DP for *<sup>p</sup>*<sup>1</sup> <sup>∈</sup> [0.260:0.002:0.280] surrounding the presumed critical value *<sup>p</sup>*<sup>c</sup> <sup>1</sup>, with the exponents corresponding to DP for *D* = 1.

Size effects already alluded to above are illustrated in Figure 12.

**Figure 12.** Size effects when *p*<sup>2</sup> = 0. (**Left**): Raw data showing late exponential decay and progressive prevalence of power-law decay as *NB* grows. (**Right**): Rescaled data. According to scaling theory, the appropriate scale for the number of steps (time *t*) is *NB <sup>z</sup>*/*D*; hence, *<sup>t</sup>* → *<sup>t</sup>*/*Nz*/*<sup>D</sup> <sup>B</sup>* , with *z* ≈ 1.58 when *D* = 1, while the turbulent fraction has to be compensated for decay as *F*<sup>t</sup> × *t <sup>α</sup>* with *α* = 0.159. The collapse of traces illustrates universality with respect to 1D-DP.

Displaying the turbulent fraction as a function of the number of steps for linear size *NB* from small systems to relatively large ones (*NB* = 64 up to 768) in lin-log scale, the left panel illustrates the late stage of decay right at criticality as obtained from the previous study summarized in Figure 11. It is seen that, in the time-window considered (0, 105) the exponential dependence observed at small sizes is progressively replaced by the power-law behavior expected at criticality at infinite size. Size effects are also ruled by scaling theory; see, e.g., [2] for DP. They relate to correlations in physical space that are associated with exponent *<sup>ν</sup>*⊥. The ratio *<sup>z</sup>* = *<sup>ν</sup>*/*ν*<sup>⊥</sup> is called the dynamical exponent and theory predicts that, for finite size systems, scaling functions depend on time with the number of sites as *tD*/*z*/*N* where *N* is the total number of sites. In the (quasi-)one-dimensional regime we are interested in, *D* = 1, *N* is just *NB* and *z* = 1.58 [2]. The right panel of Figure 12 indeed shows extremely good collapse of the traces corresponding to those in the left panel, once the number of steps is rescaled as

*t*/*N*1.58 *<sup>B</sup>* and the turbulent fraction is compensated for decay as *F*t(*t*) × *t* 0.159, both exponents taking on the 1D-DP values already mentioned.

Of interest in the context of channel flow decay, the crossover from 2D behavior for *p*<sup>2</sup> sizable (e.g., *p*<sup>2</sup> = 0.1, Figures 9 and 10) to 1D behavior for *p*<sup>2</sup> = 0 is of interest since *p*<sup>2</sup> is associated with the progressive importance of off-aligned longitudinal splitting as Re increases. A series of values of *p*2, decreasing to zero roughly exponentially, has been considered and the corresponding DP threshold has been determined as given in Table 1 and shown in Figure 13 (left).


**Table 1.** Values of *p*<sup>1</sup> at criticality at given *p*<sup>2</sup> (*p*<sup>5</sup> = 0.9 and *p* <sup>4</sup> = 0).

Except for *p*<sup>2</sup> = 0 determined as explained above (Figure 11), these values have been obtained in domains 192 × 192 with averaging over 10 independent experiments. Figure 13 (right) displays the averaged time-series of the turbulent fraction at criticality for each of these values of *p*2, once compensated for decay according to 2D-DP (*α* = 0.451).

**Figure 13.** Crossover *p*<sup>2</sup> → 0. (**Left**): Criticality condition separating the sustained active regime from the absorbing regime. (**Right**): A DP-like process governs the decay when the line is crossed, the characteristics of which can be understood from the asymptotic power-law decrease of the turbulent fraction as a function of time, here compensated by *t <sup>α</sup>*, with *α* = *α*DP(*D* = 2).

The results for *p*<sup>2</sup> = 0, evolving as *t <sup>α</sup>*2*D*−*α*1*<sup>D</sup>* are marked with (∗) and (∗∗) are obtained in the 192 × 192 domain and in the 3072 × 12 quasi-1D domain, respectively. In the time span considered here, the latter is free from finite-size effects which is not the case of the former with the corresponding compensated data decaying exponentially at the largest times. It is easily seen that, except for *p*<sup>2</sup> = 0, the compensated time series display a wide plateau indicating that 2D behavior holds for a certain amount of time. Whereas traces for *p*<sup>2</sup> = 0.1 and *p*<sup>2</sup> = 0.05 cannot be distinguished, for smaller values of *p*<sup>2</sup> the plateau regime starts at larger and larger times and develops after having followed the 1D trace for longer and longer durations, clearly indicating the influence of the anisotropy controlling the effective dimensional reduction. A similar consequence of the crossover affects the decrease of the mean turbulent fraction with the distance to threshold but, apart from this qualitative observation, no reliable information can be obtained on exponent *β* owing to the difficulty to reach the relevant critical regime.

We shall not document the case when *p*<sup>1</sup> = 0 and *p*<sup>2</sup> varies. This situation is not observed in the simulations since off-aligned longitudinal splitting is conspicuous only sufficiently above Reg, in the vicinity of which decay is fully accounted for by in-line longitudinal splitting modeled by a variable *p*<sup>1</sup> = 0, but the possibility remains, at least conceptually. The decay when *p*<sup>1</sup> = 0 happens to follow the same 1D-DP scenario though the argument is slightly less immediate. It relies on the observation that no growth is possible in the propagation direction of a given LTB species, whereas off-aligned longitudinal splitting (*p*<sup>2</sup> = 0) permits growth and diffusion in the transverse direction. Under the combined effects of transversal diffusion (*p*<sup>2</sup> small) and propagation (*p*<sup>5</sup> large), near-threshold, the sustained turbulent regime is made of quasi-1D clusters that are aligned with and drift along the diagonal of the lattice, i.e., the stream-wise direction, and get thinner and thinner when decaying, supporting the reduction to a "*D* = 1" scenario. Here, the trick used for *p*<sup>2</sup> = 0 does not work, and simulations in square domains are necessary with no escape for size effects which hinders the observation of the critical regime. Nevertheless, *p*<sup>c</sup> <sup>2</sup> when *<sup>p</sup>*<sup>1</sup> = 0 seems close to *<sup>p</sup>*<sup>c</sup> <sup>1</sup> when *p*<sup>2</sup> = 0, suggesting some symmetry between *p*<sup>1</sup> and *p*2.

The relevance of the results with *p* <sup>4</sup> ≡ 0 to transitional channel flow will be discussed in the concluding section. We now turn to the general two-sided case with transversal collisions and splittings.

#### **4. Beyond Onset of Transversal Splitting,** *P-* **<sup>4</sup>** *>* **0**

In statistical thermodynamics systems, critical properties at a second order phase transition leads to define a full set of exponents governing the variation of macroscopic observables close to criticality [26]. The concept of universality was introduced to support the observation that these systems can be classified according to the value of their exponents depending on a few qualitative characteristics, the most prominent ones being the symmetries of the order parameter and the dimension of physical space. This viewpoint can be extended to far-from-equilibrium systems such as coupled map lattices (CMLs) displaying nontrivial collective behavior. The associated ordering properties present many characteristics of thermodynamical critical phenomena at equilibrium. Universality classes beyond those known from equilibrium thermodynamics have been shown to exist with different sets of exponents. An additional criterion, the synchronous or asynchronous nature of the dynamics, has been found relevant to distinguish among them [27]. In the context of the present model, as soon as probability *p* <sup>4</sup> grows from zero, fully one-sided configurations previously reached after the termination of a possibly long transient are now unstable against the presence of states with the complementary color. The stationary regime that develops in the long term can be, either ordered, i.e., one-sided with one dominant active state (*B* or *R*), or disordered, i.e., two-sided with statistically equal fractions of each active state (*B* and *R*). Furthermore, a transition at some critical value *p* 4 <sup>c</sup> is expected to take place on general grounds. This gives us the motivation to study the response of the model to the variation of *p* <sup>4</sup> as a critical phenomenon beyond the mean-field expectations of Section 2.3.

The results of the mean-field approach, system (11) and (12), rephrased from [14], are depicted in Figure 14 (left). Upon variation of parameter *c* representing *p* <sup>4</sup> up to an unknown rescaling factor, all along the one-sided regime (*c* < *c*cr), the total turbulent fraction is seen to decrease while the order parameter measuring the lack of symmetry similarly decreases to zero according to the usual Landau square-root law. Obviously symmetrical, the two-sided regime (*c* > *c*cr) is then characterized by a total turbulent fraction that regularly grows due to the contribution of splitting, whatever the type of active state.

From now on, we shall simply refer to the turbulent fractions and other statistical quantities as their time average over a sufficiently long duration, up to 2 × <sup>10</sup><sup>6</sup> simulation steps, after elimination of an appropriate transient, up to 105 steps, the largest values being necessary close to the transition point owing to the well-known critical slowing down. On the one hand, the total turbulent fraction is obviously defined as *F*t(*B*) + *F*t(*R*), where the over-bar denotes the time averaging operation. (Later on, we shall omit this over-bar when no ambiguity arises between the instantaneous value of a quantity and its time average, especially for the axis labelling in figures.) On the other hand, the lack of symmetry can be measured by the signed difference averaged over time *F*t(*B*) − *F*t(*R*), able to distinguish global *B* orientation from its *R* counterpart, or rather its absolute value *F*t(*B*) − *F*t(*R*) since we are only interested in the amplitude of the asymmetry (called '*A*' in [14]) and not in which

orientation is dominant, the two being equivalent a priori for symmetry reasons. However, due to the finite size of the system, in the symmetry-broken regime close to threshold, orientation reversals can be observed as illustrated later (Figure 15), so that blind statistics in the very long durations are no longer representative of the actual ordering. Like in thermal systems [28] or their non-equilibrium counterparts [27], it is thus preferable to define the order parameter through the mean of the unsigned difference: |*F*t(*B*) − *F*t(*R*)| . Corresponding simulation results are displayed in Figure 14 (right) for a system of size (256 × 256). The general agreement between the two diagrams is remarkable, up to an unknown multiplicative factor translating *c* into *p* <sup>4</sup>, as discussed earlier. One can notice that the order parameter is minimal but not zero in the two-sided regime, which is due to fluctuations and the fact that the two operations of averaging over time and taking the absolute value do not commute. Finite-size effects are also apparent as a rounding of the graph at the location of the would-be critical point in the thermodynamic limit.

**Figure 14.** (**Left**): Bifurcation diagram of system (11) and (12) after [14]. The total turbulent fraction is *F*t(*B*) + *F*t(*R*) and the order parameter characterizing the transition is abs(*F*t(*B*) − *F*t(*R*)). A standard supercritical bifurcation is expected for this quantity with abs(*F*t(*B*) <sup>−</sup> *<sup>F</sup>*t(*R*)) <sup>∝</sup> (*c*<sup>c</sup> <sup>−</sup> *<sup>c</sup>*)1/2 in the one-sided regime, whereas *F*t(*B*) = *F*t(*R*) in the two-sided regime. (**Right**): Time average of turbulent fractions as functions of the control parameter *p* <sup>4</sup> after elimination of an appropriate transient as obtained from simulations of the stochastic model.

The current justification for taking the absolute value is that the time between orientation reversals diverges with the system size and the phase transition only takes place once we have taken the thermodynamic limit of infinitely large systems studied over asymptotically long durations [28]. Accordingly, very long well-oriented intermissions can be considered as representative of the symmetry-broken regime. The problem is illustrated in Figure 15 displaying the time series of *F*t(*B*) − *F*t(*R*) and histograms of |*F*t(*B*) − *F*t(*R*)| for *p* <sup>4</sup> = 0.0121, still in the one-sided but already alternating regime, next for *p* <sup>4</sup> = 0.0125 and 0.0126, where one can notice a change in the shape of the histogram, and finally for *p* <sup>4</sup> = 0.0140, sufficiently deep inside the two-sided regime where the histogram displays a sharp maximum at the origin. On this basis one could use the histograms of the "order parameter" and determine the threshold from the position of its most probable value, whether non-zero in the symmetry-broken state or at the origin when symmetry is restored. This procedure would give *p* 4 <sup>c</sup> <sup>≈</sup> 0.01255.

The symmetry-breaking bifurcation can now be studied beyond the mean-field description as other collective phenomena studied in equilibrium and far-from-equilibrium statistical physics: In addition to the order parameter, the variation of which leads to the definition exponent *β* in the ordered regime, another observable of interest is the susceptibility measuring the response to an applied field conjugate to the order parameter, vis. *M* = *χH* with the magnetization *M* coupled to magnetic field *H* in the case of magnets. The susceptibility diverges near the critical point, with leads to the definition of two exponents *γ* and *γ* in the disordered and ordered regime, respectively. Universality implies *γ* = *γ* , as can already easily be derived in the mean-field framework. When a conjugate field cannot be defined, one uses the property that fluctuations take the instantaneous value of the order parameter away from its average value, which can be understood as resulting from the response to a conjugate field. This helps one to relate the susceptibility to the variance of fluctuations of the order parameter. The identification is up to a multiplication by the "volume" of the system that has to be introduced in order to compare the results from systems with different sizes. This is what will be done here; hence, *χ* = *NBNR* × var(|*F*t(*B*) − *F*t(*R*)|). As shown in Figure 16 (top), this quantity displays a sharp maximum, indicative of the singularity expected at the thermodynamic limit. In a finite-size but large system, the critical point is then estimated from the position of the maximum of the susceptibility. Here, this gives *p* 4 <sup>c</sup> <sup>≈</sup> 0.0123 slightly smaller but compatible with the value obtained above from the examination of the histograms. Unfortunately, this discrepancy due to size-effects forbids us to determine exponents *β* and *γ* with some confidence.

**Figure 15.** Time series of the instantaneous mean orientation as measured by *F*t(*B*) − *F*t(*R*) in a 256 × 256 domain for *p*<sup>1</sup> = *p*<sup>2</sup> = 0.1 and *p*<sup>5</sup> = 0.9, at *p* <sup>4</sup> = 0.0121 below the onset of the one-sided regime, at *p* <sup>4</sup> = 0.0125 and 0.0126 near the bifurcation point, and at *p* <sup>4</sup> = 0.0140 in the two-sided regime. Bottom line: Corresponding histograms of |*Ft*(*B*) − *Ft*(*R*)|. The histograms were all built using 75 bins and contain the same number of points for 10<sup>5</sup> <sup>&</sup>lt; *<sup>t</sup>* <sup>&</sup>lt; <sup>2</sup> <sup>×</sup> <sup>10</sup>6, but the vertical scales are not identical.

Having in mind results of the mean-field approach, namely, *β* = 1/2 and *γ* = 1, we can, however, estimate the range where stochastic fluctuations have nontrivial effects. The bottom-left panel of Figure 16 displays the variation of the order parameter already shown in Figure 14 (right), but now squared in order to show that, far from the critical point, the system fulfils the mean-field square-root prediction to an excellent approximation, with an extrapolated threshold *p* 4 MF 0.0133, shifted upwards with respect to the estimates obtained from the simulations *p* 4 c. 0.0123–0.0125.

In the same way, the divergence of the susceptibility with exponents *γ* = *γ* = 1 expected from the mean-field argument shows up upon retreating the data already given in Figure 16 (top) and plotting 1/*χ* as a function of *p* <sup>4</sup>. This is done in Figure 16 (bottom-right) showing the same linear variation of 1/*χ* below and above the transition point, in agreement with the theory. The extrapolation of the linear fits on both sides of the transition yield *p* 4 <sup>c</sup> <sup>≈</sup> 0.0127 in reasonable agreement with the value obtained from the order parameter variation in the same conditions and definitely larger than the empirical values. Clearly, deviations seen in the boxed parts of these two figures warrant further scrutiny, motivating our current approach via finite-size scaling theory [28] in search for universality. On going work attempts at a full characterization of the critical regime through exponents determination. Though local agents do not behave as Ising spins, symmetries are basically identical, so that the equilibrium 2D Ising universality class or its non-equilibrium extension [27] might be relevant. We shall discuss this further below.

**Figure 16.** (**Top**): Variation of the susceptibility *χ* as a function of *p* <sup>4</sup>. (**Bottom**): Evidence of mean-field behavior away from the critical regime: The order parameter squared (**left**) and the inverse of the susceptibility (**right**) both vary linearly with *p* <sup>4</sup> sufficiently far from the mean-field extrapolated threshold.
