{'fraction_non_alphanumeric': 0.04736951983298539, 'fraction_numerical': 0.016263048016701462, 'mean_word_length': 4.109985065073608, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Emergence and seismological implications of phase transition and universality in a system with interaction between thermal pressurization and dilatancy

A dynamic earthquake source process is modeled by assuming interaction among frictional heat, fluid pressure, and inelastic porosity. In particular, fluid pressure increase due to frictional heating (thermal pressurization effect) and fluid pressure decrease due to inelastic porosity increase (dilatancy effect) play important roles in this process. Two nullclines become exactly the same in the system of governing equations, which generates non-isolated fixed points in the phase space. These lead to a type of phase transition, which produces a universality described by the power law between the initial value of one variable and the final value of the other variable. The universal critical exponent is found to be 1/2, which is independent of the details of the porosity evolution law. We can regard the dynamic earthquake slip process as a phase transition by considering the final porosity or slip as the order parameter. Physical prediction of phase emergence is difficult because the porosity evolution law has uncertainties, and the final slip amount is difficult to predict because of the universality. Finally, nonlinear mathematical application of the result is also discussed.


I. INTRODUCTION
Dynamic earthquake slip processes are frictional slips inside the Earth. They show various qualitatively different behaviors, and discontinuous transitions between them are widely observed. For example, several earthquakes show pulse-like slips, whereas some generate crack-like slips as their slip profiles [1]. The existence of slow earthquakes with slip velocity and rupture velocity (the propagation velocity of the fault tip) much smaller than those for ordinary earthquakes is widely known [2]. There are also a variety of stress drops, which is defined as the difference between the shear stress acting on the fault plane and the residual shear stress. Almost all earthquakes have scale-independent stress drops; however, some are considered to have enormously large ones [3]. Whether intermediate cases exist or not is a controversial problem. If we refer to each aspect as a phase, these behaviors may be understood in terms of phase transitions in a unified manner.
To explain such transitions, it is insufficient to assume cracks in classical elastic bodies. Several researchers have noticed several aspects of the interior of the Earth, such as fault rock melting [4] and chemical effects [5]. For example, the slow slip velocity and rupture velocity for slow earthquakes are impossible to reproduce with the classical crack model; generation of slow earthquakes is considered to be promoted by certain factors, e.g., migration of the fluid (water) [6,7]. Some studies have treated interactions among these effects. For example, frictional heating and fluid pressure are considered to interact as thermal pressurization [8]. This interaction describes elevation of fluid pressure due to frictional heating. Such fluid pressure elevation induces a decrease in the normal stress acting on the fault plane, leading to frictional stress decrease and slip acceleration. Many studies, in- * t-suzuki@phys.aoyama.ac.jp cluding a sequence of studies by the author, have considered another interaction between thermal pressurization and dilatancy (referred to as ITPD below) [6,7,[9][10][11][12][13]. The dilatancy is slip-induced inelastic pore creation, leading to fluid pressure decrease and slip deceleration. This interaction enables us to treat slip acceleration and deceleration in a single framework, with which we can understand both ordinary and slow earthquakes.
However, some problems remain unsolved in the treatment of ITPD. For example, it should be emphasized that the porosity evolution laws are not yet firmly understood [14][15][16][17][18]. Laboratory experiments under conditions deep in the Earth are so hard to perform that exact reproduction of porosity evolution is considered difficult. A unified treatment independent of the details of the law is therefore required to understand the behavior of the system in the presence of ITPD. Additionally, geophysical studies have solely focused on explaining geophysical phenomena, and mathematically and physically important characteristics have not been investigated, as indicated below.
The studies of ITPD mentioned above employed nonlinear governing equation systems. In fact, the behaviors of solution orbits for nonlinear equation systems have been studied widely. The competitive Lotka-Volterra (LV) equation system, a model describing competition between two species, is an example of such a nonlinear equation system and has been a frequently treated topic recently [19][20][21]. For such a system, isolated fixed points on the phase space have been found, and the features of solution orbits crossing the points are well understood. The points are attractors, a saddle node and a repeller: the attractors describe extinction of one species, the saddle node corresponds to coexistence of both species, and the origin is the repeller. Discontinuous change of the solution behaviors depending on the initial value is observed there, similar to the discontinuous behaviors observed in the dynamic earthquake slip process. Nonetheless, two characteristics need to be emphasized in the ITPD model. First, continuous non-isolated fixed points appear for the ITPD model as attractors, as shown in this paper. Second, the initial values of the variables construct the continuous geometry in the ITPD model; the group of the initial values becomes a line in the phase space. We can therefore conclude that the initial and final values construct the continuous geometries in the ITPD model. We can expect a universal relation between these continuous geometries, although constructing the framework to treat such a relation has not been achieved.
This paper is organized as follows. The model setup and governing equation system are clarified in Sec. II. The slip velocity and inelastic porosity are the variables governing system behavior. Geometrically different attractors generated by the nullclines common to both variables are found mathematically in Sec. III. Discontinuity of the solution behaviors can be regarded as phase transition behavior. The criticality in the vicinity of the phase transition point is found in Sec. IV. The criticality is universal, and is unaffected by the assumed details of the porosity evolution law. Physical and seismological application of the results obtained is carried out in Sec. V. Dynamic earthquake slip behavior is concluded to be the phase transition phenomenon. In particular, predicting which phase emerges and the extent of the final slip amount is difficult. The paper is summarized, and nonlinear mathematical application is performed in Sec. VI.

II. MODEL SETUP
We consider a system consisting of a homogeneous and isotropic thermoporoelastic medium, i.e. the medium has pores whose volume ratio to the whole volume (porosity) is initially homogeneous. The initial porosity is referred to as elastic porosity, φ e . The pores are assumed to be filled with fluid (water). We also assume that the thermal pressurization and dilatancy effects emerge in the slip zone located at −w h /2 < y < w h /2 along the x-axis (Fig.  1); the relative movement between opposite surfaces is assumed to be accommodated entirely within the slip zone, which is regarded macroscopically as the one-dimensional (1D) mode III slip plane. The slip zone can be considered as a boundary of the medium from a macroscopic viewpoint. Thermal pressurization describes the elevation in fluid pressure based on the frictional heating associated with the dynamic fault slip [8]. If frictional slip occurs, frictional heating increases the temperature, inducing expansion of the solid and fluid phases. However, because it is easier to expand the fluid phase than the solid phase, the fluid pressure rises. On the other hand, dilatancy represents the inelastic porosity increase due to the fracturing of fault rocks by the fault slip, which reduces the fluid pressure [7]. It should be emphasized that when the fluid pressure increases (decreases), the frictional stress decreases (increases) due to a reduction (increment) in the effective normal stress acting on the fault plane, inducing the slip velocity increase (decrease). The competition between thermal pressurization and dilatancy induces complex feedback in the slip behavior, which can explain many aspects of the dynamic earthquake slip process (e.g., Suzuki and Yamashita (2014) [6], referred to as SY14 below). However, this system is not yet researched from the viewpoint of nonlinear mathematics, particularly with regard to the behaviors of attractors. This research therefore primarily aims at understanding the system with regard to the behaviors of attractors.
We have the following system of governing equations within the slip zone (SY14): Table I summarizes the meanings of the parameters. We neglected the advection effect and adiabatic expansion of the solid phase; the validity of this assumption has been demonstrated in our previous study [13]. We also neglected the diffusions of the fluid and heat, such that the Heaviside unit step function used in SY14 is not required and variables do not depend on y because the slip zone becomes an isolated system with such assumptions. We now consider the characteristic fluid diffusion time, t f ch , and roughly estimate it here to show that neglecting fluid diffusion can be valid in the interior of the Earth.
For the estimation, we should first mention that the slip zone thickness is known to be on the order of µm-cm [9,[22][23][24][25][26][27][28][29][30]. For subduction thrust faults, the order of mm-cm seems to be reasonable [25][26][27][28]. On the other hand, some researchers insist that the zones have width on the order of µm [29,30]. However, significant alongstrike variability may be observed in the localized zone thickness [30], and the fault zones accommodating displacements are considered to have complex geometrical structures [31]. We therefore assume an order of a few cm or less in this study. Additionally, the permeability k is known to be on the order of 10 −14 -10 −21 m 2 [32]. Let us estimate t f ch from the relation w h = kM t f ch /η, where η is the fluid phase viscosity. This relation is obtained because fluid diffusion can be described by adding a term (k/η)∂ 2 p D /∂y 2 to the right hand side of Eq. (1). When k is 10 −21 m 2 , w h = 3 mm and 3 cm give t f ch ∼ 8.54 × 10 s and 8.54 × 10 3 , respectively, whereas when k is 10 −14 m 2 , w h = 3 mm and 3 cm lead to t f ch ∼ 8.54 × 10 −6 s and 8.54×10 −4 , respectively (the values of M and η are listed in Table I).
We then consider the characteristic heat diffusion time, where λ is the thermal conductivity of the medium. This definition is reasonable because the term λ∂ 2 T D /∂y 2 added to the right hand side of Eq. (2) describes the heat diffusion.
Using the values shown in Table I and λ ∼ 1 J/mKs [33], we have t h ch ∼ 2.39 × 10 s and 2.39 × 10 3 s for w h = 3 mm and 3 cm, respectively.
We can show the conditions for neglecting fluid and heat diffusions based on the results of t f ch and t h ch . If we consider the time scale ≤ 10 s, neglecting heat diffusion is considered to be valid because t h ch > 10 s. Additionally, if we consider events with this time scale, we can conclude that when k is near the lower limit, fluid diffusion can be neglected (i.e., the system is undrained), whereas when k is near the upper limit, the diffusion effect should be taken into account (the system is drained). We can thus insist that the undrained assumption corresponds to that of lower k, and such low k is assumed henceforth. To summarize, the approximation of neglecting fluid and heat diffusions has been found to be valid for the slip duration ≤ 10 s with k near the lower limit. A duration lower than 10 s is characteristic time scale of earthquake duration with a moment magnitude ≤ 7, hence earthquakes with a moment magnitude ≤ 7 will be considered below. However, note that the model here developed is 1D; therefore, the moment magnitude cannot be exactly defined. This will be stated again in Sec. V.
Equations (1) and (2) govern the temporal evolution of the fluid pressure and temperature, respectively. Equation (3) provides a solution for the equation of motion (EOM) of the medium, with a boundary condition that difference between the applied shear stress and the frictional stress (stress drop, ∆σ) is given by σ 0 s −σ res [34]. In SY14, the displacement appearing in the EOM was that of the solid phase; nonetheless, the displacements of the solid and fluid phases are exactly the same in the present system because fluid diffusion is neglected. Equation (4) provides a definition for σ res . Note that the normal stress acting on the fault σ 0 n is negative because the compression stress is defined as a negative value. In particular, the first and second terms of the right hand side of Eq.
The framework here is constructed based on SY14. In fact, researchers have shown no agreement on the mathematical treatment of the motion of the thermoporoelastic medium, even without the dilatancy effect [35][36][37][38]. Despite the lack of agreement, the differences emerging in each framework do not cause any qualitative change in the system behavior because only the analytical forms of the coefficients appearing in the governing equations differ from one another. We can use the framework here as a general framework for treating ITPD.
To close the governing equation system, we also need the equation governing inelastic porosity evolution, which will be referred to as a porosity evolution law. Some forms of the law have been suggested based on many laboratory experiments, although there exists no agreement concerning its analytical form [14][15][16][17][18]. Therefore, we should derive robust results independent of the details of the porosity evolution law. For consistency with previous studies, we only assume that φ d is a function in terms of the slip, u D , and is initially zero, i.e., φ d = φ d (u D ) and φ d (0) = 0. However, we require some conditions on the form of φ d from the physical viewpoint. First, the function φ d (u D ) is assumed to monotonically increase with increasing u D since we do not consider pore healing in an entire single slip event; this implies that ∂φ d /∂u D is always nonnegative. Second, we also assume that lim uD →∞ ∂φ d /∂t = 0, because the relation 0 ≤ φ d ≤ 1 must be satisfied and φ d must have an upper limit, which we call φ UL henceforth. Note that φ UL need not be equal to unity. From these statements, we have the porosity evolution law for which the condition lim We can obtain the normalized governing equation system from Eqs. (1)(2)(3)(4)(5) where v (0 ≤ v ≤ 1) and u are normalized slip velocity and slip, respectively, φ is normalized inelastic porosity (0 ≤ φ ≤ 1), β is a positive constant, and f (u) is a function of u. The overdot stands for differentiation with respect to normalized time, τ . To derive Eqs. (6) and (7), note that v D and p D are linearly related (see Eqs. (3) and (4)), allowing us to rewrite p D in terms of v D .
The forms of u, v, φ, β, f (u) and τ are given by respectively. Based on the assumption for ∂φ d /∂u D , f (u) is a nonnegative function. We also have the condition lim u→∞ f (u) = 0, owing to the condition lim uD →∞ ∂φ d /∂u D = 0. Moreover, the parameter β describes the contribution of inelastic porosity increase to the slip velocity change, because from Eqs. (6) and (7) one In fact, the temporal evolution equation for the normalized temperature can be given bẏ where Nonetheless, this will not be handled in the investigation below because the governing equation system (6) and (7) is closed in terms of v and φ (note that u = vdτ ), and the temperature can be determined by We now describe φ in terms of u for convenience during later analytical treatments. From Eq. (7), we have where F −1 is an inverse function of F . Because f (u) is a nonnegative function, F (u) = f (u)du is a monotonically increasing function in terms of u. Therefore, F (u) clearly has its inverse function. Using Eq. (15), we have the following governing equationṡ Equations (16) and (17) form the general framework treating ITPD in the 1D model, and our analytical treatment will be based on these equations henceforth. In particular, we will derive a kind of phase transition and universal criticality emerging near the phase transition point.

A. Geometrically different attractors
We mathematically demonstrate that two geometrically different attractors emerge within the present framework including v and φ. To show this, we consider the qualitative behavior of the solution orbit in φ − v space (Fig. 2). We first consider nullclines, which are obtained by the conditionsv = 0 andφ = 0. Forv = 0, the straight line v = 0 and the curve v = 1 − βf (F −1 (φ))(≡ g(φ)) are nullclines from Eq. (16). The curve v = g(φ) on the φ − v space will be referred to as C crit henceforth. Forφ = 0, the straight line v = 0 and the curve f (F −1 (φ))(= (1 − g(φ))/β) = 0 are nullclines. Clearly, the line v = 0 is the nullcline for both equations. Moreover, we show here that the curve f (F −1 (φ)) = 0 can be described by the straight line φ = 1 in φ − v space. Note that the condition f (F −1 (φ)) = 0 corresponds to F −1 (φ) → ∞ based on the assumption lim u→∞ f (u) = 0; therefore, φ = lim u→∞ F (u) must be satisfied on such a nullcline. It should also be emphasized that (a) φ is a monotonically increasing continuous function in terms of u, and (b) φ is normalized to have a maximum value of unity. These statements suggest that the nullcline f (F −1 (φ)) = 0 is given by the straight line φ = lim u→∞ F (u) = 1. From these statements, we can also conclude that C crit is absorbed into point (1, 1). Additionally, we utilize an important condition here for C crit . We assume the relation This assumption and the statement that C crit is absorbed into point (1, 1) allow us to conclude that C crit crosses the φ-axis in the region 0 < φ < 1, i.e., at least one positive φ c satisfying is assumed to exist (see Fig. 2). The same condition has been treated in several previous studies (SY14, [9]). We require another important assumption for the solution orbit, namely that the orbits cross the v-axis in the region 0 < v < 1; that is, we use the condition 0 < v 0 < 1, where v 0 is the value of v at φ = 0, for the sake of physical and seismological applications in the latter part of this paper (the physically meaningful regions for φ and v are 0 ≤ φ ≤ 1 and 0 ≤ v ≤ 1, respectively). However, because the mathematical treatment is performed in the present subsection, the orbit can pass through the region v < 0.
First, we assume that g ′ (φ) ≥ 0 is satisfied for 0 ≤ φ ≤ 1, where the prime denotes differentiation with respect to φ. We also assume that if points satisfying g ′ (φ) = 0 exist, they are isolated. With these assumptions, Eq. (19) has a single positive root (which can be a multiple one). In this case, C crit is clearly right-upward and crosses the φ-axis once (Fig. 2). The sign of the gradient at the given point on the solution orbit depends on whether the point is on the upside or underside of C crit because dv/dφ = 0 on C crit . If the point is on the upside of C crit , the gradient is negative, whereas if it is on the underside of C crit , the gradient is positive. The orbit becomes horizontal at the point crossing C crit . Moreover, it should be noted that the orbit is neither horizontal nor vertical at the point crossing the line v = 0, even though the line is a nullcline. This occurs because v = 0 is a nullcline for both equations; the conditionv =φ = 0 is satisfied on the line v = 0, enabling dv/dφ to be nonzero and finite there. From these statements and Fig. 2, the orbit is found to connect the point (0, v 0 ) with (1, 1). In fact, the point (1, 0) is also a fixed one for Eqs. (16) and (17), and no orbit connecting (0, v 0 ) with (1, 0) exists, as shown in Appendix A. Figure 2 also shows arrows indicating the directions of evolution of the solutions with increasing time. These arrows can easily be obtained by the relationφ = f (F −1 (φ))v, because f (F −1 (φ)) is always positive except on the nullcline φ = 1, as noted in Sec. II, and the direction of the solution is determined only by the sign of v. If v > 0, the solution moves rightward, whereas if v < 0, it moves leftward with increasing time. The solution orbits and the arrows shown in Fig. 2 illustrate that we have an attractor and repeller on the φ-axis; where φ a and φ r are real numbers. In particular, note that φ a and φ r take continuous values. These non-isolated fixed points appear because the line v = 0 is a common nullcline for both equations, and this is a noteworthy behavior of the present system. In addition, the point (1, 1) is also an attractor because C crit and all the orbits are absorbed into the point (1, 1), and the solutions move rightward with increasing time where v > 0. We can therefore categorize the attractors into two geometrically different groups, which are given by the line {(φ a , 0)| 0 ≤ φ a ≤ φ c } or the point (1, 1) (see the green line and point in Fig.  2). We will refer to the former and latter attractors as line and point attractors, respectively. Exact analytical investigation determining which attractor appears will be carried out in Sec. IV.
In Fig. 2, we assumed that g ′ (φ) ≥ 0 for 0 ≤ φ ≤ 1 and that the curve C crit was right-upward. Here, we allow g ′ (φ) to be negative at a certain value of φ. Transitions from right-upward (right-downward) to right-downward (right-upward) with increasing φ (referred to as up−down (down−up) transitions henceforth) emerge with this condition. These cases correspond to the condition that g(φ) has maximal (for up−down transition) and minimal (for down−up transition) values. First, let us consider C crit with a single up−down transition and a single down−up transition. With this condition, the case where C crit crosses the φ-axis once is shown in Fig. 3(a). In Fig.  3(a), although an orbit crossing C crit three times appears, the attractors are the same as those observed in Fig. 2. Additionally, the case where C crit crosses the φaxis more than once is shown in Fig. 3(b). When there exist two or more φ c , we refer to the mth smallest one as φ m c , where m is a positive integer. When φ c is a multiple root, it is regarded as a single φ m c . In this case, the area where v < g(φ) (i.e., the underside of C crit ) cannot include attractors for the following reasons. The solution must move away from the φ-axis there because (a) dv/dφ > 0 is satisfied and (b) the solution describing infinitesimal perturbation from the φ-axis to the positive (negative) v direction moves rightward (leftward) with increasing time. Therefore, there exist two line attractors on the φ-axis, and the point (1, 1) is also an attractor. We can insist that even if C crit is not monotonically rising with increasing φ, two geometrically different attractors emerge and only the number of line attractors changes.
We should investigate the behavior of line attractors near transition points in detail. Let us first consider the case where C crit has a down−up transition. If the transition occurs under the φ-axis, C crit crosses the φ-axis and the line attractor appears (Fig. 4(a)). However, if the φ-axis becomes tangent to C crit , the attractor converges to a point ( Fig. 4(b)). This converged line attractor is a special case for the line attractor, and we do not refer to it as a point attractor. If the down−up transition occurs above the φ-axis, the attractor vanishes (Fig. 4(c)). We now treat the case in which C crit has an up−down transition. If such a transition occurs above the φ-axis, C crit crosses the φ-axis and two line attractors divided by C crit emerge in the vicinity of the transition point (Fig. 4(d)). However, if the φ-axis becomes tangential to C crit , the line attractors coalesce and turn into a single line attractor (Fig. 4(e)). If the up−down transition occurs under the φ-axis, a single line attractor will be observed ( Fig.  4(f)).
We can summarize the behaviors of the line attractors in terms of g and its derivatives. Let ∀j, g (j) (φ m c ) = 0 and g (n) (φ m c ) = 0, where g (j) stands for the jth derivative of g with respect to φ, n is a positive integer and j is a nonnegative integer satisfying 0 ≤ j ≤ n − 1 (note that g (0) = g). If n is odd and g (n) (φ m c ) < 0, the point (φ m c , 0) is the left end point of a line attractor, as φ 2 c in Fig. 3(b) and φ m c in Fig. 5(a). If n is odd and g (n) (φ m c ) > 0, the point (φ m c , 0) is the right end point of a line attractor, as φ 1 c and φ 3 c in Fig. 3(b) and φ m c in Fig. 5(b). If n is even and g  , 0), the point is neither the left nor the right end point. Finally, we can also emphasize that the emergence of two geometrically different attractors does not depend on the details of g(φ)(= 1 − βf (F −1 (φ))), i.e., the porosity evolution law. Only the number of line attractors varies owing to this detail.

B. Important concepts from the physical viewpoint
Here, we introduce some physically meaningful concepts associated with the current mathematical treatment. First, we introduce "phases" of the attractors for later discussion about phase transition. The line and point attractors physically correspond to lim τ →∞ v = 0 and 1, respectively. Hence, we can refer to the former and latter cases as the cessation and high-speed phases, respectively, from the physical viewpoint. The physical elementary process realizing these phases will be explained in Sec. V.
We now clarify the physically meaningful solution orbits. For the cases shown in Fig. 2, the orbit starts on the v-axis (φ = 0) physically. Within the orbits beginning with v 0 = v 01 or v 02 in Fig. 2, only the parts before absorption into the green line are physical, and the other parts are mathematical and unphysical. Additionally, for the cases shown in Fig. 3(b), the solution orbits crossing the green lines are physical before absorption, which is the same as in Fig. 2. However, a physically prohibited region (PPR) can emerge near φ i+1 left (Fig. 7). We define the point (φ i PPR , 0) as being where the orbit crossing the point (φ i right , 0) again crosses the φ-axis. If such an orbit is absorbed into the point attractor, we define φ i PPR = 1. The solutions cannot be absorbed into the region {(φ, 0)| φ i+1 left ≤ φ ≤ φ i PPR } with a physically meaningful initial condition and infinitely long time, and this region will be referred to as the PPR. The left end points of the line attractors other than the origin (φ 1 left , 0) must be accompanied by the PPR. If φ i+1 right < φ i PPR , the solutions cannot be absorbed into the i + 1th line attractor. To summarize, although all green lines in Fig. 6 actually represent line attractors in a mathematical sense, not all points on the attractors describe physically meaningful solutions.

IV. PHASE TRANSITION AND UNIVERSALITY
From the investigation above, we can conclude that a system including ITPD generates a phase transition between the cessation and the high-speed phases. We now show that universal (scale-independent) behavior emerges in the vicinity of the phase transition point in the present model by considering the perturbation in v 0 , which is similar to other phase transitions, e.g., susceptibility in phase transition of the second kind, and the spanning probability of percolation [39,40]. Such universality has many implications for the behavior of the final slip amount, which plays an important role in our understanding of the dynamic earthquake slip process, as shown in Sec. V.
To derive the universality, see the area near the origin in Fig. 2, which is enlarged in Fig. 8. In this case, C crit is assumed to cross the φ-axis once at (φ 1 c , 0) and g ′ (φ 1 c ) > 0, and we have the single solution φ c = φ 1 c = φ 1 right for Eq. (19). Note that we have a manifold dividing the point attractor and the line attractor (i.e., the basin boundary), which is drawn by the light blue curve in Fig. 8. We refer to the manifold as the critical manifold, and it begins at (0, v 1 right ) and ends at (φ 1 right , 0), where (0, v i right ) is a point where the manifold crossing (φ i right , 0) crosses the v-axis. If v 0 > v 1 right , the point attractor emerges, whereas if v 0 < v 1 right , the line attractor is realized. Let us consider here the manifolds in the neighborhood of the critical manifold. We assume that these manifolds are absorbed into the line attractor, such that the assumption v 0 < v 1 right is used here. In the φ − v space, if the initiation point of the manifold is given by (0, v 1 right − δv 1 + ), its ending point is expected to be described by (φ 1 right − δφ 1 + , 0), where δv 1 + and δφ 1 + are positive amounts that satisfy δv 1 + ≪ 1 and δφ 1 + ≪ 1. We will show that δφ 1 + and δv 1 + are related via a simple power law, and obtain the value of the critical exponent α in the following part of this section. Moreover, we will show that the details of the porosity evolution law do not affect the critical exponent.

A. Derivation of the power law
To show the power law (20), note that Eqs. (6) and (7) yield From Eq. (21), we have the solution for v based on the method of variation of constants in terms of φ, which leads to where at φ| τ =0 = 0 is used. Equation (22) gives the solution orbit in φ − v space. We then obtain the value of φ 1 right from Eq. (19). This equation reads as where g −1 is the inverse function of g. The value g −1 (0) can be uniquely defined within the current framework. Note that the character φ 1 right remains in the analytical treatment below to simplify the description because using g −1 (0) makes the representation complex. If we do not use φ 1 right explicitly, we should employ solution (23). Next, we obtain v 1 right from Eq. (22). This value can be derived from the condition v = 0 with v 0 = v 1 right and φ = φ 1 right in Eq. (22): which reads as   given by because we have the relations and (see also Eq. (23)). Equations (26) and (27) lead to ). Additionally, we have another important expansion in terms of δφ 1 + : Using Eqs. (30) and (31) and expanding the exponential function in Eq. (30), we find that the terms of orders (δφ 1 + ) 0 and (δφ 1 + ) 1 vanish, leaving us with (32) Neglecting the terms of order (δφ 1 + ) 3 and higher, we obtain a simple power law between δφ 1 + and δv 1 + : Note that we have assumed g ′ (φ 1 right ) > 0. The character φ 1 right vanishes by using Eq. (23), which leads to a power law relation described in terms of β and g(φ). We should emphasize that the power 1/2 is universal and does not depend on either β or the details of g(φ).
We also have the other right end points, φ i1 right (i 1 ≥ 2). Furthermore, we should treat the case g ′ (φ m c ) ≤ 0. However, as noted in Sec. III B, the end points other than φ 1 right and φ 1 left can be in the PPR depending on the form of C crit , and physically natural solutions may not approach their vicinity. In addition, the possibility that the φ-axis actually becomes tangential to C crit (i.e., g ′ (φ m c ) = 0) in natural fault is considered to be negligibly low. We will consider the region only near (φ 1 right , 0) from the physical and seismological viewpoints in Sec. V, and mathematical discussions about the other φ m c will be performed in Sec. VI.

B. Other important suggestions about the phase transition and the power law
We have another important conclusion about the parameters governing the phase transition based on the result obtained in this section. Note that the sign of v 0 − v 1 right is concluded to be important for the phase transition. If it is positive, the high-speed phase emerges, whereas if it is negative, the cessation phase is realized. Based on Eq. (25), we have the relation which includes the parameter β and the functions A(φ 1 right ) and B(φ 1 right ). These functions can be described by g(φ) by definition, and φ 1 right can be written in terms of g −1 (see Eq. (23)). The parameters β and v 0 as well as the function g(φ) are found to govern the phase transition, and the governing function is given by the right hand side of Eq. (34).
We can also obtain a relation similar to Eq. (33) by considering the value of u = vdτ . We consider the region near (φ 1 right , 0) here. Let us assume that u ∞ = u 1 right − δu 1 + for the cessation phase, where u ∞ ≡ lim τ →0 u is the final slip amount, u 1 right = u ∞ | v0=v 1 right and δu 1 + is the positive amount satisfying δu 1 + ≪ 1. The relation δφ 1 + = (1 − g(φ 1 right ))δu 1 + /β = δu 1 + /β (derived from Eqs. (17) and (23)) and Eq. (33) give the simple power law between δu 1 + and δv 1 + : Thus, we have a power law with the same critical exponent 1/2 as observed in Eq. (33). Equations (33) and (35) derived from the mathematical viewpoint imply that, on the phase space, the values on two line segments, one for the initial values and the other for the final values, are related via a simple power law from a physical viewpoint. Physical and seismological implications associated with this statement are given in Sec. V. Note that we only consider the region in the vicinity of the point (φ 1 right , 0), and write φ c , v c , u c , δv, and δu instead of φ 1 right , v 1 right , u 1 right , δv 1 right , and δu 1 right , respectively, below.

V. APPLICATION TO NATURAL FAULTS
Dynamic earthquake slip processes show phase transitions by considering φ ∞ or u ∞ as the order parameter from the conclusion in Sec. IV. This has some implications for the diversity observed in natural dynamic earthquake slip behavior. For example, the dependence of stress drops on earthquake size can be explained. The stress drop ∆σ is the difference between the applied shear stress acting on the fault plane and the residual frictional stress, as mentioned in Sec. II. Some researchers insist that large earthquakes sometimes have larger dynamic stress drops than other earthquakes [3]. This mechanism can be understood by the framework here. For such large earthquakes, the point attractor (high-speed phase) may be realized because ∆σ ≡ σ 0 s − σ res = µv/2β v = σ 0 s v (see Eqs. (3) and (9)) remains nonzero. This corresponds to the case wherein the acceleration by the fluid pressure increment due to the thermal pressurization effect completely governs the system behavior, and the shear stress acting on the fault plane is completely released owing to thermal pressurization. The fluid pressure approaches −σ 0 n in this case. On the other hand, other ordinary earthquakes may be realizations of the line attractor (cessation phase), because we clearly have ∆σ = 0 in the cessation phase. This behavior corresponds to a situation in which deceleration due to fluid pressure reduction induced by the dilatancy effect completely governs the slip behavior, and spontaneous slip cessation (v = 0) is realized. Note that the 1D system utilized here is an approximated one, and ∆σ remains nonzero near the fault tips for real three-dimensional systems. However, the near-tip area becomes negligibly small compared with the whole fault area with propagation of the tip, and the 1D approximation is expected to work well for natural faults.
The slip behavior from the onset to the attainment of two phases can be interpreted in terms of two physical processes, the thermal pressurization and dilatancy effects. Note that the dominant physical processes can exchange during the slip. If the solution orbit is on the upside (underside) of C crit , the gradient of the orbit is negative (positive) and v decreases (increases) with increasing time. The deceleration (acceleration) of the slip is physically observed, and the dilatancy (thermal pressurization) effect is dominant. The curve C crit , i.e., the function g(φ), determines which physical process is dominant. However, it should be emphasized that the highspeed and cessation phases are completely governed by the thermal pressurization and dilatancy effects, respectively, and an intermediate state at τ → ∞ does not exist, as mentioned above. Although this phase transition was also suggested in SY14, only a single form for the porosity evolution law was assumed, and the order parameter (φ ∞ or u ∞ ) was not clarified there.
The question here is whether we can predict which phase appears based on the physical viewpoint. As noted in Sec. IV B, β, v 0 , and g(φ) govern which phase emerges via the sign of relation (34). We can completely predict slip behavior mathematically. Nonetheless, depends upon the porosity evolution law, which has not been firmly understood, as mentioned in Sec. I. Thus, the law is so uncertain that we cannot predict which phase emerges from the physical viewpoint.
Though predicting which phase emerges is difficult, it is physically meaningful to assume that the cessation phase emerges because earthquakes with enormously large stress drops are rare [3]. This observational result implies that v 0 − v 1 right < 0 is satisfied for many earthquakes. The question arising from the assumption of the cessation phase is whether we can predict the final slip amount. Note here that the value of the final slip amount u ∞ = u c − δu is a measure of earthquake magnitude (which cannot be defined exactly here because the model is 1D), and that studying the behavior of u ∞ is important for understanding the dynamic earthquake slip process. We adopt the model of SY14 as an example, and perform numerical calculations. We will derive implications independent of β and g(φ).
The governing equations in SY14 are given bẏ where S u and T a are nondimensional positive parameters. Applying the notation in the present study to this equation system, we have f (F −1 (φ)) = T a (1 − φ), g(φ) = 1 − S u (1 − φ), and β = S u /T a . With these conditions and Eqs. (15), (23), and (25), we can write v c , φ c and u c for SY14: respectively. Finally, we can apply Eq. (35) to the model of SY14 to obtain The values of S u and T a to be utilized are those which are concluded to be appropriate for ordinary earthquakes in SY14. The Runge-Kutta method with the fourth-order accuracy is adopted. The calculated u was found to approach u ∞ within the numerical accuracy in finite time. Figure 9(a) shows the numerically obtained relation between δv and δu, which agrees with the result of Eq. (41). In addition, Figs. 9(b) and (c) show that the relative disturbance in the final slip amount, δu/u c , is several to 10 times larger than that in the initial slip velocity, δv/v c . This occurs because we have the universal critical exponent 1/2. Since the uncertainty in δv is significantly amplified in δu, predicting δu is almost impossible, and thus the final slip amount of the earthquake is hard to predict. This non-predictability is newly suggested here, and was not studied in previous studies including SY14.
Equation (35) shows that δu is proportional to δv 1/2 , and that the factor of proportionality depends on β and g(φ). Therefore, we should emphasize the nonpredictability as a universal conclusion independent of uncertainties of g(φ), and the non-predictability is not specific to the model of SY14.
The parameter δv = v c − v 0 depends on β, g(φ) and v 0 (see Eq. (34)), such that it is hard to evaluate it from the physical viewpoint. Moreover, δv has another uncertainty. Note that v 0 = 1 + µ slid (σ 0 n + p D0 )/σ 0 s , where p D0 is the initial fluid pressure. Among the parameters µ slid , σ 0 s , σ 0 n , and p D0 , note that p D0 is considered to be the most susceptible to the surrounding environment. The change of the fluid pressure within the fault rocks may be caused by some chemical processes such as the dehydration of hydrous minerals [41,42], hence its quantitative evaluation becomes very arduous, because it is a microscopic phenomenon. It can be concluded that the final slip amount may reflect the fluid pressure profile, and that its exact prediction is significantly difficult.

VI. DISCUSSION AND CONCLUSIONS
The present study provides a unified framework to treat thermal pressurization and dilatancy effects simultaneously. In particular, the nullclines common to the two equations are important in this framework. Such a system generates geometrically different attractors: the point attractor (high-speed phase in a physical sense) and the line attractor (cessation phase). The transition behavior between the two phases is observed near the basin boundary, and the universality in the vicinity of the transition point is represented by the power law between δφ and δv, regardless of the details of the porosity evolution law: δφ ∝ δv 1/2 . Dynamic earthquake slip processes can be regarded as phase transition phenomena, with the final inelastic porosity or the final slip amount as the order parameter. The prediction of the emergent phase is completely performed mathematically, whereas predicting it physically is difficult, mainly because the porosity evolution law has uncertainties. Moreover, even if we assume that the cessation phase emerges, the universality suggests the non-predictability of the final slip amounts of natural earthquakes.
If a complete framework for the porosity evolution law could be developed, we may predict the phase emergence by identifying parameters used to construct v 0 = 1 + µ slid (σ 0 A more detailed study of the porosity evolution law is important for understanding the phase emergence. However, it should be noted that as mentioned in Sec. V, evaluating v 0 is difficult. A stochastic approach could eventually make it easier, and the prediction would have a stochastic character. We should also emphasize that even in the case of high probability for the cessationphase emergence, the final slip amount could not be predicted because of the universality independent of β and g(φ).
Nonetheless, it is meaningful to investigate the effect of β because it is the single parameter emerging in the governing equation system. As mentioned in Sec. II, β is a measure of the contribution of the inelastic porosity increase to the slip velocity change. Considering the region 0 ≤ φ ≤ 1 and 0 ≤ v ≤ 1, and assuming that f is completely understood, it can be found that larger β generates larger area where v > g within such a region since g = 1 − βf ; in other words C crit shifts downward, but it keeps crossing the point (1, 1). Therefore, it can be concluded that larger β (e.g., larger M or smaller σ 0 s ) is more likely to generate the cessation phase. This means that if the porosity evolution law is understood, the tendency by which the cessation phase is likely to emerge can be estimated from β.
We can perform nonlinear mathematical applications of the framework constructed in this paper. In particular, we now discuss the case of g ′ (φ m c ) = 0 mathematically. Let us assume again that ∀j, g (j) (φ m c ) = 0, show the relation (see details in Appendix C). In this case, the universal critical power value is 1/(n + 1), which decreases with increasing n. In particular, if n is even and g (n) (φ m c ) < 0, φ m c −φ ∞ and v m c −v 0 can take both positive and negative values because the point (φ m c , 0) is on a line attractor and not its end point. In this case, Eq. (42) predicts that the region on the φ-axis near the point (φ m c , 0) is harder for the solution to approach with larger n, because the disturbance in δv (m) is enlarged beyond that in δφ (m) (note that δv (m) < 1 and δφ (m) < 1) and the enlargement is stronger, even though the region is on an attractor. These treatments are important from the viewpoint of nonlinear mathematics.
Finally, we can show that the treatments performed in this study can be extended to systems such as competition relations between two species (the LV model). The competitive LV model is a simple model of the population dynamics of species competing for some common resource. The framework constructed in the present article represents a special case for the system. To show this, first, note that the competitive LV model (in the absence of the diffusion terms) is given by the following equations: where x 1 (x 2 ) is the population size of species 1 (2), r 1 (r 2 ) is the inherent per-capita growth rate of 1 (2), K 1 (K 2 ) is the carrying capacity of 1 (2), and a 12 (a 21 ) represents the effect that species 2 (1) has on the population of species 1 (2). Let us consider the system under conditions that (A) the growth rate for species 2 is negligibly smaller than that for species 1, and (B) when species 1 consumes species 2, species 1 also dies, e.g., species 2 is poisonous. If we consider the limit r 2 → 0 while maintaining r 2 a 21 /K 2 constant, this system can be depicted and the governing equation system is given bẏ where X 1 ≡ x 1 /K 1 and X 2 ≡ x 2 /K 2 , and temporal differentiation is performed with respect to τ LV ≡ r 1 t. The system can be described by the equation system in the same manner as Eqs. (16) and (17). If we introduce the variables v ≡ X 1 and φ ≡ 1 − X 2 , Eqs. (45) and (46) are exactly the same as Eqs. (36) and (37), respectively, by replacing S u and T a with a 12 K 2 /K 1 and r 2 a 21 K 1 /r 1 K 2 , respectively. Which species survives is determined by two important values. The first one is v LV 0 , which is the v value at the point where the solution orbit crossing the point (φ init , v init ) crosses the v-axis on the phase space, where v init and φ init are the initial values for v and φ, respectively (note that v init and v LV 0 are different). The variables (v, φ, v 0 ) = (v init , φ init , v LV 0 ) must satisfy Eq. (22), and we have the relation Solving this equation for v LV 0 gives us The other value is v LV c , which is the v value at the point where the critical manifold crosses the v-axis. Replacing S u and T a with a 12 K 2 /K 1 and r 2 a 21 K 1 /r 1 K 2 , respectively, in Eq. (38) gives the value, which leads to We have expanded the function 1/h(φ * ) around φ * = 1, and C i is the coefficient of the ith order of 1 − φ * for the expansion. Moreover, the term e −β(A(φ)−A(0)) (v 0 − 1) clearly vanishes with the limit φ → 1, such that we can conclude that orbit (22) always crosses the point (1, 1) and does not cross (1, 0) when lim φ→1 A(φ) = ∞.
We now consider the condition where lim orbit crossing (1, v 1 ) describes the situation in which the porosity approaches the upper limit (φ = 1) within finite time. Then, the orbit ascends the straight line φ = 1 vertically, and eventually approaches the point (1, 1) (see Fig. 10). This describes the case where C crit is a tangent to the line φ = 1, and this case can be included in the investigation performed in this paper. We have shown that there exist no orbits connecting the points (0, v 0 ) and (1, 0) for any forms of A(φ). because terms lower than (δφ (m) ) n vanish. If n is odd and g (n) (φ m c ) < 0, the point (φ m c , 0) is the left end point and δφ (m) = −(φ m c − φ ∞ ) and δv (m) = −(v m c − v 0 ). If n is odd and g (n) (φ m c ) > 0, the point (φ m c , 0) is the right end point and δφ (m) = φ m c − φ ∞ and δv (m) = v m c − v 0 . If n is even and g (n) (φ m c ) > 0, a converged line attractor is observed at (φ m c , 0) and relation (42) does not appear because the region near the point cannot be a part of the line attractor. If n is even and g (n) (φ m c ) < 0, φ m c − φ ∞ and v m c −v 0 can accept both positive and negative values. We have obtained the relation between δφ (m) and δv (m) even in the case of n ≥ 2.