{'fraction_non_alphanumeric': 0.03926393124480177, 'fraction_numerical': 0.01098558358746881, 'mean_word_length': 4.500762485703393, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatio-temporal chaos of one-dimensional thin elastic layer with the rate-and-state friction law

Independent of specific local features, global spatio-temporal structures in diverse phenomena around bifurcation points are described by the complex Ginzburg-Landau equation (CGLE) derived using the reductive perturbation method, which includes prediction of spatio-temporal chaos. The generality in the CGLE scheme includes oscillatory instability in slip behavior between stable and unstable regimes. Such slip transitions accompanying spatio-temporal chaos is expected for frictional interfaces of a thin elastic layer made of soft solids, such as rubber or gel, where especially chaotic behavior may be easily discovered due to their compliance. Slow earthquakes observed in the aseismic-to-seismogenic transition zone along a subducting plate are also potential candidates. This article focuses on the common properties of slip oscillatory instability from the viewpoint of a CGLE approach by introducing a drastically simplified model of an elastic body with a thin layer, whose local expression in space and time allows us to employ conventional reduction methods. Special attention is paid to incorporate a rate-and-state friction law supported by microscopic mechanisms beyond the Coulomb friction law. We discuss similarities and discrepancies in the oscillatory instability observed or predicted in soft matter or a slow earthquake.


I. INTRODUCTION
Earthquakes are fast shear slips accompanied by rock fracturing, and they frequently occur at the upper boundary of subducting tectonic plates. In addition to these regular earthquakes, other types of shear slip events, called slow earthquakes, have been discovered [1]. Slow earthquakes have a much longer duration and smaller radiation of seismic energy than what would be expected from their seismic moment. This suggests that they obey scaling laws different from those for regular earthquakes [2]. Moreover, the cumulative distribution of the seismic energy rate is close to an exponential distribution [3][4][5], which is quite different from the power-law distribution expected from the Gutenberg-Richter law [6].
Similar slow slips are repeated many times quasi-periodically [7]. Because slow earthquakes release stress stored around a plate boundary, a simple description of their dynamics will be helpful to understand the spatio-temporal dynamics of slip events in subduction zones.
An indication of how to make a simple description is the fact that slow earthquakes occur in the transition zone between a stably sliding aseismic zone and a seismic zone hosting regular earthquakes [8] [see Fig. 1 (a),(b)]. A shift in the sliding behavior along the subducting tectonic plate invokes a bifurcation, where stable sliding become unstable [9,10]. According to weakly nonlinear expansions near the bifurcation point, the equation of motion for each point is reduced to the Stuart-Landau equation as small-amplitude oscillators with a cubic polynomial [11,12]. Furthermore, the local in-field coupling incorporates a diffusion term with a complex coefficient representing the interaction between adjacent oscillators into a governing equation, so that the reduced equation turns into the complex Ginzbrug-Landau equation (CGLE) [11,12]. This equation is simple, but it captures a rich variety of phenomena near a bifurcation or transition point, such as transitions found in superconductivity, liquid crystals, Bernard convection, and Taylor vortices [11,13]. Given the broad range of applicability of the CGLE, a natural question here arises: does the CGLE describe some features of slow earthquakes? If yes, a viewpoint from the CGLE may contribute to progress in the study of the repeating nature of slow earthquakes and their predictability.
In this paper, we investigate slow earthquake analytically and numerically by focusing on an analog between the transition zone and the bifurcation. The analyses are based on the reductive perturbation method [11], which simplifies the original governing equation. This simplification is more than an approximation because it extracts the universal features that do not rely on specific properties. In sec. II, we first introduce a fault system modeled by a thin elastic layer under the rate-and-state friction (RSF) law. In sec. III, linear and weakly non-linear analyses are applied near the bifurcation point to show that the slow and global structures are represented by the CGLE, which may display spatio-temporal chaos due to Benjamin-Feir (BF) instability. Section IV demonstrates a numerical simulation that obeys the original governing equation before the reductive perturbation method, and discusses the equation and method together with the analytical consequences. An interesting point is that the size distribution of slip events numerically obtained for the chaotic regime is exponential, invoking the cumulative distribution of the seismic energy rate reported for slow earthquakes. In sec. V, we argue that applying the CGLE to slow earthquakes is valid.
Section VI concludes this study.

II. THIN LAYER MODEL WITH THE RSF LAW
We begin by combining the original governing equation and the RSF law. A cross section of the subduction zone is illustrated in Fig. 1 (a), where the transition zone is between the seismogenic and stable sliding zones. Figure 1 (b) shows the interface between the plates from a different angle. For simplicity, we consider it as an antiplane fault (mode III), where the displacement is perpendicular to thex-axis. In addition, although the plates bend on a large scale, the upper plate is modeled as an infinite planar layer, with the coordinate system embedded as shown in Fig. 1 (c). The layer is sheared and displaced along theŷ-axis. Note that displacements denoted byû(x,t) are uniform for theŷ-component [see Fig. 1 (d)], and thus, we can deal with the system as a one-dimensional thin elastic layer spanning thexdirection. The bottom boundary is fixed in the coordinate system, and the top boundary is dragged with constant loading velocityv l along theŷ-axis, as shown in Fig. 1 (d), while the lower plate has an relative movement in the down-dip direction [ Fig. 1 (b)]. When the thickness of the elastic layerĥ is small enough, imposing the boundary condition on Navier's equation simplifies the equation of motion for the thin layer [14] aŝ wheret,ρ,Ĝ, µ,σ, andv ≡ dû/dt denote the time, mass density, shear modulus for the elastic body, friction coefficient, normal stress, and velocity along theŷ-axis, respectively.
On the right side, the first two terms of the elastic force are derived from the continuum body model of the thin layer, where the interaction range is finite and on the order of the thickness of ≃ĥ. The restoring force for the global inhomogeneous deformation along thexaxis is replaced with the operator (∂ 2 /∂x 2 ), which serves as the local coupling. In addition, the top boundary condition emerges as loading in the second term. The last term is the friction that obeys the RSF law [15][16][17], which acts on the bottom boundary interface with magnitude µ(x,t)σ. The RSF law is commonly exploited to describe rock friction [18,19], and is also known to be readily applicable to other materials, as seen in the literature [20]. In addition, it is often used to describe regular [21] and slow [22] earthquakes. The constitutive law for the RSF is written as wherev * ≡D c /θ * , andθ denotes a variable that represents the interface state history. As an evolution law, this paper adopts the slip law [10,23], because it provides a better description of a fairly symmetric response with the characteristic distanceD c to a discontinuous increase or decrease in stress [24]. The RSF law gives the steady-state friction as meaning that velocity weakening occurs for a − b < 0 while velocity strengthening occurs for a − b > 0. A Hopf bifurcation is observed below around a − b ≃ 0 and a − b < 0, which are the main regions of interest for slow earthquakes.
In preparation for the reductive perturbation method, we recast eqs. where with three dimensionless variables: Note that, for compact notation, the units for length, time, and stress are introduced, respectively, as follows:D We use the dimensionless variables t =t/t u , In addition, parameters associated with the physical properties are listed below: where G represents the coupling with the driving factor provided by the upper plate, β is the distance from the velocity weakening point (a − b = 0), and c is the sound velocity in the idealized plate. Bear in mind that the first two parameters identify the sliding stability, although the last one turns out to be irrelevant. In the following discussion, we restrict ourselves to the region where ξ 2 + v l = v remains positive, and thus, sgn(ξ 2 + v l ) = 1. In the positive condition, we may circumvent the nonanalytic point ξ 2 + v l = v = 0 that might cause a numerical instability (see the appendix).

A. Hopf Bifurcation
We first determine a fixed point through a linear stability analysis without spatial coupling by setting c 2 ∂ 2 ξ 1 /∂x 2 = 0 in Eq. (5). The steady sliding is given by the solution ξ s = t (ξ 1s , ξ 2s , ξ 3s ) to the stationary state dξ i /dt = F i (ξ) = 0: We look at small deviations δξ ≡ ξ − ξ s around the steady-state solution to linearize Eqs.
The eigenvalues λ of the matrix L satisfy the characteristic equation det (L − λI) = 0. Under the assumption that the roots of the cubic equation are one real root α and two complex roots ǫ ± iω, the real root is obtained as α = −Gv l /(ǫ 2 + ω 2 ) < 0. This ensures that the directions irrelevant to the Hopf bifurcation are stable (see the appendix). In contrast, the bifurcation is signified by ǫ. At v l = v c , when ǫ turns from a negative to a positive value, a supercritical Hopf bifurcation appears with the following velocity and angular frequency: which agrees with previous results [10]. It is worth noting that only two crucial parameters β and G enter into the equation to determine the stability behavior. The parameter β is the distance from the velocity-weakening point, while G represents the effective spring constant.
For Eq. (14) to represent real values, G − β > 0 and β > 0, or 0 < β < G, is required. As long as this condition persists, the Hopf bifurcation point v c = (G − β)/β always exists, and the non-oscillatory solution for v l > v c is unstable.

B. CGLE
We move on to a weakly nonlinear analysis of Eqs. where The matrix L is expanded in a series of ε together with the 0-th order L 0 given by replac- ing v l with v c in the matrix of Eq. (13). Here, we employ notation such as (N 0 uvw) i = (N 0 ) ijkl u j v k w l , with summations running over repeated indices. In accordance with expansions in the matrices, the eigenvalue λ is replaced with A supercritical Hopf bifurcation is observed for Re[λ 1 ] > 0, which Eq. (55) satisfies because we now just focus on v l > v c ⇔ ǫ > 0. In addition, the diffusion (local coupling) term in the present system is specified by (D) ij = c 2 δ i2 δ j2 .
Taking the term up to the ε 3/2 order with an appropriate procedure [11], we arrive at the CGLE with the slow time variable τ = εt: Note that the complex-valued W (τ, x) is associated with the original governing Eq. (5) through where the right eigenvector for L 0 with eigenvalue λ 0 =iω c is given by and the coefficients with complex values are The time scale between t and τ is separated enough for the slow variable W (τ, x) to undergo an adiabatic change. Equation (90) [13], where conservative nonlinear wave phenomena are found without a dissipation mechanism. The present system is between these regimes and identified with a combination of the complex coefficients g and d.
Before including spatial coupling, let us discuss the general role of the complex coefficient g of the cubic term through a homogeneous example with d = 0, or equivalently, the one reduced to the Stuart-Landau equation. Qualitative features are grasped in the amplitudephase representation by rephrasing Eq. (90) with W = Re iΘ : where amplitude R and phase Θ are real numbers. Eq. (23) gives the steady solution , which implies that the solution to original Eq. (5) oscillates with a finite amplitude. The cubic term in Eq. (90) provides the solution not obtained from the linear analysis with Eq. (52). In addition, a close inspection shows that a leading part of the angular frequency is determined by ω c , but modified by dΘ/dτ [Eq. (24)]. Notably, the oscillation amplitude R affects the angular frequency through a higher-order perturbation term.

C. Benjamin-Feir Instability
When (ii) spatial coupling participates in Eq. (90), spatio-temporal chaos, referred to as the BF instability, emerges for where the spatial gradient of the phase distribution is enhanced. The BF instability appears when the gradient in the phase modifies the amplitude R, and it can be intensified through Eq. (24). Indeed, in the original Eq. (5), we find that the coupling term c 2 ∂ 2 ξ 1 /∂x 2 appears in the time derivative of ξ 2 , not in ξ 1 . The coupling term does not smooth the phase but shrinks the amplitude when the phase advances more than those of neighboring locations; consequently, the BF instability may appear. From the positive nature of β, G, and G − β, the condition for the BF instability is reduced to which separates the BF-stable and BF-unstable regions in the G − β/G plane, as shown in Fig. 2 (b). We also use Eq. (26) to distinguish the BF stability from BF instability by coloring the boundary curves with red and black, respectively, on the phase diagram in Fig. 3.

IV. NUMERICAL RESULTS
We performed a numerical simulation using the original governing Eqs. (5)-(8) to verify the reductive perturbation predictions. Note that an auxiliary viscous friction linearly proportional to the relative velocity between adjacent segments ν∂ 2 ξ 2 /∂x 2 was incorporated as the second component of Eq. (5) to suppress the numerical instability. Indeed, the supplemental material [25] provides positive evidence that the appropriate viscous friction was not strong enough to alter the qualitative behavior.
An explicit scheme under a periodic boundary condition was employed. Bear in mind that, as mentioned earlier, the simulation was halted when the velocity took on an invalid Fig. 3 (a)]. The simulation started at t = 0, when the initial condition was prepared from the homogeneous steady solution of Eqs. (6)- (8) with the addition of small spatio-temporal disturbances. Because a strong artifact from the initial random conditions is unfavorable, we waited for the system to settle, and then the numerical samples were collected (see the appendix for detailed conditions).

A. Phase Diagram
Typical phase diagrams are given in Fig. 3 (a), (b), where Hopf-bifurcation points sweep along solid curves determined with Eq. (14). In fact, persistent oscillation was not seen below these curves. The black and red curves signify whether the BF instability is absent  or present, respectively.
Yellow squares around the black curve in Fig. 3 (a) indicate regular oscillations, whose dynamic behavior is exemplified in Fig. 4 (a-i). The pattern that appeared was not propagation but rather a standing wave. Indeed, as in Fig. 4 (a-i), it looked like the wave propagated faster than the sound speed indicated by the white triangular area at the lower left corner.
Incidentally, the translational symmetry of the pattern was broken as a consequence of the initial randomness. Additionally, Fig. 4 (a-i) shows the velocity cross section taken from the density diagram along the red broken dashed line. The velocity profile looks symmetric about the loading velocity (v l = 2). More quantitative analyses are shown in Fig. 4 (a-iii) (bottom).
The duration of time spent around the maxima or minima of v was the longest, which means that p(v) has peaks near v = 0 or v = 4.
In contrast, irregular oscillations were observed (green crosses) near the red curves in Fig. 4 (b). The typical behavior of irregular oscillation has instability at short wavelengths, as in Fig. 4 (b-i), (b-ii). The wave pattern propagated at the velocity indicated by the white triangular area, which corresponds to the sound speed. This irregular oscillation meets with an example of spatio-temporal chaos due to the BF instability, because the phase distributions have random spikes in Fig. 3 (b-ii). A velocity profile was also taken from the spatio-temporal density map along the red dashed line. Although the oscillations are centered at v = 2, amplitudes became smaller than those of Fig. 3 (i). This point is made clear by looking at p(v) in Fig. 3 (iii), where bimodal distributions are found. In addition, the duration time spent around the lowest velocity was longest for this condition. Thus, a peak around v ≃ 1 is found, and p(v) turns out to be asymmetric.
The other phase is a slip pulse discovered at the points marked by blue asterisks in Fig. 3 (b). In this parameter region, spatially localized domains with finite slip velocity propagated at the sound speed, as those have a slope comparable to that of the white triangular area at the lower left corner. The slip velocity distributions look like a power law, where a large number of events on the distribution was found at small v. The qualitative condition for pulse occurrence is consistent with that in previous studies, where the slip pulse was reported for the model with velocity weakening friction [24,26].
Focusing attention on the irregular oscillation associated with the BF instability among these four types of characteristic dynamics, we look for common features shared with slow earthquakes, such as tectonic tremors, or slow slip events. One of the noteworthy points is the event size distribution. The simulation was performed in a larger space and with longer duration to collect statistical data [see Fig. 5 (a) and the appendix for detailed conditions].
The slip behavior is exemplified in Fig. 5 (a), where an irregular oscillation is observed. is supposed to be associated with the contour fractal curves characterized by the Hurst exponent H through ζ N = 1 − H/2 [27]. In particular, this was discussed for the fractal contours drawn with fractional Brownian motion [28]. This is not, however, true in the present system. As a matter of fact, neither the appearance of Fig. 5 (c) nor the distribution in Fig. 5 (d) seem to provide exponents H and ζ N .
The probability density function p(e s ), with e s estimated from binarized plots like those in Fig. 5 (c), is shown in Fig. 5 (e). As v th changes, the exponential decreases along the event size e s are maintained. We are also aware that it is unlike the power law. Such exponential distributions invoke the cumulative distributions of seismic energy rates observed in tectonic tremors [4]. The exponential distributions are meant to have their own characteristic sizes.
Fitting the exponential decay p(e s ) ∼ exp (−e s /e s0 ) on Fig. 5 (e), we estimated the characteristic event sizes e s0 , which are plotted in Fig. 5 (f). In addition, Fig. 5 (e) is rescaled with e s0 , as in Fig. 5 (g).

V. DISCUSSION
Numerical and analytical aspects are discussed in secs. III and IV. Let us first verify the agreement between them from the viewpoint of qualitative and quantitative consistencies.
The qualitative point is the distribution shape. Recall that the numerical event size displays the exponential distribution, where the temporal plots on the spatio-temporal plane are fairly periodic [see the right side profile of fig. 4 (b-i)], but the spatial plots are disordered [see fig. 4 (b-ii)]. To make the point clear, we address both the temporal and spatial points. Although the temporal period varies over the long term with spatio-temporal chaos, a periodic pattern is found, implying that the characteristic time may correspond to the temporal period T c . In contrast, the spatial profiles are disordered rather than periodic because the neighboring differences in phase are intensified due to the BF instability mechanism. Especially, the spatial correlations get lost in the finite distance Λ c , suggesting that the correlations decay in an exponential manner. Therefore, the characteristic time and length may be defined using the temporal period and correlation length, respectively. Because the event size is defined as its area (time × length) on the spatio-temporal plane, it is a natural consequence if an event has the characteristic event size e sc ∼ T c Λ c , for which the size distribution is exponential.
Let us next review the quantitative applicability. As shown in sec. III, the characteristic time should be comparable to T c = 2π/ω c = 2π/ √ G − β. The characteristic length Λ c is estimated from the wavelength of the most unstable mode on the propagating wave solution (see the appendix): So far, we have seen a correspondence between the analyses and the numerical simulations. We now move on to the main topic -slow earthquakes -in light of the CGLE. In recent seismic observations [4], exponential distributions have been reported, and the GGLE around the BF instability may qualitatively reproduce this exponential feature. However, a quantitative agreement has not been achieved. Indeed, when plausible parameters [21,[29][30][31][32] are plugged into the CGLE, we encounter unrealistically huge characteristic length and time scales ofΛ c ≃ 10 19 m andT c ≃ 10 7 s (see table I in the appendix). We should not, however, rush to the conclusion that the CGLE is not completely appropriate in the study of slow earthquakes, because the CGLE is a very general approach independent of the specific structure of the system. One possible reason for the unrealistic result is the inhomogeneity of the system. The spatial heterogeneity of frictional properties in slow earthquake zone is suggested from the diversity of tectonic tremor activity [4]. Another reason is the boundary condition, especially along the dip direction, which is not taken into account in the present model. Indeed, the shallower part is strongly locked while the deeper part is steadily sliding.
Further, additional variables may play relevant roles. Recent observations reveal that fluid migration and precipitation of dissolved chemicals may control slow earthquakes [33]. Such fluid migration is modeled to reproduce slow as well as regular earthquake [34].
To proceed further with the CGLE, analog experiments are inevitable because it may be feasible to check the CGLE validity by preparing a homogeneous system. If the CGLE is applied to analog laboratory experiments with silicone gels [35], the estimated values for λ c ≃ 10 0 m andT c ≃ 10 −1 s, which is accessible range. The corresponding phenomena have not been seen yet at the laboratory scale, but it is a necessary step to explore slow earthquakes with the CGLE.

VI. CONCLUDING REMARKS
We have discussed slow earthquakes in terms of the CGLE. The location of slow earthquakes, between stable sliding and seismogenic zones, naturally leads us to a bifurcation analysis. The CGLE has a long history of being employed as a successful framework close to the bifurcation (transition) points for various phenomena. Our study is motivated by this simple idea.
One side of a slow earthquake may be captured by the CGLE, but the other side may not.
The cumulative distribution of the seismic energy rate reported for slow earthquakes is rather exponential, which is reminiscent of the numerical results for the size distribution based on the original governing equation. The exponential feature is also a natural conclusion deduced from the CGLE on the BF instability. However, the order estimate of the spatial size is far much larger than the actual size. The procedure to derive the CGLE from the governing equation is quite general. Thus, we should go back to the original governing equation to look for a reason for the discrepancy in the order estimate. The local coupling expressed by the second spatial derivative might be too rough, because the inhomogeneity of the plate in the subduction zone is totally disregarded, although the thin elastic sheet model is justified from the analyses under homogeneous rigidity, as discussed elsewhere by one of the authors [36]. Moreover, even if the governing equation is modified, the conventional reductive perturbation method to derive the CGLE seems not to be applicable to complicated inhomogeneous systems in its current form.
It would be interesting if efforts to use the concept to extract the universal features around the bifurcation discover that slow earthquakes share characteristics in common with other scalable phenomena. To proceed further along these lines, an elaborate investigation is clearly needed to see whether the CGLE concept is workable. It should be noted that, when ξ 2 + v l = v = 0, Eqs. (7) and (8) are not analytic. We tried several modifications of Eqs. (7) and (8) to avoid the non-analytic behavior near v = 0, especially to stabilize the numerical calculation. However, violent oscillation appeared, even after modifications, when v changed its sign. Resolving this situation is beyond the scope of this study. Thus, we analytically and numerically restricted ourselves to ξ 2 + v l = v > 0.
If ξ 2 + v l = v became 0 or negative, the numerical simulation was halted.
The parameter region β > 0 and β > G finds the critical velocity and angular frequency, respectively, in Eq. (14), as in the main text.
In contrast, if β < 0, oscillation does not appear in any sliding velocity v l . Also, if β > G, oscillation emerges for any positive sliding velocity v l . The reductive perturbation method cannot be applied near the bifurcation point in either case, but they are not of interest here.

C. Numerical Setup
An explicit scheme under a periodic boundary condition was employed. The system was discretized with sizes ∆t = 0.0002 and ∆x = 0.05. In addition, the dimensionless parameters Each simulation started at t = 0 from a random initial condition. Because a strong artifact due to the initial condition is unfavorable, we waited for the system to settle. In and The solution exists for k < k max or λ > λ min with k max = 2π/λ min , where Moreover, for the propagating wave to be stable against infinitesimal perturbations, we require an additional condition: with the most unstable wavelength given by Eq. (27). Eventually, we arrive at the analytical estimate for e sc = T c Λ c with T c = 2π/ω c , Eq. (14), and Eq. (35).
When the value used for the numerical simulation in Fig. 5 is inserted, T c = 4π, Λ c = 3π, and e s0 = 12π 2 ≃ 100.

E. Order estimation of parameters
We estimated the order ofΛ c andT c for geologically plausible conditions, as well as laboratory experiments using soft materials. By referring to the parameters relevant to actual earthquakes [21,[29][30][31][32] or laboratory experiments with gels [35,37,38], we estimated the order of parameters listed in Table I using In addition, we applied the following assumptions for (i) the tectonic plate and (ii) the gels: (i) The thicknessĥ was chosen to satisfŷ  The derivation is essentially follows [11].
We first consider weakly nonlinear analysis of given equation without spatial coupling.
Starting from following derivatives: and the steady state solutions: the following expressions are obtained at the steady state: Loading velocity v l is set close to Hopf bifurcation point as v l = v c (1 + ε). Expanding Eq. (43)-(46), by noting up to second order of ε, we have The fluctuation around the steady state, δξ obeys where We define λ satisfy λ = λ 0 + ελ 1 + ε 2 λ 2 , where Ξ and Ξ * are right and left eigenvector for L 0 .
Explicit description of L is and ω c = √ G − β. Apparently, α c < 0 indicating the corresponding direction in phase space is stable to perturbation. By contrast, ±iω c have 0 real part. Indeed, the eigenvalue space for ±iω c , become unstable at the Hopf bifurcation point.
The right eigenvector Ξ corresponding for iω c is , while the left eigenvector Ξ * for iω c is where p is a normalization factor to satisfy Ξ * Ξ = 1 and Using these expressions, one can obtain the leading term due to the perturbation ε, Following results are relevant in the latter discussion, and Re [λ 1 ] > 0 indicates that the steady state solution is not stable when ε > 0. Furthermore, Re [λ 1 ] represents the growth rate of oscillation amplitude.
The higher order contribution of ε is necessary to converge the amplitude of destablized oscillation, and is included in the coefficient g: where For the calculation, one should use the following matrices to calculate g: /G E and For ε=0, ∂ 2 F 2 ∂ξ 2 and From G > 0, β > 0 and G − β > 0, the condition for Benjamin-Feir instability is or, we may rewrite